Contents lists available at ScienceDirect

Physics Letters B

www.elsevier.com/locate/physletb

Work distribution of an expanding gas and transverse energy production in relativistic heavy ion collisions

Bin Zhang *, Jay P. Mayfield

Department of Chemistry and Physics, Arkansas State University, P.O. Box 419, State University, AR 72467-0419, USA

CrossMark

ARTICLE INFO

ABSTRACT

Article history:

Received 13 September 2013 Received in revised form 21 January 2014 Accepted 4 March 2014 Available online 12 March 2014 Editor: W. Haxton

Keywords: Work distribution Transverse energy Relativistic heavy ion collisions

The work distribution of an expanding extreme relativistic gas is shown to be a gamma distribution with a different shape parameter as compared with its non-relativistic counterpart. This implies that the shape of the transverse energy distribution in relativistic heavy ion collisions depends on the particle contents during the evolution of the hot and dense matter. Therefore, transverse energy fluctuations provide additional insights into the Quark-Gluon Plasma produced in these collisions.

© 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/3.0/). Funded by SCOAP3.

1. Introduction

Transverse energy is an important characteristic of relativistic heavy ion collisions aiming at creating conditions similar to those existed in the early Universe [1-6]. The hot and dense matter produced immediately after the collision of a nucleus going in the longitudinal direction and the other in the opposite direction can be considered as being composed of transverse fluid slices undergoing longitudinal expansion following the two receding nuclei [7]. In the nucleon-nucleon center-of-mass frame, a slice that is closer to a nucleus has a bigger longitudinal speed. Transverse energy is longitudinally boost-invariant and therefore directly reflects the condition of the local rest frame (slice) irrespective of its longitudinal flow speed. It is sensitive to the longitudinal work between adjacent slices and thus carries information about the evolution of the hot and dense matter produced in heavy ion collisions [8-10].

Experimentally, the transverse energy distribution is approximately a gamma distribution [11]. Interestingly, the work distribution for the adiabatic compression or expansion of a dilute and interacting classical gas has been calculated by Crooks and Jarzyn-ski [12], and it is also a gamma distribution. The analogy between a longitudinally expanding Quark-Gluon Plasma and the adiabatic expansion of a classical gas in a cylinder prompts us to look at the latter more carefully. Crooks and Jarzynski's calculation is for a non-relativistic gas. Relativistic effects can be important for quarks and gluons in a Quark-Gluon Plasma. In the following, we will

Corresponding author.

E-mail address: bzhang@astate.edu (B. Zhang).

show that for an extreme relativistic gas, the work distribution is also a gamma distribution. But the extreme relativistic work distribution has a different shape parameter relative to the non-relativistic one. The transverse energy distributions can also be calculated, and they are gamma distributions similar to the work distributions. Because of the simplicity of the model, the two parameters of the transverse energy distribution have clear physical meanings. They are shown to reflect important properties of the evolution of the system.

2. Work distribution and transverse energy production

The work distribution for a non-relativistic, dilute, interacting, classical gas undergoing adiabatic compression or expansion is given by [12]

p(W ) = ^(ßWW

r | œ | V, a

-ßW )0(aW).

Here W is the work on the system. It is positive for compression and negative for expansion. In three dimensions, a is related to the initial volume V0 and final volume V\ by

a =1 - 1.

a is positive for compression and negative for expansion. The unit step function 6(-) in Eq. (1) ensures that W and a always have the same sign. /) is the inverse of the initial fundamental temperature. In three dimensions, k is related to the total number of particles N

http://dx.doi.org/10.1016/j.physletb.2014.03.007

0370-2693/© 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/). Funded by SCOAP3.

via k = 3N/2. Unless stated otherwise, we will use the natural unit system, in which the reduced Planck constant h, the speed of light in vacuum c, and the Boltzmann constant kB are set to 1.

The distribution of the magnitude of work then acquires the form

B i №1\k-1

P(|W |) =

\a\T(k) \ |a| )

, ß W |\

exp| - -w)

It is a gamma distribution described by shape k = 3N/2 and scale s = \a\/B.

In the following, we will derive the work distribution for the extreme relativistic case and make some comparisons. A good starting point is the number of energy states with the energy of the gas less than E. It can be described by the asymptotic formula [13]

$(E; V) =

1 VN (8n)NE

