Harmonic flow correlations in Au+Au reactions at 1.23 AGeV: A new testing ground for the Equation-of-State and expansion geometry

Correlations between the harmonic flow coefficients $v_1$, $v_2$, $v_3$ and $v_4$ of nucleons in semi-peripheral Au+Au collisions at a beam energy of 1.23~AGeV are investigated within the hadronic transport approach Ultra-relativistic Quantum Molecular Dynamics (UrQMD). In contrast to ultra-relativistic collision energies (where the flow coefficients are evaluated with respect to the respective event plane), we predict strong correlations between the flow harmonics with respect to the reaction plane. Based on an event-by-event selection of the midrapidity final state elliptic flow of protons we show that as a function of rapidity, I) the sign of the triangular flow changes, II) that the shape of $v_4$ changes from convex to concave, and III) that $v_3\propto v_1v_2$ and $v_4\propto v_2^2$ for all different event classes, indicating strong correlations between all investigated harmonic flow coefficients.


I. INTRODUCTION
The exploration of the properties of dense and hot nuclear matter is an ongoing endeavor since more than 40 years. Nowadays, the main focus lies in two areas: I) The systematic investigation of the phase diagram of Quantum-Chromo-Dynamics to understand the onset of deconfinement [1] and II) the extraction of the nuclear matter Equation-of-State (EoS) at moderate temperatures and its relation to astrophysical objects [2]. On the experimental side, the Large Hadron Collider (LHC) probes the strong interaction at the highest available collision energies. Here, a deconfined and nearly net-baryon density free system at very high temperatures is created, and such an environment is theoretically well accessible by calculations of Quantum-Chromo-Dynamics (QCD) on the lattice. At this high-temperature frontier the Fourier decomposition of the azimuthal angle distribution of emitted particles can shed light on the expansion dynamics of the created matter. Especially, the elliptic flow (usually called v 2 ) is a fabric of pressure gradients in the transverse direction and allows to extract the viscosity (and other transport coefficients) of the Quark-Gluon-Plasma (QGP). The studies at the LHC extend and complement previous and ongoing measurements at Relativistic Heavy Ion Collider (RHIC) [3,4], which have pioneered the extraction of the viscosity of a Quark-Gluon-Plasma (QGP) and have lead to the strongly coupled perfect liquid hypothesis of the QGP [5][6][7][8]. Further studies at RHIC and at LHC on the triangular (v 3 ) and quadrangular (v 4 ) flow have tied both coefficients to the initial state fluctuations [9][10][11] which can provide further insights on the parton structure of the impinging nuclei [12].
With decreasing energy, the high-density frontier is probed.
Here prominent examples are the RHIC beam energy scan program (RHIC-BES) [13] and the future FAIR [14] and NICA [15] facilities. In this energy domain, the onset of deconfinement is expected and a gradual transition (with increasing beam energy) to a QGP is expected. At even lower beam energies, e.g. currently actively explored by the GSI facility, the nuclear Equationof-State (EoS) of hadronic matter [16] is explored and will allow to bridge the gap to binary neutron star mergers, which might provide complementary information on the EoS via the detection of gravitational waves [17]. In this energy regime, the Fourier decomposition of the azimuthal angle distribution of the created particles is driven by an intricate interplay between the EoS [18,19], time dependent expansion dynamics/geometry [20] and viscous corrections [21]. Precise measurements up to the sixth flow coefficient of protons and light clusters (d, t) have already been successfully accomplished by the HADES collaboration at GSI [22].
A novel tool to explore the properties and geometry of the created QCD matter, namely harmonic flow correlations, was recently suggested and measured in ultrarelativistic nucleus-nucleus reactions at RHIC and at LHC [23,24]. At both energies an anti-correlation between event-by-event fluctuations of v 2 and v 3 was observed, while the event-by-event fluctuations of v 2 and v 4 were found to be correlated. These results were also confirmed by hydrodynamic [25] and transport simulations [26,27]. Generally, the (anti-)correlations were tracedback to the initial state eccentricities and it was concluded that the investigation of the correlations between different flow harmonics have a substantially higher sensitivity to the details of theoretical calculations (transport coefficients, equation-of-state, initial state modeling) in comparison to individual v n coefficients.
In this article we explore for the first time the eventby-event correlations between the first four harmonic flow coefficients (v 1 to v 4 ) to extract further information on the properties and geometries of the matter created at low beam energies. We propose to test these predictions with the currently running HADES experiment.

