# BHAC-QGP: three-dimensional MHD simulations of relativistic heavy-ion
collisions
II. Application to Au-Au collisions

Markus Mayer<EMAIL_ADDRESS>Goethe University Frankfurt, Max-
von-Laue-Straße 1, 60438 Frankfurt am Main, Germany  Luciano Rezzolla  Goethe
University Frankfurt, Max-von-Laue-Straße 1, 60438 Frankfurt am Main, Germany
Frankfurt Institute for Advanced Studies, Ruth-Moufang-Straße 1, 60438
Frankfurt am Main, Germany  School of Mathematics, Trinity College, Dublin 2,
Ireland  Hannah Elfner  GSI Helmholtzzentrum für Schwerionenforschung,
Darmstadt, Germany  Goethe University Frankfurt, Max-von-Laue-Straße 1, 60438
Frankfurt am Main, Germany  Frankfurt Institute for Advanced Studies, Ruth-
Moufang-Straße 1, 60438 Frankfurt am Main, Germany  Helmholtz Research Academy
Hesse for FAIR, Max-von-Laue-Straße 12, 60438 Frankfurt am Main, Germany
Gabriele Inghirami  GSI Helmholtzzentrum für Schwerionenforschung, Darmstadt,
Germany  Dirk H. Rischke  Goethe University Frankfurt, Max-von-Laue-Straße 1,
60438 Frankfurt am Main, Germany  Helmholtz Research Academy Hesse for FAIR,
Max-von-Laue-Straße 12, 60438 Frankfurt am Main, Germany

###### Abstract

We present BHAC-QGP, a new numerical code to simulate the evolution of matter
created in heavy-ion collisions. BHAC-QGP is based on the Black Hole Accretion
Code (BHAC), which has been designed to model astrophysical processes through
the solution of the equations of general-relativistic magnetohydrodynamics.
Like the mother code, BHAC-QGP uses Adaptive Mesh Refinement (AMR), which
allows for a dynamic adjustment of the resolution in regions of the
computational domain where a particularly high accuracy is needed. We here
discuss a number of applications of BHAC-QGP to Au-Au collisions at
Relativistic Heavy-Ion Collider (RHIC) energies and show that the code is able
to reproduce results of other simulations of these scenarios, but with much
higher accuracy.

##  I Background

In non-central relativistic heavy-ion collisions, the presence of charged
protons leads to the generation of immensely strong electromagnetic fields [1,
2, 3, 4, 5, 6, 7, 8, 9, 10, 5, 11]. Their magnitude in the laboratory frame
can be naively estimated with the help of the Liénard-Wiechert potentials as:

$eB\approx\dfrac{e^{2}}{4\pi}Zv\gamma\dfrac{8}{b^{2}}\;.$ (1)

Here, $Z$ represents the charge number and $R$ the radius of the nucleus,
$e\approx 0.303$ denotes the electric charge, and
$\gamma=\left(1-v^{2}\right){}^{-1/2}$ is the Lorentz factor of the nucleus in
motion at velocity $v$. This Lorentz factor is related to the center-of-mass
energy $\sqrt{s_{\text{NN}}}$ via $\gamma=\sqrt{s_{\text{NN}}}/(2m_{N})$, with
$m_{N}$ being the mass of the nucleon. Furthermore, $b$ in Eq. (1) corresponds
to the impact parameter, which is defined as the magnitude of the vector
connecting the centers of the two colliding nuclei. Typically, the impact
parameter establishes the direction of the $x$-axis. Together with the beam
direction, which is usually taken along the $z$-axis, the impact parameter
forms the so-called reaction plane. Notably, the generated magnetic field is
oriented mainly in the direction perpendicular to the reaction plane.

The Relativistic Heavy-Ion Collider (RHIC) at the Brookhaven National
Laboratory (BNL) collides Au nuclei with an energy of
$\sqrt{s_{\text{NN}}}=200~{}\text{GeV}$. According to Eq. (1), for an impact
parameter of approximately $b\sim 10~{}\text{fm}$ such a collision produces a
magnetic field of the order of $eB\approx 10~{}m_{\pi}^{2}\approx 6\times
10^{18}~{}\text{G}$, where $m_{\pi}\approx 0.14~{}\text{GeV}$ is the mass of a
pion [12]. Hence, the magnetic field strength exceeds not only the squared
pion mass but also the squared masses of up and down quarks. Such strong
fields can induce substantial quantum effects in the produced quark-gluon
plasma (QGP). Notably, due to the linear dependence on the center-of-mass
energy, the magnetic fields generated in Pb-Pb collisions at the Large Hadron
Collider (LHC) at CERN are even stronger. Moreover, the electric field
generated by the relativistically fast moving charges is also enhanced by a
Lorentz gamma factor and therefore equally strong as the magnetic field (1).
These enormously strong electromagnetic fields are believed to have a
considerable influence on the particles within the QGP [13, 14]. Given the
extreme densities and temperatures within the QGP, individual quarks are no
longer confined within hadrons and can move freely and interact with each
other. At such conditions, the masses of light quarks become negligible,
resulting in an approximate chiral symmetry in the Lagrangian of Quantum
Chromodynamics (QCD), the theory of strong interactions. Lattice-QCD
calculations have shown that this chiral-symmetry restoration occurs at
temperatures $T>170~{}\mathrm{MeV}$ [15, 16]. Topological fluctuations can
lead to a chiral imbalance in the medium [17], which combined with strong
external electromagnetic fields lead to a wide range of novel effects,
including the Chiral Magnetic Effect (CME), Chiral Separation Effect (CSE),
and Chiral Magnetic Wave (CMW) (see, for instance, Refs. [18, 19, 20, 12]).
The CME is caused by an imbalance of left- and right-handed charges and
results in a vector charge current of opposite polarization charges along the
direction of the external magnetic field [21, 22]. The CSE exists in a medium
with a nonzero vector charge, resulting in an axial-vector current along the
direction of the magnetic field. The interaction between CME and CSE causes
the CMW [23]. Another effect which occurs in a globally rotating system of
chiral fermions (i.e., in non-central heavy-ion collisions) is the Chiral
Vortical Effect (CVE), which is the generation of an axial-vector current
along the rotation axis. Such an axial-vector current is converted into a spin
polarization of quarks in the QGP, with fermions having their spins preferably
aligned with the vorticity [24, 25, 26, 27, 28].

The key to an experimental observation of the CME, CSE, and CMW lies in the
strength of the electromagnetic fields. If these fields are too weak prior to
the thermalization of the fireball, they may not lead to significant effects.
In fact, the lifetime of the electromagnetic field represents a major
uncertainty. In the simplest way it can be estimated as $t_{B}\sim R/(\gamma
v)$ [6, 12]. Since the spectators quickly leave the collision zone, the
remaining fireball plays an important role because its composition can
significantly prolong the decay of the electromagnetic fields. If a charged
QGP forms rapidly, it will respond to the fields and modify their evolution.
This phenomenon is attributable to Faraday induction, where an external
magnetic field induces an electric current within the medium, thereby
generating a magnetic field that counteracts the decaying external field.

In order to understand the dynamics of the electromagnetic fields in heavy-ion
collisions and their effects on chiral phenomena, a comprehensive
understanding of the entire evolution of the produced matter is essential.
Relativistic hydrodynamics serves as a successful framework for modeling the
evolution of the QGP near thermodynamical equilibrium. Second-order
dissipative relativistic hydrodynamics has proven to be capable of reproducing
key observables, such as the transverse-momentum spectra and the elliptic
flow. These measurements suggest that the medium is strongly interacting,
leading to a fast thermalization of the fireball and a strongly interacting
nature [29, 30, 31, 32, 33]. Incorporating electromagnetic fields into
hydrodynamics leads to a theory of magnetohydrodynamics (MHD).

With the help of numerical MHD simulations, it is possible to calculate the
(3+1)-dimensional evolution of the fireball under the influence of these
electromagnetic fields, enabling the exploration of uncertainties and chiral
effects. This paper presents the results of numerical simulations of heavy-ion
collisions performed with BHAC-QGP, a new code capable of describing the
space-time evolution of the QGP under the influence of magnetic fields.

This paper is organized as follows: We begin by providing a brief summary of
the numerical formalism of BHAC-QGP. Next, we present the procedures
implemented for initializing heavy-ion collisions within BHAC-QGP,
specifically explaining the optical Glauber model and the equations used for
calculating initial magnetic fields. Afterwards, we explain the application of
the Cooper-Frye prescription, which we employ to evaluate the freeze-out of
particles from the fluid evolution. From the Cooper-Frye prescription, one can
then compute particle spectra and anisotropic flow coefficients. Following
this introduction, we present the results derived from simulations performed
with BHAC-QGP.

Throughout this article we use natural Heaviside-Lorentz units,
$\hbar=c=k_{B}=\epsilon_{0}=\mu_{0}=1$. In these units the electric charge is
dimensionless $e:=\sqrt{4\pi\alpha\hbar c}\approx 0.303$, where $\alpha\approx
1/137$ is the fine-structure constant. As signature of the metric tensor, we
choose $\left(-,+,+,+\right)$. We use Greek indices to indicate the components
of a four-vector, while Latin indices range from 1 to 3 for its spatial
components. Bold letters indicate three-vectors. The scalar and vector
products between three-vectors are as usual denoted by $\cdot$ and $\times$,
respectively. Furthermore, we will use the Levi-Cività tensor with
$\epsilon^{0123}=1$. In this article we will work with two coordinate systems,
Cartesian coordinates $(t,x,y,z)$ and Milne coordinates $(\tau,x,y,\eta_{S})$,
where $\tau$ is the proper time and $\eta_{S}$ the so-called space-time
rapidity.

