Dynamical evolution of critical fluctuations and its observation in heavy ion collisions

TL;DR

This paper investigates how critical fluctuations evolve over time near the QCD critical point in heavy ion collisions, proposing that their rapidity dependence can serve as a robust experimental signature.

Contribution

It introduces a stochastic diffusion model to describe the evolution of critical fluctuations and identifies non-monotonic rapidity dependence as a key observable.

Findings

Non-monotonic rapidity dependence of cumulants signals critical fluctuations

Diffusion property allows probing early-time fluctuations

Robust experimental signatures for QCD critical point

Abstract

We study time evolution of critical fluctuations of conserved charges near the QCD critical point in the context of relativistic heavy ion collisions. A stochastic diffusion equation is employed in order to describe the diffusion property of the critical fluctuation arising from the coupling of the order parameter field to conserved charges. We show that the diffusion property gives rise to a possibility of probing the early time fluctuations through the rapidity window dependence of the second-order cumulant and correlation function of conserved charges. It is pointed out that their non-monotonic behaviors as functions of the rapidity interval are robust experimental signals for the existence of the critical enhancement around the QCD critical point.

Equations

$$Q_{\Delta y}(\tau )={\int }_{-\Delta y/2}^{\Delta y/2}dyn(y,\tau ).$$

$$⟨Q_{\Delta y}(\tau {)}^2{⟩}_{\mathrm{c}}$$

$$={\int }_{-\Delta y/2}^{\Delta y/2}d{y}_{1}d{y}_2⟨\delta n(y_1,\tau )\delta n(y_2,\tau )⟩,$$

$$⟨Q_{\Delta y}(\tau {)}^2{⟩}_{\mathrm{c}}={\int }_{-\Delta y}^{\Delta y}d\bar{y}(\Delta y-|\bar{y}|)⟨\delta n(\bar{y},\tau )\delta n(0,\tau )⟩.$$

$$\underset{\Delta y\to 0}{\lim }\frac{d}{d\Delta y}\frac{⟨Q_{\Delta y}(\tau {)}^2{⟩}_{\mathrm{c}}}{\Delta y}=\underset{\bar{y}\to 0}{\lim }⟨\delta n(\bar{y},\tau )\delta n(0,\tau )⟩.$$

$$\frac{\partial }{\partial \tau }\delta n(y,\tau )=D_y(\tau )\frac{{\partial }^2}{{\partial }^{2}y}\delta n(y,\tau )+\frac{\partial }{\partial y}\xi (y,\tau ),$$

$$⟨\xi (y_1,{\tau }_1)\xi (y_2,{\tau }_2)⟩$$

$$n(q,\tau )$$

$$\displaystyle d(\tau_{1},\tau_{2})=\Big{[}2\int_{\tau_{1}}^{\tau_{2}}d\tau^{\prime}D(\tau^{\prime})\Big{]}^{1/2}$$

$$⟨\delta n(q_1,\tau )\delta n(q_2,\tau )⟩$$

$$=$$

$$+{\int }_{{\tau }_0}^{\tau }d{\tau }_{1}d{\tau }_2⟨i{q}_1\xi (q_1,{\tau }_1)i{q}_2\xi (q_2,{\tau }_2)⟩$$

$$\times e^{-q_1^2[d({\tau }_1,\tau ){]}^2/2}{e}^{-q_2^2[d({\tau }_2,\tau ){]}^2/2},$$

$$⟨\delta n(y_1,{\tau }_0)\delta n(y_2,{\tau }_0)⟩=\chi ({\tau }_0)\delta (y_1-y_2).$$

$$⟨\delta n(q_1,\tau )\delta n(q_2,\tau )⟩$$

$$=$$

$$\displaystyle+2q_{1}^{2}\int_{\tau_{0}}^{\tau}d\tau^{\prime}\chi(\tau^{\prime})D(\tau^{\prime})e^{-q_{1}^{2}[d(\tau^{\prime},\tau)]^{2}}\bigg{)}.$$

$$⟨\delta n(y_1,\tau )\delta n(y_2,\tau )⟩$$

$$=$$

$$\displaystyle+\int_{\tau_{0}}^{\tau}d\tau^{\prime}\chi(\tau^{\prime})\frac{d}{d\tau^{\prime}}G\big{(}y_{1}-y_{2};2d(\tau^{\prime},\tau)\big{)}$$

$$=$$

$$\displaystyle-\int_{\tau_{0}}^{\tau}d\tau^{\prime}\chi^{\prime}(\tau^{\prime})G\big{(}y_{1}-y_{2};2d(\tau^{\prime},\tau)\big{)},$$

$$G(\bar{y};d)=\frac{1}{\sqrt{\pi }d}{e}^{-{\bar{y}}^2/d^2}.$$

$$\displaystyle\frac{\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}}{\Delta y}=\chi(\tau)-\int_{\tau_{0}}^{\tau}d\tau^{\prime}\chi^{\prime}(\tau^{\prime})F\Big{(}\frac{\Delta y}{2d(\tau^{\prime},\tau)}\Big{)},$$

$$F(X)$$

$$=\mathrm{erf}(X)+\frac{e^{-X^2}-1}{\sqrt{\pi }X}$$

$$\underset{X\to 0}{\lim }F(X)=0,\underset{X\to \infty }{\lim }F(X)=1.$$

$$\frac{⟨Q_{\Delta y}(\tau {)}^2{⟩}_{\mathrm{c}}}{\Delta y}\underset{\Delta y\to 0}{\to }$$

$$\frac{⟨Q_{\Delta y}(\tau {)}^2{⟩}_{\mathrm{c}}}{\Delta y}\underset{\Delta y\to \infty }{\to }$$

$$=$$

