# Collective flow at SIS energies within a hadronic transport approach:
Influence of light nuclei formation and equation of state

J.-Mohs1,2, M. Ege2, H. Elfner3,2,1 and M. Mayer2 1Frankfurt Institute for
Advanced Studies, Ruth-Moufang-Strasse 1, 60438 Frankfurt am Main, Germany
2Institute for Theoretical Physics, Goethe University, Max-von-Laue-Strasse 1,
60438 Frankfurt am Main, Germany  3GSI Helmholtzzentrum für
Schwerionenforschung, Planckstr. 1, 64291 Darmstadt, Germany

###### Abstract

Collective flow observables are known to be a sensitive tool to gain insights
on the equation of state of nuclear matter from heavy-ion collision
observations. Towards more quantitative constraints one has to carefully
assess other influences on the collective behaviour. In this work a hadronic
transport approach SMASH (Simulating Many Accelerated Strongly-interacting
Hadrons) is applied to study the first four anisotropic flow coefficients in
Au+Au collisions at $E_{\rm lab}=1.23A$ GeV in the context of the recently
measured data by the HADES collaboration. In particular, the formation of
light nuclei is important in this energy regime. Two different approaches are
contrasted to each other: A clustering algorithm inspired by coalescence as
well as microscopic formation of deuterons via explicit cross-sections. The
sensitivity of directed and elliptic flow observables to the strength of the
Skyrme mean field is explored. In addition, it is demonstrated that the
rapidity-odd $v_{3}$ coefficient is practically zero in this energy regime and
the ratio of $v_{4}/v_{2}^{2}$ is close to the value of 0.5 expected from
hydrodynamic behaviour. This study establishes the current understanding of
collective behaviour within the SMASH approach and lays the ground for future
more quantitative constraints on the equation of state of nuclear matter
within improved mean field calculations.

Collective flow observables, equation of state, light nuclei formation, Monte-
Carlo simulations, Heavy-ion collisions

##  I Introduction

In heavy-ion collisions the anisotropy of the particle production in the
transverse plane is described in terms of flow coefficients. In transport
calculations flow observables were shown to be very sensitive to nuclear
potentials and the equation of state (EoS) Scheid  _et al._ (1974);
Danielewicz  _et al._ (2002). Hence one can put constraints on the EoS by
comparing hadronic transport calculations with experimental flow data.

Dependent on the active degrees of freedom in the different transport
approaches different conclusions about the strength of the mean-field
interactions are drawn. Approaches mainly based on nucleons and pions with
sophisticated potentials including momentum dependent interactions as
described in Gale  _et al._ (1987) come to the best agreement with the wealth
of existing FOPI data Reisdorf  _et al._ (2012), when mean fields
corresponding to a ’soft’ equation of state are employed Aichelin  _et al._
(1987); Isse  _et al._ (2005); Le Fèvre  _et al._ (2016, 2018). On the other
hand, approaches like UrQMD with many resonance states arrive at the same
quality of agreement with experimental data incorporating a mean field
corresponding to a ’hard’ equation of state Petersen  _et al._ (2006);
Hillmann  _et al._ (2018, 2020). Modifying the equation of state via
modifications of the collision term offers an orthogonal avenue to study
collective behaviour within transport approaches Nara  _et al._ (2017).

The nuclear equation of state at high densities is relevant also for the
dynamical description of neutron star mergers. After the first detection of
gravitational wave signals Abbott  _et al._ (2016), there is a large interest
in sophisticated theoretical calculations of their detailed evolution which is
relevant for the nucleosynthesis including different assumptions about the
equation of state of QCD matter at high densities (see Bauswein  _et al._
(2019); Most  _et al._ (2018) for recent examples). Despite the isospin
difference that is encoded in the symmetry energy that can also be extracted
from heavy-ion reactions Wang  _et al._ (2020); Tsang  _et al._ (2012),
further knowledge on the equation of state is important. Before a quantitative
extraction including uncertainty quantification using Bayesian techniques is
sensible, all possible systematic uncertainties in the transport approaches
have to be understood. One of those is the formation of light clusters that
incorporates about half of the nucleons in heavy-ion reactions at few GeV per
nucleon beam energy. Coalescence models based on phase-space separation Butler
and Pearson (1961); Schwarzschild and Zupancic (1963); Gutbrod  _et al._
(1976) have been devised Sombun  _et al._ (2019) as well as more comprehensive
dynamical approaches Oh  _et al._ (2009); Aichelin  _et al._ (2020). Using
coalescence, protons and deuterons are predicted to follow a constituent
number scaling Hillmann  _et al._ (2020).

In this work we assess how SMASH-1.8 (Simulating Many Accelerated Strongly-
interacting Hadrons) Weil  _et al._ (2016); Oliinychenko  _et al._ (2020a)
performs compared to recent measurements of flow coefficients at SIS-18
enegies Kardan (2019); Adamczewski-Musch  _et al._ (2020) with a rather simple
Skyrme parametrization of the EoS Xu  _et al._ (2016) and investigate the
sensitivity of the results to the stiffness of the EoS. This is meant as a
baseline for further calculations employing more sophisticated mean field
implementations such as the one in Sorensen and Koch (2020). The main emphasis
in this work is on exploring different treatments for light nuclei production,
which constitutes one of the major uncertainties. We concentrate on the most
abundant lightest one namely the deuteron and compare two different ways of
taking the light nuclei formation into account. Coalescence of nucleons in the
final state as described in Sombun  _et al._ (2019) is contrasted to the
dynamical treatment by producing and propagating deuterons throughout the
evolution of a heavy-ion collision as described in Oliinychenko  _et al._
(2019).