##  II Formalism

This section is dedicated to our numerical setup of heavy-ion collisions. Our
goal is to simulate the evolution of the matter produced in relativistic
heavy-ion collisions under the influence of external magnetic fields. BHAC-QGP
is capable to solve the equations of ideal general-relativistic MHD with any
time-independent metric. For the purpose of calculating heavy-ion collisions,
we implemented time-dependent Milne coordinates, because these are best suited
to describe the longitudinal expansion of a fluid, cf. Fig. 1. BHAC-QGP
benefits from the underlying parallel Adaptive Mesh Refinement (AMR) approach,
that allows for a dynamic adjustment of the numerical resolution. Because the
magnetic fields in heavy-ion collisions also exist outside the collision zone,
BHAC-QGP has to deal with regions where the magnetic pressure is much stronger
than the fluid pressure. One-dimensional (1D) and two-dimensional (2D)
inversion schemes are likely to fail in such cases. Fur such scenarios, BHAC-
QGP is equipped with a so-called “entropy switch”, such that in highly
magnetized regions, BHAC-QGP can use the entropy evolution equation to find a
valid solution. More details on the set of equations solved and on the
numerical techniques employed for their solution can be found in a companion
paper, hereafter referred to as Paper I.

Figure 1: (1+1)-dimensional Milne diagram of a relativistic heavy-ion
collision. At $t=0~{}\text{fm}$ the two ions collide at $z=0~{}\text{fm}$ and
produce an extremely hot fluid in the collision region. This region approaches
thermodynamical equilibrium after $\tau_{0}\sim 1~{}\text{fm}$, so that it can
be described with hydrodynamics. This description is possible up to a certain
freeze-out time, $\tau_{\text{f.o.}}\sim 10~{}\text{fm}$, after which the
fluid becomes a gas of free-streaming particles.

###  II.1 Initial energy density - Glauber model

The initial energy-density distribution of the matter produced in heavy-ion
collisions is computed using the optical Glauber model [34, 35, 36, 37, 38,
39, 40, 41, 42], which is also used in many other codes (see, for example,
Refs. [38, 39, 40, 41, 34, 42]). The basis of this model are two ions,
projectile $A$ and target $B$, colliding at relativistic velocities and
arbitrary impact parameter. The Glauber model establishes a relationship
between the impact parameter $b$, the number of participants
$N_{\text{part}}$, and the number of collisions $N_{\text{coll}}$. A small
impact parameter corresponds to a central collision, while collisions with a
large impact parameter are called non-central or peripheral.

For estimating of $N_{\text{part}}$ and $N_{\text{coll}}$, the model requires
remarkably only two inputs: an inelastic nucleon-nucleon cross section
$\sigma_{\text{inel}}$ and an assumption on the shape of the colliding nuclei.
The inelastic cross section is usually taken from experimental measurements,
while for the nucleus one commonly assumes a smooth mass density described by
a Woods-Saxon distribution (see, for instance, Ref. [41]),

$\rho:=\dfrac{\rho_{0}}{1+\text{exp}\\!\left[\frac{\sqrt{x^{2}+y^{2}+z^{2}}-R}{d}\right]}\;,$
(2)

where $R$ is the radius of the nucleus, which is related to the mass number
$A$ via $R\approx(1.12~{}\text{fm})\,A^{1/3}-0.86\,A^{-1/3})$ [43], while $d$
is the surface thickness of the nucleus. The parameter $\rho_{0}$ represents a
normalization factor defined to ensure $\int\rho\\!\left(r\right)d^{3}r=A$.
With the Woods-Saxon distribution (2) it is then possible to evaluate the so-
called nuclear profile function, which corresponds to the probability per unit
transverse area of a nucleon being located in a projectile flux tube,

$\hat{T}_{\text{A,B}}\\!\left(x,y\right):=\int\limits_{-\infty}^{\infty}\hat{\rho}_{\text{A,B}}\\!\left(x,y,z\right)dz\;.$
(3)

Here, $\hat{\rho}_{\text{A,B}}:=\rho_{\text{A,B}}/\rho_{\text{A,B};0}$ is the
normalized Woods-Saxon distribution of either projectile $A$ or target $B$.
Because the impact-parameter vector $\boldsymbol{b}$ defines the $x$-axis,
$\boldsymbol{b}:=\left(b,0\right)^{T}$, the nuclear profile function actually
depends only on the magnitude of the impact parameter,
$\hat{T}_{\text{A,B}}\\!\left(x,y\right)=\hat{T}_{\text{A,B}}\\!\left(b\right)$.
The integration over the joint probability then defines the so-called nuclear
thickness function,

$\hat{T}_{\text{AB}}:=\int\limits_{-\infty}^{\infty}\\!\int\limits_{-\infty}^{\infty}\hat{T}_{\text{A}}\\!\left(x+\dfrac{b}{2},y\right)\hat{T}_{\text{B}}\\!\left(x-\dfrac{b}{2},y\right)dx\,dy\;.$
(4)

The nuclear thickness function $\hat{T}_{\text{AB}}$ can be interpreted as
effective overlap area for which a specific nucleon in projectile $A$
interacts with a given nucleon in target $B$. The probability of a collision
between two particles in the transverse plane, that is the number density of
collisions, is then simply the product of nuclear thickness function and
inelastic nucleon-nucleon cross section,

$n_{\text{coll}}=\sigma_{\text{inel}}\,\hat{T}_{\text{AB}}\;,$ (5)

so that the total number of collisions has a strong dependence on the impact
parameter. In addition to $n_{\text{coll}}$, the number density of
participating nucleons $n_{\text{part}}$ can also be estimated as:

$\displaystyle n_{\text{A};\text{part}}\left(x,y;b\right)=$ $\displaystyle
A\,\hat{T}_{\text{A}}\\!\left(x+\frac{b}{2},y\right)\left\\{1-\left[1-\hat{T}_{\text{B}}\\!\left(x-\frac{b}{2},y\right)\sigma_{\text{inel}}\right]^{B}\right\\}\;,$
(6) $\displaystyle n_{\text{B};\text{part}}\left(x,y;b\right)=$ $\displaystyle
B\,\hat{T}_{\text{B}}\\!\left(x-\frac{b}{2},y\right)\left\\{1-\left[1-\hat{T}_{\text{A}}\\!\left(x+\frac{b}{2},y\right)\sigma_{\text{inel}}\right]^{A}\right\\}\;.$
(7)

The total density of participants is then $n_{\text{part}}\
=n_{\text{A};\text{part}}+n_{\text{B};\text{part}}$. Since experimental
measurements indicate that the deposited energy density scales with the number
density of participants [44, 45], we follow Refs. [38, 40] and parametrize the
initial energy density in the transverse plane (at $z=\eta_{S}=0$) according
to the two-component model as:

$\displaystyle e_{\mathrm{2D}}\\!\left(x,y;b\right)=$ (8) $\displaystyle
e_{0}\Bigg{[}\dfrac{\alpha_{H}n_{\text{coll}}\\!\left(x,y;b\right)+\left(1-\alpha_{H}\right)n_{\text{part}}\\!\left(x,y;b\right)}{\alpha_{H}n_{\text{coll}}\\!\left(0,0;0\right)+\left(1-\alpha_{H}\right)n_{\text{part}}\\!\left(0,0;0\right)}\Bigg{]}\;,$

where $e_{0}:=e\\!\left(0,0,0\right)$, i.e., the energy density $e$ at the
origin for $b=0$. The parameter $\alpha_{H}\in\left[0,1\right]$ is the so-
called collision hardness and allows for a scaling between hard and soft
processes. Hard processes are high-energy interactions between quarks and
gluons that produce particles such as jets, leptons, or pions with large
transverse momenta, $p_{T}>1.0~{}\text{GeV}$, which are neglected when
$\alpha_{H}=0$. Soft processes are characterized by longer time scales and
smaller momentum transfers, which occur, for example, in hadronization or
parton energy loss. Experimental measurements indicate that soft processes
dominate in heavy-ion collisions [44]. This means that most of the
interactions involve little transverse-momentum transfer, such that
$\alpha_{H}$ is usually relatively small. In summary, $e_{0}$ and $\alpha_{H}$
are the two tunable parameters of the two-dimensional Glauber model (8), which
can be fixed by comparing the particle distributions from the hydrodynamical
model with experimental data.

Three-dimensional simulations of heavy-ion collisions have to take into
account also the longitudinal dependence of the initial energy density. A
widely used parameterization [40, 38, 45, 46, 47] is based on a flat rapidity
region around $\eta_{S}=0$ and a Gaussian damping in forward and backward
direction,

$\\!\\!H\\!\left(\eta_{S}\right)=\text{exp}\\!\left[-\Theta\\!\left(|\eta_{S}|-\frac{\Delta_{\eta_{S}}}{2}\right)\dfrac{\left(|\eta_{S}|-\Delta_{\eta_{S}}/2\right)^{2}}{2\sigma_{\eta_{S}}^{2}}\right].$
(9)

where the tunable parameter $\Delta_{\eta_{S}}$ controls the width of the flat
rapidity region, while the parameter $\sigma_{\eta_{S}}$ sets the damping of
the Gaussian fall-off. The full parameterization of the energy density reads
then as:

