From Harmonic Oscillators to Multiplicity Distributions in Particle Production Processes

Introduction

Many particles are produced in a high-energy collision event, and the total number of produced particles fluctuates from one event to another. The number of produced particles form a distribution $P(N)$ that (upon proper normalization) represents the probability that a given collision (event) produces exactly $N$ particles. For particles emitted from independent and identically distributed sources, the Poisson distribution provides a baseline for multiplicity fluctuations. However, It has been observed that the multiplicity distribution $P(N)$ obtained from pp, pA, and AA collisions are represented by distributions similar to the Negative Binomial Distribution (NBD) [1] that is connected with significantly larger fluctuations in particle multiplicity than that showed by the Poisson distribution.

In this article, we study the multiplicity distribution from the information theory approach [2, 3] using harmonic oscillators whose energies $\omega_i$ may correspond to the central values of momentum bins representing particle momenta. Hence, populating each bin amounts to changing the harmonic oscillator quantum numbers $n_i$, and the total number of particles produced in an event is given by $\sum_i n_i = N$.

With the help of the probability $P(N)$, we define an entropy that we maximize subject to a minimal (and physically relevant) set of constraints and find out the equilibrium set of probabilities $P(N)$. Using this approach for our oscillator model, we study the origin of the ‘Poissonian and beyond-Poissonian’ distributions. First, we study a single-mode harmonic oscillator as the building block of our model, a perspective also utilized recently in Ref. [4]. Eventually we generalize our results to multi-mode oscillators. The harmonic oscillator approach offers analytical tractability while capturing the essential physical scenario.

Eqs. (8), (21), and (24) highlight the main findings of the paper.

• Starting from the Shannon entropy, we show that the equilibrium probability distribution $P(N)$ for a system of degenerate harmonic oscillators is represented by the NBD (Eq. (8)). Also, for infinitely many such oscillators, this distribution approaches the Poisson distribution when the mean multiplicity is kept fixed.

• The NBD derived from the Boltzmann-Gibbs statistics assumes factorization of probabilities of micro-states, and we generalize this scenario by introducing the Tsallis entropy. We derive a generalized negative binomial distribution (gNBD) (Eq. (21)) considering degenerate harmonic oscillators following Tsallis statistics. We observe that, for the chosen datasets, the gNBD describes particle multiplicity data at high multiplicity region where a single NBD starts deviating from experimental data [3]. This indicates that the gNBD may provide an alternative to the multi-component distributions used in the literature [5, 6, 7, 8, 9]. However, more comprehensive analysis is required to ascertain this conclusion.

• The generating function of the probability distribution is also a quantity of interest. While studying the generating function of the multi-mode non-degenerate oscillators within the Tsallis statistics framework, we encounter

a special function that we call the ‘Lerch-type generalization’ of the Barnes zeta function (Eq. (24)). Some special cases of this function have been studied in physics [10] and mathematics [11]. However, to the best of our knowledge, the generalized Lerch Barnes zeta function reported in this article has not been considered earlier in the context of nonadditive statistics and multiplicity generating studies. It would be interesting to study the analytical structure of this function.

The rest of the paper is devoted to a detailed description of the model and their implications. In the following section, we consider the Boltzmann-Gibbs (BG) statistics and Shannon entropy. Later, we generalize the calculations by involving the Tsallis statistics (section III) that introduces correlation among the oscillator modes. We detail phenomenological implications in section IV, and we summarize and conclude in § 7.

Boltzmann-Gibbs Statistics, Oscillators, and Multiplicity Distributions

A single-mode oscillator

Below we present the results that are already known by studying a single-mode oscillator in Boltzmann-Gibbs statistics. By extremizing the Shannon entropy (see Appendix § 8 A), subject to probability normalization and fixed mean energy constraints (expressed in terms of the energy of a micro-state $\epsilon$), the equilibrium set of probabilities is given by

For a single-mode harmonic oscillator, probability for $N$-th state with (normal ordered) energy $\epsilon= \omega N$ is

where the partition function is given by

Hence, the probability distribution in Eq. (2)

is the geometric distribution with mean

and variance $\operatorname{Var}(N) = \langle N\rangle(1+\langle N\rangle)$.