The manuscript is structured as follows: First the SMASH transport with the
current state of nuclear potentials present in the calculation and the
different ways of treating light nuclei formation is described in Section II.
Then the effect of different Skyrme parameters and light nuclei treatments on
the directed flow in Section III and elliptic flow in Section IV of nucleons
and deuterons are evaluated. Afterwards the evolution of the directed and
elliptic flow with time is studied in Section V. Last we present some higher
flow coefficients and show the $v_{3}$ computed with the scalar product method
and the ratio of $v_{4}/v_{2}^{2}$ in Section VI. Finally we conclude and
present an outlook in Section VII.

##  II SMASH with potentials and light nuclei production

###  II.1 Hadronic transport with mean field

In this work, we apply the transport approach SMASH version 1.8 Weil  _et al._
(2016); Oliinychenko  _et al._ (2020a) to calculate flow observables. The
relativistic Boltzmann equation is effectively solved in a testparticle ansatz
with hadronic degrees of freedom. Starting with nulcei sampled from a Woods-
Saxon distribution hadrons are explicitly propagated between interactions. At
low collision energies most inelastic collisions form resonances with
properties adopted from Tanabashi  _et al._ (2018). For the resonances the
spectral functions are described by relativistic Breit-Wigner distributions
without medium modification. The mass dependence of the decay width is
included according to Manley and Saleski (1992).

For the propagation between the interactions nuclear potentials are optionally
taken into account. Let us note here, that for the current work we stick with
a simplified treatment, since this is the same setup that has been used in
other recent works in comparable approaches Hillmann  _et al._ (2018) and the
default parameters have been determined in a comparison between multiple
different transport codes Xu  _et al._ (2016).

At this point the Skyrme and symmetry potential are present in the calculation

$U=U_{\rm Sk}+U_{\rm Sym}.$ (1)

The contributions to the potential energy is given in terms of the density as

$U_{\rm
Sk}=A\left(\frac{\rho_{B}}{\rho_{0}}\right)+B\left(\frac{\rho_{B}}{\rho_{0}}\right)^{\tau},$
(2)

where $\rho_{B}$ is the net-baryon density, $\rho_{0}=0.168\,\rm fm^{-3}$ is
the nuclear ground state density and $A$, $B$ and $\tau$ are parameters.
Furthermore

$U_{\rm Sym}=\pm 2S_{\rm pot}\frac{\rho_{I_{3}}}{\rho_{0}},$ (3)

where $\rho_{I_{3}}$ is the density of the relative isospin projection
$I_{3}/I$, the sign is positive for positive isospin and negative for negative
isospin of the particle of interest.

Since the potentials are not written in a covariant form, they have to be
evaluated in the local rest frame for a Lorentz invariant treatment. In
practice, the calculation is performed in an arbitrary calculation frame such
as the fixed target or center of mass frame. Therefore a boost needs to be
introduced which results in the following form for the force acting on a
particle in the calculation frame:

$\begin{split}\vec{F}&=\frac{\partial
U_{\mathrm{Sk}}}{\partial\rho_{B}}\left[-\left(\vec{\nabla}\rho_{B}+\partial_{t}\vec{j}_{B}\right)+\dot{\vec{x}}\times\left(\vec{\nabla}\times\vec{j}_{B}\right)\right]\\\
&+\frac{\partial
U_{\mathrm{Sym}}}{\partial\rho_{I_{3}}}\left[-\left(\vec{\nabla}\rho_{I_{3}}+\partial_{t}\vec{j}_{I_{3}}\right)+\dot{\vec{x}}\times\left(\vec{\nabla}\times\vec{j}_{I_{3}}\right)\right]\end{split}$
(4)

$\vec{j}_{B}$ and $\vec{j}_{I_{3}}$ are the net-baryon current and the $I_{3}$
current respectively and $\dot{\vec{x}}$ is the velocity of the particle of
interest. Given the force in the calculation frame, the spacial components
$p_{i}$ of the particle momentum are updated at each time step

$p_{i}\rightarrow p_{i}+F_{i}\Delta t\,.$ (5)

The calculation therefore relies on small time steps.

Since the potentials and the equations of motion depend on the densities and
their spacial derivatives, a high resolution in the density is necessary. In
the calculations with potentials within this work, each particle is
represented by 20 test-particles for a more precise estimate of the densities
and their derivatives.

Different sets for the parameters of the Skyrme potential are described in
Kruse  _et al._ (1985). They are labeled soft and hard corresponding to their
stiffness. In addition, we include the default parameter set in SMASH Xu  _et
al._ (2016), which lies between the soft and hard one in terms of stiffness.
The parameter sets are explicitly given in Table 1.

In the future one should aim for a more continuous variation of parameters to
find the best values based on experimental data within a Bayesian multi-
parameter analysis. The momentum dependence will be added and the
implementation of relativistic mean fields based on density functional theory
is work in progress Sorensen and Koch (2020). In the current work, we do not
aim yet at a quantitative extraction of the nuclear equation of state, but
would like to provide a baseline for further studies and investigate the
influence of light cluster formation.

