Fluctuations of conserved charges in relativistic heavy ion collisions: An introduction

# Fluctuations of conserved charges in relativistic heavy ion collisions: An introduction

Masayuki Asakawa and Masakiyo Kitazawa
Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan
###### Abstract

Bulk fluctuations of conserved charges measured by event-by-event analysis in relativistic heavy ion collisions are observables which are believed to carry significant amount of information on the hot medium created by the collisions. Active studies have been done recently experimentally, theoretically, and on the lattice. In particular, non-Gaussianity of the fluctuations has acquired much attention recently. In this review, we give a pedagogical introduction to these issues, and survey recent developments in this field of research. Starting from the definition of cumulants, basic concepts in fluctuation physics, such as thermal fluctuations in statistical mechanics and time evolution of fluctuations in diffusive systems, are described. Phenomena which are expected to occur in finite temperature and/or density QCD matter and their measurement by event-by-event analyses are also elucidated.

## 1 Introduction

### 1.1 Background

The medium described by quantum chromodynamics (QCD) is expected to have various phase transitions with variations of external thermodynamic parameters such as temperature . Although the basic degrees of freedom of QCD, quarks and gluons, are confined into hadrons in the vacuum, they are expected to be liberated at extremely high temperature and form a new state of the matter called the quark-gluon plasma (QGP). It is also known that the chiral symmetry, which is spontaneously broken in vacuum, is restored at extremely hot and/or dense environment. These phase transitions at vanishing baryon chemical potential () are investigated with lattice QCD Monte Carlo simulations. The numerical analyses show that the phase transition is a smooth crossover [1, 2]. On the other hand, various models predict that there exists a discontinuous first order phase transition at nonzero . The existence of the endpoint of the first order transition, the QCD critical point [3], and possibly multiple critical points [4], are anticipated in the phase diagram of QCD on - plane [5, 6, 7].

After the advent of the relativistic heavy ion collisions, the quark-gluon plasma has come to be created and investigated on the Earth. At the Relativistic Heavy Ion Collider (RHIC) [8] and the Large Hadron Collider (LHC) [9], active experimental studies on the QGP have been being performed. The discovery of the strongly-coupled property of the QGP near the crossover region [8, 9] is one of the highlights of these experiments. With the top RHIC energy GeV and LHC energy TeV, hot medium with almost vanishing is created [10, 11]. On the other hand, the chemical freezeout picture for particle abundances [12] suggests that the net-baryon number density and of the hot medium increase as is lowered down to GeV [11, 13]. The relativistic heavy ion collisions, therefore, can investigate various regions of the QCD phase diagram on - plane by changing the collision energy . Such an experimental program is now ongoing at RHIC, which is called the beam-energy scan (BES) program [14]. The upgraded stage of the BES called the BES-II is planned to start in 2019 [15]. The future experiments prepared at FAIR [16], NICA [17] and J-PARC will also contribute to the study of the medium with large . The searches of the QCD critical point [18] and the first-order phase transition are among the most interesting subjects in this program.

In relativistic heavy ion collisions, after the formation of the QGP the medium undergoes confinement transition before they arrive at the detector. During rescatterings in the hadronic stage, the signals formed in the QGP tend to be blurred. In order to study the properties of QGP in these experiments, therefore, it is important to choose observables which are sensitive to the medium property in the early stage.

Recently, as unique hadronic observables which reflect thermal property of the primordial medium created by relativistic heavy ion collisions, the bulk fluctuations have acquired much attentions [19]. Although these observables are hadronic ones, it is believed that they reflect the thermal property in the early stage [20]. They are believed to be good observables in investigating the deconfinement transition [21, 22, 23] and finding the location of the QCD critical point [18, 24, 25]. Active experimental studies have been carried out [26, 27, 28, 29, 30, 31, 32, 33] as well as analyses on the numerical simulations on the lattice [34, 35]. In particular, fluctuations of conserved charges and their higher order cumulants representing non-Gaussianity [23, 36, 37] are actively studied recently. The purpose of this review is to give a basic introduction to the physics of fluctuations in relativistic heavy ion collisions, and give an overview of the recent experimental and theoretical progress in this field of research.

### 1.2 Fluctuations

Before starting the discussion of relativistic heavy ion collisions, we first give a general review on fluctuations briefly. When one measures an observable in some system, the result of the measurements would take different values for different measurements, even if the measurement is performed with an ideal detector with an infinitesimal resolution. This distribution of the result of measurements is referred to as fluctuations. In typical thermal systems, the fluctuations are predominantly attributed to thermal effects, which can be calculated in statistical mechanics. Quantum effects also give rise to fluctuations.

In contrast to standard observables, fluctuations are sometimes regarded as the noise associated with the measurement and thus are obstacles. As expressed by Landauer as “The noise is the signal” [38], however, the fluctuations sometimes can become invaluable physical observables in spite of their obstacle characters. Here, in order to spur the motivations of the readers we list three examples of the physics in which fluctuations play a crucial role.

1. Brownian motion: The first example is a historical one on Brownian motion. As first discovered by Brown in 1827, small objects, such as pollens, floating on water show a quick and random motion. Due to this motion, the position of the pollen after several time duration fluctuates even if the initial position is fixed. The origin of this motion was first revealed by Einstein. In his historical paper in 1905 [39], Einstein pointed out that the Brownian motion is attributed to the thermal motion of water molecules. This prediction was confirmed by Perrin, who calculated the Avogadro constant based on this picture [40]. In this era, the existence of atoms had not been confirmed, yet. These studies served as a piece of the earliest evidence for the existence of molecules and atoms. In other words, human beings saw atoms for the first time behind fluctuations.

This example tells us that fluctuations are powerful tools to diagnose microscopic physics. One century after Einstein’s era, now that we know substructures of atoms, hadrons, and quarks and gluons, it seems a natural idea to utilize fluctuations in relativistic heavy ion collisions in exploring subnuclear physics. In this review, a Brownian particle model for diffusion of fluctuations will be discussed in Sec. 5.4.

