Multiplicity fluctuations near the QCD critical point

TL;DR

This paper develops a framework to calculate statistical moments of particle multiplicities in heavy-ion collisions, accounting for experimental limitations and background effects, to better identify signals of the QCD critical point.

Contribution

It introduces a comprehensive method for evaluating multiplicity moments considering correlations, resonance decays, acceptance, and efficiency, aiding the search for the QCD critical point.

Findings

Framework accounts for critical and spurious correlations.

Includes effects of resonance decay and detector efficiency.

Focuses on second-order moments with straightforward extension to higher orders.

Abstract

Statistical moments of particle multiplicities in heavy-ion collision experiments are an important probe in the exploration of the phase diagram of strongly interacting matter and, particularly, in the search for the QCD critical end point. In order to appropriately interpret experimental measures of these moments, however, it is necessary to understand the role of experimental limitations, as well as background contributions, providing expectations on how critical behavior should be affected by them. We here present a framework for calculating moments of particle multiplicities in the presence of correlations of both critical and spurious origins. We also include effects from resonance decay and a limited acceptance window, as well as detector efficiency. Although we focus on second-order moments, for simplicity, an extension to higher-order moments is straightforward.

Equations

$$\mathcal{P}[\sigma ]\sim e^{-\Omega [\sigma ]/T},$$

$$\Omega [\sigma ]\approx \int d^{3}x\left\{\frac{1}{2}(\nabla \sigma {)}^2+\frac{1}{2}{m}_{\sigma }^2{\sigma }^2\right\},$$

$$P({\sigma }_0)\sim e^{-\frac{V}{T}\frac{m_{\sigma }^2}{2}{\sigma }_0^2}.$$

$${\mathcal{L}}_{int}=G{\sigma }_0\vec{\pi }\cdot \vec{\pi }+g{\sigma }_0{\bar{\psi }}_p{\psi }_p,$$

$$\frac{d\xi }{dt}=A{\left(\frac{\xi }{{\xi }_0}\right)}^{2-z}\left(\frac{{\xi }_0}{\xi }-\frac{{\xi }_0}{{\xi }_{eq}(t)}\right),$$

$$\begin{array}{rl}V(b,R_N)& =C\mathcal{A}(b,R_N)\\ & =2C{R}_N^2{\cos }^{-1}\left(\frac{b}{2R_N}\right)-b\sqrt{R_N^2-\frac{b^2}{4}},\end{array}$$

$$\begin{array}{rl}\omega & =\sqrt{p^2+m_0^2+\delta m^2}\\ & \approx {\omega }_0\left[1+\frac{1}{2}\frac{\delta m^2}{{\omega }_0^2}-\frac{1}{8}\frac{(\delta m^2{)}^2}{{\omega }_0^4}+\cdots \right]\\ & \approx {\omega }_0\left[1+\frac{m\delta m}{{\omega }_0^2}+\frac{1}{2}\left(1-\frac{m^2}{{\omega }_0^2}\right)\frac{(\delta m{)}^2}{{\omega }_0^2}+\cdots \right],\end{array}$$

$$p_i=$$

$$\begin{array}{rl}{\omega }_i=& \sqrt{p_i^2+m^2}\approx {\omega }_{0i}\{1-{\left(\frac{p_{0i}}{{\omega }_{0i}}\right)}^2\frac{\delta R}{R}+\\ & +\frac{1}{2}\left[3{\left(\frac{p_{0i}}{{\omega }_{0i}}\right)}^2-{\left(\frac{p_{0i}}{{\omega }_{0i}}\right)}^4\right]{\left(\frac{\delta R}{R}\right)}^2+\cdots \}.\end{array}$$

$$\frac{{\omega }_{0i}-\mu }{T+\delta T}=\frac{{\omega }_{0i}-\mu }{T}\left(1-\frac{\delta T}{T}+\frac{\delta T^2}{T^2}+\cdots \right),$$

$$\delta {\omega }_i^T=({\omega }_{0i}-\mu )\left(-\frac{\delta T}{T}+\frac{\delta T^2}{T^2}+\cdots \right).$$

$${\omega }_{0i}-(\mu +\delta \mu )={\omega }_{0i}+\delta {\omega }_i^{\mu }-\mu ,$$

$$\displaystyle\begin{split}&\langle n_{i}\rangle=f(\zeta_{i}):=(e^{-\zeta_{i}}\mp 1)^{-1}\,,\\ &\langle(\Delta n_{i})^{2}\rangle=f^{\prime}(\zeta_{i})=f(\zeta_{i})\;\big{(}1\pm f(\zeta_{i})\big{)}\,,\\ &\langle(\Delta n_{i})^{3}\rangle=f^{\prime\prime}(\zeta_{i})=f^{\prime}(\zeta_{i})\;\big{(}1\pm 2\,f(\zeta_{i})\big{)}\,,\\ &\langle(\Delta n_{i})^{4}\rangle=f^{\prime\prime\prime}(\zeta_{i})=f^{\prime}(\zeta_{i})\big{(}1\pm 6\,f^{\prime}(\zeta_{i})\big{)}\,,\end{split}$$

$$\begin{array}{rl}⟨n_i⟩\approx & f({\zeta }_{0i})-\beta f^{\prime }({\zeta }_{0i})\delta {\omega }_i+\\ & +\frac{{\beta }^2}{2}{f}^{\prime \prime }({\zeta }_{0i})(\delta {\omega }_i{)}^2+\cdots ,\end{array}$$

$$\begin{array}{rl}\overline{⟨n_i⟩}\approx & f({\zeta }_{0i})-\beta f^{\prime }({\zeta }_{0i})\overline{\delta {\omega }_i}+\\ & +\frac{{\beta }^2}{2}{f}^{\prime \prime }({\zeta }_{0i})\overline{(\delta {\omega }_i{)}^2}+\cdots ,\end{array}$$

$$\begin{array}{rl}\overline{⟨\Delta n_i\Delta n_j⟩}:=& \overline{⟨n_i{n}_j⟩}-\overline{⟨n_i⟩}\overline{⟨n_j⟩}\\ =& \overline{⟨n_i⟩⟨n_j⟩}-\overline{⟨n_i⟩}\overline{⟨n_j⟩}+{\delta }_{ij}\overline{f^{\prime }(-\beta {\omega }_i)}.\end{array}$$