| soft | default | hard
---|---|---|---
$A$ |  $-365$ MeV |  $-209.2$ MeV |  $-124$ MeV
$B$ |  $303$ MeV |  $156.4$ MeV |  $71$ MeV
$\tau$ | $1.17$ | $1.35$ | $2.0$
$K$ |  $200$ MeV |  $240$ MeV |  $375$ MeV
Table 1: Parameter sets for Skyrme potential with corresponding
compressibilities $K$.

###  II.2 Clustering and light nuclei

At low beam energies, many protons are bound into light nuclei Reisdorf  _et
al._ (2010). Therefore light nuclei formation needs to be taken into account
from the theory side to be comparable with experiments. One option to do so is
to employ a coalescence model. Here one can directly calculate the spectrum of
nuclei given the proton spectra Butler and Pearson (1961); Schwarzschild and
Zupancic (1963); Gutbrod  _et al._ (1976) by assuming nucleons to form
clusters once they are sufficiently close in phase space. A similar approach
is to perform clustering on a microscopic basis in each event separately also
using the distance in phase space as a criterion for coalescence (see Sombun
_et al._ (2019)). Finally, one can explicitly form nuclei through production
cross sections from hadrons and treat the nuclei dynamically as active degrees
of freedom as described in Section II.3.

In the default SMASH setup, light nuclei are not treated explicitly as degrees
of freedom. Instead they are made up of individual nucleons that interact via
nuclear potentials. Despite the potentials, nucleons do not form stable bound
states in the final state of a heavy ion collision. We therefore follow the
method of identifying light nuclei as described in Sombun  _et al._ (2019),
where the condition for two particles to form a cluster is a small distance in
both momentum and coordinate space in the center of mass frame of the two
particles at the time where the last interaction occurs in which one of the
two particles took part. Compared to the simple clustering algorithm used in
Weil  _et al._ (2016) the main difference is the time at which the coordinate
space distance of particles is evaluated. The threshold for the momentum and
position distances cannot be adopted from Hillmann  _et al._ (2020), where the
improved method was used for a similar study, since the particles are not
represented by multiple test-particles in that work and the average distance
in phase-space is therefore not comparable. Instead, new thresholds $r_{0}$
and $p_{0}$ in coordinate and momentum space are found based on the proton
rapidity spectrum at $1.23A$ GeV using Bayesian parameter estimation where the
values $r_{0}=0.87^{+0.03}_{-0.03}$ fm and $p_{0}=0.43^{+0.03}_{-0.02}$ GeV
were found.

###  II.3 Deuterons in SMASH

A second option to consider the formation of light nuclei was introduced in
Oliinychenko  _et al._ (2019) and successfully applied to describe the
centrality dependence of deuteron formation at LHC energies Oliinychenko
(2020) as well as the energy dependence Oliinychenko  _et al._ (2020b). The
idea is to explicitly implement cross sections for nuclei production and treat
the nuclei as active degrees of freedom as previously done in Danielewicz and
Bertsch (1991); Oh  _et al._ (2009). The relevant processes for deuteron
formation including protons $p$, neutrons $n$, pions $\pi$ and generic
nucleons $N$ are

$\displaystyle pnN$ $\displaystyle\leftrightarrow$ $\displaystyle dN$ (6)
$\displaystyle pn\bar{N}$ $\displaystyle\leftrightarrow$ $\displaystyle
d\bar{N}$ (7) $\displaystyle pn\pi$ $\displaystyle\leftrightarrow$
$\displaystyle d\pi$ (8) $\displaystyle NN$ $\displaystyle\leftrightarrow$
$\displaystyle d\pi.$ (9)

The first two are $3\leftrightarrow 2$ processes that are in practice realized
with $2\leftrightarrow 1$ and $2\leftrightarrow 2$ processes since
multiparticle interactions are not yet present in SMASH. For the intermediate
step, a fictional dibaryon resonance $d^{\prime}\leftrightarrow pn$ is used
that can react in $2\leftrightarrow 2$ processes to a deuteron like
$Nd^{\prime}\leftrightarrow Nd$ or $\pi d^{\prime}\leftrightarrow\pi d$.

To investigate the relevant processes for deuteron formation Figure 1 shows
the relative number of reaction partners of the $d^{\prime}$ resonance in
percent in gold-gold collisions at $E_{\rm kin}=1.23A$ GeV. The numbers do not
sum up to 100% because in many cases the $d^{\prime}$ decays without a
collision partner. The nucleons clearly dominate the reaction partners, so the
$pnN\leftrightarrow dN$ process is as expected the most important at this
energy Danielewicz and Bertsch (1991).

Figure 1: Relative number of reaction partners of the fictional $d^{\prime}$
resonance in gold-gold collisions at $E_{\rm kin}=1.23A$ GeV. The $d^{\prime}$
often decays into $pn$ without a collision partner, therefore the numbers
don’t sum up to 100%.

At higher energies, the most important process is the one involving a pion
Oliinychenko  _et al._ (2019). This difference is caused by the fact that at
higher energies more pions are produced while at low energies a dense medium
of nuclear matter is produced.

In this work, we compare the simple coalescence approach with the explicit
deuteron formation to investigate which description is more favored by the
experimental data on flow measurements. The two approaches differ in the sense
that the coalescence takes the formation of arbitrary clusters into account
while at this point the explicit formation of only deuterons is implemented.
It is possible to extend the model also to larger nuclei but the complexity of
the problem increases with the mass number.

In Figure 2 a comparison of the two methods in terms of the $p_{T}$-spectra of
nucleons is shown. The curve labeled ”no clustering” refers to a calculation
without dynamic deuterons and without clustering.

