Nonequilibrium time evolution of higher order cumulants of conserved charges and event-by-event analysis

We investigate the time evolution of higher order cumulants of conserved charges in a volume with the diffusion master equation. Applying the result to the diffusion of non-Gaussian fluctuations in the hadronic stage of relativistic heavy ion collisions, we show that the fourth-order cumulant of net-electric charge at LHC energy is suppressed compared with the recently observed second-order cumulant at ALICE, if the higher order cumulants at hadronization are suppressed compared with their values in the hadron phase in equilibrium. The significance of the experimental information on the rapidity window dependence of various cumulants in investigating the history of the dynamical evolution of the hot medium created in relativistic heavy ion collisions is emphasized.


INTRODUCTION
Statistical mechanics tells us that observables are fluctuating even in equilibrated medium. Because the fluctuations are determined by the microscopic nature of the medium and sensitive to critical phenomena, they can be exploited to reveal and characterize properties of the medium. In experimental attempts to map the global nature of QCD phase transition at nonzero baryon density in relativistic heavy ion collisions, fluctuation observables, especially those of conserved charges, are believed to be promising observables to diagnose the property of the hot medium [1][2][3][4][5][6][7][8]. Active investigation in heavy ion collisions by event-by-event analyses has recently been performed at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) [9][10][11]. Numerical analyses of higher order cumulants in equilibrium have been also carried out in lattice QCD Monte Carlo simulations [12].
As an experiment in which fluctuations are measured, heavy ion collisions have several notable features. First, higher order cumulants, which characterize the non-Gaussian nature of fluctuations, have been experimentally observed with good statistics up to the fourth order [9,10]. The measurement is possible because the system observed in the experiments is not large; the observed particle number is at most of order 10 3 , while the event-by-event statistics exceed 10 7 . Second, the fluctuations observed in experiments, especially those of conserved charges, are not necessarily the same as those in an equilibrated medium, because of the dynamical nature of the hot medium created by heavy ion collisions. Because of these properties, an appropriate description of the dynamical evolution of non-Gaussian fluctuations is required in order to understand the experimental results on higher order cumulants.
Concerning the second point, we remark that the recent experimental result on the net-electric charge fluc-tuation, (δN (net) Q ) 2 , by ALICE Collaboration at LHC [11] supports the non-thermal nature of the observed fluctuation. The value of (δN (net) Q ) 2 in this experiment is suppressed compared with the one in the equilibrated hadronic medium which has been calculated by lattice QCD simulations [12] and the hadron resonance gas (HRG) model [13]. Moreover, the dependence of (δN (net) Q ) 2 on the size of the rapidity window, ∆η, shows that the suppression of (δN (net) Q ) 2 is more pronounced for larger ∆η. These experimental results are reasonably explained if one attributes the suppression to the survival of fluctuations generated in the primordial deconfined medium [1,3,4,14]. The charge fluctuation normalized by the total number of charged particles, (δN (net) Q ) 2 / N ch , is known to take 2 ∼ 3 times smaller value in the equilibrated deconfined medium than the hadronic one [3,4]. After hadronization, this small fluctuation approaches the equilibrated value in the hadronic medium. Since the variation of the local density of a conserved charge is achieved only through diffusion, the approach of the fluctuation to the equilibrated value becomes slower as the volume where the charge is counted becomes larger. (δN (net) Q ) 2 thus takes a smaller value as ∆η becomes larger, which is consistent with the experimental result at ALICE.
On the other hand, the value of (δN (net) Q ) 2 observed at RHIC energy is consistent with the one in the equilibrated hadronic medium [21,22]. The difference between RHIC and LHC energies indicates that the evolution of the fluctuation in the hot medium is qualitatively different between these energies.
In order to confirm the validity of the above explanation on the fluctuation measured at ALICE and clarify the origin of the qualitative difference between the experimental results at RHIC and LHC energies, it should be instructive to measure the ∆η dependences of various other fluctuation observables in addition to (δN ) 2 , is an experimentallyobservable conserved charge fluctuation [15,16], although neutral baryons are not directly observable. Because the diffusion of the baryon number in the hadronic phase is slower than that of the electric charge due to the large mass of its carriers, baryons, if the origin of the suppression of (δN In event-by-event analyses, one can also measure higher order cumulants of conserved charges [9,10] such as the fourth-order ones (N ) 4 c . The experimental analysis of these observables as functions of ∆η can obviously provide us more information on the time evolution of fluctuations in the hot medium. So far, however, systematic studies on the dynamical evolution of higher order cumulants in heavy ion collisions, whose results can be compared with their experimental observation, have not been carried out to the best of the authors' knowledge. The purpose of the present Letter is to make the first investigation on this issue using a simple but theoretically lucid model, and to make a prediction on the ∆η dependence of higher order cumulants in relativistic heavy ion collisions.

STOCHASTIC FORMALISM TO DESCRIBE DIFFUSIVE SYSTEMS
In relativistic heavy ion collisions with sufficiently large collision energy per nucleon, √ s NN , the hot medium created at mid-rapidity has an approximate boost invariance. Useful coordinates to describe such a system are the coordinate-space rapidity η and proper time τ . We denote the net number of a conserved charge per unit coordinate-space rapidity as n(η, τ ). In a class of experiments, event-by-event fluctuations of the charge at kinetic freezeout in a phase space corresponding to the rapidity window determined by the experiment are observed. The phase space approximately corresponds to a finite coordinate-space rapidity interval [1]. Assuming that the kinetic freezeout takes place at a certain proper time τ fo , the experimentally-observed conserved-charge number at mid-rapidity at RHIC and LHC is given by at τ = τ fo with the rapidity window of the detector ∆η.
In the following, we investigate the time evolution of higher order cumulants of Q(τ ). In a sufficiently large space-time scale where hydrodynamic equations at first order are applicable, the average of n(η, τ ) follows the diffusion equation where D is the diffusion constant in this coordinate system. In a boost-invariantly expanding system, D receives a factor τ −2 compared with the diffusion constant in Cartesian coordinate. In order to describe fluctuations around the solution of Eq.
(2), one may employ a stochastic model, in which the time evolution of the deterministic part satisfies Eq. (2). A choice of such stochastic models is the theory of hydrodynamic fluctuations [17,18], in which the hydrodynamic equations are promoted to Langevin equations with stochastic terms representing fast random forces arising from microscopic interactions. In the equation corresponding to the conservation law of a charge, Eq. (2), the derivative of the stochastic force, ∂ η ξ(η, τ ), is added to the right-hand side of Eq. (2) [14]. The equation is referred to as stochastic diffusion equation. Up to Gaussian fluctuations, property of ξ(η, τ ) is completely determined by the fluctuation-dissipation relation, which is obtained from the locality of ξ(η, τ ) and large time behavior of n(η, τ ) [17]. It is known that the stochastic equation determined in this way well describes Gaussian fluctuations in fluids [17]. However, extension of this formalism to treat higher order fluctuations is nontrivial. There is no unique generalization of the fluctuationdissipation relation to higher orders, or no a priori justification of such extensions.
Concerning the difficulty in the description of non-Gaussian fluctuations, it is worthwhile to note a theorem on Markov process, which states that stochastic forces in a Langevin equation for Markov process are of Gaussian when the equation describes stochastic variables which are continuous and vary continuously 1 [19]. Since the standard theory of hydrodynamic fluctuations describes a Markov process and the hydrodynamic variables are continuous, the theorem demands that ξ(η, τ ) be of Gaussian; to allow for nonzero higher order correlations of ξ(η, τ ), one must relax at least one of the two conditions, i.e. Markovian and the continuity.
Without the higher order correlation of ξ(η, τ ), all higher order cumulants vanish in equilibrium. Even in such a formalism, the relaxation of non-Gaussianity starting from a particular initial condition can be described. In the physics of fluctuations in relativistic heavy ion collisions, however, nonzero higher order cumulants in equilibrium play a crucial role. First, the experimental results obtained so far report nonzero higher order cumulants near the equilibrated values [9][10][11]. Second, higher order cumulants are expected to increase toward the equilibrated values in the hadronic medium [3,4]. To reproduce these features, the stochastic model must obviously has nonzero higher order cumulants in equilibrium.
In the present study, instead of directly extending the theory of hydrodynamic fluctuations, we investigate the time evolution of higher order cumulants starting from a microscopic model. In this exploratory analysis, as such a model we consider a simple one-dimensional system composed of Brownian particles. Instead of tracking the motion of each Brownian particle separately, however, we represent the system as follows (See, Fig. 1). First, the coordinate η is divided into discrete cells with an equal length a. Second, we consider a single species of particle for the moment, and denote the number of particles in each cell, labeled by an integer m, as n m , and the probability that each cell contains n m particles as P (n, τ ) with n = (· · · , n m−1 , n m , n m+1 , · · · ). Finally, we assume that each particle moves to adjacent cells with a probability γ per unit proper time, as a result of microscopic interactions and random motion. The probability P (n, τ ) then follows the differential equation which is referred to as diffusion master equation in literature [19], where e m is the vector that all components are zero except for mth one, which takes unity. We will see later that in the continuum limit, a → 0, the average density and Gaussian fluctuation of n(η, τ ) in Eq. (3) agree with those in the stochastic diffusion equation with D = γa 2 . Before solving Eq. (3), some remarks are in order here. First, Eq. (3) is often solved with an approximation that n m are sufficiently large so that they can be regarded as continuous [19]. For the present purpose, however, the equation should be solved without this approximation, because once n m are set continuous, using an argument similar to the proof of the above-mentioned theorem on Markov processes one can show that all higher order cumulants become zero in equilibrium [19]. Discrete nature of n m should be kept to give rise to non-Gaussianity in equilibrium. Second, particles described by Eq. (3) behave as Brownian particles [20] without correlations with one another. From this property, it is immediately concluded that the distribution of Q in the τ → ∞ limit becomes of Poissonian, in which all cumulants take the same value Q n c = Q , provided that ∆η is sufficiently narrower than the total length of the system. Note that this behavior of cumulants is consistent with the hadronic fluctuations in the HRG model, owing to the non-interacting nature of hadrons in this model. Finally, since Eq. (3) describes a Markov process, stochastic effects in this model do not have temporal correlations. This model, therefore, would not be suitable to describe non-Gaussian fluctuations enhanced near the critical point, where the stochastic forces would have strong temporal correlation and thus the process becomes a non-Markov one due to critical slowing down. In the experimental results obtained so far, however, significant enhancements of cumulants which could be attributed to critical phenomena are not observed [9][10][11]. This result implies that critical phenomena do not come into play in the diffusion in the hot medium at least in the late stage in heavy ion collisions. In spite of the simplicity of the model Eq. (3), it is thus expected that the model qualitatively well describes the time evolution of Q n c in the hadronic stage.

SOLVING DIFFUSION MASTER EQUATION
Now, let us determine the time evolution of cumulants for the stochastic process Eq. (3). We first consider the time evolution of the probability P (n, τ ) with a fixed initial condition i.e. the initial particle numbers are fixed as n m = M m for all m without fluctuation. By introducing the factorial generating function, one obtains Equation (6) is a first-order partial differential equation, and solved with the method of characteristics. The solution with the initial condition Eq. (4) is given by with where N denotes the total number of cells and i is the imaginary unit. The factorial cumulants of the Fourier transform of n m ,ñ k = m n m e −2πkmi/N , are given by with K f (s, τ ) = log G f (s, τ ). Using Eqs. (10) and (7), the factorial cumulants up to the fourth order are calculated to be withM k = m M m e −2πkmi/N . The cumulants of n m are given by n m1 n m2 · · · n m l c = ∂ l K ∂θ 1 · · · ∂θ l θ=0 .
Since we are interested in the solution of Eq. (3) in the continuum limit, a → 0, here we introduce the shorthand notation for this limit: the particle density per unit rapidity n(η) = n m /a, with the rapidity of the mth cell η = ma, ω p = γa 2 p 2 with the conjugate momentum p = 2πk/N a. Also, the probability P (n, τ ) is promoted to a functional, which we denote as P [n(η), τ ]. Note, however, that these notations are conceptual; In actual applications, the functional P [n(η), τ ] is understood as the limit of the function P (n, τ ) with small but finite a. One finds from Eq. (11) that the continuum limit has to be taken with fixed D = γa 2 , so that the deterministic part of Eq. (3) follows Eq. (2). The factorial cumulants ofñ(p) in the continuum limit are obtained by simply replacing k with p in Eqs. (11) - (14).
In the following, we consider an infinitely long system without boundaries. With the fixed initial condition n(η, 0) = M (η), the cumulants of Eq. (1) at proper time τ are calculated to be with H (1) and I X (z) = where ∆η and τ dependences are encoded in the dimensionless parameter Next, we extend this result to general initial conditions containing fluctuations. We also extend the result to the system with two particle species with densities n 1 (η, τ ) and n 2 (η, τ ), and consider cumulants of the difference in order to compare the results with the cumulants of the net-electric charge and baryon numbers, which are given by the difference of particle numbers. When the two particle species separately follow Eq. (3), the probability with the initial condition P [n 1 (η), n 2 (η), 0] = F [M 1 (η), M 2 (η)] is given by the superposition of the solutions of fixed initial conditions using the superposition formula given in Refs. [23,24]. The general results will be presented in Ref. [24]. In the present Letter, we focus on the results for initial conditions satisfying spatial uniformity and locality, i.e.
where this equation defines [M i1 M i2 · · · M i l ] c on the right-hand side, which are generalized susceptibilities to higher orders and off-diagonal components. The condition Eq. (27) is satisfied, for example, in free gas, as well as systems well described by hydrodynamic equations, in equilibrium. With this initial condition, cumulants of Q (net) (τ ) are given by and M (net),(tot) (η) = M 1 (η) ∓ M 2 (η), respectively. F From Eq. (29), one also finds that the time evolution of the Gaussian fluctuation with the initial condition Eq. (27) is equivalent with the one in the stochastic diffusion equation [14] with the fluctuation in equilibrium Eq. (34). The details of the procedure to deal with Eq. (3) omitted in this Letter will be elucidated in the forthcoming publication [24].

DIFFUSION IN HADRONIC MEDIUM
Next, let us investigate the diffusion of higher order cumulants in the hadronic medium in relativistic heavy ion collisions using the above results. To make the argument simple, we assume that a boost invariant system with local equilibration is realized just above the critical temperature of the deconfinement transition. Then, the fluctuations at this time satisfy the locality condition Eq. (27). The hot medium then undergoes hadronization and chemical freezeout, which take place at almost the same time, τ = τ 0 , for large √ s NN . Because of the local charge conservations, fluctuations of local densities of the conserved charges are unchanged during these processes [1]. The fluctuations of conserved charges just after hadronization thus also satisfy Eq. (27) to a good approximation, with cumulants in the equilibrated deconfined medium, which are significantly smaller than the ones in the equilibrated hadronic medium [3,4,6]. We take this configuration at τ = τ 0 as the initial condition.
Due to the diffusion in the hadronic phase, the fluctuations approach the equilibrated values in the hadronic medium until kinetic freezeout at τ = τ fo . Provided that the τ dependence of the diffusion constant D in Eq. (2) is weak, this diffusion process is approximately described by Eq. , as conserved charges. Since experimental measurements of odd order cumulants of these charges are difficult at LHC energy because of their smallness, we also limit our attention to the second-and fourth-order cumulants. Because these charges in the hadron phase are predominantly carried by charged pions and nucleons, respectively, the cumulants of N at kinetic freezeout are approximately given by Eqs. (29) and (31) with τ = τ fo − τ 0 , where n 1 (η, τ ) and n 2 (η, τ ) are the densities of positive and negative pions, and nucleons and anti-nucleons, respectively.
To see the behavior of the cumulants determined by Eqs. (29) and (31), one must fix the parameters for the initial condition in these equations. As discussed above, the cumulants of the conserved charges at τ = τ 0 are suppressed compared with the equilibrated values. (tot) ] c in Eq. (31) is not constrained by the conservation laws and strongly depends on the hadronization mechanism. In the following, we treat this quantity as a parameter that characterizes the hadronization mechanism, and propose to utilize this parameter to constrain the ambiguity in it.
Since the value of [M 2 (tot) ] c is not constrained by the local charge conservation at hadronization as we discussed above, we regard the ratio as a free parameter and vary c in the range 0 ≤ c ≤ 1.5. The result on the second order, Q 2 (net) c , reproduces the one obtained by the stochastic diffusion equation [14].
In the figure, one finds that the behaviors of Q 4 (net) c and Q 2 (net) c as functions of 1/X are qualitatively different. While the behavior of Q 4 (net) c depends on the value of c, Q 4 (net) c is suppressed compared with Q 2 (net) c in the parameter range c 1.5. Note that Q 4 (net) c can become negative in a range of 1/X, although the fourthorder cumulants at the initial time and in equilibrium are both non-negative. It is interesting that the sign of the fourth-order cumulant can be flipped owing to the non-equilibrium effect.
Since 1/X is proportional to ∆η, the result in ) 2 c which has been already measured [11]. This statement is, however, altered for c 2. The same conclusion is also anticipated for the relation between the baryon number cumulants, (N ) 4 c . The suppression of the fourth-order cumulant compared with the second-order one can be intuitively understood as follows. In our analysis, we consider time evolution with the initial condition with small fluctuation. With this initial condition, the probability distribution becomes wider as τ becomes larger, and thus longer time is required until the probability in the tail is equilibrated. Because the higher order cumulants are sensitive to the tail behavior of the probability distribution, the approach of the cumulants to the equilibrated values is slower for higher order ones.
The result in Fig. 2 is obtained with the idealized initial condition Eq. (35), where fluctuations of the conserved charges completely vanish at τ = τ 0 . Next, we see how this result is modified when these cumulants have small but nonzero values at τ = τ 0 . In Fig. 3 we show 1/X dependences of Q 2 (net) c and Q 4 (net) c with the initial condition, The figure shows that the suppression of Q 4 (net) c compared with Q 2 (net) c is not as clear as the result in Fig. 2, but Q 4 (net) c < Q 2 (net) c is still satisfied for a rather wide range of c. The suppression of Q 4 (net) c thus is a rather robust feature reflecting the small initial fluctuations. Therefore, this behavior of the cumulants at LHC energy can be used as an experimental probe to confirm the suppression of the fluctuations at hadronization.
The results in Figs. 2 and 3 also show that the cumulants Q 2 (net) c and Q 4 (net) c have characteristic behaviors as functions of ∆η depending on the initial condition. These results indicate that experimental measurements of not only the magnitudes of various cumulants at a fixed ∆η but also their ∆η dependence enable us to explore various information on the time evolution of the hot medium and the hadronization mechanism in the experiments. In particular, these analyses would enable us to estimate the magnitude of the parameter c. Because this parameter is sensitive to the hadronization mechanism, such experimental information would be used as an important clue to it.