$$\overline{\langle\Delta n_{i}\,\Delta n_{j}\rangle}\approx f^{\prime}(\zeta_{0\,i})\,f^{\prime}(\zeta_{0\,j})\;\beta^{2}\,\overline{\delta\omega_{i}\delta\omega_{j}}+\\ +\delta_{ij}\bigg{(}f^{\prime}(\zeta_{0\,i})-\beta\;f^{\prime\prime}(\zeta_{0\,i})\;\overline{\delta\omega_{i}}+\\ +\dfrac{\beta^{2}}{2}f^{\prime\prime\prime}(\zeta_{0\,i})\;\overline{(\delta\omega_{i})^{2}}\bigg{)}\,,$$

$$\begin{gathered} \overline{\delta {\omega }_i}\approx \frac{{\omega }_{0i}}{2}\left[3{\left(\frac{p_{0i}}{{\omega }_{0i}}\right)}^2-{\left(\frac{p_{0i}}{{\omega }_{0i}}\right)}^4\right]\overline{{\left(\frac{\delta R}{R}\right)}^2}+ \\ +\frac{g^2}{2{\omega }_{0i}}\left(1-\frac{m^2}{{\omega }_{0i}^2}\right)\overline{(\delta \sigma {)}^2}+({\omega }_{0i}-\mu )\overline{{\left(\frac{\delta T}{T}\right)}^2}, \end{gathered}$$

$$\begin{gathered} \overline{\delta {\omega }_i\delta {\omega }_j}\approx \frac{g^2{m}^2}{{\omega }_{0i}{\omega }_{0j}}\overline{(\delta \sigma {)}^2}-\frac{(p_{0i}{p}_{0j}{)}^2}{{\omega }_{0i}{\omega }_{0j}}\overline{{\left(\frac{\delta R}{R}\right)}^2} \\ +({\omega }_{0i}-\mu )({\omega }_{0j}-\mu )\overline{{\left(\frac{\delta T}{T}\right)}^2}+\overline{(\delta \mu {)}^2}. \end{gathered}$$

$$\overline{⟨N⟩}=\sum _i\overline{⟨n_i⟩},$$

$$\overline{⟨(\Delta N{)}^2⟩}=\sum _{i,j}\overline{⟨\Delta n_i\Delta n_j⟩}.$$

$$\Phi (R,\theta ,\varphi )=0.$$

$$\begin{array}{rl}\overline{⟨(\Delta N{)}^2⟩}=& \sum _{\begin{array}{c}{\ell }_1,m_1,i_1\\ {\ell }_2,m_2,i_2\end{array}}\overline{⟨\Delta n_{{\ell }_1,i_1}\Delta n_{{\ell }_2,i_2}⟩}\\ =& \sum _{\begin{array}{c}{\ell }_1,i_1\\ {\ell }_2,i_2\end{array}}(2{\ell }_1+1)(2{\ell }_2+1)A_{({\ell }_1,i_1),({\ell }_2,i_2)}+\\ & +\sum _{{\ell }_1,i_1}(2\ell +1)B_{(\ell ,i)},\end{array}$$

$$\begin{gathered} A_{({\ell }_1,i_1),({\ell }_2,i_2)}=f^{\prime }({\zeta }_{0({\ell }_1,i_1)})f^{\prime }({\zeta }_{0({\ell }_2,i_2)})\times \\ \times {\beta }^2\overline{\delta {\omega }_{({\ell }_1,i_1)}\delta {\omega }_{({\ell }_2,i_2)}}, \end{gathered}$$

$$\begin{gathered} B_{(\ell ,i)}=f^{\prime }({\zeta }_{0({\ell }_1,i_1)})-\beta f^{\prime \prime }({\zeta }_{0({\ell }_1,i_1)})\overline{\delta {\omega }_{({\ell }_1,i_1)}}+ \\ +\frac{{\beta }^2}{2}{f}^{\prime \prime \prime }({\zeta }_{0({\ell }_1,i_1)})\overline{(\delta {\omega }_{({\ell }_1,i_1)}{)}^2}. \end{gathered}$$

$$\begin{gathered} A_{({\ell }_1,i_1),({\ell }_2,i_2)}=f_A^{\prime }({\zeta }_{0({\ell }_1,i_1)}^A)f_B^{\prime }({\zeta }_{0({\ell }_2,i_2)}^A)\times \\ \times {\beta }^2\overline{\delta {\omega }_{({\ell }_1,i_1)}^A\delta {\omega }_{({\ell }_2,i_2)}^B}, \end{gathered}$$

$$B_{(\ell ,i)}=0,$$

$$\begin{array}{r}{u}^2\ge 1-{\left(\frac{p_{+}}{p}\right)}^2,\end{array}\begin{array}{r}{u}^2\le 1-{\left(\frac{p_{-}}{p}\right)}^2.\end{array}$$

$$e^{-2{\eta }_C}\le \frac{1+u}{1-u}\le e^{2{\eta }_C},$$

$$e^{-2{\eta }_C}-1\le u(1+e^{-2{\eta }_C})u(1+e^{2{\eta }_C})\le e^{2{\eta }_C}-1,$$

$$|u|\le \tanh {\eta }_C.$$

Full text

I Introduction

The QCD phase diagram in the temperature versus baryonic chemical potential plane is marked by very distinct regimes, characterized by different (approximate) symmetries and degrees of freedom. However, the nature of these regimes and the transitions between them is still not settled. While first-principle calculations face the breakdown of perturbation theory at low and intermediate energy scales and lattice calculations remain limited to relatively low baryon densities, where an analytic crossover is found Aoki et al. (2006); Borsanyi et al. (2010), effective models suggest a first-order phase transition at high densities. The first-order coexistence line, related to the restoration of chiral symmetry, should end at lower densities in a second-order critical end point of unknown location Stephanov (2006, 2004); Rajagopal (1995).

The possibility of experimentally probing this point with current technology and the very peculiar singular behavior associated with a second-order phase transition raise very appealing hopes of accessing reliable information on the main features of the phase diagram of strongly interacting matter, which have motivated considerable experimental and theoretical efforts. The search for the QCD critical point is currently ongoing at the RHIC Beam Energy Scan program Luo (2016, 2015); Xu (2016).