II. MODEL SETUP AND FLOW EXTRACTION
For the present study we use the Ultra-relativistic Quantum Molecular Dynamics (UrQMD) model [28][29][30] in its most recent version (v3.5). UrQMD is a dynamical microscopic transport simulation based on the explicit propagation of hadrons in phase-space. The imaginary part of the interactions is modeled via binary elastic and inelastic collisions, leading to resonance excitations and decays or color flux-tube formation and their fragmentation. The real part of the interaction potential is implemented via different equations of state (following the usual notion of a hard and soft equation-of-state), alternative equations-of-state, e.g. a chiral mean field EoS can also be introduced, see e.g. [31]. In its current version, UrQMD includes a large body of baryonic and mesonic resonances up to masses of 4 GeV. The model is well established in the GSI energy regime. For recent studies of the bulk dynamics, we refer the reader to [32,33]. For the analysis of the integrated harmonic flows at GSI energies see [19,34].
The flow coefficients are identified with the Fourier coefficients in the series expansion of the azimuthal angle distribution which can be written as in which v n is n-th order flow coefficient of the even (cosine) term,ṽ n is the n-th order flow coefficient of the odd (sine) term, φ is the azimuthal angle and Ψ RP is the angle of the reaction plane. The HADES experiment uses a forward wall to reconstruct the event plane from the spectator nucleons [22,35]. In the simulation, the reaction plane is fixed because the impact parameter is known and thus Ψ RP = 0 is used for the present analysis of the simulation. However, the first order spectator event plane can still fluctuate on an event by event basis which is reflected in general in nonzero sine terms for a single event. But by taking the average over all events (in a given event class) they will be zero due to symmetry. We hence restrict our investigation to the cosine terms (i.e. the projection of the flow vector of each particle onto the known reaction plane). The flow coefficients are then calculated as v n = cos(n(φ − Ψ RP )) , where the average · is taken over all nucleons in a fixed rapidity or transverse momentum range in a given event.
Let us stress that the extraction of the flow coefficents in the HADES experiment is different from the methods used at higher energies. At high energies, e.g. at RHIC and LHC, the flow coefficients are usually extracted with respect to the n-th order event plane which is constructed via the flow vector Q [36,37]. Another possibility is to extract the flow from the 2-and 4-particle cumulants [38] or to use the Lee-Yang zeroes method [39]. Generally such methods allow to obtain the magnitude of the flow coefficients and work best (especially in case of the Lee-Yang zero method), if the multiplicites are sufficiently high and the flow harmonics are not too small [40]. At low energies the situation is different, the multiplicities are limited, the flow coefficients might be rather small and the sign of the elliptic flow is important (because of the change from in-plane to out-of-plane emission (squeeze-out) due to the blocking by the spectators towards lower energies). Therefore, the HADES experiment uses a different method and extracts the coefficients with respect to the reaction plane. In the simulation this can be done exactly, because the reaction plane is known. However in the experiment the reaction plane has to be reconstructed. The HADES collaboration uses the first order spectator event plane as a good estimator for the reaction plane with high resolution [22,35].
Previously, the flow correlations have been quantified using Pearson correlation functions [25]. However, for the present study we focus mainly on different event classes straightforwardly and show the correlation directly, because showing the full distributions allows for a more direct interpretation of the results than compressing the distribution into a single number. However, at the end of the paper, we also provide the Pearson coefficients to allow for a comparison to the data at higher energies.

III. RESULTS
All results were obtained by simulating 20-30% peripheral Au+Au collisions at E beam = 1.23 AGeV kinetic beam energy with the UrQMD (v3.5) model. The centrality is selected via impact parameter cuts following previous Monte-Carlo Glauber simulations [41]. We employ the model mainly with a hard equation-of-state as it was shown to yield the best description of the measured HADES data (cf. Ref. [19,34]). We focus our analysis on participating nucleons and also exclude nucleons that are bound in light clusters. It has been shown [34] that both effects need to be taken into account to reliably describe the measured data [22].