Figure 2: Transverse momentum spectra of nucleons (full lines) and deuterons
(dashed lines) at mid-rapidity in 20%-30% central gold-gold collisions at
$E_{\rm kin}=1.23A$ GeV. For all calculations the hard equation of state is
used.

Naturally, the number of nucleons is the largest in this case. The clustering
result has the lowest number of nucleons because the formation of different
light nuclei can be taken into account as compared to the dynamic deuteron
calculation. The clustering mainly affects the low momentum part of the
spectrum, where the density in phase space is large. Comparing the dynamic
treatment of deuterons to the other curves, one can see that the deuteron
production influences the dynamics of nucleons significantly as nucleons are
shifted to larger $p_{T}$. This is related to interactions with deuterons that
also have large transverse momentum on average.

###  II.4 Flow calculations

In the following we present results for the flow coefficients of protons and
deuterons and compare to experimental data Kardan (2019); Adamczewski-Musch
_et al._ (2020). To be consistent with the measurement, the flow coefficients
are evaluated with respect to the reaction plane which is fixed to the x-z-
plane in the calculation setup. Hence, using the angle
$\phi=\arctan(p_{y}/p_{x})$ the $n$-th order flow coefficient can be evaluated
as an average over particles

$v_{n}=\langle\cos(n\phi)\rangle.$ (10)

The centrality of the events is selected by constraining the impact parameter
to a specific range. We focus on the 20% to 30% most central gold-gold
collisions, which correspond to the impact parameter range
$6.6\,\mathrm{fm}<b<8.1\,\mathrm{fm}$ Adamczewski-Musch  _et al._ (2018).

##  III Directed flow of protons and deuterons

We begin with the first order flow coefficient of protons and deuterons in
gold-gold collisions at $E_{\rm kin}=1.23A$ GeV. The selection of rapidity and
transverse momentum bins is chosen to cover a wide region in phase space.
Since flow coefficients in general are sensitive to the strength of the
potentials, results from calculations using different equations of state as
introduced in Section II.1 are compared. Figure 3 shows the $v_{1}$ of
nucleons as a function of rapidity for different bins in transverse momentum.
The shape reproduces the experimental measurements very well but only the
calculation with the hard equation of state gives the correct magnitude.

Figure 3: Directed flow of protons as a function of rapidity in 20%-30%
central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$ for
different $p_{T}$ bins compared to experimental data points Kardan (2019). The
full lines are obtained with a hard equation of state while for the dashed
lines the default equation of state is used.

Two calculations consider the formation of light nuclei by clustering in the
final state as described in Section II.2. In addition to the two equations of
state, a third calculation is shown, where no clustering is performed but the
deuterons are explicitly propagated as degrees of freedom in the calculation
as described in Section II.3. For that calculation, also the hard equation of
state is used.

Figure 4 shows the directed flow as a function of transverse momentum for
different rapidity bins. Also as a function of $p_{T}$, the hard equation of
state describes the experimental data best. The difference between the two
options to account for nuclei formation is clearly visible at low transverse
momentum but vanishes at larger $p_{T}$. The importance of clustering in the
low momentum region can be explained by the large phase-space density that
makes the formation of nuclei more likely.

Figure 4: Directed flow of protons as a function of transverse momentum in
20%-30% central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$
for different rapidity bins compared to experimental data points Kardan
(2019). The full lines are obtained with a hard equation of state while for
the dashed lines the default equation of state is used. The curve labeled ”No
clustering” takes the formation of light nuclei into account by explicitly
producing deuterons during the collision as described in Section II.3.

Compared to the data, clustering gives the best description and produces a
small kink in the low $p_{T}$ region, while the explicit deuterons give a
reasonable description of the data and the curve is more smooth but does not
follow the data as closely.

From the calculation where deuterons are explicitly propagated, the directed
flow of deuterons themselves is extracted. The results shown in Figure 5
compare the deuteron flow with and without nuclear potentials. For the
calculation with potentials, the hard equation of state is used.

Figure 5: Directed flow of deuterons as a function of transverse momentum in
20%-30% central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$
for different rapidity bins compared to experimental data points Kardan
(2019). A hard equation of state was employed here and deuterons were
dynamically treated as particles in this calculation.

Considering that the model was originally designed for high energy collisions
where the composition of the system is very different and nuclear potentials
are negligible, it is not obvious that the directed flow of deuterons can be
described. A surprisingly good agreement with the experimental data is
observed when the potentials are included.

To conclude the findings from the comparison to the measured $v_{1}$, the hard
equation of state is clearly favored in this setup. Overall a good agreement
with the data is observed, where the clustering setup follows a bit closer the
transverse momentum dependent $v_{1}$ of nucleons, while the deuteron flow is
matched extremely well in the calculation with explicit deuteron formation.

##  IV Elliptic flow of protons and deuterons

Continuing with the second order flow harmonic, Figure 6 shows the $v_{2}$ as
a function of rapidity in semi central gold-gold collisions. Again, two
calculations with clustering and different equations of state and a
calculation with explicit deuterons employing a hard equation of state are
compared.

Figure 6: Elliptic flow of protons as a function of rapidity in 20%-30%
central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$ for
different $p_{T}$ bins compared to experimental data points Kardan (2019). The
full lines are obtained with a hard equation of state while for the dashed
lines the default equation of state is used.

