Differential HBT Applied to Relativistic Fluid Dynamics
NNew method to detect rotation in high energy heavy ion collisions
L.P. Csernai , S. Velle , and D.J. Wang , Institute of Physics and Technology, University of Bergen, Allegaten 55, 5007 Bergen, Norway Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China. (Dated: November 10, 2018)With increasing beam energies the angular momentum of the fireball in peripheral heavy ioncollisions is increasing, and the proposed Differential Hanbury Brown and Twiss analysis is able toestimate this angular momentum quantitatively. The method detects specific space-time correlationpatterns, which are connected to rotation.
PACS numbers: 25.75.-q, 25.75.Gz, 25.75.Gz, 25.75.Nq
I. INTRODUCTION
In high energy peripheral heavy ion collisions the highangular momentum is realized in rotating flow, large ve-locity shear, vorticity and circulation. Viscous, explo-sive expansion leads to the decrease of vorticity and cir-culation with time, however, with small viscosity thevorticity remains significant at the final freeze-out (FO)stages. The proposed Differential Hanbury Brown andTwiss (HBT) analysis, a combination of standard twoparticle correlation functions, is adequate to analyze ro-tating systems. At the present collision energies the an-gular momentum and rotation is becoming a dominantfeature of reaction dynamics, and up to now the rotationof the system was never analysed, neither with the HBTmethod nor in any other way. We present and analysethis method and its results in a high resolution, ParticleIn Cell fluid dynamics model. Fluid dynamics is provento be the best theoretical method to describe collectiveflow phenomena. The same model was used to predictthe rotation in peripheral ultra-relativistic reactions [1],to point out the possibility of Kelvin Helmholtz Instabil-ity (KHI) [2], flow vorticity [3] and polarization arisingfrom local rotation, i.e. vorticity [4]. The model wasalso tested for its numerical viscosity and the resultingentropy production [5]. The formation of KHI was alsoobserved recently in AdS/CFT holography, where the in-stability is even more pronounced in peripheral reactions,although the time scale is sufficiently short only at highquark chemical potentials as at FAIR, NICA and RHIC-BES [6].The total angular momentum of the fireball is maximalat b = 0 . b max [7], while the angular momentum per netbaryon charge is maximal around b = (0 . − . b max .At ultra-peripheral collisions fluctuations dominate col-lective effects. According to the present analysis the Dif-ferential HBT method is indicating rotation via particlesat collective momenta, p t ≈ (0 . −
2) GeV/c the best, andthe magnitude of the introduced Differential CorrelationFunction is monotonically increasing with the angularmomentum.
II. CORRELATION FUNCTION FROM FLUIDDYNAMICS
The pion correlation function is defined as the inclusivetwo-particle distribution divided by the product of theinclusive one-particle distributions, such that [8]: C ( p , p ) = P ( p , p ) P ( p ) P ( p ) , (1)where p and p are the 4-momenta of pions. We intro-duce the center-of-mass momentum : k = ( p + p ) , and the relative momentum q = p − p , where from themass-shell condition [8] q = kq /k . We use a methodfor moving sources presented in Ref. [9]. C ( k, q ) = 1 + R ( k, q ) (cid:12)(cid:12)(cid:82) d x S ( x, k ) (cid:12)(cid:12) , (2)where R ( k, q ) = (cid:90) d x d x cos[ q ( x − x )] × S ( x , k + q/ S ( x , k − q/ . (3)Using the emission function S ( x, k ), discussed in refs.[10], here R ( k, q ) can be calculated [9] via the function J ( k, q ) = (cid:90) d x S ( x, k + q/
2) exp( iqx ) , (4)and we obtain the R ( k, q ) function as R ( k, q ) = Re [ J ( k, q ) J ( k, − q )].We estimate the local pion density by the specific en-tropy, σ ( x ), as n π ( x ) ∝ n ( x ) σ ( x ), where n ( x ) is the The vector k is the wavenumber vector, k = p / (cid:126) so for numericalcalculations we have to use that (cid:126) c = 197.327 MeV fm., The sameapplies to q . a r X i v : . [ nu c l - t h ] M a r proper net baryon charge density. Thus the local in-variant pion density is given by the J¨uttner distributionas f J ( x, p ) = n ( x ) σ ( x ) C π exp (cid:18) − p µ u µ ( x ) T ( x ) (cid:19) , (5)where C π = 4 πm π T K ( m π /T ), at temperature T , and K is a modified Bessel function.We assume that the single particle distributions, f ( x, p ), in the source functions are J¨uttner distributions,which depend on the local velocity, u µ ( x ), and we usethe notation u = u ( x ) = u µ ( x ).By using the Cooper-Frye (CF) freeze out descriptionwe can connect the Source function, S , to the phase spacedistribution function on the freeze out hypersurface. Letthe space-time points of the hyper-surface be given inparametric form x F O = x F O ( x ), which can be givenby the freeze out condition (e.g t =const., τ =const., T =const. or other). In the source function formalismthis corresponds to a 4-volume integral (cid:90) d x S ( x, p ) = (cid:90) d x f J ( x, p ) P ( x, p ) = (cid:90) d x f J ( x, p ) δ ( x − x F O ) p µ ˆ σ µ , where the emission probability is [11] P ( x, p ) = δ ( x − x F O ) p µ ˆ σ µ . This CF freeze out treatment is the mostfrequent in fluid dynamical models. This sudden freezeout assumption can be relaxed by assuming an extendedfreeze out layer in the space time via replacing theDirac delta function with a freeze out profile functionin P ( x, p ), e.g.: P ( x, p ) = δ ( x − x F O ) p µ ˆ σ µ −→ √ ∆ π exp (cid:18) − ( s − s F O ) ∆ (cid:19) p µ ˆ σ µ , where s is a local coordinate in the direction of ˆ σ µ , andthe local width of the freeze out layer is ∆ = ∆( x ) (whichshould tend to zero to get the Dirac delta function forthe emission probability). This description is then com-pletely general, with the only assumption that the emis-sion probability has a Gaussian profile. (This could alsobe relaxed.)If we assume that the two coincident particles origi-nate from two points, x and x , the expression of the At the latest times presented here, t = 3 .
56 fm/c, ( ∼ correlation function, Eq. (3) will be become [10] R ( k, q ) = (cid:90) d x d x S ( x , k ) S ( x , k ) cos[ q ( x − x )] × exp (cid:20) − q · (cid:18) u ( x ) T ( x ) − u ( x ) T ( x ) (cid:19)(cid:21) , (6)and the corresponding J ( k, q )-function is J ( k, q ) = (cid:90) d x S ( x, k ) exp (cid:20) − q · u ( x )2 T ( x ) (cid:21) exp( iqx ) , (7)In Ref. [10] different one, two and four source sys-tems were tested with and without rotation. Here westudy only the case where the emission is asymmetric and dominated by the fluid elements facing the detector.In numerical fluid dynamical studies of symmetric(A+A) nuclear collision the initial state is symmetricaround the center of mass (c.m.) of the system, and (ifwe do not consider random fluctuations) this symmetryis preserved during the fluid dynamical evolution.Let us consider the usual conventions, z is the beamaxis, and the positive z -direction is the direction of theprojectile beam. The impact parameter vector pointsinto the positive x -direction, i.e. towards the projectile.Finally the y -axis is orthogonal to both.The fluid dynamical system, without fluctuations canbe considered as a set of symmetric pairs of fluid cells.The emission probabilities from the two fluid cells of asource pair are not equal.If we have several fluid cell sources, s , with Gaussianspace-time (ST) emission profiles, then the source func-tion in J¨uttner approximation is (cid:90) d x S ( x, k ) = (cid:88) s (cid:90) s d x s dt s S ( x s , k ) =(2 πR ) / (cid:88) s γ s n s ( x ) ( k µ ˆ σ µs ) C s exp (cid:20) − k · u s T s (cid:21) , (8)where n s = n π , and the spatial integral over a cell volumeis, V cell = (2 πR ) / while the time integral is normalizedto unity. Similarly the J -function is J ( k, q ) = (cid:88) s e − q · usTs e iqx s (cid:90) s d x S s ( x, k ) e iqx . (9)We then assume that the FO layer is relatively narrowcompared to the spatial spread of the fluid cells, so thatthe peak emission times, t s , of all fluid cells are thesame. Then the exp( iq t s ) factor drops out from the J ( k, q ) J ( k, − q ) product. This FO simplification is jus-tified for rapid and simultaneous hadronization and FO If the emission is happening through a layer with time-like nor-mal, but the peak is not at constant t s , but rather at constant τ s , then we can adapt the coordinate system accordingly, i.e. wecan use the τ, η coordinates instead of t, z , see e.g. [11]. from the plasma. For dilute and transparent matter thecorrelations from the time dependence of FO should behandled the same way as the spatial dependence.Due to mirror symmetry with respect to the [ x, z ], re-action plane, it is sufficient to describe the cells on thepositive side of the y -axis. The other side is the mirrorimage of the positive side. Then we can evaluate thecorrelation function the same way as in Ref. [10].Thus we define the quantities: Q c = (cid:0) πR (cid:1) / exp (cid:20) − R q (cid:21) ,P s = γ s n s C s exp (cid:20) − k u s T s (cid:21) ,Q ( q ) s = exp (cid:20) − q u s T s (cid:21) ,w s = ( k µ ˆ σ µs ) exp (cid:20) − Θ s q (cid:21) , (10)where u s = γ s , the local 4-direction normal of the meanparticle emission from an ST point of the flow is ˆ σ µs (as-sumed to be time-like), R is the size (radius) of the fluidcells, and Θ s is the path length of the time integral fromthe ST point of the source, s , while assuming a Gaussianemission time profile [10]. The weights, w s arise directlyfrom the Cooper-Frye formula [11].We can reassign the summation for pairs, so that s = { i, j, k } will correspond to a pair of cells: at { i, j, k } andits reflected pair across the c.m. point at the same timeat { i ∗ , j ∗ , k ∗ } . Then the function S ( k, q ) becomes (cid:90) d xS ( x, k ) = (cid:0) πR (cid:1) / × (cid:88) s P s (cid:20) w s exp (cid:18) ku s T s (cid:19) + w ∗ s exp (cid:18) ku ∗ s T s (cid:19)(cid:21) , (11)while, the function J ( k, q ) becomes J ( k, q ) = Q c (cid:88) s P s (cid:20) Q ( q ) s w s exp (cid:20)(cid:16) k + q (cid:17) u s T s (cid:21) e i qx s + Q ( q ) s w ∗ s exp (cid:20)(cid:16) k + q (cid:17) u ∗ s T s (cid:21) e i qx ∗ s (cid:21) Only the mirror symmetry across the participant c.m.is assumed, which is always true for globally symmetric,A+A, heavy ion collisions in a non-fluctuating fluid dy-namical model calculation. Then the correlation functioncan be evaluated using Eqs. (2-4).By using few fluid cell sources for tests, in Ref. [10]it was shown that in case of a globally symmetric fluiddynamical configuration the correlation function only in-cludes cos( c ku s ) and cosh( c ku s ) terms, therefore it willnot depend on the direction of the velocity, only on itsmagnitude. The direction dependence becomes apparentin the correlation function only if we take into accountthat due to the radial expansion and the opacity of the strongly interacting QGP, the emission probability fromthe far side of the system is reduced compared to the sideof the system facing the detector.Based on the few source model results the DifferentialHBT method was introduced by evaluating the differenceof two correlation functions measured at two symmet-ric angles, forward and backward shifted in the reactionplane in the participant c.m. frame by the same angle,i.e. at η = ± const., so that the Differential CorrelationFunction (DCF) becomes∆ C ( k, q ) ≡ C ( k + , q out ) − C ( k − , q out ) . (12)For the exactly ± x -symmetric spatial configurations (i.e. k + x = k − x and k + z = − k − z ), e.g. central collisions orspherical expansion, ∆ C ( k, q ) would vanish! It wouldbecome finite if the rotation introduces an asymmetry. III. THE FREEZE-OUT
The HBT method is sensitive to the time developmentof particle emission, and well suited to transport mod-els where emission happens during the ST history of thecollision, although the emission is concentrated at a FOlayer. The fluid dynamical model is also able to describe w i t h o u t r o t a t i o n k = 0 . 2 / f m k = 5 / f m w i t h r o t a t i o n k = 0 . 2 / f m k = 5 / f m
C(k,q) q ( 1 / f m ) t = 3 . 5 6 f m / c
FIG. 1. (color online) The dependence of the standard corre-lation function in the k + direction from the collective flow, atthe final time, 3.56 fm/c after reaching local equilibrium and8.06 fm/c from the first touch, including the initial longitudi-nal expansion Yang-Mills field dynamics [19]. emission in a ST volume or layer [12, 13], or in hybridmodels where a microscopic transport module is attachedto the fluid dynamics, e.g. [14]. The determination of theFO surface normal or the emission direction from the STFO layer and the emission profile in this layer are thesubjects of present theoretical research, see [7, 15–18].This complex FO process certainly has an influence onthe HBT analysis, but our present aim is not to repro-duce exactly a given experiment.We focus on a single collective effect, the rotation, de-veloping from the angular momentum during the initialstages of the fluid dynamics. Thus we constrain our dis-cussion to the fluid dynamical stage, and adopt a rela-tively simple FO description from ref. [11], which canbe implemented in Eq. (10). This provides the emissionweight factors, w s , which depend on the local mean emis-sion direction ˆ σ µ , and the flow velocity at the emissionpoint.The detector configuration is given by the two particlesreaching a given detector in the direction of k . Thus theemission weights depend on the normal of the emissionsurface and of the emission, i.e. ˆ σ and ˆ k . Most of theparticles FO in a layer along a constant proper time hy-perbola, with a dominant flow 4-velocity normal to thishyperbola: ˆ σ µ ≈ u µ . The origin of the hyperbola is at aST point, which may precede the impact of the Lorentzcontracted nuclei [15].We assume in the actual numerical calculations that inthe expression of the weight, in Eq. (10), is the same forall surface layer elements: Q ( q ) s = Q ( q ) and Θ s = Θ, sothat w s = ( k µ ˆ σ µs ) exp( − Θ q / , where ˆ σ sµ = ( σ s , σ s ),so that k µ ˆ σ µs = k σ s + kσ s . If the emission path time-length, Θ, tends to zero, then the time modifying factorbecomes unity. With the choice ˆ σ µ = u µ , the time-likeFO normal is ˆ σ sµ = ( γ s , u s ). Then ( k µ ˆ σ µs ) = γ s k + ku s .So the weight becomes w s = ( γ s k + ku s ) exp( − Θ q / . (13)This weight is explicitly different for the mirror im-age cell at x ∗ s → − x s , where u ∗ s → − u s and then w ∗ s = ( γ s k − ku s ) exp( − Θ q / . The weight factors appear both in the nominator anddenominator of the correlator, so its normalization is bal-anced. On the other hand the role of the different factorsin the weight have an effect to determine, which cells con-tribute more, which cells contribute less to the integratedresult. w i t h o u t r o t a t i o n k = 0 . 2 / f m k = 5 / f m D C(k,q) q ( 1 / f m ) w i t h r o t a t i o n k = 0 . 2 / f m k = 5 / f mt = 3 . 5 6 f m / c
FIG. 2. (color online) The Differential Correlation Function,∆ C ( k, q ), at the final time with and without rotation. IV. RESULTS
The sensitivity of the standard correlation function onthe fluid cell velocities decreases with decreasing dis-tances among the cells. So, with a large number ofdensely placed fluid cells where all fluid cells contributeequally to the correlation function, the sensitivity on theflow velocity becomes negligibly weak.Thus, the emission probability from different ST re-gions of the system is essential in the evaluation. Thisemission asymmetry due to the local flow velocity occursalso when the FO surface or layer is isochronous or if ithappens at constant proper time.We studied the fluid dynamical patterns of the calcula-tions published in Ref. [2], where the appearance of theKHI is discussed under different conditions. We chosethe configuration, where both the rotation [1], and theKHI occurred , at b = 0 . b max with high cell resolutionand low numerical viscosity at LHC energies, where theangular momentum is large, L ≈ (cid:126) [7].In spatially symmetric few source configurations [10],the standard correlation function did not show any dif-ference if it is measured at two symmetric k and q -out angles, e.g. in the reaction, [x-z] plane at k + =( k x , , + k z ), q + = ( q x , , + q z ) and k − = ( k x , , − k z ), q − = ( q x , , + q − ), i.e. ∆ C ( k, q ) vanished. Here we havechosen two directions at η = ± .
76, that is at polar an-gles of 90 ±
40 degrees. These are measurable with theALICE TPC and at the ATLAS and CMS also.The standard correlation function is both influencedby the ST shape of the emitting source as well as itsvelocity distribution. The correlation function becomesnarrower in q with increasing time primarily due to therapid expansion of the system. At the initial configura-tion the increase of | k | leads to a small increase of thewidth of the correlation function.Nevertheless, in theoretical models we can switch offthe rotation component of the flow, and analyse how therotation influences the correlation function and especiallythe DCF, ∆ C ( k, q ).Fig. 1 compares the standard correlation functionswith and without the rotation component of the flow atthe final time moment. Here we see that the rotationleads to a small increase of the width in q for the dis-tribution at high values of | k | , while at low momentumthere is no visible difference.In Fig. 2 ∆ C ( k, q ) is shown for the configuration withand without rotation. For k = 5 / fm the rotation in-creases both the amplitude and the width of ∆ C . Thedependence on | k | is especially large at the final time.In the original K frame defined by the beam directionand the impact parameter, we can describe the vector k with coordinates, k = (cid:26) k x k z (cid:27) . In the K (cid:48) frame the samevector is then k (cid:48) ( α ) = (cid:26) k x (cid:48) k z (cid:48) (cid:27) = (cid:26) k x cos α − k z sin αk z cos α + k x sin α (cid:27) . (14) ΔC=0, by Def.Symmetry axis α = 0non-tiltednon-rotating (a) k x k z ΔC ≠ 0non-tiltedrotating (b) k z k x α Δ C≠ 0 Δ C α = 0 (c) k x k z k x ’ k z ’ rotatingtilted Symmetry axis α ≠ 0 tiltedrotating (d) α k z k x k z ’ k x ’ Δ C≠ 0 Δ C α ≠ 0 FIG. 3. (Color online) Sketch of the configuration in differentreference frames, with and without rotation of the flow. Thenon-rotating configurations have may have radial flow velocitycomponents only. The DCF, ∆ C α ( k, q ) is evaluated in a K (cid:48) reference frame rotated by and angle α in the x, z , reactionplane. We search for the angle α , where the non-rotataingconfiguration is ”symmetric”, so that it has a ”minimal” DCFas shown in Fig. 4. Then one can define the DCF,∆ C α ( k (cid:48) , q (cid:48) ) , (15)which depends on the angle α . We have to find the propersymmetry axes of the emission ellipsoid. The conven-tional way would be the standard azimuthal HBT, how-ever, we can exploit the features of the DCF. As theanalytic examples [10] show if (i) the shape is symmetricaround the x (cid:48) axis, and (ii) there is no rotation in theflow, then ∆ C α (cid:0) k (cid:48) , q (cid:48) (cid:1) = 0 . (16)Thus we can use the DFC to find the angle, α (cid:48) , of theproper frame K (cid:48) also. For a given | k | (e.g. | k | = 5 / fm),we search for the minimum of the norm of the DCF as afunction of α .This procedure is performed and the result is shown inFig. 4 We separated the effect of the rotation by findingthe symmetry angle where the rotation-less configurationyields vanishing or minimal ∆ C for a given transversemomentum k .The figure shows the result where the rotation com-ponent of the velocity field is removed. The DCF showsa minimum in its integrated value over q , for α = − α . Unfortunately this is not possi-ble experimentally, so the direction of the symmetry axes should be found with other methods, like global flow anal-ysis and/or azimuthal HBT analysis. To study the de-pendence on the angular momentum the same study wasfor lower angular momentum also, i.e. for a lower (RHIC)energy Au+Au collisions at the same impact parameterand time. We identified the angle where the rotation-lessDCF was minimal, which was α = − D C a (k,q) q ( 1 / f m ) a ( d e g r e e ) - 2 0 - 1 1 . 5 - 1 1 - 1 0 0 1 0 2 0w i t h o u t r o t a t i o nk = 5 / f mt = 3 . 5 6 f m / c FIG. 4. (Color online) The DCF at average pion wavenum-ber, k = 5 / fm and fluid dynamical evolution time, t =3 . q (in units of 1/fm). The DCFis evaluated in a frame rotated in the reaction plane, in thec.m. system by angle α . We did this for two different energies, Pb+Pb /Au+Au at √ s NN = 2 . / . α = − / − respectively. In these deflected frames we evalu-ated ∆ C for the original, rotating configurations, whichare shown in Fig. 5. This provides an excellent measureof the rotation.On the other hand the rotation-less configuration can-not be generated from experimental data in an easy way.Other methods like the Global Flow Tensor analysis, orthe azimuthal HBT analysis [20] can provide an estimatefor finding the deflection angle α . V. CONCLUSION
The analysed model calculations show that the Differ-ential HBT analysis can give a good quantitative measureof the rotation in the reaction plane of a heavy ion col-lision. These studies show that this measure is propor-tional to the beam energy or total angular momentum The negative angles are arising from the fact that our modelcalculations predict rotation, with a peak rotated forward [1]. k = 5 / f m t = 3 . 5 6 f m / c D C(k,q) q ( 1 / f m )
P b + P b @ 2 . 7 6 T e V a = - 1 1 d e g r e e s w i t h r o t a t i o n w i t h o u t r o t a t i o n A u + A u @ 2 0 0 G e V a = - 8 d e g r e e s w i t h r o t a t i o n w i t h o u t r o t a t i o n FIG. 5. (color online) The DCF with and without rotationin the reference frames, deflected by the angle α , where therotation-less DCF is vanishing or minimal. In this frame theDCF of the original, rotating configuration indicates the effectof the rotation only. The amplitude of the DCF of the originalrotating configuration doubles for the higher energy (higherangular momentum) collision. (while the polarization [4] is not). It shows the best sen-sitivity at higher collective transverse momenta.To perform the analysis in the rotation-less symmetryframe one can find the symmetry axis the best with theazimuthal HBT method, which provides even the trans-verse momentum dependence of this axis [20].It is also important to determine the precise Event byEvent c.m. position of the participants [21], and min-imize the effect of fluctuations to be able to measureaccurately the emission angles, which are crucial in thepresent ∆ C ( k, q ) studies. ACKNOWLEDGEMENTS
This work was supported in part by the HelmholtzInternational Center for FAIR. We thank F. Becattini,M. Bleicher, G. Graef, P. Huovinen and J. Manninen forcomments. [1] L.P. Csernai, V.K. Magas, H. St¨ocker, and D.D.Strottman, Phys. Rev. C , 024914 (2011).[2] L.P. Csernai, D.D. Strottman and C. Anderlik, Phys.Rev. C , 054901 (2012).[3] L.P. Csernai, V.K. Magas, and D.J. Wang, Phys. Rev. C87, 034906 (2013).[4] F. Becattini, L.P. Csernai, D.J. Wang, Phys. Rev. C ,034905 (2013).[5] Sz. Horv´at, V.K. Magas, D.D. Strottman, L.P. Csernai,Phys. Lett. B , 277 (2010).[6] B. McInnes, and E.Teo, Nucl. Phys. B , 186 (2014).[7] V. Vovchenko, D. Anchishkin, and L.P. Csernai, Phys.Rev. C , 014901 (2013).[8] W. Florkowski: Phenomenology of Ultra-relativisticheavy-Ion Collisions , World Scientific Publishing Co.,Singapore (2010).[9] A.N. Makhlin, Yu.M. Sinyukov, Z. Phys. C , 69 (1988);S.V. Akkelin, Yu.M. Sinyukov, Z. Phys. C , 501 (1996);T. Cs¨org˝o, S.V. Akkelin, Y. Hama, B. Luk´acs, and Yu.M.Sinyukov, Phys. Rev. C , 034904 (2003) .[10] L.P. Csernai, S. Velle, (2013) arXiv:1305.0385[11] T. Cs¨org˝o, Heavy Ion Phys. , 1-80, (2002); arXiv:hep-ph/0001233v3 [12] E. Moln´ar, L. P. Csernai, V. K. Magas, Zs. I. Lazar, A.Nyiri, and K. Tamosiunas, J. Phys. G , 1901 (2007).[13] E. Moln´ar, L. P. Csernai, V. K. Magas, A. Nyiri, and K.Tamosiunas, Phys. Rev. C , 024907 (2006).[14] Yu-Liang Yan, Yun Cheng, Dai-Mei Zhou, Bao-GuoDong, Xu Cai, Ben-Hao Sa, and Laszlo P Csernai, J.Phys. G , 025102 (2013).[15] D. Anchishkin, V. Vovchenko, and L.P. Csernai, Phys.Rev. C , 014906 (2013).[16] Y. Cheng, L.P. Csernai, V.K. Magas, B.R. Schlei, and D.Strottman, Phys. Rev. C , 064910 (2010).[17] P. Huovinen, H. Petersen, Eur. Phys. J. A , 171 (2012)[18] D. Anchishkin, V. Vovchenko, and S. Yezhov, Int. J.Modern Phys. E , 014901 (2001), and Nucl. Phys. A , 167-204(2002).[20] M.A. Lisa, N.N. Ajitanand, J.M. Alexander, et al., Phys.Lett. B , 1 (2000); M.A. Lisa, U. Heinz, U.A. Wiede-mann, Phys. Lett. B , 287 (2000); E. Mount, G.Graef, M. Mitrovski, M. Bleicher, M.A. Lisa, Phys. Rev.C , 014908 (2011); G. Graef, M. Bleicher, and M. Lisa,Phys. Rev. C , 014903 (2014).[21] L.P. Csernai, G. Eyyubova, V.K. Magas, Phys. Rev. C86