2. Cosmic microwave background: The second example is found in cosmology. As a remnant of Big Bang and as a result of transparent to radiation, our Universe has  K thermal radiation called cosmic microwave background (CMB) [41]. The temperature of this radiation is almost uniform in all directions in the Universe, but has a tiny fluctuation at different angles. This fluctuation is now considered as the remnant of quantum fluctuations in the primordial Universe. With this picture the power spectrum of this fluctuation tells us various properties of our Universe [42]; for example, our Universe has started with an inflational expansion billion years ago. In other words, we can see the hot primordial Universe behind the fluctuation of CMB.

This example tells us that fluctuations can be powerful tools to trace back the history of a system. It thus seems a natural idea to utilize fluctuations to investigate the early stage of the “little bang” created by relativistic heavy ion collisions. A common feature in the study of fluctuations in CMB and heavy ion collisions is that the non-Gaussian fluctuations acquire attentions. In fact, the non-Gaussianity of the CMB has been one of the hot topics in this community [43, 44], although the Planck spacecraft has not succeeded in the measurement of statistically significant non-Gaussianity thus far [45].

3. Shot noise: The final example is the fluctuations of the electric current in an electric circuit called the shot noise. The electric current at a resistor is generally fluctuating. Even without an applied voltage, the current has thermal noise proportional to , which is called the Johnson-Nyquist noise [46, 47]. On the other hand, there is a contribution of the noise which takes place when a voltage is applied and is proportional to the average current . When a circuit has a potential barrier, variance of such a noise tends to be proportional to , where is the electric charge of the elementary degrees of freedom carrying electric current. (This proportionality comes from the Poisson nature of the noise as will be clarified in Sec. 3.2.3.) This noise is called the shot noise [48]. Because of this proportionality, this fluctuation can be used to investigate the quasi-particle property. When the material undergoes the phase transition to superconductivity, for example, electrons are “confined” into Cooper pairs and the electric charge carried by the elementary degrees of freedom is doubled. This behavior is in fact observed in the measurement of the shot noise [49]. More surprisingly, in the materials in which the fractional quantum Hall effect is realized, the shot noise behaves as if there were excitations having fractional charges [50].

This example tells us that the fluctuations are powerful tools to investigate elementary degrees of freedom in the system although they are macroscopic observables. It thus seems a natural idea to utilize this property of fluctuations to explore the confinement/deconfinement property of quarks in relativistic heavy ion collisions. In fact, this is a relatively old idea in heavy ion community [21, 22, 23]. It is also notable that the non-Gaussianity of the shot noise has been observed in mesoscopic systems [51].

### 1.3 Bulk fluctuations in relativistic heavy ion collisions

In this review, among various fluctuations we focus on the bulk fluctuations of conserved charges. When one measures a charge in a phase space in some system, the amount of the charge, , fluctuates measurement by measurement. We refer to the distribution of as the bulk fluctuation (or, simply fluctuation). When we perform this measurement in a spatial volume in a thermal system, this fluctuation is called the thermal fluctuation.

The bulk fluctuations are closely related to correlation functions. The total charge in a phase space is given by the integral of the density of the charge as

 Q=∫Vdxn(x), (1)

where is the coordinate in the phase space. The variance of thus is given by

 ⟨δQ2⟩V=⟨(Q−⟨Q⟩V)2⟩V=∫Vdx1dx2⟨δn(x1)δn(x2)⟩, (2)

where . In this equation, the left-hand side is the quantity that we call (second-order) bulk fluctuation, while the integrand on the right-hand side is called correlation function. Equaiton (2) shows that can be obtained from the correlation function by taking the integral. If one knows the value of for all the correlation function can also be constructed from . In this sense, the correlation function carries the same physical information as the bulk fluctuation, and the choice of observables, bulk fluctuation or correlation function, is a matter of taste for the second-order fluctuation . (For higher orders, correlation functions contain more information than bulk fluctuations.) Phenomenological studies on the correlation functions of conserved charges in relativistic heavy ion collisions are widely performed, especially in terms of the so-called balance function [52]. The relation between the correlation functions and bulk fluctuations are also discussed in the literature [53, 54, 55]. In this review, however, we basically stick to bulk fluctuations in our discussion.

In relativistic heavy ion collisions, the bulk fluctuations are observed by the event-by-event analyses. In these analyses, the number of some charge or a species of particle observed by the detector is counted event by event. The distribution of the numbers counted in this way is called event-by-event fluctuation. As we will discuss in detail in Sec. 4, the fluctuations observed in this way are believed to carry information on the thermal fluctuation in the primordial stage. The event-by-event fluctuations, however, are not the thermal fluctuations themselves. The hot media created by collisions are dynamical systems, and the detector can only measure their final states. Moreover, fluctuations other than thermal fluctuations contribute to event-by-event fluctuations. Careful treatment and interpretation, therefore, are required in comparing the event-by-event fluctuations with theoretical analysis on thermal fluctuations. As will be discussed in Sec. 4, however, there are sufficient reasons to expect that event-by-event fluctuations can be compared with thermal fluctuations with an appropriate treatment.

In Figs. 1 and 2, we show two examples of recent experimental results on event-by-event analyses of bulk fluctuations. Figure 1 shows the experimental result obtained by ALICE Collaboration at LHC [27]. This figure shows the variance of the net-electric charge; the right vertical axis shows the quantity called D-measure, i.e. the variance normalized in such a way that the value in the equilibrated hadronic medium becomes [22]. The horizontal axis is the rapidity window to count the particle number. The figure shows that the experimental result has a nontrivial suppression, which cannot be described by the equilibrated hadronic degrees of freedom. The result thus suggests that the net-electric charge fluctuation at the LHC energy contains non-hadronic or non-thermal physics in the primordial medium. The origin of the suppression in Fig. 1 will be discussed in Sec. 3.3. The experimental result also shows that the fluctuation is more suppressed for larger . This behavior will be discussed in detail in Sec. 5.2.