$\displaystyle e_{\text{3D}}$
$\displaystyle\left(x,y,\eta_{S};b\right)=\dfrac{e_{0}}{w}\,H\\!\left(\eta_{S}\right)\,\Theta\\!\left(Y_{B}-|\eta_{S}|\right)$
$\displaystyle\times$
$\displaystyle\Bigg{[}\alpha_{H}n_{\text{coll}}\\!\left(x,y;b\right)+\left(1-\alpha_{H}\right)$
$\displaystyle\times\Bigg{(}\dfrac{Y_{B}-\eta_{S}}{Y_{B}}n^{\text{A}}_{\text{part}}\\!\left(x,y;b\right)+\dfrac{Y_{B}+\eta_{S}}{Y_{B}}n^{\text{B}}_{\text{part}}\\!\left(x,y;b\right)\Bigg{)}\Bigg{]}\;.$

The Heaviside function guarantees that initially the matter does not move
faster in longitudinal direction than with beam rapidity, i.e.,
$|\eta_{S}|\leq Y_{B}$, where the beam rapidity is defined as:

$Y_{B}:=\dfrac{1}{2}\,\text{ln}\\!\left(\dfrac{1+v}{1-v}\right)\;.$ (11)

The energy density in Eq. (II.1) is weighted by a factor

$w=\left(1-\alpha_{H}\right)n_{\text{part}}\\!\left(0,0;0\right)+\alpha_{H}n_{\text{coll}}\\!\left(0,0;0\right)\;.$
(12)

Figure 2: Initial energy-density distribution of a Au-Au collision at RHIC
energy $\sqrt{s_{\text{NN}}}=200~{}\text{GeV}$, calculated with the
(3+1)-dimensional optical Glauber model. The impact parameter of
$b=8~{}\text{fm}$ chosen here corresponds to a peripheral collision, which is
why the fireball has an almond-like shape in the $x$-$y$-plane (at
$\eta_{S}=0$).

###  II.2 Initial magnetic field

Since nucleons are positively charged, they generate an electric current when
in motion, which subsequently produces a magnetic field according to Maxwell’s
equations [9, 14]. Because heavy-ion collisions take place at almost the speed
of light, the resulting magnetic fields also reach enormous strengths.
Sufficiently strong magnetic fields may have consequences on the dynamics of
the produced QGP and on the particles itself. Therefore, to ultimately study
the effect of the electromagnetic fields, a realistic initial parametrization
of these fields produced during the collision of the two positively charged
ions is required. Earliest calculations indicated that the electromagnetic
field could decay and even disappear so quickly that it may not have a
considerable effect on the QGP (see, e.g.,Ref. [48]). However, the situation
is quite different when the medium is electrically conducting. Indeed, a QGP
with a large electric conductivity will have a significant response to the
change of the external magnetic field, delaying the decay of the
electromagnetic field. This is simply the consequence of the induced currents:
An external magnetic field will induce an electric current in the medium
which, in turn, will cause a magnetic field compensating the decaying external
magnetic field. Lattice-QCD calculations predict that the electric
conductivity of a static QGP is relatively large [49, 50, 51, 9, 11],

$\sigma_{E}=\left(5.8\pm 2.9\right)\dfrac{T}{T_{C}}~{}\text{MeV}\;,$ (13)

where $T$ is the plasma temperature and $T_{C}$ the transition temperature
between QGP and hadronic phase. However, these calculations are done in the
quenched approximation (without dynamical quarks) and assume that the medium
itself is static. Furthermore, a uniform but temperature-dependent electric
conductivity is clearly an approximation since the temperature in a collision
varies rapidly in spacetime.

As described in Ref. [9], it is possible to solve Maxwell’s equations for an
infinite homogeneous medium that has a constant electric conductivity
$\sigma_{E}$ and a constant chiral magnetic conductivity
$\sigma_{\chi}\ll\sigma_{E}$. The magnetic field created by one particle with
electrical charge $e$ moving in $z$-direction with velocity $v$ is a solution
of the Maxwell equations, which in cylindrical coordinates takes the form:

$\displaystyle B_{\rho}\\!\left(t,\boldsymbol{x}\right)$
$\displaystyle=-\sigma_{\chi}\dfrac{e}{8\pi}\dfrac{v\gamma^{2}x_{T}}{\Delta^{3/2}}\left[\gamma\left(vt-z\right)+A\sqrt{\Delta}\right]e^{A}\;,$
(14) $\displaystyle B_{\phi}\\!\left(t,\boldsymbol{x}\right)$
$\displaystyle=\dfrac{e}{4\pi}\dfrac{v\gamma
x_{T}}{\Delta^{3/2}}\left(1+\dfrac{\sigma_{E}v\gamma}{2}\sqrt{\Delta}\right)e^{A}\;,$
(15) $\displaystyle B_{z}\\!\left(t,\boldsymbol{x}\right)$
$\displaystyle=\sigma_{\chi}\dfrac{e}{8\pi}\dfrac{v\gamma}{\Delta^{3/2}}\Big{[}\gamma^{2}\left(vt-z\right)^{2}$
(16)
$\displaystyle\times\left(1+\dfrac{\sigma_{E}v\gamma}{2}\sqrt{\Delta}\right)+\Delta\left(1-\dfrac{\sigma_{E}v\gamma}{2}\sqrt{\Delta}\right)\Big{]}e^{A}\;,$

where $x_{T}^{2}:=x^{2}+y^{2}$,
$\Delta:=\gamma^{2}\left(vt-z\right)^{2}+x_{T}^{2}$ and
$A:=\left(\sigma_{E}v\gamma/2\right)\left[\gamma\left(vt-z\right)-\sqrt{\Delta}\right]$.
The numerical integration over all charged particles is explained in detail in
Ref. [3]. With the nuclear charge density $\rho=Z/(\frac{4}{3}\pi R_{A}^{3})$,
the magnetic fields of the projectile nucleus moving in $+z$-direction and the
one of the target nucleus moving in $-z$-direction are evaluated as:

$\displaystyle\boldsymbol{B}_{\text{proj}}\\!\left(x_{-},|\boldsymbol{b}_{1}|\right)$
$\displaystyle=$ $\displaystyle\int
2\,\rho\,\sqrt{R^{2}-{b^{\prime\,2}}}\,\,\boldsymbol{B}\\!\left(x_{-},|\boldsymbol{b}_{1}-\boldsymbol{b}^{\prime}|\right)\cdot\left(-\sin\psi_{1}\boldsymbol{e}_{x}+\cos\psi_{1}\boldsymbol{e}_{y}\right)d^{2}b^{\prime}\;,$
(17)
$\displaystyle\boldsymbol{B}_{\text{targ}}\\!\left(x_{+},|\boldsymbol{b}_{2}|\right)$
$\displaystyle=$ $\displaystyle\int
2\,\rho\,\sqrt{R^{2}-{b^{\prime\,2}}}\,\,\boldsymbol{B}\\!\left(x_{+},|\boldsymbol{b}_{2}-\boldsymbol{b}^{\prime}|\right)\cdot\left(-\sin\psi_{2}\boldsymbol{e}_{x}+\cos\psi_{2}\boldsymbol{e}_{y}\right)d^{2}b^{\prime}\;,$
(18)

where $x_{\mp}:=t\mp z/v$, while $\psi_{1,2}$ is the angle between the vector
$\boldsymbol{b}_{1,2}-\boldsymbol{b}^{\,\prime}$ and the $x$-axis,

$\cos\psi_{1,2}=\dfrac{b_{1,2}\cos\phi_{1,2}-b^{\prime}\cos\phi^{\prime}}{\sqrt{b_{1,2}^{2}+b^{\prime\,2}-2\,b_{1,2}\,b^{\prime}\cos\\!\left(\phi^{\prime}-\phi_{1,2}\right)}}\;.$
(19)

We note that the total magnetic field is then simply the sum of the two,
$\boldsymbol{B}_{\text{tot}}=\boldsymbol{B}_{\text{proj}}+\boldsymbol{B}_{\text{targ}}$,
and while the transformation from cylindrical to Cartesian coordinates is
straightforward, the transformation from Cartesian to Milne coordinates reads
as [14]:

$\tilde{B}_{x}=\dfrac{B_{x}}{\cosh\eta_{S}}\;,\quad\tilde{B}_{y}=\dfrac{B_{y}}{\cosh\eta_{S}}\;,\quad
B_{\eta_{S}}=\dfrac{B_{z}}{\tau}\;.$ (20)

Figure 3: Initial magnetic field in the transverse plane (at $\eta_{S}=0$) as
it is produced in a Au-Au collision at a center-of-mass energy
$\sqrt{s_{\text{NN}}}=200~{}\text{GeV}$ and an impact parameter of
$b=10~{}\text{fm}$. The top row shows the magnetic field for a medium with
$\sigma_{E}=5.80~{}\text{MeV}$, $\sigma_{\chi}=0~{}\text{MeV}$, while the
bottom row represents the magnetic field calculated for a medium with
$\sigma_{E}=5.80~{}\text{MeV}$, $\sigma_{\chi}=1.50~{}\text{MeV}$. Clearly,
the symmetry under reflection at the $y-$axis is broken in case of a
nonvanishing value of $\sigma_{\chi}$, which is a consequence of the current
that is produced due to the chiral imbalance. This current in turn generates
the symmetry-breaking magnetic field.

###  II.3 Cooper-Frye prescription