(2n)3N N! (3N)! n2NN! (3N)!' The density of states can now be calculated as

g(E- V) = — =

dE n2NN! (3N - 1)!'

This leads to the partition function Z(P, V) = J dEg(E; V) exp(-BE) =

n2NN! P3N' Energy E follows the canonical distribution

g(E' V) (BE)3N-1

P (E; B) = SL^I exp(-B E) = ,, exp(-B E).

Z(B, V) (3N - 1)!

The work during an adiabatic process is the change in the internal energy, i.e.,

W = E1 - E0.

Assuming ergodicity, $(E; V) is an adiabatic invariant, and E1 can be related to E0 via

E1 = ( tT E 0. (9)

Therefore,

W ^(V0) / - ^E0 = aE0. (10)

r- (")

is different from the non-relativistic formula given in Eq. (2). Now the work distribution is given by

p(W) = j dE0 P(E0; B)S(W - aE0) k-1

_ ß f-wy

|a|r(kA a J

ex^ - —W )°(aW),

where k = 3N and a is given by Eq. (11). The distribution of the magnitude of work is also given by Eq. (3) with k and a given by the extreme relativistic formulas above. Therefore, the extreme relativistic case has the same work distribution compared to the non-relativistic case, but it has different formulas for the parameters. In particular, the shape parameter changes from k = 3N/2 in the non-relativistic case to k = 3N in the extreme relativistic case.

0.15 -

#=4,^^=16 dashed: compression solid: expansion thin: non. rel. thick: ex. rel.

0.05 -

Fig. 1. Work magnitude distributions for the compression and expansion of a non-relativistic (non. rel.) gas and an extreme relativistic (ex. rel.) gas.

Making use of the expression for the free energy

F(ß, V) = —— lnZ—, V) = --ln. ,„,

) ß ) ß 2NN! -3N

the work distribution can be shown to satisfy the Jarzynski equality [14,15]

- ln(exp(-B W)} = -lnj dW p(W) exp(-B W)

= Nln( V0 )= ß^F,

where AF = F(B, V1) - F(B, V0) is the change in the free energy between two states with the same B and different volumes.

It is also straightforward to show that the work distribution for the forward process pF and that for the corresponding reverse process pR satisfy the Crooks fluctuation theorem [16,17]

PF (W) PR (-W)

/ v \N

i^j exp(BW) = exp(B(W - AF)). (15)

Now let us compare the non-relativistic case and the extreme relativistic case. Fig. 1 gives the work magnitude distributions. In both the non-relativistic case and the extreme relativistic case, the compression curve and the expansion curve meet at B \AF\ = 4ln(16). The extreme relativistic curves are closer to 4ln(16) than the corresponding non-relativistic curves. Consequently, the extreme relativistic case has the higher probability of having W < AF. However, the second law of thermodynamics, i.e., the average work (W) is not smaller than AF, is still valid [18]. From the fact that the average of a gamma distribution is the product of the shape and scale parameters, for the compression case, the non-relativistic average work (W)n = 3N/(2B)((V0/V 1)2/3 - 1), and the extreme relativistic average work (W)e = 3N/B((V0/V 1)1/3 - 1). This can be compared with A F = N/B ln(V0/V1). We arrive at the relation (W)n > (W)e > AF > 0. It shows the ordering of (W)n and (W)e, and that both the non-relativistic and extreme relativis-tic cases satisfy the second law of thermodynamics. Likewise, for the expansion case, 0 > (W)n > (W)e > AF.

With the same compression or expansion ratio V0/V1, the non-relativistic and the extreme relativistic cases have different a values as given by Eqs. (2) and (11). Thus the work magnitude distributions have different scale parameters s = \a\/B. During an adiabatic process, the ensemble remains canonical. Therefore, for

Fig. 2. Work magnitude distributions for a non-relativistic gas and an extreme rela-tivistic gas undergoing the same temperature change.

the non-relativistic case with initial temperature T0 and final temperature T1,

and for the extreme relativistic case,

T ' = < £f T0-

This leads to an interesting expression for the scale parameter, i.e., s = |T 1 — T0| = |AT|. This can be more relevant for relativis-tic heavy ion collisions, where the hot and dense matter can be considered as starting at some initial temperature and stopping at some freeze-out temperature. Fig. 2 shows that non-relativistic and extreme relativistic gases have very different work distributions. Since they have the same scale parameter, the difference in the means comes from different shape parameters, and the extreme relativistic one is about twice that of the non-relativistic case.