In Fig. 2, we show the experimental results on the non-Gaussian fluctuations measured by STAR Collaboration at RHIC [28, 32]. The two panels show the same quantity, the ratios of the net-proton number cumulants, as a function of the collision energy ; since the baryon chemical potential of the hot medium becomes smaller as increases [11, 13], these plots can be interpreted as the dependence of the cumulants. The right panel [32] is the updated version of the left panel [28]; the quality of the experimental analysis is ever-improving in this field [32]. In the right panel, the vertical axes are quantities which is expected to take unity in the equilibrated hadronic medium. The panel shows that these quantities take values which are close to the hadronic one but have statistically-significant deviation from those values. These deviations are believed to be important observables to explore the QCD phase structure. In this review, we will elucidate the meanings of the vertical axes in Fig. 2 and the reason why these quantities are widely discussed.

### 1.4 Contents of this review

In this review, we give a pedagogical introduction to the physics of fluctuations in relativistic heavy ion collisions. In particular, one of the objectives of the introductory part is to understand the meanings of Figs. 1 and 2. We aimed at answering, for example, the following questions in this review:

• What are the cumulants? Why should we focus on these quantities in the discussion of non-Gaussian fluctuations? Why are the cumulants sometimes called susceptibility, and what are the relation of cumulants with moments, skewness and kurtosis?

• Meanings of the vertical axes of Figs. 1 and 2. How to understand these experimental data?

• What happens in the fluctuation observables in heavy ion collisions if the hot medium undergoes a phase transition to deconfinement transition, or passes near the QCD critical point?

• Why can the event-by-event fluctuation be compared with thermal fluctuations? Why should one not directly compare the event-by-event fluctuations with thermal fluctuations?

• The concept of “equilibration of the fluctuation of conserved charges.” What is the difference of this concept from the “local equilibration,” and why should we distinguish them?

In addition to the answers to these questions, we have tried to describe recent progress in this field of research.

The outline of this review is as follows. In the next section, we give a pedagogical review on the probability distribution function, which is a basic quantity describing fluctuations. The cumulants are introduced here, and their properties, especially the extensive nature, are discussed. In Sec. 3, we discuss the thermal fluctuations, i.e. the fluctuations in an equilibrated medium. The behaviors of cumulants in the QCD phase diagram is also considered. Sections 46 are devoted to a review on the event-by-event fluctuations in experimental analyses. In Sec. 4, we summarize general properties of the event-by-event fluctuations. Various cautions in the interpretation of these quantities are given in this section. In Sec. 5, we focus on the non-equilibrium property of the event-by-event fluctuations, by describing the time evolution of fluctuations using stochastic formalisms. In Sec. 6, we consider a model for probability distribution functions which treats the efficiency problem in the observation of fluctuations. The difference between net-baryon and net-proton number cumulants are also discussed here. We then give a short summary in Sec. 7.

## 2 General introduction to probability distribution function

Because fluctuation is a distribution of some observable, it is mathematically represented by probability distribution functions. For example, if one repeats a measurement of an observable in an equilibrated medium many times, the result of the measurement would fluctuate measurement by measurement. The distribution of the result of the measurement is represented by a histogram. After accumulating the results of many measurements, the histogram with an appropriate normalization can be regarded as the probability distribution function. This distribution is nothing other than the fluctuation. In many contexts, the width of the distribution is particularly called fluctuations.

In this section, as preliminaries of later sections we give a pedagogical introduction to basic concepts in probability. We introduce moments and cumulants as quantities characterizing probability distribution functions. Advantages of cumulants, especially for the description of non-Gaussianity of distribution functions, are elucidated. We also discuss properties of specific distribution functions, Poisson, Skellam, binomial and Gauss distributions. The properties of these distribution functions play important roles in later sections for interpreting fluctuation observables in relativistic heavy ion collisions.

### 2.1 Moments and cumulants

We start from a probability distribution function satisfying for an integer stochastic variable . One of the set of quantities which characterizes is the moments. The -th order moment is defined by

 ⟨mn⟩=∑mmnP(m), (3)

where the bracket on the left-hand side represents the statistical average with . If the moments for all exists, they carry all information encoded in . For a probability distribution function for a continuous stochastic variable , the moments are defined by

 ⟨xn⟩=∫dxxnP(x), (4)

where the integral is taken over the range of .

To calculate the moments for a given probability distribution , it is convenient to introduce the moment generating function,

 G(θ)=∑memθP(m)=⟨emθ⟩. (5)

Moments are then given by the derivatives of as

 ⟨mn⟩=dndθnG(θ)∣∣∣θ=0. (6)

For the continuous case the generating function is defined by

 G(θ)=∫dxexθP(x). (7)

For many practical purposes, it is more convenient to use cumulants rather than moments for characterizing a probability distribution. To define the cumulants, we start from the cumulant generating function,

 K(θ)=lnG(θ). (8)

The cumulants of are then defined by

 ⟨mn⟩c=dndθnK(θ)∣∣∣θ=0. (9)

As we will see below, cumulants have several useful features for describing fluctuations, especially their non-Gaussianity.

Before discussing the advantages of cumulants, let us clarify the relation between moments and cumulants. These relations are obtained straightforwardly from their definitions Eqs. (6) and (9). For example, to write cumulants in terms of moments, we calculate as follows:

 ⟨m⟩c =ddθlnG(θ)∣∣∣θ=0=G(1)(0)G(0)=⟨m⟩, (10) ⟨m2⟩c =d2dθ2lnG(θ)∣∣∣θ=0=G(2)(0)G(0)−(G(1)(0))2(G(0))2=⟨m2⟩−⟨m⟩2=⟨δm2⟩, (11)

