Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan

TL;DR

This paper reviews the current efforts and findings in searching for phase transitions, the critical point, and anomalous transport phenomena in Quantum Chromodynamics using the Beam Energy Scan program at RHIC.

Contribution

It provides a comprehensive overview of experimental observables, results, and future directions in the quest to map QCD phases and critical phenomena.

Findings

Experimental results suggest possible signals of the QCD critical point.

Fluctuation and correlation measurements are key to identifying phase transitions.

Further data and theoretical work are needed to confirm the existence of the critical point.

Abstract

We review the present status of the search for a phase transition and critical point as well as anomalous transport phenomena in Quantum Chromodynamics (QCD), with an emphasis on the Beam Energy Scan program at the Relativistic Heavy Ion Collider at Brookhaven National Laboratory. We present the conceptual framework and discuss the observables deemed most sensitive to a phase transition, QCD critical point, and anomalous transport, focusing on fluctuation and correlation measurements. Selected experimental results for these observables together with those characterizing the global properties of the systems created in heavy ion collisions are presented. We then discuss what can be already learned from the currently available data about the QCD critical point and anomalous transport as well as what additional measurements and theoretical developments are needed in order to discover these phenomena.

Tables

Table 1. Table 1: RHIC beam energy scan schedule (as originally approved by BNL management), event statistics, collision energies as well as the corresponding values of the chemical freeze-out temperature and baryonic chemical potential. These values are for the most central Au+Au collisions at RHIC. The 2019 BES-II runs not shown in this Table include a short RHIC electron-cooling test run (3M events at s N ​ N = 7.7 ​ GeV subscript 𝑠 𝑁 𝑁 7.7 GeV \sqrt{s_{NN}}=7.7\rm GeV ) as well as a series of fixed-target mode collisions (50M events at 7.7 ​ GeV 7.7 GeV 7.7\rm GeV , 50M events at 3.9 ​ GeV 3.9 GeV 3.9\rm GeV , 200M events at 3.2 ​ GeV 3.2 GeV 3.2\rm GeV and 300M events at 3.0 ​ GeV 3.0 GeV 3.0\rm GeV ). In addition to the originally planned energies, a new request for collisions with s N ​ N subscript 𝑠 𝑁 𝑁 \sqrt{s_{NN}} at about 17 ​ G ​ e ​ V 17 G e V 17\rm GeV (possibly taking 250 million events in year 2021) has also been put forward.

BES-II	STAR	BES-I	STAR	PHENIX	$\sqrt{s_{NN}}$	${\mu }_B$	$T_{ch}$
	Events (106)		Events (106)	Events (106)	(GeV)	(MeV)	(MeV)
		2010	238	1681	200	25	166
		2010	45	474	62.4	73	165
		2017	1200		54.4	92	165
		2010	86	154	39	112	164
		2011	32	21	27	156	162
2019	580	2011	15	6	19.6	206	160
2019	320	2014	13		14.5	264	156
2020-21	230	2010	7		11.5	315	152
2020-21	160	2008	0.3		9.2	355	140
2020-21	100	2010	3	2	7.7	420	139

Equations

$${\mathcal{L}}_{QCD}={\bar{\psi }}_f\left(i{\gamma }^{\mu }{D}_{\mu }-m\right){\psi }_f-\frac{1}{4}{G}_{\mu \nu }^a{G}_a^{\mu \nu }$$

$${\mathcal{V}}_{\mu }^i$$

$${\mathcal{A}}_{\mu }^i$$

$$⟨0|{\mathcal{A}}_{\mu }^j(x)|{\pi }^k(q)⟩=i{f}_{\pi }{q}_{\mu }{\delta }^{jk}{e}^{-iq\cdot x}$$

$$⟨0|{\partial }^{\mu }{\mathcal{A}}_{\mu }^j(x)|{\pi }^k(q)⟩=-f_{\pi }{q}^2{\delta }^{jk}{e}^{-iq\cdot x}=-f_{\pi }{m}_{\pi }^2{\delta }^{jk}{e}^{-iq\cdot x}$$

$$\delta {\mathcal{L}}_{quarks}={\bar{\psi }}_f\left[i\left({\gamma }_{\mu }{D}^{\mu }\right)-m\right]{\psi }_f$$

$${\Lambda }_V:$$

$${\Lambda }_A:$$

$${\psi }_R=\frac{1+{\gamma }_5}{2}\psi ,{\psi }_L=\frac{1-{\gamma }_5}{2}\psi .$$

$$J_R^{\mu ,i}=\frac{1}{2}\left(V^{\mu ,i}+A^{\mu ,i}\right)={\bar{\psi }}_R{\gamma }^{\mu }{\tau }^i{\psi }_R,{\partial }_{\mu }{J}_R^{\mu ,i}=0$$

$$J_L^{\mu ,i}=\frac{1}{2}\left(V^{\mu ,i}-A^{\mu ,i}\right)={\bar{\psi }}_L{\gamma }^{\mu }{\tau }^i{\psi }_L,{\partial }_{\mu }{J}_L^{\mu ,i}=0$$

$$U(1):$$

$$U(1{)}_A:$$

$${\partial }_{\mu }{J}_5^{\mu }=\left(\frac{1}{2{\pi }^2}\right)(qe\mathbf{E})\cdot (qe\mathbf{B})$$

$$\frac{d{N}_5}{dt}={\int }_V{d}^3\mathbf{r}\frac{\partial J_5^0}{\partial t}=\left(\frac{1}{2{\pi }^2}\right){\int }_V{d}^3\mathbf{r}(qe\mathbf{E})\cdot (qe\mathbf{B})=C_A{\int }_V{d}^3\mathbf{r}(qe\mathbf{E})\cdot (qe\mathbf{B}).$$

$${\partial }_{\mu }{J}_R^{\mu }=+\frac{C_A}{2}(qe\mathbf{E})\cdot (qe\mathbf{B}),{\partial }_{\mu }{J}_L^{\mu }=-\frac{C_A}{2}(qe\mathbf{E})\cdot (qe\mathbf{B})$$

$${\partial }_{\mu }{J}_5^{\mu }=-\frac{g^2}{16{\pi }^2}{G}_a^{\mu \nu }{\tilde{G}}_{\mu \nu }^a=\frac{1}{4{\pi }^2}(g{\mathbf{E}}^a)\cdot (g{\mathbf{B}}_a)=\frac{C_A}{2}(g{\mathbf{E}}^a)\cdot (g{\mathbf{B}}_a)$$

$$N_5(t\to +\infty )-N_5(t\to -\infty )=-\frac{g^2}{16{\pi }^2}\int dt{d}^3\mathbf{r}{G}_a^{\mu \nu }{\tilde{G}}_{\mu \nu }^a=C_A\int d^{4}x{\partial }^{\mu }{K}_{\mu },$$

$$Q_w=\frac{1}{2\pi }{\int }_{\Sigma }d\vec{l}\cdot (e\vec{A})$$

$$Q_w=-\frac{1}{32{\pi }^2}\int d^{4}x(g{G}_a^{\mu \nu })\cdot (g{\tilde{G}}_{\mu \nu }^a)=\frac{C_A}{2}{\int }_{boundary}d{\Sigma }^{\mu }\cdot K_{\mu },$$

$$N_R-N_L=N_5=2Q_w.$$

$$M=\frac{1}{N}\sum _{i=1}^N{s}_i$$

$$P(M)\sim \exp \left\{-\frac{N}{T}\left(F(M)-hM\right)\right\}.$$

$$\left(\frac{\partial F}{\partial M}\right)_{T}=h\quad\mbox{at $M=\langle M\rangle$}.$$

$$F(M)=\frac{1}{2}a(T)M^2+\frac{1}{4}b(T)M^4,\text{mbox}(mean-fieldfreeenergy)$$

$$⟨(\delta M{)}^2⟩=\frac{T}{N}{\left(\frac{{\partial }^{2}F}{\partial M^2}\right)}_T^{-1}=\frac{T}{N}{\left(\frac{\partial ⟨M⟩}{\partial h}\right)}_T.$$

$$⟨\delta \sigma (\mathbf{x})\delta \sigma (\mathbf{y})⟩=\int \mathcal{D}\sigma \mathcal{P}[\sigma ]\delta \sigma (\mathbf{x})\delta \sigma (\mathbf{y})\sim \exp (-|\mathbf{x}-\mathbf{y}|/\xi ),\mathrm{for}|\mathbf{x}-\mathbf{y}|\to \infty .$$

$$P[\sigma ]\sim \exp \left\{-\Omega [\sigma ]/T\right\},$$

$$\Omega [\sigma ]=\int d^3\mathbf{x}\left[U(\sigma )+\frac{1}{2}(\nabla \sigma {)}^2+\dots \right].$$

$$U(\sigma )=U_0+\frac{1}{2}{m}_{\sigma }^2(\sigma -{\sigma }_0{)}^2+\dots$$

$$\delta \sigma =\sigma -⟨\sigma ⟩,$$

$$⟨\delta \sigma (\mathbf{x})\delta \sigma (\mathbf{y})⟩=\frac{T}{4\pi |\mathbf{x}-\mathbf{y}|}\exp \left(-\frac{|\mathbf{x}-\mathbf{y}|}{\xi }\right),\text{mbox}(Gaussianapproximation)$$

Full text

††journal: Physics Reports

1. 1 Introduction

2. 2 QCD and Chiral Symmetry

3. 2.1 Chiral symmetry

4. 2.2 Anomaly

5. 2.3 Topology

6. 3 Phase Diagram of QCD

7. 3.1 Chiral symmetry restoration and deconfinement

8. 3.2 Phase transition at finite baryon density

9. 3.3 First-order transition and the critical point

10. 3.4 High baryon density regime

11. 4 Theory and Phenomenology of the Critical Point

12. 4.1 Critical phenomena

13. 4.2 The role of fluctuations

14. 4.3 Long-range correlations

15. 4.4 Mean-field critical behavior

16. 4.5 Higher-order moments and cumulants of the fluctuations

17. 4.6 The cumulants and the equation of state in parametric representation

18. 4.7 Mapping to QCD equation of state

19. 4.8 Baryon number cumulants near the critical point

20. 4.9 Fluctuation cumulants in heavy-ion collisions

21. 4.10 Finite time dynamics and limitations on the correlation length in heavy-ion collisions

22. 4.11 Experimentally observable fluctuations and correlations

23. 4.12 Translating critical point correlations into observable momentum correlations

24. 4.13 Higher-order cumulants

25. 5 Theory and Phenomenology of Anomalous Chiral Transport

26. 5.1 The many-body physics of the chiral anomaly

27. 5.1.1 Chiral Magnetic and Vortical Effects

28. 5.1.2 Chiral collective excitations

29. 5.1.3 Fluid dynamics with chiral anomaly

30. 5.1.4 Chiral kinetic theory

31. 5.2 Anomalous chiral transport in heavy ion collisions

32. 5.3 Experimental signals of anomalous chiral transport

33. 5.4 The beam energy dependence of anomalous chiral transport

34. 6 The Beam Energy Scan Program at RHIC

35. 6.1 Introduction

36. 6.2 Bulk properties and collectivity

37. 6.2.1 Particle spectra and ratios

38. 6.2.2 Collective flow and azimuthal asymmetries

39. 6.2.3 HBT and femtoscopic correlations

40. 6.3 Criticality

41. 6.3.1 Experimental data on proton, anti-proton, net-proton and net-charge high order cumulants

42. 6.3.2 Efficiency corrections

43. 6.4 Chirality

44. 7 Discussions on BES-I Results

45. 7.1 Discussion of the STAR data on proton number cumulants

46. 7.1.1 Proton factorial cumulants

47. 7.1.2 Baryon number conservation

48. 7.1.3 Volume fluctuation

49. 7.1.4 Rapidity dependence

50. 7.1.5 Bimodal distribution

51. 7.2 Discussions on the measurements for anomalous chiral transport

52. 7.2.1 The issues with background correlations

53. 7.2.2 Quantitative modeling of anomalous chiral transport

54. 7.2.3 Developing new measurement methods

55. 7.2.4 The isobaric collisions

56. 8 Summary and Outlook

57. 8.1 Present status

58. 8.2 Open questions

59. 8.3 Theory efforts

60. 8.4 Future Upgrades and Physics Program at RHIC and the global picture

61. A Introduction to Cumulants

1 Introduction

The understanding of matter and its phases has always been of fundamental interest in science. Matter in our physical world is governed by a hierarchical microscopic structure, reflecting various length scales: “normal” matter such as water or crystals is made out of atoms and molecules which by way of their interaction arrange themselves in various ways, such as liquid or gas or crystalline structures. At a shorter length scale, molecules consist of atoms and again the interactions among them result in a variety of configurations ranging from diatomic molecules to complex polymers. The atoms themselves are made of a positively charged nuclei surrounded by negatively charged electrons. The nucleus itself consists of protons and neutrons, which themselves are made of even more fundamental particles, quarks and gluons which are fundamental degrees of freedom of the Standard Model. In this review we are concerned with matter which is governed by the strong interactions. It therefore involves nucleons or more generally hadrons as well as quarks and gluons. The properties of strongly interacting matter is not only of fundamental interest but it also important for our understanding of the Universe [1, 2]. Indeed a few milliseconds after the “Big Bang” the entire Universe was occupied by strongly interacting matter of quarks and gluons. Also, highly compressed strongly interacting matter at densities of several times nuclear saturation density ($\rho_{0}\sim 0.16/{\rm fm}^{3}$) is expected to exist inside the compact stars and thus affects their properties [3, 4].

Matter can occur in various states, or phases. Most well known are probably the liquid, gas and solid phases of water. So the obvious question arises, if strongly interacting matter can also exhibit different phases. And indeed, soon after the discovery of Quantum Chromodynamics (QCD) [5], the theory of the strong interactions, and following the realization that QCD exhibits asymptotic freedom [6, 7], it was recognized that QCD likely predicts a transition from a hadronic phase to one dominated by quarks and gluons, the so called quark-gluon plasma (QGP) [8, 9, 10, 11, 12, 13]. However, already before the advent of QCD the possibility of a limiting temperature of the hadron gas was discussed and a quantitative prediction of $T_{limit}\simeq 170{\rm\,MeV}$ was obtained in the statistical bootstrap model of Hagedorn [14]. The existence of a new phase was confirmed in the first calculations using the lattice formulation of QCD, initially for pure $SU(2)$ gauge theory [15, 16, 17].

These results soon sparked ideas to create and study the quark-gluon plasma in high energy heavy ion collisions (see e.g. [11, 18, 19] ) extending early considerations of creating thermodynamically equilibrated matter in high energy hadronic collisions by Fermi [20], Landau [21], and Hagedorn [14]. These early initiatives have meanwhile led to very active and large research programs at the Bevatron of LBL, the Alternating Gradient Synchrotron (AGS) and the Relativistic Heavy Ion Collider (RHIC) of Brookhaven as well as the Super Proton Synchrotron (SPS) and the Large Hadron Collider (LHC) of CERN. The central discoveries at RHIC, and later confirmed and refined at the LHC, include a large azimuthal asymmetry, known as elliptic flow $v_{2}$, in the particle yields [22], as well as a strong suppression of high energy jets and heavy quarks [23]. The observed elliptic flow, which is consistent with predictions from nearly ideal hydrodynamics, together with the strong suppression of jets suggested that the system created in these collisions is strongly coupled. These findings required a paradigm shift, away from the expectations, based on asymptotic freedom, of a weakly coupled quark-gluon plasma [24, 25, 26, 27] (for a recent review see, e.g., [2]). Of course at very high temperatures and/or baryon density QCD will eventually become weakly coupled due to asymptotic freedom. To which extent these conditions can be reached in experiments at even higher energies or in the interior of neutron stars is at present not known.

Together with the advances in experiment, important developments occurred on the theoretical side. In particular, lattice QCD methods became sufficiently refined and powerful to perform continuum extrapolated calculations in full QCD with physical quark masses. These simulations for example showed that the transition at vanishing baryon density is a crossover [28] with a pseudo-critical temperature of $T\simeq 150{\rm\,MeV}$ [29, 30, 31, 32]. They also showed that the (approximate) chiral symmetry of QCD, which is spontaneously broken in the vacuum, is restored at high temperatures.

Due to the fermion sign problem111In the presence of a baryon-number chemical potential, the fermion determinant in the path integral of the partition function becomes complex, making importance sampling difficult if not impossible., lattice QCD methods are restricted to the region of vanishing or small baryon density. Also most experiments both at the highest energies at RHIC as well as at the LHC explore this region of nearly vanishing baryon density. Therefore, very little is known about the properties of strongly interacting matter at large baryon density. Many model calculations, see, e.g. Refs. [2, 33, 34, 1, 35, 36] for recent reviews, predict a first order phase transition at large baryon density or equivalently baryon number chemical potential. If true, this phase transition will end at a critical point since we know the transition to be a crossover at small baryon chemical potential. These model ideas are typically summarized in a (conjectured) phase diagram of QCD matter shown in Fig. 1. Clearly such a phase diagram has to be rather schematic. The only regions where we have firm knowledge are: (a) at low temperature and baryon chemical potential, where we have a dilute gas of hadrons, which are predominantly pions, where interactions are small corrections which can be systematically described in chiral perturbation theory and the experimental knowledge of hadronic interactions. (b) At small values of the baryon number chemical potential and finite temperatures ($T\gtrsim 130{\rm\,MeV}$) from lattice QCD, and, (c), for small temperatures close to the nuclear matter saturation density from the extrapolation of well tested nuclear forces and experiments of nuclear fragmentation [37, 38, 39] .

In order to experimentally explore the QCD phase diagram at finite net-baryon density, one needs to create systems with finite net-baryon density in heavy ion collisions. Since baryon number is conserved, the only way to increase the net-baryon density is to ensure that some of the baryons from the colliding nuclei are transported to the mid-rapidity region. This can be achieved by lowering the beam energy, and available particle production systematics [41] confirm that this strategy indeed works. Therefore, a systematic scan of heavy ion measurements over a range of beam energies enables the exploration of the high baryon density region of the QCD phase diagram and the search for the existence of a first order phase transition and its associated critical point. Such a beam energy scan (BES) program has been started at RHIC in 2010 (see e.g. [42]) and its next phase with improved beam quality, such as increased luminosity and smaller beam packages, and detector capability has just started in 2019. Also, several experiments at other facilities such as NA61/SHINE and HADES are able to measure some of the same observables at energies even lower than those achievable at RHIC. In addition, new experiments extending the reach of the RHIC beam energies towards lower energies are planned at the FAIR facility in Darmstadt (CBM), at NICA in Dubna (MPD), as well as at the CSR in Lanzhou (CEE). Finally, a fixed-target experimental program enabling much lower energy collisions at RHIC is underway, which allows measurements at center of mass energies down to $\sqrt{s_{\rm NN}}=3{\rm\,GeV}$, although with somewhat limited acceptance. These current and future programs provide unique opportunities for exploring and mapping the phases of QCD across a wide range of conditions in the laboratory.

The first set of measurements resulting from the RHIC beam energy scan made a number of intriguing observations, such as a non-monotonic dependence on the beam energy of some of the key observables and the disappearance at low energy of certain key signals observed at high energy. These observations underline the discovery potential for locating the Critical Point (CP) as a landmark of QCD phase diagram, as well as for unambiguously observing the anomalous chiral transport effects thus experimentally verifying chiral symmetry restoration. Therefore, and in view of the second phase of the beam energy scan, it appears to be appropriate to provide a review of the key physics motivations together with the experimental results of the RHIC beam energy scan so far. Special emphasis in this review will be given to observables which are sensitive to the QCD critical point, namely fluctuations of conserved charges, and to chiral restoration, namely various correlation functions which measure the induced currents from the chiral anomaly.

The rest of this review is organized as follows. In the next section we provide an overview of those aspects of QCD pertinent to the topics discussed in this review. Next we provide a more detailed discussion of the QCD phase diagram followed by a section on the physics of the critical point. We then turn to the dynamics induced by the chiral anomaly and how it can be used to extract information about chiral symmetry restoration. After that we present an overview of the various measurements, the implications of which will then be discussed in the context of the QCD phase diagram and chiral restoration. Finally we will summarize the current status and provide an outlook towards future experiments and developments.

2 QCD and Chiral Symmetry

In this section we will discuss the properties of QCD, the theory of the strong interaction, which are relevant to this review. The Lagrangian of QCD is given by

[TABLE]

Here $D_{\mu}=\partial_{\mu}-igA_{\mu}^{a}\lambda_{a}$ is the co-variant derivative, and $A_{\mu}^{a}$ is the gluon field of $SU(3)$ color symmetry with $a$ the color index, $\lambda_{a}$ the Gell-Mann matrices, and $G_{\mu\nu}^{a}$ the corresponding non-Abelian gauge field strength tensor [43]. In the following we will consider only the light (up and down) quarks and assume their mass, $m$, to be the same. The quark fields $\psi_{f}$ are then two-component spinors in flavor space, denoted by the index $f$, and we have suppressed the additional color indices. One of the most widely known properties of QCD is asymptotic freedom, i.e., the fact that the strength of the strong coupling constant, $\alpha_{s}=\frac{g^{2}}{4\pi}$, decreases at short distances or, equivalently, at large momentum scales. This property, first pointed out by Gross, Wilczek, and Politzer [6, 7] has been widely tested and confirmed in experiment and has been awarded the Nobel prize in 2004. Asymptotic freedom allows for a perturbative treatment of processes at high energies, such as deep inelastic scattering etc. Asymptotic freedom also implies an emergent energy scale $\Lambda_{QCD}\simeq 200\ \rm MeV$ at which the QCD running coupling becomes large and the theory gets strongly coupled. Equivalently this implies an emergent distance scale $1/\Lambda_{QCD}\simeq 1\ \rm fm$ that is empirically found to set the size of all hadrons containing light flavor valence quarks. At such low energy scales or long distances, QCD is highly non-perturbative, and its treatment requires either numerical methods such as lattice QCD or effective field theories which respect the symmetries of QCD. It is the low energy regime of QCD which is most relevant for the proprieties of QCD probed in heavy ion collisions, such as the phase structure, its collective behavior etc. Especially relevant for the topics discussed in this review are chiral symmetry, the axial anomaly as well as the topological structure of the QCD vacuum, which we will briefly discuss in the following. 222We note that there are other highly important non-perturbative features of QCD such as confinement, which are not much discussed in the present review. We refer interested readers to e.g. [44]. Given the quark fields one can define the (iso-vector) vector and axial vector currents, ${\mathcal{V}}_{\mu}^{i}$ and ${\mathcal{A}}_{\mu}^{i}$, given by:

[TABLE]

which play an important role in the context of chiral symmetry which we will discuss below.

2.1 Chiral symmetry

Chiral symmetry is a symmetry of QCD in the limit of vanishing (light) quark masses. In reality, however, the current masses of the up and down quarks are finite, $m_{u,d}\simeq 10{\rm\,MeV}$, but very small compared to the scale of QCD $\Lambda_{QCD}\simeq 200{\rm\,MeV}$. Therefore, chiral symmetry may be considered an approximate symmetry of the strong interactions. Before we provide a proper definition of chiral symmetry in QCD let us briefly review its history.

Long before QCD was known to be the theory of the strong interactions, phenomenological indications for the existence of chiral symmetry came from the study of the nuclear beta decay. There one finds that the weak coupling constants for the vector and axial-vector hadronic-currents, $C_{V}$ and $C_{A}$, did not (in case of $C_{V}$) or only slightly ($25\%$ in case of $C_{A}$) differ from those for the leptonic counterparts. Consequently strong interaction ‘radiative’ corrections to the weak vector and axial vector ‘charge’ are essentially absent. This absence of radiative corrections indicates the presence of conserved currents and thus symmetries, similar to the familiar case of the electric charge, which is protected from radiative corrections as consequence of electric charge conservation. In case of the vector current, the underlying symmetry is the well known isospin symmetry of the strong interactions and the hadronic vector current is identified with the isospin current. The identification of the axial current, on the other hand is not so straightforward. Contrary to the vector or iso-spin symmetry, which is clearly visible in the hadronic mass spectrum – degeneracy of proton and neutron or the three pion states, for example – no such symmetry exists for states related by axial symmetry. For instance, the difference between the masses of the $\rho$ and $a_{1}$ mesons, $m_{a_{1}}-m_{\rho}\simeq 500{\rm\,MeV}$, is of the order of QCD scale, $\Lambda_{QCD}$, rather than the electromagnetic or weak scale.

However there is additional phenomenological evidence that the axial current is (almost) conserved. Since the pion is a pseudo-scalar, the weak decay of the pion with four-momentum $q$ involves the axial current $\mathcal{A}_{\mu}^{j}$ in Eq. (3). Since $\mathcal{A}_{\mu}^{j}$ is a Lorentz vector and since we are considering a pion with given four-momentum, the matrix element has to be proportional to $q_{\mu}e^{-iq\cdot x}$. Thus the matrix element can be written as

[TABLE]

where $f_{\pi}=93{\rm\,MeV}$ is the pion decay constant and $j,k$ are isospin indices. From its divergence

[TABLE]

we find that $\partial^{\mu}\mathcal{A}_{\mu}^{j}(x)\sim m_{\pi}^{2}$. Given that the mass of the pion is very small, the divergence of the axial current is also small, and thus the axial current is “almost” conserved. This observation led to the so called partially conserved axial current (PCAC) hypothesis, which typically refers to the relation (5). The above relations, Eqs.(4) and (5), suggest that the pion contributes to the axial current and that the divergence of the axial current may serve as an interpolating pion field. This observation led to further developments, such as the Goldberger-Treiman relation [45], which relates the strong pion nucleon coupling constant $g_{\pi}$ with the weak pion decay constant, $f_{\pi}$, a nice example of the power of symmetries. Also, the so-called Weinberg sum rules [46] which are concerned with moments of the spectral functions of the vector and axial currents, predict the mass difference of the $\rho$ and $a_{1}$ mesons to reasonable accuracy. The mass splitting is simply a result of the contribution of the pion to the axial channel, which has no equivalent in the vector channel.

The above examples, which were obtained prior to the advent of QCD, show that there is indeed a (partial) axial symmetry in the strong interaction. But why is this symmetry not visible in the mass spectrum? The answer to this is the mechamism of “spontaneous broken symmetry” first discovered by Y. Nambu [47, 48] in the context of the conserved axial current of the strong interaction, and he was awarded the 2008 Nobel prize for this work. The essential concept of spontaneous symmetry breaking is the following: while the Hamiltonian possesses the symmetry, its ground state does not. A familiar example is a ferromagnet, where the interaction between the spins is invariant under rotation, but the ground state, which has a magnetic moment pointing in a certain direction, obviously breaks rotation symmetry. An important consequence of the spontaneous symmetry breaking, as expressed in the Goldstone theorem [49], is the existence of a massless mode, the so-called Goldstone mode. While a more formal proof of the Goldstone theorem may be found in textbooks such as [43], its essence may be understood in the simple example shown in Fig. 2. In both (a) and (b) we see a potential (Hamiltonian) which is symmetric under rotations. However, in (a) the ground state indicated by the little ball, is in the middle, so that both the potential and ground state are invariant under rotations. In (b), on the other hand, the ground state, i.e., the position of the little ball obviously breaks the rotational symmetry. Thus the symmetry is spontaneously broken. However because of the rotational symmetry of the potential, rotational excitations, or small excitations in the $(y,\pi)$ direction, do not cost any energy. Or in other words, we have a “massless” mode. Excitations in the radial or $(x,\sigma)$ direction, on the other hand, cost energy and are thus “massive”. Compare this with the situation in (a). There excitations in both the $(x,\sigma)$ and $(y,\pi)$ direction are massive. And, due to the symmetry of the potential, they are degenerate. Thus spontaneous breaking leads to a non-degenerate excitation spectrum and the symmetry reveals itself in the appearance of a massless Goldstone mode. Indeed, the very existence of a “flat direction” of the potential, and the resulting Goldstone mode, in the spontaneously broken vacuum is guaranteed by the original symmetry of the Hamiltonian. In case of the strong interaction, the Goldstone bosons are the pions. The massive radial excitation is often referred to as the $\sigma$ meson, although such a meson has not been identified as a well defined state in experiment. The reason is that the $\sigma$, which carries the quantum number of the vacuum, is likely to mix with many other states, and thus broadens considerably due to, e.g., decay into two or more pions. 333For example in the pion-pion phase shift, a broad feature occurs at an invariant mass of about $500\ \rm MeV$ that could be plausibly due to the $\sigma$ meson. Finally, similar to the ferromagnet, where the magnetization disappears above the Curie temperature, one expects the chiral symmetry to be restored at sufficiently high temperatures.

After this historical introduction, we next turn to QCD. Chiral symmetry, is a symmetry of the quark and flavor sector of QCD and thus it is sufficient to consider only the fermionic part of the QCD Lagrangian (1),

[TABLE]

where we restrict ourselves again to the light (up and down) flavors.

Chiral symmetry is related to the following vector and axial-vector transformations of the quark fields

[TABLE]

where $\vec{\tau}=(\tau_{1},\tau_{2},\tau_{3})$ represents the isospin generators, which are Pauli matrices acting upon the flavor indices of the quark fields, and $\vec{\Theta}=(\theta_{1},\theta_{2},\theta_{3})$ are the corresponding transformation angles. The above transformations do not involve the gluon fields, and it is easy to see that, in the limit of massless quarks, $m=0$, often referred to as the chiral limit, the QCD Lagrangian is invariant under both $\Lambda_{V}$ and $\Lambda_{A}$. The conserved Noether currents (see, e.g., Ref. [43]) associated with these symmetries are the axial and vector currents, Eqs. (2) and (3).

Before we discuss chiral symmetry further let us step back and examine the special case of massless quarks or, generally, massless fermions a bit closer. In the absence of any fermion masses, the Dirac Hamiltonian reduces to $\hat{H}\to\gamma^{0}\vec{\gamma}\cdot\mathbf{p}$ and it commutes with $\gamma_{5}$. Therefore, we can characterize the eigenstates of $\hat{H}$ in terms of eigenstates of $\gamma_{5}$. This is done simply by the following decomposition of a given spinor into its left-handed and right-handed component

[TABLE]

with $\gamma^{5}\Psi_{R/L}=(+/-)\Psi_{R/L}$. So instead of classifying the fermions for example by spin-up and spin-down, in case of massless fermions it is more appropriate to sort them into left- and right-handed states, as they are eigenstates of both Hamiltonian and $\gamma_{5}$. For free fermions, where one frequently works with helicity states, the helicity of a left/right handed fermion is indeed left/right handed, which explains the origin of the concept of handedness or chirality in Greek. And as long as we are dealing with massless fermions the handedness or chirality is conserved classically. That is, a right-handed (left-handed) fermion with its spin parallel (antiparallel) to its momentum remains right-handed (left-handed).

Returning to QCD, in the limit of vanishing quark masses, both vector and axial vector are conserved and thus we can form left and right handed combinations which are both conserved

[TABLE]

Thus, both left and right handed quarks form individually conserved (iso-vector) currents. The combined symmetries for right handed and left handed quarks is, therefore, referred to as chiral symmetry, which is governed by the group $SU(2)_{L}\times SU(2)_{R}$.

It is easy to see that the scalar combination of the quark fields, $\bar{\psi}\psi$, which comes with the mass term in the Lagrangian, is not invariant under the axial vector transformation. The vector transformation, $\Lambda_{V}$, on the on the other hand, leaves the mass term unchanged. Therefore, for real QCD, where even the light quarks have a finite but small mass, the axial vector symmetry is (slightly) broken, whereas the vector symmetry remains intact. However, as already discussed, the dominant effect is the spontaneous breaking of the chiral symmetry, where the ground state does not posses the symmetries of the Hamiltonian. In QCD this happens due to the dynamics in the gluon fields. In particular the so-called instantons, which characterize certain topological fluctuation of the gluon field, are believed to be the main source for the spontaneous chiral symmetry breaking [50]. Therefore, even in case of vanishing quark masses chiral symmetry would be (spontaneously) broken, and the explicit breaking, due to the small current quark masses is only a small correction. This is evidenced by the fact that in the chiral limit the pion mass would be zero while in nature $m_{\pi}\simeq 139.6{\rm\,MeV}$, which is still very small compared to that of other hadrons. As a result of the spontaneous symmetry breaking, the vacuum expectation value of the scalar bi-linear quark operator, $\left\langle\bar{q}q\right\rangle$ obtains a finite value. The chiral condensate, $\left\langle\bar{q}q\right\rangle$, is not invariant under the chiral transformation, and may serve as an order parameter for the spontaneous symmetry breaking.

As one increases the temperature, thermal fluctuations suppress the instantons and as a result the chiral condensate $\left\langle\bar{q}q\right\rangle$ will vanish and chiral symmetry is restored, in analogy to the a ferro-magnet, where above the Curie temperature the (spontaneous) magnetization vanishes. That this is indeed the case has been established by lattice QCD calculations [31, 51] and in Fig. 3 we show the temperature dependence of the chiral condensate.444What is plotted is actually the ratio of the chiral condensate at finite temperature over that at $T=0$. In addition some subtractions had to be made in order to properly renormalize this quantity. For details see Ref. [51]. We see that at a temperature of $T\gtrsim 180{\rm\,MeV}$ the chiral condensate has essentially vanished and thus chiral symmetry is restored. As a consequence the chirality of the light quarks becomes well defined, which is an essential feature for the possible occurrence of the so-called chiral magnetic effect which we will discuss in detail in Section 5.

The system created in heavy ion collisions at RHIC and LHC is expected to reach the chiral transition temperature and it would be interesting and important to devise a measurement which would confirm the restoration of chiral symmetry. Chiral symmetry predicts that the vector and axial-vector spectral function should be the same. The vector spectral function is accessible via dilepton measurements. These have been carried out over a wide range of collision energies, from $\sqrt{s_{NN}}\simeq 2.2{\rm\,GeV}$ up to $\sqrt{s_{NN}}=200{\rm\,GeV}$. The analysis of these data [52, 53] indicate that the spectral function flattens, i.e. the peak associated with the $\rho$-meson disappears. Unfortunately there is no easy access to the axial spectral function and one has to rely on model estimates to see if chiral symmetry is restored. One such attempt [54] finds chiral symmetry would be restored and that both the vector and axial vector spectral functions become featureless, similar to those of free quarks.

Another approach is to look for the so-called Chiral Magnetic Effect (CME), an anomalous chiral transport phenomenon whose presence requires that the chirality of the quarks in the hot matter is well defined. The correlation between the spin and momentum of a chiral fermion, which would be completely spoiled by the vacuum chiral condensate mixing up right-/left-handed sectors, is an essential ingredient of CME. Thus experimental proof for its existence would immediately imply that chiral symmetry is restored in these collisions. This idea and the present status of experimental searches and their interpretation is discussed in detail in Sections 5 and 7.2.

So far we have concentrated on the vector and axial vector transformations $\Lambda_{V}$ and $\Lambda_{A}$ in the flavor $SU(2)$ iso-vector channels, Eq. (7), since they constitute what is commonly referred to as chiral symmetry. However QCD is also invariant under the analogous $U(1)$ transformations:

[TABLE]

The former, $U(1)$, which holds even in case of finite (and non-equal) quark-masses, gives rise to the conservation of baryon number (or equivalently quark number). The latter, $U(1)_{A}$, which appears to be a symmetry of QCD in the chiral limit, is actually broken by quantum effects, the so-called anomaly which we discuss next.

2.2 Anomaly

In the classical limit, the axial $U(1)_{A}$ symmetry would lead to a conserved current (i.e. its corresponding Noether current) $J^{\mu}_{5}=\bar{\psi}\gamma^{\mu}\gamma_{5}\psi$ in the case of massless quarks. In quantum field theory, however, the anomaly gives rise to yet another type of symmetry breaking, in addition to the explicit or the spontaneous breaking discussed previously. Even if the quarks are massless, upon the quantization of the Dirac field when coupled to any (Abelian or non-Abelian) gauge theory, the axial symmetry is explicitly broken and the axial current is thus not conserved. This phenomenon, which is different from the aforementioned spontaneous symmetry breaking, is the well-known chiral anomaly, often also referred to as axial anomaly or triangle anomaly. It was first discovered in the 1960’s in the effort to try to understand the neutral pion decay, $\pi^{0}\to 2\gamma$, with pioneering works from Adler as well as from Bell and Jackiw (thus sometimes also called ABJ anomaly) in the QED context [55, 56]. It was shown later in an insightful paper by Fujikawa in the path integral formalism that the chiral anomaly is an essential and general feature for the quantum description of chiral fermions with gauge interaction [57]. The chiral anomaly has now become an indispensable and integral ingredient in many basic aspects of quantum field theories and string theories.