On the other hand, notwithstanding its very characteristic features in an ideal scenario, the observation or exclusion of the critical point in the noisy and dynamic context of heavy-ion collision experiments remains an extremely challenging task. Whereas the main probe in the pursuit for such a point is the increase of long-range fluctuations caused by a diverging correlation length in its neighborhood Stephanov et al. (1999); Stephanov (2009); Athanasiou et al. (2010), the short lifetime and small size of the formed system will considerably limit its growth Berdnikov and Rajagopal (2000); Braun et al. (2006); Kiriyama et al. (2006); Palhares et al. (2011); Fraga et al. (2011); Mukherjee et al. (2015, 2016); Hippert et al. (2016). Moreover, collisions are not perfectly controlled and are affected by initial-state fluctuations, so that spurious sources of fluctuations could possibly hide or mimic the expected signatures Stephanov et al. (1999); Hippert et al. (2016); Gorenstein (2015); Sahoo et al. (2013); Herold et al. (2016); Luo et al. (2013). Finally, the medium which is formed in experiments undergoes very fast expansion and cooling, followed by freeze-out and rescattering, which may also have a strong influence upon critical behavior Herold et al. (2016); Antoniou et al. (2007, 2008); Stephanov (2010). It is thus clear that care should be taken when comparing naive theoretical expectations with actual measurements. Such difficulties call for the use of observables that can serve as robust probes of criticality. The search for reliable signatures demands a correct understanding of both critical fluctuations and relevant background contributions, as well as their dependence on center-of-mass energy, collision centrality and kinematic cuts Sahoo et al. (2013); Fraga et al. (2011); Bzdak and Koch (2012); Ling and Stephanov (2016).

In this paper we introduce a simple, yet reasonably general framework for understanding how both critical and background fluctuations affect different fluctuation measures. For that purpose, we adopt a simplified model of long-range fluctuations Stephanov (2006); Stephanov et al. (1999); Hippert et al. (2016). We also show how kinematic cuts and effects from resonance decay can be included and how they affect our results. Since we restrict ourselves to equilibrium and isotropy assumptions, we focus on multiplicity fluctuations. This simple framework provides a useful tool for understanding multiplicity fluctuations near the critical point.

The outline of the paper is as follows. In Sec. II, we explain our models for critical and spurious fluctuations and discuss the growth of the correlation length near the critical point. In Sec. III, we describe a method for analytically calculating multiplicity fluctuation measures in an equilibrated gas of quasi-free particles. Sections IV and V discuss effects from a limited acceptance window and resonance decay, respectively. Finally, Sec. VI summarizes our results, while Sec. VII presents final remarks.

II Fluctuations near the critical point

II.1 Critical fluctuations

Inspired by Refs. Stephanov (2006); Stephanov et al. (1999), we describe critical fluctuations as fluctuations of an order parameter $\sigma(x)$, with a probability distribution

[TABLE]

and, for small fluctuations, we approximate the effective free energy $\Omega[\sigma]$ as

[TABLE]

where we have imposed $\langle\sigma\rangle=0$. The mass $m_{\sigma}$ can be shown to be the inverse of the correlation length $\xi$ and vanishes at the critical point, rendering long-wavelength fluctuations essentially costless.

Thus, as the critical point is approached, long-range fluctuations get increasingly strong until, for very small values of $m_{\sigma}=\xi^{-1}$, they get too large and Eq. (2) loses validity altogether. For the sake of simplicity, since long wavelengths dominate, we look exclusively at fluctuations of $\sigma_{0}:=\int d^{3}x\,\sigma(x)/V$. Integrating over shorter wavelengths in (1) allows us to define a probability distribution for $\sigma_{0}$,

[TABLE]

Non-Gaussian fluctuations, which are important probes in the search for the critical point, can be described by including in (2) terms of higher order in $\sigma$, which would also extend the validity of the model Stephanov (2006).

Equation (3) describes global fluctuations of the order parameter in equilibrium in some region of the phase diagram around the critical end point where fluctuations are not too large and a classical treatment for the order parameter is valid. It neglects dynamical effects which might make global correlations unattainable in the relevant time scales and provides only a crude model. However, this simplification enables us to make both simulations and analytical calculations in a relatively simple fashion, even in the case of non-Gaussian fluctuations, providing an optimistic, yet valuable tool for understanding how different effects affect critical fluctuations. The off-equilibrium evolution of the order parameter is tackled, for instance, in Refs. Berdnikov and Rajagopal (2000); Mukherjee et al. (2015, 2016); Herold et al. (2016).

A connection with observables requires coupling the order parameter to particle fields. We consider linear couplings in $\sigma_{0}$ which enter the Lagrangian as mass corrections for both protons and pions,

[TABLE]

where the value of $G$ can be estimated as $\sim 300$ MeV near the critical point by symmetry arguments Stephanov et al. (1999), but only the vacuum value of $g\sim 10$ is known, from the linear sigma model Gell-Mann and Levy (1960).

In the absence of reliable information on its value as a function of collision energy or temperature and chemical potential, we treat $\xi=m_{\sigma}^{-1}$ as an input to our calculations. This correlation length, whose size determines the strength of the fluctuations Stephanov (2006); Stephanov et al. (1999), is not determined solely by an equilibrium equation of state or static universality class exponents, but is also influenced by both finite-size effects and dynamics, with the latter being expected to dominate Berdnikov and Rajagopal (2000). We consider limitations to the growth of $\xi$ from critical slowing down as discussed in Ref. Berdnikov and Rajagopal (2000), with the evolution of $\xi(t)$ in the proper-time $t$ given by the ansatz equation

[TABLE]

where $\xi_{eq}(t)=\xi_{0}\;|{t}/{\tau}|^{-\nu/\beta\delta}$ and $\alpha=0.11$, $\nu=0.63$, $z=2+\alpha/\nu$, $\beta=0.326$ and $\delta=4.80$, coming from universality class arguments Zinn-Justin ; Guida and Zinn-Justin (1997); Hohenberg and Halperin (1977). The nonuniversal dimensionless parameter $A$ cannot be directly estimated, but is limited by causality, depending on the combination $x=\tau/\xi_{0}$, where $\tau$ defines a cooling time scale and $\xi_{0}$ is the typical value of the correlation length at which universal behavior sets in Hippert et al. (2016). However, we believe our previous estimate for the time spent in the critical region, $\tau=5.5$ fm, to be overly optimistic, comparable to the lifetime of the system Hippert et al. (2016). This parameter strongly influences the maximum value of $\xi/\xi_{0}$, which depends on the combination $A\,x\leq 1.3\,x^{1.33}$. Taking the more conservative, but still optimistic value of $\tau=1$ fm yields, for $\xi_{0}=1.6$ fm, a maximum value of $\xi/\xi_{0}=1.3$, instead of the previous maximum of $\xi/\xi_{0}=1.9$.