where represents the -th derivative of and we have used . In the last equality we defined . By repeating a similar manipulation, one can extend the relation to an arbitrary order. The results for third- and fourth-orders are given by

 ⟨m3⟩c =⟨δm3⟩, (12) ⟨m4⟩c =⟨δm4⟩−3⟨δm2⟩2. (13)

Note that the first-order cumulant is equal to the first-order moment, or the expectation value. The second- and third-order cumulants are given by the central moments,

 ⟨δmn⟩=⟨(m−⟨m⟩)n⟩. (14)

In particular, the second-order cumulant corresponds to the variance. This quantity is sometimes called simply fluctuation, because for many purposes the cumulants higher than the second-order are not physically significant. The cumulants for are given by nontrivial combinations of central moments with -th and lower orders.

All cumulants except for the first-order one are represented by central moments and do not depend on the average . To prove this statement, we consider a probability distribution function in which the distribution is shifted by compared with . The cumulant generating function of is calculated to be

 K′(θ) =ln∑memθP′(m)=ln∑memθP(m−m0)=ln∑me(m+m0)θP(m) =ln∑memθP(m)+m0θ=K(θ)+m0θ, (15)

where is the cumulant generating function of . Equation (15) shows that the difference between and is a term . The derivatives of and higher than the first-order thus are equivalent. Therefore, the cumulants higher than the first-order do not depend on , and they are represented only by central moments.

The expressions of moments in terms of cumulants are similarly obtained as follows:

 ⟨m⟩ =ddθeK(θ)∣∣∣θ=0=K(1)(0)eK(0)=⟨m⟩c, (16) ⟨m2⟩ =d2dθ2eK(θ)∣∣∣θ=0=(K(2)(0)+(K(1)(0))2)eK(0)=⟨m2⟩c+⟨m⟩2c, (17)

and so forth. Up to the fourth-order, one obtains

 ⟨m3⟩ =⟨m3⟩c+3⟨m2⟩c⟨m⟩c+⟨m⟩3c, ⟨m4⟩ =⟨m4⟩c+4⟨m3⟩c⟨m⟩c+3⟨m2⟩2c+6⟨m2⟩c⟨m⟩2c+⟨m⟩4c. (18)

### 2.2 Sum of two stochastic variables

An important property of cumulants becomes apparent when one considers the sum of two stochastic variables. Let us consider two integer stochastic variables and which respectively obey probability distribution functions and which are not correlated. Then, the probability distribution of the sum of two stochastic variables, , is given by

 P(m)=∑m1,m2δm,m1+m2P(m1)P(m2). (19)

(To understand Eq. (19), one may, for example, imagine the probability distribution of the sum of the numbers of two dices.) The moment and cumulant generating functions for are calculated to be

 G(θ) =∑memθP(m)=∑memθ∑m1,m2δm,m1+m2P1(m1)P2(m2) =∑m1em1θP1(m1)∑m2em2θP2(m2)=G1(θ)G2(θ), (20) K(θ) =lnG(θ)=K1(θ)+K2(θ), (21)

where and are the moment and cumulant generating functions of , respectively, for and . By taking derivatives of the both sides of Eq. (21), one finds

 ⟨mn⟩c=⟨mn1⟩c+⟨mn2⟩c. (22)

This result shows that the cumulants of the probability distribution for the sum of two independent stochastic variables are simply given by the sum of the cumulants. (This is the reason why the cumulants are called in this way.) Note that this result is obtained for two independent stochastic variables; when the distributions of and are correlated, Eq. (22) no longer holds.

### 2.3 Cumulants in statistical mechanics

In statistical mechanics, results of measurement of observables in a volume are fluctuating, and one can define their cumulants from the distribution of the results. From Eq. (22) one can argue an important property of the cumulants in statistical mechanics that the cumulants of extensive variables in grand canonical ensemble are extensive variables.

To see this, let us consider the number of a conserved charge in a volume in grand canonical ensemble. From the distribution of the result of measurements, one can define the cumulants of the charge. Next, let us consider the cumulants of the particle number in a twice larger volume, . This system can be separated into two subsystems with an equal volume . In statistical mechanics, it is usually assumed that the subsystems are uncorrelated when the volume is sufficiently large, and the property of the system does not depend on the shape of . Therefore, the particle number in the total system is regarded as the sum of the two independent particle numbers in the two subsystems. From Eq. (22), thus are represented as

 ⟨Nn⟩c,2V=2⟨Nn⟩c,V. (23)

By similar arguments one obtains,

 ⟨Nn⟩c,λV=λ⟨Nn⟩c,V (24)

for an arbitrary number . Equation (24) shows that the cumulants of in statistical mechanics are extensive variables. As special cases of this property, the average particle number and the variance in statistical mechanics are extensive variables.

From Eq. (24), the cumulants in volume can be written as

 ⟨Nn⟩c,V=χnV. (25)

Here, is the density of the particle, and is the quantity which is referred to as susceptibility because of the linear response relation discussed in Sec. 3.1.3. We call for as generalized susceptibilities.

Remarks on the extensive nature of the cumulants are in order. First, the above argument is valid only when the volume of the system is sufficiently large. When the spatial extent of the volume is not large enough, the correlation between two adjacent volumes becomes non-negligible. This would happen when the spatial extent of the volume is comparable with the microscopic correlation lengths. Next, in the above argument we have implicitly assumed grand canonical ensemble in addition to the equilibration. The argument, for example, is not applicable to subvolumes in canonical ensemble in which the number of in the total system is fixed. In this case, the fixed total number gives rise to correlation between subvolumes; because the total number is fixed, if the particle number in a subvolume is large the particle number in the other subvolume tends to be suppressed. This correlation violates the assumption of independence between the particle numbers in subvolumes unless the subvolume is small enough compared with the total volume.