Before calculating the transverse energy distribution, we will look at the final energy distribution. The final energy can be related to the initial energy by E1 = (V0/V 1){2/3>1/3} = qE0. Unless stated otherwise, the first choice in the braces is for a non-relativistic gas, and the second is for the extreme relativistic one. Now the final energy distribution

p(E 1) = dE0 P(E0; P)S(E 1 - qE0)

= ; a = p(e 1; p-q \q J V q

is a gamma distribution with shape k = {3N/2, 3N} and scale s = q/fi = T0(V0/V 1)<2/3,1/3} = T1. This is expected as the ensemble remains canonical during the adiabatic process.

In order to get the transverse energy distribution, we need to approximate the sum over particles by an integral. For the non-relativistic case,

E1 EU i=1

= f d3p C exp T-r^) p-J \ 2mT1 / 2m

= C I dp p2 exp (--P— 1 — I d cos 9 1 dé

' 2mTi)2mJ J é

= D J d cos 9 = 2 D. 1

In the above, E1ti is the final energy of particle i. C and D are constants. 9 and é are the polar and azimuthal angles in the spherical coordinate system where the polar axis goes along the longitudinal direction. The final transverse energy

EE1X.i = /

En = > .En; = d3pC exp

2mT1 2m

= CI dp p2 exp(--p—)— I sin2 9 d cos 9 I dé

J 2mT1)2mJ J é

= Dj sin2 9 d cos 9 = 3 D. 1

Therefore, E 1x = f E1.

For the extreme relativistic case,

E1 = ^E1,i = J d3 pC ex^ - ^ )p 1

= D' j d cos 9 = 2 D' -1

where C' and D' are constants.

I>x,i = , i=1

-1±,i = I d3pC'exp( -^ )p±

D j sin 9 d cos 9 = — D'. -1

Hence, E1x = §■ E1.

The non-relativistic and the extreme relativistic cases can be summarized into one formula E 1x = aE1, where a = {2/3, n/4}. Now the final transverse energy distribution

d(E 1±) = I dE1 p(E 1)S(En- aE1) = 1