Apart from the different value of $\tau$, the treatment described here is the same as in Ref. Hippert et al. (2016), where some of its aspects are discussed in more detail and the distribution (3) is used to build a simple Monte Carlo implementation of critical fluctuations to which background contributions can be superposed.

II.2 Background fluctuations

One of the most challenging limitations in the search for fluctuations which might serve as a probe of criticality is understanding the role and behavior of background contributions. These contributions may come from any source of undesired statistical correlation, including quantum statistics, fluctuations of imperfectly controlled thermodynamic parameters, such as temperature and system size, and correlations from initial conditions Stephanov et al. (1999); Hippert et al. (2016); Luo and Xu (2017); Gorenstein (2015); Sahoo et al. (2013); Herold et al. (2016); Luo et al. (2013). An experimental discovery of the critical point can hardly be claimed before background is either controlled, well described in terms of beam energy dependence or eliminated by a good choice of signatures.

We use the same model for spurious fluctuations as in Ref. Hippert et al. (2016). For temperature fluctuations we simply take a Gaussian distribution of $5\%$ width. For volume fluctuations, we take a probability distribution $\mathcal{P}_{b}(b)\propto b$ for the impact parameter $b$, and suppose the final volume $V$ to be proportional to the initial overlap area $\mathcal{A}$ between the two nuclear disks of radius $R_{N}$, with a proportionality constant $C$,

[TABLE]

in which we take $R_{N}=r_{0}=6.38$ fm for gold nuclei (see Refs. Vries et al. (1987); Jager et al. (1974)) and $C$ can be fixed by the average radius, or average volume, of the system at freeze-out. This model enables us to find the system-volume distribution in a given centrality range, while disregarding initial-state fluctuations and other effects Luo et al. (2013).

These are very crude models for background fluctuations, but more complete models can be very straightforwardly included within the present framework. Volume fluctuations have previously been considered in Refs. Luo and Xu (2017); Skokov et al. (2013); Gorenstein and Gazdzicki (2011); Luo et al. (2013), as well as the methods for partially correcting them in data analysis Luo and Xu (2017); Luo et al. (2013). Certain fluctuation measures, called strongly intensive, are, by construction, insensitive to such fluctuations Gorenstein and Gazdzicki (2011); Gorenstein (2015).

III Computing fluctuations from energy-level shifts

III.1 Methodology

Besides simulating fluctuations using Monte Carlo techniques Hippert et al. (2016), it is possible to derive analytical expressions for signatures of the critical point within this scenario. We now develop a scheme to calculate observables in the presence of spurious and critical fluctuations. It has the advantage of naturally including finite-size effects.

Inspired by Sec. II.1, we note that a fluctuation of the order parameter $\sigma_{0}$, produces a shift in each single-particle energy level $\omega$,

[TABLE]

where $\delta m=g\;\delta\sigma$ is the corresponding mass correction and $\omega_{0}$ is the original one-particle energy.

Spurious fluctuations of parameters such as volume, temperature and chemical potential can also be described as energy-level shifts. Suppose we have a gas of free particles confined to a typical length scale $R$. For both cubic and spherically symmetric boundary conditions, the allowed momentum levels will be distributed depending on the value of $R$, $p_{i}={\alpha_{i}}/{R}$, where $\alpha_{i}$ represents constants which depend on the boundary conditions. A small change $\delta R$ in the size of the system, hence, modifies the momenta, $p_{i}$, and the one-particle energies, $\omega_{i}$, of the allowed modes,

[TABLE]

Temperature and chemical potential fluctuations can be included in a similar fashion. All the quantities we calculate in this work depend on temperature, chemical potential and energy levels only through the Boltzmann factor $e^{-(\omega_{i}-\mu)/T}$. As a consequence, fluctuations of both temperature and chemical potential are equivalent to fluctuations of the energy levels. For a fluctuation $\delta T$ in the temperature, for instance,

[TABLE]

so it is equivalent to a fluctuation in $\omega_{i}$ of the form

[TABLE]

For a fluctuation $\delta\mu$ in the chemical potential,

[TABLE]

with $\delta\omega^{\mu}_{i}=-\delta\mu$.

In order to calculate the effect of these fluctuating energy shifts, an ensemble of different $(\delta\sigma,\delta R,\delta T,\delta\mu,\cdots)$ and occupation numbers, $\{n_{i}\}$, can be considered. Furthermore, configurations in this ensemble can be organized according to the values of $\delta R$ and $\delta\sigma$, so that averages can be calculated by first averaging over a regular grand-canonical ensemble for a fixed $(\delta\sigma,\delta R,\delta T,\delta\mu,\cdots)$ (denoted by $\langle\cdots\rangle$) and then averaging over $(\delta\sigma,\delta R,\delta T,\delta\mu,\cdots)$ (denoted by $\overline{\cdots}$).

We proceed by expanding the first averages, $\langle\cdots\rangle$, as power series in $\delta\sigma,\delta R,\cdots$ before averaging over their different values. Defining $\overline{\delta\sigma}=\overline{\delta R}=\overline{\delta T}=\overline{\delta\mu}=0$, an expansion to at least order $(\delta\omega)^{2}$ in the single-particle energy corrections is necessary for consistency.

III.2 Second-order expansion

We now turn to the actual calculations. The lowest-order contributions to observables come from the variance of each fluctuating quantity, $\overline{\delta\sigma^{2}}$, $\overline{\delta R^{2}}$ and so on, corresponding to a Gaussian approximation. Before we go further, it is useful to introduce some relations, which can be derived from $\partial^{n}(\beta F)/\partial\zeta_{i}^{n}=\langle(\Delta n_{i})^{n}\rangle$:

[TABLE]