### 2.4 Examples of distribution functions

Now, we see some specific distribution functions, which play important roles in relativistic heavy ion collisions.

#### 2.4.1 Binomial distribution

The binomial distribution function is defined by the number of “successes” of independent trials, each of which yields a success with probability . The binomial distribution function is given by

 Bp,N(m)=NCmpm(1−p)N−m, (26)

where

 NCm=N!m!(N−m)! (27)

is the binomial coefficient. It is easy to show that using the binomial theorem. The moment and cumulant generating functions of the binomial distribution are calculated to be

 GB(θ) =∑memθBp,N(m)=∑mNCm(eθp)m(1−p)N−m =(1−p+eθp)N, (28) KB(θ) =Nln(1−p+eθp). (29)

From Eq. (29), the cumulants of the binomial distribution function are given by

 ⟨mn⟩c=ξnN (30)

with explicit forms of up to the fourth-order

 ξ1=p,ξ2=p(1−p),ξ3=p(1−p)(1−2p),ξ4=p(1−p)(1−6p+6p2). (31)

Equation (30) shows that the cumulants of the binomial distribution are proportional to , which is a reasonable result from the extensive nature of cumulants.

The sum of two stochastic variables obeying independent binomial distribution functions with an equal probability is again distributed with a binomial one. This can be shown by explicitly deriving the identity,

 Bp,N1+N2(m)=∑m1,m2δm,m1+m2Bp,N1(m1)Bp,N2(m2). (32)

An alternative way to prove this statement is to use Eqs. (22) and (30). Suppose that and obey the binomial distribution and , respectively. From Eq. (22), one finds that the cumulants of the sum are given by , which are nothing but the cumulants of the binomial distribution . Because all cumulants are those of the binomial distribution, the distribution of is given by the binomial one.

#### 2.4.2 Poisson distribution

The Poisson distribution function is defined by the number of successes of independent trials, each of which yields a success with infinitesimal probability. The Poisson distribution is thus obtained by taking the limit of the binomial distribution function with fixed . Replacing the binomial coefficients in Eq. (26) as , which is valid in this limit, and using the definition of Napier’s number in Eq. (26), one obtains

 Pλ(m)=λmm!e−λ. (33)

The cumulant generating function of Eq. (33) is obtained as

 Kλ(θ)=λ(eθ−1). (34)

By taking derivatives of Eq. (34), one finds that all the cumulants of the Poisson distribution are the same,

 ⟨mn⟩c=λ~{}~{}(for any n≥1). (35)

This property will be used frequently in later sections. This result also shows that the Poisson distribution is characterized by a single parameter , while the binomial distribution function is characterized by the two parameters, and .

The sum of two stochastic variables obeying two independent Poissonian obeys Poissonian,

 Pλ1+λ2(m)=∑m1,m2δm,m1+m2Pλ1(m1)Pλ2(m2), (36)

as one can explicitly show easily. Similarly to the case of the binomial distribution, however, a much easier way to show this is to use Eqs. (22) and (35).

The Poisson distribution is one of the most fundamental distribution function, as it naturally appears in various contexts. One interesting example among them is the classical ideal gas. One can show that the distribution of the number of a classical free particles in a volume in grand canonical ensemble obeys the Poisson distribution. To show this, the following intuitive argument suffices. First, consider a canonical ensemble for a volume with a fixed particle number . The grand canonical ensemble is defined by a subvolume in this system in the limit . Next, consider a particle in this system. The probability that this particle exists in the subvolume is given by . Moreover, because we consider classical free particles, the probability that each particle is in the subvolume is uncorrelated. The probability distribution function of the particle number in thus is given by the binomial distribution with probability . Because the grand canonical ensemble is defined by limit with fixed , the distribution in this limit obeys Poissonian. In Sec. 3.2 we will explicitly calculate the cumulants in the classical ideal gas and show that they satisfy Eq. (35). When the effect of quantum statistics, Bose-Einstein or Fermi-Dirac, shows up, this argument does not hold owing to the quantum correlation even for ideal gas.

#### 2.4.3 Skellam distribution

Although the sum of stochastic variables obeying independent Poisson distributions again obeys Poissonian, the difference of two stochastic variables obeying independent Poisson distributions is not given by a Poisson distribution. This is obvious from the fact that the difference can take a negative value, while the Poisson distribution take nonzero values only for positive . The difference,

 Sλ1,λ2(m)=∑m1,m2δm,m1−m2Pλ1(m1)Pλ2(m2), (37)

is called the Skellam distribution. The generating functions of the Skellam distribution are calculated to be

 G(θ) =∑memθ∑m1,m2δm,m1−m2Pλ1(m1)Pλ2(m2) (38) =∑m1em1θPλ1(m1)∑m2e−m2θPλ2(m2) (39) =Gλ1(θ)Gλ2(−θ), (40) K(θ) =Kλ1(θ)+Kλ2(−θ), (41)

where generating functions and are those of Poisson distribution. By taking derivatives of this result with Eq. (34), one obtains

 ⟨mn⟩c=⟨mn1⟩c+(−1)n⟨mn2⟩c=λ1+(−1)nλ2 (42)

for the Skellam distribution. This result shows that all even cumulants take a common value , while the odd cumulants take . This result also shows that the Skellam distribution is characterized by two parameters, and , or alternatively odd and even cumulants.

The explicit analytic form of is given by

 Sλ1,λ2(m)=e−(λ1+λ2)(λ1λ2)Im(2√λ1λ2), (43)

where is the modified Bessel function of the first kind.

The Skellam distribution plays an important role in later sections, because they describe the distribution of net particle number, i.e. the difference of the numbers of particles and anti-particles, in the classical ideal gas. The fluctuation of net-baryon number in hadron resonance gas, for example, obeys the Skellam distribution.