Multi-mode oscillators, NBD, and Poisson

Probability distribution

For a collection of $k$ harmonic oscillators of different frequencies $\omega_1$, $\omega_2$, .......$\omega_k$ (with indices $n_1$, $n_2$, .......$n_k$) in the Boltzmann-Gibbs statistics, the set of probabilities $P(\{n_i\})$ is given by

where $\mathcal{Z}$ is the partition function. From the above equation, the multiplicity distribution $P(N)$ can be equated to

Considering degenerate oscillators having $\omega_i\approx\omega$ and expressing $p\equiv\exp(-\beta\omega)$, the distribution becomes as follows

Hence, we observe that our model involving degenerate oscillators in the BG statistics gives rise to the negative binomial distribution that describes multiplicity data obtained in many experiments.

Generating function

For the sake of completeness, we also study the generating function given below

where $z$ is a dummy variable. In the BG oscillator model

For degenerate oscillators we obtain

The function $G(z)$ is the generating function for the NBD and is related to $P(N)$ by

It is also worth mentioning that when $k\to\infty$, the NBD approaches the Poisson distribution when the mean multiplicity is kept fixed.

Nonadditive statistics, oscillators, and multiplicity distributions

In the BG statistics, we assume probability micro-states to be independent. However, for systems produced in high-energy collisions, fluctuations, long-range correlations, and memory effects may be prevalent, and the assumption of factorization may be relaxed. Hence, in this section we explore the statistical mechanics of oscillators in a generalized setting inspired by nonadditive statistics.

A single-mode oscillator

By extremizing the Tsallis entropy, subject to probability normalization and fixed mean energy constraints (expressed in terms of the energy of a micro-state $\epsilon$), the equilibrium set of probabilities is given by

where (for $q > 0$, and $x > 0$)

The normalization factor is given by the Hurwitz zeta function $\zeta$

where $\zeta(s,a)=\sum_{n=0}^{\infty}(n+a)^{-s}$. In stead of a geometric distribution (4), now we obtain the Zipf-Mandelbrot distribution [12]

The mean multiplicity is given by

For several more analytical results involving two and three oscillators and nonadditive statistics, please refer to Appendix.

Multimode oscillators

Probability distribution

Following (§ 2), we now study multiplicity distributions using multimode oscillators in nonadditive statistics. The set of nonadditive probabilities are given by

For a fixed $N=\sum_i n_i$, the distribution $P(N)$ is given by the following constrained summation

The partition function in the denominator can be expressed in terms of the Barnes zeta function $\zeta_B$ [10] as follows

In the above equation, $\omega_{\mathrm{ref}}$ is a reference energy, $\{a_i\} \equiv(a_1, a_2, a_3, \ldots) \equiv\{\omega_i/\omega_{\mathrm{ref}}\}$, and $\lambda= \{(q-1)\beta\omega_{\mathrm{ref}}\}^{-1}$. If we assume degenerate oscillators such that $\omega_i \approx\omega$ $\forall$ $i$, and $\sum_i n_i=N$, Eq. (19) becomes

The partition function $Z_q$ can be evaluated as

Generating function

Since the generating function $\mathcal{G}(z)$ admits a sum of $N$ ranging from zero to infinity, we obtain the following

The summation in the above equation can be thought of as a Lerch-type generalization [13] of the Barnes zeta function that we can denote as $\Phi_B(z,1/(q-1),\lambda|\{a_i\})$. The generating function for non-degenerate oscillators is obtained as follows

To the best of our knowledge, combining a Barnes structure with a Lerch-like generalization has not been explicitly considered in the context of nonadditive statistics. Some special cases have been considered in physics ($z=1$ case in Ref. [10]) and mathematics ($\{a_i\}=\{1\}$ case in Ref. [11]). It would be interesting to explore the analytical structure of such a generalization and other physical scenarios where it may arise.

For degenerate oscillators, the generating function in Eq. (23) becomes

The above formula can be expressed in terms of the Lerch’s $\Phi$ function [13] (see Appendix § 8 D)

In the limit $q \to1$, the above equation converges to Eq. (11)