where $\Delta n_{i}:=n_{i}-\langle n_{i}\rangle$, $F$ is the free-energy, $\beta:=1/T$, $\zeta_{i}:=-\beta(\omega_{i}-\mu)$ and the upper (lower) sign stands for bosons (fermions).

We start by the average occupation number of the $i$th energy level. Expanding the ordinary expression for the grand-canonical ensemble up to $\delta\omega_{i}^{2}$, for noninteracting particles,

[TABLE]

and taking the average over fluctuations of the induced shifts,

[TABLE]

where $\zeta_{0\,i}:=-\beta(\omega_{0\,i}-\mu_{0})$.

A calculation of the microscopic correlator111When in an overlined expression, $\Delta n_{i}:=n_{i}-\overline{\langle n_{i}\rangle}$. Stephanov et al. (1999) is also possible,

[TABLE]

From Eq. (15), since $\overline{\delta\omega_{i}}\,\overline{\delta\omega_{j}}\sim{\cal O}(\delta\sigma^{4},\delta R^{4},\cdots)$,

[TABLE]

where derivatives of $f(\zeta_{i})$ can be extracted from Eq. (13) and we have used that $\langle\Delta n_{i}\Delta n_{j}\rangle=\delta_{ij}\,\langle(\Delta n_{i})^{2}\rangle$ for free particles.

Finally, from Eqs. (7), (9), (11) and (12), neglecting terms of order $(\delta\sigma\,\delta R)^{2}$, for instance, and assuming independent fluctuations,222It is, of course, very straightforward to relax this assumption, provided the relevant correlations.

[TABLE]

[TABLE]

The average and variance of the particle multiplicity can then be written as

[TABLE]

[TABLE]

The results presented above can be generalized to other observables and higher-order moments, by including higher orders of $\delta\omega$. Since volume fluctuations modify the momentum modes, the calculation of fluctuation measures involving the momenta requires the computation of terms such as $\overline{p_{i}\,\langle\Delta n_{i}\,\Delta n_{j}\rangle}$ and $\overline{p_{i}\,p_{j}\,\langle\Delta n_{i}\,\Delta n_{j}\rangle}$, for instance.

III.3 Hard-sphere boundary conditions

The finite size of the system requires the specification of boundary conditions in order to use results from Sec. III.2. We choose to work with Dirichlet boundary conditions on a sphere of radius $R$,

[TABLE]

Although they are not very realistic, these boundary conditions are more natural than cubic ones and yet fairly simple.

The eigenstates of the squared momentum $p^{2}$ in a spherical geometry are given by the spherical Bessel functions of the first kind $j_{\ell}(p\,r)$, where the orbital angular momentum $\ell$ is a positive integer. The allowed eigenvalues of $p$ are given by $p^{(\ell)}_{i}=\alpha_{i}^{(\ell)}/R$, where $\alpha^{(\ell)}_{i}$ is the $i$th root of $j_{\ell}(x)$, $j_{\ell}(\alpha^{(\ell)}_{i})=0$. For each pair $(\ell,i)$ there are $2\,\ell+1$ linearly independent one-particle quantum states, corresponding to eigenvalues of the $z$ component of the angular momentum $-\ell\leq m\leq\ell$.

In the absence of kinematic cuts,

[TABLE]

where, from Eq. (17),

[TABLE]

[TABLE]

It is also possible to calculate moments involving different kinds of particles, $A$ and $B$, such as $\overline{\langle\Delta N^{A}\Delta N^{B}\rangle}$, for instance. In this case, the previous calculations apply with small modifications and we have, in Eq. (23),

[TABLE]

[TABLE]

where Eq. (19) can be used with adaptations for the chemical potentials, masses and couplings.

IV Acceptance cuts

Expression (23) is valid for the total multiplicity of particles under the chosen boundary conditions. Experimentally, however, only the particles produced in a certain range of transverse momentum and pseudorapidity are detected and used in data analysis. The behavior of fluctuation measures under variations of these acceptance cuts has been discussed, for example, in Refs. Ling and Stephanov (2016); Karsch et al. (2016); Bzdak and Koch (2012); Garg et al. (2013); Pruneau et al. (2002). Here we show how they can be included in our framework.

The quantum numbers $\ell$ and $p$, chosen in Sec. III.3, give us no information about the separate longitudinal and transverse parts of the momentum. However, it is possible to calculate the fraction of particles of a given momentum $p$ which will be in the selected window of rapidity and transverse momentum, assuming their momenta are distributed isotropically and without angular correlation. Because of quantum statistics, this assumption is not obviously true for indistinguishable particles, but should hold as a reasonable approximation.

The upper and lower transverse momentum cuts, $p_{T}\leq p_{+}$ and $p_{T}\geq p_{-}$, may be written as

[TABLE]

where $u:=\cos\theta$, $p_{T}:=p\,|\sin\theta|$ and $\theta$ is the angle between the momentum of the particle and the beam axis.

We also select a window for the pseudorapidity $\eta=\ln[(p+p_{z})/(p-p_{z})]/2$. Since $\eta=\ln[(1+u)/(1-u)]/2$, $|\eta|\leq\eta_{C}$ yields

[TABLE]

resulting in

[TABLE]

[TABLE]

Combining both rapidity and transverse momentum cuts, we have

[TABLE]

Since we assume an uncorrelated and uniform angular distribution of particles and $|u|$ is uniformly distributed in the interval $(0,1]$, the resulting fraction of accepted particles is

[TABLE]

where

[TABLE]

A plot of $F(p)$ for different acceptance windows is shown in Fig. 1.

Event-by-event cumulants of observables may now be calculated. Attention should be drawn to the fact that the acceptance of particles is probabilistic and follows a binomial distribution, which also contributes to fluctuations. Hence, denoting as $\langle\cdots\rangle_{\textrm{acc}}$ an average over accepted particles only, for a single value of $p$,

[TABLE]

[TABLE]

with which,

[TABLE]

where we denote $\Delta n_{p}:=n_{p}-\langle n_{p}\rangle_{\textrm{acc}}$ or $\Delta n_{p}:=n_{p}-\overline{\langle n_{p}\rangle}_{\textrm{acc}}$ in the corresponding averages over accepted particles.