Let us examine the chiral anomaly in a more concrete situation, with a single species of massless Dirac fermions of electric charge $qe$ coupled to an electromagnetic field $A_{\mu}$. The anomaly relation reads:

[TABLE]

where $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic fields. The constant in front of the right-hand side is the anomaly coefficient $C_{A}=\frac{1}{2\pi^{2}}$. Clearly this implies that the axial charge, $N_{5}\equiv\int_{V}d^{3}\mathbf{r}\,J_{5}^{0}$, is no longer conserved. This is evident upon integrating the above equation over a certain spatial volume $V$ (with vanishing current on the boundary): 555We note that in condensed matter physics a nonzero axial charge is generated in samples of the so-called Dirac or Weyl semimetals [58] by simply applying parallel electric and magnetic fields.

[TABLE]

It deserves mentioning that the anomaly relation may be viewed as occurring for the left-handed and right handed quarks separately:

[TABLE]

Each sector “suffers” half of the axial symmetry violation but they together maintain vector current conservation $\partial_{\mu}J^{\mu}=0$ which of course is required by gauge invariance.

A unique feature of the violation of the axial symmetry via the anomaly, is that the anomaly coefficient $C_{A}$ is universal, i.e. independent of the coupling strength and other dynamical details of the system. For example, in the context of perturbative computation of the anomaly relation, it first arises at one-loop level via the famous triangle diagram where one immediately obtains the full value of $C_{A}$. Naively one would expect that there will be all kinds — in fact an infinite number — of higher order diagrams that would contribute and thus modify the coefficient into a series of expansion ordered by the increasing power of the coupling constant (as is typically the case for a perturbative calculation). Interestingly, all the higher order contributions simply add to zero, thus leaving the coefficient $C_{A}$ intact. That is, the anomaly coefficient, $C_{A}$, is dictated by the one-loop contribution. Calculations based on nonperturbative methods show the same value for $C_{A}$ as well. As we will see later, such a universal value of the anomaly coefficient has profound consequences for its macroscopic manifestation.

Let us now turn to the anomaly in QCD in the chiral limit, $m_{q}\simeq 0$. In this case the non-Abelian anomaly occurs similarly through the triangle diagram and it is analogous to the Abelian anomaly Eq. (11). One simply replaces the QED coupling with the QCD coupling and electromagnetic fields with the gluon fields as well as sums over color index. Thus one obtains the following anomaly relation (for each flavor):

[TABLE]

with $g$ the color charge and $\mathbf{E}^{a},\mathbf{B}_{a}$ the chromo electric and magnetic fields.

Clearly this implies that the axial charge (per flavor) in a certain spatial volume could be changed by gluon field configurations according to $\frac{dN_{5}}{dt}=\frac{C_{A}}{2}\,\int_{V}d^{3}\mathbf{r}\,(g\mathbf{E}^{a})\cdot(g\mathbf{B}_{a})$. That is, chromo electric and magnetic fields in parallel (or anti-parallel) configurations will induce the generation of nonzero axial charges. As it turns out, such configurations naturally occur during the very initial states, often referred as the glasma, in relativistic heavy ion collisions [59, 60]. The strong chromo electric and magnetic fields in the glasma organize themselves into many flux tubes with collinear $\mathbf{E}^{a}$ and $\mathbf{B}_{a}$ fields extending along the beam direction and with a typical transverse size determined by the so-called saturation scale $Q_{s}^{2}=(1\sim 2){\rm\,GeV}^{2}$, which characterizes the gluon wave function of a nucleus [59, 60]. Consequently nonzero axial charge will be generated inside each flux tube with fluctuating signs depending upon whether the $\mathbf{E}^{a},\mathbf{B}_{a}$ are parallel or anti-parallel. Rough estimates as well as recent real-time lattice simulations have suggested that sizable axial charge fluctuations can be generated through the glasma flux tubes in heavy ion collisions [61, 62, 63].

A much more nontrivial question is whether axial charge could be generated globally. To appreciate this issue, one may perform integration of Eq. (14) over the whole spacetime and notice that its right-hand side can actually be rewritten as a full derivative:

[TABLE]

with the topological current $K_{\mu}\equiv\frac{1}{8}\epsilon_{\mu\nu\rho\sigma}\left(\tilde{A}_{a}^{\nu}\partial^{\rho}\tilde{A}_{a}^{\sigma}-\frac{1}{3}f_{abc}\tilde{A}_{a}^{\nu}\tilde{A}_{b}^{\rho}\tilde{A}_{c}^{\sigma}\right)$ for a given gluon configuration $\tilde{A}_{a}^{\mu}=gA_{a}^{\mu}$. This integral would evaluate to zero, unless there is a nonzero contribution from the spacetime boundary due to gluonic configurations with highly nontrivial boundary topology. These do exist in QCD, known as the instantons and sphalerons,666We note that these configurations are fluctuations of the QCD vacuum such that $\left\langle N_{5}\right\rangle=0$, i.e. the CP symmetry is not violated by QCD. which we shall discuss next.

2.3 Topology

Topology is a mathematical concept pertaining to the properties of shapes (in its general yet abstract sense) that would not change under arbitrary continuous deformations. It has wide applications in many branches of physics. In the context of field theories, there exist configurations with nontrivial topological structures in various spacetime dimensions. In particular, such nontrivial topological configurations are found to play important roles in non-Abelian gauge theories including QCD.

Let us first use a simple example, the vortex in two spatial dimension, to illustrate the emergence of topological configurations. This is shown in Fig. 4. Let us consider possible configurations of a $U(1)$ gauge field $A_{\mu}$ on the boundary ${\Sigma}$ at spatial infinity which is a circle or $S(1)$ (represented by the blue solid circle in Fig. 4). The possible values of the field $A_{\mu}$ on ${\Sigma}$ must be pure gauge 777By “pure gauge” one means field configurations $A_{\mu}$ with zero field strength $F_{\mu\nu}=0$ and thus carrying no energy or action. in order to ensure a finite action. In other words, $A_{\mu}$ should be zero up to a gauge transformation characterized by an internal $U(1)$ phase angle $\theta(x)$ (represented by the green arrows in Fig. 4). Therefore at each point on ${\Sigma}$, the field $A_{\mu}$ takes a value from a group of gauge-equivalent $A_{\mu}$ labeled by the internal angle $\theta$. The space of gauge-equivalent $A_{\mu}$ or angles, which is represented by the dashed red circles at the boundary in Fig. 4 is often referred to as the functional space $\mathcal{F}$. In other words, we have a map from the configuration space, $\Sigma$ to the functional space, $\mathcal{F}$, which in our simple example are both circles. Now let start from an arbitrary point P (with and initial value for the field $\theta_{0}$) on the boundary ${\Sigma}$ and go around for a full circle back to the point P while watching the change of the internal phase angle $\theta$ along the way. Upon returning to P, $\theta$ must return to its initial value $\theta_{0}$: this however could be realized in many ways, e.g. allowing $\theta$ to stay constant888It actually does not have to be constant, it may fluctuate around the initial angles while going around the blue circle. along $\Sigma$ (left configuration of Fig. 4), or changing for a full $2\pi$ along $\Sigma$ (middle configuration of Fig. 4), or changing for a full $4\pi$ along $\Sigma$ (right configuration of Fig. 4), etc. Clearly, in these distinctive configurations, the gauge field winds around its internal space by zero, one, or two times while going around the spatial boundary once. This brings the concept of an integer number $Q_{w}$ whose magnitude simply counts the number of times the gauge field winds around. One can further assign a positive or negative sign to $Q_{w}$, depending on whether one is winding in a counterclockwise or clockwise direction (while assuming always going counterclockwise on the spatial boundary $\Sigma$). This intuitive observation can be made more precise by the following mathematical definition: 999We note that if one considers this 2D system as actually embedded in a 3D space, then this definition is counting the quantized magnetic flux through the vortex.

[TABLE]

which is a line-integral of the gauge field over the spatial boundary $\Sigma$. Let us briefly show that this definition actually works. As already said, on the boundary the gauge field is pure gauge, given by the local gauge phase angle $\theta(x)$ via $\vec{A}=\frac{1}{e}\nabla\theta(x)$. It is most transparent if one performs the above integral with polar coordinates $(r,\phi)$, in which case $d\vec{l}\to Rd\phi$ (with $R$ the radius of the boundary $\Sigma$) while $(e\vec{A})\to\frac{1}{R}\partial_{\phi}(\theta)$. It shall now become obvious that, going around $\Sigma$ from point P back to P, the above integral yields $Q_{w}=[\theta(P)_{final}-\theta_{0}]/(2\pi)$ where $\theta(P)_{final}$ must end up as $(\theta_{0}+2\pi\,n)$ (where $n$ is any integer) for periodicity. Thus the so-defined $Q_{w}$ indeed counts the number of windings.

It is not difficult to recognize that configurations with the same $Q_{w}$ could be continuously deformed into each other while configurations with different $Q_{w}$ could never be. That is precisely a statement about topology: configurations with the same $Q_{w}$ are topologically equivalent, while configurations with different $Q_{w}$ are topologically distinctive. The number $Q_{w}$ is the topological number labeling various topological sectors. Mathematically speaking, one defines a mapping from the spacetime manifold $\Sigma$ to the field space manifold $\mathcal{F}$. Such a mapping could be topologically trivial or nontrivial depending upon the specific $\Sigma$ and $\mathcal{F}$. There are a number of famous examples for topologically nontrivial configurations in field theories, including the kink (domain wall) in one spatial dimension, the vortex in two spatial dimensions discussed above, the monopole in three spatial dimensions, as well as the instantons and sphalerons in four spacetime dimensions. Interested readers are referred to a number of excellent and pedagogical review articles on this topic, see e.g. [50, 64, 65].

Focusing our discussions on QCD, the relevant topological configurations are the instantons and sphalerons arising from the gluon field. It is most convenient to discuss them in the four dimensional Euclidean spacetime. In this case, the boundary is the $O(3)$ hypersphere at infinity in 4D (Euclidean) spacetime. As in our example, any gluon configuration $A_{\mu}(x)$ with finite action can only take values on the boundary which are pure gauge. Therefore the field space $\mathcal{F}$ would be the collection of all pure gauge field configurations connected with each other through gauge transformations which in the case of QCD would be $\mathrm{SU(3)}$. This mapping may be topologically nontrivial, and gluon field configurations can be labelled by a topological winding number $Q_{w}$ (in analogy to the winding number $Q_{w}$ of the 2D vortex example):

[TABLE]

where $g$ is QCD gauge coupling, $G$ and $\tilde{G}$ are gluon field strength tensor and dual tensor, and $a=1,2,...8$ is the adjoint color index to be summed. Explicit solutions with nonzero integer $Q_{w}$, i.e. the instantons, were discovered a while ago [66]. They were found to play vital roles in understanding many fundamental aspects of QCD from vacuum to finite temperature and density: see reviews in [50, 67].

Next, let us highlight a few important features of such topological configurations with nonzero integer $Q_{w}$. We first note that, given Eq. (17), the gluon field is strong and thus non-perturbative, in particular if the gauge coupling $g$ is weak. Indeed, Eq. (17) suggests that the gluon field strength scales as $G\sim\sqrt{Q_{w}}/g$ and its contribution to the action is $S\sim G^{2}\sim Q_{w}/g^{2}$. A precise calculation reveals the action of such a configuration, which is a solution to the classical equation of motion, is given by $Q_{w}(8\pi^{2}/g^{2})$. This has nontrivial implications: as the gauge coupling becomes stronger, the action “cost” of such configurations becomes less. As a result these configurations become more important at strong coupling or equivalently at the low energy scale in QCD. Furthermore the exact solution shows that the total action of such configurations, instead of spreading widely over large spacetime volume, is more concentrated across a domain of certain size, typically on the order of a fraction of $1\,\rm fm$ in QCD.101010The scale is actually set by the running of the coupling constant as the action $S\sim 1/g^{2}$. Furthermore, let us emphasize that the above topological winding number of gluon field is ${\hat{\mathcal{P}}}$-odd and ${\hat{\mathcal{C}}}{\hat{\mathcal{P}}}$-odd, which perhaps becomes apparent upon the rewriting $G^{\mu\nu}_{a}\tilde{G}_{\mu\nu}^{a}\to\mathbf{E}^{a}\cdot\mathbf{B}_{a}$. In short, the gluonic topological transitions are non-perturbative phenomena naturally occurring in QCD and are actually important in the regime where the gauge coupling is strong. Accompanying such transitions are the creation of ${\hat{\mathcal{P}}}$-odd and ${\hat{\mathcal{C}}}{\hat{\mathcal{P}}}$-odd local domains. Needless to say, an experimental verification of such topological transitions would be of fundamental interest.

This however is nontrivial as the gluons do not carry the kind of quantum numbers (such as electric charge or baryon number) that could be readily measured. The gluonic topological fluctuations could be better “seen” if their unique features would translate into the quark sector. And indeed, such a “translation” is readily provided by the chiral anomaly, as discussed in the previous section. By comparing Eq. (15) with Eq. (17) one immediately recognizes a key relation between the change of chirality in the quark sector and the gluonic field topology:

[TABLE]

The above relation bears a deep connection to mathematics, being essentially a special case of the Atiyah-Singer index theorem [68] for the Dirac operator in the instanton background fields. Therefore, every gluonic topological transition induces a corresponding fluctuation in the quark chirality imbalance. Consequently, topological transitions in the gluon sector are always accompanied by considerable fluctuations of the net chirality in local domains of matter, and an experimental measurement of the latter would be a unique and direct means for probing the former. To which extent this is possible we will discuss in Section 5.

3 Phase Diagram of QCD

3.1 Chiral symmetry restoration and deconfinement

The phase diagram of QCD is a map of thermodynamic states of QCD in the plane of temperature $T$ and baryon chemical potential $\mu_{B}$. It is shaped by the fundamental properties of QCD: spontaneous chiral symmetry breaking and confinement. If quark masses were zero, the QCD Lagrangian would possess exact chiral symmetry, which would be spontaneously broken in the vacuum state at $T=\mu_{B}=0$ and restored at sufficiently high temperatures. On the other hand, if there were no light quarks at all, the color charges would be linearly confined in the vacuum by color flux tubes or QCD strings, while the string tension would vanish at sufficiently high temperatures leading to deconfinement.

Strictly speaking, in our exquisitely imperfect world, neither of these two fundamental properties are realized exactly. The chiral symmetry is, of course, explicitly broken by non-zero quark masses. Meanwhile, confinement, in the rigorous sense of the infinite energy required to separate two opposite color sources, does not hold exactly unless quarks are absent (or infinitely heavy): string breaking via creation of quark-anti-quark pairs limits the range of the linear growth of the confining potential.

An alternative characterization of confinement in the vacuum state is based on the fact that the physical spectrum of QCD contains only color-singlet states, i.e. the various hadrons. This definition, however, is not directly useful for studying the phase diagram. Rather than being a consequence of dynamics, for a non-abelian gauge theory such as QCD, the Hilbert space of physical states must contain only color-singlet states by definition, as a necessary condition for quantization of the theory.  111111 For many subtle facets of confinement, especially in a theory like QCD with dynamical quarks, see, e.g., Ref. [44].

Both idealized limits discussed above – massless quarks and infinitely massive (or absent) quarks have their own order parameters. In the massless limit, the chiral condensate (as discussed in Sec. 2.1) is a local order parameter identifying different phases at finite temperature and density: being zero in the symmetric phase and being non-zero in the symmetry-breaking phase. Somewhat similarly in the limit without dynamical quarks (i.e. in a pure non-Abelian gauge theory), one can define an order parameter by considering finite-temperature theory as a 4-dimensional theory in Euclidean space (with compactified imaginary time direction): Polyakov loop L\equiv\left\langle\frac{1}{N_{c}}{\rm Tr}\left\{\hat{P}\ \exp\left(i\int_{0}^{1/T}dx_{4}A_{4}\right)\right\}\right\rangle\ where $N_{c}=3$ is the number of colors and $\hat{P}$ signifies path-ordering for the integration of gauge field component $A_{4}$ along temporal direction $x_{4}$. As shown in [69] there is a $Z_{N_{c}}$ symmetry of the theory for which $L$ is an order parameter. This symmetry is present at zero temperature (vacuum) but is broken at high temperature signifying deconfinement transition in pure gauge theory. Thus the change of the Polyakov loop from zero (in vacuum) to nonzero value at finite temperature signifies deconfinement transition in a pure-gauge theory.

The $Z_{3}$ symmetry is, however, explicitly broken in QCD because of quarks and the Polyakov loop is nonzero both in vacuum and at high-temperatures and thus cannot be used as an order parameter for rigorous definition of deconfinement in QCD. Likewise the chiral symmetry is also explicitly broken by nonzero physical quark masses in QCD so the chiral condensate is not an order parameter in QCD. Nevertheless, examining the behavior of the chiral condensate and the Polyakov loop at finite temperature and density is found to be very useful for studying QCD phase diagram, as we shall discuss further.

Notwithstanding chiral symmetry’s approximate nature, the smallness of quark masses relative to the relevant QCD scale makes this symmetry a useful guide to the states of QCD at finite temperature and density. Lattice calculations at zero baryon density show the rapid drop of the measure of the spontaneous breaking of the chiral symmetry – the chiral condensate – as a function of temperature (see Fig. 3). This behavior is expected due to the restoration of the chiral symmetry at high temperature. However, consistent with the approximate nature of the symmetry, the chiral condensate does not drop to zero and no thermodynamic singularity or discontinuity has been found in lattice calculations (at $\mu_{B}=0$) [70, 28].

As a function of increasing temperature QCD undergoes the transition from a thermal state of hadron gas, which is a pion-dominated, dilute and nearly ideal gas for sufficiently low temperatures, to a much denser state of quark-gluon plasma, which is nevertheless also nearly ideal at sufficiently high temperatures due to asymptotic freedom. The temperature of the transition region is determined largely by the confinement scale $\Lambda_{\rm QCD}\sim 200$ MeV (inverse of the confinement radius). The fact that this transition is associated with deconfinement can be seen from the first-principle lattice calculations of thermodynamic functions. The main feature, largely due to liberation of colored gluons, is a dramatic rise of the effective number of the degrees of freedom, measured, e.g., by entropy density in units of temperature, $s/T^{3}$ as seen in Fig. 5.

Lattice calculations show that this color liberation transition is not abrupt, and is not accompanied by a thermodynamic singularity. Rather, it is a smooth crossover [70, 28]. This is in agreement with the observation that the confinement-deconfinement transition is not well-defined in the presence of quarks, i.e., there is no order parameter, such as string tension or the Polyakov loop expectation value, whose non-analytic behavior would signal a true phase transition [69]. Similarly, because of the explicit breaking of the chiral symmetry by the quark masses, there is no order parameter for the chiral symmetry breaking which would have to be strictly zero when the system is in the state where the symmetry is restored.

In contrast, in the idealized chiral limit $m_{q}\to 0$, the thermodynamic realization of the chiral symmetry (in this case exact on the level of the QCD Lagrangian) qualitatively changes at a certain transition temperature (from being broken to being restored). This, by a standard textbook argument [73], has to lead to a true thermodynamic singularity – a phase transition. The chiral symmetry plays a central role in this argument as it enforces the vanishing of the order parameter, the chiral condensate, in the high-temperature phase. An analytic function cannot remain constant on any finite interval unless it is constant everywhere. Therefore, the high-temperature phase must be separated from the low-temperature phase by a singularity.

Using symmetry, universality and renormalization group considerations, it can be further argued [74] that due to the SU(3)$\times$ SU(3) chiral symmetry and the U(1)A anomaly in an idealized world with three lightest quarks ($u$, $d$ and $s$) being massless, the finite-temperature transition must be of the first order. By continuity, this must remain true for a range of quark masses near the point $m_{u}=m_{d}=m_{s}=0$. The size and shape of this parameter region, usually drawn on the $m_{s}$ vs $m_{u,d}\equiv m_{u}=m_{d}$ diagram known as “Columbia plot” [70], has been and still remains a subject of lattice QCD studies, with the main problem currently being strong cutoff dependence of the result [75]. All lattice results, however, are in agreement on the conclusion that the region of the first-order transition (at $\mu_{B}=0$) on the Columbia plot does not extend to the physical point given by the real-world values of quark masses $m_{s}$ and $m_{u,d}$. In other words, the finite temperature transition in QCD with physical quark masses is not a first-order transition but an analytic crossover.

3.2 Phase transition at finite baryon density

As we increase the baryon chemical potential $\mu_{B}$ the nature of the finite-temperature transition from hadron gas to quark-gluon plasma can change from the continuous (second-order transition or crossover) to discontinuous (first-order) transition. In the idealized world with massless $u$ and $d$ quarks and sufficiently heavy $s$ quark, where the transition is second order at $\mu_{B}=0$, the change of the order of the transition at finite $\mu_{B}$ would correspond to a tricritical point on the phase diagram found in models of QCD [76, 77, 78, 33]. In the world where $u$ and $d$ quarks have small masses and the second-order transition at $\mu_{B}=0$ is replaced by a crossover, the first-order transition would have to begin at a critical point as shown in the sketch of the QCD phase diagram in Fig. 6. Many models of QCD show such a phenomenon (see, e.g., Refs. [33, 34] for a review). The experimental search for the signatures of the QCD critical point and associated first-order transition is underway and is a main subject of this review.

First-principle lattice QCD calculations provide a powerful method for studying QCD thermodynamics at $\mu_{B}=0$ [79]. However, reliable lattice calculations at finite chemical potentials are impeded by the notorious sign problem (see, e.g., Ref. [80] for a review). Estimates based on various techniques circumventing the sign problem such as reweighting or Taylor expansion, have yielded results which are consistent with the existence of the critical point at $\mu_{B}\gtrsim 200-300$ MeV, while the imaginary $\mu_{B}$ extrapolation does not show a critical point at least until $\mu_{B}\gtrsim 500$ MeV. Taken together at face value, these results disfavor the existence of the critical point for $\mu_{B}\lesssim 200$ MeV. It should be emphasized that some of these results are obtained with lattice spacing in the range where the cutoff dependence is still very strong and makes the extrapolation to continuum limit uncertain. There has been also much effort and progress in addressing the sign problem recently which will have to be reviewed elsewhere. This review will focus on the experimental search for the QCD critical point and the corresponding first order phase transition.

3.3 First-order transition and the critical point

Many physical properties of the critical point can be better understood as a consequence of the fact that this point terminates the coexistence curve in the $T$ vs $\mu_{B}$ plane. On this curve the two coexisting phases are separated by a first-order phase transition and the densities of conserved quantities, such as energy or baryon density, as well as other quantities, such as the chiral condensate, are discontinuous. The discontinuities vanish at the critical point, where the two coexisting phases become one and the same. The increasing similarity between the two phases, and the diminishing thermodynamic barrier between them are at the origin of the large fluctuations at the critical point. These qualitative features are universal, i.e., independent of the physical nature of the underlying system or substance with the critical point.

The discontinuities across the first-order phase transition line reflect the difference in the physics of the phases. Mechanical, thermal and chemical equilibrium on the coexistence line require continuity of pressure $p$, $T$ and $\mu_{B}$ respectively. However, the entropy density, $s$, is higher on the QGP side, i.e., $\Delta s>0$, due to the liberation of color degrees of freedom from confinement. Similarly, the baryon number density $n_{B}$ is higher on the QGP side, i.e., $\Delta n_{B}>0$, because the degrees of freedom carrying the baryon number – quarks – are lighter due to restoration of the chiral symmetry. These discontinuities across the first-order phase transition curve are related by the Clasius-Clapeyron relationship to the slope of the curve: $dT/d\mu_{B}=-\Delta n_{B}/\Delta s$, which follows from equilibrium conditions $\Delta p=\Delta T=\Delta\mu_{B}=0$ and $dp=sdT+n_{B}d\mu_{B}$. The negative slope $dT/d\mu_{B}<0$ of the transition line on the phase diagram is consistent with the physics associated with the discontinuity: $\Delta s>0$ and $\Delta n_{B}>0$ [77].

As we shall discuss in more detail in the next section, a critical point terminating the first-order phase transition is a ubiquitous phenomenon which occurs in physically very different systems such as ferromagnets and liquids. Despite completely different microscopic nature, the collective effects near such critical points are remarkably universal. After a certain remapping the leading singularity of the equation of state near all such critical points (within the same universality class) is the same. This is a consequence of the fact that critical phenomena are caused by the long-range correlations which emerge in the thermal fluctuations of the system. Qualitatively, these fluctuations are related to the transitions between two coexisting phases of the system which become increasingly similar, i.e., close to each other as the first-order transition approaches its terminal (critical) point. In the limit of infinite correlation length such long-range fluctuations are described by a conformal (scale-invariant) field theory defined as the infrared fixed point of a renormalization group evolution (coarse-graining) of the system. The universality of critical phenomena directly follows from the independence of the fixed-point theory on the initial point of the evolution (within a universality class). The scaling and universality of the critical phenomena is the key to our ability to make predictions about the QCD critical point and its signatures in heavy-ion collisions despite having only incomplete information about the QCD equation of state. We shall discuss this in more detail in Section 4.

3.4 High baryon density regime

At still larger baryon chemical potentials the QCD phase diagram is believed to be very rich. Much of the physics at asymptotically large $\mu_{B}$ can be addressed perturbatively, where color-flavor locking determines the ground state of QCD matter [81]. The most interesting regime of intermediate $\mu_{B}$ relevant to neutron stars [3, 4] and, possibly, low-energy heavy-ion collisions is a subject of current theoretical research. Many interesting possible phases and regimes have been proposed, with predictions often depending on the interplay between non-perturbative effects and dynamical assumptions made in calculations. The equation of state of QCD in this regime is, unfortunately, inaccessible to present lattice first-principle calculations due to the sign problem.

It should also be mentioned that gravitational wave astronomy and the recent discovery of neutron star mergers have opened a new exciting avenue for exploring the QCD equation of state in this domain. We shall leave these interesting topics to surveys elsewhere. The present review will focus on the investigation of the QCD phase diagram in the regime readily accessible by heavy-ion collisions.

4 Theory and Phenomenology of the Critical Point

4.1 Critical phenomena

The critical phenomena which we shall focus on occur at an end-point of a first-order transition in a thermodynamic system. The first-order transition corresponds to a situation when a thermodynamic system under given external conditions (such as $T$ and $\mu$, for example) can be in equilibrium in two distinct thermodynamic states. Such a two-phase coexistence can occur only for special values of external parameters, typically, on a manifold of one less dimension than the space of external parameters. E.g., in the $T-\mu$ plane this manifold is a first-order transition line. One of the two states is thermodynamically stable on one side of the first-order phase transition, and the other – on the other side. By adjusting parameters along the phase-coexistence line one could arrive at a special point where the difference between the two coexisting phases disappears. This is a critical point, also known as a second-order phase transition. This point is characterized by critical phenomena which manifest themselves in singular thermodynamic and hydrodynamic properties.

The two most common examples of such critical points are the end point of the liquid-gas coexistence curve and the Curie point in a (uniaxial) ferromagnet. Although the two systems in which these two examples occur are different on a fundamental, microscopic level, the physics near the critical point is remarkably similar on qualitative as well as quantitative level. This observation is the basis of the concept of universality of the second-order phase transitions.

The uniaxial, or Ising, ferromagnet is the simplest of such systems. It can be modeled by a lattice of spins $s_{i}=\pm 1$, or two-state systems, with local (e.g., nearest neighbor) interaction favoring the alignment of spins in the same direction. There are two ground states, with all the spins pointing in one of the two possible directions. The degeneracy is lifted if one applies a magnetic, or ordering, field $h$, which changes the energy of the spins by $h\sum_{i}s_{i}$. The two ordered states are distinguished by the value of the magnetization

[TABLE]

which equals $+1$ or $-1$ depending on the sign of $h$, or more precisely, by its thermal average $\langle M\rangle$. At finite, low enough temperature the ordering persists and $\langle M\rangle$ plays the role of the order parameter which flips sign at $h=0$. The two ordered phases coexist on the line $h=0$ in the $T-h$ plane as shown in Fig. 7

The magnetization $M$ along the coexistence line, $h=0$, decreases with increasing temperature due to thermal fluctuations. At the Curie temperature, $T_{c}$, the magnetization completely vanishes and remains zero for all higher temperatures. The coexistence line (the first-order phase transition) ends at $T=T_{c}$ – the critical point. There is only one phase at and above the Curie point temperature.

Similarly, liquids (e.g., water) coexists with their vapour at given pressure $p$ at the boiling temperature $T$, which defines a line in the $T$ vs $p$ plane. At any of the coexistence points on this line the molecules making up the substance can be arranged in two possible ways, or phases, with the same pressure and temperature. The liquid phase has higher density, but the kinetic energy of molecules in it is smaller due to the molecular attraction. The resulting pressure is the same as that of the gas phase of a lower density.

As the temperature and pressure increase along the coexistence curve the density difference $\Delta\rho\equiv\rho_{\mathrm{liquid}}-\rho_{\mathrm{gas}}$ decreases and vanishes at the critical point, where the coexistence line ends.

In order to enhance the similarity between the two examples of the phase transitions, it is instructive to think of the line of vanishing ordering (magnetic) field, $h=0$, between $T=0$ and $T_{\mathrm{c}}$ on the $T$ vs $h$ plane of an Ising ferromagnet as a coexistence line where two possible arrangements of spins with equal magnitude and opposite sign of magnetization coexist. The discontinuity of the magnetization $\Delta M=2M$ is then the Ising ferromagnet analog of the density discontinuity $\Delta\rho$ in the liquid-gas transition. Even though $\rho$ is not a measure of “order” in the same sense as $M$ is, the similarity is enough to justify common terminology “order parameter” in reference to density $\rho$.

In the following we will review the most important properties of criticality pertinent to the experimental search for the QCD critical point, such as the divergences of the fluctuations and their relation with the correlation length $\xi$. For many other aspects of critical behavior of the phase transitions we refer the reader to standard textbooks and reviews on the subject (e.g., Refs. [82, 83, 84]).

4.2 The role of fluctuations

A useful way to understand the state of a system at a critical point is to think of it as a point where the difference between two coexisting states, or phases, has just vanished. These states coexist along the first-order transition which terminates at the critical point. In the Landau-Ginzburg picture [73] of the transition each of the states corresponds to a minimum of the Landau-Ginzburg free energy (as a function of an order parameter). At any point on the first order transition line the two alternative (or coexisting) states are equally deep and are separated by the Landau-Ginzburg free-energy barrier. The path from one state to the other lies through non-equilibrium states (or more precisely, through inhomogeneous states). At the critical point this barrier shrinks and disappears (the lowest free energy cost, which is proportional to interface tension, vanishes).

The approach to the critical point is characterized by an increase of fluctuations, which is the manifestation of the ease of changing the state of the system from one state to the other. Away from the critical point, on the first order line, these fluctuations are suppressed by the barrier of non-equilibrium states, while on the crossover side there is simply only one state.

More quantitatively, one can describe the fluctuations by the probability distribution $P$ of the quantities characterizing a macroscopic state of the system. For the Ising model – magnetization $M$ is a natural choice, since the coexisting states differ by the value of $M$. For a large, but finite, system, $1\ll N<\infty$, the probability distribution $P(M)$ at finite temperature $T$ and given $h$ is a two-peak function, with the two peaks reaching equal height at the coexistence point $h=0$, as shown schematically in Fig. 8. Above $T_{c}$ the probability distribution $P(M)$ has only one peak. The critical point is characterized by the probability distribution $P$, which has a single peak, with a very flat (i.e., zero curvature) top – the two coexisting states (peaks) have just merged together.

The probability distribution of $M$ is related to the Helmholtz free energy $F(M)$:

[TABLE]

Indeed, for large $N$ the average (equilibrium) value of $M$, $\langle M\rangle$, given by maximum of $P(M)$, satisfies

[TABLE]

A qualitatively useful but, due to fluctuations, quantitatively incorrect expression for $F(M)$ is given by the mean-field (or Landau-Ginzburg) free energy:

[TABLE]

where $b(T)>0$ at all temperatures while $a(T)>0$ for $T>T_{c}$ and $a(T)<0$ for $T<T_{c}$. Thus $F$ has two minima at $T<T_{c}$, and correspondingly, $P(M)$ has two peaks (see Fig. 8).

The critical point is characterized by the vanishing of the second derivative of $F$ at its minimum, which is related to the magnitude of the fluctuations of $M$ around the equilibrium, $\delta M\equiv M-\langle M\rangle$, as well as isothermal magnetic susceptibility $(\partial\langle M\rangle/\partial h)_{T}$:

[TABLE]

Thus the divergence of the susceptibility, i.e., flatness of $F$ as a function of $M$ at the critical point, is directly related to the fact that the fluctuations are large at the critical point (see Fig. 8).

The magnitude of fluctuations is of order $\delta M\sim 1/\sqrt{N}$. The subject of thermodynamic fluctuations is usually peripheral in the study of condensed matter systems because the typical number of degrees of freedom $N$ in such systems is huge (e.g., $10^{24}$). There is one notable exception. Near the critical point the correlation length $\xi$ of fluctuations grows arbitrarily large. The actual measure of the importance of the fluctuations is the ratio of the volume of the system $V$ to the typical correlation volume, i.e., the cube of the correlation length, $\xi^{3}$. Although it is still impractical to fine-tune the system to achieve correlation volumes of macroscopic sizes, in the phenomenon of critical opalescence the relevant size (for light scattering) is only the wavelength of visible light. This length is few orders of magnitude longer than the typical correlation length, but near the critical point the correlation length can achieve such large values, resulting in the observable loss of transparency due to scattering of light on the density fluctuations.

Contrary to condensed matter systems, fluctuations in heavy-ion collision experiments do play a significant role even away from a critical point. These systems are large enough to be treated thermodynamically (or hydrodynamically), and yet not too large ($10^{2-4}$ particles) for the fluctuations to become unimportant. As in condensed matter systems, the fluctuations should be even larger near the critical point.

Another fruitful example of a system where fluctuations are important is provided by the statistical lattice models and lattice field theory simulated on computers, since they require the number of degrees of freedom to be not too large for practical reasons.

4.3 Long-range correlations

It is important to keep in mind that the divergence of fluctuations at the critical point is not a microscopically local phenomenon. Rather, it is a collective phenomenon – the result of the correlation between fluctuations of many degrees of freedom. One can say that it is related to the divergence of the number of the degrees of freedom in a correlation volume $\xi^{3}$, not to the magnitude of the fluctuations of each individual degree of freedom (as is obvious in the case of the Ising model, where $s_{i}=\pm 1$).

This can be quantified by considering the correlation function of an order parameter. To make the discussion unified, we shall call this order parameter ${\sigma}(\bm{x})$. For the ferromagnet ${\sigma}(\bm{x})$ is $M(\bm{x})$ – magnetization of a small patch of the lattice centered around point $\bm{x}$, and for the liquid-gas transition it could be chosen as density $\rho(\bm{x})$.

To describe spatial correlations we consider the probability distribution as a functional, ${\cal P}[{\sigma}(\bm{x})]$, for a spatially varying parameter ${\sigma}(\bm{x})$. To determine the correlation length $\xi$ we measure the thermal expectation value of $\delta{\sigma}(\bm{x})\delta{\sigma}(\bm{y})$, where $\delta{\sigma}={\sigma}-\langle{\sigma}\rangle$:

[TABLE]

where the functional integral is taken over all configurations of the fluctuating field $\sigma$ with weights given by $\mathcal{P}[\sigma]$. The logarithm of the probability density ${\cal P}$ defines the corresponding Ginzburg-Landau free energy functional, $\Omega[{\sigma}]$:

[TABLE]

also known as effective Hamiltonian or effective action, depending on the context, which we shall expand in the gradients of $\sigma$:

[TABLE]