The fireball produced in a relativistic heavy-ion collision will approach
local thermodynamical equilibrium through multiple interactions of the
constituent particles. Consequently, a thermodynamical-equilibrium pressure
with respective gradients is established, which drives the expansion of the
fluid into the surrounding vacuum, accompanied by cooling. Once the
temperature falls below the transition temperature between the QGP and the
hadronic phase, the system hadronizes. Upon further cooling, the hadrons cease
to interact and propagate as free-streaming particles. Hydrodynamics is then
no longer applicable, and the system freezes out. The freeze-out of particles
from the fluid is commonly accomplished using the widely adopted Cooper-Frye
prescription [52]. The Cooper-Frye method essentially counts the number of
particle worldlines passing through the freeze-out hypersurface, which is
commonly taken to be a surface of constant temperature or energy density. The
freeze-out hypersurface is a three-dimensional surface
$\Sigma\\!\left(x\right)$ in four-dimensional space-time. This hypersurface is
subdivided into infinitesimal elements $d^{3}\sigma$, each associated with an
outward-pointing four-vector $d^{3}\sigma_{\mu}\\!\left(x\right)$
perpendicular to $\Sigma$ at point $x$. The number of particle worldlines
passing through this surface is determined by the expression:

$\displaystyle N_{i}$ $\displaystyle=$
$\displaystyle\dfrac{g_{i}}{\left(2\pi\hbar\right)^{3}}\int\limits_{\Sigma}\\!\int\dfrac{d^{3}p}{p^{0}}\,f_{i}\\!\left(x,p\right)p^{\mu}d^{3}\sigma_{\mu}\;,$
(21)

where $g_{i}$ is the degeneracy factor that accounts for the internal number
of degrees of freedom (i.e., spin, isospin, etc.) of particles with mass
$m_{i}$ and chemical potential $\mu_{i}$, while $p^{0}\equiv
E_{i}:=\sqrt{\mathbf{p}^{2}+m_{i}^{2}}$ is the on-shell energy of these
particles. If the fireball is in local thermodynamical equilibrium, the
Lorentz-invariant local-equilibrium distribution $f_{i}\\!\left(x,p\right)$ is
given by:

$f_{i}\\!\left(x,p\right)=\dfrac{1}{e^{\beta\left(x\right)\big{[}p^{\nu}u_{\nu}\left(x\right)-\mu_{i}\left(x\right)\big{]}}+a}\;.$
(22)

Here, $a$ accounts for the proper quantum statistics of particles species $i$:
$a=-1$ for bosons, $a=+1$ for fermions, and $a=0$ for Boltzmann particles. The
factor $p^{\nu}u_{\nu}\\!\left(x\right)$ represents the energy of the particle
in the local rest frame, boosted to the observer’s frame by the fluid four-
velocity $u^{\mu}$, while $\beta(x):=1/T(x)$ is the inverse temperature. The
fluid four-velocity, the temperature, and the chemical potential on the
freeze-out hypersurface are obtained from the hydrodynamical calculation.
Dividing Eq. (21) by the Lorentz-invariant momentum-space measure
$d^{3}p/E_{i}$ yields the Cooper-Frye formula for particle species $i$:

$E_{i}\,\dfrac{dN_{i}}{d^{3}p}=\dfrac{g_{i}}{\left(2\pi\hbar\right)^{3}}\int\limits_{\Sigma}f_{i}\\!\left(x,p\right)p^{\mu}d^{3}\sigma_{\mu}\;.$
(23)

Evaluating Eq. (23) requires finding the location of the hypersurface $\Sigma$
and its normal on a discrete grid, which can be achieved by a simple
comparison: A cell $i$ contains part of the hypersurface if its temperature
$T_{i}$ fulfills the following inequality:

$\left(T_{i+1}-T_{\text{f.o.}}\right)\left(T_{\text{f.o.}}-T_{i}\right)\geq 0$
(24)

Here, $T_{\text{f.o.}}$ is the freeze-out value of the temperature, while
$T_{i+1}$ is the temperature in one of the neighboring cells. Once we have
found the hypersurface, we employ the Cornelius method [53] to calculate the
discretized normal vector $\Delta\sigma_{\mu}$ of a two- or three-dimensional
isosurface element in a three- or four-dimensional discrete grid. This then
allows to evaluate a discretized version of the Cooper-Frye formula (23),
namely

$E_{i}\,\dfrac{d^{3}N_{i}}{dp^{3}}\approx\dfrac{g_{i}}{\left(2\pi\hbar\right)^{3}}\sum_{\sigma}\Delta\sigma_{\mu}p^{\mu}f_{i}\\!\left(x,p\right)\;.$
(25)

The four-momenta of the outgoing particles can then be obtained by a rejection
sampling of Eq. (25). It is also possible to use the output provided by BHAC-
QGP for any particle sampler, such as frzout [54], to generate particles.

Figure 4: Total-energy density distribution $\tau=3.4~{}\text{fm}$
($\tau_{0}=0.4~{}\text{fm}$) resulting from a Au-Au collision at
$\sqrt{s_{\text{NN}}}=200~{}\text{GeV}$, and an impact parameter of
$b=8~{}\text{fm}$. The simulation was performed using BHAC-QGP with a freeze-
out energy density of $e_{\text{f.o.}}=0.3~{}\text{GeV}/\text{fm}^{3}$, which
is indicated by the black contour line. The blue arrows correspond to the
normal vectors of the hypersurface, which are evaluated with the Cornelius
method [53]. For visibility, only every fourth normal vector is plotted.

###  II.4 Anisotropic flow