$$\displaystyle\chi^{\prime}(\tau)\left\{\begin{array}[]{cc}\geq 0\\ \leq 0\end{array}\right.\quad\Rightarrow\quad\frac{d~{}}{d\Delta y}\frac{\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}}{\Delta y}\left\{\begin{array}[]{cc}\leq 0\\ \geq 0\end{array}\right..$$

$$\displaystyle\boxed{\begin{array}[]{c}\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y\\ \mbox{has a local maximum}\\ \mbox{as a function of $\Delta y$}\end{array}}\Rightarrow\boxed{\begin{array}[]{c}\chi(\tau)\\ \mbox{has a local maximum}\\ \mbox{as a function of $\tau$}\end{array}}$$

$$\displaystyle\chi^{\prime}(\tau)\left\{\begin{array}[]{cc}\geq 0\\ \leq 0\end{array}\right.\quad\Rightarrow\quad\frac{d~{}}{d\bar{y}}\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle\left\{\begin{array}[]{cc}\geq 0\\ \leq 0\end{array}\right.,$$

$$\displaystyle\boxed{\begin{array}[]{c}\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle\\ \mbox{has a local minimum}\\ \mbox{as a function of $\bar{y}$}\end{array}}\Rightarrow\boxed{\begin{array}[]{c}\chi(\tau)\\ \mbox{has a local maximum}\\ \mbox{as a function of $\tau$}\end{array}}$$

$$\chi (T)$$

$$\frac{1}{D(T)}$$

$$M(R,\theta )$$

Full text

††preprint: J-PARC-TH-0088

I Introduction

The search for the QCD critical point is one of the most intriguing topics in the physics of QCD in medium. The existence of the critical point(s) in the QCD phase diagram in the temperature ($T$) and baryon chemical potential ($\mu$) plane is suggested by effective models Asakawa:1989bq ; Barducci:1989wi ; Halasz:1998qr ; Berges:1998rc ; Scavenius:2000qd ; Hatta:2002sj and lattice QCD Monte Carlo simulations Fodor:2001pe ; deForcrand:2002hgr ; Endrodi:2011gv . However, the number Kitazawa:2002bc ; Hatsuda:2006ps , locations, and even existence itself are still controversial. The experimental search for the critical point is one of the fundamental purposes of the relativistic heavy ion collision experiments Luo:2017faz . Because the medium created by these experiments pursues different trajectories in the $T$–$\mu$ plane depending on collision energies $\sqrt{s_{{}_{NN}}}$, this search would be achieved by comparing the characteristics of collision events at various $\sqrt{s_{{}_{NN}}}$. For this purpose, active experimental studies in the Beam-Energy Scan (BES) program at the Relativistic Heavy-Ion Collider (RHIC) are ongoing BESII . Future heavy-ion program at J-PARC J-PARC-HI , FAIR FAIR and NICA NICA will also contribute to this subject.

In the experimental search for the QCD critical point, fluctuation observables, such as the cumulants of the net-baryon number, are believed to be natural and promising observables Asakawa:2015ybt . At the critical point in the thermodynamic limit, equilibrated fluctuations of various observables diverge reflecting the softening of the effective potential. Such singular behavior is expected to be found in event-by-event analyses in heavy ion collisions Stephanov:1999zu ; Jeon:2003gk ; Stephanov:2008qz ; Asakawa:2009aj . The experimental analyses of the fluctuation observables, especially focusing on non-Gaussianity, have been performed actively recently Adamczyk:2013dal ; Adamczyk:2014fia ; Luo:2017faz .

In interpreting the experimentally-observed fluctuations, however, it should be remembered that the fluctuations observed in these experiments are not those in an equilibrated medium near the critical point. First, the more closely the system approaches the critical point, the longer becomes the relaxation time toward equilibrium because of the critical slowing down Berdnikov:1999ph . Owing to this effect, the enhancement of fluctuations is limited even if the medium passes right through the QCD critical point. Second, the critical fluctuations are to be distorted during the subsequent evolution until the detection Kitazawa:2013bta ; Sakaida:2014pya . Remember that particles are measured only in the final state. To extract the information on the critical point from the experimental data, proper understanding and description of these effects are indispensable, in addition to various restrictions from the real experimental settings Asakawa:2015ybt .

The time evolution of critical fluctuations has been discussed in the literature Berdnikov:1999ph ; Nonaka:2004pg ; Kapusta:2012zb ; Mukherjee:2015swa ; Herold:2016uvv . In Refs. Berdnikov:1999ph ; Nonaka:2004pg ; Mukherjee:2015swa , the effect of the critical slowing down is discussed for fluctuation of a uniform order parameter field $\sigma=\langle\bar{q}q\rangle$. It is known, however, that the critical mode of the QCD critical point is not the pure $\sigma$ mode, but rather is given by a linear mixing of the $\sigma$ mode and conserved charges Fujii:2003bz ; Fujii:2004jt ; Son:2004iv ; Minami:2011un . The soft mode of the critical point thus is a diffusion mode. As we will see below, the time evolution of a diffusion mode depends on the length scale, and this property is crucial for the description of its dynamics. One of the main goals of the present study is to reveal the dynamical evolution of the critical fluctuation with the diffusion property.

Another aim of this paper is to make a connection between the critical fluctuation of thermodynamics and experimental observation. This connection will be clearest for the fluctuations of conserved charges, especially that of the net-baryon number (see Ref. Asakawa:2015ybt for detailed discussions). In order to understand the evolution of the conserved-charge fluctuations, we first need to know that of the charge density itself in a collision event. This problem has been discussed for a hadronic medium in a simple setup Kitazawa:2013bta ; Sakaida:2014pya , but to the best of our knowledge it has not yet been studied for the case where a critical enhancement is encountered during the evolution. Carrying out this subject is the second goal of this study.

In this study, in order to address these issues we employ a simple stochastic diffusion equation (SDE). In this approach, the singularity associated with the critical point is encoded in the time-dependent susceptibility and diffusion coefficient. The SDE is a counterpart of the stochastic hydrodynamics Kapusta:2011gt , and is a good model to describe fluctuations of conserved charges. Moreover, as shown in Refs. Fujii:2003bz ; Fujii:2004za ; Fujii:2004jt ; Son:2004iv , the critical mode at the QCD critical point has the diffusion character and should be described by the SDE.

In this study, we first write down a formal solution for the density fluctuation within rapidity window $\Delta y$. Next, with a phenomenological parametrization for the susceptibility and diffusion coefficient, we analyze numerically the time evolution of the fluctuation. We find that the time evolution of the fluctuations strongly depends on the size of the rapidity window $\Delta y$. We also show that the second-order cumulant can have a non-monotonic $\Delta y$ dependence only when the medium fluctuation undergoes a critical enhancement during the time evolution. It is argued that this non-monotonic behavior serves as a robust experimental signal of the critical enhancement.

In the present study we restrict our attention only to the second-order cumulant, and leave the analysis of third and still higher order cumulants for future study, because even at the second order non-trivial outcomes are obtained from the analysis. We also deal with the correlation functions of conserved charges. We show that the same argument on the non-monotonic behavior holds for this function with rapidity window replaced by rapidity separation. The relation between the cumulant and the correlation function is also studied in detail.

This paper is organized as follows. In the next section we introduce the SDE, and argue that this model is suitable for describing the dynamics of conserved-charge fluctuations including the critical ones. We then solve this equation analytically, and discuss general properties in Sec. III. In Sec. IV, we introduce a phenomenological model, which is numerically solved in Sec. V. The last section is devoted to discussions and a summary. In Appendix A, we give a brief review of Refs. Fujii:2003bz ; Fujii:2004jt ; Son:2004iv , which clarified the diffusion property of the critical fluctuation. In Appendix B, we compare the cumulant and correlation function and discuss the condition for the appearance of the non-monotonic behaviors in these functions in detail.

II Model

II.1 Second-order cumulant and correlation function

In this study, we investigate the time evolution of the second-order cumulant and the correlation function of conserved charges in the hot medium produced in heavy ion collisions. We assume a boost-invariant Bjorken model for the event evolution throughout this paper in order to simplify the analysis, and adopt the Milne coordinates, the spacetime rapidity $y=\tanh^{-1}(z/t)$ and proper time $\tau=\sqrt{t^{2}-z^{2}}$.

Let us consider a conserved-charge density per unit rapidity $n(y,\tau)$, where transverse coordinates have been integrated out. As the conserved charge, we can take specifically the net-baryon number Kitazawa:2011wh or net-electric charge in heavy ion collisions. The amount of the charge in a finite rapidity interval $\Delta y$ at mid-rapidity at proper time $\tau$ is given by

[TABLE]

The second-order cumulant, or variance, of $Q_{\Delta y}(\tau)$ is written as

[TABLE]

where $\langle\delta n(y_{1},\tau)\delta n(y_{2},\tau)\rangle$ is the correlation function with $\delta n(y,\tau)=n(y,\tau)-\langle n(y,\tau)\rangle$ and $\langle\cdot\rangle$ stands for event average. In a boost invariant system, the correlation function depends only on the rapidity difference $\bar{y}=y_{1}-y_{2}$, and Eq. (2) is rewritten as

[TABLE]

This reveals a close relation between $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}$ and $\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle$. We investigate both of these functions in this paper. From Eq. (3), one finds that

[TABLE]

II.2 Stochastic diffusion equation

We consider the time evolution of the conserved charge density $n(y,\tau)$ at long time and length scales. This is well described by the stochastic diffusion equation Kapusta:2011gt ; Asakawa:2015ybt , which is written in the $\tau$–$y$ coordinates as

[TABLE]

where the diffusion coefficient $D_{y}(\tau)$ is related to the Cartesian one $D_{\rm C}(\tau)$ as $D_{y}(\tau)=D_{\rm C}(\tau)\tau^{-2}$. The noise term $\xi(y,\tau)$ represents the coupling with “short-time” fluctuations, whose average should vanish $\langle\xi(y,\tau)\rangle=0$. The noise term appears with the $y$-derivative in Eq. (5) so as to satisfy the conservation constraint.

When the noise correlation is local in time and space, the fluctuation-dissipation relation specifies its value to be Asakawa:2015ybt

[TABLE]

Here, $\chi_{y}(\tau)$ denotes the susceptibility of the conserved charge per unit rapidity and is related to the one in the Cartesian coordinates $\chi_{\rm C}(\tau)$ as $\chi_{y}(\tau)/\tau=\chi_{\rm C}(\tau)$. The susceptibility $\chi_{y}(\tau)$ is related to the second-order cumulant in equilibrium as $\langle Q_{\Delta y}^{2}\rangle_{\rm c,eq}=\chi_{y}\Delta y$. Furthermore, the noise is independent of the value of the density at earlier times: $\left<n(y,\tau)\xi(y^{\prime},\tau^{\prime})\right>=0$ for $\tau\leq\tau^{\prime}$. In the rest of this paper, we suppress the subscripts of the diffusion coefficient and susceptibility, and denote $D(\tau)=D_{y}(\tau)$ and $\chi(\tau)=\chi_{y}(\tau)$, respectively.

The SDE (5) is not only a suitable equation for describing the time evolution of conserved charges in a non-critical system, but also a good phenomenological equation to deal with the slow dynamics near the QCD critical point. This is because there the critical mode is identified as a linear combination of the $\sigma$ mode and conserved charges, and therefore its evolution must be consistent with the conservation law Fujii:2003bz ; Fujii:2004jt ; Son:2004iv ; Kapusta:2012zb . As a consequence, the equation for the critical mode is given by the same form as in Eq. (5). See Appendix A for more detailed discussion. Hence, in the present study we use Eq. (5) solely to describe the time evolution throughout the event trajectories passing near and away from the critical point.

In our study, the critical enhancement and slowing down of fluctuations are represented by the $\tau$-dependent susceptibility $\chi(\tau)$ and diffusion coefficient $D(\tau)$: Near the critical point, $\chi(\tau)$ grows sharply and $D(\tau)$ becomes vanishingly small reflecting, respectively, the large fluctuation and critical slowing down.

In the next section, we derive the formal solution of Eq. (5) with $\tau$ dependent $\chi(\tau)$ and $D(\tau)$, and discuss its general property. We then model the trajectories of collision events in terms of possible $\tau$-dependence of the susceptibility $\chi(\tau)$ and diffusion coefficient $D(\tau)$ including the effect of the critical point, and analyze the time evolution numerically in Secs. IV and V.

II.3 Comment on critical slowing down

Before closing this section, we comment on the difference of the treatment of the critical fluctuation in the present study from previous ones Berdnikov:1999ph ; Nonaka:2004pg ; Mukherjee:2015swa . In Refs. Berdnikov:1999ph ; Nonaka:2004pg ; Mukherjee:2015swa , the time evolution of the uniform $\sigma$ field was analyzed as the slowest mode of the system111 In Refs. Berdnikov:1999ph ; Nonaka:2004pg , the evolution of the “correlation length” of the $\sigma$ mode is studied, which is equivalent to treating the spatially uniform $\sigma$ mode, as shown in Ref. Mukherjee:2015swa .. Since the $\sigma$ mode is non-conserved, it can relax locally in rapidity space and follows a relaxation equation without $\partial/\partial y$.

This assumption is in contrast to the conserved charge fluctuation discussed in the present study, which can relax only through diffusion. As discussed already, the evolution equation of the critical fluctuation has to be consistent with the conservation law Fujii:2004za ; Fujii:2004jt . In order to respect this property, the evolution has to be dealt with the equations like the SDE consistent with the continuity equation222 Strictly speaking, neglected in Eq. (5) is the coupling of the soft mode with the momentum-density modes, which plays an important role in describing the critical dynamics more precisely Minami:2011un . This effect is discussed in Ref. Kapusta:2012zb . . As we will see in the next section, the time evolution of the critical mode then depends on their length scale $\Delta y$. This fact makes the problem complicated because the critical mode can no longer be regarded as a spatially uniform mode. At the same time, however, this diffusion property opens a possibility to study the critical fluctuation through $\Delta y$ dependence in experiments, as we will see later.

We also remark that our model can describe the fluctuations throughout the time evolution of the hot medium including the critical region as well as late stages. This property is advantageous in understanding dynamics behind experimental observables.

III Analytic properties

In this section, we formally solve the SDE (5) and study general properties of the second-order cumulant and correlation function analytically.

III.1 Solution of SDE

Defining the Fourier transform of $n(y,\tau)$ via $n(q,\tau)=\int dy\ e^{-iqy}n(y,\tau)$, the formal solution of Eq. (5) with the initial condition $n(q,\tau_{0})$ at $\tau=\tau_{0}$ is obtained as

[TABLE]

where

[TABLE]

denotes the diffusion “length” in rapidity space from $\tau_{1}$ to $\tau_{2}$ with $\tau_{1}\leq\tau_{2}$. The diffusion length $d(\tau_{1},\tau_{2})$ is a monotonically increasing (decreasing) function of $\tau_{2}$ ($\tau_{1}$), satisfying the boundary condition $d(\tau,\tau)=0$. The correlation function at proper time $\tau$ is obtained by taking the average of the product of Eq. (7) as333 This procedure to solve the stochastic equation corresponds to Stratonovich integral Gardiner . There is alternative method called Ito stochastic integral. These two stochastic integrals give the same result for Eq. (5).

[TABLE]

where we have used $\langle n(y_{1},\tau_{0})\xi(y_{2},\tau)\rangle=0$ for $\tau_{0}\leq\tau$.

To proceed further, we assume that the initial fluctuation satisfies the locality condition,

[TABLE]

Indeed, this condition should hold in a thermal medium at the length scale at which the extensive property of thermodynamic functions is satisfied Asakawa:2015ybt . We then obtain $\langle\delta n(q_{1},\tau_{0})\delta n(q_{2},\tau_{0})\rangle=2\pi\delta(q_{1}+q_{2})\chi(\tau_{0})$, and Eq. (III.1) is calculated to be

[TABLE]

The correlation function in $y$ space is obtained from Eq. (III.1) as

[TABLE]

where $\chi^{\prime}(\tau)=d\chi(\tau)/d\tau$ and we have defined the normalized Gauss distribution

[TABLE]

We note that the correlation function depends on $D(\tau)$ only through the diffusion length $d(\tau^{\prime},\tau)$.

By substituting Eq. (13) into Eq. (3), the second-order cumulant is calculated to be

[TABLE]

where

[TABLE]

with the error function ${\rm erf}(x)=\tfrac{2}{\sqrt{\pi}}\int_{0}^{x}dze^{-z^{2}}$. The behavior of $F(X)$ is shown in Fig. 1. As is shown in the figure, $F(X)$ is a monotonically increasing function satisfying

[TABLE]

III.2 Properties of fluctuation observables

From Eqs. (15) and (13), we find several notable features in the rapidity dependences of the cumulant and the correlation function.

First, we consider the behavior of $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}$ in the small and large $\Delta y$ limits. Using Eq. (17), the cumulant in these limits is easily calculated to be

[TABLE]

These results show that $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ takes the local-equilibrium value $\chi(\tau)$ in the small $\Delta y$ limit, while it recovers the initial value in the opposite limit. This shows that the relaxation toward the equilibrium state is sufficiently fast as $\Delta y$ becomes smaller, but it becomes arbitrarily slow with increasing $\Delta y$. The latter means that equilibration of the conserved-charge fluctuation over the whole system cannot be achieved within a finite time, because it takes infinite time to transport a charge from one end to the other. The analysis here also implies that $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ with smaller (larger) $\Delta y$ bears the information of $\chi(\tau)$ at later (earlier) $\tau$. This suggests that one can study the $\tau$ dependence of $\chi(\tau)$ from the $\Delta y$ dependence of $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ Kitazawa:2013bta .

Second, when $\chi(\tau)$ increases (decreases) monotonically in $\tau$, $\chi^{\prime}(\tau)\geq 0$ ($\leq 0$), then $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ for a given $\tau$ is a monotonically decreasing (increasing) function of $\Delta y$:

[TABLE]

This can be easily shown from Eq. (15) and the fact that $F(X)$ is a monotonically increasing function. Taking the contraposition of Eq. (24), one concludes that $\chi(\tau)$ must have at least one extremum when $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ is non-monotonic as a function of $\Delta y$. In particular,

[TABLE]

The same argument also applies to the correlation function. From the fact that $G(\bar{y},d)$ monotonically decreases as $\bar{y}$ increases, it is again easy to show that

[TABLE]

for $\bar{y}>0$. From the contraposition of Eq. (36), one obtains

[TABLE]

The properties (31) and (43) are quite useful in extracting the $\tau$ dependence of $\chi(\tau)$ in relativistic heavy ion collisions. If the experimental results of $\langle Q_{\Delta y}(\tau)^{2}\rangle/\Delta y$ and/or $\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle$ show non-monotonic behavior as a function of rapidity, this immediately confirms the existence of non-monotonicity in $\chi(\tau)$ as a function of $\tau$. It is known that the susceptibilities of baryon number and electric charge have a peak structure along the phase boundary around the QCD critical point Hatta:2002sj ; Asakawa:2009aj . The peak in $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ or $\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle$ serves as an experimental signal for this critical enhancement. This is the most important conclusion of this paper.

It should be kept in mind that the inverses of Eqs. (31) and (43) do not necessarily hold. That is, even if $\chi(\tau)$ is a non-monotonic function of $\tau$, there is a possibility that $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ and $\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle$ are monotonic. Therefore, from the monotonic behavior of $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ and/or $\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle$ one cannot conclude anything about the $\tau$ dependence of $\chi(\tau)$. In Appendix B, we discuss the condition for the appearance of the non-monotonic behaviors in these functions in more detail.

III.3 Comment on higher order cumulant

Using SDE (5), it is possible to calculate the time evolution of third and still higher order cumulants and correlation functions. As is easily shown, however, the higher order correlation functions $\langle\delta n(y_{1},\tau)\delta n(y_{2},\tau)\cdots\delta n(y_{N},\tau)\rangle$, and accordingly the higher order cumulants $\langle Q_{\Delta y}(\tau)^{N}\rangle_{\rm c}$, too, vanish in the $\tau\to\infty$ limit for $N\geq 3$ Asakawa:2015ybt ; the fluctuations of $n(y,\tau)$ in equilibrium described by Eq. (5) obey the Gaussian distribution.

Because of this property, SDE (5) is not capable of describing the relaxation of higher order fluctuations toward nonzero non-Gaussian equilibrium values. In relativistic heavy ion collisions, observed higher order cumulants take values close to their nonzero equilibrium values Asakawa:2015ybt . If this result is a consequence of the relaxation process in the hadronic medium, the relaxation cannot be described by the SDE (5). This is one of the reasons why we limit our attention to the second-order cumulant and correlation function in the present study, in spite of the useful properties of the higher order cumulants Ejiri:2005wq ; Stephanov:2008qz ; Asakawa:2009aj ; Friman:2011pf . To describe the relaxation of higher-order cumulants toward nonzero non-Gaussianity, different approaches are needed. In Ref. Kitazawa:2013bta , for example, the non-interacting Brownian particle model is employed to describe this process.

IV Model of collision evolution

In the previous section we showed that non-monotonic behavior of $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ and/or $\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle$, if observed, is a direct experimental evidence for the existence of a peak structure in susceptibility $\chi(\tau)$. In the rest of this paper, we demonstrate the appearance of the non-monotonic behavior by studying the behavior of $\langle Q_{\Delta y}(\tau)^{2}\rangle_{\rm c}/\Delta y$ and $\langle\delta n(\bar{y},\tau)\delta n(0,\tau)\rangle$ with a phenomenological parametrization of $\chi(\tau)$ and $D(\tau)$ for a collision event evolution passing near and away from the QCD critical point. In this section we first introduce the model for $\chi(\tau)$ and $D(\tau)$. Then, the time evolution of fluctuation is studied in the next section.

In this study, we write the susceptibility and the diffusion coefficient at temperature $T$ as a sum of their singular and regular contributions:

[TABLE]

where $\chi^{\rm cr}(T)$ and $\chi^{\rm reg}(T)$ denote the singular and regular parts of susceptibility per unit rapidity, respectively. We also define the singular and regular parts of diffusion coefficients $D_{\rm C}^{\rm cr}(T)$ and $D_{\rm C}^{\rm reg}(T)$ in Cartesian coordinate. We then parametrize the map of the evolution time to the temperature $T=T(\tau)$ to fix the $\tau$ dependences.

IV.1 Singular part

First, we discuss the singular parts $\chi^{\rm cr}(T)$ and $D^{\rm cr}(T)$. It is known that the QCD critical point belongs to the same static universality class as the 3D Ising model. The magnetization $M$ of the Ising model as a function of the reduced temperature $r$ and the dimensionless magnetic field $H$ near the critical point is parametrized with the two variables $R\geq 0$ and $\theta$ in the linear parametric model Guida:1996ep ; Schofield:1969zza as

[TABLE]

where $r$ and $H$ are expressed as

[TABLE]

The critical point is located at $r=H=0$. The crossover $(r>0,H=0)$ and first-order transition $(r<0,H=0)$ lines correspond to $\theta=0$ and $|\theta|=\sqrt{3/2}$ with $R>0$, respectively. We adopt approximate values $\beta=1/3$ and $\delta=5$ for the Ising critical exponents Guida:1996ep . From Eq. (46), one can calculate the magnetic susceptibility as

[TABLE]

As the susceptibility of a conserved charge $\chi$ near the QCD critical point should share the same critical behavior as $\chi_{M}(r,H)$, we set Stephanov:2011pb ; Mukherjee:2015swa ; Bzdak:2016sxg

[TABLE]

with a dimensionless proportionality constant $c_{\rm c}$. The susceptibility in the hadronic medium $\chi^{\rm H}$ will be defined in Sec. IV.3. We fix the normalization constants $m_{0}$ and $h_{0}$ by imposing $M(r=-1,H=0^{+})=1$ and $M(r=0,H=1)=1$.

In reality, the finite system size effect in heavy ion collisions prevents the divergence of $\chi$ Stephanov:1999zu . Nevertheless we ignore this effect because the growth of fluctuation would be limited more severely by the finiteness of the evolution time owing to the critical slowing down Berdnikov:1999ph .

For determining $D_{\rm C}^{\rm cr}$, we employ the dynamic universality argument Hohenberg:1977ym . Since the QCD critical point belongs to the model H Son:2004iv in the classification of Ref. Hohenberg:1977ym , the singular part $D_{\rm C}^{\rm cr}$ scales with the correlation length $\xi$ as $D_{\rm C}^{\rm cr}\sim\xi^{-2+\chi_{\eta}+\chi_{\lambda}}$, where the exponents $\chi_{\eta}$ and $\chi_{\lambda}$ for model H are obtained by the renormalization group calculation as $\chi_{\eta}\simeq 0.04$ and $\chi_{\lambda}\simeq 0.916$ Hohenberg:1977ym . As the correlation length $\xi$ and the susceptibility $\chi^{\rm cr}$ are related as $\chi^{\rm cr}\sim\xi^{2-\chi_{\eta}}$, the singular part $D_{\rm C}^{\rm cr}$ can be expressed in terms of $\chi^{\rm cr}$:

[TABLE]

with a proportionality constant $d_{\rm c}$ having the dimension of diffusion coefficient. We set $d_{\rm c}=1$ fm in what follows. Note that $D_{\rm C}^{\rm cr}$ vanishes at the critical point, reflecting the critical slowing down.

In our model, we leave the strengths of the critical component $c_{\rm c}$ as a free parameter. The reduced temperature $r$ controls the distance of the trajectory from the critical point. We will vary $r$ to simulate the change of the collision energy in the next section.

IV.2 Parametrizing the medium evolution

To utilize the above universality argument for describing heavy ion collision events, we need a map between the Ising variables $(r,H)$ and the physical variables $(T,\mu)$ in QCD, in addition to a map from the proper time $\tau$ to $(T,\mu)$ at a given collision energy. To skip these mappings, we follow a simple approach adopted in Refs. Berdnikov:1999ph ; Mukherjee:2015swa : We assume that $(T,\mu)$ in QCD are linearly mapped to the Ising variables $(r,H)$ around the critical point, and that only $H$ changes while $r$ is fixed during the time evolution. We write the linear relation between $H$ and $T$ as

[TABLE]

with $T_{\rm c}$ being the critical temperature. The ratio $\Delta T/\Delta H$ controls the width of the critical region in the QCD phase diagram. To relate $T$ and $\tau$, we assume the one-dimensional Bjorken expansion and conservation of total entropy. The relation is then obtained as Mukherjee:2015swa

[TABLE]

where $c_{\rm s}^{2}$ is the sound velocity and $T_{0}$ is the initial temperature of the system at the initial proper time $\tau_{0}$. We ignore possible entropy production in the critical region in this work.

In our calculation we set the initial temperature $T_{0}=220$ MeV at proper time $\tau_{0}=1.0$ fm, the critical temperature of the QCD critical point $T_{\rm c}=160$ MeV, and the kinetic freeze-out temperature $T_{\rm f}=100$ MeV Kumar:2014tca , where we stop the evolution. The parameters for the critical region in Eq. (52) are set to $\Delta T/\Delta H=10$ MeV. For the sound velocity $c_{\rm s}^{2}$, we adopt $c_{\rm s}^{2}=0.15$, which is indicated in a lattice calculation in the transition region at $\mu=0$ Bazavov:2014pvz . In Eq. (53), the effect of transverse expansion is not taken into account. Transverse expansion makes the duration of the hadron phase shorter. Thus, the following calculation is likely to estimate the effect of diffusion stronger than the actual one. As we will see, the suppression of the diffusion will be advantageous to detecting the critical point in experiments.

IV.3 Regular $+$ singular

We assume that the susceptibility per unit rapidity $\chi(T)$ approaches a constant value $\chi^{\rm Q}$ ($\chi^{\rm H}$) at high (low) temperature, which we call quark-gluon plasma (hadronic) value. We note that the value of $\chi^{\rm Q}$ depends on the trajectory in the $(T,\mu)$ plane as well as the thermal property Asakawa:2000wh ; Jeon:2000wg ; see also Refs. Koch:2008ia ; Asakawa:2015ybt . We also note that the susceptibility per unit rapidity approaches a constant value in the late stage in heavy ion collisions because the particle abundances are fixed after the chemical freeze-out. In the present study, we use the value estimated in Ref. Asakawa:2000wh assuming entropy conservation,

[TABLE]

For the diffusion coefficient, we assume that the coefficient in the Cartesian coordinates approaches constant values, $D_{\rm C}^{\rm Q}$ and $D_{\rm C}^{\rm H}$, at high and low temperatures, respectively. We take $D_{\rm C}^{\rm Q}=2.0$ fm from an estimate in the lattice QCD calculation Aarts:2014nba and $D_{\rm C}^{\rm H}=0.6$ fm from Ref. Sakaida:2014pya .

The regular parts $\chi^{\rm reg}(T)$ and $D_{\rm C}^{\rm reg}(T)$ are then constructed by smoothly interpolating these values at high and low temperatures,

[TABLE]

with

[TABLE]

Here $\chi_{0}^{\rm Q,H}$ and $D_{0}^{\rm Q,H}$ are determined so that $\chi(T)$ and $D_{\rm C}(T)$ coincide with the presumed values $\chi^{\rm Q,H}$ and $D_{\rm C}^{\rm Q,H}$ at $T=T_{0,\rm f}$, respectively. We set the width of the crossover region, $\delta T=10$ MeV.

In Fig. 2, we plot the susceptibility $\chi(T)/\chi^{\rm H}$ and the diffusion coefficient $D_{\rm C}(T)$ as a function of $T$ for several values of $r$ and $c_{\rm c}$ with $d_{\rm c}=1$. The upper panel of Fig. 2 shows that $\chi(T)/\chi^{\rm H}$ for $r=0$ diverges at $T=T_{\rm c}$. The sharp peak around $T_{\rm c}$ remains even for $r=1$ with $c_{\rm c}=4$. The regular part Eq. (55), labeled “reg”, is also shown for comparison. The lower panel of Fig. 2 shows that $D_{\rm C}(T)$ with the singular part vanishes at $T=T_{\rm c}$ for $r=0$, which is a manifestation of the critical slowing down.

V Effects of criticality on observables

Now, we analyze the time evolution of the fluctuation and correlation function using the parametrization obtained in the previous section, and study the effect of the QCD critical point on observables.

One remark concerned with the experimentally-observed fluctuations is that the diffusion described by the SDE (5) proceeds in coordinate space, but the experimental measurements are performed in momentum space. The imperfect correlation between the two rapidities owing to thermal motion gives rise to the “thermal blurring” effect Ohnishi:2016bdf ; Asakawa:2015ybt . For nucleons, this effect increases the apparent diffusion length in rapidity space by about $0.25$ at and after the thermal freeze-out Ohnishi:2016bdf . We take account of this effect in the subsequent analyses.

In this section, we show the numerical results of the cumulant and correlation function in the following normalized forms,

[TABLE]

where $\langle Q_{\Delta y}^{2}\rangle_{\rm c,H}=\chi^{\rm H}\Delta y$ is the cumulant in the equilibrated hadronic medium.

V.1 Non-critical trajectory

First, we study the case without the singular parts by setting $\chi(\tau)=\chi^{\rm reg}(T(\tau))$ and $D(\tau)=D^{\rm reg}(T(\tau))$. This corresponds to the collision events which pursue a trajectory away from the critical point in the crossover region (or, in the QCD phase diagram without a first order phase transition). As in Fig. 2, $\chi^{\rm reg}(T)$ behaves monotonically as a function of $T$ in this case.

In Fig. 3, we show the results of $K(\Delta y)$ and $C(\bar{y})$ for several values of $T$ from the initial temperature $T_{0}=220$ MeV to the kinetic freeze-out $T_{\rm f}=100$ MeV, together with the result after the thermal blurring. Note that the result only after thermal blurring can be compared with experimental results. The other results are shown to understand the time evolution of these quantities.

At the initial time with $T=220$ MeV, $K(\Delta y)$ is given by a constant, while $C(\bar{y})$ vanishes, in accord with the locality condition Eq. (10). As $T$ is lowered, nontrivial structures emerge in these functions. As discussed in Sec. II, $K(\Delta y)$ at $\Delta y=0$ is equal to its thermal value, i.e. $\chi(T)/\chi^{\rm H}$, which increases monotonically with time in this non-critical case. $K(\Delta y)$ at nonzero $\Delta y$ follows this trend but the increase is slower because of the finite diffusion time. As a result, the cumulant of a conserved charge depends on $\Delta y$ strongly. We also note that $K(\Delta y)$ decreases monotonically with increasing $\Delta y$, which is consistent with the statement, Eq. (24), which tells us that $K(\Delta y)$ should be monotonic when $\chi(\tau)$ is monotonic. $K(\Delta y)$ approaches its initial value $\chi^{\rm Q}/\chi^{\rm H}$ as $\Delta y$ increases. Notice that the behavior of $K(\Delta y)$ after thermal blurring is qualitatively consistent with the experimental result of the second-order cumulant of net-electric charge observed at the Large Hadron Collider ALICE .

The $\Delta y$ dependence of the second-order cumulant has been studied in Refs. Shuryak:2000pd ; Kitazawa:2013bta ; Sakaida:2014pya . The analysis of these studies corresponds to the parameter choice $\delta T=0$ in Eq. (55), i.e. $\chi(\tau)$ jumps discontinuously at $T_{\rm c}$. Since our result is qualitatively unchanged from the previous one, the nonzero width $\delta T$ seems not crucial for the argument here.

The lower panel of Fig. 3 shows $C(\bar{y})$. One finds that this function also behaves monotonically as a function of $\bar{y}$, which is consistent with Eq. (36). We also notice that $C(\bar{y})$ always takes a negative value. This is directly confirmed by substituting $\chi^{\prime}(\tau)>0$ into Eq. (13).

V.2 Trajectory passing through the critical point

Let us examine the case where the trajectory in heavy ion collisions passes right through the critical point ($r=0$). In Fig. 4, we show the evolution of $K(\Delta y)$ and $C(\bar{y})$ along the critical trajectory ($r=0$) for several values of $T$ and after thermal blurring with $c_{\rm c}=4$.

In the upper panel of Fig. 4, $K(\Delta y)$ at $T=T_{\rm c}=160$ MeV shows a remarkable enhancement, which comes from the divergence of $\chi(\tau)$ at $T=T_{\rm c}$ as shown in Fig. 2. This figure shows, however, that the cumulant stays finite for nonzero $\Delta y$ even at the critical point. This result is a manifestation of the critical slowing down for conserved charges. We remark that the effect of the critical slowing down is dependent on $\Delta y$ as discussed in Sec. III.2. After passing through the critical point, the value of $K(\Delta y)$ at $\Delta y=0$ decreases rapidly in accordance with the suppression of $\chi(\tau)$, while the decrease at nonzero $\Delta y$ is slower because of the slower diffusion for the larger $\Delta y$. As a consequence, a non-monotonic structure appears in $K(\Delta y)$. In Fig. 4, the non-monotonic behavior continues to exist in $K(\Delta y)$ until the kinetic freeze-out time and even survives the thermal blurring. Therefore, the non-monotonic behavior of $K(\Delta y)$ can be observed experimentally in this case. As discussed in Sec. III.2, this non-monotonic behavior, if observed, is a direct signal for the existence of the critical enhancement of $\chi(\tau)$.

An important lesson to learn from this result is that the non-monotonic behavior of $K(\Delta y)$ can survive whereas the magnitude of fluctuation itself is almost smeared to the equilibrated hadronic value $K(\Delta y)=1$; the maximum value of $K(\Delta y)$ after thermal blurring is $K(\Delta y)\simeq 1.2$ at $\Delta y=0.75$. This result suggests that the study of the non-monotonicity of $K(\Delta y)$ is advantageous for the search of the critical enhancement than the value of $K(\Delta y)$ with fixed $\Delta y$. Therefore, it is quite an interesting experimental subject to analyze its $\Delta y$ dependence.

From the lower panel of Fig. 4, one can draw the same conclusion on the $\bar{y}$ dependence of $C(\bar{y})$: $C(\bar{y})$ at $\bar{y}\to 0$ changes from negative to positive around $T_{\rm c}$. Triggered by this behavior, the non-monotonic $\bar{y}$ dependence of $C(\bar{y})$ manifests itself. The non-monotonicity again survives thermal blurring, suggesting that it can be measured experimentally.

We note that similar non-monotonic behaviors of correlation functions are also observed in Ref. Kapusta:2012zb and that of the mixed correlation function in Ref. Pratt:2012dz . The appearance of the non-monotonicity in these studies is understood completely the same way as that in Sec. III.

Next, we consider the case of a weaker critical enhancement by setting $c_{\rm c}=1$, but still keeping $r=0$. We show the results in Fig. 5. Although results above $T=155$ MeV look almost the same as those for $c_{\rm c}=4$, the non-monotonicity of $K(\Delta y)$ disappears already at $T=100$ MeV. As the growth of the susceptibility with $c_{\rm c}=1$ is weaker, the signal is drowned out by the diffusion in the hadronic phase. This exemplifies that, as discussed in Sec. III.2, the absence of the non-monotonicity in $K(\Delta y)$ does not necessarily mean the absence of the peak structure in $\chi(\tau)$.

The lower panel of Fig. 5 shows the result for the correlation function. The figure shows that the non-monotonic behavior of $C(\bar{y})$ generated at the critical point survives at $T=100$ MeV and even the thermal blurring. This result suggests that the non-monotonic signal in $C(\bar{y})$ is observable even when it disappears in $K(\Delta y)$. In fact, as we will discuss in Appendix. B, the non-monotonicity in $C(\bar{y})$ is more sustainable than that in $K(\Delta y)$.

V.3 Trajectory passing near the critical point

Finally, let us study the time evolution of fluctuations for $r>0$, which corresponds to the case where the system undergoes a crossover transition. Shown in Fig. 6 are the results of $K(\Delta y)$ and $C(\bar{y})$ for $r=1$ and $c_{\rm c}=4$. The results are qualitatively the same as those in Fig. 4. By closely comparing these results, one finds that the non-monotonic signal with $r=1$ for $T\leq 150$ MeV is much clearer than that in Fig. 4, although $\chi(T)$ does not diverge with $r=1$.

There are two reasons behind this result. First, $D(\tau)$ for $r=1$ does not vanish because the trajectory does not pass right through the critical point (see, Fig. 2). Therefore, the critical slowing down for $r=1$ is less important than that for $r=0$, and the fluctuations can grow faster around $T_{\rm c}$. Second, $\chi(T)/\chi^{\rm H}$ for $r=1$ and $c_{\rm c}=4$ is larger than that for $r=0$ at $T\lesssim 155$ MeV in our parametrization, as seen in Fig. 2. Therefore, $\chi(T)$ for $r=1$ approaches $\chi^{\rm H}$ more slowly. This behavior makes the non-monotonic peaks in $K(\Delta y)$ and $C(\bar{y})$ more prominent. Note, however, that the second observation may be dependent on the parametrization of $\chi(T)/\chi^{\rm H}$.

This argument suggests that the non-monotonic signals can be observed even when the trajectory does not pass right through the critical point. Moreover, it is possible that the trajectory off the critical point is more favorable for the emergence of the non-monotonicity. However, the signal of the critical enhancement, of course, weakens and finally disappears as the trajectory departs further off the critical point. In our analysis with $c_{\rm c}=4$, the non-monotonic behavior of $K(\Delta y)$ and $C(\bar{y})$ disappears at $r\simeq 5$ and $8$, respectively.

VI Discussions and a summary

The most important conclusion of the present study is Eqs. (31) and (43), i.e. the non-monotonic behaviors of $K(\Delta y)$ and/or $C(\bar{y})$, if observed, are direct experimental signals of the critical enhancement in the susceptibility. Now, let us consider the application of this conclusion to real heavy ion collisions. Throughout this study we assumed the Bjorken space-time evolution. This assumption, however, is violated in lower energy collisions. In particular, in the energy range of the BES program at RHIC, one has to consider the effect of the violation severely. For lower energy collisions, the effects of global charge conservation Sakaida:2014pya arising from the finite system size would also affect fluctuation observables seriously. It is expected, however, that these effects do not alter our conclusions on non-monotonicity as long as the system size in rapidity space is sufficiently large compared to the diffusion length. First, the violation of the Bjorken picture makes the correspondence between the coordinate and momentum space rapidities worse. This makes the thermal blurring effect stronger Ohnishi:2016bdf ; Asakawa:2015ybt , and the experimental measurement of non-monotonicity would become difficult. Even in this case, however, the relations Eqs. (31) and (43) should hold, because the thermal blurring effect only acts to enhance the diffusion length Ohnishi:2016bdf . From the analysis in Ref. Sakaida:2014pya , it is expected that the fluctuation around mid-rapidity is not affected by the effect of global charge conservation for sufficiently high energy collisions as long as the system size is sufficiently large. However, our conclusion may need to be altered for collisions at BES energies and lower energies such as those at FAIR, NICA, and J-PARC.

In this study, we described the critical fluctuation by the stochastic diffusion equation (5). Although this model well describes sufficiently long and slow fluctuations, the fluctuations in heavy ion collisions may not be slow enough compared to the time scale of the medium evolution. To take account of these effects, the SDE (5) has to be modified. One direction is to include higher order derivative terms. Another interesting extension is to include the $\sigma$ field as a dynamical field, and solve the coupled equation of $n$ and $\sigma$. Near the critical point, the coupling of $\sigma$ with momentum density would also be important Kapusta:2012zb . To deal with these subjects, numerical simulations of medium evolution will be needed with adopting a certain model, for example, the chiral fluid model Herold:2016uvv .

In this study, we investigated the time evolution of the second order cumulant and the correlation function of conserved charges in heavy ion collisions which pass through or near the critical point focusing on the effects of critical slowing down near the critical point and dissipation in the late stages. We adopted the stochastic diffusion equation with critical nature being encoded in the time-dependent susceptibility and diffusion coefficient. This model can describe the dynamics of the critical mode respecting its diffusion property, which was not considered in previous studies on the critical slowing down. We have pointed out that the critical enhancement in susceptibility can be observed as the non-monotonic behaviors in the second-order cumulant and correlation function. Our numerical results suggest that these non-monotonic behaviors are a more robust experimental signal than the value of these functions themselves. It is, therefore, quite an interesting experimental subjects to analyze the rapidity dependences of these functions in heavy ion collisions.

M. A. and M. K. thank V. Koch and M. Lisa for inviting them to INT workshop “Exploring the QCD Phase Diagram through Energy Scans”, Sep. 19 - Oct. 14, 2016, Seattle, USA, and stimulating discussions. The authors thank M. Nahrgang and M. Bluhm for fruitful discussions. They also thank A. Bzdak, K. Redlich, and M. Stephanov for useful conversations. This work was supported in part by JSPS KAKENHI Grant Numbers 16J01314, 26400272, and 16K05343.

Appendix A Time evolution of the soft mode near the QCD critical point

In this appendix, we explain that the appropriate equation to describe the time evolution of the soft mode near the critical point is the SDE (5), on the basis of Refs. Hohenberg:1977ym ; Fujii:2004jt ; Son:2004iv .

At sufficiently long distance and time scales, the dynamics of a finite temperature system near equilibrium is described by hydrodynamic theory, which only contains the modes whose excitation energies vanish in the long wavelength limit – “hydrodynamic modes”. Near the critical point, the hydrodynamic variables are given by the fluctuations of the order parameter and the densities of conserved charges Fujii:2004jt ; Son:2004iv . The QCD critical point shares the same dynamical universality class with the model H in the classification of Ref. Hohenberg:1977ym , and the chiral order parameter field $\sigma=\langle\bar{q}q\rangle$ has nonzero couplings with baryon number density $n$, and energy and momentum densities. In this appendix, we neglect the energy-momentum density for a simple illustration.

We start from the Ginzburg-Landau free energy functional

[TABLE]

where the coefficients $A,\ B$, and $C$ are functions of temperature $T$ and baryon chemical potential $\mu$. The neglected terms include higher order terms in $\delta\sigma$ and $\delta n$, and derivative terms, which are not important to describe slow modes. Here, $B\neq 0$ because the coupling between $\sigma$ and $n$ is allowed at the QCD critical point because of the finite quark masses and finite baryon density Fujii:2004jt .

Deviation of $\sigma$ and $n$ from the equilibrium values gives rise to relaxation of the system. The evolutions of $\sigma$ and $n$ are given by the following stochastic hydrodynamic equations:

[TABLE]

where the noise correlators are local,

[TABLE]

with $i,j=\sigma,n$. From Onsager’s principle, we have $\gamma_{\sigma n}=\gamma_{n\sigma}$. In the small momentum limit, $\gamma_{\sigma\sigma}$ is given by a nonzero constant, while the coefficients for $\dot{n}$, $\gamma_{n\sigma}$ and $\gamma_{nn}$, are proportional to space-derivative squared because of the conservation law, parity invariance, and analyticity. We thus have in the Fourier space

[TABLE]

with $q^{2}$ being the momentum squared. Inserting Eqs. (71) and (A) into Eq. (69), we obtain the hydrodynamic equation to leading order in $q^{2}$ as

[TABLE]

Solving Eq. (78), we obtain two eigenfrequencies,

[TABLE]

with $\Delta=AC-B^{2}$ and the corresponding eigenmodes,

[TABLE]

which decay with $|\omega_{1}|$ and $|\omega_{2}|$, respectively. At the critical point, the energy functional develops a flat direction, $\Delta=0$, and the susceptibilities of $\sigma$ and $n$ become divergent. The diffusion coefficient

[TABLE]

goes to zero at the critical point, which represents the critical slowing down.

This result tells us that a small fluctuation of $\sigma$ and $n$ in the system relaxes with two distinct time scales. The first slow mode $v_{1}$ is just the conserved charge $n$, whose time evolution is described by the SDE

[TABLE]

Note that $\xi_{n}$ should be proportional to space-derivative as in Eq. (5) so that this equation is consistent with the continuity equation. On the other hand, $v_{2}$ is a relaxation mode with nonvanishing relaxation time scale $(\gamma_{\sigma\sigma}A)^{-1}$. This time scale is fast compared to that of $v_{1}=\delta n$. In the faster time scale, $\sigma$ alone relaxes to the value $\delta\sigma=-(B/A)\delta n$, which minimize $F[\sigma,n]$ with $n$ being fixed to a given value. After that, with the longer time scale $(Dq^{2})^{-1}$ the mode $v_{1}$ relaxes to $\delta n=0$. In this stage, $\sigma$ simply traces the profile of $n$, and the time evolution of the slow mode is described by the SDE (85). The effect of the critical point is encoded in vanishing of $D$ in the SDE.

In the above discussion we neglected the energy-momentum density. In model H, with which the QCD critical point shares the same dynamic universality, the hydrodynamic slow modes in the long wavelength limit are in fact given by the baryon number diffusion and the diffusion of two transverse momentum components Son:2004iv having a nonlinear coupling via the nonzero Poisson bracket Hohenberg:1977ym . The analysis incorporating the nonlinear coupling is left for future works.

Appendix B Conditions for the appearance of non-monotonicity

In Sec. III, we showed that the non-monotonic behaviors of $K(\Delta y)$ and $C(\bar{y})$ in Eqs. (58) and (59) serve as direct experimental evidence for the existence of a peak structure in $\chi(\tau)$. In this appendix, we take a much closer look at the conditions for the appearance of the non-monotonic behaviors in these functions and discuss which function is better in sustaining the non-monotonicity.

To simplify the problem, in this appendix we consider the functional form of $\chi(\tau)$ which has only one maximum as a function of $\tau$. Then, $\chi^{\prime}(\tau)$ changes its sign only once from positive to negative. In this case, $K(\Delta y)$ ($C(\bar{y})$) can have only one local maximum (minimum). Thus, the necessary and sufficient conditions for the non-monotonic behaviors of $K(\Delta y)$ and $C(\bar{y})$ are given by

[TABLE]

respectively.

From Eq. (15), the $\Delta y$ derivative of $K(\Delta y)$ is given by

[TABLE]

where, $F^{\prime}(X)=dF(X)/dX$. Using

[TABLE]

one obtains

[TABLE]

Here, $d(\tau^{\prime},\tau_{\rm f})$ is a monotonically decreasing function of $\tau^{\prime}$ with $d(\tau_{\rm f},\tau_{\rm f})=0$. The integrand in Eq. (90) is $\chi^{\prime}(\tau^{\prime})$ with a weight $1/d(\tau^{\prime},\tau_{\rm f})$, which takes larger value for larger $\tau^{\prime}$. The sign of $\chi^{\prime}(\tau^{\prime})$ with later $\tau^{\prime}$ is more strongly reflected to the sign of Eq. (90). On the other hand, in Eq. (91) $\chi^{\prime}(\tau^{\prime})$ is integrated with a weight $d(\tau^{\prime},\tau_{\rm f})$ taking larger value for earlier $\tau^{\prime}$. The sign of $\chi^{\prime}(\tau^{\prime})$ with earlier $\tau^{\prime}$ is more responsible for that of Eq. (91).

Next, the $\bar{y}$ derivative of $C(\bar{y})$ is calculated to be

[TABLE]

By taking the small $\bar{y}$ limit, we obtain

[TABLE]

Equation (93) shows that the sign of $dC(\bar{y})/d\bar{y}$ in the small $\bar{y}$ limit is determined by the integral of $\chi^{\prime}(\tau^{\prime})$ with a weight $1/d(\tau^{\prime},\tau_{\rm f})^{3}$. In the large $\bar{y}$ limit, the weight is given by $e^{-\bar{y}^{2}/d(\tau^{\prime},\tau_{\rm f})^{2}}/d(\tau^{\prime},\tau_{\rm f})^{3}$, which concentrates at the initial time $\tau^{\prime}=\tau_{0}$ in the large $\bar{y}$ limit. The sign in this limit thus is determined only by $\chi^{\prime}(\tau_{0})$, which is positive in the present situation.

Now, let us compare the conditions Eqs. (86) and (87). As discussed above, the second condition in Eq. (87) is always satisfied, while that in Eq. (86) is not necessarily true but dependent on the functional form of $\chi(\tau)$. Next, the first conditions in Eqs. (86) and (87) are not always satisfied, but the latter is more favored, because the weight of the latter, $d(\tau^{\prime},\tau_{\rm f})^{-3}$, is more concentrated at later $\tau^{\prime}$ than the former, $d(\tau^{\prime},\tau_{\rm f})^{-1}$. From these observations, one concludes that the appearance of the non-monotonicity is more likely to appear in $C(\bar{y})$ than in $K(\Delta y)$.

Although the manifestation of non-monotonicity is more robust in $C(\bar{y})$ than $K(\Delta y)$, in experimental analyses it is meaningful to analyze both of these functions. In the above argument, the position of the local extremum is not determined. The ranges of $\Delta y$ and $\bar{y}$ which can be measured in experiments are limited owing to the coverage of the detector, and the manifestation of the extremum in this range depends on the functional form. Therefore, by analyzing both $K(\Delta y)$ and $C(\bar{y})$, the chance to find the non-monotonic behaviors is enhanced.