= ipfe = p(e 1±;£

aq aq aq

is a gamma distribution with shape k = {3N/2, 3N} and scale s = (aq)/fi ={2/3, n/4}T1. This tells us that the shape parameter is very sensitive to the particle contents during the evolution while the scale parameter is slightly sensitive and is mainly determined by the final temperature.

3. Summary and discussions

The work distribution for an extreme relativistic gas undergoing an adiabatic process is shown to be a gamma distribution with a shape parameter twice as large as that for the non-relativistic gas. Both cases have a scale parameter that can be related to

the change in the system temperature. The corresponding transverse energy distributions are also gamma distributions. In both the extreme relativistic and the non-relativistic cases, the shape parameter is the same as that for the work distribution, and the scale parameter is related to the final temperature. This gives insights into relativistic heavy ion collisions where transverse energy distributions are approximate gamma distributions. In particular, a change in the particle mass can lead to a change in the shape parameter, and the scale parameter is mainly determined by the freeze-out temperature and depends slightly on the particle mass. This may lead to deviations from the superposition of nucleon-nucleon collisions as has been calculated in [11]. On the other hand, the transverse energy distribution in proton-proton collisions is also approximately a gamma distribution [11]. Thus the parameters may reflect details of proton-proton collisions.

In general, the density of states cannot be written as a power of the energy. In particular, Appendix A shows that

g(E; V) =

(4n m3)

(2n)3NN! (N - 1)!m

--+ N - V^ cosh tk

ЬП/dtj cosh(2tj)

x 0|--+ N cosh fyj,

t* = ln - + 1 +JI - + 1 - 1

The energy distribution can be calculated numerically by making use of Eq. (24), but it is not a simple gamma distribution, and it is unlikely that the work distribution can be expressed as a gamma distribution. It is possible that the work distribution can be approximated by a gamma distribution interpolating between the non-relativistic case and the extreme relativistic case. If so, qualitative changes due to the change of particle mass are also expected to show up strongly in the effective shape parameter.

In addition to the finite mass corrections, other factors come into play in relativistic heavy ion collisions. The longitudinally expanding gas with local thermal equilibrium can only give some general guidance for the effects of final state interactions. Even for the extreme relativistic case, as the kinetic energy is much larger than the rest mass energy, particle production and annihilation can happen, and entropy is expected to change accordingly. At this moment, the role of particle number changing processes is still under intense investigation [19-26]. As particle production during the expansion process leads to more cooling, the entropy increase is expected to be small [27]. It is not clear whether and by how much particle number changing processes affect the proportion relationship between the initial and final energies. If the relation is significantly modified, there could be large deviations from the gamma distribution for the transverse energy distribution. Other than the particle number changing processes, the changing particle contents, the transverse expansion, and the differential freeze-out all contribute to the evolution of the transverse energy distribution. Studies with dynamical models will be necessary to sort out the details.

The above derivation of the transverse energy distribution depends on local thermal equilibrium. Deviations from equilibrium can lead to deviations from the gamma distribution. For example,

initial conditions based on the Glauber model [28] or various saturation models [29,30] can be very different from thermal initial conditions. The difference may lead to observable deviations from the gamma distribution when higher order moments are studied. In other words, longitudinal flow results can complement the widely investigated transverse flow analyses [31]. During the late stage of the evolution, the viscosity is large, and the system cannot maintain equilibrium. The longitudinal work will be significantly reduced. At the other extreme from equilibrium, no longitudinal work is expected for the free streaming case, and the transverse energy distribution does not change. Therefore, as a first approximation, the change of the transverse energy distribution due to the late stage can be neglected. To be a little more precise, if the average kinetic energy can be used as a measure of the temperature in the non-equilibrium case, the freeze-out temperature that goes into the transverse energy distribution estimate should be higher than the non-equilibrium "temperature".

The gamma distribution is also limited to a classical gas. Recently, there are some discussions of the possibility of forming a Bose-Einstein condensate in the early stage of a relativistic heavy ion collision [32,33]. If a Bose-Einstein condensate is formed, the transverse energy distribution may have some noticeable difference from a gamma distribution. One can examine the key elements in the derivation of the gamma distribution and see the difference. One thing that does not change in the derivation is the relation between the work during an adiabatic process and the initial energy. The reason is that the relation between the energy of a single particle in a box and the volume is independent of whether the particle is a Boson or a Fermion. If interactions are negligible, the relation between the total energy of the system and the volume is independent of whether the particles are Bosons or Fermions. However, the number of states with energy smaller than a given value depends on the number of ways to partition an integer [34,35]. Thus the density of states depends on whether the system is composed of Bosons or Fermions. Only in the non-degenerate limit does the work distribution become a gamma distribution.

It is interesting to see what experimental data can teach us. Our focus will be on relativistic heavy ion collisions where hydrodynamics is successful in describing various experimental observations. In particular, preliminary data from the PHENIX Collaboration at the Relativistic Heavy Ion Collider [36] will be used. The transverse energy here is the transverse energy measured in the electromagnetic calorimeter. If г is the average transverse energy, and a2 is the variance, they can be related to the shape parameter k and the scale parameter s via a/г = 1/Vfc and a2/г = s. For Au + Au collisions at the nucleon-nucleon center-of-mass frame energy */sNN = 200 GeV, the 0-5% centrality bin has a/г = 0.116 and a2/г = 1.60 GeV. Therefore, the shape parameter k = 74, and the scale parameter s = 1.60 GeV. They appear to be outside the ranges expected from the commonly accepted initial particle number and final temperature values. However, each centrality bin has a distribution of particle numbers [37]. Selection according to the number of particles may improve on the situation.

The experimental data may provide more information. If using the average number of particles for a centrality bin only introduces a common rescaling factor for all centralities, and if the particle contents remain the same over different centralities, the ratio of a/г to 1/VN is expected to be independent of centrality. Assuming the charge particle pseudo-rapidity distribution dNch/dn is 2/3 of the particle pseudo-rapidity distribution, the centrality dependence of a/¡i/(1/t/dNch/dn) should reflect the particle contents. Fig. 3 shows this ratio for „JsNN = 62.4 GeV and 200 GeV. If other complications are not important, for the 200 GeV case, as dNch/dn increases toward 100, the particles responsible for the longitudi-

4 3.5 3 2.5

O û o

open: 62.4GeV solid: 200GeV

0 100 200 300 400 500 600 700 dNJdr,

Fig. 3. Centrality dependence of the ratio of the relative width o/i to 1/ydNch/dn for ysNN = 62.4 GeV and 200 GeV.

nal work become lighter (more relativistic). The particles become heavier (less relativistic) as dNch/dn increases further from 100. One should bear in mind that the error bars are on the order of 5-10%, and there can be other factors affecting the centrality dependence. It also appears that 200 GeV has heavier particles for many centrality bins compared with the 62.4 GeV case.

In a recent preprint [38], the PHENIX Collaboration compared the experimental transverse energy distributions to results from some models based on the superposition of nucleon-nucleon collisions. They demonstrated that the experimental data favor the number-of-constituent-quark-participant model. However, the superposition of nucleon-nucleon collisions cannot generate collective flow (longitudinal, radial, or elliptic). As longitudinal flow reduces the event transverse energy, even the number-of-constituent-quark-participant model can give very different transverse energy distributions when the hydrodynamic evolution of the hot and dense matter is taken into account. On the other hand, if no dynamical models can give satisfactory explanation of the measured transverse energy distributions, it will be a real challenge to reconcile different flow phenomena.

Acknowledgements

We thank G.E. Crooks, P. Danielewicz, C. Jarzynski, M. Murray for helpful discussions. This work was supported by the U.S. National Science Foundation under grant No. PHY-0970104.

Appendix A. Derivation of the density of states formula

The number of states with energy less than E is

$(E; V) =

(2n)3NN This leads to the density of states

-J d3Np ^ E - .

g(E; V) = — =

d E (2n)3NN

-J d3Np ^ E - .

The integral

f (E) = j d3Np ^E E;

can be simplified by looking at its Laplace transform

F(ß) = j dE exp(-ßE) j d3Np si E E 0

= J d3Np exp^-ß ^^E

= H f d3Pi exp(-ßEi), i=i J

where Re ß > 0. For each particle, J d3 pi exp(-ßEi)

= 4n j dpi pi exp^-ß^\jp2 + m2 - mY) 0