A. Elliptic flow fluctuations and event class selection
We start our investigation with an analysis of the event-by-event distribution of the final state elliptic flow. ons in the rapidity window 1 |y cm | ≤ 0.5 in 20-30% central Au+Au collisions at kinetic beam energy of 1.23 AGeV from UrQMD with a hard EoS. We observe that the event-by-event value of the elliptic flow fluctuates substantially around the mean value of v 2 |ycm|≤0.5 = -0.05 (the mean value is in line with the HADES data [22] and also consistent with older data in a similar energy region [42][43][44]). The width of the v 2 distribution is rather broad and has a FWHM of 0.15. Thus, single events can even show an overall positive v 2 , as well as highly negative v 2 values. This finding naturally gives rise to the idea that heavy-ion collision events can be categorized into different classes based on their final elliptic flow around midrapidity, while keeping everything else fixed. Such a categorization is well known from ultra-relativistic collisions under the name "event shape engineering" [45] connecting event wise final flow to the initial geometric (spatial) configuration. Although the correlation between initial spatial fluctuations and final elliptic flow is less obvious at low energies than at ultra-relativistic energies, due to the intricate emission dynamics between spectator blocking and expansion, the observation of strong final state elliptic flow fluctuations opens new possibilities for quantitative tests of the properties of high density QCD matter and the dynamical models used for its investigation.

B. Flow in different event classes
Let us now select event classes based on the eventwise elliptic flow at midrapidity, but keeping everything else (collision system, energy and centrality) fixed. To this aim, we split the events into 14 classes defined by the event wise integrated elliptic flow at midrapidity, v 2 |ycm|≤0.5 with the specific values for the (momentumspace-)ellipticity given in Tab. I. event class ellipticity For each of these event classes, we now investigate the first to fourth order flow coefficients extracted in each specified event class as shown in Fig. 2  Let us begin the discussion with the directed flow v 1 . We observe that the directed flow is only marginally affected by the selection of the event class. This is expected, because the directed flow is mainly driven by the bounce-off of the impinging nuclei. Nevertheless, a small dependence on the v 2 event class can be observed with a positive correlation leading to a larger magnitude of v 1 at forward/backward rapidities (i.e. an increase of dv 1 /dy cm | ycm=0 ) with the selection of higher v 2 values. Turning to the elliptic flow, we observe that the event class selection via the final state v 2 shifts the magnitude of the elliptic flow, while its qualitative rapidity dependence stays unaffected. The major effects are however found when exploring the triangular flow and its correlation to v 2 . The magnitude and rapidity dependence of v 3 vary drastically with different v 2 event classes. While for averaged (i.e. unselected) events one finds a slightly negative v 2 and also a v 3 with a slightly negative slope at midrapidity, for positive elliptic flow event classes the triangular flow develops a strongly positive slope. Also for the opposite event class (strongly negative v 2 ) we observe a correlation indicated by a steeper negative slope of v 3 for such events. This clearly shows the strong correlation between v 2 and v 3 in the HADES energy regime. Finally, we address the rapidity dependence of the quadrangular flow coefficient v 4 . Also v 4 shows an interesting and unexpected behavior: While events in the negative v 2 event class show a concave shape of v 4 as function of rapidity, event classes with a positive v 2 tend to develop a convex shape in rapidity. This indicates that also v 2 and v 4 are intertwined even on an event-by-event basis and show a strong anti-correlation (larger v 2 leads to smaller v 4 at midrapidity). One should note that the correlations among the flow coefficients are strikingly different at low energies in comparison to high energies and also much stronger pronounced at the low energies in- during the overlap phase leading to negative v 2 , followed by an expansion in impact parameter direction resulting in a positive v 2 ) pattern which are fundamentally different than at high energies where the expansion in the impact parameter direction dominates already from the beginning.