DISCUSSIONS
In this study, we employed a simple model, Eq. (3), to investigate the time evolution of non-Gaussian fluctuations. The model is expected to describe well the qualitative feature of the diffusion of higher order cumulants in the hadronic stage as already discussed. In order to explore the diffusion in each stage in heavy ion collisions more quantitatively, however, one must take various effects into account. In the present model, for example, fluctuations in equilibrium are given by the Poisson or Skellam distribution as a result of the absence of interactions between each Brownian particle. The model thus is not suitable to describe systems where equilibrated fluctuations do not follow either of these distributions. The easiest way to treat non-Poissonian fluctuations is to consider a system composed of several particle species having different charges which separately follow Eq. (3). The fluctuations of total charge of all particle species then becomes neither Poissonian nor Skellam distribution. This modification would be used as a first approximation to model the non-Skellam behavior of the net-electric charge fluctuations. When the non-Poissonian behavior comes from interaction between particles, terms describing the interaction should be introduced in Eq. (3). Next, the model Eq. (3) describes a Markov process and the stochastic effect does not have temporal correlations. To take account of nonzero temporal correlations, the model should be extended to non-Markov ones. Such a treatment would be required in, for example, dealing with the non-Gaussianity associated with critical phenomena and considering phenomena of the time scale of microscopic interaction. In the present study, it is also assumed that the numbers of two particle species n 1 (η, τ ) and n 2 (η, τ ) which carry opposite charges are separately conserved. While this assumption would be well justified for baryon numbers after chemical freezeout, effects of pair-creation and annihilation of electric charges would modify the time evolution of the net-electric charge fluctuations qualitatively. Finally, whereas the system investigated in the present study corresponds to the one with fixed diffusion constant D, τ dependence of this parameter should also be taken into account especially when one investigates the expanding system [18].
Besides these technical issues on the model, several assumptions also have been introduced on the geometry of the system and initial condition. To simplify the arguments, in the present Letter we limited our attention to the system which is infinitely long and the event-by-event configuration of the initial distribution is uniform and local, and thus satisfy Eq. (27). The hot medium created by heavy ion collisions, on the other hand, is a finite system, and the effects of global charge conservation and the violation of boost invariance modify the values of the fluctuation observables [1] and render them dependent on the position of the rapidity window. The validity of locality condition in Eq. (27) should also be investigated carefully. To investigate these effects, the model Eq. (3) has to be solved with appropriate boundary conditions and various initial conditions. Part of these issues will be investigated in the forthcoming paper [24].
In the present study, we limited our attention to sufficiently large √ s NN , where the measurement of odd order cumulants are difficult owing to their smallness. On the other hand, in the range of √ s NN where the beam energy scan program at RHIC is exploring, the third order cumulants are also observable [9,10] and will carry significant information on the QCD phase transition [7]. The time evolution and ∆η dependence of the third-order cumulants will also be discussed in Ref. [24].
In the present Letter, we investigated the time evolution of non-Gaussian fluctuations of conserved charges in a volume through diffusion using the diffusion master equation Eq. (3), which is intended to model the time evolution of higher order cumulants of conserved charges N ) 2 c at ALICE [11] is a consequence of the survival of the smallness of the fluctuation at the primordial stage. The same conclusion is also anticipated for baryon number cumulants. It should, however, be noted that this result qualitatively depends on the parameter c in Eq. (36), which is not constrained only by the local charge conservation and sensitive to hadronization mechanism. We emphasize that the fourth-order cumulants of net-electric charge and baryon numbers, (N ) 4 c , can be observable as functions of ∆η in the mid-rapidity region in heavy ion collisions at LHC and RHIC. The comparison of the experimental data on various cumulants with the present analysis will reveal the dynamical evolution of fluctuations and hadronization mechanism in heavy ion collisions.