= exp(ßm)2n j dm} m}

fdy cosh y exp(-ßm± cosh y)

= exp(ßm)4n y dm±m}Ki(ßm±)

= exp (ßm) —— K 2 (ßm), ß

where Ki(x) and K2(x) are modified Bessel functions. Therefore,

(4n m2 \N

F(fi) = exp(Nfim){—~K2(fim)j .

Now f (E) can be calculated from F(fi) by using the inverse Laplace transform.

fi'+iTO

f (E) = ¿i / F (fi) exp (fi E) dfi

fi'—¡TO fi'+iTO

1 f (4nm2 \N

= 2n J exp(N ßm)y—ß— K2(ßm)\ exp(ß E ) dß,

ß'-iœ

where ß' > 0. The modified Bessel function can be expressed as

K2(ßm) = /dt cosh(2t)exp(-ßm cosht). (A.8)

This leads to

f (E) = (4nm2)N^ f dtj cosh(2tj) j=1 J

e'+i f

1 f exp(ß(E + Nm - m ¿J^ cosh tk))

2n J ßN

fi'-iœ

= (4n mifnfdtj j=i0

cosh(2tj)

(E + Nm - mjjt=\ cosh tk) (N -1)!

/ N \ 01 E + Nm - m^^coshtl I

(4nm3)N N

(N - 1)!m

mU J dtj cosh(2tj)

x I E + N - V^ cosh tk I m k=1 k \

x 0| — + N -X^ cosh tH, m = /

where the upper bound t* is the positive solution of

+ 1 - cosh t* = 0.

(A.10)

References

[1] K. Adcox, et al., PHENIX Collaboration, Phys. Rev. Lett. 87 (2001) 052301, arXiv:nucl-ex/0104015.

[2] J. Adams, et al., STAR Collaboration, Phys. Rev. C 70 (2004) 054907, arXiv: nucl-ex/0407003.

[3] K. Adcox, et al., PHENIX Collaboration, Nucl. Phys. A 757 (2005) 184, arXiv: nucl-ex/0410003.

[4] J. Adams, et al., STAR Collaboration, Nucl. Phys. A 757 (2005) 102, arXiv: nucl-ex/0501009.

[5] A. Toia, J. Phys. G 38 (2011) 124007, arXiv:1107.1973 [nucl-ex].

[6] M. Malek, CMS Collaboration, Nucl. Phys. A 904-905 (2013) 385c.

[7] J.D. Bjorken, Phys. Rev. D 27 (1983) 140.

[8] P.V. Ruuskanen, Phys. Lett. B 147 (1984) 465.

[9] K.J. Eskola, M. Gyulassy, Phys. Rev. C 47 (1993) 2329. [10] M. Gyulassy, Y. Pang, B. Zhang, Nucl. Phys. A 626 (1997) 999, arXiv:nucl-th/ 9709025.

11] A.K. Chaudhuri, Phys. Rev. C 47 (1993) 2875.

12] G.E. Crooks, C. Jarzynski, Phys. Rev. E 75 (2007) 021116.

13] R.K. Pathria, Statistical Mechanics, second ed., Elsevier ButterworthHeinemann, Oxford, U.K., 1996, p. 41.

14] C. Jarzynski, Phys. Rev. Lett. 78 (1997) 2690.