C. Event class selected scaling of flow coefficients in rapidity
With this splitting in different event classes in hand we investigate now the multiplicative relations between the lower flow harmonics with the higher order harmonics that were previously found, namely v 3 ≈ v 1 · v 2 and v 4 ≈ 0.5v 2 2 (we refer to [46] for a deeper discussion of the v 4 scaling relation). Especially the scaling of v 3 is interesting, because at ultra-relativistic collision energies the triangular flow has been shown to be mainly sensitive to the initial state [9][10][11], while at the low energies investigated here where flow is extracted with respect to the first order event plane, triangular flow is mostly attributed to geometry and the intricate space and time dependent emission pattern of hadrons [19,47]. Fig. 3 shows the rapidity dependence of triangular flow v 3 (upper left) in comparison to the multiplicative relation v 1 · v 2 (upper right) and the rapidity dependence of quadrangular flow v 4 (lower left) in comparison to the scaling relation 0.5v 2 2 (lower right) in the different event classes (from the most positive v 2 class in red to most negative v 2 class in purple, see legend). The integrated flow coefficient, i.e. without event class selection is shown as a black line. All calculations are done for 20-30% peripheral Au+Au collisions with hard EoS at 1.23 AGeV. For the triangular flow, a nearly perfect matching with the product of directed and elliptic flow in all event classes can be observed. A similar observation can be made for the scaling of the quadrangular flow with the square of the elliptic flow. One should note that this scaling is remarkable because the elliptic flow changes sign as a function of event class, but still its square stays proportional to the v 4 flow harmonic and even the change of shape in rapidity from convex to concave is reproduced. Nevertheless a systematic numerical scaling factor is necessary which can be quantified as v 4 /v 2 2 ≈ 0.5. How can we interpret these results? First of all, the results demonstrate that higher flow harmonics (n > 2) are strongly intertwined with the first and second flow harmonics emphasizing the influence of geometry. In addition the quadrangular flow and its scaling with the elliptic flow was suggested to provide specific information on the applicability of ideal fluid dynamics [48]. In [48] the authors argued, using a hydrodynamic model, that a scaling relation of v 4 = 0.5v 2 2 , as found in the present studies (see also [34,46]), suggests applicability of ideal hydrodynamics, which means local equilibrium seems to be achieved. It should however be noted that such a conclusion is at variance with a study of the viscosity at HADES energies suggesting a substantial time dependent viscosity over entropy ratio [21] and also not supported by the analysis of Ref. [46].

D. Event class selected analysis of flow shapes at midrapidity
Let us finally use our event class analysis to obtain more information on the Equation-of-State. We start from the correlation of the v 2 shape in rapidity in dependence of the v 2 trigger. Specifically we propose to use the curvature of the elliptic flow around midrapidity, i.e. d 2 v 2 /dy 2 cm | ycm=0 . Numerically, the curvature is obtained by fitting the rapidity dependence of v 2 with a quadratic polynomial and extracting the curvature from the fit. Fig. 4 shows the curvature of the elliptic flow component v 2 of nucleons at y cm = 0 as a function of the final state elliptic flow v 2 (colored full symbols) integrated over −0.5 ≤ y cm ≤ 0.5 for a hard EoS (circles), a soft EoS (triangles) and in cascade mode (squares) from 20-30% peripheral Au+Au collisions at 1.23 AGeV from UrQMD. The curvature of the elliptic flow is increasing with increasing trigger v 2 for the hard and soft equationsof-state while the curvature in case of the cascade simulation is decreasing with increasing v 2 event class. This suggests that the curvature of the elliptic flow is strongly sensitive to the nuclear EoS.
Next, we turn to the triangular flow shape and its correlation with the v 2 trigger. We specifically propose to use the correlation between the final v 2 and the slope dv 3 /dy cm | ycm=0 at midrapidity for the different event classes. Given our scaling assumption one would assume in first approximation that dv 3 /dy cm | ycm=0 = dv 1 /dy cm | ycm=0 · v 2 . Thus we expect to observe a linear dependence of dv 3 /dy cm | ycm=0 on v 2 . In Fig. 5 we show dv 3 /dy cm | ycm=0 as a function of the v 2 event class for the nucleons in semi-peripheral Au+Au collisions at 1.23 AGeV for cascade calculations (no potential, squares), a soft EoS (triangles) and a hard EoS (circles) as full colored symbols. For all studied equations-of-state a linear dependence between v 2 and dv 3 /dy cm | ycm=0 is observed. The slope of the correlation depends on the stiffness of the EoS, a stiffer EoS shows a stronger incline. This correlation observable allows therefore also to pin down the equation-of-state more precisely than is usually possible.
Last, we turn to the quadrangular flow. Motivated by measurements of the ATLAS collaboration [49], we analyze the quadrangular flow at midrapidity as a function of the elliptic flow trigger. Fig. 6 shows the quadrangular flow component v 4 of nucleons at y cm = 0 as a function of the final state elliptic flow v 2 (colored full symbols) integrated over −0.5 ≤ y cm ≤ 0.5 for a hard EoS (circles), a soft EoS (triangles) and in cascade mode (squares) from 20-30% peripheral Au+Au collisions at 1.23 AGeV from UrQMD. The lines depict least-squares quadratic fits. We observe a quadratic dependence of v 4 on the v 2 trigger as seen by the quadratic polynomial fit functions. The coefficients of the quadratic terms are 0.521 (hard EoS), 0.521 (soft EoS) and 0.522 (cascade) which again confirm the scaling relation v 4 ∝ 0.5v 2 2 . However, the value of v 4 alone is not sensitive to the EoS. Therefore, we perform a similar analysis as for the elliptic flow and extract the v 4 shape in rapidity as a function of the v 2 trigger. Again, we use a quadratic polynomial and extract the curvature from the fit, i.e. d 2 v 4 /dy 2 cm | ycm=0 . Fig. 7 shows the curvature of the quadrangular flow of nucleons at y cm = 0 as a function of the final state elliptic flow v 2 (colored full symbols) integrated over −0.5 ≤ y cm ≤ 0.5 for a hard EoS (circles), a soft EoS (triangles) and in cascade mode (squares) from 20-30% peripheral Au+Au collisions at 1.23 AGeV from UrQMD. We observe a strong splitting of the quadrangular flow shapes towards negative trigger v 2 , which indicates a strong dependence on the nuclear equation of state.