Naturally, the second order flow coefficient is more difficult to describe
than the first order one. In general, the magnitude of $v_{2}$ is too small at
large transverse momenta. However the flow in the intermediate $p_{T}$ region,
where most nucleons are located, is comparable to the experimental data. Same
as for $v_{1}$, the hard equation of state produces a stronger flow signal and
is therefore preferred by the data. In addition to that, the magnitude of
$v_{2}$ is larger when the deuterons are treated explicitly in the
calculation. That calculation agrees best with the data. Similar to the
observation for $v_{1}$, in the low $p_{T}$ region (here the
$0.4\,\mathrm{GeV}<p_{T}<0.45\,\mathrm{GeV}$ bin) the results are very
sensitive to how the formation of light nuclei is taken into account.

The elliptic flow of nucleons as a function of the transverse momentum is
shown in Figure 7. Here it is again rather obvious that the elliptic flow is
not well described, especially in the high $p_{T}$ region. In the low $p_{T}$
region, the clustering performs better, while the calculation with explicit
deuteron production works better at intermediate transverse momentum.

Figure 7: Elliptic flow of protons as a function of transverse momentum in
20%-30% central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$
for different rapidity bins compared to experimental data points Kardan
(2019). The full lines are obtained with a hard equation of state while for
the dashed lines the default equation of state is used.

The difficulties in describing the $v_{2}$ at high $p_{T}$ possibly comes from
the lack of momentum dependence of the potentials Welke  _et al._ (1988).

Figure 8 shows the elliptic flow of deuterons for the calculation, where they
are treated as active degrees of freedom. Same as for the directed flow, the
calculation matches the data pretty well.

Figure 8: Elliptic flow of deuterons as a function of rapidity in 20%-30%
central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$ for
different transverse momentum bins compared to experimental data points Kardan
(2019). A hard equation of state was employed here and deuterons were
dynamically treated as particles in this calculation.

Also looking at the $p_{T}$ differential elliptic flow of deuterons in Figure
9, the agreement with the data is good. At forward rapidity the calculation is
slightly off but at mid-rapidity the data is described almost perfectly.

Figure 9: Elliptic flow of deuterons as a function of transverse momentum in
20%-30% central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$
for different rapidity bins compared to experimental data points Kardan
(2019). A hard equation of state was employed here and deuterons were
dynamically treated as particles in this calculation.

To conclude this section, the elliptic flow is observed to be very sensitive
on how the light nuclei formation is taken into account, which leads to some
uncertainty for the extraction of the EoS. The elliptic flow of nucleons at
large transverse momenta is underestimated in the SMASH calculations but the
best agreement is found with the hard equation of state and treating the
deuterons as active degrees of freedom.

##  V Evolution of flow coefficients

In this section we focus on the evolution of flow with time to see when the
anisotropy is developed and what the most important stages are. We concentrate
on the setting that worked best in comparison to experimental data, namely the
calculation with explicit deuteron formation and a hard equation of state.

Figure 10: Nucleon and deuteron multiplicities, slope of the first order flow
coefficient at mid-rapidity and $v_{2}$ of nucleons and deuterons as a
function of time in 20% - 30% central gold-gold collisions at
$E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$. The flow coefficients are evaluated
for particles that satisfy $0.75\,\mathrm{GeV}<p_{T}<0.8\,\mathrm{GeV}$ while
the number of particles refers to the total multiplicity.

The first two panels of Figure 10 show the multiplicity of nucleons and
deuterons respectively as a function of time. The initially present nucleons
from the two gold nuclei start interacting, forming resonances and deuterons,
so that their number decreases. If no deuterons are present in the
calculation, the number of nucleons increases again towards the end of the
calculation when all resonances decay. The number of deuterons rises
throughout the evolution of the heavy ion collision. In the third panel one
can see the slope of the directed flow at mid-rapidity as a proxy for the
magnitude of $v_{1}$ of nucleons and deuterons. To evaluate the slope, at each
point in time the $v_{1}(y)$ is fit with a function $v_{1}(y)=ay+by^{3}$ where
the cubic term is added to take the curvature observed in Figure 3 into
account. The fitting is performed using Bayesian parameter estimation with
Markov chain Monte Carlo sampling. The inner band gives the 68% and the outer
band the 95% credible interval. At early times the uncertainty is very large
since not enough particles have interacted yet and are not located in the
relevant phase-space region and deuterons are not yet produced. Afterwards the
flow starts building up from zero to a maximum after $\approx 15$ fm. This is
the time the nuclei take to pass through each other. Afterwards one observes
for both nucleons and deuterons that the flow signal slightly weakens before
saturating but the effect is more pronounced for deuterons.

The final panel of Figure 10 shows a comparison of the elliptic flow at mid-
rapidity between nucleons and deuterons. Once enough particles are produced
and located in the investigated kinematic region the elliptic flow of nucleons
and deuterons look almost exactly the same. This is not intuitively clear
since no mass number scaling is applied here. Looking at the evolution of
$v_{2}$ in more detail one can see that at first a positive elliptic flow
signal starts to build up but soon the $v_{2}$ drops below zero where it stays
until the end of the evolution. Positive elliptic flow is typically associated
with pressure gradients in the initial state. At low collision energies a
competing effect is the squeeze out that results in a negative $v_{2}$ signal
when the slow spectator nuclei are blocking the path and push particles out of
the reaction plane Gutbrod  _et al._ (1989). Both effects contribute to the
observed elliptic flow.

It is interesting to note that especially for nucleons the $v_{2}$ shows a
very different behavior as a function of time than $v_{1}$. The difficulties
in describing the elliptic flow and directed flow with the same parameter set
might arise from the two being sensitive to different stages of the evolution.