For a given field $\sigma$, the coefficient of the gradient term $(\bm{\nabla}\sigma)^{2}$ will not necessarily have the canonical value $1/2$. We can, however, always redefine the field $\sigma$, multiplying it by an appropriate constant, to reduce the gradient term to the canonical form.

Note that for a uniform field ($\sigma=M$) $\Omega[\sigma]$ is directly related to the free energy $F-hM$ according to Eq. (20).

4.4 Mean-field critical behavior

The Gaussian approximation (also known as saddle-point or mean-field approximation) to the path integral in Eq. (24) amounts to expanding $\Omega[{\sigma}(\bm{x})]$ around its minimum, $\sigma_{0}$, or

[TABLE]

and dropping all terms beyond quadratic. The Gaussian path integral is then easy to take. One finds $\langle\sigma\rangle=\sigma_{0}$ and, for fluctuation

[TABLE]

the correlator

[TABLE]

where the correlation length $\xi$ is related to the curvature of the potential $U(\sigma)$ at its minimum:

[TABLE]

At the critical point, two coexisting minima of $U$ continuously merge into one, and the curvature $U^{\prime\prime}(\sigma_{0})$ vanishes. Thus $\xi=(U^{\prime\prime})^{-1/2}$ diverges. The divergence of the correlation length Eq. (30) is the source of the non-analyticity of all thermodynamic functions at the critical point.

Indeed, let us define the volume integral of the fluctuation of the order parameter field

[TABLE]

This extensive quantity fluctuates around its expectation value $\langle\sigma_{V}\rangle=0$. The (square) magnitude of these fluctuations is measured by the second moment $\langle\sigma_{V}^{2}\rangle$ which is related to the integral of the correlator of $\sigma(\bm{x})$

[TABLE]

This quantity is also extensive, i.e., $\kappa_{2}\sim V$, provided the correlator $\langle\sigma\sigma\rangle$ is local, i.e., $\xi$ is finite, or more precisely, $V\gg\xi^{3}$. In the mean-field approximation

[TABLE]

This result is easy to understand as the zero-momentum value of the propagator of the field $\sigma$, $(q^{2}+m_{\sigma}^{2})^{-1}$. It is also useful to view the divergence $\xi^{2}$ as a result of the divergence of the integral of $1/|\bm{x}|$ being cut off at $|\bm{x}|\sim\xi$ according to Eq. (29). I.e., the divergence of the magnitude of the fluctuations is due to the increase of the range of the correlations, not their local magnitude.

Going beyond the mean-field approximation is a non-trivial task requiring the application of renormalization group (see, e.g., [83]), outside of the scope of this review. This typically results in small corrections to the power laws in Eqs. (29) and (33): $\langle{\sigma}{\sigma}\rangle\sim 1/|\bm{x}-\bm{y}|^{1+\eta}$ for $|\bm{x}-\bm{y}|\ll\xi$ and correspondingly $\chi\sim\xi^{2-\eta}$, with very small $\eta\approx 0.04$. 121212 The smallness can be understood as a consequence of the fact the corrections to Gaussian approximation for $\eta$ appear only at two-loop order. In terms of $\epsilon=4-d$ expansion: $\eta=\mathcal{O}(\epsilon^{2})$. There are other power laws, and corresponding exponents, which differ from the mean-field ($\epsilon=0$) values more significantly, i.e., at one-loop order. Most notably, the correlation length diverges as $\xi\sim|T-T_{c}|^{-\nu}$, where $\nu=1/2$ in mean-field approximation, but, in fact, is closer to $2/3$ in $d=3$ dimensions. This means that the specific heat exponent $\alpha$ in $C_{V}\sim|T-T_{c}|^{-\alpha}$ is close to zero $\alpha=2-d\nu\approx 0.11$ in $d=3$ (note that it is exactly zero in $d=2$ and $d=4$), i.e., the specific heat is only weakly divergent.

These corrections become visible only for very large correlation lengths, and for finite non-static systems such as heavy-ion collisions are negligible compared to other more important effects.

4.5 Higher-order moments and cumulants of the fluctuations

The divergence of the correlation length at the critical point has even stronger effect on the measures of the “shape” of the probability distribution $\mathcal{P}[\sigma]$, as pointed out in Ref. [85]. As the critical point is approached the probability distribution becomes less Gaussian, i.e., flatter at the top if the symmetry is maintained as in the case $h=0$ (see Fig. 8) and, in general, less symmetric around a maximum for a given $h\neq 0$. These non-Gaussian properties can be quantified using higher-order cumulants of fluctuations. Extending the expansion of the potential $U(\sigma)$ to quartic order in $\sigma$

[TABLE]

and performing the path integrals in the saddle point approximation one finds

[TABLE]

where in the last line the subscript ‘$c$’ means ‘cumulant’ or ‘connected’. Diagrammatically this calculation is represented by tree graphs in Fig. 9, where each propagator of the $\sigma$ field (wavy line) contributes a factor $\xi^{2}$.

Going beyond the tree-level (saddle-point) approximation, which is necessary near the critical point due to the importance of fluctuations, requires resummation of infinitely many graphs of all loop orders, which is what the renormalization group achieves. The result is that the couplings “run”, i.e., the dependence of the non-Gaussian cumulants $\kappa_{3}$ and $\kappa_{4}$ on the correlation length $\xi$ appears not only explicitly, but also via the coefficients $\lambda_{3}$ and $\lambda_{4}$ which, similarly to $m_{\sigma}^{2}$, depend on $\xi$. Near the critical point this dependence follows from the invariance of the Landau-Ginzburg free energy $\Omega$ under rescaling $x\to x/b$, $\xi\to\xi/b$, $\sigma\to b^{1/2}\sigma$. One can check that such rescaling leaves the gradient term in Eq. (26) invariant.131313Fluctuations slightly modify the power $1/2$ into $1/2-\eta/4$, but for our purposes it is sufficient to approximate $\eta\approx 0$ as we already noted at the end of previous subsection. The invariance of the term $\int d^{3}\bm{x}U(\sigma)$ requires $m_{\sigma}^{2}\sim\xi^{-2}$, $\lambda_{3}\sim\xi^{-3/2}$ and $\lambda_{4}\sim\xi^{-1}$. The scaling invariance is one of the most important properties of the critical phenomena and is a consequence of the fact that in this regime the correlation length $\xi$ is much larger than any other microscopic scale, and thus becomes the only relevant length scale in the system. Therefore

[TABLE]

where we introduced dimensionless couplings $\tilde{\lambda}_{3}$ and $\tilde{\lambda}_{4}$ via $\lambda_{3}=\tilde{\lambda}_{3}T(T\xi)^{-3/2}$ and $\lambda_{4}=\tilde{\lambda}_{4}(T\xi)^{-1}$. These rescaled couplings are scale-invariant and universal.141414Up to very slow “running” due to small $\eta$ which, as we already noted, we can neglect. Eqs. (36) and (37) demonstrate that the dependence on the correlation length is much stronger in higher-order cumulants, which measure the shape of $\mathcal{P}[\sigma]$, compared to the second order cumulant in Eq. (33), which measures the width of $\mathcal{P}[\sigma]$. Furthermore, the higher-order cumulants in Eqs. (36) and (37) may change sign, depending on $T-T_{c}$ and $h$ [86].

This is a good place to emphasize that the results such as the ones shown in Eqs. (33), (36) and (37) are universal in the sense that the scaling with certain powers of $\xi$ is independent of the microscopic details of the theory. This universality is the key to our ability to search for the QCD critical point even despite having only incomplete knowledge of the QCD equation of state and the microscopic dynamics.

4.6 The cumulants and the equation of state in parametric representation

Cumulants of fluctuations can be determined directly from the equation of state, i.e., the free energy $F(M)$ in Eq. (20). Indeed, identifying

[TABLE]

we can write, using Eq. (23),

[TABLE]

In order to generalize this to higher-order cumulants we will introduce the Gibbs free energy as a Legendre transform of $F(M)$ with respect to $M$:

[TABLE]

Since $\langle M\rangle=-(\partial G/\partial h)_{T}$ we can rewrite Eq. (39) as $\kappa_{2}=-NT(\partial^{2}G/\partial h^{2})_{T}$. This relation can be easily generalized to cumulants of all orders:

[TABLE]

We can thus directly determine the dependence of $\kappa_{n}$’s on $T$ and $h$ from the equation of state $\langle M\rangle(h)$, or from $G(h)$, without knowing the values of $\xi$ and $\lambda_{n}$’s at given $T$ and $h$.

Even for the mean-field free energy in Eq. (22) the explicit expression for $G(h)$ of Eq. (40) or $\langle M\rangle(h)$ looks cumbersome. A tangible expression of the equation of state in this case can be written in the form $h=h(\langle M\rangle)$ where:

[TABLE]

This equation of state is singular at $T=T_{c}$ and $h=0$, since the susceptibility $(\partial M/\partial h)_{T}=1/a$ diverges because $a(T_{c})=0$. The dominant part of this singularity can be described by approximating $a(T)$ and $b(T)$ by the leading term of their Taylor expansion around $T=T_{c}$:

[TABLE]

where $a^{\prime}=da/dT|_{T=T_{c}}$ and $b=b(T_{c})$ and we introduced a “reduced temperature” $r=T-T_{c}$. The most notable feature of the equation of state is invariance under rescaling

[TABLE]

with specific values of $\beta=1/2$ and $\delta=3$, known as mean-field critical exponents. The scaling is observed near Curie critical points in ferromagnets and near liquid-gas critical points. It is characterized by the universal values of $\beta$ and $\delta$ for all systems in the Ising model universality class, which includes uniaxial ferromagnets, liquid-gas and binary fluid transitions.151515This is a ubiquitous critical universality class because it occurs in systems with a one-component (singlet) order parameter. For example, the Heisenberg ferromagnet universality class requires an exact O$(3)$ symmetry and order parameter which is a triplet (a vector). However, the universal values of the critical exponents observed ($\beta\approx 1/3$ and $\delta\approx 5$) deviate significantly from the mean-field values. The role of fluctuations which are responsible for these deviations from the mean-field values can be understood within the renormalization group approach to critical phenomena which is a subject of many classic textbooks (see, e.g., [82, 83]).

Thus the mean-field equation of state, while providing a simple and intuitive description of the phase transition, is not adequate for describing the leading singular behavior of the critical equation of state quantitatively. A quantitative description can be achieved by using the parametric representation [87, 83], where, instead of using the variables $h$, $r$ and $M$ directly, as in Eq. (42), one introduces two auxiliary variables $R$ and $\theta$, roughly parameterizing the “distance” from the critical point and the “angle” relative to the crossover direction (i.e., $h=0$ for $T>T_{c}$) respectively. In terms of $R$ and $\theta$ the variables of the Ising model $M$, $h$ and $r$ are given by

[TABLE]

The correct scaling Eq. (43) is built into this representation and corresponds to $R\to\lambda R$ with $\theta$ being scale invariant.

All the information about the equation of state $M(r,h)$ is in the values of critical exponents $\beta$ and $\delta$ and the universal scaling function $H(\theta)$ which can be calculated order by order in the expansion around dimension $d=4$, i.e., $\varepsilon$ expansion where $\varepsilon=4-d$. To an approximation sufficient for our purposes it is given by

[TABLE]

Note that when scaling exponents are assigned their mean-field values $\beta=1/2$ and $\delta=3$ the mean-field equation of state Eq. (42) emerges from the parametric representation given by Eqs. (45) and (46), with parameters $a^{\prime}=3$ and $b=1$. At $d=4$ (i.e., $\varepsilon=0$) the exponents $\beta$ and $\delta$ are given by their mean-field values. At order $\varepsilon$ the exponents differ from mean-field values 161616More precisely, $\beta=1/2-\varepsilon/6+\mathcal{O}(\varepsilon^{2})$ while $\beta\delta=3/2+\mathcal{O}(\varepsilon^{2})$., while Eq. (46) is correct to order $\varepsilon^{2}$ [83].

One can see from Eq. (45) that parameter $R$ is the measure of the distance from the critical point $r=h=0$. The correlation length diverges as $\xi\sim R^{-\nu}$ when $R\to 0$. The mean-field value of the exponent $\nu$ is $1/2$, as follows from Eq. (30) and $m_{\sigma}^{2}\sim r$. However, the actual value of the exponent $\nu$ (in $d=3$) is significantly different: $\nu\approx 2/3$. One can check that this value (together with $\beta\approx 1/3$ and $\delta\approx 5$) is consistent, via Eq. (41), with the scaling of the cumulants in Eqs. (36) and  (37).

The parameter $\theta$ is a measure of the direction in the $r$ vs $h$ plane (or more precisely, $r^{\beta\delta}$ vs $h$), with $\theta=0$ along the crossover line, $\theta=1$ along $r=0$, and $\theta=\sqrt{3/2}$ along the coexistence (first-order transition) line. The sign of the cumulants depends on the value of the scale-invariant parameter $\theta$. As an example, the corresponding dependence of $\kappa_{4}$ on $r$ and $h$ is shown in Fig. 10. The curves along which $\kappa_{4}$ changes sign correspond to $\theta\approx\pm 0.32$. The negative sign of $\kappa_{4}$ along the crossover line is the reflection of the fact that the distribution is “flatter” than Gaussian when one approaches the critical point, as seen in Fig. 8.

4.7 Mapping to QCD equation of state

The QCD equation of state with a critical point is an essential ingredient for hydrodynamic simulations of the heavy-ion collision fireball evolution at collision energies relevant to the beam-energy scan program. Since first-principle lattice simulations, hindered by the sign problem, cannot yet produce the needed equation of state at nonzero $\mu_{B}$, the following approach has been proposed in Ref. [88]. The available information from lattice (at $\mu_{B}=0$) can be combined with the known universal behavior of the equation of state (pressure vs $T$ and $\mu_{B}$) near the critical point. Of course, this information is not sufficient to completely fix the equation of state. However, it allows to construct a family of equations of state which depends on a relatively small number of parameters (e.g., position of the critical point) while satisfying known constrains. Comparisons of hydrodynamic simulations with such a flexible parametric equation of state and experimental data could help detect and identify the signatures of the QCD critical point.

In order to use the universality of the critical phenomena to describe the leading singularity of the equation of state at the critical point we need to map the coordinates $T$ and $\mu_{B}$ on the QCD phase diagram to the coordinates $r$ and $h$ on the phase diagram of the Ising model (see Fig. 11). 171717 In this section $T$ denotes the temperature in QCD (or, more generally, the temperature in a liquid-gas system under consideration), while $r$ is the reduced temperature of the Ising model $r=T_{\rm Ising}-T_{c,\rm Ising}$, as before. The mapping is linear, i.e.,

[TABLE]

Nonlinearities in the mapping have no effect on the leading singularity of the equation of state. The coefficients $r_{T}$, etc., determine directions of the $r$ and $h$ axes in the $T$ vs $\mu$ plane shown in Fig. 11. In particular, the slope of the coexistence line at the critical point is given by setting $h=0$ in Eq. (47): $(dT/d\mu)_{h=0}=-h_{\mu}/h_{T}$. The slope of the $h$-axis ($r=0$) is non-universal and corresponds to the leading asymmetric correction to scaling [89, 90]. Taking the mixing in Eq. (47) into account one can write the QCD pressure as:

[TABLE]

where $G(r,h)$ is the leading singular contribution to Gibbs free energy of the Ising model. The relation between the QCD pressure (or the pressure in any system in the vicinity of the critical point in the same universality class – that of the Ising model) in Eq. (48) is the consequence of the universality of critical phenomena which follows from renormalization group theory. This relation allows us to make certain predictions about the behavior of QCD near the critical point despite having only incomplete knowledge about the QCD equation of state, thus enabling an experimental search for this phenomenon.

As we discussed in the previous section an explicit expression for $G$ is cumbersome, but for all practical purposes it can be found from Eq. (40) with $F$ given by integrating $h=(\partial F/\partial M)_{r}$ using the equation of state $h(M)$ given by the parametric representation Eqs. (45) and (46). The remaining part of the pressure, $p_{bg}(T,\mu)$, is the ‘background’ which contains less singular (corrections to scaling) and non-singular contributions.

The parameters $r_{T}$, $r_{\mu}$, $h_{T}$ and $h_{\mu}$ are nonuniversal. If the information about the equation of state is available, e.g., from a lattice calculation or beam energy scan experiment, one can attempt to determine these parameters by fitting the parametric equation of state such as given in Eq. (48). The parameters which need to be fitted are the location $T_{c}$ and $\mu_{c}$ of the critical point, the slope $-h_{\mu}/h_{T}$ of the coexistence line ($h=0$), and the slope $-r_{\mu}/r_{T}$ of the $r=0$ line at the critical point. To set the scale, one needs to supply also $r_{T}$ and $h_{\mu}$. The overall scale of the singular part is not an additional independent parameter because of the scaling property of the leading singularity of the equation of state:

[TABLE]

for an arbitrary $\lambda$.

The background part of the pressure, $p_{bg}(T,\mu)$, can be chosen to smoothly match the equation of state at $\mu=0$ known from the lattice [88]. In this form the equation of state can incorporate the information reliably known from lattice QCD calculations as well as the correct leading singular behavior at the critical point, while being flexible enough to accommodate a critical point in the range of $T$ and $\mu$ accessible by heavy-ion collisions.

4.8 Baryon number cumulants near the critical point

In order to understand better the equation of state near the QCD critical point described by Eq. (48) it is helpful to study the behavior of baryon number cumulants. Due to the relationship between the pressure and the partition function of QCD

[TABLE]

the cumulants of the fluctuations of the baryon charge $N$ are related to the derivatives of the pressure:181818In the context of lattice calculations the susceptibilities are often defined as dimensionless quantities, i.e., $\chi_{k}^{\rm lattice}=\partial^{k}(p/T^{4})/\partial(\mu/T)^{k}$.

[TABLE]

where

[TABLE]

Fig. 12 shows a close-up of the QCD phase diagram in Fig. 6 near the critical point. Using the universal equation of state given by the mapping in Eqs. (48) and (47) (where we choose $r_{T}=0$ for simplicity), we illustrate the behavior of susceptibilities $\chi_{k}$. It is instructive to follow $\chi_{k}$ along lines of fixed $T$. Three such lines are shown in Fig. 12 (top row): two isothermal lines traverse the crossover region above the critical point and the behavior of $\chi_{k}$ along these lines is shown in the second row and one isothermal line traverses the first-order coexistence line with the corresponding $\chi_{k}$ shown in the third row.

As we traverse the crossover region (panel (a) in Fig. 12) the density increases continuously with a steeper slope for the case where the isothermal line is closer to the critical point. When we cross the first order line, the baryon density, $n$, jumps, as expected. The baryon number cumulants, or susceptibilities, $\chi_{k}$, being derivatives of the density (see Eq. (51)), will be sensitive to the proximity of the critical point in the crossover region as the change of the density $n$ becomes steeper. This is illustrated in the panels (b) through (d) in Fig.  12, where we show the second to fourth order susceptibilities. We see that, not surprisingly, the steeper increase in the density when traversing the pseudo-critical region closer to the critical point is reflected in larger values of the cumulants. This difference gets more pronounced the higher the order of the susceptibility or cumulant. Furthermore, the sign changes of the various cumulants shown in the contour plots can be easily understood as simply changes in the slope (for $\chi_{2}$), curvature (for $\chi_{3}$) and higher derivatives of the density $n$ in the first column. Finally, when crossing the first-order line (third row) we find that away from the critical line the cumulants are only modestly changed. On the critical line, of course, they are undefined due to a discontinuity. 191919The absence of visible discontinuity in even cumulants in Fig. 12 is a consequence of our simplification $r_{T}=0$.

This simple example qualitatively explains what happens near the critical point as discussed in section 4.5: The higher the order of the cumulant the stronger is its dependence on the correlation length. As we get closer to the critical point, where correlation length diverges, the transition gets sharper and the cumulants also diverge at the critical point.

As we have seen, the high-order cumulants show nontrivial dependence on $T$ and $\mu$ in the crossover region. This observation suggests that the measurement of net-baryon cumulants may also provide an avenue to establish the existence of a cross-over transition at $\mu_{B}=0$, as predicted by lattice QCD [28]. As discussed in [91, 92, 93, 94] in the context of model as well as lattice QCD calculations, a cross-over transition results in negative sixth and eighth order cumulants at the freezeout temperature, $\kappa_{6}/\kappa_{2}<0$ and $\kappa_{8}/\kappa_{2}<0$. Therefore, the measurement of these cumulant ratios could provide experimental evidence that the systems created in high energy heavy ion collisions freeze out close to the cross-over transition.

4.9 Fluctuation cumulants in heavy-ion collisions

The baryon number cumulants, or susceptibilities, are not directly measurable in heavy-ion collision experiments which detect charge particles, leaving neutrons out of the acceptance. However, the fluctuations near the critical point affect fluctuations of charged particles as well as the neutral ones because the coupling of the critical mode is isospin blind. Thus cumulants of the fluctuations of proton number (or net proton number) show a similar pattern near the critical point. In Section 4.12 we shall describe how to relate the critical mode fluctuations with the observable fluctuations of the particle multiplicities.

The experiments also do not scan the phase diagram along fixed $T$ lines as in Fig. 12. The scanning parameter, such as $\sqrt{s_{NN}}$, affects both $T$ and $\mu$ of the freezeout. A typical freezeout trajectory along which $T$ and $\mu$ are varied is shown in Fig. 13 superimposed on the density plot of the quartic cumulant of a critical order parameter, such as, e.g., baryon density. The position of the freezeout point on the curve depends on the collision energy $\sqrt{s_{NN}}$ and can be determined experimentally using the chemical freezeout systematics extracted from the measured particle yields [95]. The systematics provides an estimate of $T$ and $\mu$ of the system at freezeout at given $\sqrt{s_{NN}}$.

As we shall see below in Section 4.12 the critical contribution to the cumulants of the fluctuations of the observed particle multiplicities are proportional to the cumulants of the fluctuations of the critical order parameter [96, 85, 86, 97]. Thus the deviation of the particle multiplicity cumulant from its noncritical baseline will depend on collision energy (or on $\mu_{B}$ along the freezeout curve) as shown qualitatively in Figure 14 .

4.10 Finite time dynamics and limitations on the correlation length in heavy-ion collisions

Phase transitions and associated thermodynamic singularities are properties of static systems in the limit of infinite volume. In a finite system the phase transition singularities are smeared out. The maximum correlation length achievable on a finite system is limited by the system size. In the case of heavy-ion collision this limit is on the order of 10 fm.

However, a more stringent constraint comes from the fact that the system does not spend enough time in the critical region for the correlation length to build up to its equilibrium value [96, 99]. The amount of time that it takes to equilibrate diverges near the critical point

[TABLE]

where $z\approx 3$ is the dynamic critical exponent [100]. Thus, for a system evolving through the critical region with a characteristic time scale $\tau$ the maximal achievable correlation length will be proportional to $\tau^{1/z}$ [96, 99]. The typical evolution time is on the order of $7-10$ fm. Estimating the coefficient of proportionality requires dynamical modeling. The simplest models of this type predict the maximum correlation length on the order of $2-3$ fm [99] (see also Ref. [101] for a more recent discussion of the physics involved in Eq.(53)).

We note that the maximum correlation length due to critical slowing down limits the applicability of finite size scaling arguments [102, 103, 104]. The finite-size scaling describes the behavior of systems so close to criticality that $\xi$ is as large as the system itself. Such arguments do not apply to systems created in heavy-ion collisions because finite-time dynamics limits $\xi$ to be several times smaller ($2-3$ fm) than the size of the system ($7-10$ fm).

A more comprehensive description of the effect of space-time evolution on the fluctuation signatures of the critical point requires embedding the fluctuations into hydrodynamic evolution code and is work in progress. The physics of the critical slowing down and of critical fluctuations are essentially the same physics, described by the evolution of equal-time correlators of the slowest mode (entropy per baryon) towards their equilibrium value. Such an extension of hydrodynamics (called “Hydro+”) incorporating additional critically slow degrees of freedom (correlators of the slow mode) has been introduced in Ref. [105]. The first application of this approach to heavy-ion collisions, albeit in a very simplified setup, has been recently reported in Ref. [106]. Generalization of this approach to non-Gaussian cumulants is still a challenge, but is expected to describe interesting “memory” effects, as in, e.g., Ref.[107, 108, 109].

In order to describe the effects of the critical fluctuations on the observable particle fluctuations and correlations, below we shall consider a simplified model scenario [96, 86, 110], where we parameterize the effects of the evolution near the critical point by the value of the correlation length $\xi$, which determines the cumulants of the fluctuations via scaling relations such as Eqs. (33), (36) and (37).

The fact that the correlation length $\xi$ may not reach large values makes higher-order cumulants important signals of the criticality in heavy-ion collisions because their dependence on $\xi$ is much stronger than that of the quadratic cumulant.

4.11 Experimentally observable fluctuations and correlations

Heavy ion collision experiments do not directly measure the values of the order parameter field or its fluctuations. What quantities should then serve as signatures for the QCD critical point?

Fluctuations in heavy-ion collisions are studied by collecting data on event-by-event basis and looking at event-by-event fluctuations. The data for each event is a set of particle momenta and other quantum numbers such as charge, baryon number, spin, etc. Let us denote the set of variables characterizing a given particle by $A$, such that

[TABLE]

We can think of $A$ being a point in a generalized phase space, which includes also (discrete) variables such as charge $q_{A}$, spin $s_{A}$, etc. and the data set for each event as a collection of such points.

What experiments cannot measure is the position of a particle at the time of freezeout. The critical fluctuations we discussed are related to spatial correlations (see, e.g., Eq. (29)) at that time. Let us consider how these spatial correlation translate into the correlations between momenta of the particles.

Does the correlation length $\xi$ determine the range of correlations in momentum space? The answer to a first approximation is ‘no’ [110]. The range of correlations is set by the thermal distribution of particles within the correlated spatial volume. That range is of order of the average momentum difference of the particles and is determined by temperature (e.g., for nonrelativstic particles of mass $m$, it is given by $\Delta p\sim\sqrt{mT}$). In particular, for protons, this corresponds to a typical range of rapidities of order unity. In principle, spatial correlations can be translated into longitudinal momentum correlations by strong longitudinal expansion, as in Bjorken flow. However, for typical correlation length $\xi\sim 0.5-1$ fm, or even near the critical point $\xi\sim 2-3$ fm, and freezout time of order $\tau\sim 10$ fm, the corresponding Bjorken rapidity correlation range $\Delta\eta\approx\xi/\tau\sim 0.05-0.3$ is significantly smaller than the range, $\Delta y_{\rm corr}=\mathcal{O}(1)$, due to thermal momentum spread.

Though the range of momentum correlations is not sensitive to $\xi$, their magnitude is directly related to the correlation length $\xi$. Indeed, the number of particles correlated with any given one grows, as these are the particles which come from the region of space which grows with $\xi$. Below we shall estimate that effect using a simple model.

First, however, we shall define the quantities in terms of which these momentum correlations can be described. At freezeout the particles move freely and the state of the system can be described by the phase space distribution function $f_{A}(x)$, where we denoted by $A$ the set of momentum and other (conserved) quantum numbers characterizing the particle, as in Eq. (54). The space-time coordinates $x$ of particles at freezeout are not directly measurable. The quantity related to $f_{A}(x)$ but directly measurable in experiment is a single-particle momentum distribution function, which is simply the integral of $f_{A}(x)$ over the fireball at freezeout (over the space-like freezeout hypersurface):

[TABLE]

The event average of this distribution is commonly used in experimental analysis

[TABLE]

The distribution function $f$ fluctuates on event-by-event basis. These fluctuations can be described in terms of the phase-space event-averaged product $\langle\delta f_{A}\delta f_{B}\rangle$, where $\delta f=f-\langle f\rangle$. In the absence of correlations this quantity is only nonzero when $A=B$, in which case it equals $\langle(\delta f)^{2}\rangle=\langle f\rangle$ according to the Poisson distribution.

Again, experiments can observe only the momentum space correlations given by the integral

[TABLE]

The quantity often used in experiments to analyze two-particle correlations is a two-particle momentum space density:

[TABLE]

which counts the (event averaged) number of pairs of particles in an element of 2-particle momentum space $d^{3}\bm{p}_{A}d^{3}\bm{p}_{B}$. Without correlations $\rho_{2}$ factorizes: $\rho_{2}(A,B)\to\rho_{1}(A)\rho_{1}(B)$. Thus correlations are characterized by

[TABLE]

This quantity is related to $\langle\delta n_{A}\delta n_{B}\rangle$ by

[TABLE]

where $\delta_{AB}=\delta_{q_{A},q_{B}}\delta^{3}\left(\frac{\bm{p}_{A}-\bm{p}_{B}}{2\pi}\right)$ and the last term makes sure that we do not count a particle in the same momentum cell twice as if it was a pair.

It is easy to express fluctuations of inclusive quantities, such as, e.g., total charge in terms of the correlators we defined:

[TABLE]

Alternatively, in terms of the correlation function $C_{2}$ [111]:

[TABLE]

Here and in what follows we use a shorthand notation for the integral over the momentum space within the detector’s acceptance window as well as particle quantum numbers (cf. Eq. (54)):

[TABLE]

Acceptance cuts also apply to discrete quantum numbers. For example, if some particles are not detected, such as electrically neutral ones, they are not to be included in the sum in Eq. (63). The relevant example is the case when $Q$ is the baryon number: neutrons carry nonzero baryon number $q_{A}=1$ in this case, but should not be included in $\int_{A}$ if they are not detected.

4.12 Translating critical point correlations into observable

momentum correlations

Clearly, the knowledge of the momentum space correlator (57) or, equivalently, $C_{2}$, would allow for quantitative predictions of fluctuation and correlation measurements. There are many sources contributing to $C_{2}$: memory of initial fluctuations, flow-induced correlations, charge conservation, hadron interactions (resonance decays) and jets, to name just a few [112, 113, 114, 115]. We will not attempt to cover all of them here, but some of the effects will be discussed in Section 7.1 . Our focus is on determining the effect of fluctuations associated with the critical point on this experimentally measurable correlator. In other words, we want to know how fluctuations and spatial correlations of the critical mode $\sigma$ affect fluctuations described by momentum-space correlator $\langle\delta n_{A}\delta n_{B}\rangle$. The main point to keep in mind is that critical point effects are localized to the critical region of the phase diagram, i.e., a relatively small interval of beam energies $\sqrt{s_{NN}}$, compared to other contributions unrelated to critical point, which typically persist across the whole range of $\sqrt{s_{NN}}$ and/or vary monotonously with $\sqrt{s_{NN}}$. This is the rationale put forward in Ref. [78] for the beam energy scan.

For the basic estimate and to illustrate the physics at work we shall consider the idealization of thermal equilibrium, which will allow us to apply the results of the earlier subsections for the fluctuations of $\sigma$. Another idealization we shall use exploits the fact that the typical relevant size of the fireball (transverse radius or Bjorken proper time at freezeout) $R\sim 7-10$ fm is large compared even with the maximum correlation length expected in heavy-ion collisions near the critical point. Under assumption $R\gg\xi$ we can consider the spatial correlations as almost local compared to the smoothly varying (in position space) distribution functions $\langle f\rangle$. We can thus describe the spatial correlator of $\sigma$ as a delta-function normalized by matching the space integral of the correlator (see Eqs. (29), (32) and (33)):

[TABLE]

This can be easily generalized to 3-point or 4-point correlation functions [110].

A simple model [85, 86] to describe the correlations between particles induced by the fluctuations of the $\sigma$ field considers the particles in question in equilibrium in the presence of a slowly varying 202020This requires the momenta of the particles to be large compared to $1/\xi$, which is a reasonable approximation near the critical point background of the field $\sigma$. We should expect the properties of the particles, i.e., the energy-momentum (dispersion) relation $E_{A}({\bm{p}};\sigma)$, to depend on $\sigma$. Then the equilibrium distribution $f_{A}$ of the particles depends on $\sigma$ and, as a result, will fluctuate, contributing to the fluctuation of the distribution function in addition to the trivial statistical fluctuations in a free particle gas:

[TABLE]

For example, we can assume, for simplicity (and the present lack of better knowledge of $E_{A}$ dependence on $\sigma$), that the mass $m_{A}$ of the particles depends on the local value of $\sigma$, which is the usual assumption in the sigma-model. 212121 This approach is inspired by the text-book sigma-model where the $\sigma$-field is an order parameter of the chiral symmetry breaking and in which the masses of hadrons depend on the expectation value of $\sigma$. Near the QCD critical point the field which we denote by $\sigma$ must be understood as a certain combination of scalar quatities such as energy density, baryon density and chiral condensate. In this case

[TABLE]

(plus/minus for bosons/fermions), where we defined a coupling $g_{A}\equiv dm_{A}(\sigma)/d\sigma$ and denoted $\gamma_{A}=(dE_{A}/dm_{A})^{-1}=E_{A}/m_{A}$ - the relativistic gamma-factor for a particle with momentum $\bm{p}_{A}$.

The coupling $g_{A}$ plays an important role in determining the magnitude of the critical fluctuation effects. The existing rough estimates of its value for various particle species are based on an assumption that $\sigma$ couples with strength proportional to the particle mass, as in the sigma-model. If this assumption is correct the strength of the critical point signatures in cumulants of the protons should be much larger than in the similar cumulants for pions [97].

The results obtained in this simple model can be reproduced in a more rigorous diagrammatic approach [116]. The same results can also be obtained in a dynamical model of the evolution of the particle gas coupled to a scalar field $\sigma$, which can also describe effects of finite-time evolution [117].

The fluctuations of $\sigma$ via Eqs. (65) and (64) give rise to an additional contribution to the correlator (at freezeout):

[TABLE]

This correlator is local in coordinate space (on distance scales much larger than $\xi$) but is nonlocal in momentum space. This nonlocality is important for understanding the acceptance dependence of the fluctuation measures.

Using equation (57) we find

[TABLE]

where the first term is the trivial statistical Poisson contribution (see also Eq.(60)).