IV. FLOW CORRELATION FUNCTIONS
To bridge the gap to previous studies of flow fluctuations and their correlations [25][26][27] we finally investigate the correlation functions among the flow harmonics. To this aim, we use the linear correlation function corr(v n , v n ) (also known as the Pearson coefficient 2 ) between the first four flow harmonics calculated as Here, the standard deviation σ vi = v 2 i − v i 2 is used to normalize the covariance. In Fig. 8 we show the Pearson correlation function corr(v n , v m ) (full symbols) between the first four flow harmonics of nucleons as a function of rapidity for a hard EoS (circles), a soft EoS (triangles) and in cascade mode (squares) from 20-30% peripheral Au+Au collisions at 1.23 AGeV from UrQMD.
We observe a strongly pronounced rapidity dependent correlation between the first and second, the second and third, and the third and forth flow harmonic which is point-symmetric around y cm = 0, while the correlation between the first and third, and second and fourth flow harmonic is symmetric around y cm = 0 but its value is negative and smaller, even vanishing for the cascade simulation. The correlation of the first and fourth flow harmonic is negligibly small. Especially, the previously observed scaling relation of quadrangular flow (v 4 ∝ v 2 2 ) is interesting to discuss. The correlation is clearly observed in our calculations and has also been extensively studied in various publications [34,46,48], but due to its quadratic dependence the Pearson correlation shows no significant signal, but only a maximal value of -0.06. Although the Pearson correlation allows at most to draw conclusions about linear dependence, the results demonstrate nonetheless a dependence of the Pearson coefficient on the employed Equation-of-State: the stiffer the EoS the stronger the correlation becomes for all harmonic flow combinations. Together with the slope of triangular flow and the curvature of elliptic and quadrangular flow shown above this composes a neat possibility to pin down the nuclear Equation-of-State with high precision.

V. CONCLUSION
We have employed the Ultra-relativistic Quantum-Molecular-Dynamics model (UrQMD v3.5) to study semi-peripheral Au+Au collisions (20-30% centrality) at a beam energy of 1.23 AGeV. We found that the final state v 2 of a single event fluctuates strongly around its mean value allowing to identify classes of events with selected elliptic flow. We employ these event classes as a trigger to extract the rapidity dependence of the first to fourth order harmonic flow coefficients. Using this trigger, the directed flow is only slightly affected in its magnitude, while the elliptic flow acquires (as expected) a linear shift with positive correlation. The higher order flow components however reveal interesting new features, namely as a function of rapidity the triangular flow changes sign and the quadrangular flow changes shape from concave to convex when going from negative to positive v 2 event classes. All event classes are found to fulfill the scaling relation v 3 ∝ v 1 v 2 nearly perfectly, also the quadrangular flow exhibits the ideal fluid scaling v 4 = 0.5v 2 2 around midrapidity in all event classes. We further demonstrated that, correlation between the slope of the triangular flow and the curvature of the elliptic and quadrangular flow at midrapidity as a function of the selected elliptic flow are highly sensitive to the EoS. We thus propose these features as novel tools to obtain further information on the nuclear Equation-of-State. Lastly, we investigated the Pearson coefficients which underline the correlations analyzed with the event class selection and allow to compare the current analysis to calculations at higher collision energies.