#### 2.4.4 Gauss distribution

So far we have considered stochastic variables taking integer values. An example of a distribution with a continuous stochastic variable is the Gauss distribution, which is defined by

 PG(x)=1σ√2πexp[−(x−x0)22σ2]. (44)

The normalization factor is required to satisfy . The generating functions are calculated to be

 G(θ) =∫dxeθx1σ√2πexp[−(x−x0)22σ2]=exp[x0θ+12σ2θ2], K(θ) =x0θ+12σ2θ2. (45)

We thus have

 ⟨x⟩ =x0,  ⟨x2⟩c=σ2, (46)

and

 ⟨xn⟩c =0for n≥3. (47)

The results in Eqs. (46) and (47) can, of course, also be obtained by explicitly calculating and , and so forth.

Equations (46) and (47) show that the cumulants higher than the second-order vanish for the Gauss distribution. In other words, nonzero higher order cumulants characterize deviations from the Gauss distribution function. This is the reason why the cumulants are used as quantities representing non-Gaussianity.

### 2.5 Variance, skewness and kurtosis

Till now, we have discussed cumulants as quantities characterizing distribution functions. When one wants to describe the deviation from the Gauss distribution, it is sometimes convenient to use the quantities called skewness and kurtosis [56]. These quantities are defined as

 S=⟨x3⟩c⟨x2⟩3/2c=⟨x3⟩cσ3,κ=⟨x4⟩c⟨x2⟩2c=⟨x4⟩cσ4, (48)

where is the variance defined by . The skewness and kurtosis are interpreted as the third- and fourth-order cumulants of the renormalized stochastic variable satisfying ,

 S=⟨~x3⟩c,κ=⟨~x4⟩c. (49)

In Fig. 3, we show typical distribution functions having nonzero skewness and kurtosis. All distribution functions shown in the figure satisfy and . As shown in the left panel, the skewness represents the asymmetry of the distribution function. The kurtosis, on the other hand, typically describes the “sharpness” of the distribution compared with the Gaussian one as in the right panel.

When the non-Gaussianity of a distribution function is discussed, the set of variables to be used, and , or the third and fourth-order cumulants, should be chosen depending on the problem. When the distribution is expected to obey some specific one of which the cumulants are known, the cumulants are much more convenient. For example, when the distribution is expected to obey the Poisson one and the difference from this distribution is concerned, one may focus on the cumulants normalized by the average, , which become unity in the Poisson distribution. The deviation from the Poisson distribution is then characterized by the difference of the ratio from unity. As we will see in Sec. 3, for example, the fluctuation of net-baryon number in the equilibrated hadronic medium is well described by the Skellam distribution. If the fluctuation carries some physics which cannot be described by hadronic degrees of freedom in equilibrium, it may deviate from the Skellam one. In order to find this deviation, the best observable to be used is the ratios of cumulants between even or odd orders [23], which become unity for the Skellam distribution as discussed in Sec. 2.4.3. Another example is the distribution of the topological charge, , in QCD. A model called the dilute instanton gas model for the topological sector in QCD predicts that the distribution of the topological charge is given by the Skellam one with [57]. The proximity of the ratios of the even order cumulants to unity thus is a useful measure to judge the validity of this model. Recently, the measurements of the topological cumulants are actively performed on the lattice [58, 59, 60].

Contrary to the ratios of cumulants, the magnitudes of skewness and kurtosis in the Poisson distribution depend on the average as

 S=λ−1/2,κ=λ−1. (50)

These quantities become arbitrary small as becomes larger. This behavior is related to the central limit theorem in statistics, which states that the sum of independent events approaches a Gauss distribution as the number of events to be summed is increased. Because the Poisson distribution can be interpreted as the result of the sum of independent events, it approaches the Gauss distribution for large . The and , the quantities characterizing the deviation from the Gauss distribution, become arbitrarily small in this limit. Similar tendency is also expected in the large volume limit of fluctuation observables in statistical mechanics. Fluctuations in macroscopic systems are well described by Gauss distributions, and and become irrelevant in the large volume limit. On the other hand, the higher order cumulants take nonzero values even in this limit. There, however, are practical difficulties in the measurement of higher order cumulants in large systems. In addition to the difficulty in the exact measurement of observables in large systems, there is a fundamental problem that the statistical error of higher order cumulants grows as the volume increases, as discussed in the next subsection.

### 2.6 Error of the cumulants

Related to the above discussion on non-Gaussianity in large systems, we give some remarks on the statistical error associated with the measurement of higher order cumulants.

The statistical error of an observable is estimated as

 (Δ^O)2=⟨(^O−⟨^O⟩)2⟩Nstat−1=⟨δ^O2⟩Nstat−1, (51)

where is the number of statistics, i.e. the number of the measurements of , and the expectation value is taken over this statistical ensemble.

Let us consider an extensive observable in statistical mechanics, such as particle number. As discussed already, cumulants of are extensive observables proportional to the volume . Now, we start from the first-order cumulant. Usually, the density of this quantity is concerned rather than . By substituting into in Eq. (51), one obtains

 (Δρ)2=⟨δN2⟩V2Nstat=χ2VNstat, (52)

where in the last equality we used the extensive nature of the second-order cumulant Eq. (25). We have also assumed that and used an approximation in the denominator. The result Eq. (52) shows that the error of is proportional to , which is a well-known result. With a fixed , the error of becomes smaller as becomes larger.

Next, we consider the statistical error of the second-order cumulant of . By substituting in , the error is calculated to be

 (Δ(δN2))2 =⟨((δN2)−⟨δN2⟩)2⟩Nstat=⟨δN4⟩−⟨δN2⟩2Nstat=⟨N4⟩c+2⟨N2⟩2cNstat =χ4V+2(χ2V)2Nstat. (53)