For the sake of simplicity, we neglect fluctuations of $F(p)$. For midrapidity, $F(p)$ is only sensitive to $p$ in the limits of the acceptance region, as seen in Fig. 1, and fluctuations of $p$ should only be relevant for the acceptance probability in these narrow regions. In this approximation,

[TABLE]

where $p_{0}$ is the unperturbed value of the momentum mode in question. A similar result is found for the correlator of different momentum modes, while lacking the second term in Eq. (39).

Hence, Eq. (23) is modified to

[TABLE]

where $A_{(\ell_{1},i_{1}),(\ell_{2},i_{2})}$ and $B_{(\ell,i)}$ are given by Eqs. (24) and (25) and $C_{(\ell,i)}$ is found from Eq. (15),

[TABLE]

As in Ref. Bzdak and Koch (2012), we use a probabilistic description for acceptance effects. However, we take into account the momentum dependence in the relevant probabilities. A similar treatment will be applied to resonance decay under acceptance cuts in the next section. Note that a finite detection efficiency $e_{0}(p)$ can be introduced by making $F(p)\rightarrow e_{0}(p)\cdot F(p)$.

The implementation of kinematic cuts in our Monte Carlo simulations is far simpler and can be done by independently sampling the momentum direction for each of the produced particles and applying the cuts, also relying on the assumption of uncorrelated, isotropic emission. It is also straightforward to include the effects of a finite efficiency.

V Resonance decay effects

So far, we have only considered direct particles from a thermal distribution. However, a non-negligible fraction of the final-state particles in a collision come from the decay of unstable particles. In this section, we consider the role of resonance decay in multiplicity fluctuations. The impact of resonance decay in fluctuation measures was previously explored in Refs. Bluhm et al. (2017); Mishra et al. (2016); Nahrgang et al. (2015); Sahoo et al. (2013); Garg et al. (2013); Begun et al. (2006); Jeon and Koch (1999), among others.

We consider particles coming from resonance decay to be affected by spurious fluctuations, but not by critical ones, such that they are expected to dilute signatures of criticality. Moreover, both the decay of unstable particles and the detection of its products in the relevant acceptance window are probabilistic processes, which affects fluctuations of the particle multiplicities.

Resonance decay contributions to particle multiplicity fluctuations can be incorporated to our calculations by using the probability of having each decay product in the acceptance window. These probabilities are computed in Sec. V.1 and later used in Sec. V.2 to calculate the desired contributions. Unlike most previous calculations of resonance decay effects, we apply acceptance cuts to the decay products themselves, not to the resonances. For simplicity, we neglect the widths of the resonances. As an example, we apply our methods to the decay of a $\rho$ meson into two pions.

V.1 Probabilistic treatment

We aim to calculate the distribution of decay products found within the acceptance window. Information about both the average fraction of accepted particles and the corresponding fluctuations is necessary. We limit ourselves to two-particle decays and start with branching ratios of $100\%$. We use the associated phase-space volume as a measure of probability for the kinematic variables.

We denote the momentum of the resonance by ${\bf p}_{in}$. Since we assume spherical symmetry, there is no preferred direction for it. The products of the decay have momenta ${\bf p}_{1}$ and ${\bf p}_{2}$, with orientations given by polar and azimuthal angles of $\theta_{1}$, $\theta_{2}$ and $\phi_{1}$, $\phi_{2}$ , respectively. Imposing energy and momentum conservation, the two-particle phase-space volume differential reads

[TABLE]

where $E_{i}:=\sqrt{p^{2}_{i}+m^{2}_{i}}$ and we integrate the direction of the resonance momentum over the whole sphere.

Changing variables from $\theta_{2},\phi_{2}$ to $\theta^{(1)}_{2},\phi^{(1)}_{2}$, defined with their zenith in ${\bf p}_{1}/p_{1}$,

[TABLE]

with $\cos\theta^{*}:=({p_{in}^{2}-p_{1}^{2}-p_{2}^{2}})/2\,p_{1}p_{2}$. Hence, using $dE_{i}/dp_{i}=p_{i}/E_{i}$,

[TABLE]

If we interpret $d\Phi_{2P}$ as proportional to a probability density, Eq. (44) suggests $\phi^{(1)}_{2}$, $\phi_{1}$, $\cos\theta_{1}$ and $E_{1}$ to be uniformly distributed among their full range of values, while $E_{2}$ and $\theta^{(1)}_{2}$, the angle between ${\bf p}_{1}$ and ${\bf p}_{2}$, are determined by the former. Notice that the conditions that $E_{2}^{*}:=E_{in}-E_{1}$ and $\cos\theta^{*}$ fall within the integration regions for $E_{2}$ and $\cos\theta^{(1)}_{2}$ implicitly limit the possible values of $E_{1}$. The constraint $\cos^{2}\theta^{*}|_{p_{2}=p_{2}^{*}}\leq 1$, for example, yields

[TABLE]

with

[TABLE]

[TABLE]

It is now possible to sample the number of decay products inside the acceptance window for a resonance of momentum $p_{in}$. We sample $\phi^{(1)}_{2}$, $\phi_{1}$, $\cos\theta_{1}$ and $E_{1}$ uniformly and calculate $E_{2}$ and $\cos\theta^{(1)}_{2}$, counting the rate of acceptance within kinematic cuts to find the probability of accepting each particle.333The change of variables $\theta_{2},\phi_{2}\leftrightarrows\theta^{(1)}_{2},\phi^{(1)}_{2}$ can be undone by a simple rotation. Results are exemplified in Fig. 2. Condition (45), together with the form of $F(p)$ (see Fig. 1), is responsible for the acceptance probabilities in Fig. 2 vanishing both above $p_{res}\approx 25$ GeV and bellow $p_{res}\approx 80$ MeV.444We thank V. Koch for pointing out that these probabilities should vanish at high $p_{res}$.

In Monte Carlo simulations, the momenta of the decay products can be sampled according to the above scheme for each single resonance.

V.2 Fluctuations from resonance decays