##  VI Higher flow coefficients

###  VI.1 Scalar product $v_{3}$

In this section we would like to quantify the triangularity of heavy ion
collisions to extend the excitation function from Karpenko  _et al._ (2015)
down to $1.23A$ GeV. Unlike for the first and second harmonic, the $v_{3}$
signal does not simply emerge from the geometry of the colliding nuclei but
rather by fluctuations in the overlap region. To really quantify the
triangularity of an event it is therefore not sufficient to calculate the flow
coefficient with respect to the reaction plane. Hence for evaluating the
triangular flow in this section the scalar product method Adler  _et al._
(2002) is employed because even though the HADES collaboration only measures
$v_{3}(\Psi_{2})$ we would like to explore, if there is a finite triangular
flow in our calculations at low beam energies (see Appendix A for higher order
flow coefficients with respect to the reaction plane).

In the scalar product method the $n$th flow coefficient is calculated as
follows:

$v_{n}=\frac{\left\langle\vec{u}\cdot\frac{\vec{Q}_{n}}{N}\right\rangle}{\sqrt{\left\langle\frac{\vec{Q}_{n}^{A}}{N_{A}}\cdot\frac{\vec{Q}_{n}^{B}}{N_{B}}\right\rangle_{E}}}\,,$
(11)

where $\vec{u}$ is the momentum unit vector of a particle in the transverse
plane, $\langle...\rangle$ denotes the average over all particles of interest
while $\langle...\rangle_{E}$ is an average over events. $N$ is the total
number of particles and $N_{A}$ and $N_{B}$ are the number of particles in the
sub-events A and B respectively. The same labeling applies for the flow-
vectors $\vec{Q}^{A}_{n}$, $\vec{Q}^{B}_{n}$ and $\vec{Q}_{n}$ defined as

$\vec{Q}_{n}=\begin{pmatrix}\sum_{i}w_{i}\cos(n\phi_{i})\\\
\sum_{i}w_{i}\sin(n\phi_{i})\end{pmatrix}\,,$ (12)

where $\phi_{i}$ is the azimuthal angle of particle $i$, $w_{i}$ is the weight
which is for the $v_{3}$ calculation $p_{T}^{3}$. The sum runs over the
particles in the (sub-)event but for the scalar product in the numerator of
Equation 11 the particle of interest needs to be excluded from the flow vector
calculation to avoid auto correlations. Each event is divided into sub-events
according to the pseudo-rapidity $\eta$ of the particles. One sub-event
contains all particles with $\eta>0.1$ while the other one has particles that
satisfy $\eta<-0.1$.

Applying the scalar product method to extract the triangular flow of nucleons
one obtains the results shown in Figure 11.

Figure 11: Triangular flow of nucleons as a function of rapidity in 20%-30%
central gold-gold collisions. The flow coefficient is obtained using the
scalar product method.

Looking at the scale of the axis the triangular flow signal in collisions at
low energies almost vanishes. That agrees with the findings from Karpenko  _et
al._ (2015) where the triangular flow signal vanishes at low collision
energies. With the same cuts ($0.2\,\mathrm{GeV}<p_{T}<2.0\,\mathrm{GeV}$ and
$|\eta|<1$) applied as in Karpenko  _et al._ (2015) we obtain for nucleons and
deuterons $v_{3}=0.00081\pm 0.00006$ and $v_{3}=-0.0014\pm 0.0002$
respectively.

###  VI.2 Quadrangular flow

In this section we show the fourth order flow coefficient of nucleons divided
by the squared elliptic flow because this quantity is suggested to be a probe
of ideal hydrodynamic behavior of the system if $v_{4}/v_{2}^{2}\approx 0.5$
Borghini and Ollitrault (2006). $v_{4}/v_{2}^{2}$ is calculated at mid-
rapidity as a function of $p_{T}$ and the result is shown in Figure 12.

Figure 12: Quadrangular flow of nucleons divided by the squared elliptic flow
at mid-rapidity as a function of transverse momentum in 20% - 30% central
gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$.

One can see that $v_{4}/v_{2}^{2}$ is close to 0.5 in the $p_{T}$ region above
0.75 GeV. In general the shape of the curve is similar to the ideal three
dimensional hydro prediction at RHIC energies from Borghini and Ollitrault
(2006) so the pure transport calculation as presented in this work shows in
terms of this observable a close to hydrodynamic behavior. This is a sign that
in low energy collisions the response to the initial geometry within a
transport approach is translated as efficiently to final state observables as
expected from a hydrodynamic picture.

##  VII Conclusions and Outlook

In this paper we compared the double differential directed and elliptic flow
of protons and deuterons to experimental flow measurements in Au+Au collisions
at $E_{\rm lab}=1.23A$ GeV as provided by the HADES collaboration. Different
parameter sets for the Skyrme potential, that each correspond to a different
equation of state, are employed. In our calculation for the flow of both
nucleons and deuterons a hard equation of state is preferred by the data. The
overall agreement for the elliptic flow of nucleons is still not satisfactory.
Therefore more sophisticated potentials including a momentum dependence will
be necessary to extract the equation of state from flow measurements.