In this derivation, we have used the central moments as basic quantities rather than moments. This choice suppresses the correlation between first- and second-order moments. The central moments are converted to cumulants in the third equality. From Eq. (53), the error of the susceptibility is obtained by dividing both sides by as

 (Δχ2)2=χ4V−1+2χ22Nstat=2χ22Nstat+O(V−1). (54)

This result shows that the error of does not have dependence for sufficiently large . Contrary to the first-order case, the increase of does not reduce the statistical error of .

Similarly, the error of the higher order cumulants can be estimated by substituting the definition of the cumulants in Eq. (51). The higher order cumulants are represented by central moments in general as shown in sec. 2.1. (To obtain the explicit forms, one may substitute in Eqs. (10)–(13) and so forth.) It is therefore instructive to first see the error of the central moments. The -th order central moment is represented by the sum of the product of cumulants

 ⟨mn1⟩c⟨mn2⟩c⋯⟨mni⟩c, (55)

with and . Moreover, one can easily show that the coefficients of all terms in this decomposition are positive and nonvanishing. Now, because the cumulants are proportional to , the term which is leading in in large volume is the one containing the product of the largest number of the cumulants. For even , it is and one obtains the behavior of the central moments in large as

 ⟨(δN)n⟩∼(⟨N2⟩c)n/2(1+O(V−1))∼(χ2V)n/2+O(Vn/2−1). (56)

For odd , one obtains .

Using these dependences on , the statistical error of the central moments is given by

 (Δ(δNn))2=⟨δN2n⟩−⟨δNn⟩2Nstat∼(χ2V)nNstat, (57)

where in the last step we have used Eq. (56) and the fact that the terms proportional to never cancel out between and . Subleading terms in are neglected on the far right-hand side. The error of the -th order central moment thus is proportional to for large and .

Now, we come back to the statistical error of higher order cumulants. As discussed already, the cumulants can be represented by the central moments. Substituting the cumulants in into Eq. (51) and representing them in terms of central moments, the error of the cumulants is given by the sum of the product of the central moments. It is then concluded that the error of the cumulants in the leading order in should be same as Eq. (57) unless the highest order terms in cancel out. The error of the generalized susceptibility thus is expected to behave as

 Δχn∼√χn2Vn−2Nstat, (58)

for large and . Eq. (58) shows that the error of for grows as becomes larger and the measurement of non-Gaussian cumulants becomes more and more difficult as the spatial volume becomes larger. This dependence is highly contrasted to the error of standard observables Eq. (52), which becomes smaller as becomes larger. The proportionality coefficients in Eq. (58) up to the fourth-order are presented in Ref. [61].

In relativistic heavy ion collisions, non-Gaussian cumulants have been observed up to the fourth-order. These measurements are possible because the size of the system observed by detectors is not large; the particle number in each event is at most of order . The growth of the statistical error of higher order cumulants is a well known feature, and discussed in various contexts; see for example Refs. [62, 63, 64].

### 2.7 Cumulants for multiple variables

Next, we discuss probability distribution functions for multiple stochastic variables and their moments and cumulants.

Let us consider a probability distribution function for integer stochastic variables and . The moments of this distribution are defined by

 ⟨mn11mn22⟩=∑m1,m2mn11mn22P(m1,m2). (59)

By defining the moment generating function as

 G(θ1,θ2)=∑m1,m2eθ1m1eθ2m2P(m1,m2)=⟨eθ1m1eθ2m2⟩, (60)

the moments are given by

 ⟨mn11mn22⟩=∂n1∂θn11∂n2∂θn22G(θ1,θ2)∣∣ ∣∣θ1=θ2=0 (61)

similarly to the case of a single variable. The cumulants for are similarly defined with the cumulant generating function as

 ⟨mn11mn22⟩c=∂n1∂θn21∂n2∂θn22K(θ1,θ2)∣∣ ∣∣θ1=θ2=0. (62)

It is easy to extend the argument in Sec. 2.1 to relate the moments and cumulants to this case. The relation of the cumulants with moments for or is equal to the previous case. For the mixed cumulants, it is, for example, calculated to be

 ⟨m1m2⟩c =∂∂θ1∂∂θ2lnG(θ1,θ2)∣∣∣θ1=θ2=0=∂∂θ1(∂G(θ1,θ2)∂θ2G(θ1,θ2)−1)∣∣∣θ1=θ2=0 =(G(θ1,θ2)−1∂2G(θ1,θ2)∂θ1∂θ2−G(θ1,θ2)−2∂G(θ1,θ2)∂θ1∂G(θ1,θ2)∂θ2)∣∣∣θ1=θ2=0 =⟨δm1δm2⟩ (63)

and so forth. Equation (63) shows that the mixed second-order cumulant is given by the mixed central moment, or correlation.

For a probability distribution function for stochastic variables , the generating functions are defined by

 G(θ1,⋯,θk) =∑m1,m2,⋯,mk(k∏i=1eθimi)P(m1,⋯,mk), (64) K(θ1,⋯,θk) =lnG(θ1,⋯,θk). (65)

With these generating functions, the moments and cumulants are defined similarly. From these definitions, it is obvious that the cumulants for multi-variable distribution functions are extensive variables in grand canonical ensembles.

#### 2.8.1 Cumulant expansion

From the definition of cumulants Eq. (9), the cumulant generating function is expanded as

 K(θ)=∞∑n=1θnn!⟨mn⟩c. (66)

By substituting in Eq. (8), one obtains

 K(1)=ln⟨em⟩=ln∑memP(m)=∞∑n=1⟨mn⟩cn!. (67)

Equation (67) is called the cumulant expansion, and plays effective roles in obtaining various properties of higher order cumulants; examples are found in Appendix A and Ref. [65].

A remark here is that the cumulant expansion Eq. (67) has the same structure as the linked cluster theorem in field theory [66]. In fact, if one regards the “connected part” of correlation functions as the cumulant, the theorem is completely equivalent with the cumulant expansion.