We now apply the results from Sec. V.1 in the context of particle multiplicity fluctuations. Consider the decay of a resonance $R$ into two distinct particles, $A$ and $B$, with probabilities $P_{A}$, $P_{B}$, $P_{2}$ and $P_{0}$ that only particle $A$, only particle $B$, both particles or neither of them are within the acceptance window, respectively. These probabilities provide all the relevant information for the study of multiplicity fluctuations. Then, we calculate averages involving the number of accepted particles, $n_{l}$, where $l=A,B$. For the decay of a single resonance with momentum $p$, for instance, $n_{l}$ is either [math] or $1$ and555For simplicity of notation, we choose to hide the dependence on $p$ of $P_{l}$, $P_{2}$, $P_{0}$ and kinematic averages.

[TABLE]

[TABLE]

where $m\in\mathbb{N}$ and $\langle\cdots\rangle_{\textrm{kin}}$ denotes an average over the angular distribution of particles within kinematic cuts. For two different decays, $i$ and $j$, of resonances with the same momentum $p$,

[TABLE]

Considering now the independent decays of $n^{R}_{p}$ resonances of momentum $p$, summing over $p$ and taking averages over thermodynamic, critical and spurious fluctuations yields for the contribution from correlations among resonance decay products:

[TABLE]

[TABLE]

[TABLE]

As in Sec. IV, we neglect the influence of fluctuations on the fraction of accepted particles, taking the acceptance probability for the unperturbed resonance momentum.666In this case, however, these probabilities are also affected by changes in the mass of the particles, driven by critical fluctuations. Such effects are present in the full simulations. For the contributions from correlations between resonance decay products and direct particles:

[TABLE]

Contributions from correlations between decays of different resonances can also be calculated in a straightforward fashion.

A similar treatment can be applied to decays into two identical particles, with $P_{A}=P_{B}=P_{1}$, in which case, for a single decay, $n_{l}$ is either [math], $1$ or $2$ and

[TABLE]

[TABLE]

[TABLE]

and Eqs. (52) and (53) are combined into a single result,

[TABLE]

A generalization of the results above to decays into more than two particles should be possible using similar arguments. A branching ratio $b_{r}$ can be easily introduced by making $P_{l}\rightarrow b_{r}\,P_{l}$.

VI Results

Calculations within the above scheme can be directly implemented in a computer code, which enables us to calculate the variances and covariances of different particle multiplicities.

Our aim here is neither a direct comparison with experimental data nor a test of different signatures with realistic background contributions. We rather examine the effect of different kinds of limitations to simple fluctuation measures and compare analytical and simulation results. Our background models and resonance contributions are not meant to be complete, but exemplify how different effects may be implemented and how they affect fluctuations. For that reason, only the decay of $\rho$ mesons into pions is considered.

Since the mass and the coupling of pions are better constrained, we concentrate on these particles. A general discussion on the limitations and advantages of our methods and results is presented in Sec. VII.

VI.1 Signatures from pions

We now turn to signatures coming exclusively from charged pions and examine their dependence on the chiral correlation length $\xi$. As in Ref. Hippert et al. (2016), we use $T=130$ MeV for the temperature of the system and $R_{p}=6.8$ fm for its average radius. We choose to display results for the increase in the ratio between the variance and the average of the multiplicity of charged pions, relative to its value at $\xi_{r}$, as a function of $\xi$:

[TABLE]

where $M_{\pi_{ch}}=\overline{\langle N_{\pi_{ch}}\rangle}$ and $V_{\pi_{ch}}=\overline{\langle(\Delta N_{\pi_{ch}})^{2}\rangle}$ are the mean and variance of the charged pion multiplicity and, unlike in Ref. Hippert et al. (2016), we use $\xi_{r}=0.4$ fm $\sim 1/m_{\sigma}^{\textrm{v}ac}$ as a baseline value for the correlation length. Figure 3 displays results for the effects of critical and background contributions as well as the decay of $\rho$ mesons into two pions. These results suggest that the chosen signature would hardly reach $\sim 1\%$. Figures 4, 5, 6 and 7 show results for the signal at $\xi=5\,\xi_{r}=2.0$ fm, $S_{5}$, as a function of the lower transverse momentum cut $p_{-}$, the pseudorapidity cut $\eta_{C}$, a constant detector efficiency and the branching ratio for the decay of $\rho$ mesons into two pions, respectively.

Figure 4 also presents simulation results, showing good agreement with analytical ones for different transverse momentum cuts. This validates the approximation, made in Sec. IV, that fluctuations of $F(p)$ are negligible. The decrease of $S_{5}$ with an increasing $p_{-}$ is due to a larger concentration of pions at low transverse momentum, which also makes the value of the superior transverse momentum cut $p_{+}$ less relevant.

Attention should be drawn to the fact that the signature in Eq. (59) is also affected by particle number fluctuations which are inherent to a grand-canonical ensemble. These fluctuations alone yield a nearly Poissonian contribution, V_{\pi_{ch}}/M_{\pi_{ch}}\big{|}_{GC}\approx 1, which dominates the variance of the particle multiplicity. Since our calculations only have meaning within a grand-canonical ensemble, this contribution is present in all of our results, even the ones for “pure” critical contributions. This is partly the reason why our signals grow more slowly in $\xi$ than one would naively expect. While subtracting a Poissonian contribution from our signature would result in much stronger relative signals, such signals would correspond to much less impressive values in absolute numbers, signaling to small increases in the particle number fluctuations. On the other hand, calculating the signal with respect to the full background contributions, including the ones which are intrinsic to a grand-canonical ensemble, has the advantage of providing a direct comparison to a number we know is of order $\sim{\cal O}(1)$. The signals obtained after subtracting a Poissonian background are related to our results by an approximately constant ratio, which is different for each case in Fig. 3. Without the additional spurious fluctuations, this ratio is of roughly $\sim 80$, but it drops to $\sim 2.2$ in their presence. If the decay of $\rho$ mesons is included, these values drop to $\sim 62$ and $\sim 1.7$, respectively. This shows that, despite growing more slowly with $\xi$, our naive signature is a lot more robust to other background contributions.

The behavior for different pseudorapidity cuts in Fig. 5 can be understood by arguments similar to the ones in Ref. Ling and Stephanov (2016), although, there, expansion effects are considered and background contributions are not contemplated. Because of the isotropy assumption, as a larger rapidity window is used, a larger number of pairs of momentum modes, roughly proportional to the square of the window in the cosine of the polar angle $[\Delta(\cos\theta)]^{2}$, is correlated by critical contributions. Since the results are normalized by the average multiplicity of charged pions ($\propto\Delta(\cos\theta)$), they scale roughly as $\propto\Delta(\cos\theta)\propto F(p)\propto\tanh\eta_{C}$ if $p_{T}$ is unrestricted, as can also be seen from Eq. (40) and is confirmed by the very successful $a\tanh(\eta_{C})$ fit to the results. Since the expansion of the system was not considered, the interpretation of these results is purely geometric.