Pressure gradients inside the fireball will result in a collective flow of
matter of these particles, resulting in a radially symmetric expansion if the
fireball is azimuthally symmetric in the transverse plane. However, in non-
central collisions, the fireball takes on an almond-like shape due to the
nonvanishing impact parameter, which will lead to stronger pressure gradients
in the $x$-direction (i.e., in the reaction plane) than in the $y$-direction
(i.e., perpendicular to the reaction plane, leading to an anisotropic flow of
matter, which is stronger in $x$\- than in $y$-direction. This anisotropic
flow is therefore an observable that reflects the geometry of the collision
zone, describing the correlation between the direction of the flow and the
distribution of matter. The anisotropic flow is reflected in the momentum
distribution of the particles emitted from the system. In order to quantify
this, one decomposes the azimuthal momentum distribution into a Fourier
series,

$\dfrac{d^{3}N}{p_{T}\,dp_{T}dYd\phi}=\dfrac{d^{2}N}{2\pi\,p_{T}\,dp_{T}dY}\left[1+2\sum\limits_{n=1}^{\infty}v_{n}\cos\\!\left(n\,\phi\right)\right]\;,$
(26)

where the Fourier coefficient $v_{n}$ of the $n$-th harmonic is defined as:

$v_{n}:=\langle\cos\\!\left(n\phi\right)\rangle\;,$ (27)

while $\langle...\rangle$ denotes an average over azimuthal angle , i.e.,

$v_{n}=\dfrac{\int\limits_{0}^{2\pi}\dfrac{d^{3}N}{p_{T}\,dp_{T}dYd\phi}\,\cos\\!\left(n\phi\right)d\phi}{\int\limits_{0}^{2\pi}\dfrac{d^{3}N}{p_{T}\,dp_{T}dYd\phi}d\phi}\;.$
(28)

The first of these harmonics, $v_{1}$, is referred to as directed flow. It
describes the deflection of particles away from the beam axis, primarily
caused by the bounce-off effect from the fireball matter. Hence, $v_{1}$
serves as a measure for quantifying the total amount of transverse flow. The
far more prominent observable is the second coefficient, $v_{2}$, or elliptic
flow, whose definition is

$v_{2}=\langle\cos\\!\left(2\phi\right)\rangle=\left\langle\dfrac{p_{x}^{2}-p_{y}^{2}}{p_{x}^{2}+p_{y}^{2}}\right\rangle\;,$
(29)

and which is most closely related to the initial anisotropy of the colliding
nuclei. As explained above, the initial energy-density distribution exhibits
azimuthal asymmetry in the transverse plane, leading to stronger pressure
gradients along the $x$-direction compared to the $y$-direction. Consequently,
the fluid will have a higher velocity in the $x$-direction than in
$y$-direction, resulting in a non-zero and positive value for $v_{2}$.
Therefore, the elliptic flow is a useful observable for gaining an
understanding of the initial anisotropy of the fireball matter.

##  III Results of Numerical Simulations for Heavy-Ion Collisions

In this section, we present some results of relativistic Au-Au collisions
calculated with BHAC-QGP. Since we are particularly interested in using BHAC-
QGP to study the effects of magnetic fields on the QGP, we first start with an
investigation of the lifetime of the magnetic field, comparing the numerical
results with those of ECHO-QGP [55]. Finally, we showcase particle spectra for
different Au-Au collisions at RHIC energies, specifically at
$\sqrt{s_{\text{NN}}}=200~{}\text{GeV}$. These results are not meant to
reproduce experimental results, but are intended to highlight the diverse
capabilities of BHAC-QGP.

If not stated otherwise, we adopt a constant inelastistic nucleus-nucleus
cross section of $\sigma_{\text{inel}}=4.2~{}\text{fm}^{2}$ and Milne
coordinates, which are defined as

$\left(\tau,x,y,\eta_{S}\right):=\left(\sqrt{t^{2}-z^{2}},x,y,\dfrac{1}{2}\text{ln}\left(\dfrac{t+z}{t-z}\right)\right)\;.$
(30)

We assume that the fireball approaches local thermodynamical equilibrium
rather fast and set as initial time $\tau_{0}=0.4~{}\text{fm}$. The system of
hydrodynamic equations is closed with the equation of state for a massless gas
of classical particles, $P=e/3$. For the (2+1)-dimensional simulations we
define a computational domain of $[-15~{}\mathrm{fm},+15~{}\mathrm{fm}]^{2}$,
with 150 cells in each direction, while the (3+1)-dimensional calculations are
performed on a grid with size $[-20~{}\mathrm{fm},+20~{}\mathrm{fm}]^{3}$ and
200 cells in each direction. For the initialization of the energy density we
set $e_{0}=55~{}\text{GeV}/\text{fm}^{3}$, while the longitudinal expansion is
computed with $\Delta_{\eta_{S}}=5.9$ and $\sigma_{\eta_{S}}=0.4$.

###  III.1 Decay of magnetic field

Since the lifetime of the magnetic field is a major uncertainty in heavy-ion
collisions, we investigate the temporal evolution of the magnetic field
strength squared $\boldsymbol{B}^{2}$ for different cases. We compare the
decay of the magnetic field strength at the center of the grid with that for
generalized Bjorken flow and various other field configurations. The reference
magnetic field corresponds to that of a peripheral Au-Au collision
$(b=10~{}\text{fm})$ at RHIC energies of
$\sqrt{s_{\text{NN}}}=200~{}\text{GeV}$, while the energy density is
initialized according to the optical Glauber model, see Eq. (II.1), with
$\alpha_{H}=0.15$ and $\sigma_{\text{inel}}=4.2~{}\text{fm}^{2}$. For the QGP,
we assume an electric conductivity of $\sigma_{E}=5.8~{}\mathrm{MeV}$ and a
chiral conductivity of $\sigma_{\chi}=1.5~{}\mathrm{MeV}$. Both the fluid and
the magnetic field contribute to the total pressure,
$P_{\text{tot}}=P+B^{2}/2$, causing the matter to expand. In Fig. 5, we
compare the evolution of this magnetic field (solid blue line in the figure)
with that for generalized Bjorken flow, corresponding to the evolution of a
uniform pressure and magnetic field. As Refs. [56, 57] have shown, in such a
scenario the magnetic field evolves according to $B/B_{0}=\tau_{0}/\tau$
(solid black line in the figure). In contrast, the magnetic field of a heavy-
ion collision is not uniform, but forms additional pressure gradients, which
are responsible for a faster decay. The stronger the magnetic field is, the
stronger the pressure gradients that are formed, which drive matter faster
away from the center (see dashed red line in the figure for which each
component of the magnetic field was enhanced by a factor of 50). As an
additional comparison, we have also explored the case where the fluid pressure
is uniform $(P=1~{}\mathrm{GeV}/~{}\mathrm{fm}^{3})$, while the magnetic field
mirrors that of a heavy-ion collision (green filled circles in Fig. 5). In
this case, the evolution of the magnetic field is more similar to the one of
generalized Bjorken flow, with small differences, resulting from the magnetic-
pressure gradients that are absent for generalized Bjorken flow. Also in this
case, increasing the magnetic field would cause a faster time decay. If, on
the other hand, the fluid pressure is set up according to the optical Glauber
model and a Gaussian distribution for $B_{x}$ and $B_{y}$ is chosen, the total
pressure gradients are stronger and the decay is correspondingly stronger.
However, there is no noteworthy difference between a symmetric
$(\sigma_{x}=3~{}\mathrm{fm},\sigma_{y}=3~{}\mathrm{fm})$ and an asymmetric
$(\sigma_{x}=3~{}\mathrm{fm},\sigma_{y}=9~{}\mathrm{fm})$ distribution of the
magnetic field components (solid orange line and green diamonds in the figure,
respectively).

Figure 5: Time evolution of the strength of the magnetic field at the center
of the grid $(x=y=0~{}\mathrm{fm})$. Various magnetic-field configurations are
compared with the magnetic field of a heavy-ion collision (HIC) and the one of
Bjorken’s solution. The magnetic field of the HIC in BHAC-QGP is initialized
with an electric conductivity of $\sigma_{E}=5.8~{}\text{MeV}$ and a chiral
conductivity of $\sigma_{\chi}=1.5~{}\text{MeV}$.

###  III.2 Influence of magnetic field on the QGP evolution

The magnetic field produced in a relativistic heavy-ion collision is
responsible for a pressure that, in particular outside the collision zone, is
much stronger than the fluid pressure. Indeed, the ratio between the magnetic
pressure $P_{\text{mag}}:=B^{2}/2$ and the fluid pressure $P$, the so-called
inverse plasma-$\beta$ parameter, is likely to exceed $100$ outside the
collision zone. Such regions are very difficult to handle for MHD codes, since
the conversion between conserved and primitive variables can quickly fail.
Using the equation of state $P=e/3$, in the ideal-MHD regime, BHAC-QGP can
exploit an additional evolution equation for the entropy, as described in
Paper I. Nevertheless, to keep the initial inverse plasma-$\beta$ parameter
controlled, it has proven to be advantageous to establish a relatively high
minimum value for the energy density of the fluid, denoted as
$e_{\text{min}}$. Unless otherwise specified, we use
$e_{\text{min}}=10^{-3}~{}\mathrm{GeV}/\mathrm{fm}^{3}$, which is a negligible
value compared to the freeze-out energy density, $e_{\text{f.o.}}\sim
0.5~{}\mathrm{GeV}/\mathrm{fm}^{3}$, and allows us to perform any
ultrarelativistic Au-Au collision without numerical problems.

Figure 6: Evolution of the fluid pressure (top row), the magnetic field
strength (second row), and the inverse plasma-$\beta$ parameter (third row) in
a Au-Au collision with $b=4~{}\text{fm}$. Further parameters are listed in
Tab. 2. Since this is a rather central collision, the fluid pressure $P$ is
almost isotropic and subsequently expands smoothly in all directions. Compared
to the fluid pressure, the magnetic field is only noticeably strong outside
the fireball. It has no noteworthy influence on the evolution of the fireball.
Albeit further away from the center some regions are highly magnetized, BHAC-
QGP can simulate these regions with the help of the entropy evolution
equation.

As a result of the strong magnetic pressure, the magnetic field affects the
expansion of the QGP even in the ideal-MHD regime. To understand this
influence, we compared and studied the evolution of a central collision
$(b=4~{}\text{fm})$ with that of a peripheral collision $(b=10~{}\text{fm})$.
We show snapshots of the fluid pressure, magnetic field strength, and inverse
plasma-$\beta$ parameter in Fig. 6 (central collision) and Fig. 7 (peripheral
collision), respectively. In both figures, the first column refers to the
initial pressure distribution of the heavy-ion collision at
$\tau_{0}=0.4~{}\text{fm}$. The second column shows the respective quantities
at $\tau=2.4~{}\text{fm}$, while the third column refers to
$\tau=4.4~{}\text{fm}$. Clearly, the evolution of the pressure in case of the
central collision is relatively uniform and smooth, as can be seen in the
first row of Fig. 6. The magnetic field has no significant influence on its
evolution, since the fluid pressure dominates over the magnetic pressure,
especially in the collision zone. The total pressure gradient is then
responsible for the relatively isotropic expansion of the magnetic field,
which can be seen in the middle row of Fig. 6. The expanding fireball also
drives the magnetic field relatively isotropically away from the center. The
broken symmetry under reflection at the $y-$axis is a consequence of the
chiral medium, since the magnetic field is dominant only far outside the
collision zone. In fact, the magnetic pressure outside the fireball can be
more than 100 times stronger than the fluid pressure, as can be seen from the
third row of Fig. 6. However, as the snapshots prove, BHAC-QGP has no problem
calculating these highly magnetized areas.

The dynamics is quite different for a peripheral collision, whose snapshots
are presented in Fig. 7, and where the evolution of the fluid pressure is
strongly influenced by the magnetic field. The reason for this are the two
regions of lower magnetic pressure at $x=\pm 8~{}\text{fm}$ and
$y=0~{}\text{fm}$. As a result, the magnetic field exerts a magnetic pressure
on these areas, while, at the same time, the fireball expands towards these
areas. This leads to an interaction from which shock waves emerge that will
influence the evolution of the whole system. However, the chosen magnetic
fields are not strong enough to noticeably affect the hypersurface of the
energy density, at least not at the RHIC energies studied in this work. Note
that peripheral collisions lead to large regions where
$\beta^{-1}=P_{\text{mag}}/P>10^{2}$, which BHAC-QGP can handle satisfactorily
well.

Figure 7: Evolution of the fluid pressure (top row), the magnetic field
strength (second row), and the inverse plasma-$\beta$ parameter (third row) in
a Au-Au collision with $b=10~{}\text{fm}$. Further parameters are listed in
Tab. 2. Since it is a peripheral collision, the initial pressure distribution
has an almond-like shape, which results in a non-isotropic expansion. The
magnetic field is particularly strong outside the collision zone, with a
region of weaker magnetic field in the vicinity of each of the original ion
centers. Subsequently, the magnetic field expands into these areas from all
sides, resulting in interactions between them. This interaction influences in
particular the expansion at later times, as one can see in the plots of the
third column $(\tau=4.4~{}\mathrm{fm})$. As can be seen in the third row, the
area outside the collision zone is highly magnetized, but BHAC-QGP is able to
handle these regions using the entropy evolution equation.

###  III.3 Comparisons of 2D and 3D simulations

In order to check the correctness of the BHAC-QGP simulations, we compare the
(2+1)-dimensional results with those of (3+1)-dimensional simulations. This is
crucial as we can use the (2+1)-dimensional results as a comparison and
benchmark for (3+1)-dimensional simulations. The parameters for this
consistency check are reported in Tab. 1. For a consistent comparison, the
energy-density profile for the (2+1)-dimensional simulation must be the same
as for the (3+1)-dimensional calculation. Hence, for the (2+1)-dimensional
initialization we have to use Eq. (II.1) with $\eta_{S}=0$.

Parameter | Value
---|---
Collision hardness $\alpha_{H}$ | $0.05$
Inelastic cross section $\sigma_{\text{inel}}$ | $4.0~{}\text{fm}^{2}$
Electric conductivity $\sigma_{E}$ | $5.8~{}\text{MeV}$
Chiral conductivity $\sigma_{\chi}$ | $1.5~{}\text{MeV}$
Minimum of initial energy density $e_{\text{min}}$ | $10^{-4}~{}\text{GeV}/\text{fm}^{3}$
Table 1: Parameters used for the (3+1)-dimensional Au-Au simulation to study
the influence of the magnetic field on the evolution of the fireball.

As a particularly important measure of the dynamics and a useful quantity for
the comparison we investigated the evolution of the elliptic flow. Figure 8
displays the elliptic flow $v_{2}$ of all pions as function of the transverse
momentum $p_{T}=(p_{x}^{2}+p_{y}^{2})^{1/2}$ and shows that BHAC-QGP is able
to yield results for (2+1)-dimensional and (3+1)-dimensional simulations that
are clearly consistent.

Figure 8: Elliptic flow of pions as function of the transverse momentum
obtained from a (2+1)-dimensional calculation compared with the one obtained
from a (3+1)-dimensional calculation. The spectra pertain to a Au-Au collision
with parameters listed in Tab. 1. As the figure proves, BHAC-QGP is able to
yield consistent results for rather central $(b=4~{}\mathrm{fm})$ and
peripheral $(b=8~{}\mathrm{fm})$ collisions.

Another validation of BHAC-QGP is offered by a comparison with the ECHO-QGP
code [55]. We recall that BHAC-QGP solves the GRMHD equations using finite-
volume methods, while ECHO-QGP uses finite-difference methods. In addition,
ECHO-QGP uses a different method for the evaluation of the Cooper-Frye
formula. Therefore, the comparison with ECHO-QGP allows us to further validate
the correct implementation of the freeze-out procedure. We compare the results
for a Au-Au collision whose parameters are again those listed in Tab. 1.
Figure 9 displays the transverse-momentum distribution of positively charged
pions, while Fig. 10 shows the corresponding elliptic flow for a Au-Au
collision at $b=10~{}\text{fm}$. As can be seen from the plots, both codes
yield results that are in good agreement. One observes that the transverse-
momentum distributions have the same slope, but are offset by a factor $\sim
2$. This is most likely due to the different numerical evaluation of the
Cooper-Frye formula employed in BHAC-QGP (which uses the Cornelius routine)
vs. ECHO-QGP (which uses a different method, see Ref. [38]).

Figure 9: Comparison between the transverse-momentum distribution of
positively charged pions for a simulation without an external magnetic field
and a simulation with an external magnetic field. One notices that the
magnetic field does not provide sufficient magnetic energy for an enhanced
pion production.  Figure 10: Comparison between the elliptic flow of
positively charged pions for a simulation without an external magnetic field
and a simulation with an external magnetic field. The magnetic field has no
visible influence on the elliptic flow because it does not significantly
change the shape of the initial fireball. Furthermore plotted is a direct
comparison with results obtained from ECHO-QGP. Overall, BHAC-QGP and ECHO-QGP
show a very good agreement on the elliptic flow.

###  III.4 Particle spectra

Having presented the consistency checks between (2+1)-dimensional and
(3+1)-dimensional simulations in the previous section we now focus entirely on
(3+1)-dimensional simulations performed in Milne coordinates. In particular,
we present various particle spectra obtained from typical Au-Au collisions at
RHIC energies of $\sqrt{s_{\text{NN}}}=200~{}\text{GeV}$, where all parameters
used are summarized in Tab. 2. The impact parameter $b$ and the freeze-out
temperature $T_{\mathrm{f.o.}}$ will be varied and quoted separately.

Parameter | Value
---|---
Initial time $\tau_{0}$ | $0.4~{}\text{fm}$
Center-of-mass energy $\sqrt{s_{\text{NN}}}$ | $200~{}\text{GeV}$
Collision hardness $\alpha_{H}$ | $0.15$
Inelastic cross section $\sigma_{\text{inel}}$ | $4.2~{}\text{fm}^{2}$
Electric conductivity $\sigma_{E}$ | $5.8~{}\text{MeV}$
Chiral conductivity $\sigma_{\chi}$ | $1.5~{}\text{MeV}$
Width $\Delta_{\eta_{S}}$ | $5.9~{}\text{MeV}$
Gaussian width $\sigma_{\eta_{S}}$ | $0.4~{}\text{MeV}$
Minimum of initial energy density $e_{\text{min}}$ | $10^{-3}~{}\text{GeV}/\text{fm}^{3}$
Equation of state | $P=e/3$
Table 2: Parameters used for the (3+1)-dimensional Au-Au collisions in the
calculation of particle spectra.

First, as simple check, we study collisions for different values of the impact
parameter $b$ and a freeze-out temperature of
$T_{\text{f.o.}}=130~{}\text{MeV}$. Figure 11 shows the total number of pions
as a function of the transverse momentum for different impact parameters.
According to the optical Glauber model, a small impact parameter corresponds
to more participants, so that more energy is deposited in the collision zone.
Accordingly, in collisions with small impact parameters, more pions can be
produced, which is clearly visible in Fig. 11. Furthermore, a peripheral
collision leads to a pronounced anisotropy, and, indeed, in these cases the
fireball shows an almond-shaped form so that anisotropic pressure gradients
are formed which are much stronger in the $x$-direction than in the
$y$-direction. The consequence is a strongly pronounced positive elliptic flow
$v_{2}$. In contrast, the more central the collision, the more uniform is the
fireball and its expansion, so that the corresponding pressure anisotropies
and subsequently the elliptic flow are weaker. All of this is seen in Fig. 12.

Figure 11: Comparison between the transverse-momentum distribution of pions
($\pi^{\pm},\pi^{0}$) for different impact parameters. The more central the
collision, the larger the fireball and the higher its energy density. A higher
energy density allows for the generation of more particles. The spectra are
obtained from Au-Au collisions with parameters listed in Tab. 2. Figure 12:
Elliptic flow $v_{2}$ of pions ($\pi^{\pm},\pi^{0}$) as function of the
transverse momentum $p_{T}$ for different impact parameters. The more
peripheral the collision, the more pronounced is the initial anisotropy of the
fireball. Since the elliptic flow is a measure of this eccentricity, larger
values of the impact parameter will lead to a more pronounced $v_{2}$. The
spectra are obtained from Au-Au collisions with parameters listed in Tab. 2.

The effect of the initial anisotropy can also be investigated when studying
pions as a function of the particle’s rapidity
$Y:=1/2~{}\mathrm{ln}\\!\left((E+p_{z})/(E-p_{z})\right)$. As can be deduced
from Fig. 13, the elliptic flow $v_{2}$ is positive as a consequence of the
stronger pressure gradients in the $x$-direction. While the two ions collide
and the spectators leave the collision zone in $z$-direction, pions are
produced in the hot and dense fireball. Figure 14 emphasizes again that much
more energy is available for particle production in central collisions.

Figure 13: Elliptic flow of pions ($\pi^{\pm},\pi^{0}$) as a function of the
rapidity $Y$ for different impact parameters. Larger impact parameters lead to
stronger pressure anisotropies, resulting in a much more pronounced elliptic
flow. The freeze-out temperature was chosen as
$T_{\text{f.o.}}=130~{}\text{MeV}$. Figure 14: Number of produced pions
($\pi^{\pm},\pi^{0}$) as a function of the rapidity $Y$. The freeze-out
temperature was chosen as $T_{\text{f.o.}}=130~{}\text{MeV}$.

In addition to the impact parameter, BHAC-QGP also allows us to study any
freeze-out temperature and the number of produced pions as a function of the
rapidity for different freeze-out temperatures. The precise relationship
between the number of produced pions and the freeze-out temperature is
determined by the dynamics of the collision and the properties of the
interacting particles. In general, a lower freeze-out temperature results in a
lower number of pions. For this reason, the freeze-out temperature has a
strong influence on the elliptic flow, as can be seen in Fig. 15. Furthermore,
at lower freeze-out temperatures, the system has a longer evolution time,
allowing for stronger interactions and collective behavior among the
particles. This leads to an enhanced elliptic flow of pions as can be seen in
Fig. 16.

Figure 15: Number of produced pions ($\pi^{\pm},\pi^{0}$) as a function of the
rapidity $Y$ for different freeze-out temperatures. The Au-Au collisions were
computed with an impact parameter of $b=8~{}\text{fm}$.  Figure 16: Elliptic
flow of pions ($\pi^{\pm},\pi^{0}$) as a function of the rapidity $Y$ for
different freeze-out temperatures. The Au-Au collisions were computed with an
impact parameter of $b=8~{}\text{fm}$.

Investigating the elliptic flow is particularly useful as it is related to the
initial anisotropy, as well as it is sensitive to the equation of state, which
governs the evolution of the fluid. In order to find out if the strong
external magnetic field also has an influence on the momentum anisotropy, we
calculated $v_{2}$ for heavy-ion collisions with and without external magnetic
fields. Although the magnetic fields are strong, especially outside the
fireball, they seem to have no significant effect for the energy normally
relevant for heavy-ion collisions. This was already suggested by Fig. 10 and
confirmed by the results obtained via ECHO-QGP [55]. Nevertheless, the
magnetic field influences the expansion of the fireball, especially for
peripheral collisions as discussed with Fig. 7. In particular, a sufficiently
strong magnetic field contributes a significant portion to the total energy
density of the fireball and as a result of this additional energy density,
many more particles can be produced. This can be seen in Fig. 17, where we
have increased each component of the magnetic field by a factor of either two,
three, or four. In this figure, lines of different color refer to different
amplifications compared to the standard value of the magnetic field. However,
the magnetic field is not strong enough to add sufficient magnetic energy
density, which scales quadratically with the magnetic field, for the
production of pions. On the other hand, a higher magnetic field strength also
leads to stronger magnetic pressures, and total pressure gradients, resulting
in a stronger expansion, and, ultimately to particles with higher transverse
momentum. Since these pressure gradients are also not isotropic, they further
amplify the anisotropy. This ultimately leads to an enhanced elliptic flow, as
highlighted by Fig. 18.

Figure 17: Effect of the magnetic field on the total number of produced pions.
The red line refers to a simulation where each component of the magnetic field
was enhanced by a factor of 2, the green line to a simulation where each
component of the magnetic field was enhanced by a factor of 3, while in case
of the yellow line each component was amplified by a factor of 4. In all
cases, the magnetic field is does not contribute sufficient energy to enhance
the production of pions. Figure 18: Effect of the magnetic field on the
elliptic flow of pions. The red line refers to a simulation where each
component of the magnetic field was enhanced by a factor of 2, the green line
to a simulation where each component of the magnetic field was enhanced by a
factor of 3, while in case of the yellow line each component was amplified by
a factor of 4. Stronger magnetic fields result in stronger pressure gradients
and a subsequent amplification of the elliptic flow.

##  IV Conclusions and Outlook

We presented results obtained from a novel (3+1)-dimensional relativistic
magnetohydrodynamics code, BHAC-QGP. BHAC-QGP is an extension of the Black
Hole Accretion Code (BHAC), which is capable of evolving the matter produced
in relativistic heavy-ion collisions. In order to perform such an evolution,
we have implemented Milne coordinates and the ultrarelativistic equation of
state $(P=e/3)$ for a generic investigation of heavy-ion collisions. We plan
to implement more equations of state in the future. Moreover, BHAC-QGP can
resort to the entropy evolution equation in the ideal-MHD regime, which allows
to simulate also the strong magnetic fields outside the collision zone. This
allows us to study novel chiral effects, such as the Chiral Magnetic Effect
(CME), especially because the novel code is capable of simulating a QGP with a
non-zero electric and a non-zero chiral conductivity.

BHAC-QGP uses the optical Glauber model to initialize the energy density that
is produced during the collision of two heavy ions. Therefore, all types of
heavy-ion collisions with arbitrary impact parameters can be calculated.
Additionally, our ongoing work involves implementing the option to use initial
energy-density data extracted from transport codes like SMASH [58]. This will
allow us to simulate a more randomly distributed and realistic energy-density
profile, which is particularly crucial for the investigation of higher flow-
harmonic coefficients, such as the triangular flow $v_{3}$.

BHAC-QGP relies on the Cornelius routine to evaluate the Cooper-Frye freeze-
out formula. Subsequently, the resulting data can be used to generate particle
spectra. With this implementation, BHAC-QGP is able to reproduce spectra from
ECHO-QGP, although ECHO-QGP uses different routines.

Another interesting topic is the investigation of a rotating QGP. A QGP with a
non-zero angular momentum can give rise to the chiral vortical effect (CVE).
However, its investigation requires very high resolutions. BHAC-QGP, with its
AMR framework, is excellently positioned for this.

## Acknowledgments

The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the CRC-TR 211 “Strong-interaction matter
under extreme conditions” – project number 315477589 – TRR 211. The work is
supported by the State of Hesse within the Research Cluster ELEMENTS (Project
ID 500/10.006). Computational resources have been provided by the Center for
Scientific Computing (CSC) at the Goethe University. M.M. would like to thank
H. Olivares, M. Chabanov, N. Kuebler, and J. Sammet for fruitful discussions.

## References

  * Skokov  _et al._ [2009] V. Skokov, A. Y. Illarionov, and V. Toneev, Estimate of the magnetic field strength in heavy-ion collisions, Int. J. Mod. Phys. A 24, 5925 (2009), arXiv:0907.1396 [nucl-th] .
  * Tuchin [2014] K. Tuchin, Electromagnetic fields in high energy heavy-ion collisions, Int. J. Mod. Phys. E 23, 1430001 (2014).
  * Tuchin [2013a] K. Tuchin, Time and space dependence of the electromagnetic field in relativistic heavy-ion collisions, Phys. Rev. C 88, 024911 (2013a), arXiv:1305.5806 [hep-ph] .
  * Tuchin [2013b] K. Tuchin, Particle production in strong electromagnetic fields in relativistic heavy-ion collisions, Adv. High Energy Phys. 2013, 490495 (2013b), arXiv:1301.0099 [hep-ph] .
  * Bzdak and Skokov [2012] A. Bzdak and V. Skokov, Event-by-event fluctuations of magnetic and electric fields in heavy ion collisions, Phys. Lett. B 710, 171 (2012), arXiv:1111.1949 [hep-ph] .
  * Deng and Huang [2012] W.-T. Deng and X.-G. Huang, Event-by-event generation of electromagnetic fields in heavy-ion collisions, Phys. Rev. C 85, 044907 (2012), arXiv:1201.5108 [nucl-th] .
  * McLerran and Skokov [2014] L. McLerran and V. Skokov, Comments About the Electromagnetic Field in Heavy-Ion Collisions, Nucl. Phys. A 929, 184 (2014), arXiv:1305.0774 [hep-ph] .
  * Voronyuk  _et al._ [2011] V. Voronyuk, V. D. Toneev, W. Cassing, E. L. Bratkovskaya, V. P. Konchakovski, and S. A. Voloshin, (Electro-)Magnetic field evolution in relativistic heavy-ion collisions, Phys. Rev. C 83, 054911 (2011), arXiv:1103.4239 [nucl-th] .
  * Li  _et al._ [2016] H. Li, X.-l. Sheng, and Q. Wang, Electromagnetic fields with electric and chiral magnetic conductivities in heavy ion collisions, Phys. Rev. C 94, 044903 (2016), arXiv:1602.02223 [nucl-th] .
  * Gursoy  _et al._ [2014] U. Gursoy, D. Kharzeev, and K. Rajagopal, Magnetohydrodynamics, charged currents and directed flow in heavy ion collisions, Phys. Rev. C 89, 054905 (2014), arXiv:1401.3805 [hep-ph] .
  * Pasechnik and Šumbera [2017] R. Pasechnik and M. Šumbera, Phenomenological Review on Quark–Gluon Plasma: Concepts vs. Observations, Universe 3, 7 (2017), arXiv:1611.01533 [hep-ph] .
  * Huang [2016] X.-G. Huang, Electromagnetic fields and anomalous transports in heavy-ion collisions — A pedagogical review, Rept. Prog. Phys. 79, 076302 (2016), arXiv:1509.04073 [nucl-th] .
  * Roy  _et al._ [2017] V. Roy, S. Pu, L. Rezzolla, and D. H. Rischke, Effect of intense magnetic fields on reduced-MHD evolution in $\sqrt{s_{\rm NN}}$ = 200 GeV Au+Au collisions, Phys. Rev. C 96, 054909 (2017), arXiv:1706.05326 [nucl-th] .
  * Inghirami  _et al._ [2020] G. Inghirami, M. Mace, Y. Hirono, L. Del Zanna, D. E. Kharzeev, and M. Bleicher, Magnetic fields in heavy ion collisions: flow and charge transport, Eur. Phys. J. C 80, 293 (2020), arXiv:1908.07605 [hep-ph] .
  * Borsanyi  _et al._ [2010] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo (Wuppertal-Budapest), Is there still any $T_{c}$ mystery in lattice QCD? Results with physical masses in the continuum limit III, JHEP 09, 073, arXiv:1005.3508 [hep-lat] .
  * Bazavov  _et al._ [2012] A. Bazavov _et al._ , The chiral and deconfinement aspects of the QCD transition, Phys. Rev. D 85, 054503 (2012), arXiv:1111.1710 [hep-lat] .
  * Kharzeev  _et al._ [1998] D. Kharzeev, R. D. Pisarski, and M. H. G. Tytgat, Possibility of spontaneous parity violation in hot QCD, Phys. Rev. Lett. 81, 512 (1998), arXiv:hep-ph/9804221 .
  * Kharzeev  _et al._ [2008] D. E. Kharzeev, L. D. McLerran, and H. J. Warringa, The Effects of topological charge change in heavy ion collisions: ’Event by event P and CP violation’, Nucl. Phys. A 803, 227 (2008), arXiv:0711.0950 [hep-ph] .
  * Pu  _et al._ [2015] S. Pu, S.-Y. Wu, and D.-L. Yang, Chiral Hall Effect and Chiral Electric Waves, Phys. Rev. D 91, 025011 (2015), arXiv:1407.3168 [hep-th] .
  * Kharzeev  _et al._ [2016] D. E. Kharzeev, J. Liao, S. A. Voloshin, and G. Wang, Chiral magnetic and vortical effects in high-energy nuclear collisions—A status report, Prog. Part. Nucl. Phys. 88, 1 (2016), arXiv:1511.04050 [hep-ph] .
  * Fukushima  _et al._ [2008] K. Fukushima, D. E. Kharzeev, and H. J. Warringa, The Chiral Magnetic Effect, Phys. Rev. D 78, 074033 (2008), arXiv:0808.3382 [hep-ph] .
  * Sadofyev and Isachenkov [2011] A. V. Sadofyev and M. V. Isachenkov, The Chiral magnetic effect in hydrodynamical approach, Phys. Lett. B 697, 404 (2011), arXiv:1010.1550 [hep-th] .
  * Kharzeev and Yee [2011] D. E. Kharzeev and H.-U. Yee, Chiral Magnetic Wave, Phys. Rev. D 83, 085007 (2011), arXiv:1012.6026 [hep-th] .
  * Kharzeev and Son [2011] D. E. Kharzeev and D. T. Son, Testing the chiral magnetic and chiral vortical effects in heavy ion collisions, Phys. Rev. Lett. 106, 062301 (2011), arXiv:1010.0038 [hep-ph] .
  * Jiang  _et al._ [2015] Y. Jiang, X.-G. Huang, and J. Liao, Chiral vortical wave and induced flavor charge transport in a rotating quark-gluon plasma, Phys. Rev. D 92, 071501 (2015), arXiv:1504.03201 [hep-ph] .
  * Son and Surowka [2009] D. T. Son and P. Surowka, Hydrodynamics with Triangle Anomalies, Phys. Rev. Lett. 103, 191601 (2009), arXiv:0906.5044 [hep-th] .
  * Asakawa  _et al._ [2010] M. Asakawa, A. Majumder, and B. Muller, Electric Charge Separation in Strong Transient Magnetic Fields, Phys. Rev. C 81, 064912 (2010), arXiv:1003.2436 [hep-ph] .
  * Stephanov and Yin [2012] M. A. Stephanov and Y. Yin, Chiral Kinetic Theory, Phys. Rev. Lett. 109, 162001 (2012), arXiv:1207.0747 [hep-th] .
  * Dumitru and Rischke [1999] A. Dumitru and D. H. Rischke, Collective dynamics in highly relativistic heavy ion collisions, Phys. Rev. C 59, 354 (1999), arXiv:nucl-th/9806003 .
  * Molnar and Huovinen [2005] D. Molnar and P. Huovinen, Dissipation and elliptic flow at RHIC, Phys. Rev. Lett. 94, 012302 (2005), arXiv:nucl-th/0404065 .
  * Petersen  _et al._ [2015] H. Petersen, J. Steinheimer, J. Auvinen, and M. Bleicher, The beam energy dependence of collective flow in heavy ion collisions, PoS CPOD2014, 024 (2015), arXiv:1503.03363 [nucl-th] .
  * Petersen [2014] H. Petersen, Anisotropic flow in transport + hydrodynamics hybrid approaches, J. Phys. G 41, 124005 (2014), arXiv:1404.1763 [nucl-th] .
  * Heinz and Snellings [2013] U. Heinz and R. Snellings, Collective flow and viscosity in relativistic heavy-ion collisions, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th] .
  * Vogt [2007] R. Vogt,  _Ultrarelativistic heavy-ion collisions_ (Elsevier, Amsterdam, 2007).
  * Jaiswal and Roy [2016] A. Jaiswal and V. Roy, Relativistic hydrodynamics in heavy-ion collisions: general aspects and recent developments, Adv. High Energy Phys. 2016, 9623034 (2016), arXiv:1605.08694 [nucl-th] .
  * d’Enterria and Loizides [2021] D. d’Enterria and C. Loizides, Progress in the Glauber Model at Collider Energies, Ann. Rev. Nucl. Part. Sci. 71, 315 (2021), arXiv:2011.14909 [hep-ph] .
  * Kolb  _et al._ [2001] P. F. Kolb, U. W. Heinz, P. Huovinen, K. J. Eskola, and K. Tuominen, Centrality dependence of multiplicity, transverse energy, and elliptic flow from hydrodynamics, Nucl. Phys. A 696, 197 (2001), arXiv:hep-ph/0103234 .
  * Del Zanna  _et al._ [2013] L. Del Zanna, V. Chandra, G. Inghirami, V. Rolando, A. Beraudo, A. De Pace, G. Pagliara, A. Drago, and F. Becattini, Relativistic viscous hydrodynamics for heavy-ion collisions with ECHO-QGP, Eur. Phys. J. C 73, 2524 (2013), arXiv:1305.7052 [nucl-th] .
  * Karpenko  _et al._ [2014] I. Karpenko, P. Huovinen, and M. Bleicher, A 3+1 dimensional viscous hydrodynamic code for relativistic heavy ion collisions, Comput. Phys. Commun. 185, 3016 (2014), arXiv:1312.4160 [nucl-th] .
  * Schenke  _et al._ [2010] B. Schenke, S. Jeon, and C. Gale, (3+1)D hydrodynamic simulation of relativistic heavy-ion collisions, Phys. Rev. C 82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Miller  _et al._ [2007] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Glauber modeling in high energy nuclear collisions, Ann. Rev. Nucl. Part. Sci. 57, 205 (2007), arXiv:nucl-ex/0701025 .
  * Florkowski [2010] W. Florkowski,  _Phenomenology of Ultra-Relativistic Heavy-Ion Collisions_ (2010).
  * Bohr and Mottelson [1998] A. N. Bohr and B. R. Mottelson,  _Nuclear Structure (in 2 volumes)_ (World Scientific Publishing Company, 1998).
  * Kharzeev and Nardi [2001] D. Kharzeev and M. Nardi, Hadron production in nuclear collisions at RHIC and high density QCD, Phys. Lett. B 507, 121 (2001), arXiv:nucl-th/0012025 .
  * Hirano [2002] T. Hirano, Is early thermalization achieved only near mid-rapidity at RHIC?, Phys. Rev. C 65, 011901 (2002), arXiv:nucl-th/0108004 .
  * Hirano and Tsuda [2002] T. Hirano and K. Tsuda, Collective flow and two pion correlations from a relativistic hydrodynamic model with early chemical freezeout, Phys. Rev. C 66, 054905 (2002), arXiv:nucl-th/0205043 .
  * Nonaka and Bass [2007] C. Nonaka and S. A. Bass, Space-time evolution of bulk QCD matter, Phys. Rev. C 75, 014902 (2007), arXiv:nucl-th/0607018 .
  * Oliva [2020] L. Oliva, Electromagnetic fields and directed flow in large and small colliding systems at ultrarelativistic energies, Eur. Phys. J. A 56, 255 (2020), arXiv:2007.00560 [nucl-th] .
  * Ding  _et al._ [2011] H. T. Ding, A. Francis, O. Kaczmarek, F. Karsch, E. Laermann, and W. Soeldner, Thermal dilepton rate and electrical conductivity: An analysis of vector current correlation functions in quenched lattice QCD, Phys. Rev. D 83, 034504 (2011), arXiv:1012.4963 [hep-lat] .
  * Amato  _et al._ [2013] A. Amato, G. Aarts, C. Allton, P. Giudice, S. Hands, and J.-I. Skullerud, Electrical conductivity of the quark-gluon plasma across the deconfinement transition, Phys. Rev. Lett. 111, 172001 (2013), arXiv:1307.6763 [hep-lat] .
  * Aarts  _et al._ [2015] G. Aarts, C. Allton, A. Amato, P. Giudice, S. Hands, and J.-I. Skullerud, Electrical conductivity and charge diffusion in thermal QCD from the lattice, JHEP 02, 186, arXiv:1412.6411 [hep-lat] .
  * Cooper and Frye [1974] F. Cooper and G. Frye, Comment on the Single Particle Distribution in the Hydrodynamic and Statistical Thermodynamic Models of Multiparticle Production, Phys. Rev. D 10, 186 (1974).
  * Huovinen and Petersen [2012] P. Huovinen and H. Petersen, Particlization in hybrid models, Eur. Phys. J. A 48, 171 (2012), arXiv:1206.3371 [nucl-th] .
  * [54] http://qcd.phy.duke.edu/frzout/index.html.
  * Inghirami  _et al._ [2016] G. Inghirami, L. Del Zanna, A. Beraudo, M. H. Moghaddam, F. Becattini, and M. Bleicher, Numerical magneto-hydrodynamics for relativistic nuclear collisions, Eur. Phys. J. C 76, 659 (2016), arXiv:1609.03042 [hep-ph] .
  * Pu  _et al._ [2016] S. Pu, V. Roy, L. Rezzolla, and D. H. Rischke, Bjorken flow in one-dimensional relativistic magnetohydrodynamics with magnetization, Phys. Rev. D 93, 074022 (2016), arXiv:1602.04953 [nucl-th] .
  * Roy  _et al._ [2015] V. Roy, S. Pu, L. Rezzolla, and D. Rischke, Analytic Bjorken flow in one-dimensional relativistic magnetohydrodynamics, Phys. Lett. B 750, 45 (2015), arXiv:1506.06620 [nucl-th] .
  * Weil  _et al._ [2016] J. Weil _et al._ , Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions, Phys. Rev. C 94, 054905 (2016), arXiv:1606.06642 [nucl-th] .