#### 2.8.2 Factorial moments and factorial cumulants

Up to now we have discussed moments and cumulants as quantities characterizing probability distribution functions. One of other sets of such quantities is factorial moments and factorial cumulants. The factorial moments are defined as

 ⟨mn⟩f=⟨m(m−1)⋯(m−n+1)⟩=dndsnGf(s)∣∣∣s=1 (68)

with the factorial moment generating function

 Gf(s)=∑msmP(m)=G(lns). (69)

The factorial cumulants are then defined by the factorial cumulant generating function,

 Kf(s)=lnGf(s)=K(lns), (70)

as

 ⟨mn⟩fc=dndsnKf(s). (71)

In the discussion of physical property of fluctuations, the standard moments and cumulants tend to be more useful than the factorial moments and cumulants. For example, as we will see in the next section the cumulants of conserved charges are directly related to the partition function and have more apparent physical meanings owing to the linear response relation. For some analytic procedure and theoretical purposes, however, the factorial moments and cumulants make analyses more concise; see, for example, Refs. [67, 68, 69, 70].

To relate the moments and cumulants with factorial ones, one may use relations between the two generating functions, such as or : By taking derivatives with respect to or , their relations are obtained similarly to the procedure to obtain the relations between moments and cumulants presented in Sec. 2.1.

## 3 Bulk fluctuations in equilibrium

In this section we consider the properties of thermal fluctuations, i.e. the fluctuations in equilibrated medium. After describing a general property of fluctuations in quantum statistical mechanics, we calculate the fluctuations in ideal gas. This simple analysis allows us to understand many interesting properties of thermal fluctuations, such as the magnitude of cumulants in the hadronic medium [23]. We also review the linear response relations, which connect cumulants of conserved charges with the corresponding susceptibilities, i.e. magnitude of the response of the system against external force. The linear response relations for higher order cumulants enable us to introduce physical interpretation to the behavior of cumulants in the QCD phase diagrams [37]. Thermal fluctuations in the hadron resonance gas model and anomalous behavior of fluctuations expected to take place near the QCD critical point are also reviewed. We also give a brief review on the recent progress in the study of fluctuations in lattice QCD numerical simulations.

The anomalous behaviors of fluctuation observables in an equilibrated medium discussed in this section serves as motivations of experimental analyses of fluctuations. It, however, should be remembered that the fluctuations discussed in this section are not directly observed in relativistic heavy ion collisions, as will be discussed in detail in Secs. 4, 5 and 6.

### 3.1 Cumulants in quantum statistical mechanics

In this subsection, we first take a look at general properties of fluctuations, in particular, the cumulants of conserved charges, in quantum statistical mechanics.

#### 3.1.1 Cumulants of conserved charges

We consider a system described by a Hamiltonian in a volume and assume that this system has an observable which is a conserved charge. Because is conserved, commutes with ,

 [H,^N]=0. (72)

The grand canonical ensemble of this system with temperature and chemical potential for is characterized by the density matrix

 ρ=1Ze−(H−μ^N)/T, (73)

with the grand partition function

 Z=tr[e−(H−μ^N)/T], (74)

where the trace is taken over all quantum states. The expectation value of an observable is given by

 ⟨^O⟩=tr[^Oρ]. (75)

As in the previous section, one can define the moments and cumulants of in quantum statistical mechanics. The moments are simply given by the expectation values of powers of . The cumulants are then defined from the moments and the relations such as Eqs. (10) - (13), which relate moments and cumulants.

We note that the moments and cumulants defined in this way are interpreted as those for the distribution of the particle number in a volume . To be more specific, imagine that you were able to count the particle number in in the equilibrated medium exactly at some time. The obtained particle numbers then would fluctuate measurement by measurement around the average. The moments and cumulants are those of this fluctuation.

The cumulants of the conserved charge in quantum statistical mechanics are given by derivatives of the grand potential with respect to . The first-order cumulant, i.e. the expectation value, is given by

 ∂(−Ω/T)∂(μ/T)=1Ztr[^Ne−(H−μ^N)/T]=⟨N⟩. (76)

The second derivative gives the second-order cumulant as

 ∂2(−Ω/T)∂(μ/T)2 =∂∂(μ/T)(1Ztr[^Ne−(H−μ^N)/T])=1Ztr[^N2e−(H−μ^N)/T]−(1Ztr[^Ne−(H−μ^N)/T])2, =⟨δ^N2⟩=⟨^N2⟩c (77)

with . Similar manipulations lead to

 ⟨^Nn⟩c=∂n(−Ω/T)∂(μ/T)n. (78)

To show this relation, one may use the fact that is the moment generating function Eq. (5) up to normalization constant,

 ⟨Nn⟩=1Ztr[Nne−(H−μ^N)/T]=1Z∂nZ∂(μ/T)n. (79)

The normalization of Eq. (5) affects the definition of the cumulant generating function in Eq. (8) only by a constant. Since the constant does not alter derivatives, derivatives of the logarithm of give the cumulants.

As already discussed in Sec. 2.3, the cumulants of are extensive variables. This property can easily be shown using the fact that the grand potential is an extensive variable, and thus can be written using the grand potential per unit volume, , as

 Ω=ωV. (80)

The cumulants are thus given by

 ⟨^Nn⟩c=∂n(−ω/T)V∂(μ/T)n≡χnV, (81)

where on the far right-hand side we introduced the cumulant per unit volume,

 χn=∂n(−ω/T)∂(μ/T)n. (82)

The quantities defined here is called susceptibilities, as the reason will be explained in Sec. 3.1.3.

The extensive property of cumulants Eq. (81) gives a constraint on the correlation function of particle number density. The particle number in a volume is related to the particle density as

 ^N=∫Vdxn(x). (83)

The extensive property Eq. (81) then implies that

 ⟨n(x