The main focus of the current work is the formation of light nuclei in low
energy heavy-ion collisions that cannot be neglected even considering just the
flow of nucleons. Two different ways of taking the formation of deuterons into
account are implemented. Within this model we observe that dynamically forming
deuterons and propagating them as particles throughout the calculation
performs better with respect to the HADES data than forming deuterons via
coalescence in the final state of a heavy-ion collision. In the setup with the
dynamic treatment of deuterons the evolution of the multiplicity and flow of
nucleons and deuterons over time is studied. Here one can see that the flow of
nucleons and deuterons show a very similar time dependence, while the directed
and elliptic flow coefficients are sensitive to different stages in the
evolution.

Finally we show for the first time that the scalar product triangular flow of
nucleons almost vanishes in collisions at low energy and confirm that
$v_{4}/v_{2}^{2}\approx 0.5$ outside the low-$p_{T}$ region within the SMASH
transport approach.

In the future a more advanced description of the nuclear potentials Sorensen
and Koch (2020) will allow for a more detailed comparison to a much broader
set of experimental measurememts. That comparison can be performed more
systematically using Bayesian methods to put reliable constraints on the
equation of state.

###### Acknowledgements.

Discussions and provision of experimental data is acknowledged to Behruz
Kardan and Christoph Blume. J. M. acknowledges funding by GSI F&E program.
Computational resources have been provided by the Center for Scientific
Computing (CSC) at the Goethe University and the GSI GreenCube. M. Mayer
acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the
grant CRC-TR 211 ”Strong-interaction matter under extreme conditions”. M.
Mayer would also like to thank Mathilde Ziegler-Himmelreich for fruitful
discussions.

##  Appendix A Reaction plane $v_{3}$ and $v_{4}$

For completeness we provide the triangular and quadrangular flow with respect
to the reaction plane in the same $p_{T}$ bin as shown in Adamczewski-Musch
_et al._ (2020). The hard equation of state is employed for all calculations
in this appendix. Figure 13 shows the first four flow coefficients of nucleons
as a function of rapidity. Compared is a calculation with explicit deuteron
formation as described in Section II.3 with final state clustering of nucleons
II.2.

Figure 13: $v_{1}$ to $v_{4}$ of nucleons as a function of rapidity in 20%-30%
central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$ for
$1.0\,\mathrm{GeV}<p_{T}<1.5\,\mathrm{GeV}$. All lines are obtained with a
hard equation of state. Compared is a calculation with explicit deuteron
formation and clustering.

We also provide the four lowest order flow coefficients calculation for
deuterons in Figure 14. The deuterons are explicitly produced and propagated
throughout the calculation.

Figure 14: $v_{1}$ to $v_{4}$ of deuterons as a function of rapidity in
20%-30% central gold-gold collisions at $E_{\mathrm{kin}}=1.23A\,\mathrm{GeV}$
for $1.0\,\mathrm{GeV}<p_{T}<1.5\,\mathrm{GeV}$. All lines are obtained with a
hard equation of state.