Following Refs. [116, 110], let us note that the dominant dependence on momentum in $\partial f_{A}/\partial\sigma$, comes from factor $f_{A}$ (see, for example, Eq.(66), taking into account that typically $f_{A}\ll 1$). Correspondingly, the dominant contribution to the integrand in Eq. (68) comes from the factor $f_{A}(x)f_{B}(x)$. The product $f_{A}(x)f_{B}(x)$ is the number of pairs of particles with momenta $A$ and $B$ from the same cell (of size $\xi$) located at $x$. In the boost-invariant Bjorken expansion scenario the local rest frame of the fireball moves with rapidity equal to its Bjorken coordinate $\eta$. Thus, most of the particles with momentum rapidity $y_{A}$ come from the cell at $\eta=y_{A}$ (see Fig. 15). The spread of the $f_{A}(\eta)$ around rapidity $y_{A}$ is of order the typical thermal rapidity $|\eta-y_{A}|\lesssim 1$. Therefore, the product $f_{A}f_{B}$ is strongly suppressed if $|y_{A}-y_{B}|$ exceeds typical thermal rapidity: a pair of particles with larger rapidity separation could not come from the same spatially correlated cell. Thus the range of the correlations in the rapidity space, $\Delta y_{\rm corr}$, is set by thermal rapidity (typically of order

1. of the particles (and not by the correlation length $\xi$, which at freezeout corresponds to very short interval of Bjorken rapidity $\Delta\eta\sim\xi/\tau\sim 0.1$) [110].

Using Eq. (61) one can also obtain the corresponding cumulant of a given charge $Q$:

[TABLE]

The value of the cumulant depends on the experimental acceptance (the momentum space region over which $\int_{A}$ and $\int_{B}$ are evaluated). This dependence is directly related to the range of the correlator $\langle\delta n_{A}\delta n_{B}\rangle$. In the Bjorken boost-invariant expansion scenario the first term on the RHS of Eq. (69), i.e., the Poisson term, is simply proportional to the volume of integration, $\int dy_{A}=\Delta y_{\rm acc}$, since the integrand is rapidity independent. The rapidity dependence of the second, correlation term in Eq. (69) is different in two regimes: small and large acceptance. For small acceptance, $\Delta y_{\rm acc}\ll\Delta y_{\rm corr}\sim 1$, each rapidity integral $\int_{A}$ is proportional to $\Delta y_{\rm acc}$ and the correlation term is proportional to $\Delta y_{\rm acc}^{2}$. In the opposite regime222222In practice, this regime is limited by the total rapidity width of the Bjorken plateau, which depends on the collision energy., $\Delta y_{\rm acc}\gg\Delta y_{\rm corr}$, the double integral $\int_{A}\int_{B}$ has support only for $|y_{A}-y_{B}|\lesssim\Delta y_{\rm corr}$ (the diagonal strip in $y_{A}-y_{B}$ space), and thus the cumulant scales as $\Delta y_{\rm acc}$ [110].

It should be kept in mind that $\int_{A}$, in addition to momentum integration, also contains summation over other quantum numbers of the particles, including $q_{A}$ (see Eq. (63)). Different particles can carry different quanta $q_{A}$ of the given charge $Q$. For example, if $Q$ is net baryon charge – protons and antiprotons have $q_{A}=\pm 1$. If their average numbers and couplings to the sigma mode, i.e., $\partial f_{A}/\partial\sigma$, are close to each other (e.g., at large $\sqrt{s_{NN}}$, when baryon/antibaryon asymmetry is small), the positive and negative contributions to $\int_{A}q_{A}(\partial f_{A}/\partial\sigma)$ will tend to cancel each other. Therefore, the critical point contribution to cumulants of positive-definite quantities (e.g., number of protons in acceptance) are typically larger than to the cumulants of such quantities as net-charge or net-baryon number [118, 97].

4.13 Higher-order cumulants

This analysis can be generalized to higher-order cumulants [85, 110]. For example,

[TABLE]

[TABLE]

One can also generalize this analysis to correlations between different charges and species [97].

Comparing to Eqs. (35) one can see that the magnitude of the critical contribution is proportional to the cumulants of the order parameter $\sigma$. I.e., critical fluctuations of the order parameter are “imprinted” on the fluctuations of the observed particles. One can view the critical contribution diagrammatically as shown in Fig.16, which also illustrates that it contributes to the irreducible correlation function.

As already noted, the cumulants contain non-critical contributions, such as the trivial statistical contribution or Poisson term. Some of these non-critical contributions can be removed by subtracting the trivial statistical contribution. Even better, by subtracting contributions from all lower order correlations, one can construct a factorial cumulant which represents irreducible correlations [110], as we shall discuss in Section 7.1.

Similarly to the quadratic cumulant in Eq. (69), the critical contribution to the higher order cumulants for small $\Delta y_{\rm acc}\ll\Delta y_{\rm corr}$ scale as $(\Delta y_{\rm acc})^{k}$, where $k$ is the order of the cumulant $\kappa_{k}$ and scale linearly with $\Delta y_{\rm acc}$ in the (idealized) limit of $\Delta y_{\rm acc}\gg\Delta y_{\rm corr}$. Since in realistic heavy-ion collisions $\Delta y_{\rm acc}\lesssim\Delta y_{\rm corr}\sim 1$, one should be able to observe the beginning of the transition from the fast $\Delta y_{\rm acc}^{k}$ rise to a slower $\Delta y_{\rm acc}^{1}$ rise. To which extent this is possible with the presently available data will be discussed in Section 7.1.4. It is often convenient to consider the ratio of cumulants such as $\kappa_{k}/\kappa_{2}$, or $\omega_{k}=\kappa_{k}/N$. In this case the $\Delta y_{\rm acc}^{k-1}$ growth would begin to saturate at a constant value at larger $\Delta y_{\rm acc}$. This behavior can be used to determine the rapidity range of an observed correlation. Since the rapidity range is different for the critical correlations ($\Delta y_{\rm corr}\sim 1$) and for, e.g., initial state correlations ($\Delta y_{\rm corr}\gg 1$), the range of correlations could be used to separate critical point signatures from certain background effects.

From the experimental perspective it is important to note [85] that, even though the $\xi$ dependence of $\kappa_{4}$ is stronger than that of $\kappa_{3}$ or $\kappa_{2}$, its measurement involves subtraction of two contributions, such as $\langle(\delta N)^{4}\rangle-3\langle(\delta N)^{2}\rangle^{2}$, each of which is order $N$ times larger than their difference ($\mathcal{O}(N^{2})$ vs $\mathcal{O}(N)$), which typically dilutes the signal-to-noise ratio in experimental measurement. This issue becomes more problematic as the order of the cumulant grows (e.g., $\langle(\delta N)^{6}\rangle_{c}=\mathcal{O}(N)$ is obtained by cancellation of contributions of order $\mathcal{O}(N^{3})$.) Therefore, even though the higher-order cumulants are even more sensitive to the divergence of the correlation length and equations such as (71) can be generalized to higher-order cumulants, the measurement of such cumulants would become impractical at too high orders.

Our analysis in this section was focused on integral measures, such as fluctuation cumulants in Eq. (69). It is clear, however, that more information is contained in the correlator, such as in Eq. (68). Some of this information, such as the range of the correlations, can be extracted by studying acceptance dependence, but some information is lost during the integration. The correlator in Eq. (68), and similar higher-point correlators $\langle\delta n_{A}\delta n_{B}\delta n_{C}\ldots\rangle$ are, of course, directly measurable, albeit with more difficulty.

It is important to note again that other sources may and do contribute to fluctuation measures: remnants of initial fluctuations, flow, jets – to name just a few obvious contributors. A quantitative study of these effects may be necessary to unambiguously identify the critical point signal. This serves to emphasize that the energy scan of the QCD phase diagram is needed to separate such background contributions from the genuine critical point effect, the latter being non-monotonous function of the initial collision energy $\sqrt{s_{NN}}$ [78, 96] (or accepted particle rapidity [119]) as the critical point is approached and then passed. The fact that non-Gaussian cumulants have stronger dependence on $\xi$ than, e.g., quadratic moments, makes those higher-order cumulants more sensitive signatures of the critical point [85].

5 Theory and Phenomenology of Anomalous Chiral Transport

5.1 The many-body physics of the chiral anomaly

Macroscopic properties of matter and symmetries of the microscopic dynamics are often in deep and direct correspondence. For example, if a physical theory is invariant under space and time translation, then energy and momentum must be conserved in any individual interaction process. Correspondingly, this mandates that the many-body system of this theory must also conserve energy momentum at the macroscopic level. The hydrodynamic equations for the energy-momentum tensor are the direct consequences of these conservation laws. Similarly the phase invariance of an underlying theory leads to the conserved charge, which in turn gives rise to the current conservation equation in hydrodynamics. These observations raise an interesting question with regards to the chiral anomaly discussed in Section 2. In a way, the axial symmetry is like a “semi-symmetry”, in that it is classically preserved but quantum mechanically broken. One naturally wonders: how would the microscopic chiral anomaly manifest itself in the macroscopic properties of chiral matter, i.e a many-body system containing chiral fermions? In the past decade, significant progress has been achieved in understanding this question, strongly driven by motivations and efforts from multiple branches of physics, notably from heavy ion collisions as well as from topological semimetals (see recent reviews in e.g. [120, 121, 122, 123]). The short answer is that the anomaly will induce a number of novel anomalous chiral transport effects, such as the Chiral Magnetic Effect (CME), Chiral Vortical Effect (CVE) and Chiral Magnetic Wave (CMW). The many-body theoretical frameworks for the quantitative description of these effects are also under active investigation. In the following, we briefly discuss these developments.

5.1.1 Chiral Magnetic and Vortical Effects

Let us consider a certain type of chiral materials which could be, e.g., the quark-gluon plasma (QGP) with (approximately) chiral quarks or the Dirac/Weyl semimetals with chiral quasi-particles. The macroscopic transport properties, which typically involve collective motion (i.e. momentum correlations) of underlying particles, become markedly different from normal materials because of the nontrivial interplay between the particles’ momentum and spin due to their chirality and the anomaly. In particular, certain transport processes that are forbidden in normal environment become possible (and necessary) in such chiral materials. As a typical example of a normal transport process, let us consider electric conductivity in usual materials: an electric current is generated in response to an applied electric field. In those materials, no current could be generated by applying a constant magnetic field. But in the chiral materials discussed here, an electric current can be generated in response to an applied magnetic field. This is the so-called Chiral Magnetic Effect (CME) [124, 125], the most widely studied example among a class of anomalous transport effects in chiral materials. Specifically the CME predicts the following relation:

[TABLE]

where $\mathbf{J}_{e}$ is the electric current, $\mathbf{J}$ is the underlying fermion vector current 232323Note that the vector current associated with a given type of fermion field is $J^{\mu}=\bar{\psi}\gamma^{\mu}\psi$, while the electric current or baryonic current associated with it can be constructed from this vector current based on the fermion’s electric or baryonic charge accordingly, e.g. $J^{\mu}_{e}=(qe)J^{\mu}$. and $\mathbf{B}$ is the magnetic field. The coefficient $\sigma_{5}$ represents the chiral magnetic conductivity for each species of Dirac fermions with electric charge $qe$. It is expressed in terms of the anomaly coefficient $C_{A}=1/(2\pi^{2})$ and the chiral chemical potential $\mu_{5}$ which quantifies the imbalance in the densities of right-handed (RH) and left-handed (LH) fermions (see Eq. (8)). For QCD with $N_{c}$ colors and $N_{f}$ light quarks of charges $Q_{f}e$, one needs to sum over the quark colors and flavors, so that $(qe)^{2}\to\left(N_{c}\sum_{f=1,...,N_{f}}Q_{f}^{2}e^{2}\right)$.

In order to provide an intuitive understanding of the CME, let us consider a system of chiral fermions with positive charge, as illustrated in Fig. 17. The magnetic field aligns the fermions’ spin preferentially along its direction, i.e. $\langle\mathbf{s}\rangle\propto(qe)\mathbf{B}$. At the same time the momentum and spin of a chiral fermion is correlated so that the momenta of RH particles are preferentially parallel and those of LH particles anti-parallel to the direction of the magnetic field, respectively. Therefore, if the system has a chirality imbalance, say we have more RH than LH particles ($\mu_{5}>0$), the net momentum (and spin) would be parallel to the magnetic field, $\langle\mathbf{p}\rangle\propto\mu_{5}\langle\mathbf{s}\rangle$. Since the particles all carry a positive charge, this would result in an electric vector current along the magnetic field, $\mathbf{J}\propto\langle\mathbf{p}\rangle\propto\mu_{5}(qe)\mathbf{B}$.

As one can see from this intuitive illustration of CME, an essential ingredient for the CME transport current to occur, is the correlation between spin and momentum for a massless chiral fermion of specific RH or LH chirality. A finite fermion mass, either found in the Lagrangian or dynamically generated by spontaneous symmetry breaking, as in the case of the QCD vacuum, would spoil this correlation. More precisely, in a massive fermion theory the chirality imbalance will relax to zero at a rate proportional to mass squared (see, e.g., Refs. [126, 127]). In order for chirality to be maintained, its relaxation rate needs to be considerably smaller than the thermal kinetic relaxation rate, which would require the mass to be small compared with the relevant matter scale like temperature. Therefore, while the effect of a small current quark mass on CME is negligible, the effect of large constituent mass due to chiral symmetry breaking is not. For this reason, the chiral symmetry restoration is important for observation of the CME.

It is instructive to contrast this phenomenon with the usual conductivity driven by an electric field $\mathbf{J}_{e}=\sigma_{e}\mathbf{E}$. Contrary to the electric conductivity $\sigma_{e}$ which depends on the properties of the material such as mean free path etc, the chiral magnetic conductivity $\sigma_{5}$ is “universal” in that it only depends on the chiral anomaly coefficient $C_{A}$, $\mu_{5}$ and particle charge, and is not affected by the interaction among the particles [125]. This reflects the fact that CME is essentially the macroscopic manifestation of the fundamental quantum anomaly in a many-body setting. Furthermore, the chiral magnetic conductivity, $\sigma_{5}$, is ${\hat{\mathcal{P}}}$-odd and ${\hat{\mathcal{C}}}{\hat{\mathcal{P}}}$-odd, due to the ${\hat{\mathcal{P}}}$ and ${\hat{\mathcal{C}}}{\hat{\mathcal{P}}}$ transformation properties of the chiral chemical potential $\mu_{5}$. This explains why the CME may occur only in chiral matter with macroscopic chirality (i.e. nonzero $\mu_{5}$). Finally the $\sigma_{5}$ is time-reversal even [128] which implies the non-dissipative nature of the CME transport current. In contrast, the electric conductivity $\sigma_{e}$ is time-reversal odd, thus involving dissipative transport.

Analogous to the anomaly relations of Eq. (13), the CME-current may be written in terms of the RH/LH sectors:

[TABLE]

Here $\mu_{R/L}$ are the chemical potentials for RH/LH fermions, respectively, which satisfy $\mu_{R}-\mu_{L}=2\mu_{5}$. The above RH/LH formulation also allows one to recognize the so-called Chiral Separation Effect (CSE) [129, 130] in a straightforward way. By noting that $\mathbf{J}_{5}=\mathbf{J}_{R}-\mathbf{J}_{L}$, one obtains from (73):

[TABLE]

where $\mu_{V}=(\mu_{R}+\mu_{L})/2$. It predicts the generation of an axial current along the external magnetic field when there is a nonzero vector charge density (i.e. with its corresponding chemical potential $\mu_{V}\neq 0$).

In addition to a magnetic field $\mathbf{B}$, the rotation of the fluid could also induce anomalous transport in chiral matter. This can be understood by the analogy between magnetic field and fluid rotation, such as the similarity between Lorentz force and Coriolis force in classical mechanics and that between magnetic flux and orbital angular momentum in quantum mechanics. The linear motion of a fluid is described by a velocity field $\mathbf{v}$ in analogy to the vector field $\mathbf{A}$ in electromagnetism. The rotational motion of a fluid, on the other hand, is described by a vorticity field $\bm{\omega}=\frac{1}{2}\nabla\times\mathbf{v}$, which plays the role similar to magnetic field $\mathbf{B}$. In the covariant formulation of relativistic fluid dynamics, one has correspondingly the velocity vector $u_{\mu}$ and the vorticity tensor $\omega_{\mu\nu}=\frac{1}{2}(\partial_{\nu}u_{\mu}-\partial_{\mu}u_{\nu})$. The counterpart of the CME in chiral matter under fluid rotation, is the so-called Chiral Vortical Effect (CVE) [131, 132, 133, 134]. The CVE predicts the generation of the underlying fermion vector current in response to the fluid rotation:

[TABLE]

In comparison to the CME in Eq. (72), the anomaly coefficient $C_{A}$ and the $\mu_{5}$ are common to both, while the vorticity $\bm{\omega}$ replaces the magnetic field $\mathbf{B}$ as the driving “force”. Finally, let us mention in passing that there are additional transport effects such as the Chiral Electric Separation Effect (CESE) [135, 136], characterizing the generation of axial current in chiral matter in response to an external electric field. More detailed discussions may be found in Refs. [120, 122, 137].

5.1.2 Chiral collective excitations

Collective excitations are important “degrees of freedom” in many-body systems. They often emerge in the long-time, large-distance behavior of a system and play key roles in its transport properties. A prime example is the sound wave, a gapless excitation arising from the coupled evolution between the fluctuations of pressure and energy density.

A system of chiral fermions in the presence of an external magnetic field $\mathbf{B}$ exhibits new collective excitations as a result of the aforementioned anomalous transport effects. Physically this is not hard to imagine: Consider such a system with a neutral background (i.e. without any net vector or axial density for the underlying fermions) and suppose an axial density fluctuation occurs somewhere. This will generate a vector current via the CME, which transports vector charges, and thus charge fluctuations, from the initial location along the direction of the magnetic field $\mathbf{B}$. These charge fluctuations, in turn trigger an axial current via the CSE, which transports axial charge fluctuations further away along the magnetic field direction, which again induce further vector charge fluctuation via the CME and so on. Consequently, as a result of the anomalous transport in the presence of a magnetic field, the fluctuations of vector and axial charges get entangled and propagate along the $\mathbf{B}$ field direction.

In order to describe this phenomenon mathematically, it is convenient to use the chiral basis i.e. to examine the fluctuations in the RH and LH sectors. We consider an equilibrated neutral background chiral matter under the presence of a constant magnetic field $\mathbf{B}$ and for simplicity assume no electric field. In this case, the small charge fluctuations of the RH and LH sectors satisfy the following continuity equations as well as the CME relations:

[TABLE]

where the label $\chi=+1$ for RH sector and $\chi=-1$ for LH sector. Furthermore the fluctuations of charge density and its corresponding chemical potential can be connected via the susceptibilities, i.e. $(\delta J^{0}_{\chi})=c_{\chi}\,(\delta\mu_{\chi})$ with $c_{\chi}$ the thermal susceptibility evaluated from the background matter. Combining all these together, we arrive at the following equation for the chiral density fluctuations:

[TABLE]

which takes the form of a wave equation. Here $v_{B}=\frac{C_{A}(qe)|\mathbf{B}|}{2c_{\chi}}$ and $\hat{\mathbf{B}}=\mathbf{B}/|\mathbf{B}|$ is the unit vector along the magnetic field direction. Decomposing the density fluctuation $(\delta J^{0}_{\chi})$ into its Fourier modes $(\delta J^{0}_{\chi})\sim e^{-(i\omega t-\mathbf{k}\cdot\mathbf{x})}$ with frequency $\omega$ and wave-vector $k\hat{\mathbf{B}}$, we find the dispersion relation

[TABLE]

where the $v_{B}$ plays the role of the wave-speed. This is the Chiral Magnetic Wave (CMW), a gapless collective excitation in chiral matter [138, 139]. It has two independent modes (with $\chi=\pm 1$ respectively): the RH wave that propagates in parallel to the magnetic field direction $\hat{\mathbf{B}}$ while the LH wave propagates in the opposite direction to $\hat{\mathbf{B}}$. Just like sound waves that can transport energy and momentum in the usual medium, the CMW can transport vector and axial charges in chiral matter. Different from the CME for which a chirality imbalance (i.e. nonzero axial density) is necessary, the existence of CMW does not require any particular background density and can be simply triggered by density fluctuations in chiral matter as long as there is a magnetic field.

The above analysis can be further generalized into the situation of collective modes on top of any non-neutral background matter under the presence of both electric and magnetic fields [135]. Furthermore it is natural to expect that collective excitations similar to the CMW should be present also under fluid rotation. Indeed, chiral collective modes for RH and LH density fluctuations, referred to as the Chiral Vortical Wave (CVW), were identified in [140] and further analyzed in [141, 142, 143].

5.1.3 Fluid dynamics with chiral anomaly

Fluid dynamics is a theoretical framework that provides a universal effective description of the long wavelength behavior of many-body systems [144], and it is successfully applied to describe a large variety of physical systems, ranging from the smallest fluid droplet (e.g. the QGP on a femto scale) to the largest fluid (e.g. the cosmic evolution). In the large wave length limit the system may be considered to be close to local thermal equilibrium and its behavior is dictated by conservation laws such as energy-momentum and charge conservation. If a system is in strict local equilibrium, it is governed by ideal fluid dynamics. Deviations from local thermal equilibrium can be systematically incorporated via an order-by-order “gradient expansion” which leads to viscous fluid dynamics.

In order to illustrate the gradient expansion let us consider a charged current $J^{\mu}$ (for a single degree of freedom) in fluid dynamics. For such a conserved current, the fluid dynamics equation is simply the continuity equation, $\partial_{\mu}J^{\mu}=0$ where the current is expressed in terms of fluid dynamical variables such as velocity field $u^{\mu}(x)$, charge density field $n(x)$, temperature field $T(x)$, etc. Such an expression for $J^{\mu}$ is the so-called constituent relation. For ideal fluid dynamics (0th order in gradient expansion), one has $J^{\mu,(0)}=nu^{\mu}$ which simply says that the charges move along with the local fluid flow. However a fluid cell cannot be totally isolated from its surroundings and the behavior of one cell will be affected by neighboring cells. Such a “gradient” in the fluid leads to the viscous transport currents. Diffusion is a prime example, which simply says that if the charge density is higher at one location compared to other, nearby locations, then the charges will develop a diffusion current and spread out from the dense spot to the dilute spots. This is a typical dissipative process that increases the total entropy of the system and tends to bring the system toward global thermal equilibrium. Turning on an external electric field has the equivalent effect of creating a gradient in the charge chemical potential along the electric field direction. As it turns out, the second law of thermodynamics (i.e. entropy must not decrease with time) completely fixes the possible form of such leading-order viscous corrections, denoted by $J^{\mu,(1)}$ and known as Navier-Stokes current, as follows:

[TABLE]

In the $J_{viscous}^{\mu,(1)}$, the first term is the diffusion current and the second term is the electric conducting current. Here, T is temperature, $\mu$ is the chemical potential corresponding to the charge density $n$ in thermal equilibrium, $q$ is the electric charge, and the projection operator $\Delta^{\mu\nu}=(g^{\mu\nu}-u^{\mu}u^{\nu})$ with the usual Minkowski metric tensor $g^{\mu\nu}=Diag(+,-,-,-)$. The electric conductivity, $\sigma$, is a diffusion coefficient which is a characteristic transport property of a given material.

We now turn to the discussion of a chiral fluid (i.e. a fluid containing chiral fermions), with a chiral current $J^{\mu}_{\chi}$ (where $\chi=\pm 1$ is a chirality label) that could be from either the right-handed (RH, with $\chi=+1$) sector or the left-handed (LH, with $\chi=-1$) sector. Such a chiral current is associated with the underlying “semi-symmetry” and subject to the chiral anomaly. As a consequence the continuity equation needs to be modified in order to account for the anomaly relation Eq. (13). Furthermore, the CME and CVE discussed before would hint at additional anomalous transport currents that are absent in usual fluid dynamics. As found in [132], the second law of thermodynamics is still able to fix completely the first order gradient expansion terms of the chiral current and in fact mandates the presence of anomalous transport currents in the form of CME and CVE. Such fluid dynamics which incorporates the chiral anomaly, up to the first order gradient expansion, can be summarized as follows:

[TABLE]

Here $E^{\mu}=F^{\mu\nu}u_{\nu}$ and $B^{\mu}=\frac{1}{2}\epsilon^{\mu\nu\alpha\beta}u_{\nu}F_{\alpha\beta}$ are external electromagnetic fields in the fluid local rest frame. The covariant fluid vorticity is given by $\omega^{\mu}=\frac{1}{2}\epsilon^{\mu\nu\alpha\beta}u_{\nu}\partial_{\alpha}u_{\beta}$. While $0^{th}$ and $1^{st}$ order currents, $J^{\mu,(0)}$ and $J^{\mu,(1)}_{viscous}$, remain the same as in usual fluid dynamics of Eq. (80), we find two new terms in $J^{\mu,(1)}_{anomalous}$ which correspond precisely to the CME and CVE respectively. The two coefficients $\xi_{B}$ and $\xi_{\omega}$ are entirely fixed by chiral anomaly and thermodynamics, and (in appropriate frame) reproduce exactly the coefficients in Eqs. (72) and (75). Different from the viscous currents, these anomalous currents are universal for various chiral fluids. Therefore, this new or extended type of fluid dynamics incorporates macroscopic quantum transport currents arising from microscopic anomaly and distinguishes the “left” from the “right”. It also provides a useful framework for the quantitative modeling of anomalous chiral transport in heavy ion collisions [145, 146, 147, 148, 149, 150], to be discussed later in Sec. 7.

5.1.4 Chiral kinetic theory

While the preceding discussions address the macroscopic effects of the anomaly for chiral matter in or near thermal equilibrium, there remains the question about the behavior of such system in an out-of-equilibrium situation. The natural framework to address this situation is kinetic theory based on transport equations for the phase space distribution function of such a system. Different from classical kinetic theory, a proper description of the chiral fermions must account for intrinsic quantum and relativistic effects. During the last several years, such a kinetic chiral kinetic theory has been developed using a variety of approaches, see e.g. [151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161]. Also, several phenomenological attempts to study out-of-equilibrium anomalous chiral transport have been proposed [63, 162, 163, 164]. In the following we will outline the key elements of the chiral kinetic theory.

Let us start by considering the classical kinetic theory, which describes the time evolution of the classical phase space distribution $f^{(c)}(t,{\mathbf{x}},{\mathbf{p}})$ for a given type of particles. Various physical quantities of interest (e.g. physical density and current, energy-momentum tensor, etc) may be constructed by properly integrating the distribution. The corresponding transport equation (under the presence of external electromagnetic fields $\mathbf{E}$ and $\mathbf{B}$) reads [165]:

[TABLE]

In the above, $q$ is the charge of the particle and $E_{\mathbf{p}}$ is the single particle energy for given momentum. Notably, the second line (84) gives the classical equation of motion: individual particles simply follow classical trajectories in between sequential collisions with other particles. The $\mathcal{C}[f^{(c)}]$ schematically represents the collision term, and its precise form would depend on the details of scattering processes under consideration (see e.g. [165]).

Let us next consider the kinetic description of chiral fermions. In this case, the spin degrees of freedom cannot be ignored and thus quantum mechanical effects have to be accounted for. For example, a spin-$\frac{1}{2}$ charged fermion will have a nonzero magnetic moment which interacts with an external magnetic field, in addition to the classical Coulomb and Lorentz forces. More importantly, a fermion with fixed chirality has its spin and momentum states “locked up”, for example a RH fermion has its spin and momentum always aligned. Consequently a change of its momentum (e.g. when the particle moves in an external magnetic field) will necessarily be accompanied by the appropriate change in its spin direction and vice versa. To properly treat the spin degrees of freedom and thus chiral fermions in kinetic theory one therefore has to adopt a quantum transport theory such as the Wigner function formalism discussed in  [166, 167]. It turns out that a semi-classical treatment which takes into account quantum effects up to the first order in $\hbar$ appears to be sufficient for describing the anomaly-induced transport [151, 152, 153]. In such a chiral kinetic theory (CKT), the corresponding transport equation for the semi-classical phase space distribution $f^{(q)}$ can be written as:

[TABLE]

with the Jacobian factor $\sqrt{G}=\left(1+\hbar q\mathbf{b}_{\chi}\cdot\mathbf{B}\right)$ and $\widetilde{\mathbf{E}}=\mathbf{E}-\frac{1}{q}\triangledown_{\mathbf{x}}E_{\mathbf{p}}$. The single particle energy is given by $E_{\mathbf{p}}=|\mathbf{p}|\left(1-\hbar q\mathbf{B}\cdot\bf{b}_{\chi}\right)$ and the corresponding velocity is $\widetilde{\mathbf{v}}=\frac{\partial E_{\mathbf{p}}}{\partial\mathbf{p}}=\hat{\mathbf{p}}\left(1+2\hbar q\mathbf{B}\cdot\mathbf{b}_{\chi}\right)-\hbar qb_{\chi}\mathbf{B}$ with $\hat{\mathbf{p}}=\mathbf{p}/|\mathbf{p}|$. Here the factor $\mathbf{b}_{\chi}=\chi\frac{\mathbf{p}}{2|\mathbf{p}|^{3}}$ is the Berry curvature, with $\chi=\pm 1$ for RH/LH fermions respectively. The collision term also needs to be carefully modified in an appropriate way (see, e.g., Ref. [155]).

In comparison with the classical kinetic theory, in the chiral kinetic theory the order-$\hat{O}(\hbar)$ level quantum effects appear in a number of places. First of all, the single particle energy $E_{\mathbf{p}}$ receives a shift by an order-$\hat{O}(\hbar)$ level magnetization energy which accounts for the interaction between the particle’s magnetic moment and the magnetic field. This also results in a quantum correction to the corresponding particle velocity $\widetilde{\mathbf{v}}$. Further quantum corrections at the same order in ${\mathbf{G}}_{\mathbf{x}}$ and ${\mathbf{G}}_{\mathbf{p}}$ arise from the so-called Berry phase, which originates from the adiabatic progression of the particle spin state in order to “catch up” with the change in momentum. Additionally, the phase space integration measure gets modified by the factor $\sqrt{G}$. Such a CKT formalism has been shown to correctly reproduce the chiral anomaly relation as well as the CME and CVE transport currents in or near the thermal equilibrium [151, 152, 153]. Also, it allows the study of such anomalous transport even far from equilibrium, see e.g. recent examples in [168, 163, 164]. Currently the exploration of the transport theory for chiral fermions is under intensive development, aiming to fully develop its conceptual foundation as well as various applications.

5.2 Anomalous chiral transport in heavy ion collisions

As discussed above, the anomalous chiral transport effects (such as the CME, CVE and CMW) stem from highly nontrivial interplay between chirality, topology and anomaly, which all concern fundamental aspects of QCD. A conclusive observation of them in heavy ion collisions would be of paramount significance, providing the first tantalizing experimental verification for the high-temperature restoration of the spontaneously broken chiral symmetry of QCD as well as a direct manifestation of the elusive gluonic topological fluctuations. In the following let us discuss the necessary conditions for these effects to occur in heavy ion collisions.

As previously discussed, both CME and CVE would occur only under the presence of macroscopic chirality, i.e. requiring a nonzero chiral chemical potential $\mu_{5}$. But how does such a situation occur in the quark-gluon plasma created by heavy ion collisions? The answer lies in the gluonic topological fluctuations. These fluctuations are strong due to nonpertubative gauge coupling and they induce the chirality imbalance of the light quarks via Eq. (15), therefore creating a sizable nonzero $\mu_{5}$ (or more precisely $\sqrt{\mu_{5}^{\,2}}$) in the quark-gluon plasma on an event-by-event basis. The CMW, on the other hand, does not require a nonzero chirality imbalance to begin with. Given the initial conditions, the propagation of CMW leads to specific charge distributions at the end, to be discussed later.

Furthermore the necessary “driving forces” for these effects, namely the fluid vorticity and magnetic field, are available in a typical off-central collision, as illustrated in Fig. 18 (left). In the reaction plane spanned by the impact parameter direction $\hat{x}$ and beam direction $\hat{z}$ the two incident nuclei have their center-of-mass separated by impact parameter $b$, see Fig. 18 (right). They both carry a longitudinal momentum $P_{z}=\pm A\sqrt{s_{NN}}/2$ (with opposite direction) where $A$ is the number of nucleons in each nucleus and $\sqrt{s_{NN}}$ is the center-of-mass energy per nucleon pair of the collision. Clearly such a colliding system carries a considerable orbital angular momentum, $L=bA\sqrt{s_{NN}}/2$. For example, a collision of AuAu at a center of mass energy of $\sqrt{s_{NN}}=200{\rm\,GeV}$ with $b=8$ fm has $L\simeq 1.6\times 10^{5}\hbar$. For comparison, the highest spin states of atomic nuclei (which are typically heavy ones with a spatial volume similar to the fireball in a heavy ion collision) have their total angular momentum values only at the order of $10^{1\sim 2}\hbar$ [171, 172]. The direction of this angular momentum is perpendicular to the reaction plane $\vec{L}\sim\hat{y}$. Given that each fast-moving nucleus also carries a large positive electric charge (e.g. $Z=79$ for Au), a collision configuration as illustrated in Fig. 18 (right) would also imply a strong magnetic field along the same direction as the angular momentum.

The majority of the initial angular momentum is carried away by the spectator nucleons, i.e. those outside the overlapping zone and not actively colliding, see Fig. 18 (left). Nevertheless, a fraction of the angular momentum, about $10\sim 20\%$ [173], remains to be carried by the hot matter created in the overlapping zone. Despite the subsequent complicated evolution of the hot fireball, angular momentum is conserved and therefore will be distributed across the system’s constituents, affecting both their orbital motion and their individual spins. As a consequence, as the fireball undergoes collective expansion in a fluid dynamical way, rotational orbital motion develops, which is characterized by nonzero fluid vorticity $\bm{\omega}$. Owing to the compressible nature of this fluid, it will develop a nontrivial distribution of local fluid vorticity which has a nonzero average $\langle\bm{\omega}\rangle$ over the whole fluid along the out-of-plane direction [173, 174, 175, 176, 177, 178]. In Fig. 19 (left) we show the results from AMPT model [173, 174] for the fluid velocity and vorticity distribution patterns, with the time evolution of averaged vorticity shown in Fig. 19 (right) for a wide range of collision energies. We note that the averaged vorticity, $\langle\bm{\omega}_{y}\rangle$, increases rapidly with decreasing collision energy, $\sqrt{s_{NN}}$. Thus the energy range explored by the RHIC beam energy scan (BES) appears to be much preferred for the search of any vorticity-driven effects. Furthermore the nonzero global angular momentum and average vorticity would also imply a preference of their spins to be aligned along the same direction due to the spin-orbital coupling [179, 180, 181, 182, 183, 178, 173, 174, 184]. A robust evidence of such global polarization effect has been reported recently by the STAR Collaboration, via sophisticated measurements of the hyperon spin orientation [185].

We next examine the strong magnetic fields from the two colliding ions which are similar to two oppositely-running strong electric currents, as depicted in Fig. 18 (right). An example of the spatial distribution of the magnetic field orientation and strength, calculated for Au-Au collisions at $b=8$ fm and $\sqrt{s_{NN}}=200$ GeV, is shown in Fig. 20 (left). In the overlapping zone the magnetic fields from each current add and become very strong, pointing approximately along the angular momentum (i.e. the out-of-plane) direction. A simple estimate suggests the strength of such field to be $(eB)\sim\frac{Z\alpha_{EM}\,\gamma}{R_{A}^{2}}$ where $\alpha_{EM}=1/137$, $R_{A}\simeq 1.1\,A^{1/3}$ fm is the nuclear radius and $\gamma=\sqrt{s_{NN}}/(2M_{N})$ is the Lorentz factor for the moving nucleons of mass $M_{N}\simeq 0.939$ GeV. For example a AuAu collision at $\sqrt{s_{NN}}=200$ GeV would have $(eB)\sim 3m_{\pi}^{2}$ or about $10^{17}$ Gauss, which is significantly stronger than any other known sources of magnetic fields. For comparison the magnetic field of a Magnetar, a neutron star which carries extremely strong magnetic field, is estimated to be on the order of $10^{14\sim 15}$ Guass [190]. Quantitative computations of such magnetic field strength, using event-by-event simulations, have been done at various beam energies and centralities as well as for different colliding systems [191, 192, 189, 193], with an example shown in Fig. 20 (middle) [189]. While many details of the magnetic fields are known at the beginning of a heavy ion collision, its subsequent time evolution is, however, less understood at the moment. The dominant contribution to the initial magnetic fields comes from the spectators which quickly fly away from the fireball. Thus, naively one would expect $eB$ in the collision zone to decay away rapidly, on a time scale as short as $\sim 2R_{A}/\gamma$. However, it has been argued that the partonic medium containing charged quarks/antiquarks could prevent the decrease of the magnetic field through the induction mechanism just as in the Lenz’s Law [194, 187, 188, 186]. Fig. 20 (right) shows the widespread model results and parameterizations of the magnetic field time dependence, illustrating the significant uncertainty associated with this quantity. To further constrain the time evolution of the magnetic field is the most important challenge for the quantitative modeling of magnetic field driven effects.

5.3 Experimental signals of anomalous chiral transport

Given that the conditions for anomalous chiral transport effects are in principle realized in the environment created in heavy ion collision, it is of great interest to examine their potential signals in order to verify these phenomena in experiment. Here, we will focus on the magnetic field driven effects, the CME and CMW, for which intensive experimental efforts have been carried out.

The CME leads to an electric current $\bf{J}_{e}$ along the magnetic field $\mathbf{B}$ direction which is approximately perpendicular to the reaction plane, as indicated by the grey plane in the left panel of Fig. 18 . This current implies transport of excessive positive charges toward one tip of the fireball (e.g. above the reaction plane) and transport of excessive negative charges toward the other tip of the fireball (e.g. below the reaction plane). As a result, a charge separation across the reaction plane emerges [124, 195], as illustrated in the left panel of Fig. 18. Such a spatial separation of positive/negative charges will convert into a charge dipole in the momentum space distribution of final state hadrons, and the strong outward or radial expansion of the fireball helps amplify this effect. Taking the configuration shown in Fig. 18 (left) as an example, the positive charges above the reaction plane will have their momenta preferably pointing upward while the negative point mostly downward. This results in a charge asymmetry in the positive/negative hadrons’ distributions in azimuthal angle $\phi$ with respect to the reaction plane (RP):

[TABLE]

There is however a subtlety: the CME current could be either parallel or antiparallel to the $\mathbf{B}$ field, depending upon the sign of chirality imbalance $\mu_{5}$ arising from event-wise fluctuations. Consequently, it is equally probable for the electric dipole in momentum space to be aligned or anti-aligned with the magnetic field or angular momentum. Therefore, the coefficient $a_{1}$ in Eq. (87) also fluctuates between positive/negative values from event to event and averages to zero, $\langle a_{1}\rangle\propto\langle\mu_{5}\rangle=0$. This simply reflects the fact that the parity symmetry is respected by QCD and is not violated in heavy ion collisions. What can be measured, however, is the variance of this charge separation coefficient, $\sqrt{\langle a_{1}^{2}\rangle}\propto\sqrt{\langle\mu_{5}^{2}\rangle}\neq 0$. To measure the variance of $a_{1}$, the following azimuthal correlator $\gamma$ for a pair of charged hadrons with either same-sign (SS) or opposite-sign (OS) was proposed [196]:

[TABLE]

where $\Psi_{2}$ is the 2nd harmonic event plane as an experimental proxy for the reaction plane. The correlation of pairs with same electric charge, $\gamma^{SS}$, has $\{\alpha\beta\}\to\{++\}$ or $\{--\}$ while that for opposite charged pairs, $\gamma^{OS}$, has $\{\alpha\beta\}\to\{+-\}$ or $\{-+\}$. A CME-induced charge dipole distribution in (87) will contribute to the above correlation: $\gamma^{SS}_{CME}\to-\langle a_{1}^{2}\rangle$ and $\gamma^{OS}_{CME}\to+\langle a_{1}^{2}\rangle$. To maximize the signal and reduce the backgrounds, one can further examine the difference between the correlation of same and opposite charged pairs, $\left(\gamma^{OS}-\gamma^{SS}\right)\sim 2\langle a_{1}^{2}\rangle$. Taking such difference would also make the charge-independent backgrounds cancel out. Extensive measurements of the $\gamma$ correlations were done for a variety of colliding systems (AuAu, CuCu, UU, PbPb, pPb, dAu) by STAR, ALICE and CMS over a wide span of beam energies [169, 197, 198, 199, 200].

The $\gamma$ correlator by itself, however, does not unambiguously determine if a pair is moving in parallel or anti-parallel to the reaction plane. For example: a pair of particles co-moving along the out-of-plane direction (which would be the case for CME-induced same-sign pairs) would have $\gamma<0$, but a pair of particles moving back-to-back along the in-plane direction would also have $\gamma<0$; likewise a pair of particles moving back-to-back along the out-of-plane direction (which would be the case for CME-induced opposite-sign pairs) and a pair of particles co-moving along the in-plane direction would both have $\gamma>0$. This ambiguity can be resolved by considering in addition the following correlator [201, 202]

[TABLE]

In this case, a pair of particles co-moving along the out-of-plane direction would have $\gamma<0$ and $\delta>0$, while a pair of particles moving back-to-back along the in-plane direction would also have $\gamma<0$ and $\delta<0$; likewise a pair of particles moving back-to-back along the out-of-plane direction would have $\gamma>0$ and $\delta<0$, while a pair of particles co-moving along the in-plane direction would have $\gamma>0$ and $\delta>0$.

A careful analysis of $\gamma$ and $\delta$ together indeed pointed to a dominance of background contributions [201, 202]. Sources of those non-CME backgrounds include e.g. transverse momentum conservation, local charge conservation, resonance decay and clusters, etc [201, 203, 204, 205, 206, 207, 208, 209, 210]. Note that the $\gamma$ correlator measures a difference between the in-plane and out-of-plane projected azimuthal correlations. As a result, the background contributions to $\gamma$ are controlled by the elliptic anisotropy coefficient $v_{2}$ that quantifies the difference between the in-plane versus out-of-plane collective expansion of the bulk matter. The substantial backgrounds pose a significant challenge for the experimental search of CME (see e.g. recent reviews [120, 202]). In order to remove or at least suppress backgrounds associated with the elliptic anisotropy, a two-component decomposition strategy was proposed [202]:

[TABLE]

where the $H$ represents the CME-like out-of-plane correlations while the $F$ represents the various background correlations whose contributions to $\gamma$ are modulated by the elliptic anisotropy $v_{2}$ as well as a kinematic factor $\kappa$. While one may expect $\kappa$ to be of the order unity, the precise value of $\kappa$ would sensitively depend on detector acceptance and analysis kinematic cut, as illustrated by the example of transverse momentum conservation in [204]. By measuring $\gamma$, $\delta$ and $v_{2}$, one may extract the $H$-correlation under certain plausible assumption of the $\kappa$ factor. Analysis based on this method was successfully carried out by STAR [211]. Clearly the $\gamma$ or $H$ correlations are not perfect due to strong background contamination, and a firm conclusion would require efforts to accurately calculate the signal, to get those backgrounds under full control, and to develop additional observables sensitive to the CME. We defer a detailed treatment of these issues to subsections 6.4 and 7.2 where the current status of measurements, their interpretations as well as pertinent background contributions will be discussed.

Another interesting experimental signal is related to the Chiral Magnetic Wave (CMW). As previously discussed, small fluctuations of vector and axial (or equivalently the RH and LH) charge densities propagate as CMW under the presence of a magnetic field, just like sound waves. In particular the RH wave propagates along the $\mathbf{B}$ field direction while the LH wave propagates in the opposite direction. If the fireball created in heavy ion collisions has a nonzero electric charge density (which is an equal mix of RH and LH) in its initial condition, then a simple calculation in terms of CMW solutions suggests that the RH wave moves toward one tip of the fireball along the out-of-plane direction while the LH wave moves toward the opposite tip. As a result, the nonzero initial (say, positive) electric charges will be transported toward the two tips along the azimuthal direction $(\phi-\Psi_{RP})=\pm\frac{\pi}{2}$. Due to charge conservation, this also leaves less positive charges around the “equator” region of the fireball along the azimuthal direction $(\phi-\Psi_{RP})=0$ or $\pi$. Such a spatial distribution pattern of the electric charge on the transverse plane has a nonzero quadruple moment with respect to the reaction plane, as first predicted and estimated in [139]. Again, the strong radial flow will “blow” all the charges outward and correlate their momenta with positions. As a consequence, the spatial quadruple pattern of the electric charge distribution will translate into a quadruple pattern in the final state charged hadrons’ azimuthal distributions with respect to the reaction plane. That is, there will be a little more in-plane charged hadrons of one sign while a little more out-of-plane charged hadrons of the opposite sign. Such a pattern implies a splitting in the elliptic flow coefficient $v_{2}$ for the positive/negative hadrons. The magnitude of this effect is proportional to the CMW-induced quadruple moment which in turn is determined by the amount of initial electric charge density. If the initial density is positive (negative), more positive (negative) charges will be transported by CMW to the out-of-plane tips of the fireball, leading to a smaller $v_{2}$ for the positive (negative) charged hadrons as compared with hadrons of the opposite sign. This initial density mainly comes from random “stopping” of the initial charges carried by the colliding nuclei. One could use the observed charge asymmetry $A_{ch}=\frac{N_{+}-N_{-}}{N_{+}+N_{-}}$ of an event (that is, the difference in the measured numbers, $N_{\pm}$ of positive/negative charged hadrons in this event) as a proxy which is correlated with the initial nonzero charge density of an event. It was predicted [139] that a specific splitting of the elliptic flow for positively and negatively charged pions should be observed:

[TABLE]

In other words, the CMW predicts a linear dependence of the splitting on the overall charge asymmetry $A_{ch}$, with the slope parameter $r$ computed from the CMW. The intercept $\Delta^{base}$ is independent of $A_{ch}$ and unrelated to CMW. This prediction was experimentally verified later by the STAR collaboration at RHIC [212]. Similar measurements were carried out by ALICE and CMS at the LHC as well [213, 214]. Certain potential background contributions were proposed and studied, such as the local charge conservation [215] and viscous transport [216]. At the moment the CMW appears to be the viable interpretation for the measured flow splitting signal at RHIC energies.

5.4 The beam energy dependence of anomalous chiral transport

In this last subsection, we discuss how the anomalous chiral transport may vary with the collision energy, which is an important question for the BES experimental program. The answer relies upon a number of key ingredients, to be discussed below.

First, the necessary condition for anomalous transport to occur is a sufficiently hot environment such that a substantial amount, in terms of both spatial volume and lifetime, of quark-gluon plasma with restored chiral symmetry is created 242424 As discussed in subsection 5.1.1., a vacuum condensate would spoil chirality of quarks and the restoration of chiral symmetry in high temperature phase is a prerequisite for CME. . One would then expect that with increasing collision energy this condition is satisfied better and better. Another relevant factor is the topological transition rate which controls the initial chirality imbalance for CME and which in general increases with the relevant medium scale, be it the pre-thermal saturation scale for the early stage or the temperature scale in the quark-gluon plasma. In either case, such a rate is expected to increase with the collision energy. 252525At present there is not yet a reliable way of precisely determining such a topological transition rate, which remains a major theoretical uncertainty. In fact, the main early motivation for studying CME was actually to help manifest and extract the topological transition rate. Thus it is natural to expect a “threshold energy” below which a QGP with sufficient topological transitions would ceases to exist in the fireball. A number of measurements from BES-I (such as the constituent quark scaling of elliptic flow, jet energy loss, directed flow, net proton fluctuations, light cluster production, etc) appear to hint at a qualitative change in the observed properties of the created bulk matter occurring around the beam energy range of $\sqrt{s_{NN}}\simeq(10\sim 20)$ GeV, as we shall discuss in Sec. 6 in more details. Therefore, one should expect any signal from anomalous chiral transport to decrease towards the low enough beam energy regime and to eventually turn off when the collision energy drops below the QGP production threshold.

The other necessary condition is the presence of the “driving forces”, i.e. the fluid vorticity $\bm{\omega}$ and the magnetic field $\mathbf{B}$. As already shown in Fig. 19, the vorticity $\bm{\omega}$ decreases with increasing beam energy and becomes negligible at the high energy end. The situation for the magnetic field is more complicated, due to an important difference between angular momentum and magnetic field: the former is conserved during the system evolution while the latter is strongly dependent on time. The peak strength of $\mathbf{B}$ scales with the nucleon Lorentz factor $\gamma\propto\sqrt{s_{NN}}$ and thus increases with increasing beam energy. On the other hand, without any medium feedback the time duration of this strong initial magnetic field scales inversely with $\gamma$ and thus decreases rapidly with increasing beam energy. This time scale can be estimated as $\tau_{B}\sim\frac{2R}{\gamma}\sim\frac{4RM_{N}}{\sqrt{s_{NN}}}$ with $R$ the nuclear radius (e.g. $R\sim 7$ fm for Au or Pb) and $M_{N}\simeq 0.939$ GeV the nucleon mass. Another important time scale is the formation time of the quarks after the initial collisions. Estimates of this timescale $\tau_{f}$ based on the glasma picture for the early stage in heavy ion collisions would suggest that $\tau_{f}\sim\frac{1}{Q_{s}}$, where $Q_{s}$ is the so-called saturation scale [60, 59]. $Q_{s}$ is expected to scale with beam energy in a specific way, ${Q_{s}}\sim{Q_{0}}\left(\frac{\sqrt{s_{NN}}}{E_{0}}\right)^{\frac{0.3}{2}}$ with $Q_{0}$ the saturation scale at a reference energy scale $E_{0}$. We may use the RHIC values $Q_{0}\simeq 1.5$ GeV and $E_{0}\simeq 200$ GeV for a quick order-of-magnitude estimate [60, 59]. Let us now compare the two time scales $\tau_{B}$ versus $\tau_{f}$, with the former decreasing rapidly with beam energy $\sqrt{s_{NN}}$ while the latter only decreasing mildly. Note that it is the quarks that are needed for both the anomalous chiral transport and for the medium induction mechanism which has the potential of prolonging lifetime of the magnetic field. One may therefore expect that at certain high enough beam energy where $\tau_{B}\ll\tau_{f}$, the initial magnetic field will exist only like an extremely short pulse before the formation of any quark or antiquark medium and would not cause any anomalous transport. That is, the signals of anomalous chiral transport effects may be expected to eventually disappear at the high beam energy end. Let us give a concrete estimate for beam energies relevant to RHIC and the LHC. At RHIC top energy $\sqrt{s_{NN}}=200$ GeV, one has $\tau_{B}\sim 0.13$ fm/c and $\tau_{f}\sim 0.13$ fm/c, with the two scales comparable $\tau_{B}\sim\tau_{f}$. At LHC energy e.g. $\sqrt{s_{NN}}=5020$ GeV, $\tau_{B}\sim 0.005$ fm/c significantly decreases from RHIC while the $\tau_{f}\sim 0.08$ fm/c only slightly decreases, resulting in a situation where $\tau_{B}\ll\tau_{f}$.

To conclude this brief discussion, it is quite plausible that potential signals from anomalous chiral transport effects would have nontrivial dependence on the collision beam energy, possibly disappearing in collisions both at very low energy end (e.g. below $\sim\hat{O}(10)\,\rm GeV$) and at very high energy end (e.g. beyond $\sim\hat{O}(1)\,\rm TeV$). From the estimate above, it appears quite likely that the optimal beam energy window may be in the range of $\hat{O}(10\sim 100)\,\rm GeV$. Due to this non-monotonic trend, the beam energy scan program at RHIC provides the unique opportunity to look for the signals from anomalous chiral transport. In fact, such a pattern of beam energy dependence could in itself be considered as a characteristic signature for the search of such effects.

6 The Beam Energy Scan Program at RHIC

6.1 Introduction

Soon after the observation that the matter created at the top RHIC energy collision behaves like a highly opaque, nearly perfect fluid [24, 25, 26, 27], subsequently dubbed as the strongly coupled Quark Gluon Plasma (sQGP), a natural question arose: how would such properties change as one increases or lowers the collision energy? One way to look at this question is from the perspective of QCD phase diagram as shown in Fig. 1. Generally speaking, an increase/decrease of the collision energy would create systems with higher/lower temperature and less/more net-baryon density. If nearly equilibrated QCD matter in the cross-over region at nearly vanishing baryon number density is being created at RHIC, the expectation would be that at even higher energies, available at the LHC, the created system would be in the same region and its basic thermodynamic properties should not change. This expectation has meanwhile been confirmed by measurements at the LHC [217, 218, 219, 220, 221]. Lowering the beam energy, and thus probing the baryon-rich region of the phase diagram, one would expect qualitative changes. Based on many models predictions, it is commonly expected that at sufficiently high baryon density the crossover turns into a first-order phase transition with a critical point marking the position in the phase diagram where the crossover ends and the first order transition begins. If the critical point and/or first order transition happen to be within reach of low energy collisions, one would expect certain significant changes in the properties of created matter. If, on the other hand, there is no true phase transition or the transition is at too high a value of the baryon density, one would at least expect the disappearance of the hot quark-gluon plasma at sufficiently low collision energy. In either case, it is plausible and tempting to expect qualitative changes in many observables as one lowers the collision energy. It were these considerations which gave rise to the RHIC beam energy scan (BES) program, which started its first phase in 2010.

The first phase of the beam energy scan program (BES-I) at RHIC has been completed with great success. The collision energies and event statistics for the BES-I as well as the planned second phase (BES-II) are listed in Table 1. The program covers the beam energies of $\sqrt{s_{NN}}=200,62.4,54.4,39,27,19.6,14.5,11.5,7.7\,{\rm\,GeV}$ 262626 By switching to fixed-target mode, the RHIC BES program has been extended further toward even lower collision energies: see caption of Table. 1. corresponding to a range of chemical potentials $25\lesssim\mu_{B}\lesssim 420{\rm\,MeV}$. 272727As we shall discuss below, at a given energy, especially at the lower beam energy, the value of the chemical potential varies with collision centrality. These values were extracted for the most central (0-5%) Au+Au collisions. Comprehensive experimental measurements of observables have been carried out by both the PHENIX and STAR collaborations, aiming to characterize the system’s bulk properties and its collective expansion as well as to search for evidences that are sensitive to a possible phase transition and/or critical point as well as the aforementioned anomalous chiral transport effects.

While in the following we concentrate on the results from the RHIC Beam Energy Scan, we should mention that there are other experiments addressing the same physics. There are on-ongoing efforts at SPS energies by the NA61/SHINE Collaboration and its predecessor NA49. Their main focus is the search for possible structures in the QCD phase diagram with various fluctuation and correlation observables for various energies and colliding systems [222]. In particular, NA49 and NA61/SHINE obtained results on the scaled variance [223, 224, 225, 226], strongly intensive quantities [227, 223, 228, 229, 230, 231, 232], intermittency [233, 234, 235], and many others [236, 237, 238, 239, 240]. Also, the HADES experiment at GSI Darmstadt, which measures at very low collision energies, has some capability to address fluctuation observables. In addition, there are several planned and approved experiments, namely: the Compressed Matter Experiment (CBM) at the future FAIR facility in Darmstadt, Germany, the Multi Purpose Detector (MPD) experiment at the NICA facility in Dubna, Russia, the CSR-External target Experiment (CEE) at HIAF in Hushou, China. Furthermore, there is the possibility to carry out an experiment at J-PARC in Tokai, Japan.

Let us now turn to reviewing the major experimental results from BES-I, with an emphasis on the beam energy dependence of the various observables.

6.2 Bulk properties and collectivity

In order to characterize the system created in a heavy ion reaction one typically measures the spectra for various hadrons from which one then obtains the hadron yields and extracts information about possible collective expansion, also referred to as flow. Further information about the size of the system and its lifetime is deduced from interferometry, or HBT measurements (see [241] for a review). In this subsection we briefly review the collision energy dependence of these measurements as obtained from the first phase of the RHIC beam energy scan (BES-I).

6.2.1 Particle spectra and ratios

As part of BES-I, both the STAR and PHENIX experiments have measured the collision energy as well as centrality dependence of the transverse momentum spectra [242, 243, 244, 245, 246] for the following hadron species: $\pi^{\pm}$, $K^{\pm}$, $p$, $\bar{p}$, $\phi$, $\Lambda$, $\bar{\Lambda}$, $\Xi^{\pm}$, $\Omega^{\pm}$ at the center of mass energies of 7.7, 11.5, 14.5, 19.6, 27, and 39 GeV. As an example, the results of the acceptance and efficiency corrected spectra for $\pi^{\pm}$, $K^{\pm}$, $p$, and $\bar{p}$ taken at $\sqrt{s_{NN}}=19.6{\rm\,GeV}$, are shown in Fig.21 [246]. These spectra are measured at mid-rapidity, $|y|\leq 0.1$, and over a transverse momentum range of $0.2\leq p_{T}\leq 2{\rm\,GeV}/\rm c$. The hadron yields at mid-rapidity are then obtained by extrapolating the measured $p_{T}$ spectra to $p_{T}=0$ and integrating them. In the case of the STAR measurements, this extrapolation gives rise to systematic uncertainties, which are of the order of 9-11%, 8-10% and 11-13% for pions, kaons and protons, respectively. Being a collider detector, for a given type of hadron, the acceptance of the STAR detector does not change as a function of collision energy. Therefore the systematic errors do not change much in these BES energies as compared with fixed target experiments. Given the hadron spectra and yields one can then proceed to extract some basic properties of the created system.

First, one observes that the shape of the spectra depends on the mass of the particles. This can be understood as the result of collective radial expansion of the system. As first discussed in [247], such a collective expansion leads to a larger apparent temperature similar to the well known blue shift for photons. In addition, for massive particles the spectrum flattens at low transverse momentum, and this flattening is more prominent for particles with higher mass. In order to test whether such a collective expansion is the plausible cause for the difference in the spectral shape of the various hadrons, as well as to extract the collective flow velocity and intrinsic temperature, one employs the so-called blast-wave model [248]. This model assumes a Hubble type radial expansion with a velocity profile, $\beta(r)=\langle\beta_{T}\rangle\left(\frac{r}{R}\right)^{\alpha}$, where $R$ is the radius of the source, $\langle\beta_{T}\rangle$ is the collective flow velocity at the surface, and $\alpha$ is some parameter which is of the order of unity. The blast-wave model works remarkably well in describing spectra of produced particles in heavy ion collisions at all energies. In Fig. 21, we show the fit for the most central collisions as dashed lines. Given the success of the blast-wave model, one can extract the flow velocity $\langle\beta_{T}\rangle$ and the temperature, $T_{fo}$, of the source at the moment of the so-called kinetic freeze-out when the particles cease to scatter with each other [249]. The kinetic freeze-out temperature and collective velocity parameters, $T_{fo}$ and $\langle\beta_{T}\rangle$, extracted from the blast-wave fits to the BES-I hadron spectra, are plotted in the right-panel of Fig. 22. For comparison, the results from 200 GeV Au+Au collisions at RHIC and 2.76 TeV Pb+Pb collisions at LHC are also shown [218]. The general trend seen in these data indicates that with decreasing beam energy the system freezes out at a higher temperature with smaller collective velocity. This may be understood as the consequence of a lower initial energy density and thus lower initial pressure of the system as well as a shorter fireball lifetime at lower collision energy. Similarly the smaller systems produced in the more peripheral collisions appear to have a higher freeze-out temperature and weaker collective expansion. The labels 1,2,3 in the plot indicate the most peripheral collisions of beam energy 7.7 GeV (1), 200GeV (2) and 2.76TeV (3) respectively, showing a stronger collectivity for the most peripheral collisions at higher collision energy.

The blast-wave analysis of the spectra provides information about the kinetic freeze out, i.e the moment the system seizes to interact entirely. The particle yields, on the other hand, allow us to deduce the conditions when the system freezes out chemically, i.e. when inelastic, or particle number changing processes, such as, e.g., $4\pi\leftrightarrow 2\pi$,282828We note that resonant (quasi) elastic processes, such as $n+\pi\rightarrow\Delta\rightarrow n+\pi$ do not change the number of final state nucleons or pions. are no longer effective. To extract the conditions, i.e. the temperature $T_{ch}$ and baryonic chemical potential $\mu_{B}$ at the chemical freeze-out, one compares the measured particle yield with the expectation of a hadron resonance gas (HRG) at certain given $T_{ch}$ and $\mu_{B}$. The HRG model is simply an ideal gas of all known hadronic resonances, with their abundances given by their thermal weights $N\sim\int d^{3}p\,n(p)$, where $n(p)$ is the appropriate thermal Fermi-Dirac or Bose-Einstein distribution functions for a given hadron resonance. The final state particles are then determined by taking into account all the decays of the unstable resonances (for details see e.g. [41, 255]). This approach, often referred as the thermal model for particle yields, also works remarkably well. In the left panel of Fig. 22 we show the extracted values for $T_{ch}$ and $\mu_{B}$ for the systems created in BES-I. The filled-circles and open-squares are from central (top 5%) and peripheral (60-80%) Au+Au collisions, respectively. The yellow-line represents the systematics of hadron resonance gas (HRG) fit of chemical freeze-out [253, 254] over a whole range of collision energies ranging from $2{\rm\,GeV}\lesssim\sqrt{s_{NN}}\lesssim 200{\rm\,GeV}$. The green-band represents the region of the cross-over transition as determined from lattice QCD [31, 252]. As one can see, the available collision energies of BES-I cover a wide range of baryonic chemical potential, $20{\rm\,MeV}\lesssim\mu_{B}\lesssim 420{\rm\,MeV}$, with the lowest energy corresponding to the highest value of $\mu_{B}$. While at a given collision energy the temperature $T_{ch}$ hardly changes with centrality, we observe a clear centrality dependence of the chemical potential, namely the more central collision tends to have larger chemical potential. The collision energy as well as the centrality dependence of the baryon number chemical potential can be easily understood. At lower energies and more central collisions, more of the incoming nucleons get stopped and thus increase the net-baryon number at mid-rapidity. We also observe, that at low baryon number chemical potential, $\mu_{B}\lesssim 300{\rm\,MeV}$, the chemical freeze temperature, $T_{ch}$, stays almost constant at a value close to the crossover transition temperature determined by lattice QCD. Above $\mu_{B}\gtrsim 400{\rm\,MeV}$, on the other hand, $T_{ch}$ starts to decrease rather rapidly according to the global fit (yellow line).

The observed higher baryon density at lower collision energy originates from the stronger baryon stopping. The effect of such baryon stopping is prominently seen in the energy dependence of the kaon to pion ratio, which is shown in Fig. 23. The BES-I results from the STAR experiment are plotted as circles in the figure. While the $K^{-}/\pi^{-}$ ratio (open symbols) shows a smooth energy dependence, the $K^{+}/\pi^{+}$ ratio (filled symbols) exhibits a maximum around $\sqrt{s_{NN}}\sim 8{\rm\,GeV}$. This behavior is well reproduced by the HRG-model [257] as seen by the yellow and green shaded areas. We note that the $\Lambda/\pi^{-}$ ratio (not shown here) also has a maximum at the same energy which is equally well reproduced by the HRG model [261]. The occurrence of this maximum can be understood by realizing that the threshold energy for associate production, $N+N\rightarrow N+\Lambda+K^{+}$ is $E^{associate}_{threshold}=1.58{\rm\,GeV}$ while that for kaon pair production is $E^{pair}_{threshold}=2.5{\rm\,GeV}$. Therefore at low energy, associate production dominates and gives rise to the steep increase of $K^{+}$ production compared to $K^{-}$ production, which are mostly produced via pair-production. With increasing collision energy very few baryons are stopped at mid-rapidity. Consequently, the contribution to the kaon yield from associate production drops so that it is dominated by kaon pair production. Therefore, the maximum of the $K^{+}/\pi^{+}$ ratio indicates the collision energy where the maximum baryon density at chemical freeze out is reached [262, 261]. For completeness we note that other model studies [263, 264, 265] have claimed that the peak might be an indication of the QCD phase transition. Indeed the existence of such a maximum was first proposed in Ref. [266] based on arguments concerning the onset of de-confinement. In Refs. [267, 268], the authors attempted to relate the peak in the ratio of $K^{+}/\pi^{+}$ to the restoration of chiral symmetry, which however relied on a rather complicated scheme in explaining why the chiral symmetry restoration does not affect $K^{-}/\pi^{-}$. At the moment, it appears that the simple HRG model provides a very successful and perhaps the most natural and “economic” explanation for the energy dependence of these ratios which serves as a useful indicator of the baryon-rich environment at low collision energy.

To summarize this subsubsection, the single particle observables such as spectra and hadron yields already characterize the bulk properties of the system very well. They provide information about the chemical and kinetic freeze-out conditions as well as the collective radial flow and baryon stopping. In order to gain further insight about bulk expansion dynamics, one needs to look at correlations, the most prominent of which are the azimuthal asymmetries often refereed to as directed, elliptic and higher order flow [22, 269, 270, 271, 272, 273, 274].

6.2.2 Collective flow and azimuthal asymmetries

The system created in off-center heavy ion collision initially will not be azimuthally symmetric in configuration space, but rather has roughly an elliptic shape with the long axis pointing perpendicular to the reaction plane. As a result, in the hydrodynamics framework, the pressure gradient and thus the expansion will be stronger in the in-plane direction. As a consequence, one expects an (elliptic) azimuthal asymmetry in the final state momentum distribution [275]. The measurement of this asymmetry, therefore, provides information about the pressure evolution of the system which in turn constrains the equation of state as well as the various transport coefficients such as the shear viscosity [276, 277]. Given the azimuthal direction of the reaction plane, $\Psi_{R}$, the momentum distribution in the transverse direction may be generally written as a Fourier series in the azimuthal angle $\phi$

[TABLE]

where the Fourier coefficients, $v_{n}(p_{T})$, characterize the azimuthal distribution for a given transverse momentum bin. Integrating the above expression over the transverse momentum and dividing by the number of particles, $N$, one arrives at the so-called integrated $v_{n}=(1/N)\int dp_{T}\frac{dN}{dp_{T}}v_{n}(p_{T})$. Due to the weighting with the transverse momentum distribution, the integrated $v_{n}$ are also sensitive to the aforementioned radial flow, or blue shift. Since for each event the orientation of the reaction plane $\Psi_{R}$ in the laboratory is different, the determination of the Fourier coefficients requires the measurement of at least a certain two-particle correlation. And indeed the $v_{n}$ can be extracted from two particle and higher order correlation functions, for example

[TABLE]

which also holds for each momentum bin separately. Here, the argument $\{2\}$ denotes that we use a two-particle correlation function. However, the above formula holds only under the assumption that no contributions from other intrinsic particle correlations are present. However, in reality hadronic decays give rise to azimuthal two particle correlations, so that a two-particle measurement, Eq. (93), will receive substantial corrections. In addition event-by-event fluctuations of the initial shape of the fireball introduce additional subtleties in determining the flow coefficients. Therefore, as discussed in detail in [278, 273], the extraction of the $v_{n}$ from higher order correlations, if possible, is preferred. However, for the discussion at hand, where we are more concerned with the energy dependence of the various $v_{n}$, the $v_{n}^{2}{\{2\}}$ are sufficient as the contamination for hadron decays are likely similar for all energies, and therefore, will not be responsible for a possible non-monotonic behavior. Let us now proceed to the energy dependence of the first two harmonics, $v_{1}$ and $v_{2}$.

Directed flow

The first harmonic, $v_{1}$, measures a dipole deformation of the azimuthal momentum distribution. Consider an off-central collision of a symmetric system such as Au+Au. If there is repulsion, forward going particle will be pushed to one side (along the impact parameter direction) and backward going particle will be pushed to the other side. Consequently $v_{1}$ will be similar in magnitude but of opposite signs for particles with positive and negative center of mass rapidities. Therefore, in order to measure the strength of this “bounce off” one considers the slope of $v_{1}$ at mid-rapidity, which is often referred to as directed flow. It has been proposed [268, 281] that the directed flow is sensitive to the equation of state and its softest point, i.e. the point where the speed of sound is minimum. It was predicted that the softest point, should reveal itself as a minimum in the excitation function of the net-proton directed flow, where its value is predicted to be negative. In Fig. 24 we show the result of a systematic study by the STAR collaboration on the energy dependence of the directed flow of net-protons as well as strange hadrons including Kaons and $\Lambda$ for Au+Au collisions at 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and 200 GeV at RHIC [279, 280]. There is a clear minimum in both the net-proton and net-lambda directed flow whereas the kaons show a monotonic decrease with collision energy. However, so far none of the available dynamical calculations involving a first order phase transition are able to even get a qualitative agreement with the STAR data [282, 283, 284]. Also, the slope of $v_{1}$ at mid-rapidity, i.e. the directed flow, for protons but also other baryons, is sensitive to the stopping of the incoming nuclei [285]. This may be the reason why the produced particles such as kaons do not show a minimum while the net protons do. 292929Indeed, when looking at protons and antiprotons separately, one finds that the energy dependence of the antiprotons is quite similar to all produced particles, whereas the protons exhibit a clear minimum at $\sqrt{s_{NN}}\simeq 15{\rm\,GeV}$. In addition, at lower energies the repulsive nuclear force among nucleons is expected to eventually result in a positive directed flow [286]. And indeed one way to obtain a directed flow consistent with the data is by including a repulsive interaction for the nucleons in the hadronic phase [284]. Whatever the detailed mechanism is, the different behavior for kaons and baryons clearly indicates that the bulk expansion dynamics in such low energy region is much more complex than that of a single fluid being pushed out as it appears to be the case in high energy collisions. The present status certainly calls for much better theoretical calculations, either hydro-based or hydro-transport-hybrid, for collisions in the BES energy range. Such calculations need to carefully model the initial conditions of hydrodynamics by fully accounting for the baryon stopping dynamics and for the pre-hydro baryon transport. Substantial progress has been made and developments continue addressing such needs [287, 288, 289, 290, 291, 292, 293, 294].

Elliptic and higher order flow

The second harmonic, $v_{2}$ corresponds to a quadrupole deformation of the azimuthal distribution and is expected to be symmetric around mid-rapidity (for symmetric collisions). As already mentioned, for off-central collisions the difference in the pressure gradients for in-plane and out-of-plane directions are expected to give rise to such a quadrupole deformation in momentum space. And indeed a finite and initially unexpected large $v_{2}$ has been observed. This, together with the observed strong jet-quenching [23, 295] have been the major steps towards the discovery of the so-called strongly coupled quark-gluon plasma. Furthermore it was found that the $v_{2}$ for baryons and mesons behaves such as if the baryons and mesons are produced via coalescence of constituents quarks, for example $1/2v^{\pi}_{2}(p_{T}/2)\simeq 1/3v_{2}^{proton}(p_{T}/3)$. This systematic pattern, often referred to as number of constituent quark (NCQ) scaling, works remarkable well for essentially all observed hadrons at top RHIC and LHC energies. With decreasing energy and increasing baryon number chemical potentials one observes an increasing difference in the elliptic flow of particles and anti-particles. However, the number of quark scaling still holds among particles and antiparticles separately [296, 297].

In Fig. 25 we present the second harmonic, $v_{2}(p_{T})$, for pions, $K^{0}_{S}$, protons, and deuterons for collision energies ranging from 7.7 GeV to 200 GeV at 0-80% centrality [298]. The dependence on particle mass of the observed $v_{2}$ at relatively small $p_{T}$ has been well reproduced by various hydrodynamical models as well as simple blast-wave model calculations, and can be understood as an effect of the previously discussed radial flow. We further observe that the magnitude of $v_{2}(p_{T})$ increases rather mildly with increasing collision energy and seems to level off for collision energies above $\sqrt{s_{NN}}\gtrsim 20{\rm\,GeV}$, contrary to the directed flow, $v_{1}$, for net-protons (see Fig. 24), the increase of $v_{2}(p_{T})$ is monotonic as a function of beam energy. This can also be seen in the beam energy dependence of the integrated $v_{2}$ for all charged particles, shown in Fig. 26.303030The change of sign of the integrated $v_{2}$ at $\sqrt{s_{NN}}\sim 4{\rm\,GeV}$ is understood to come from the shielding of the in-plane transverse expansion by the spectators. The increase of $v_{2}$ at even lower energies, on the other hand, reflects the energy dependence of the nuclear force which turns attractive at low energies. We note that the integrated $v_{2}$ receives also contributions from the radial flow. Therefore, we do not observe the leveling off in the integrated $v_{2}$ since the radial expansion velocity steadily increases with beam energy, see Fig. 22. We finally note, that higher order flow coefficients such as the triangular flow $v_{3}$ have also been measured over a wide span of beam energies [299], and the general trends of these higher order flow coefficients with changing beam energy resemble that of the elliptic flow.

Flow in small systems

Recently it has been observed that flow-like multi-particle correlations also exist in small systems, such as the ones created in high energy p+p, d+Au and p+Pb collisions. While it is still being debated if the origin of these correlations is hydrodynamic evolution as it is in nucleus-nucleus collisions [300, 301, 302], it is interesting to see if these flow-like correlations disappear at lower collision energies. Such a measurement has been recently carried out by the PHENIX collaboration [303] and in Fig. 27 we show their results for the beam energy dependence of $v_{2}(p_{T})$ for deuteron-gold reactions for collision energies ranging from $\sqrt{s_{NN}}=19.6{\rm\,GeV}$ to $\sqrt{s_{NN}}=200{\rm\,GeV}$. Similar to the $v_{2}$ in heavy-ion collisions shown in Fig. 25, the elliptic flow hardly changes over this range of collision energies.

To summarize this subsubsection, the measured collective flow coefficients demonstrate interesting patterns with changing collision energy. The elliptic and higher order flows appear to be mainly driven by anisotropic pressure gradients arising from the initial geometric shape which does not change much with beam energy and which occurs for both large and small colliding systems. The directed flow for baryons, on the other hand, shows highly nontrivial, non-monotonic collision energy dependence, with a minimum structure around $15{\rm\,GeV}$. The excitation function for directed flow of mesons, however, is monotonic. In addition a difference in the particle and antiparticle elliptic flow also becomes important in the low energy region. These observations suggest the presence of additional interaction dynamics in the bulk collectivity that pertains to the baryon-rich matter created at low collision energy and that is absent at high energy. Aside from an enhanced role of baryon stopping, this may provide a first hint that one probes the region in the phase diagram where the equation of state changes rapidly with the control parameters, as one would expect in the vicinity of a phase transition.

6.2.3 HBT and femtoscopic correlations

At very low relative momentum, correlations between a pair of identical particles are sensitive to the quantum interference resulting from the symmetry or anti-symmetry of the two-particle wave-function for bosons and fermions, respectively. Analogous to the well known Hanbury-Brown Twiss interferometry measurements of the sizes of stellar objects [304], one can utilize the quantum interference to determine the size of the systems produced in heavy ion collisions. These measurements are often referred to as HBT or femtoscopic correlation analysis, and are mostly done with pion, kaon and proton pairs. These measurements require appropriate corrections for the Coulomb interaction as well (in case of protons) correction for the s-wave scattering amplitude, see, e.g., the review [241] for details.

To extract the size of a colliding system, one examines the two-particle correlation function

[TABLE]

for a given total pair momentum $\vec{P}=(\vec{p_{1}}+\vec{p_{2}})$ as a function of the relative momentum $\vec{q}=(\vec{p_{1}}-\vec{p_{2}})$ and parameterizes it in the following form (in case of pions)

[TABLE]

Here, $N$ is an overall normalization, $\lambda$ represents the correlation strength and $F_{c}$ is the Coulomb correction factor [305]. In the function $G(\vec{q})$, we have decomposed the relative momentum $\vec{q}$ into three components: ${q}_{long}$ which points along the beam axis, and the transverse part $\vec{q}_{T}$ is further split into the component ${q}_{out}$ which points along the transverse pair-momentum, $\vec{P}_{T}$, and the component perpendicular to it, ${q}_{side}$. While the above formula is an approximation – corrections are discussed in detail, e.g., in Ref. [241] – it captures the relevant physics we wish to discuss here. In the above formula, one assumes that the source in configuration space is a three dimensional Gaussian. The parameters $R_{long}$, $R_{out}$, and $R_{side}$ then represent the spatial size along the beam axis, $R_{long}$, and the transverse sizes along the pair momentum, $R_{long}$, and perpendicular to it, $R_{side}$.

Suppose we have a system, where particles are slowly emitted from the surface in the transverse direction. In this case the apparent size in the “out” direction appears larger than that in the “side” direction, since the first emitted particle has already traveled a long distance before the last emitted particle emerges. In other words, the difference between $R_{out}$ and $R_{side}$ is sensitive to the emission duration, $\Delta\tau^{2}\sim(R_{out}^{2}-R_{side}^{2})$ [306, 307, 308]. Such a slow emission could, for example happen, if the speed of sound is very small due to a first order phase transition or a critical point [309, 310]. Therefore, a maximum in $R_{out}^{2}-R_{side}^{2}$ as a function of the collision energy could indicate the energy where the system gets closest to the “softest point” in the QCD equation of state. This is actually seen in the experimental data shown in the left panel of Fig. 28 [311]. In addition to the emission duration $\Delta\tau$ one can also attempt to estimate the expansion velocity from the HBT-correlation measurements. For a system which expands longitudinally in a nearly boost invariant fashion [312] the longitudinal size of a system $R_{long}$ provides a measure for the total lifetime of the system, $\tau_{life}\sim R_{long}$. The spatial expansion in the transverse direction, on the other hand, can be estimated by difference between the final side-ward radius as determined by $R_{side}$ and the initial size of the fireball $R_{init}=\sqrt{2}\bar{R}$, where $\bar{R}=\frac{1}{\sqrt{1/\sigma_{x}^{2}+1/\sigma_{y}^{2}}}$ with the Gaussian widths $\sigma_{x}$, $\sigma_{y}$ of the fireball density distribution determined by a Glauber model calculation. Therefore, the combination $(R_{side}-\sqrt{2}\bar{R})/R_{long}\sim\Delta R/\tau_{life}\simeq v_{expand}$ provides an estimator for the expansion velocity, $v_{expand}$. The beam energy dependence of this quantity is shown in the right panel of Fig. 28. We observe a minimum at about the same energy where the emission time exhibits a maximum. This combined behavior is rather intriguing and may be a hint about the softest point of the equation of state. However, a note of caution is in order here. First, the blast-wave analysis discussed above (see Fig. 22) shows no minimum of the expansion velocity. Second, both quantities shown in Fig. 28 change as much when going from the LHC energy (blue points) to the top RHIC energy, as they change from the top RHIC energy to the maximum/minimum at around $\sqrt{s_{NN}}\simeq 20{\rm\,GeV}$. However, for both the LHC and top RHIC energy collisions one probes the equation of state at approximately $\mu_{B}\simeq 0$. Therefore, in addition to a possible effect of the equation of state there is obviously some other mechanism which changes these two observables in between top RHIC energy and LHC energy. Unless this other mechanism has been identified and subtracted, it is too premature to conclude that the softest point of the equation of state has been found.

To summarize this subsection on bulk properties and collectivity, we find that the kinetic and chemical freeze-out temperatures as well as the radial and elliptic flow exhibit a rather smooth dependence on the collision energy. However, some observables such as the $K/\pi$ ratio, directed flow as well as the emission time extracted from HBT show a clear non-monotonic behavior, with maxima/minima occurring in a collision energy range of $10{\rm\,GeV}\lesssim\sqrt{s_{NN}}\lesssim 20{\rm\,GeV}$, corresponding to a baryon-rich region. To which extend these maxima/minima are related to the equation of state or are more a reflection of the change from baryon to meson dominance of the system cannot be easily concluded based on flow and global observables alone. To make further progress one needs to study observables that are more sensitive to a possible phase transition, such as fluctuations, which we will turn to next.

6.3 Criticality

6.3.1 Experimental data on proton, anti-proton, net-proton and net-charge high order cumulants

In the previous section, Sec. 4, we demonstrated that fluctuations, including fluctuations of conserved charges, are expected to be sensitive observables for the existence of a potential QCD critical point and/or first order phase transition. We further discussed that the correlation length, which would be infinite in a static infinitely large thermal system, will be limited to $\xi\simeq 2-3\rm\,fm$, because the systems created in heavy-ion collisions have a finite size and a finite lifetime. Therefore, it is advantageous to study cumulants and correlations of higher order, since they scale with the higher powers of the correlations length (see Sec. 4.5 and [85, 317, 318]). Meanwhile all collider heavy-ion experiments ALICE, PHENIX and STAR have carried out the measurement of cumulants of various orders and for several particle species [319, 320, 321, 316]. In the context of the RHIC BES-I, cumulants of net-proton, net-Kaon and net-charge multiplicity distributions have been measured as part of the search for the QCD critical region [320, 321, 316]. In the following we will discuss selected results from these measurements with an emphasis on net-protons, as they are believed to couple most strongly to the fluctuating $\sigma$ field near the critical point, as we discussed in Sec. 4.12.

Let us start with the net-electric charge distribution. In Figure 29 we show the beam energy dependence of the cumulant ratio $\kappa\sigma^{2}\equiv\kappa_{4}^{Q}/\kappa_{2}^{Q}$ for the net-charge distribution.313131 See the Appendix for discussions on quantities like kurtosis, skewness etc. and their relations to the cumulants. Panel (a) shows the results obtained by the STAR collaboration for two centralities (0-5%, filled-circles, and 70-80% open-squares) at collision energies $\sqrt{s_{NN}}=$ 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and 200 GeV [313]. The results from the PHENIX experiment, shown in panel (b), are for the top 5% central Au+Au collisions at energies $\sqrt{s_{NN}}=$ 7.7, 19.6, 27, 39, 62.4 and 200 GeV [316]. Within statistical and systematic uncertainties, both data sets show no significant features as a function of the collision energies over the range of $7.7{\rm\,GeV}\leq\sqrt{s_{NN}}\leq 200{\rm\,GeV}$, corresponding to a range of the baryon-chemical potential $420{\rm\,MeV}\gtrsim\mu_{B}\gtrsim 20{\rm\,MeV}$. We note that the statistical errors for the STAR measurement are significantly larger than those for the PHENIX data. The reason for this difference is essentially the larger acceptance of the STAR detector. It can be shown [322] that the statistical error for the measurement for the cumulant ratio $\kappa\sigma^{2}$ scales like323232While this scaling is true in general, there are distributions such as the bi-modal one discussed in Section 7.1.5 which are remarkably statistics friendly. $error(\kappa\sigma^{2})\propto\sigma^{2}/\left(\sqrt{N}\epsilon^{4}\right)$, where $N$ is the number of events, $\epsilon$ the detection efficiency, and $\sigma$ is the width of the observed net-charge distribution. Even though PHENIX recorded less events at $7.7{\rm\,GeV}$ and $19.6{\rm\,GeV}$ than STAR, see Table 1, the width of the *observed * net-charge distribution, which, to a good approximation is given by the total number of *observed * charged particles, $\sigma\sim N_{ch}$, is significantly larger for STAR. This is simply a result of the larger acceptance in the STAR detector. To be specific, the STAR events were collected with full azimuthal acceptance within $|\eta|<0.5$ and $0.2\leq p_{T}\leq 2.0{\rm\,GeV}/\rm c$ while the PHENIX events were taken with its two central arm spectrometers with $|\eta|<0.35$ and $0.3\leq p_{T}\leq 2.0{\rm\,GeV}/\rm c$, and each of the two arms cover $\pi/2$ radians in azimuth.

As mentioned in Section 4.12, the coupling between the fluctuating sigma field with baryons is expected to be stronger than that for pions, which are the most abundant charged particles. Therefore, the rather featureless beam energy dependence of the cumulant ratio $\kappa\sigma^{2}=\kappa_{4}/\kappa_{2}$ 333333In the following we will use the notation for the cumulant ratios $\kappa_{3}/\kappa_{2}=S\sigma$ and $\kappa_{4}/\kappa_{2}=\kappa\sigma^{2}$ interchangeably. for the net-charge distribution may very well reflect this weaker coupling. Since one expects a stronger signal in the baryon sector, for the rest of this Section we will focus our discussions on the energy dependence of the proton and anti-proton distributions.

In Fig. 30 we present the preliminary results from the STAR collaboration for the cumulant ratio $\kappa_{4}/\kappa_{2}=\kappa\sigma^{2}$ of net-protons (filled-circles), anti-protons (open-triangles) and protons (open-squares) for the 5% most central Au+Au collisions. We also show, as yellow band, the results obtained with the UrQMD event generator. We see that with decreasing collision energy, and thus increasing net-baryon density, the proton and net-proton $\kappa\sigma^{2}$ decreases until about $\sqrt{s_{NN}}\sim 20{\rm\,GeV}$, followed by a sharp increase above unity for collision energies below 10 GeV. As discussed in more detail in Section 4, such a non-monotonic behavior was predicted in [86] for a scenario where the freezeout line traverses the critical region, as shown in Figs. 13 and 14. The UrQMD prediction, on the other hand, exhibits a continuous monotonic decrease, which is mostly due to baryon-number conservation [112](see Section 7.1.2). As we shall discuss in Section 7.1.1, a more detailed analysis of all cumulants shows that the observed non-monotonic behavior could result from the interplay of negative two-particle correlations, which are likely the result of baryon number conservation, and rapidly increasing four-particle correlations, which then dominate the signal at the lowest measured energy of $\sqrt{s_{NN}}=7.7{\rm\,GeV}$.

The non-monotonic behavior of the cumulant ratio $\kappa_{4}/\kappa_{2}=\kappa\sigma^{2}$ is most pronounced for the most central collisions and gradually disappears as the collisions become less central. At the most peripheral centrality, 70-80%, $\kappa\sigma^{2}$ for net-protons is flat in the measured energy range $\sqrt{s_{NN}}=7.7-200{\rm\,GeV}$. This can be seen in Fig. 31 where in the left panel we show the STAR preliminary data [323] for $\kappa_{4}/\kappa_{2}=\kappa\sigma^{2}$ for several centralities. The right panel of Fig. 31 shows the third order cumulant ratio, $\kappa_{3}/\kappa_{2}\equiv S\sigma$, for the same centralities. Here the situation is different. While a non-monotonic behavior is seen in the most central collisions (filled-circles) all data from other centrality bins show a smooth increase as the energy decreases. In addition for the lowest two energies we see a reduction of $S\sigma$ with increasing centrality while for $\kappa\sigma^{2}$ we see an enhancement. As we shall discuss in Section 7.1.1, this decrease can be attributed to negative three particle correlations. At present the statistical uncertainties for $\kappa\sigma^{2}$ are rather sizable, especially for the collision energies below 20 GeV. This should be remedied in the planed second phase of the RHIC beam energy scan at RHIC, as we shall discuss in Section 8.4.

Note that when $\sqrt{s_{NN}}\lesssim 15{\rm\,GeV}$ all data points for $\kappa\sigma^{2}$ are above unity, for which a natural possible interpretation could be clustering of protons due to certain attractive interactions among them [325]. However, all known transport models fail to reproduce the observed enhancement of $\kappa_{4}/\kappa_{2}$ in the high baryon density region $\mu_{B}\gtrsim 300{\rm\,MeV}$ even when “attractive interaction” was turned on [326]. A recent work [325] suggested an enhanced attraction in the vicinity of the critical point. Using molecular dynamics the authors showed that this enhanced attraction lead to cluster formation, which in turn was able to reproduce the observed enhancement of $\kappa 4/\kappa_{2}$. However, this approach failed to reproduce the observed suppression of the third-order cumulant ratio, $\kappa_{3}/\kappa_{2}$, which is a general problem of cluster formation models, as we shall discuss in more detail in Sec. 7.1.1. As discussed earlier in this paper, Section 6.2.2, a similar sharp increase of net-proton directed-flow slope parameter ($dv_{1}/dy|_{y=0}$) has also been observed in the same energy region [279, 280], see Fig. 24. An increase in the directed flow, however, would rather suggest repulsive interactions. Also, as we shall discuss in Section 7.1.1, the third order factorial cumulant is negative, which would also naively suggest repulsion instead of attraction. Obviously, a simple picture of nucleon interactions does not suffice to explain the observed trends and more refined calculations are needed in order to understand the nature of the energy dependence of the observed correlation pattern.

Let us next turn to the dependence of the cumulant ratios $\kappa\sigma^{2}$ and $S\sigma$ on the size of the rapidity and transverse momentum space window. This will provide insight on the correlation range in momentum space. As discussed in Section 4.12 and in Ref. [110], one expects that the momentum correlations due to critical fluctuations should be similar to that of the thermal system, which translate into a rapidity correlation length of $\Delta y_{\rm corr}\simeq 1$. In Fig. 32 we present the preliminary STAR results for various transverse momentum windows (left-plots) and rapidity windows (right-plots) for net-proton $\kappa\sigma^{2}$ (top panels). In the bottom panels we show the same dependence for the third order cumulant ratio, $S\sigma$/Skellam, which is normalized by the expectation for an uncorrelated system which follows a Skellam distribution. As indicated in the Figure the data are for the top 5% Au+Au collisions. In case of the fourth order cumulant ratio $\kappa\sigma^{2}$, the maximum non-monotonicity is observed for the largest possible phase-space, i.e. $0.4\leq p_{T}\leq 2$ GeV/c and $|y|\leq 0.5$. When the phase space is reduced, either by reducing transverse momentum window or the rapidity coverage, the observed non-monoticity is reduced and the cumulant ratio $\kappa\sigma^{2}$ move closer to the expected value for an uncorrelated system, where $\kappa\sigma^{2}=1$. A similar trend is also seen in case of the normalized third order cumulant ratio, $S\sigma$/Skellam: the maximum deviation from unity is seen for the maximum phase-space available to the STAR experiment. Once the phase-space is reduced, the cumulant ratio again approaches the expectation for an uncorrelated system, $S\sigma/\rm Skellam=1$. As we shall discuss in more detail in Section 7.1.4, these data are consistent with long-range correlations in both transverse momentum and rapidity. However, we should keep in mind that the maximum rapidity window is rather small so that the present data are consistent with the expectation that the correlation length is $\Delta y_{\rm corr}\simeq 1$, the aforementioned expectation of a thermal rapidity correlation length. Clearly, the presently available acceptance of the STAR detector is too small to determine the momentum space correlation range, which would be a very important piece of information in general and in particular with regards to possible criticality. If the momentum correlation range is indeed thermal, as suggested in [110], then one should eventually see a leveling off of the cumulant ratio with increasing size of the rapidity window. In the future, as part of its upgrade program, the STAR experiment will extend the rapidity coverage for protons significantly. The upgrade extends the proton rapidity coverage from $|y|<0.5$ to $|y|<0.8$, as shown in Figure 33. Also shown are the expected statistical errors from BES-II for both 7.7 GeV and 19.6 GeV central collisions as green-bars, while the expected statistics from BES-II are listed in Table 1. Also shown in Fig. 33 as dashed lines are results obtained with the AMPT event generator. In this case the cumulant ratios decrease with increasing rapidity window as one would expect from baryon number conservation.

6.3.2 Efficiency corrections

Let us finish this subsection with some remarks concerning detection efficiencies and detector effects. Before comparing to theoretical expectations, all measured quantities will have to be corrected for experimental efficiencies including (i) detecting efficiencies, (ii) identification probabilities and (iii) and fiducial coverages. Although the last correcting factor involves a constant number for a given measurement, both (i) and (ii) involve the stochastic nature of the detector performance under specific circumstances. Furthermore, they are all dependent on the momentum of the particle under study. For the proton and anti-proton measurements discussed in this section, this implies that one must understand all of the aspects involving proton (and anti-proton) measurements at RHIC including effects of run-to-run fluctuating luminosity, pile-up, event-multiplicity, particle momenta and so on. For single spectra measurements, the most relevant are the efficiency corrections. In that case, the averaged values of the efficiency are often extracted, either as a function of either transverse momentum or rapidity or event-multiplicity or all of them, from the careful simulations of a given colliding system. However, for the multi-particle measurements, this procedure becomes much more involved as we shall briefly explain below. For more details see, e.g., Ref. [327].

In general, if there are $N$ produced particles, say protons, the detector will, due to finite detection efficiencies register $n$ particles, with $n\leq N$. However, the fraction of particles detected fluctuates from event to event. Therefore, we may detect the same number of protons for events with different number of produced protons, and, vice versa, we may detect different numbers of protons for events which have the same number of protons. To formalize this situation, let us denote by $P(N)$ the multiplicity distribution of the produced particles and by $p(n)$ the multiplicity distribution of the observed (detected) particles. The relation between the two is then given by

[TABLE]

where $B(n,N)$ is the probability to observe $n$ particles given $N$ produced particles. In general $B(n,N)$ may depend on phase space, particle multiplicity etc, which we ignore here for simplicity. The aforementioned (single particle) detection efficiency, $\epsilon$, is given by the ratio of the number of detected particles, $\left\langle n\right\rangle=\sum_{n}np(n)$, over the average number of produced particles, $\left\langle N\right\rangle=\sum_{N}NP(N)$,

[TABLE]

In practice, the detection probability, $\epsilon$, for single particles is, to first approximation independent from the presence of other particles. In this case the $B(n,N)$ is given by a binomial distribution with the Bernoulli probability equal to the detection efficiency $\epsilon$

[TABLE]

In this case it is easy to convince oneself that the single particle efficiency is indeed $\epsilon$. Also, in this case the efficiency to detect pairs

[TABLE]

is simply the product of two single particle efficiencies. And more generally one can show that the factorial moments for the produced protons, $F_{i}$, are given by

[TABLE]

where $\langle\frac{n!}{(n-i)!}\rangle$ is a factorial moment for the measured (detected) particles. Given the above relation between the factorial moments of produced and detected particles, Eq. (100), the efficiency corrections to the cumulants then simply involve expressing the cumulants in terms of factorial moments. We note that the preliminary results for the cumulant measurements shown in this section are all based on efficiency corrections based on the assumption that $B(n,N)$ follows a binomial distribution, Eq. (98) [320, 328].

As already mentioned, in reality, the situation is more complicated. For example, $\epsilon$ usually depends on the phase-space, $y,$ $p_{t}$ and $\varphi$. Also in case of net-protons one has to deal we two particle species, protons and anti-protons. In this case more refined methods are needed as discussed e.g. in Refs. [329, 330, 322, 331, 318, 332]. It is also possible that $B(n,N)$ in Eq. (96) is not given by a binomial distribution which seriously complicates the experimental analysis, see Ref. [327]. If $B(n,N)$ is binomial then we decide for each particle independently if it is detected or not. This procedure obviously does not introduce any (artificial) correlations. If $B(n,N)$ is non-binomial, on the other hand, the detection mechanism is not independent for each particle and artificial correlations are introduced into the system by the detector. In this case one needs to invert the matrix $B(n,N)$ so that one can extract the distribution of produced particles, $P(N)$ from that of the measured ones, $p(n)$,

[TABLE]

This procedure, which is commonly referred to as unfolding, is rather complicated [333, 334, 327]. A simple example is the situation where $\epsilon$ depends on the number of produced particles, $\epsilon=\epsilon(N)$. Usually the more produced particles the smaller efficiency and this is exactly the case for the STAR detector, see for example Fig. 1 in Ref. [323]. In Ref. [327] a simple model was considered where the efficiency linearly depends on the number of produced particles

[TABLE]

where $\epsilon_{0}=\sum_{N}P(N)\epsilon(N)$ is the average efficiency. It was found that even a very small value of $\epsilon^{\prime}\sim-0.0005$ can generate quite large effects. The problem of non-binomial efficiency and its effect on the measured cumulants are currently under study by the STAR Collaboration.

6.4 Chirality

In non-central heavy ion collisions, the positively charged and fast moving nuclei not only bring in tremendous amount of angular momenta but also produce a very strong magnetic field at mid-rapidity, as already discussed at length in Sec. 5.2. In a typical collision at RHIC, the angular momentum reached could be as large as $10^{4}\,\hbar$ [173] and the strength of the magnetic field is on the order of 1017 Gauss [124], both being aligned approximately along the direction perpendicular to the event plane. In analogy to the quantum hall effect in condensed matter physics where a strong magnetic field plays the crucial role, the angular momentum and magnetic field provide the necessary extreme conditions to induce novel quantum phenomena, such as the chiral magnetic and chiral vortical effect (see Sec. 5). These effect, if seen in experiment, in turn reveal nontrivial properties of the hot and dense matter created in heavy ion collisions such as chiral symmetry restoration.

The measurement of hyperon spin polarization with respect to the global event plane by the STAR collaboration during phase-I of BES has been the first experimental observation connected with the angular momentum of a high-energy nuclear collision. The results for Au+Au collisions at RHIC are shown in Fig. 34, where the observable $P$ quantifies the net polarization of hyperons, i.e. their level of preferential spin orientation along the direction of global angular momentum, which in turn is determined from the directed flow of forward- and backward-traveling fragments via the beam-beam counters [185]. Both $\Lambda$ and $\overline{\Lambda}$ hyperons show a nonzero global polarization at a level of a few percent. Model estimates indicate the vorticity of the system to be as high as $(9\pm 1)\times 10^{21}\,sec^{-1}$ [185, 183]. As an example, the result from the AMPT transport model is shown as yellow line in Fig. 34. The estimated magnitude of vorticity of approximately $10^{21}\,sec^{-1}$ is many orders of magnitude larger than the previously known maximum vorticity of about $10^{7}\,sec^{-1}$ in superfluids [185]. The hot QCD fluid created in non-central heavy ion collisions, therefore, represents the “most vortical fluid” so far produced in the laboratory. We further observe that the global polarization signal strongly increases with decreasing collision energy, in line with the arguments given in Sec. 5.2. Obviously, the observation of the global $\Lambda$ polarization provides strong evidence that the fireball created in these collisions carries angular momentum. Since the global polarization becomes stronger with decreasing collision energy, future experiments with the STAR fixed-target program, and at NICA, FAIR and CSR are well positioned to study this new phenomenon.

While the angular momentum with its associated vorticity is necessary for the spin polarization and possible existence of the Chiral Vortical Effect (CVE), the presence of a magnetic field is required for the Chiral Magnetic Effect (CME) and Chiral Magnetic Wave (CMW) in such collisions. As already discussed in Sec. 5.2 (see, e.g., Fig. 18), an observation of nonzero angular momentum and vorticity would be an indirect indication that the initial strong magnetic field may persist for a while in the system created as a consequence of the Faraday effect [194]. What is more, the systematic difference in the polarization between lambda (filled-circles) and anti-lambda (open-squares) in Fig. 34, with the latter having a larger polarization than the former, may have already provided a direct hint at the effect of the magnetic field. Clearly a magnetic field would affect the polarization of particles and anti-particles with opposite sign [183, 335, 336, 337]. However, as shown in [338], a vortical net-baryon current, may lead to a similar splitting. In addition, recent STAR measurements on the dielectron production at very low transverse momentum and low mass regime at RHIC [339], the difference in the directed flow $v_{1}$ for charmed-hadron $D^{0}$ and $\overline{D^{0}}$ [340] at RHIC as well as the large $v_{2}$ of $J/\Psi$ at high transverse momentum at LHC [341] might provide independent evidence for the strong magnetic fields that may persist long enough to affect observables in such collisions. We now turn to the measurements related to the anomalous chiral effects driven by the magnetic field.

As discussed in Sec. 5.3, one observable consequence of the CME is the electric charge separation along the direction of the external magnetic field. Such a charge separation can be measured through the charge-dependent azimuthal correlator $\gamma^{OS}$ and $\gamma^{SS}$, which are defined in Eq. (88). These observables have been measured at both RHIC and the LHC over a wide span of collision energies [169, 197, 198, 199, 200] and are shown in Fig. 35. For collision energies of $\sqrt{s_{NN}}\gtrsim 19.6{\rm\,GeV}$ we observe a significant charge asymmetry in the correlators, that is, a difference between $\gamma^{OS}$ and $\gamma^{SS}$. This difference disappears for the lowest two energies, $7.7{\rm\,GeV}$ and $11.5{\rm\,GeV}$. Also, as discussed in Sec. 5.3, one expects for the CME that $\gamma^{OS}>\gamma^{SS}$ and that this difference increases toward more peripheral collisions. This trend can be seen in the data.

However, as discussed in detail in Sec. 5.3, at this point we should remind the reader that the interpretation of the $\gamma$-correlator is complicated by strong background contributions. These backgrounds arise from certain genuine two-particle correlations coupled with the anisotropic collective expansion, and their contributions to $\gamma$ are thus (roughly) linearly dependent on the elliptic coefficient $v_{2}$, see Eq. (90). One possible way to study and potentially remove this flow induced background is by means of the so-called event-shape engineering analysis [342, 343]. Within a given centrality (where the magnetic field is expected to approximately stay the same), one could further divide the events into subgroups according to their measured $v_{2}$ and then examine how the $\gamma$-correlator varies with respect to $v_{2}$. An example of such an analysis for Pb+Pb collisions at $2.76\rm\,TeV$ by ALICE is shown in Fig. 36 [343]. There, in a given centrality class we see a more or less linear dependence of $\gamma$ on $v_{2}$. An extrapolation of this trend toward $v_{2}=0$ indicates a possibly nonzero intercept that could be the CME signal. Combing results from the $10-50\%$ centrality range, ALICE concludes that the fraction of CME signal could be about $8\sim 10\%$ of the measured $\gamma$ correlators but this is still subject to significant systematic and statistical uncertainties.

Another approach [202], discussed in detail in Sec. 5.3, is to decompose the $\gamma$-correlator into one part which scales with $v_{2}$, and one which does not (see Eq. (90)). An analysis along this line was completed with BES-I data by STAR [211], with suitable assumptions for the kinematic $\kappa$ factor introduced in Eq. (90) and discussed therein. Fig. 37 shows the STAR results [211] for the collision energy dependence of the component $H$, which is independent of $v_{2}$ and thus may be interpreted as the signal for charge separation. These results are from the mid-central ($10-30$%) Au+Au collisions from $7.7{\rm\,GeV}$ to $200{\rm\,GeV}$. For comparison, the ALICE result [198] is also shown on the high beam energy end in the figure. One uncertainty involved in such an analysis, as discussed already in Sec. 5.3, is the choice of the $\kappa$ factor for background correlations. As shown in, e.g., Ref. [204, 202], the value of this factor is sensitive to the momentum spectra of charged particles, the differential flow coefficient $v_{2}(p_{T})$, the detector acceptance and kinematic cuts, etc. The factor $\kappa$ may be estimated by performing simulations of the bulk evolution (without the presence of CME) with e.g. event generators such AMPT or UrQMD or using hydrodynamics. Such an exercise, taking into account the STAR acceptance as well as kinematic cuts, leads to values of $\kappa$ in the range of $1\lesssim\kappa\lesssim 2$ [344]. For the analysis shown in Fig. 37, three lines are presented corresponding to the choice of $\kappa=1,1.5,2$ respectively to reflect the influence of this uncertainty. As one can see from the Figure, the so-obtained signal shows a nontrivial dependence on the collision energy, reaching a maximum around $39$ GeV while approaching zero when the energy is lower than 11 GeV. This specific trend appears in qualitative agreement with expectations from CME, as discuss in Sec. 5.4. The apparent decrease of the H-correlator toward the very high energy end deserves a special note. The ALICE collaboration [343] combined a similar two-component decomposition strategy with the event-shape engineering analysis to constrain the level of a possible CME signal at $2.76\rm\,TeV$ in the measured $\gamma$-correlator and the results are consistent with the above H-correlator results. Furthermore, recently the CMS collaboration reported its results on charge dependent azimuthal particle correlations with respect to the second-order event plane in both p+Pb and Pb+Pb collisions at 5.02 TeV [199]. At the same multiplicities, a similar behavior in p+Pb and Pb+Pb collisions is observed. A further analysis by CMS carefully examined such correlations in Pb+Pb collisions with respect to both the second and third harmonic flow event planes and attempted to extract the background contributions [200]. These analyses put a very stringent upper limit, no more than about $4\%$ for Pb+Pb according to CMS, on any potential CME fraction in the $\gamma$-correlator at 5.02 TeV collisions [200]. The ALICE and CMS measurements at the LHC, therefore, consistently suggest a gradual decrease of the potential CME signal with increasing beam energy and its possible disappearance at the very high energy end (e.g. $\sim 5$ TeV and beyond).

As discussed in Secs. 5.1 and 5.3, in addition to the CME, the magnetic field can also induce the Chiral Magnetic Wave (CMW). The CMW predicts a specific splitting between the elliptic flow, $v_{2}$, of $\pi^{-}$ and $\pi^{+}$, as given in Eq. (91) of Sec. 5.3, with a linear dependence on the charge asymmetry $A_{ch}=\frac{N_{+}-N_{-}}{N_{+}+N_{-}}$. Such a feature is indeed confirmed by experimental data, as shown in Fig. 38 (left). In the CMW scenario, the slope parameter $r$ extracted from this linear dependence is directly related to the CMW-induced quadruple moment of the charge distribution in the QGP. The estimated values for $r$ from CMW calculations are in good agreement with STAR data [139, 212].

The STAR measurements for the collision energy dependence of the parameter $r$ are shown in Fig. 38 (right), for Au+Au collisions in two centrality bins, 40-50% and 50-60% at RHIC BES energies. As one can see in this figure, at high energies, $\sqrt{s_{NN}}>20$ GeV, the slope parameter is clearly finite, $r>0$, consistent with the CMW expectations. At energies below 15 GeV, on the other hand, $r$ is consistent with zero within the albeit large error bars. As discussed earlier in Sec. 5.3, the effect of the CMW requires the existence of QGP with chiral restoration. Therefore, the vanishing trend of the $r$ parameter at low energies could be an indication of broken chiral symmetry in the hadronic phase dominating the fireball at such energies [212]. A similar nonzero slope parameter was also measured and reported by ALICE for Pb+Pb collisions at $2.76$ TeV [213]. At the even higher energy of $5.02\rm\,TeV$, the CMS measurements [214] appear to suggest a negligible CMW signal at such energy, which is similar to and consistent with the CME measurements by CMS at the same energy.

Let us briefly summarize this subsection, which focused on measurements pertinent to the chirality aspect in high energy nuclear collisions. There are strong experimental evidences for a nonzero angular momentum and vorticity structure in the fireball created by such collisions. There also appear compelling indications for the existence of very strong magnetic fields that may last for sufficient amount of time during the fireball evolution to induce observable effects. Given such external conditions, one expects the occurrence of anomalous transport like the CME and CMW. An unambiguous observation of these effects would provide, for the first time, the long sought-after experimental evidence for a QGP with restored chiral symmetry as well as the presence of gluonic topological fluctuations. The CME is predicted to induce a charge separation that can be measured through charge dependent azimuthal correlations, while the CMW is predicted to induce a charge quadruple in the fireball that can be measured through charge dependent elliptic flow. These measurements have been done over a wide span of collision energies from RHIC to LHC, with interesting hints of CME and CMW signals that qualitatively agree with theoretical expectations. A conclusion however can not be drawn yet, due to strong background contamination that is hard to separate unambiguously or evaluate precisely. The current understanding of these measurements and relevant issues will be further discussed in the next section. However, it is important to note here that the potential signals of both CME and CMW appear to have a non-monotonic dependence on collision energy, vanishing below a certain threshold energy while also disappearing at very high energies. In light of the discussions on the beam energy dependence of anomalous chiral transport in Sec. 5.4, this nontrivial trend may actually be considered as a characteristic to be expected for CME and CMW signals. It also appears that the optimal window for the precision search of anomalous chiral transport effects would be in the RHIC BES energy region.

We end this subsection by mentioning the exciting prospect, brought by the collisions of isobars [345], for disentangling the background contamination and unambiguously identifying the CME signal. These experiments were completed in 2018 and substantial amount of data has been taken, which are being processed and analyzed at the time of this writing. A more detailed discussion of the isobaric collisions will be presented in Sections 7.2 and 8.

7 Discussions on BES-I Results

In the following we will be discussing the experimental results obtained during the first phase of the RHIC beam energy scan, which we have presented in the previous section. We will concentrate on the fluctuation measurements, which are relevant for the critical point and first-order transition search, and the results for the correlation functions which are sensitive to anomalous transport.

7.1 Discussion of the STAR data on proton number cumulants

In Section 6.3.1 we have presented the current status with regards to the measurement of particle number cumulants. Especially the beam energy dependence of the net-proton cumulants show some rather interesting features for energies below $\sqrt{s_{NN}}\lesssim 19{\rm\,GeV}$. In this Subsection we present several observations regarding the preliminary STAR data on the (net)-proton cumulants. We will focus mostly on the lowest energies, $\sqrt{s_{NN}}\leq 19$ GeV, where several experimental observables are characterized by non-monotonic behaviour, including the cumulant ratios of the proton distribution. We will further discuss effects on the cumulants due to baryon-number conservation and participant (or volume) fluctuations, which need to be understood in order to relate the measurements with most model calculations. Finally, since the number of anti-protons is negligible for energies below $\sqrt{s_{NN}}\lesssim 19{\rm\,GeV}$ in the following we will discuss mostly cumulants of the protons, which simplifies some of the developments.

7.1.1 Proton factorial cumulants

As already eluded to in Section 4 factorial cumulants provide an alternative measure of fluctuations. As discussed in A, factorial cumulants represent the integrated irreducible (genuine) correlation functions. Therefore, they are more sensitive to the underlying interactions and have no contribution from simple finite number statistics, i.e. they vanish in absence of any correlations. The factorial cumulants of the proton distribution, $\hat{C}_{n}$, are given in terms of the regular cumulants, $\kappa_{n}$, by (see Eq. (143))

[TABLE]

Based on the above relations and using the preliminary STAR data [323] measured in the acceptance region $|y|<0.5$ and $0.4<p_{t}<2.0$ GeV, we obtain the results for the factorial cumulants $\hat{C}_{2}$, $\hat{C}_{3}$ and $\hat{C}_{4}$ presented in Fig. 39. To best illustrate the contributions of the various factorial cumulants to the fourth-order cumulant $\kappa_{4}$ we plot them with the appropriate weights: $7\hat{C}_{2}$, $6\hat{C}_{3}$ and $\hat{C}_{4}$, see Eq. (143). For the most central collision points, which shall be the focus of the subsequent discussion, the average number of protons is approximately $\left\langle N\right\rangle\approx 39,$ $31$ and $25$ for energies $7.7$, $11.5$ and $19.6$ GeV, respectively [323, 318].

As can be seen in Fig. 39, at $7.7$ GeV the fourth order cumulant in central collisions is clearly dominated by the fourth-order factorial cumulant and thus there are substantial four-proton correlations present in the system. We obtain

[TABLE]

As we increase the energy the situation is clearly changing. At $11.5$ GeV, the $\hat{C}_{4}$ and thus four particle correlations are still sizable. However, since the contribution of $\hat{C}_{4}$ to the fourth order cumulant, $\kappa_{4}$, is approximately balanced by these from $\hat{C}_{2}$ and $\hat{C}_{3}$, the presence of the four-particle correlation is less apparent in the forth order cumulant. This is a nice example of the usefulness of factorial cumulants, which exhibit the correlation structure more directly than the regular cumulants, $\kappa_{n}$. At $19.6$ GeV, for the most central collisions, $\kappa_{4}$ is dominated by the (negative) second order factorial cumulant, $\hat{C}_{2}$. The results shown in Fig. 39 can be summarized as follows: As we decrease the energy from $19.6$ GeV to $7.7$ GeV, the relative importance of the four-order factorial cumulant, and therefore four-proton correlations, is increasing. A more detailed discussion, including the centrality dependence, can be found in Ref. [346].

To put the above numbers, and in particular the value of $\hat{C}_{4}$ at $7.7$ GeV into perspective, let us consider a simple model, already described in the A, see also [347]. Suppose that we have clusters of protons, distributed according to the Poisson distribution, which always decay into $m$ protons. According to Eq. (155), we obtain $\hat{C}_{k}=\left\langle N_{c}\right\rangle\frac{m!}{(m-k)!}$. Having four-proton clusters ($m=4$) we obtain $\hat{C}_{4}=24\left\langle N_{c}\right\rangle$ and in order to get $\hat{C}_{4}=170$, see Eq. (104), we need to assume $\left\langle N_{c}\right\rangle\sim 7$. Consequently, on average roughly $28$ of the observed 40 protons should originate from such clusters. For five-proton clusters ($m=5$) we have $\hat{C}_{4}=120\left\langle N_{c}\right\rangle$ thus $\left\langle N_{c}\right\rangle\sim 1.5$. Given these results, it is clear that $\hat{C}_{4}\sim 170$ is a rather large number and a strong source of four-proton correlations is required to explain the magnitude of the preliminary STAR signal. However, the simple cluster model also predicts all $\hat{C}_{k}$ to be positive in obvious disagreement with the preliminary data, where both $\hat{C}_{2}$ and $\hat{C}_{3}$ are negative. Therefore, it is not straightforward to explain the observed large correlations and their signs with baryonic clusters or any other more conventional means. Could this be the first hint for the critical point and/or first order transition? This question will be addressed to some extent in Section 7.1.5. Before we turn to this and further discussion of the data, however, let us next discuss corrections due to baryon-number conservation and participant (volume) fluctuations.

7.1.2 Baryon number conservation

Since baryon number is conserved in a heavy ion reaction, the fluctuation of the protons is different that that obtained, e.g., in a grand canonical ensemble, where most theoretical calculation are carried out. Therefore, it is essential to determine the corrections due to baryon number conservation on the various fluctuation measures such as cumulants. At low energy, $\sqrt{s_{NN}}\lesssim 20{\rm\,GeV}$, where the cumulants show the most interesting features almost all measured baryons originate from the incoming nuclei. The baryon number is simply decelerated to the mid-rapidity region. The simplest model, which does not introduce any new correlations, is to assume that baryons stop independently in a given rapidity interval. In this case the distribution of protons is given by the binomial distribution

[TABLE]

where $p$ is the probability that the initial baryon will end up in our rapidity bin as a proton. $B$ is the initial number of baryons that can potentially be observed as protons. In the following we will assume $B\approx 400$, being approximately the total baryon number in both colliding nuclei. Here we assume that all initial baryons are subjected to a maximum isospin randomization. That is, the initial proton can be stopped either as a proton or a neutron with equal probability, see Refs. [348, 349].

The generating function reads343434For one baryon, the probability to be observed as a proton in a given interval is $p$ and the probability not to be observed is $1-p$. Thus, the generating function for a single baryon equals $(1-p)z^{0}+pz^{1}$. Having $B$ independent sources, the final generating function is given by a product of $B$ one-baryon generating functions, which gives Eq. (106).

[TABLE]

and $\hat{C}_{1}=\left\langle N\right\rangle=pB$, and

[TABLE]

where $\left\langle N\right\rangle$ is the average number of measured protons, which varies with the size of the chosen rapidity interval. Taking $B=400$ and for $7.7$ GeV $\left\langle N\right\rangle=40$ we obtain: $\hat{C}_{2}=-4$, $\hat{C}_{3}=0.8$ and $\hat{C}_{4}=-0.24$. We observe that $7\hat{C}_{2}=-28$ is of the same order of magnitude and roughly a factor of $2$ larger in magnitude than the STAR data. $6\hat{C}_{3}=4.8$ is significantly smaller than in Eq. (104) and the baryon conservation driven $\hat{C}_{4}=-0.24$ is completely negligible, when compared to STAR $\hat{C}_{4}\approx 170$. Finally, we note that the above results for the factorial cumulants due baryon number conservation together with Eq. (143) give a cumulant ratio of $\kappa_{4}/\kappa_{2}=0.46$ for $7.7{\rm\,GeV}$, consistent with the UrQMD result shown in Fig. 30.

We conclude that at lower energies, the effect of baryon conservation plays some role in $\hat{C}_{2}$, is an order of magnitude too small for $\hat{C}_{3}$, and is almost three orders of magnitude too small for $\hat{C}_{4}$. This shows a clear advantage of the factorial cumulants when compared to the cumulants, which are always ”contaminated” by $\hat{C}_{2}$. A detailed discussion of the baryon-number conservation effects, including the presence of anti-baryons, can be found in Refs. [113, 350].

7.1.3 Volume fluctuation

Another source of fluctuations in heavy ion collisions arises from the fact that even with the tightest centrality cuts the impact parameter and thus the size of the created system fluctuates. These volume or participant fluctuations [351] give rise to additional non-dynamical fluctuations which need to be understood and, if large, subtracted from the experimentally determined fluctuation measures. This aspect has been addressed in, e.g., Ref. [347], where a minimal model, incorporating the volume fluctuation and baryon conservation at $\sqrt{s_{NN}}=7.7$ GeV was discussed. In general, for a given centrality class, the fluctuations of the volume of the system is due to the fluctuations of the number of participating or wounded nucleons, $N_{\text{part}}$. Following Ref.  [347] the distribution of number of participating nucleons $P(N_{\text{part}})$ for a given centrality selection can be determined from a standard Glauber model. In addition, at low energies, where the production of baryon–anti-baryon pairs is negligible, the number of wounded nucleons in a given event is, to a good approximation, identical to the total number of baryons. In absence of any correlations, each wounded nucleon will then end up with a probability $p$ as a proton a given rapidity (or acceptance) bin. Consequently, the distribution of protons at a given $N_{\text{part}}$ is described by the binomial distribution. Note, in this approach we assume that each wounded nucleon is decelerated independently. Any deviation from this assumption would introduce new correlations and we are interested in the minimal model, where only $N_{\text{part}}$ fluctuations (and baryon conservation) are present. Therefore, in this model the distribution of protons is given by

[TABLE]

with the generating function

[TABLE]

This situation is very similar to the one discussed above in the context of baryon number conservation, see Eq. (105), except that now in each collision the maximum number of baryons that can end up in our rapidity bin is not the total baryon-number, $B$ but the number of wounded nucleons, $N_{\text{part}}$, which fluctuates from event to event. A straightforward calculation then gives $\hat{C}_{1}=\left\langle N\right\rangle=p\left\langle N_{\text{part}}\right\rangle,$ with $\left\langle N_{\text{part}}\right\rangle$ being the average number of wounded nucleons in a given centrality class. $\langle N\rangle$ is the average number of observed protons. The higher order factorial cumulants read

[TABLE]

and $\hat{C}_{4}$ ($\hat{C}_{4}/\langle N\rangle^{4}$ to be more precise) can be found in Ref. [347].353535In Ref. [347] $\hat{C}_{k}/\langle N\rangle^{k}$ is denoted by $c_{k}$. As seen from Eq. (110) the fluctuation of $N_{\text{part}}$ contributes to the factorial cumulants.

Applying the centrality cuts similar to those used in the STAR measurement [323] (for details see Ref. [347]) one arrives at results presented in Fig. 40, where we plot the factorial cumulants as a function of the average number of wounded nucleons, $\langle N_{\text{part}}\rangle|_{N_{\text{ch}}}$, and the average is calculated at a given number of produced pions. The solid lines represent the full calculations and the dashed lines represent the results without $N_{\text{part}}$ fluctuations, that is, only the leading terms in Eq. (110) are included, $\hat{C}_{k}\sim\langle N\rangle^{k}/\langle N_{\text{part}}\rangle^{k-1}$. The open symbols represent the results after averaging $\hat{C}_{k}$ over bins in centrality $0-5\%$, $5-10\%$, etc. (the first five points are shown). We see that, in general volume fluctuation introduces some nontrivial features, however the signal is negligible for $\hat{C}_{4}$ when compared to the preliminary STAR data, see Eq. (104). $\hat{C}_{3}$ is visibly modified by $N_{\text{part}}$ fluctuations but again the signal is way to small to explain the STAR data. Interestingly, $\hat{C}_{2}$ is strongly modified by volume fluctuation and the signal is modified (when compared to the case without $N_{\text{part}}$ fluctuations) by roughly a factor of two.

From this we conclude that at $7.7$ GeV both $\hat{C}_{3}$ and $\hat{C}_{4}$ in central collisions cannot be explained by a simple volume fluctuation or baryon conservation.

7.1.4 Rapidity dependence

Next, we will comment on the acceptance dependence of the cumulants which has been discussed already in Sections 4.12 and 6.4. STAR has measured cumulants up to the maximum rapidity interval given by $|y|<0.5$. Therefore, the maximum rapidity distance between two measured protons is one unit of rapidity. We will demonstrate that the preliminary STAR data are consistent with a constant long-range rapidity correlation within this interval. Here, for simplicity we assume that the single-particle rapidity distribution, $dN/dy\equiv\rho(y)$, can be approximated by a constant in $|y|<0.5$, $\rho(y)=\rho_{0}=const$ which is not far from reality [352].

Long-range rapidity correlations mean that the normalized multi-particle correlation function, defined as363636In Refs. [346, 353] this object is denoted by $c_{k}(y_{1},...y_{k})$.

[TABLE]

does not depend on the rapidities of the particles, $y_{i}$,

[TABLE]

Therefore, the factorial cumulants

[TABLE]

scale with $\left\langle N\right\rangle^{k}$, where in the above equations the integrations are performed over a rapidity interval $-\Delta y/2<y_{i}<$ $\Delta y/2$. Since the single particle rapidity density is, to a good approximation, constant, $\rho(y)=\rho_{0}$, the average number of particles scales with the size of the rapidity interval, $\left\langle N\right\rangle=\rho_{0}\Delta y\sim\Delta y$, and we obtain373737We should also point out that short-range correlations lead to $\hat{C}_{k}\sim\Delta y$ and the cumulant ratios do not depend on $\Delta y$, see, e.g., Ref. [346]. For example assuming $R_{2}(y_{1},y_{2})\sim\delta(y_{1}-y_{2})$ we have $\hat{C}_{2}\sim\int dy_{1}dy_{2}\rho(y_{1})\rho(y_{2})\delta(y_{1}-y_{2})$ which gives $\hat{C}_{2}\sim$ $\int dy\rho(y)^{2}\sim\Delta y$.

[TABLE]

Following Eq. (143) we see that the cumulants $\kappa_{i}$ scale with $\left\langle N\right\rangle$ and $\Delta y$ in a rather non-trivial way. For example

[TABLE]

We note that as $\Delta y\rightarrow 0$ ($\left\langle N\right\rangle\rightarrow 0$), the cumulants are dominated by $\left\langle N\right\rangle$ and the cumulant ratios go to unity which is the Poissonian limit, even in the presence of large correlations (given here by $R_{k}^{0}$).

A similar scaling of the factorial cumulants with the particle number, $\hat{C}_{k}\sim\left\langle N\right\rangle^{k}$ has been already obtained when we discussed baryon number conservation and volume fluctuations, see Eqs. (107) and (110). This is not unexpected, since both baryon conservation and volume fluctuation are, in the simple approach presented here, rapidity independent long-range phenomena.

The maximum value of $\Delta y$ presently available at RHIC is $\Delta y=1$. Now, we can fit the values of $R_{k}^{0}$ using the data for $\Delta y=1$ and calculate the rapidity dependence of the cumulants. This is discussed carefully in Ref. [346] and demonstrated in Fig. 41 for (a) $7.7$ GeV and (b) $19.6$ GeV. We see that the preliminary STAR data are consistent with a constant rapidity correlation. Of course, since the maximum rapidity interval covered by STAR is $|y|<0.5$, in practice this means that the correlation range in rapidity is larger than one unit, $\delta y>1$.

In Fig. 41 we also demonstrate the importance of various $\hat{C}_{n}$ contributions to $\kappa_{4}/\kappa_{2}$. Clearly the signal at $7.7$ GeV is driven by $\hat{C}_{4}$ and at $19.6$ GeV by $\hat{C}_{2}$, in agreement with Fig. 39.

As explained in Ref. [353], similar analysis can be done in the transverse direction. Here we need to take into account $\rho(p_{t})$ when calculating the dependence of the cumulant ratio on the size of the transverse momentum interval. Assuming that, within the current acceptance range,

[TABLE]

we again arrive at the conclusion that the cumulants and the factorial cumulants depend on the average number of measured protons $\left\langle N\right\rangle$ (in a given transverse momentum interval) as demonstrated in Fig. 42.

To summarize, the preliminary STAR data at low energies are consistent with $R_{k}(y_{1},...,y_{k};p_{t1},...,p_{tk})=const$ within the acceptance range of the STAR detector, $|y|<0.5$ and $0.4{\rm\,GeV}\leq p_{T}\leq 2{\rm\,GeV}$. However, baryon number conservation [353] as well as the arguments given in Section 4.12 suggest that the correlations will eventually change as one increase the rapidity coverage. For the upcoming second phase of the RHIC beam energy scan the rapidity coverage of STAR detector has been increased by a factor of 1.6. Thus, it will be interesting to see if and when the scaling presented in Figs. 41 and 42 is broken.383838We add that a mild repulsive rapidity dependence, studied in Ref. [353], is also consistent with the data at $7.7$ GeV.

7.1.5 Bimodal distribution

As already discussed in this Section, the surprisingly large factorial cumulants measured at $7.7$ GeV cannot be understood with baryon conservation or volume fluctuation. One can easily generate large values of certain factorial cumulants with, e.g., baryonic clusters but the signs of the measured $\hat{C}_{n}$ are rather challenging to obtain. Recently in Ref. [354] it was demonstrated that the preliminary STAR data at $7.7$ GeV can be naturally reproduced assuming that the proton multiplicity distribution has a bimodal form

[TABLE]

where $P_{(a)}$ and $P_{(b)}$ are the proton multiplicity distributions characterized by different means $\langle N_{(a)}\rangle$ and $\langle N_{(b)}\rangle$. It is interesting that one can reproduce the STAR data assuming that both $P_{(a)}$ and $P_{(b)}$ are characterized by small or even vanishing factorial cumulants. In the above equation $\alpha$ determines the relative strength of the two distributions. The STAR data (up to $\hat{C}_{4}$) can be reproduced assuming that $\alpha\approx 0.0033$, $\langle N_{(a)}\rangle\approx 40$ and $\langle N_{(b)}\rangle\approx 25.3$, as shown in panel (a) of Fig. 43. Here $P_{(a)}$ and $P_{(b)}$ are taken as binomial ($B=350$ and $p\approx 0.144$) and Poisson distributions, respectively. In panel (b) we show the same distribution with an imposed efficiency of $0.65$, which is more relevant for the STAR environment.

The bimodal distribution given by Eq. (117) and presented in Fig. 43 arises naturally when the system fluctuates between two distinct phases, as described in Section 4.2. A (small) system close to the critical point or the first order phase-transition is described by a similar bimodal distribution. This is demonstrated in Fig. 44 (see also Fig. 8), where the multiplicity distributions at various points in the phase diagram for the van der Waals model are presented (see [354] for details). If the two modes of the bimodal distribution indeed represent two distinct phases, one would expect this to be visible also in other observables, such as spectra, flow etc. As discussed in [354], this can be tested by selecting events with small or large proton number to enhance either of the two phases. Of course, there is also a more mundane explanation for the appearance of a bi-modal distribution: A small fraction of events can be, e.g., contaminated by some detector effects, which would result in (at least) two event classes. This scenario can be carefully studied and hopefully excluded.

It is imperative to test whether the underlying proton distribution at $7.7$ GeV is indeed described by the bimodal distribution. As shown in [354] and detailed in A, the factorial cumulants are approximately given by $\hat{C}_{n}\approx-\alpha\langle(N_{(a)}-N_{(b)})\rangle^{n}$ and thus are quickly growing with $n$ with alternating sign. Moreover, as found very recently in Ref. [355], the bimodal distribution describing the STAR data is statistically friendly, i.e., can be successfully measured through the factorial cumulants of high orders with a relatively small number of events. In this case, the factorial cumulants are driven by the bimodal structure away from the tails and thus can be measured with a limited number of events. This is in contrast to the Poisson, binomial, negative binomial, etc. distributions, which are statistically very demanding. In Ref. [355] the prediction up to $\hat{C}_{9}$ (with errors) are presented based on the current STAR statistics in $0-5\%$ central Au+Au collisions. Confirmation of these results could very well indicate the observation of the first-order phase transition.

7.2 Discussions on the measurements for anomalous chiral transport

In this subsection, we discuss the current understanding and existing issues of the various measurements for anomalous chiral transport. We also discuss the theoretical and experimental progress in addressing these issues.

7.2.1 The issues with background correlations

A main observable that has been used for the search of the chiral magnetic effect (CME) is the $\gamma$-correlator, Eq. (88). As discussed in Sec. 5.3, it is sensitive to the CME but also receives background contributions. Thus, the key question is what fraction of the observed signal arises from the CME signal compared to those originating from backgrounds. In terms of Eq. (90), it is a question of quantitatively separating the F and H components.

If the CME were to be dominant contribution, its extraction would have been easy. As it turns out, however, the measured $\gamma$-correlator from STAR at RHIC energy (with data shown in Sec. 6.4) is dominated by background contributions, as already pointed out in [201] shortly after the data came out. This situation can be recognized by a joint analysis of the $\gamma$ and $\delta$ correlators as defined in Eqs. (88) and (89), both of which were measured, as discussed in Sec. 5.3. Given these two correlators, one can define in-plane and out-of-plane projected azimuthal correlations: $\langle\cos\phi^{\alpha}\cos\phi^{\beta}\rangle=(\delta^{\alpha\beta}+\gamma^{\alpha\beta})/2$ and $\langle\sin\phi^{\alpha}\sin\phi^{\beta}\rangle=(\delta^{\alpha\beta}-\gamma^{\alpha\beta})/2$. These projections are shown in Fig. 45 for the STAR data at 200 GeV. In a scenario of CME dominance, one would expect the out of plane projection to be larger than the in-plan projection. Specifically, one expects same side, out-of plane correlation for same charge pairs, i.e. $\langle\sin\phi^{\pm}\sin\phi^{\pm}\rangle>0$, and opposite side out-of plane correlations for opposite charge pairs $\langle\sin\phi^{+}\sin\phi^{-}\rangle<0$. At the same time, the in-plane correlation should vanish, $\langle\cos\phi^{\alpha}\cos\phi^{\beta}\rangle\simeq 0$ for both same and opposite charge pairs. The data obviously deviates from such a pattern, in fact demonstrating a strong presence of background contributions. Indeed, the sizable $\langle\cos\phi^{\alpha}\cos\phi^{\beta}\rangle$ can only be background correlations which most likely would also make a sizable contribution toward the out-of-plane $\langle\sin\phi^{\alpha}\sin\phi^{\beta}\rangle$ correlation. This however does not exclude the possible presence of CME, and actually may still hint at certain CME contribution. For example, let us take a closer look at the same charge correlations: there appears a sizable negative $\langle\cos\phi^{\alpha}\cos\phi^{\beta}\rangle$ correlation in contrast to the $\langle\sin\phi^{\alpha}\sin\phi^{\beta}\rangle$ correlation that appears nearly zero and becomes slightly positive in the relatively peripheral region. Such a significant “mismatch” between in-plane and out-of-plane correlations could have two possible explanations. It could either be due to a negative (i.e. back-to-back type) correlation that somehow occurs only in-plane. However, at present it is not clear what mechanism may cause such a background. Or, it could be due to a positive out-of-plane-only correlation, such as the CME signal, on top of a back-to-back bulk background, which induces comparable negative contributions to both in-plane and out-of-plane correlations. In any case, the above discussion shows that the situation is rather involved, and the challenge is to extract a small signal immersed in a dominant background.

A number of specific sources of background correlations have been identified (see e.g. reviews [202, 120] for detailed discussions). We briefly discuss them here:

Transverse momentum conservation (TMC) — This leads to a back-to-back correlations. To see this, consider the extreme case of producing only two particles. Obviously, their momenta must be balanced in the transverse direction, $\vec{\bf p}$ and $-\vec{\bf p}$, resulting in a back-to-back correlation. For a system of N particles, this correlation is “diluted” by $1/N$. A quantitative analysis of TMC can be found in [204]. Of course, if the transverse momentum conservation is confined to certain local scale, then the effect would be amplified. Nevertheless, transverse momentum conservation does not distinguish between charges of particles and thus the correlations are the same for same-sign and opposite-sign pairs. Therefore, one can remove their influence by considering the difference between same and opposite charged correlations, $\gamma^{OS}-\gamma^{SS}$.

Local charge conservation (LCC) — This effect arises from possible charge neutrality over small size spatial domains upon freeze-out. Due to strong collective expansion, particles produced in close spatial proximity tend to have their momentum directions collimated on average. The neutrality enforces a companion negative charge for every positive charge produced from the same local domain. This leads to a near-side azimuthal correlations for opposite-sign pairs, resulting in positive contributions to both in-plane and out-of-plane correlations with the former slightly stronger. Quantitative analysis [206, 207, 208, 356] suggests that the LCC could account for a substantial amount of the observed opposite-sign correlation patterns. However the precise mechanism of such local neutrality is unclear and the spatial size of such neutral domains, which is a key parameter controlling the strength of this correlation, remains poorly constrained.

Resonance decay — Hadrons produced from resonance, or more generally cluster, decay will have their momenta correlated in azimuthal angle, inheriting the parent resonance’s momentum [205]. A most relevant example for the $\gamma$-correlator background is the neutral $\rho$ decay into charged pions, i.e. $\rho^{0}\to\pi^{+}\pi^{-}$, which makes a considerable contribution to a near-side correlation for opposite-sign hadron pairs and has similar characteristics to the LCC contribution. The hadronic decay contributions could be quantified by hybrid bulk evolution models which include a hadronic cascade stage, e.g., implemented via UrQMD simulations.

All the background sources discussed roughly scale inversely with the multiplicity. Therefore, they should increase as one decreases the centrality, i.e. considers more peripheral collisions. They all induce similar angular correlations for both in-plane and out-of-plane projected components, however with the in-plane component larger in magnitude due to stronger in-plane collective expansion. Therefore their contributions to the $\gamma$-correlator, being a difference between the in-plane and out-of-plane components, are roughly proportional to $v_{2}$ as in Eq. (90) and are absorbed together into the $F$-term there. Finally let us note, that centrality dependence of the same-sign out-of-plane correlation, left panel of Fig. 45, is much weaker as compared to the other correlations.

To wrap up this discussion, the current understanding is that the measured $\gamma$-correlator, while being sensitive to CME contributions, is dominated by strong non-CME background correlations but also bears hints of nonzero CME-like signal at RHIC energy region. To convincingly extract a possible CME signal, progress needs to be made in three directions. First, one needs to get the backgrounds under control by quantitatively understanding and computing them. Above we discussed the present status but obviously more needs to be done. Second, one needs to be able to quantitatively predict the CME signal. Here major progress has been and continues to be made, as we shall discuss next. Lastly, one needs to develop new experimental approaches that may be more sensitive to the CME and may allow suppression/separation/subtraction of the backgrounds: a lot of interesting developments are ongoing along this line, which we discuss in the last two parts of this subsection.

7.2.2 Quantitative modeling of anomalous chiral transport

To address the difficulty in the experimental detection of CME signals, it is imperative to develop a quantitative modeling framework that simulates the anomalous chiral transport while accurately accounts for the realistic environment in heavy ion collisions. In the past few years a matured framework, called the Anomalous-Viscous Fluid Dynamics (AVFD), has been developed [145, 146, 357]. We note in passing that there have been phenomenological study of CME based on kinetic transport models as well [358, 359, 163, 360, 361, 362, 363, 364].

The main idea of AVFD is to extend the hydrodynamic bulk evolution [365] to properly include the evolution of the fermion currents following the anomalous fluid dynamics framework described in Sec. 5.1.3. In AVFD the evolution of fermion currents is calculated on top of a viscous hydrodynamics background. This approach accounts for both anomalous and normal viscous transport effects, however, it neglects the feedback of the anomalous current onto the bulk-evolution. The most important feature of AVFD is the evolution of the ${\bf B}$-field driven anomalous current $J^{\mu,(1)}_{anomalous}$, Eqs. (81) and (82), where the left-handed and the right-handed components contribute with opposite sign. In Fig. 46 we demonstrate the effect of such chiral transport within AVFD. Specifically we consider the evolution of up-quark densities starting from a given initial density distribution. The initial up-quark density at time $\tau_{0}=0.60$ fm/c is depicted in panel-(a). In panels (b)-(d) we show the up-quark densities after the system has evolved to the time $\tau=3.00$ fm/c. In panel-(b) we see the density in absence of a magnetic field. In this case right-handed and left handed quarks show the same distribution. In panels (c) and (d) we show the distributions with a non-zero magnetic field ${\bf B}$ along positive y-axis for right-handed (panel (c)) and left-handed (panel(d)) up-quarks. The effect of the anomalous transport is clearly visible: the right-handed quarks are predominantly transported along the direction of the magnetic field, i.e. towards the positive $y$, while the left-handed quarks are transported in the opposite direction. For the negatively charged down-quarks one would see just the opposite pattern: Right-handed d-quarks will predominantly propagate towards negative $y$ while left-handed ones propagate to positive $y$.

Within the AVFD framework in a heavy ion collision the charge separation arises then as follows. The $\vec{B}$ field points along the out-of-plane direction, which we will denote by the $y$-direction. If initially we have a nonzero axial charge density, i.e. an imbalance between right-handed and left-handed quarks, we will have a situation shown similar to those shown in panels (c) or (d) of Fig. 46. For example, if we have more right-handed then left-handed quarks, the positively-charged up-quark density would move more towards the positive $y$-direction while, at the same time, the negatively charged d-quark density would move towards the negative y-direction. As a result we have a net charge current parallel to the magnetic field. If, on the other hand we have more left-handed than right-handed quarks, the charge current would be anti-parallel to the magnetic field. These currents cause a separation of electric charge across the reaction plane which can be quantified by the coefficient in Eq. (87). AVFD, thus, allows for a quantitative understanding of the generation and evolution of a CME-induced charge separation signal during the hydrodynamic stage as well as its dependence on various theoretical ingredients (for detail see [146]). With reasonable estimates of key parameters, such as the strength and duration of the magnetic field, the results of the AVFD model for the CME-induced values of the H-correlator, defined in Eq. (90), are shown in Fig. 47 and found to be consistent with those extracted from the STAR measurements.

More recently this framework has been further developed into event-by-event simulations (EBE-AVFD) [357], which is the state-of-the-art tool for quantifying anomalous chiral transport in heavy ion collisions. The EBE-AVFD includes fluctuating initial conditions as well as a hadronic cascade stage via UrQMD after the hydrodynamic evolution which accounts for the influence of hadronic re-scatterings and resonance decays. Further improvements of the framework shall address a number of issues, including more accurate determination of axial charge initial conditions [62, 63], the fluctuation and dissipation of axial charge during the hydrodynamic evolution [366, 367, 368, 369], as well as the possible generation of the CME current in the pre-equilibrium stage [164].

7.2.3 Developing new measurement methods

Recognizing the limitations of the $\gamma$-correlator, it is natural to look for alternative measurement methods that would help extracting possible CME signals. As can be seen from the two-component decomposition in Eq. (90), a useful strategy is to “dial” certain external control conditions so as to vary one component while keep the other component stay constant, or to vary the two components in opposite ways, i.e. allowing one to increase while the other to decrease. A number of new approaches were put forward in recent years, all of which follow more or less that philosophy. As we discussed in Section 5.3, the signal is driven by the magnetic field whereas the background is driven by the elliptic flow, $v_{2}$. Therefore, in order to disentangle signal and background, it would be helpful to vary the magnetic field while keeping the elliptic flow constant or vice versa. The event-shape-engineering approach, already discussed in Sec. 6.4, is an attempt to vary the elliptic flow, while keeping the magnetic field more or less constant [342, 343]. Alternatively, one may want to keep the elliptic flow constant and change the magnetic field. This is the idea behind the isobar run which we discuss at the end of this subsection. Let us turn to the event-shape-engineering approach first. Within a given centrality class of events, the background contributions can be considerably varied by selecting the events’ observed elliptic flow coefficient $v_{2}$. This could allow a subtraction of the $v_{2}$-dependent portion from the inclusive $\gamma$-correlator and, therefore, provide a potential CME signal that is independent of $v_{2}$. The working assumption of this method is that the magnetic-field-driven CME signal is uncorrelated with the geometric shape fluctuations and stays constant across the whole centrality class. This however is not entirely true and the magnetic field component projected to the event-plane (which is the relevant component for measured CME signal) would actually change with the event selection. There is thus the danger of over-subtraction in the $\gamma$-correlator. There could also be issues with potential selection bias and event plane resolution.

A different approach is to compare the $\gamma$-correlator measurements with respect to the event plane (EP) and that with respect to the reaction plane (RP) [370, 371]. Event-by-event fluctuations in the initial conditions cause a de-correlation in azimuthal angles among event-wise elliptic anisotropy, the magnetic field orientation and the reaction plane. The event-plane, by definition is the direction maximizing the elliptic anisotropy $v_{2}$ and therefore the EP-measured $v_{2}$ value is greater than that from RP method. On the other hand, the magnetic field direction is more correlated with RP than with EP. Thus the RP-projected magnetic field value is found to be greater than the EP-projected value. Therefore, the $v_{2}$-driven background and the magnetic-field-driven signal would vary differently between EP and RP based measurements. This would allow the possibility of inferring the signal/background ratio by contrasting EP-measured and RP-measured $\gamma$-correlator. A main caveat with this approach is that reaction plane is a theoretical concept and can only be approximately inferred, e.g., from the so-called ZDC plane at RHIC. This uncertainty requires careful study.

Another category of new measurements focuses on characterizing the background components by examining the correlations in a situation where the backgrounds will be present while the CME contribution is expected to be absent. For example, one can measure the $\gamma$-correlator in a small colliding system (e.g. pPb at LHC or dAu at RHIC) in which the magnetic field is decorrelated with event plane and thus the CME does not contribute to the $\gamma$-correlator. Or, one can measure a modified $\gamma$-correlator that is not defined with respect to the elliptic event plane but to the triangular event plane: in this case again the magnetic field is decorrelated with event plane and, again, the CME does not contribute to the $\gamma$-correlator. In both cases, the various background correlations should still be present and thus can be determined from measurements. These approaches were developed by CMS and successfully used to put a rather stringent limit on the existence of CME signal in PbPb collisions at $5.02$ TeV [199, 200]. Recently, the STAR collaboration also applied these analysis methods for constraining background contributions at RHIC.

Yet another strategy aims to find observables in which the background correlations could be significantly suppressed. For example, one may apply an invariant mass selection on the $\pi^{+}\pi^{-}$ pairs which dominate the opposite-sign pairs used in the $\gamma$-correlator measurements [372, 371]. The aforementioned resonance decay contributions, e.g. $\rho^{0}\to\pi^{+}\pi^{-}$ could be suppressed by an invariant mass selection away from the $\rho^{0}$ mass region. This method helps control resonance decay contributions but can not suppress other background correlations. The CME signal would be mostly produced in the soft regime with low invariant mass, where recent measurements suggest a nonzero signal.

There is also a new proposal of the so-called $R$-correlator method [373, 374]. To construct this observable, one first examines the distribution $N(\Delta S_{\parallel})$ of the event-wise out-of-plane dipole charge separation $\Delta S_{\parallel}$ for a given set of events. The next step is to construct a ratio $C(\Delta S_{\parallel})$ via normalizing $N(\Delta S_{\parallel})$ by a baseline distribution $N_{shuffled}$ obtained from the same set of events after random reassignment of the charge of each particle in an event. This ratio is to be further contrasted with a similarly constructed ratio for the in-plane separation $C(\Delta S_{\perp})$, leading to the $R$-correlator: $R(\Delta S)=C(\Delta S_{\parallel})/C(\Delta S_{\perp})$. Since CME-induced charge separation will only influence the out-of-plane distribution, the $R$-correlator is expected to be sensitive to the CME. Indeed AVFD simulations unambiguously demonstrate such sensitivity, with the so-constructed $R(\Delta S)$ distribution becoming more and more concave-shaped when the CME signal is increased. The step of normalizing the charge separation distributions by shuffled baseline distributions also appears to effectively suppress background contributions to certain extent. A preliminary STAR analysis [375] based on this method shows that at comparable multiplicity, a rather flat distribution of the $R$-correlator for p+Au and d+Au collisions (where pure backgrounds prevail) while a strongly concave-shape distribution is seen for peripheral Au+Au collisions, providing an indication of CME signal. At present, potential issues with this method mainly concern the lack of a full understanding on the behavior (e.g. convex versus concave shape) of the $R$-correlator in response to various backgrounds [376] as well as the influence of factors like kinematic cuts, multiplicity fluctuations and event-plane resolution. Finally we note in passing a newly proposed observable to look for charge separation via signed balance function based on the idea of examining the momentum ordering of charged pairs along the out-of-plane direction [377].

Clearly each of the new methods discussed above, just like the $\gamma$-correlator, has its own advantages and drawbacks. The best strategy to move forward is perhaps to look for a consistent pattern regarding the presence of a CME signal by a global analysis of the various methods and observables. A first attempt at such a global analysis was recently carried out by the STAR collaboration and reported at the Quark Matter 2018 conference [375]. The STAR results suggest that a CME signal could likely be at the 10% level of the measured inclusive $\gamma$-correlator. However, the currently large uncertainty as well as the discrepancy among different measurement methods would not rule out a signal in the range of $0\sim 20$%. A rigorous conclusion, therefore, requires future experimental and theoretical effort to fully understand and improve each approach’s methodology, to carefully examine and address potential issues with them, and to meaningfully combine all analyses together.

7.2.4 The isobaric collisions

As can be seen from the above discussions, while new proposals for measurement methods are under active development, each of them may face difficulties of one sort or another. The theoretical uncertainties and experimental limitations have made it very challenging to quantitatively separate the possible signal from the backgrounds. In the past few years, a novel idea emerged and matured, based on contrasting collisions of the isobaric nuclei pairs [345, 378]. In one sentence, this idea could be summarized as follows: the bulk collective dynamics controls the background correlations, while the magnetic field controls the CME signal, and the strategy is to create and compare two colliding systems that have identical bulk properties but different magnetic field.

Let us explain how this idea works in some details. A pair of isobars are two different nuclei that have the same mass number A but different electric charge (i.e. proton number) Z. Specifically for our discussions, we consider Ruthenium (Ru) and Zirconium (Zr), both having A=96 nucleons while Z=44 for Ru and Z=40 for Zr. One can then perform experiments with Ru+Ru collisions and with Zr+Zr collisions at the same beam energy, and look for the difference of desired observables between the two colliding systems. The 10% difference in the electric charge implies a difference in the magnetic field of also about 10%. This can also be seen from the computed magnetic field values with event-by-event simulations in Fig. 48 (left). Since the CME current is proportional to the magnetic field, see Eq. (72), and since the $\gamma$-correlator measures the variance of the CME current, the $\gamma$-correlator scales with the square of the magnetic field $\gamma\sim B^{2}$. Therefore, a 10% difference in the magnetic field translates into a 20% difference in the CME-induced contribution to the $\gamma$-correlator. This can also be seen in the results from the AVFD simulation shown in the middle panel of Fig. 48.

The key for success of the isobar run is ensure that the background correlations will be exactly the same. This requires that one compares the $\gamma$-correlator for events from Ru+Ru and Zr+Zr, which have identical bulk properties, in particular multiplicity and elliptic flow $v_{2}$. Since the two isobars have the same number of nucleons, naively one would expect that a simple centrality (or multiplicity) selection would ensure this requirement. However, uncertainties in the initial nucleon distributions in the two isobar nuclei may actually bring percent-level difference in the elliptic flow between the two systems [379] for a given centrality. Sine the CME signal in the $\gamma$-correlator is expected to be rather small, this presents a serious problem which needs to be controlled. As shown in [357], this problem can be overcome by selecting events on both multiplicity and elliptic-flow. Event-by-event simulations [357] for both initial conditions and for the final state observables have demonstrated that the joint-multiplicity-ellipticity cut would be effective and sufficient for this purpose.

Experimentally one compares the measured $\gamma$-correlator between Ru+Ru and Zr+Zr, for which the signal differs by about 20% while the backgrounds are (ideally speaking) identical. The actual difference in the $\gamma$-correlator depends of course on the background level in $\gamma$-correlator. This is shown in the right panel of Fig. 48, where we show the results of calculations Refs. [345, 360] for the expected differentiation capability with 400-million events (vertical axis) of the contrast measurement versus the background level in $\gamma$-correlator (horizontal axis) for the isobaric collisions at RHIC. For example, if the background level were to be less than 2/3 of the $\gamma$-correlator, then the isobar experiment would allow an observation of CME signal at $5\sigma$ level or beyond with just 400-million events.

The isobaric collision experiment was successfully completed in the RHIC 2018 with more than twice the planned 400-million events collected by the STAR detector. With this improved statistics and with the latest estimate of from the STAR analysis of a CM signal of about 10%, one could hope for a signal of at least $2\sigma$ significance. Furthermore, there are now an array of new observables besides just the $\gamma$-correlator, as discussed in the preceding part. One may therefore still anticipate a fairly reasonable chance to have a statistically significant statement on the presence of CME signal based on a synthetic analysis of multiple observables.

Clearly, the isobaric collision experiment provides a unique opportunity to detect the possible presence of CME in heavy ion collisions. A conclusive observation of the CME signal would directly establish, for the first time, the existence of chirality-flipping topological transitions that stem from the vacuum structure of QCD as well as provide evidence for the restoration of chiral symmetry in the high temperature quark-gluon plasma. Therefore, given the characteristic beam energy dependence of the CME signal as discussed in Sec. 5.4, a positive CME signal from the isobar run at $200$ GeV would strongly suggest future measurements at lower beam energies. A measurement in the intermediate BES energy region would possibly detect a maximal CME signal while one in lowest BES energy region may help establish the expected disappearance of the CME signal.

8 Summary and Outlook

We close this review by summarizing the present status as well as discussing open questions. We then finish with discussions of experimental upgrades at current experimental facilities as well as future opportunities at upcoming new facilities.

8.1 Present status

The measurements taken during the first phase of the RHIC Beam Energy Scan (BES-I) have provided rather intriguing results, with hints at potentially very interesting physics.

Global observables such as elliptic flow and particle yields suggest that the systems created in the energy range of the RHIC BES are (locally) equilibrated and that their bulk dynamics is controlled by the collision geometry and collective expansion. 2. 2.

Several observables, however, exhibit quite an intriguing energy dependence below $\sqrt{s_{NN}}\lesssim 20{\rm\,GeV}$. The directed flow as well as the higher order cumulants ($\kappa_{3}$ and $\kappa_{4}$) of net-protons show a non-monotonic behavior which may indicate dramatic changes in the equation of state which, for example, may be due to a critical point or first-order phase transition. In addition femtoscopy measurements seem to indicate a maximum emission time and minimum expansion velocity at the same beam energy (about $\sqrt{s_{NN}}\simeq 20{\rm\,GeV}$), suggesting that the system might cross the softest point of the EoS at this energy. A closer analysis of the measured net-proton cumulants finds that below $\sqrt{s_{NN}}\lesssim 15{\rm\,GeV}$ the system exhibits significant four-particle correlations which increase in strength rapidly with decreasing collision energy. Model calculations, such as UrQMD, as well as general arguments based on baryon number conservation predict the cumulants (within fixed rapidity window $\Delta y=1$) to decrease toward lower beam energy. 3. 3.

The charge asymmetry of the azimuthal correlator, $\gamma$, which is sensitive to the anomalous chiral transport current, seems to show features in line with expectations from the Chiral Magnetic Effect at high collision energy while it vanishes for energies below $\sqrt{s_{NN}}\lesssim 10{\rm\,GeV}$. This may be an indication that at this energy the system created in these collisions changes from one dominantly in the phase with restored chiral symmetry to one dominantly in the phase with spontaneously broken chiral symmetry. Similarly, analyses of the charge asymmetry by means of the H-correlator as well as the slope parameter of the charge-dependent elliptic flow, which is sensitive to the Chiral Magnetic Wave, show a positive signal at high energy which disappears for energies below $\sqrt{s_{NN}}\lesssim 10{\rm\,GeV}$.

We should of course point out that the above findings are partially based on preliminary data and may depend upon model inputs in certain aspects. Also the presently available statistics, especially for the cumulant measurements as well as those for the charge dependent correlator $\gamma$, is rather limited. Much better statistics would be required in order to draw any firm conclusions. This situation will improve dramatically in the second phase of the RHIC beam energy scan (BES-II). Aside from the limited statistics, there are a number of open questions which need to be addressed.

8.2 Open questions

Criticality and net-proton cumulants:

• (a)

The critical point signatures described in Section 4 are characterized by specific patterns of non-monotonous behavior of various cumulants some of which (such as $\kappa_{4}$) appear to be intriguingly matched by experimental data described in Section 6. Before drawing conclusions, it is important to keep in mind that the predictions in Section 4 are based on an assumption of approximate static equilibrium. This assumption has to be relaxed in order to address the effects of the expansion dynamics on critical point signatures, some of which, such as memory effects, are anticipated. The hydrodynamic framework incorporating fluctuations and critical slowing down needs to be developed and compared to experimental measurements.

• (b)

The proton cumulant measurements at the lowest BES energies show that the system exhibits large and positive four-particle correlations and sizable three-particle correlations which are negative. Such a pattern cannot be explained by simply invoking attractive interactions among the protons. Similarly, the directed flow, $v_{1}$, does not support a simple picture based on attractive interactions either. This raises the question of what dynamics can explain both the net-proton cumulants and the directed flow. Can the observed patterns be related to a change in the equation of state and if not, what other dynamical effect plays a role? Clearly the steep increase of the fourth-order net-proton cumulant ratio $\kappa_{4}/\kappa_{2}$ calls for measurements at even lower energies or equivalently higher baryon densities. In addition, of course, higher statistics for all energies below $\sqrt{s_{NN}}\lesssim 20{\rm\,GeV}$ are needed.

• (c)

While not discussed in detail in this review, it would be equally interesting to experimentally explore and establish the existence of a crossover transition at $\mu_{B}\simeq 0$. As discussed in Section 4.8, to achieve this goal future high energy collider experiments should pursue the measurement of cumulants of the net-proton distribution to sixth, $\kappa_{6}/\kappa_{2}$, eigth, $\kappa_{8}/\kappa_{2}$, or even higher order. Naturally, for these measurements to be successful high statistics data sets as well as excellent control of the systematics are mandatory. 2. 2.

Chirality and chiral magnetic effect (CME):

• (a)

A crucial but missing piece of information for a quantitative understanding of the CME is the time evolution of the in-medium magnetic field. Can it be quantitatively computed including its collision energy dependence? Can it be connected with the fluid rotation in some way? Can the lifetime of the magnetic field be quantitatively connected with certain other experimental observables, such as the polarization difference between particles and antiparticles or the di-electron transverse momentum broadening?

• (b)

A major obstacle in the search for the Chiral Magnetic Effect is the characterization of non-CME physics backgrounds and not all of them are well understood. Could combining various experimental methods together help constrain, extract, or remove the backgrounds? Can background effects be quantified and implemented into dynamical models so that they are able to predict also the energy dependence of both the signal and expected backgrounds? What additional theoretical developments and measurements are required to achieve this?

• (c)

At present, the most promising experimental strategy to isolate the CME signal is the isobar run (see Sec. 7.2.4) which has just been completed at RHIC where data have been taken at a collision energy of $\sqrt{s_{NN}}=200{\rm\,GeV}$. If the analysis of the data show a positive CME signal, a definitive discovery will require additional measurements at lower beam energies. These measurements are also necessary to unambiguously establish that chiral symmetry is actually restored in the system created in these collisions.

8.3 Theory efforts

To address these questions one needs not only additional measurements but also progress in theory. In particular, a definitive answer to the questions of the existence of the QCD critical point and anomalous transport requires a coherent theory effort aiming at a quantitative understanding of the dynamical evolution of the system created in a heavy ion collision. Essential to this is the development of a quantitative dynamical framework which describes the entire evolution of the system from the initial collision to the final state. Such a framework should build on the progress achieved over the last decade which resulted in a quantitative description of heavy ion collisions at the highest RHIC and LHC energies. However to address the physics discussed in this review, several important further developments are needed:

The initial conditions for the hydrodynamic evolution are substantially different at lower energies, where the penetration time of the nuclei could be potentially as long as the equilibration time. Also the stopping of the incoming nucleons needs to be understood very well as they are the source of net-baryon density at midrapidity. Some progress in this direction has already been made [287].

• 2.

The (viscous) hydrodynamic evolution needs to be extended to propagate all conserved currents as well as the anomalous current. In addition, it needs to account for and propagate the fluctuations which occur in the vicinity of the critical point. Also, a definitive calculation of the CME will require magneto-hydrodynamics to properly determine the time evolution of the magnetic field and its impact on the anomalous currents. And, as already discussed in Section 4.7, an equation of state with critical point which is consistent with lattice QCD results at $\mu_{B}=0$ is needed [88].

Considerable progress has been made on many fronts: First results for the the propagation of conserved and anomalous currents have been obtained [287, 145, 146] (see also Sect.7.2.2). The treatment of fluctuations within stochastic hydrodynamics [380, 381, 382, 383, 384] and their coupling to critical modes [385, 386, 387] have been studied. Alternatively, the time evolution of hydrodynamic fluctuations by means of deterministic propagation of correlations functions have been explored [388, 389], including the self-consistent description of critical fluctuations and critical slowing down [105, 106]. First calculations have applied magneto-hydrodynamics to study the time evolution of the magnetic field [186, 390].

• 3.

A transformation of the hydrodynamic fields to actual particles needs to preserve the fluctuations and all conserved and anomalous currents [391]. In addition, as the critical point and phase transition affect also the hadronic phase, the kinetic evolution of the hadronic phase needs to account for this [117].

All these essential pieces still need to be put together in one common framework which can be used to perform a quantitative comparison with experimental data. Clearly there are many open conceptual as well as technical problems and an ever growing community of theorists worldwide is actively working to address these challenges.

8.4 Future Upgrades and Physics Program at RHIC and the global picture

During the years of 2019 - 2021, the STAR experiment at RHIC will undertake the second beam energy scan campaign (BES-II), with the goal to dramatically reduce both statistical and systematical uncertainties for all of the observables and to map the QCD phase diagram with highest possible precision. As shown in Table 1, BES-II will cover the energy range $7.7{\rm\,GeV}\leq\sqrt{s_{NN}}\leq 19.6{\rm\,GeV}$ corresponding to $420{\rm\,MeV}\geq\mu_{B}\geq 200{\rm\,MeV}$ [392]. The expected statistical uncertainties for the observables relevant for the search for anomalous transport and the QCD critical point are illustrated in Fig. 30 (green band) for the cumulants and in Fig. 37 (pink band) for the $\gamma$ correlator. In addition to the increased statistics STAR detector upgrades such as that of the inner TPC, the end-cap time of flight detector, and the event plane detector allow for increased rapidity coverage, better event plane resolution, and centrality selection and particle identification at forward rapidities. These new capabilities will enable STAR to address questions concerning the range of the proton correlations, more precise measurements of the $\gamma$ correlator, and a suppression of auto-correlations for the fluctuation measurements.

Both the fourth-order proton cumulants, which have the largest value at the lowest available RHIC collider energy, as well as the CME observable, $\gamma^{OS}-\gamma^{SS}$, which vanishes right at the lowest energy, call for measurements at even lower energies. To address this issue, the STAR collaboration has developed a fixed-target (FXT) program, which uses a gold-target at the one end of the TPC accompanied with a time of flight wall at the other end of the TPC (see Fig. 49). With the fixed target setup, STAR will be able to measure particle productions and correlations in Au+Au collisions in the energy range $3{\rm\,GeV}\leq\sqrt{s_{NN}}\leq 19.6{\rm\,GeV}$, extending maximum accessible value of the chemical potential from $\mu_{B}\simeq 400{\rm\,MeV}$ to $\mu_{B}\simeq 750{\rm\,MeV}$. Importantly for the control of systematics, the highest center-of-mass energy of the FXT mode overlaps with the lowest energy from colliding mode at $\sqrt{s_{NN}}=7.7{\rm\,GeV}$, while its lower range overlaps with the collision energies explored by future experiments, such as CBM at FAIR [393] and MPD at NICA.

The exciting physics presented in this review has sparked many new experimental initiatives worldwide. In addition to BESII and the STAR FXT program, there are a number of approved experiments coming online in the near future, which will be able to address the physics discussed here. They are the compressed baryonic matter (CBM) experiment at the FAIR facility, the multi-purpose-detector (MPD) at NICA and the CSR-external target experiment (CEE) at HIAF. In addition to those there are plans to have an experiment at J-PARC in Japan. And, last but not least, there is the HADES experiment at GSI which has already taken data for Au+Au collisions at $\sqrt{s_{NN}}\simeq 2.2{\rm\,GeV}$ and is presently analyzing proton cumulants which, if successful, will provide and important low energy reference point for the fluctuation measurements. The energy range together with the (expected) interaction rate of the various facilities are shown Fig. 50, which demonstrate the worldwide excitement about the physics presented in this review.

Acknowledgements

This material is based upon work partly supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract number DE-AC02-05CH11231, contract number DE-AC03-76SF00098, contract number DE-FG0201ER41195, as well as within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration. The work is also supported in part by MOST of China 973-project No. 2015CB856900 and NSFC Grant No. 11735007, by Ito Science Foundation (2017) and JSPS KAK-ENHI Grant No. 25105504, as well as by NSF Grant No. PHY-1913729. AB also acknowledges partial support by the Polish Ministry of Science and Higher Education and the National Science Centre Grant No. 2018/30/Q/ST2/00101. JL also acknowledges hospitality during sabbatical visits at CCNU and SCNU where part of this work was completed. The authors thank Dmitri Kharzeev, Swagato Mukherjee, Shuzhe Shi and Paul Sorensen for useful discussions and communications.

Appendix A Introduction to Cumulants

In this Appendix we present a general discussion of the multiplicity distribution $P(N)$ and its various characteristics. $P(N)$ is the probability to produce $N$ particles and by definition $0\leq P(N)\leq 1$, and $\sum_{N}P(N)=1$. For example $P(N)$ could describe the probability to produce $N$ protons in Au+Au collisions. In this review we discuss several distributions: the Poisson distribution, the binomial distribution and the negative binomial distribution. Here we mostly consider a simplified case where $P(N)$ depends on one variable only but we also generalize to the higher number of variables.

Factorial moments

It is convenient to define the generating function

[TABLE]

where the value of $H(z)$ at $z=1$ is determined by the normalization condition $\sum_{N}P(N)=1$. Having $H(z)$ we can readily calculate the average number of particles, the average number of pairs, triplets etc. In general, it is easy to see that

[TABLE]

where $F_{k}$ are called the factorial moments and $k\geq 1$. $F_{1}=\left\langle N\right\rangle$ is the average number of particles, $F_{2}=\left\langle N(N-1)\right\rangle$ is the average number of pairs, $F_{3}=\left\langle N(N-1)(N-2)\right\rangle$ is the average number of triplets etc. Let us add here that $F_{k}$ are directly related to typically measured multi-particle rapidity393939In this Section we will usually discuss rapidity distributions but our discussion is general and applies to any other variable or variables. (or any other variable) distributions, see, e.g., Ref. [111]

[TABLE]

where $\rho_{1}(y)$ is a single-particle rapidity distribution404040In this review we will usually denote it by $\rho(y)$., $\rho_{2}(y_{1},y_{2})$ is the two-particle density414141To be more precise, $\rho_{2}(y_{1},y_{2})dy_{1}dy_{2}$ is the number of pairs of particles such that one particle is located in a rapidity interval $y_{1}-dy_{1}/2<y<y_{1}+dy_{1}/2$ and another one is located in $y_{2}-dy_{2}/2<y<y_{2}+dy_{2}/2$. (distribution of pairs) etc.

In the case of two species of particles we have

[TABLE]

where $P(N,\bar{N})$ is, e.g., the probability to produce $N$ protons and $\bar{N}$ antiprotons. Here the factorial moments are given by

[TABLE]

where $F_{1,1}=$ $\left\langle N\bar{N}\right\rangle$, $F_{2,1}=$ $\left\langle N(N-1)\bar{N}\right\rangle$ etc. If $P(N,\bar{N})=P(N)P(\bar{N})$, that is, the number of, say, protons is independent from the number of, say, antiprotons we obtain

[TABLE]

and $F_{i,k}=F_{i}F_{k}$. As we will discuss later, the factorial moments are very useful in experimental analysis since they are easy to correct for detector efficiency [329, 318].

Factorial cumulants

Another useful characteristics of $P(N)$ is given by the factorial cumulants, which are generated from the following function

[TABLE]

and the factorial cumulants are given by

[TABLE]

It is instructive to list the first three $\hat{C}_{k}$. $\hat{C}_{1}=\left\langle N\right\rangle$ is the average number of particles and

[TABLE]

[TABLE]

where we expressed $\hat{C}_{2}$ and $\hat{C}_{3}$ through the factorial moments and $\hat{C}_{2}$. In general it is straightforward to express $\hat{C}_{k}$ through $F_{i}$ by taking appropriate derivatives of $\ln\left(H(z)\right)$, see Eq. (124), and using Eq. (119).

We note that $\hat{C}_{k}$ are directly related to the integrals of the multi-particle correlation functions

[TABLE]

where $C_{k}(y_{1},\ldots,y_{k})$ is the $k$-particle genuine correlation function. The integration is performed over the same region of phase-space where $P(N)$ is measured. For example, the two-particle rapidity density

[TABLE]

is given by a product of two single-particle distributions, $\rho(y_{1})\rho(y_{2})$, plus the two-particle correlation function $C_{2}(y_{1},y_{2})$. Taking the integral of Eq. (129) over a given interval in rapidity we obtain

[TABLE]

which is the same as Eq. (126). Analogously, for three particles we have [394] (see also, e.g., Ref. [395])

[TABLE]

where $C_{3}(y_{1},y_{2},y_{3})$ is the three-particle genuine correlation function in rapidity. Calculating the integral over a given bin in rapidity we obtain

[TABLE]

which is the same as Eq. (127). In Ref. [395] explicit formulas for $\rho_{k}(y_{1},...,y_{k})$ up to $k=6$ are presented.

The factorial cumulants are very useful in analysis of physical processes since they are directly related to correlations present in the system and the order of the factorial cumulant is related to the genuine correlation of the same order, see, e.g., [110, 346, 396]. This is not the case for the factorial moments, where, e.g., as seen from Eq. (132), $F_{3}$ is a combination of factorial cumulants of different orders.

In the case of two species of particles we have

[TABLE]

and

[TABLE]

where $C_{n+m}(...)$ is the $n+m$ genuine correlation function involving $n$ and $m$ particles of different species. For example

[TABLE]

which is zero if, e.g., the number of protons $N$ is independent from the number of antiprotons $\bar{N}$.

Again, the mixed factorial cumulants $\hat{C}_{n+m}^{(n,m)}$ are directly related to the correlations present in the system. Suppose that $P(N,\bar{N})=P(N)P(\bar{N})$, that is the number of, say, protons is independent from the number of, say, antiprotons. In this case

[TABLE]

and according to Eq. (134), the mixed factorial cumulants $\hat{C}_{n+m}^{(n,m)}=0$ for $n\geq 1$ and $m\geq 1$.

Cumulants

The cumulants of a given multiplicity distribution $P(N)$ are generated from the cumulant generating function $K(t)$ given as

[TABLE]

and the cumulants read

[TABLE]

As already mentioned in Section 6.3.1, the STAR Collaboration measured cumulants of net-proton [320], net-charge [313] and net-kaon [321] distributions in Au+Au collisions, see also Refs. [323, 318]. The definition of $K(t,\bar{t})$ for two species of particles is analogous to Eq. (133), where we sum over $P(N,\bar{N})e^{tN}e^{\bar{t}\bar{N}}$. The net-proton cumulants, for example, are given by

[TABLE]

where the corresponding $K(t,\bar{t})$ is taken at $\bar{t}=-t$. In this case the cumulants are expressed through the moments of $N-\bar{N}$ distribution. Clearly, one could define also the cumulants of the sum $N+\bar{N}$ by taking appropriate derivatives of $K(t,+t)$ at $t=0$.

We note here that if $P(N,\bar{N})=P(N)P(\bar{N})$, that is, the number of, say, protons, is independent from the number of, say, antiprotons we obtain

[TABLE]

an the cumulants of, say, net-proton number can be expressed as a sum of proton cumulants plus (even cumulants) or minus (odd cumulants) antiproton cumulants. Performing straightforward derivatives of $K(t,-t)$ we obtain the net-proton (or baryon, charge, strangeness, etc.) cumulants

[TABLE]

where $\delta[N-\bar{N}]=N-\bar{N}-\left\langle N-\bar{N}\right\rangle$. As discussed in Section 4.8, the cumulants appear naturally in statistical physics as derivatives of the partition function.

It is useful to express the cumulants through the factorial cumulants $\hat{C}_{k}$. Let us start with one species of particles. In this case $K(t)=G(x(t))$, where $x(t)=e^{t}$, see Eq. (137), and the second order cumulant reads

[TABLE]

where, according to Eq. (125), $dG(x)/dx$ and $d^{2}G(x)/dx^{2}$ at $x=1$ equal $\hat{C}_{1}=\left\langle N\right\rangle$ and $\hat{C}_{2}$, respectively. Clearly, $dx/dt=e^{t}$ which is one at $t=0$. Performing straightforward calculations we obtain424242Additionally: $\kappa_{5}=\left\langle N\right\rangle+15\hat{C}_{2}+25\hat{C}_{3}+10\hat{C}_{4}+\hat{C}_{5}$ and $\kappa_{6}=\left\langle N\right\rangle+31\hat{C}_{2}+90\hat{C}_{3}+65\hat{C}_{4}+15\hat{C}_{5}+\hat{C}_{6}$. [110, 346]

[TABLE]

For two species of particles the formulas are a bit more complicated (see the Appendix of Ref. [346])

[TABLE]

which obviously reduces to Eq. (143) if $\bar{N}=0$ since in this case only $C_{k}^{(k,0)}$ may be different than zero. Here $\hat{C}_{2}^{(1,1)}$ measures the correlations between, say, protons and antiprotons, $\hat{C}_{2}^{(0,2)}$ is the factorial cumulant for antiprotons only etc. It is clear that the cumulants for the net-protons only are not enough to determine all $\hat{C}_{i+k}^{(i,k)}$. As seen from Eqs. (143) and (144), the cumulants mix the factorial cumulants of different orders and thus their interpretation is not straightforward. We will come back to this problem shortly.

At the end of this part, let us clarify the notation frequently used in, e.g., experimental plots, see for example Fig. 31. The standard deviation squared, $\sigma^{2}$, is simply the second central moment $\langle(x-\langle x\rangle)^{2}\rangle$, where $x$ equals to $N$ or $N-\bar{N}$. Clearly $\sigma^{2}=\kappa_{2}$. The skewness, $S$, is defined as $\langle(x-\langle x\rangle)^{3}\rangle/\sigma^{3}$ and thus $S\sigma$ equals to $\langle(x-\langle x\rangle)^{3}\rangle/\sigma^{2}=\kappa_{3}/\kappa_{2}$. The kurtosis or, more precisely, the excess kurtosis, $\kappa$ is defined as $\kappa=\langle(x-\langle x\rangle)^{4}\rangle/\sigma^{4}-3$ and therefore $\kappa\sigma^{2}=\langle(x-\langle x\rangle)^{4}\rangle/\sigma^{2}-3\sigma^{2}=\kappa_{4}/\kappa_{2}$, see Eq. (141). 434343 For the normal distribution, $\langle(x-\langle x\rangle)^{4}\rangle/\sigma^{4}=3$ which provides a kind of baseline and therefore commonly used to shift the kurtosis value for defining the excess kurtosis.

Discussion

Here we discuss the factorial moments, factorial cumulants and cumulants in different physical situations. Let us start with one species of particles distributed according to the Poisson distribution

[TABLE]

In this case, the generating function, given by Eq. (118), reads

[TABLE]

and consequently, see Eq. (119), we obtain

[TABLE]

that is, the number of pairs, triplets, etc., is simply expressed by a product of $\left\langle N\right\rangle$. This result corresponds to the situation where multiparticle rapidity distribution, $\rho_{k}$ for each $k$, is given by a product of single particle distributions

[TABLE]

that is, the multiparticle genuine correlation functions $C_{k}(y_{1},...,y_{k})=0$ for each $k$. Indeed, performing integration over a given bin in acceptance, see Eq. (120), we obtain Eq. (147).444444Strictly speaking, $F_{k}=\left\langle N\right\rangle^{k}$ means that the factorial cumulants $\hat{C}_{k}$, being the integrals over $C_{k}(y_{1},...,y_{k})$, equal zero, which does not necessary mean $C_{k}(y_{1},...,y_{k})=0$.

Using Eq. (124) we obtain $G(z)=\left\langle N\right\rangle\left(z-1\right)$, and consequently $\hat{C}_{k}=0$ for $k\geq 1$. By the same token we have $K(t)=$ $\left\langle N\right\rangle\left(e^{t}-1\right)$ and $\kappa_{n}=\left\langle N\right\rangle$. The last relation can be immediately seen from Eq. (143) in the absence of $\hat{C}_{k}$. Clearly, for the Poisson distribution all cumulant ratios equals one

[TABLE]

Before we discuss more complicated and interesting examples, we want to explain why physically the factorial cumulants $\hat{C}_{k}$, being the integrals of the genuine correlation functions, are zero for the Poisson distribution. Suppose that we generate some number of particles, $N$, from a given multiplicity distribution $P(N)$, which we further assume to be the Poisson distribution. Next we randomly split the particles between the left and the right bin (for example in rapidity) with probability $1/2$, namely, for each particle we choose randomly whether it goes to the left or to the right bin. Now, we can calculate $P(N_{L},N_{R})$, being the joint probability to observe $N_{L}$ particles in the left interval and $N_{R}$ particles in the right interval. Obviously $\left\langle N_{L}\right\rangle=\left\langle N_{R}\right\rangle$. We obtain

[TABLE]

where the multiplicity distribution $P(N)$ is taken at $N=N_{L}+N_{R}$. We observe that in general $P(N_{L},N_{R})$ cannot be factorized to the product $P_{L}(N_{L})P_{R}(N_{R})$ where $P_{L}$ is the multiplicity distribution in the left bin and $P_{R}$ that for the right bin. The factorization works for the Poisson distribution, Eq. (145), since $(N_{L}+N_{R})!$ cancels and we obtain454545This result can be easily generalized to more than two bins and to arbitrary values of probabilities of particle distribution between the bins.

[TABLE]

which is the product of the two Poisson distributions. It means that for the Poisson distribution there is no correlation between the numbers of particles in the left and the right bins. Obviously any other distribution would introduce multiplicity correlations. For example, if a source is producing always (in each event), say, $10$ particles, and we know that, say, $N_{L}=2$ we immediately know that $N_{R}=8$. In this case $N_{L}$ and $N_{R}$ are maximally correlated (or more precisely anti-correlated).

To better understand the cumulants and the factorial cumulants let us consider a slightly more complicated model. Suppose that clusters (e.g., resonances) are distributed according to Poisson distribution $P_{c}(N_{c})$, where $N_{c}$ is the number of clusters. Next, each cluster decays into exactly $m$ final particles, namely, $N_{c}$ clusters result in $N=mN_{c}$ particles. In this case, the final multiplicity distribution is given by

[TABLE]

The generating function is given by

[TABLE]

where the last equality assumes $P_{c}(N_{c})$ to be the Poisson distribution with the average number of clusters $\left\langle N_{c}\right\rangle$. It follows directly from Eq. (146) where $z$ is replaced by $z^{m}$. Performing suitable calculations we obtain $F_{1}=m\left\langle N_{c}\right\rangle=\left\langle N\right\rangle$ and

[TABLE]

and more complicated expressions for the higher order factorial moments. Clearly $F_{i}>\left\langle N\right\rangle^{i}$ which indicates specific correlations.

On the other hand, the factorial cumulants, generated from $\ln(H)=\left\langle N_{c}\right\rangle\left(z^{m}-1\right)$, read $\hat{C}_{1}=m\left\langle N_{c}\right\rangle=\left\langle N\right\rangle$ and for $k\geq 2$ we have [347]

[TABLE]

In other words, the factorial cumulants vanish if the order of the factorial cumulant is larger than $m$ - the number of particle from a single cluster. This is not surprising. The factorial cumulants are directly related to correlations present in the system. Having clusters decaying exactly into $m$ particles we introduce two- and up to $m$-particle correlations but not $(m+1)$-particle correlations. Consequently $\hat{C}_{m+1}$, $\hat{C}_{m+2}$ etc. equal zero. It is worth mentioning that in our simple example clusters themselves do not bring any correlations since their multiplicity distribution is given by Poisson. However, any other distribution different than Poisson would in general result in $\hat{C}_{k}\neq 0$ also for $k>m$.464646The interpretation of factorial cumulants can at times be rather tricky. For example, one may have n-particle correlations even if one has less that $n$ particles in the system. This is best illustrated by the extreme example, where we have a source that in each event produces exactly one particle. In this case $P(N=1)=1$ and $H(z)=z$. Obviously $F_{1}=1$ and $F_{i}=0$ for $i\geq 2$. The same relations hold for the cumulants. However, the factorial cumulants are $\hat{C}_{i}=(-1)^{i-1}(i-1)!$. For example, the two particle density is zero (we have only one particle) and thus following Eq. (129), we obtain $C_{2}(y_{1},y_{2})=-\rho(y_{1})\rho(y_{2}),$ which integrates to $\hat{C}_{2}=-1$. Therefore, even though we have only one particle we have integrated n-particle correlations for all $n$. This seemingly paradoxical situation is best resolved by realizing that the factorial cumulants, $\hat{C}_{k}$, “measure” the deviation from Poisson distribution. In our example the distribution with exactly one particle deviates considerably from Poisson with $\left\langle N\right\rangle=1$. It is these deviations which are encoded in the $\hat{C}_{k}$. This example also nicely demonstrates that interactions between produced particles should not be confused with correlations since both are not necessarily related, as is evident by a source producing one particle.

As seen from Eq. (143), the cumulants $\kappa_{i}$ are influenced by the factorial cumulants up to the $i$-th order. In particular, if $m=2$, that is we have only two-particle correlations in the system, all cumulants are influenced, which makes the interpretation of the cumulants rather tricky. Performing straightforward calculations we obtain

[TABLE]

which follows directly from the cumulant generating function, which in our simple model reads $K(t)=\left\langle N_{c}\right\rangle\left(e^{mt}-1\right)$. For example, taking $m=2$ we obtain $\kappa_{4}/\kappa_{2}=4$ even though the system has only two particle correlations. This demonstrates that the higher order cumulants do not necessarily carry information about nontrivial multiparticle correlations, and that factorial cumulants seem to be better suited for this purpose.

Another example worth considering in more detail is the bi-modal or two-component distribution discussed in Sec. 7.1.5. In this case the multiplicity distribution is given by Eq. (117)

[TABLE]

so that the factorial moment generating function is given by

[TABLE]

Here $H_{(a,b)}(z)$ are the generating functions for the distributions $P_{(a,b)}(N)$. The factorial cumulant generating function is then

[TABLE]

where $G_{(a,b)}(z)=\ln[H_{(a,b)}(z)]$ are the factorial cumulant generating functions for the distributions $P_{(a,b)}(N)$. Noting, that the second term in the above equation resembles the cumulant generating function $K_{Bn}(t)$ of the Bernoulli distribution with the probability $p=\alpha$, which is given by

[TABLE]

we have

[TABLE]

Let us now assume that both, $P_{(a)}(N)$ and $P_{(b)}(N)$ are Poisson distributions so that all but the first factorial cumulants vanish,

[TABLE]

In this case the factorial cumulants for the bi-modal distribution $P(N)$ are simply given by

[TABLE]

Here, $\kappa_{Bn,k}$ denote the cumulants of the Bernoulli distribution with $\kappa_{Bn,1}=\alpha$, $\kappa_{Bn,2}=\alpha(1-\alpha)$, $\kappa_{Bn,3}=\alpha(1-\alpha)(1-2\alpha)$ and so on. For the bi-modal distribution discussed in Sec. 7.1.5, which fits the measured STAR cumulants, we have $\alpha\ll 1$ and $\left\langle N_{(a)}\right\rangle-\left\langle N_{(b)}\right\rangle\simeq 15$ so that, to a good approximation, the higher order factorial cumulants are given by

[TABLE]

Therefore, with increasing order the factorial cumulants increase in magnitude (by a factor of $\sim 15$ in case of the fit to the STAR data ) and alternate in sign. Most interestingly, we obtain large factorial cumulants even though in each event particles are produced from a single Poisson distribution, which by itself has vanishing factorial cumulants. What happens in the case of the two-component model can be easily understood following the discussion around Eq. (150). Suppose we split particles originating from $P(N)$ between the left and the right bin. Having a large number of particles in the left bin it is more likely that they are originating from $P_{(a)}(N)$, since in our example $\left\langle N_{(a)}\right\rangle>\left\langle N_{(b)}\right\rangle$. Consequently, we expect a large number of particles in the right bin as well. Similarly, heaving a small number of particles in the left bin indicates that $P_{(b)}(N)$ is active and we expect a smaller number of particles in the right bin. In other words, having a superposition of two Poisson distributions with different means results in $P(N_{L},N_{R})$ being different than $P(N_{L})P(N_{R})$, which indicates multiplicity correlations.

The final example we want to discuss is related to the multiplicity distribution of charged particles, which is typically described with the negative binomial distribution (see, e.g., Ref. [397])

[TABLE]

where $k$ measures deviation from the Poisson distribution.474747The negative binomial distribution goes to the Poisson distribution if $k\rightarrow\infty$ at fixed $\left\langle N\right\rangle$. For example $\left\langle N^{2}\right\rangle-\left\langle N\right\rangle^{2}=\left\langle N\right\rangle\left(1+\frac{\left\langle N\right\rangle}{k}\right)$ which goes to $\left\langle N^{2}\right\rangle-\left\langle N\right\rangle^{2}=\left\langle N\right\rangle$. In this case, the generating function reads

[TABLE]

The first factorial moment is obviously $F_{1}=\left\langle N\right\rangle$ and for $i\geq 2$ we have

[TABLE]

which goes to the Poisson limit, $F_{i}=\left\langle N\right\rangle^{i}$, if $k\rightarrow\infty$ at fixed $\left\langle N\right\rangle$. Taking the logarithm of Eq. (166) and using Eq. (125) we obtain the factorial cumulant

[TABLE]

which goes to zero if $k\rightarrow\infty$ at fixed $\left\langle N\right\rangle$, being the proper Poisson limit. Finally, the cumulants are given by a sum of the factorial cumulants as shown in Eq. (143). The negative binomial distribution is usually defined with $k>0$ and thus the factorial cumulants are all positive and consequently the cumulants are larger than $\left\langle N\right\rangle$, and for example

[TABLE]

which is larger then one (unless $k\rightarrow\infty$).