15] C. Jarzynski, Phys. Rev. E 56 (1997) 5018.

16] G.E. Crooks, Phys. Rev. E 60 (1999) 2721.

17] G.E. Crooks, Phys. Rev. E 61 (2000) 2361.

18] C. Jarzynski, Annu. Rev. Condens. Matter Phys. 2 (2011) 329.

19] Z. Xu, C. Greiner, Phys. Rev. Lett. 100 (2008) 172301, arXiv:0710.5719 [nucl-th].

20] J.-W. Chen, H. Dong, K. Ohnishi, Q. Wang, Phys. Lett. B 685 (2010) 277, arXiv: 0907.2486 [nucl-th].

21] B. Zhang, W.A. Wortman, Phys. Lett. B 693 (2010) 24, arXiv:1006.0270 [nucl-th].

22] J.-W. Chen, J. Deng, H. Dong, Q. Wang, Phys. Rev. D 83 (2011) 034031, arXiv:1011.4123 [hep-ph];

J.-W. Chen, J. Deng, H. Dong, Q. Wang, Phys. Rev. D 84 (2011) 039902 (Erratum).

23] B. Zhang, J. Phys. Conf. Ser. 420 (2013) 012035, arXiv:1208.1224 [nucl-th].

24] O. Fochler, J. Uphoff, Z. Xu, C. Greiner, Phys. Rev. D 88 (2013) 014018, arXiv: 1302.5250 [hep-ph].

25] X.-G. Huang, J. Liao, arXiv:1303.7214 [nucl-th].

26] B. Zhang, arXiv:1307.6234 [nucl-th].

27] J. Letessier, J. Rafelski, A. Tounsi, Phys. Rev. C 50 (1994) 406, arXiv:hep-ph/ 9711346.

28] M.L. Miller, K. Reygers, S.J. Sanders, P. Steinberg, Annu. Rev. Nucl. Part. Sci. 57 (2007) 205, arXiv:nucl-ex/0701025.

29] H.-J. Drescher, Y. Nara, Phys. Rev. C 75 (2007) 034905, arXiv:nucl-th/0611017.

30] H.-J. Drescher, Y. Nara, Phys. Rev. C 76 (2007) 041903, arXiv:0707.0249 [nucl-th].

31] U. Heinz, R. Snellings, Annu. Rev. Nucl. Part. Sci. 63 (2013) 123, arXiv:1301.2826 [nucl-th].

32] J.-P. Blaizot, F. Gelis, J. Liao, L. McLerran, R. Venugopalan, Nucl. Phys. A 873 (2012) 68, arXiv:1107.5296 [hep-ph].

33] J.-P. Blaizot, J. Liao, L. McLerran, Nucl. Phys. A 920 (2013) 58, arXiv:1305.2119 [hep-ph].

34] F.C. Auluck, D.S. Kothari, Math. Proc. Camb. Philos. Soc. 42 (1946) 272.

35] B.K. Agarwala, F.C. Auluck, Math. Proc. Camb. Philos. Soc. 47 (1951) 207.

36] R.L. Armendariz, UM1-32-59544.

37] S.S. Adler, et al., PHENIX Collaboration, Phys. Rev. C 71 (2005) 034908, arXiv:nucl-ex/0409015;

S.S. Adler, et al., Phys. Rev. C 71 (2005) 049901 (Erratum). [38] S.S. Adler, et al., PHENIX Collaboration, arXiv:1312.6676 [nucl-ex].