## References

  * Scheid  _et al._ (1974) W. Scheid, H. Muller, and W. Greiner, Phys. Rev. Lett. 32, 741 (1974).
  * Danielewicz  _et al._ (2002) P. Danielewicz, R. Lacey, and W. G. Lynch, Science 298, 1592 (2002), arXiv:nucl-th/0208016 .
  * Gale  _et al._ (1987) C. Gale, G. Bertsch, and S. Das Gupta, Phys. Rev. C 35, 1666 (1987).
  * Reisdorf  _et al._ (2012) W. Reisdorf _et al._ (FOPI), Nucl. Phys. A 876, 1 (2012), arXiv:1112.3180 [nucl-ex] .
  * Aichelin  _et al._ (1987) J. Aichelin, A. Rosenhauer, G. Peilert, H. Stoecker, and W. Greiner, Phys. Rev. Lett. 58, 1926 (1987).
  * Isse  _et al._ (2005) M. Isse, A. Ohnishi, N. Otuka, P. Sahu, and Y. Nara, Phys. Rev. C 72, 064908 (2005), arXiv:nucl-th/0502058 .
  * Le Fèvre  _et al._ (2016) A. Le Fèvre, Y. Leifels, W. Reisdorf, J. Aichelin, and C. Hartnack, Nucl. Phys. A 945, 112 (2016), arXiv:1501.05246 [nucl-ex] .
  * Le Fèvre  _et al._ (2018) A. Le Fèvre, Y. Leifels, C. Hartnack, and J. Aichelin, Phys. Rev. C 98, 034901 (2018), arXiv:1611.07500 [nucl-th] .
  * Petersen  _et al._ (2006) H. Petersen, Q. Li, X. Zhu, and M. Bleicher, Phys. Rev. C 74, 064908 (2006), arXiv:hep-ph/0608189 .
  * Hillmann  _et al._ (2018) P. Hillmann, J. Steinheimer, and M. Bleicher, J. Phys. G 45, 085101 (2018), arXiv:1802.01951 [nucl-th] .
  * Hillmann  _et al._ (2020) P. Hillmann, J. Steinheimer, T. Reichert, V. Gaebel, M. Bleicher, S. Sombun, C. Herold, and A. Limphirat, J. Phys. G 47, 055101 (2020), arXiv:1907.04571 [nucl-th] .
  * Nara  _et al._ (2017) Y. Nara, H. Niemi, J. Steinheimer, and H. Stöcker, Phys. Lett. B 769, 543 (2017), arXiv:1611.08023 [nucl-th] .
  * Abbott  _et al._ (2016) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Bauswein  _et al._ (2019) A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Phys. Rev. Lett. 122, 061102 (2019), arXiv:1809.01116 [astro-ph.HE] .
  * Most  _et al._ (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120, 261103 (2018), arXiv:1803.00549 [gr-qc] .
  * Wang  _et al._ (2020) Y. Wang, Q. Li, Y. Leifels, and A. Le Fèvre, Phys. Lett. B 802, 135249 (2020).
  * Tsang  _et al._ (2012) M. Tsang _et al._ , Phys. Rev. C 86, 015803 (2012), arXiv:1204.0466 [nucl-ex] .
  * Butler and Pearson (1961) S. Butler and C. Pearson, Phys. Rev. Lett. 7, 69 (1961).
  * Schwarzschild and Zupancic (1963) A. Schwarzschild and C. Zupancic, Phys. Rev. 129, 854 (1963).
  * Gutbrod  _et al._ (1976) H. Gutbrod, A. Sandoval, P. Johansen, A. M. Poskanzer, J. Gosset, W. Meyer, G. Westfall, and R. Stock, Phys. Rev. Lett. 37, 667 (1976).
  * Sombun  _et al._ (2019) S. Sombun, K. Tomuang, A. Limphirat, P. Hillmann, C. Herold, J. Steinheimer, Y. Yan, and M. Bleicher, Phys. Rev. C 99, 014901 (2019), arXiv:1805.11509 [nucl-th] .
  * Oh  _et al._ (2009) Y. Oh, Z.-W. Lin, and C. M. Ko, Phys. Rev. C 80, 064902 (2009), arXiv:0910.1977 [nucl-th] .
  * Aichelin  _et al._ (2020) J. Aichelin, E. Bratkovskaya, A. Le Fèvre, V. Kireyeu, V. Kolesnikov, Y. Leifels, V. Voronyuk, and G. Coci, Phys. Rev. C 101, 044905 (2020), arXiv:1907.03860 [nucl-th] .
  * Weil  _et al._ (2016) J. Weil _et al._ , Phys. Rev. C 94, 054905 (2016), arXiv:1606.06642 [nucl-th] .
  * Oliinychenko  _et al._ (2020a) D. Oliinychenko, V. Steinberg, J. Weil, M. Kretz, J. Staudenmaier, S. Ryu, A. Schäfer, J. Rothermel, J. Mohs, F. Li, H. E. (Petersen), L. Pang, D. Mitrovic, A. Goldschmidt, L. Geiger, J.-B. Rose, J. Hammelmann, and L. Prinz, “smash-transport/smash: Smash-1.8,”  (2020a).
  * Kardan (2019) B. Kardan (HADES), Nucl. Phys. A 982, 431 (2019), arXiv:1809.07821 [nucl-ex] .
  * Adamczewski-Musch  _et al._ (2020) J. Adamczewski-Musch _et al._ (HADES), (2020), arXiv:2005.12217 [nucl-ex] .
  * Xu  _et al._ (2016) J. Xu _et al._ , Phys. Rev. C 93, 044609 (2016), arXiv:1603.08149 [nucl-th] .
  * Sorensen and Koch (2020) A. Sorensen and V. Koch, (2020), arXiv:2011.06635 [nucl-th] .
  * Oliinychenko  _et al._ (2019) D. Oliinychenko, L.-G. Pang, H. Elfner, and V. Koch, Phys. Rev. C 99, 044907 (2019), arXiv:1809.03071 [hep-ph] .
  * Tanabashi  _et al._ (2018) M. Tanabashi _et al._ (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Manley and Saleski (1992) D. Manley and E. Saleski, Phys. Rev. D 45, 4002 (1992).
  * Kruse  _et al._ (1985) H. Kruse, B. Jacak, and H. Stoecker, Phys. Rev. Lett. 54, 289 (1985).
  * Reisdorf  _et al._ (2010) W. Reisdorf _et al._ (FOPI), Nucl. Phys. A 848, 366 (2010), arXiv:1005.3418 [nucl-ex] .
  * Oliinychenko (2020) D. Oliinychenko, in  _28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions_ (2020) arXiv:2003.05476 [hep-ph] .
  * Oliinychenko  _et al._ (2020b) D. Oliinychenko, C. Shen, and V. Koch, (2020b), arXiv:2009.01915 [hep-ph] .
  * Danielewicz and Bertsch (1991) P. Danielewicz and G. Bertsch, Nucl. Phys. A 533, 712 (1991).
  * Adamczewski-Musch  _et al._ (2018) J. Adamczewski-Musch _et al._ (HADES), Eur. Phys. J. A 54, 85 (2018), arXiv:1712.07993 [nucl-ex] .
  * Welke  _et al._ (1988) G. Welke, M. Prakash, T. Kuo, S. Das Gupta, and C. Gale, Phys. Rev. C 38, 2101 (1988).
  * Gutbrod  _et al._ (1989) H. Gutbrod, B. Kolb, H. Schmidt, A. M. Poskanzer, H. Ritter, and K. Kampert, Phys. Lett. B 216, 267 (1989).
  * Karpenko  _et al._ (2015) I. Karpenko, P. Huovinen, H. Petersen, and M. Bleicher, Phys. Rev. C 91, 064901 (2015), arXiv:1502.01978 [nucl-th] .
  * Adler  _et al._ (2002) C. Adler _et al._ (STAR), Phys. Rev. C 66, 034904 (2002), arXiv:nucl-ex/0206001 .
  * Borghini and Ollitrault (2006) N. Borghini and J.-Y. Ollitrault, Phys. Lett. B 642, 227 (2006), arXiv:nucl-th/0506045 .