The dependence of $S_{5}$ on the detector efficiency, shown in Fig. 6, can be modeled by assuming a detection probability proportional to the efficiency $e$, so that $V_{\pi_{ch}}/M_{\pi_{ch}}=\lambda(\xi)\,e+1$ and $S_{5}$ takes the form

[TABLE]

which is in excellent agreement with results.777Fits to the results yield $\kappa=(2.5107\pm 0.0000003)\,10^{-2},\;\gamma=(1.25693\pm 0.00001)\,10^{-4}$ without background contributions and $\kappa=(2.36801\pm 0.0000004)\,10^{-2},\;\gamma=(8.54053\pm 0.000004)\,10^{-3}$ when taking them into account.

As for the dependence of the results on the branching ratio of the decay of a $\rho$ meson into two pions, shown in Fig. 7, it can be roughly described by assuming a Poissonian distribution on the number of decay products, with a mean proportional to the branching ratio $b_{r}$. Combining contributions from decay products with critical ones while assuming statistical independence between the two and neglecting the dependence of $M_{\pi_{ch}}$ in $\xi$ yields

[TABLE]

which can fit the results quite well.888Fits to the results yield $\chi=2.45388\pm 0.002523,\;\zeta=(2.59846\pm 0.02242)\,10^{-3}$ without background contributions and $\chi=1.28827\pm 0.00122,\;\zeta=(7.27304\pm 0.02838)\,10^{-3}$ when taking them into account.

The signal in Fig. 3 scales quadratically with the correlation length, with a proportionality factor which depends mainly on the square of the coupling $G$ in Eq. (4), as well as our models for spurious contributions and the chosen acceptance window. While the estimate of $G=300$ MeV near the critical point carries large uncertainties related to medium effects and the properties of the $\sigma$ meson, it is clear that it should be much smaller than in vacuum, by a factor of $\sim 1/6$. Since this estimate considers a vacuum sigma mass of $600$ MeV, we believe it to be optimistic Stephanov et al. (1999). Regarding spurious contributions, while our models most likely overestimate the role of temperature fluctuations, they largely neglect the complex nature of heavy-ion collision experiments. Among others, the effect of indirect particles is underestimated, since we have only considered the decay of $\rho$ mesons. Finally, although we rely on the the simplified interaction Lagrangian of Eq. (4), we believe it captures most of the physics, at least for second-order fluctuations. Further caveats to our results are discussed in Sec. VII.

VII Concluding remarks

In this paper, we have outlined a simple framework for calculating multiplicity fluctuations in the neighborhood of the QCD critical point. As we have shown, it allows for the introduction of general long-range fluctuations, acceptance cuts and resonance decay in a systematic way. Acceptance window effects were treated probabilistically, giving rise to extra contributions to fluctuations, and applied to resonance decay. Finite-size effects were also included, which could be important at low collision energies. Moreover, there is no reason one should not extrapolate results to the continuum limit.

While we have calculated simple signatures from Gaussian fluctuations of charged pions, the same treatment might as well be applied to the calculation of more interesting signatures, such as higher moments of the net-proton distribution. Nevertheless, the inclusion of other particles, as well as an extension to non-Gaussian fluctuations, requires the inclusion of poorly known parameters, such as the proton mass $m_{p}$ and the proton coupling constant $g$ near the critical point, and the parameters required for an appropriate description of non-Gaussian fluctuations Stephanov (2009), a problem which may be partially solved by somehow extracting information about fluctuations in equilibrium from a provided equation of state.

Our approach presents some important limitations. It includes no dynamics, except for the estimation of the maximum value of the correlation length, and no freeze-out mechanism or rescattering. Fluctuations of the $\sigma_{0}$ and other parameters are taken to be perfectly homogeneous and particles are assumed to obey perfect symmetry and excessively restrictive boundary conditions.

Moreover, our background fluctuations are still unrefined, although they could be easily substituted for more physically and experimentally motivated models. Information on the statistical distribution of thermodynamic parameters may be extracted from the event-by-event distribution of the mean transverse momentum of particles, specially for temperature fluctuations. In this case, different particle species, with different masses and charges, might be used to separate effects from temperature, chemical potential and volume fluctuations. In the case of volume fluctuations, a more complete model could be achieved by assigning a probability distribution for the proportionality factor $C$ in Sec. II.2.

A relevant restriction is that we assume isotropy and global equilibrium, neglecting the effects of flow. This means that our results for the acceptance window dependence should not be taken at face value, especially at large pseudorapidities (compare with Ling and Stephanov (2016)) and, even so, their validity is limited to central collisions. Although these effects were not considered, they can be roughly included by using a prescribed velocity profile and writing

[TABLE]

where $\tilde{\Omega}_{\textrm{acc}}(p,{\bf x})$ is the solid angle coverage of the acceptance window when boosted with the velocity of the fluid element at $\bf x$. This would, however, provide only a crude solution.

The methods we have employed seem to be very adaptable, enabling several effects to be considered at once in reasonably simple analytical calculations and simulations which are relatively cheap in computer power. These simulations produce integer numbers of particles without introducing artificial Poissonian noise and might be used to feed some kind of rescattering algorithm.

Another nice feature of our simulations is that they can be adapted to include global conservation laws (see, for instance Refs. Bzdak et al. (2013); Begun et al. (2006, 2004)) by restricting them to configurations with fixed values for conserved quantities, discarding samples which do not follow this criterion. However, this would raise the computing time by increasing both the number of necessary samples and the number of considered momentum modes.

Acknowledgements.

The authors would like to thank V. Koch, M. Lisa, L. Moriconi, D. Parganlija, P. Sorensen, M. Tannenbaum and G. Torrieri for elucidative conversations and constructive criticism, and the organizers and participants of the INT Program “Exploring the QCD Phase Diagram through Energy Scans” for two weeks of very inspiring discussions and presentations. We thank the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work. We also acknowledge CNPq and FAPERJ for their financial support.