The Oscillator Framework and Phenomenology

In this section, we compare the results obtained using the oscillator framework with the observed multiplicity distribution data. We find that the NBD obtained from the oscillator picture describes multiplicity distribution data up to certain multiplicity values (Figure 1), and starts deviating from the experimental observations for higher $N$ values. A similar observation has been reported in the literature. For example, Refs. [5, 6] utilize a phenomenological two-component NBD for describing this deviation. However, the present work utilizes the NBD and a generalized NBD (Eqs. (8) and (21)) derived from an information theoretic model to describe experimental data. Comparison of Eqs. (8) and (21), referred to as NBD and gNBD, with observed multiplicity distributions has been made in Figure 1-Figure 4. The figures demonstrate that the while both NBD and gNBD describe low multiplicity region satisfactorily, they deviate at higher multiplicity region whereas the gNBD describes the data points. To demonstrate this deviation, we have chosen two rapidity windows, $|\eta|<0.5$ and $|\eta|<2$, at two different pp collision energies 900 GeV and 7 TeV. It is apparent that choosing a higher rapidity window helps explore higher multiplicity region where the effect of the generalized NBD derived in the manuscript is prominent. It is also noteworthy that comparison with experimental data leads to one of the fewer instances in high energy physics where the entropic parameter $q<1$.

Summary, Conclusion, and Outlook

In this work, we have proposed a harmonic oscillator-based model to describe multiplicity distributions in particle production processes. The probability distributions are derived from maximization of Boltzmann-Gibbs and generalized nonadditive entropies. We observe that while both Boltzmann-Gibbs statistics and nonadditive statistics provide good description of chosen datasets at low multiplicity region, generalized NBD derived from nonadditive statistics describe data more efficiently in higher multiplicity region.

The mathematical construction involving a system of degenerate oscillators provides a connection with the phenomenological NBD used in the literature. The phenomenological NBD is expressed in terms of the mean multiplicity $\langle N\rangle$ and the $k$ parameter. Since the total number $N$ is equal to $\sum_i n_i$ in the oscillator picture, the mean multiplicity $\langle N\rangle$ is equal to $\langle n_1\rangle+\langle n_2\rangle+\cdots\langle n_k\rangle$. Each term $\langle n_i\rangle$ represents the Bose-Einstein single particle distribution (since we consider a bosonic oscillator) and is expressible in terms of the quantity $p=e^{-\beta\omega}$ (defined above Eq. (8)). For a system of $k$ degenerate oscillators, $\langle N\rangle=kp/(1-p)$, which implies that $\langle N\rangle/(\langle N\rangle+k)=p$. Hence, the phenomenological $k$ parameter, that signifies dispersion, may correspond to the effective number of degenerate oscillators. Also, the factor $\langle N\rangle/(\langle N\rangle+k)$ can be explained as the Boltzmann-Gibbs weight of adding one quantum to an oscillator, where $\langle N\rangle/k$ is the single particle distribution.

Analysis of experimental data using the oscillator framework reveals that the phenomenological parameter $q<1$ (which, when subtracted from 2, gives the entropic parameter $q$). A possible explanation for the sub-unity values of $q$ may be provided from the finite reservoir picture [15]. If we consider the observed system to be in contact with the unobserved part that may act as a reservoir, a nonadditive power law distribution may arise when the effective number of degrees of freedom of the bath is finite from the perspective of the system. For wider acceptance intervals (e.g. wider rapidity window), finite reservoir effects may be more prominent, potentially explaining the observations made in the paper. However, $q<1$ values introduce an upper cut-off in $N$ for real-valued probabilities, as we can see from Eq. (21).

The oscillator picture has been utilized earlier to describe transverse momentum spectra [4], and starting from the same framework, the present article describes multiplicity distributions. This framework can also be extended to explore other observed quantities such as inelasticity, higher order multiplicity cumulants, and forward-backward correlation. We reserve these studies for the future.

Appendix

Maximum entropy approach

We maximize the Shannon entropy

subject to some constraints. Among all the discrete distributions with a specified mean value $\langle n\rangle$, the maximum entropy distribution is the geometric distribution. Probability distribution $P(n)$ with mean $\langle n\rangle$ corresponds to maximum entropy when the quantity (expressed in terms of the multipliers $\lambda_1$ and $\lambda_2$)

satisfies $\partial F/\partial P(n)=0$ for all $n$. The derivative yields $-1-\ln P(n)=\lambda_1+\lambda_2n$ implying $P(n)=\exp(\lambda_1-1)\exp(\lambda_2n)$. Setting $r=\exp(\lambda_2)$, we obtain $P(n)=Cr^n$, where $C=\exp(\lambda_1-1)$. The normalization constraint $\sum_n P(n)=1$ implies $C=1-r$, and

Similarly, the Zipf-Mandelbrot distribution maximizes the Shannon entropy under the constraint [16]

Following the standard procedure, we consider the quantity

From $\partial F/\partial P(n)=0$ for each $n$, and using the substitutions $\lambda_1-1=\ln\lambda$ and $\lambda_2=1/(1-q)$, we obtain

The normalization constraint implies $\lambda^{-1}=Z=[(q-1)\beta\omega]^{\frac{1}{1-q}}\zeta\left(\frac{1}{q-1},\frac{1}{(q-1)\beta\omega}\right)$.

For Tsallis entropy

under constraint $\sum_n nP(n)=\langle n\rangle$, we consider the quantity

Setting $\partial F/\partial P(n)=0$ for each $n$ leads to

Substituting $q\to2-q$ and reorganizing equation (36) we have

where $\lambda=\beta\omega$ for $\lambda_1=\frac{1-(2-q)Z^{q-1}}{q-1}$ and $\lambda_2=\beta\omega(q-2)Z^{q-1}$.

Two oscillator model: nonadditive statistics

We consider two oscillators representing two bosonic particles and impose the constraint that the total number of particles is $N$. Hence

where $\zeta$ represents the Hurwitz zeta function. The nonadditive partition function $Z_q$ can be evaluated in terms of the two-parameter Barnes zeta function $\zeta_{\mathrm{B}}^{(2)}$ as follows.

In the above equation, $\omega_{\mathrm{ref}}$ is a reference energy, $\{a_i\}\equiv(a_1,a_2)\equiv(\omega_1/\omega_{\mathrm{ref}},\ \omega_2/\omega_{\mathrm{ref}})$, and $\lambda=\{(q-1)\beta\omega_{\mathrm{ref}}\}^{-1}$.

Three oscillator model: nonadditive statistics

We present a three oscillator model comprising two bosonic and one fermionic oscillator. Hence $P(\{n_i\})$ can be written as:

If we consider the index $n_3$ to be associated with a fermionic oscillator it assumes only two values 0 and 1. Hence, if we impose the constraint that $n_1+n_2+n_3=N$, we have either $n_2=N-n_1$ (when $n_3=0$) or $n_2=N-n_1-1$ (when $n_3=1). Hence, $P(N)$ can be written as a summation of two terms,

Following the previous section, the first sum can be evaluated as:

The second sum can be evaluated as:

$Z_q^{(1)}$ and $Z_q^{(2)}$ can be evaluated following the method outlined in the previous section.

Generating functions: degenerate oscillators in nonadditive statistics

Let us evaluate the quantity in the parentheses of Eq. (25). When $j=0$, the quantity is proportional to the Lerch’s $\Phi$ function

When $j=1$

When $j=2$

When $j=3$

By induction,

Hence the generating function corresponding to the multiplicity distribution $P(N)$, considering that multiplicity $N$ is distributed among $k$ degenerate oscillators in a nonadditive thermal bath, is given by

For the numerical purposes, expressing derivatives of the Lerch’s $\Phi$ in terms of the shifted $\Phi$ functions will be more useful. Let us re-evaluate the parentheses of Eq. (25) for the $j=1$ case.

Hence the generating function $\mathcal{G}_{\mathrm{D}}^{(2)}(z)$ for the two oscillator ($k=2$) case is

Since for the $k=2$ case $g(N)=N+1\Rightarrow a_0=a_1=1$

Acknowledgements

TB acknowledges funding from the European Union’s Horizon Europe Programme, via the ERA Fellowship Grant Agreement number 101130816.