# Multiparticle Bose-Einstein Correlations

###### Abstract

Multiparticle symmetrization effects are contributions to the spectra of Bose-symmetrized states which are not the product of pairwise correlations. Usually they are neglected in particle interferometric calculations which aim at determining the geometry of the boson emitting source from the measured momentum distributions. Based on a method introduced by Zajc and Pratt, we give a calculation of all multiparticle symmetrization effects to the one- and two-particle momentum spectra for a Gaussian phase space distribution of emission points. Our starting point is an ensemble of -particle Bose-symmetrized wavefunctions with specified phase space localization. In scenarios typical for relativistic heavy ion collisions, multiparticle effects steepen the slope of the one-particle spectrum for realistic particle phase space densities by up to 20 MeV, and they broaden the relative momentum dependence of the two-particle correlations. We discuss these modifications and their consequences in quantitative detail. Also, we explain how multiparticle effects modify the normalization of the two-particle correlator. The resulting normalization conserves event probabilities, which is not the case for the commonly used pair approximation. Finally, we propose a new method of calculating Bose-Einstein weights from the output of event generators, taking multiparticle correlations into account.

###### pacs:

PACS numbers: 25.75.+r, 07.60.ly, 52.60.+h^{†}

^{†}preprint: CU-TP-867

## I Introduction

Most hadrons are emitted in the final stage of a relativistic heavy ion collision. They do not probe directly the hot and dense intermediate stages where quarks and gluons are expected to be the relevant physical degrees of freedom for equilibration processes. The geometrical size and dynamical state of the hadronic phase space emission region, however, depends sensitively on the entire evolution of the collision. This motivates current attempts to reconstruct its spatial and dynamical state from the experimental hadron spectra and to use it as a starting point for a dynamical back extrapolation into the hot and dense intermediate stages [1]. Two-particle correlations of identical particles which are sensitive to the space-time characteristics of the collision [2], play a crucial role in this approach. The reconstruction program based on their analysis has very good prospects: due to the increasing event multiplicities and larger statistics of the CERN SPS lead beam program (and the yet better quality data expected from Relativistic Heavy Ion Collider RHIC at BNL), particle interferometric measurements start showing statistical errors on the few percent level. Also, systematic errors are increasingly better understood. Theoretical calculations should aim for a similar accuracy and control necessary approximations quantitatively.

One uncontrolled approximation used so far in almost all particle interferometric calculations is to neglect for the particle momentum spectra of Bose-Einstein symmetrized -particle states all multiparticle correlations which cannot be written in terms of simpler pairwise ones. This reduces the number of terms contributing to the two-particle correlator to a manageable sum over all particle pairs. Beyond this approximation, two approaches have been used in the literature. First, Zajc has employed Monte Carlo techniques [3] to generate events with realistic multiparticle correlations. This amounts to a shifting prescription which modifies the momentum distribution of simulated events according to unit weights (which themselves depend on the space-time structure of the source). Second, Zajc has found [3] the first steps towards a calculational scheme brought into final form by Pratt [4]: for model distributions, the -particle spectra are given by a simple algorithm involving only two types of terms: and . In practice, this reduces the sums over all permutations, typical for the calculation of momentum spectra, to sums over all partitions of .

For both these approaches, there are first numerical calculations [3, 4, 5] and related analytical attempts [6, 7, 8, 9, 10] to control multiparticle effects to the one- and two-particle spectra, but a detailed study of their momentum dependence is missing, even for simple models. This work aims at filling this gap, making quantitative statements about the extent to which the slope of the one-particle spectrum and the relative momentum dependence of the two-particle correlations are modified due to multiparticle symmetrization effects. Our investigation takes the set of phase space emission points as initial condition. For notational simplicity, we restrict the discussion to one particle species, like-charge pions say. To the set , we associate a symmetrized -particle wave function [11]

(1) |

where the sum runs over all permutations of the indices, is a shorthand for the 3-dimensional coordinates , and the functions denote single particle wavepackets centered around at initial time and propagated according to the free time-evolution. Final state interactions, which imply a structure of the -particle state different from (1), will not be considered in the present work. The wave function defines the boson emitting source. What we are interested in is the calculation of the one- and two-particle momentum spectra resulting from (1), the information they contain about the initial distribution of the ‘emission points’ , the extent to which these results modify the predictions of the pair approximation, and finally the algorithm which implements the numerical calculation of multiparticle spectra from the initial distribution .

Our work is organized as follows: Section II shortly sets up and illustrates the general formalism via which particle momentum spectra are calculated from an -particle state. In section III, we discuss the properties of Gaussian wavepackets which we choose for the single-particle states in (1). Section IV contains the main results. It gives a complete qualitative and quantitative analysis of multiparticle contributions to the one- and two-particle spectra for a boson emitting source of Gaussian phase space distribution. In section V, we discuss the implications of this model study for an algorithm which calculates the Hanbury Brown Twiss (HBT) radius parameters and momentum spectra for an arbitrary distribution of emission points. Results and related conceptual questions are summarized and discussed in the Conclusions.

## Ii The formalism

We want to determine for the -particle symmetrized wave function the detection probability for measuring identical bosons at time at the positions and momenta . We calculate first the -particle Wigner phase space density for [12]

(2a) | |||

(2b) |

Here, denotes the momentum operator acting on . The one-particle pseudo-Wigner functions provide the basic building blocks for the calculation of the -particle momentum spectrum which is obtained by integrating (2a) over all spatial coordinates,

(3a) | |||||

(3b) |

where is a normalization constant, , ensuring that the probability of detecting particles is one. For a free time evolution of , the integration over the spatial components of (3a) leads to a time-independent expression , since interactions between the particles are necessary to change the momentum distribution in time. In contrast, integrating over all momenta leads to the detection probability of the bosons at positions , which is a time-dependent quantity since free evolving bosons change their positions in time. The calculation of the -particle momentum spectrum according to (3a) involves a sum over terms. Due to the factorization of , this reduces to a sum over terms,

(4) |

In what follows, we are especially interested in the one- and two-particle momentum spectra , associated with the -particle state . These are obtained by integrating over all but one, respectively two momenta,

(5a) | |||||

(5b) | |||||

(5c) |

All particle momentum spectra are given in terms of the building blocks (which determine ) and . Once the analytical form of the single particle wavefunctions is specified, these are readily calculated. In what follows, capital letters denote measurable position and momentum coordinates, small letters , , denote the centers of wavepackets which are not directly measurable. The only exception to this is the measurable relative momentum of the two-particle correlator which we denote by a small letter.

### ii.1 The Zajc-Pratt algorithm

Dynamical correlations between particles in the source are reflected in correlations in the set of emission points . If there are no correlations, then the initial distribution of the centers of single particle wavepackets is given by a one-particle probability distribution . The -particle spectra for a set of events with multiplicity are obtained by averaging over this distribution

(6a) | |||||

(6b) |

A particular model distribution with correlations which assumes that the emission probability of a boson is increased if there is another emission in its vicinity [10], is e.g. obtained by replacing in (6a)

(7) |

Here is an averaged normalization defined below. The technical advantage of adopting (7) is that the -dependent normalization factor in the spectrum of (6a) is canceled. This allows to write without approximation all spectra in terms of the building blocks [4]

(8a) | |||||

(8b) |

The resulting Zajc-Pratt (ZP) algorithm for the calculation of one- and two-particle spectra reads [3, 4, 10]

(9) | |||

(10) | |||

(11) |

These spectra are normalized to unity. In Appendix A, we give their derivation in some combinatorical detail. Due to the ZP-algorithm, the calculation of the -particle spectra is reduced from sums over all permutations to sums over all partitions of a set of points into subsets of points, . The number of partitions of grows asymptotically like . For explicit calculations with event multiplicities in the hundreds, it is hence important to get control over the -dependence of the Pratt terms and . This is our strategy in section IV.

### ii.2 A simple example: the Zajc model

To illustrate the above formalism, we consider a normalized -particle density matrix for multiparticle states , created by repeated operation of the single-particle creation operator ,

(12a) | |||||

(12b) | |||||

This density matrix specifies in particular the one- and two-particle spectra and . For the Gaussian model distribution [3]

(13a) | |||||

(13b) |

the effects of multiparticle correlations on the HBT-radius parameters have been considered already by Zajc. His discussion however is restricted to an explicit calculation of three-particle symmetrization effects and to qualitative estimates of higher order contributions. Here, we demonstrate that the ZP-algorithm allows for a complete quantitative analysis of the Zajc model. The first step is to identify the building blocks of the particle spectra (5),

(14a) | |||||

(14b) |

The calculation of the terms reduces then to -dimensional Gaussian integrations. One finds and

(15a) | |||||

(15b) |

where denotes Chebysheff polynomials of the first kind [26]. The momentum-dependent terms read

(16a) | |||||

(16b) | |||||

(16c) | |||||

(16d) |

The main message of these involved but explicit expressions is contained in the m-dependence of the terms and . These specify the - and -dependence of the building blocks and hence of all spectra. Here, they are functions of the phase space volume only, and their behaviour can be understood by simple arguments:

The factors generically increase with increasing , see Figure 1. Especially, in a limiting case, one finds

(17) |

The reason is that Bose-Einstein symmetrization effects enhance the low momentum region of the one-particle spectrum, leading to steeper slopes. This one-particle spectrum is a linear superposition of Gaussian terms , which due to show increasingly steeper slopes.

The -dependence of governs the -dependence of the two-particle spectrum. The terms depicted in Figure 1 decrease with increasing and have the limiting values

(18) |

Zajc has concluded on the basis of this behaviour that [3] “the two-particle correlation function becomes a superposition of terms with successively broader distribution in , leading to an increasingly smaller value for the inferred radius.” The reason is that Bose-Einstein symmetrization effects enhance the probability of finding bosons closer together in configuration space and hence result in broader -distributions of the two-particle spectra.

The above arguments explain the effect of multiparticle correlations qualitatively. For a quantitative understanding, the weights of higher order terms contributing to the one- and two-particle spectra are important. These weights are governed by the terms which (up to a correction factor of order unity) essentially decrease like th powers of the inverse phase space volume. For event multiplicities in the hundreds, a quantitative analysis can then be done numerically, using the analytical expressions (15) and (16). We defer such a study to a slightly more general model in section IV where certain analytical simplifications allow for a more transparent discussion. A short comparison of the qualitative and quantitative properties of the models studied here and in section IV is given in the text following Eq. (27).

## Iii Gaussian wavepackets

In the example of section II.2, we have started from single particle creation operators whose momentum support does not carry a label . As a consequence, the -particle states considered in section II.2 are build up from single particle wavefunctions with identical phase space localization. We now adopt a more general setting in which a set of phase space points is associated with Gaussian wavepackets , centered at initial time around the points with spatial width , [13, 14, 11]

(19) |

The Fourier transform of is proportional to and can be compared to the momentum support of in (12a). The corresponding -particle symmetrized states reduce to those considered in section II.2 for and , when the momentum support function becomes independent of the particle label .

The one boson state (19) is optimally localized around in the sense that it saturates the Heisenberg uncertainty relation , with at initial time . Different choices of single particle wavefunctions can lead to different results for the calculated spectra. We have some control over the extent to which different choices matter by finally checking the -dependence of our results. The model parameter can range between . In Refs. [11, 15], it was argued that a realistic value for is for pions of the order of the Compton wavelength.

The free time evolution of is specified by the one-particle hamiltonian which acts as a multiplication operator in momentum space,

(20) | |||||

The resulting building blocks for the -particle spectra are

(21a) | |||||

(21b) |

To streamline our notation, we have neglected in a factor . The product of these factors cancels in the calculation of the Wigner function (2b) and a fortiori in all the functions derived from it. The function denotes the probability that a boson in the state is detected with momentum . This measured momentum has a Gaussian distribution around the central momentum of the wavepacket.

The functions in (5c) characterize the overlap between the wavepackets and and play an important role in the ZP-algorithm. They take a particularly simple form if all particles are emitted in a flash,

(22) | |||||

All terms contributing to or contain the same number of factors and hence, the normalization of does not matter in what follows. For notational convenience, we change it by fixing . The functions measure the distance between the phase space points and . This distance measure depends on the wavepacket width but leaves the phase space volume independent of ,

(23) |

## Iv Multiparticle correlations for a Gaussian Model

We now determine quantitatively multiparticle correlation effects for a source of identical bosons whose wavepackets of spatial width are emitted instantaneously according to a Gaussian -particle phase space distribution (7) with

(24) |

Our main aim is to study for this model the extent to which multiparticle correlations modify the slope of the one-particle spectrum and the width of the two-particle correlator. To this end, we calculate first the building blocks of the ZP-algorithm. Having explicit expressions for , and in terms of simple polynomials will simplify our discussion considerably. The corresponding first order Pratt terms are

(25a) | |||||

(25b) | |||||

All higher order Pratt terms can be calculated explicitly as averages over relative and average pair distributions. Details are given in Appendix B. The momentum-independent terms read

(26a) | |||||

(26b) | |||||

(26c) | |||||

(26d) |

Here, denotes the greatest integer not larger than (‘floor’), and the least integer not smaller than (‘ceiling’). The notational shorthands and range between 0 and 1 depending on the phase space localization of the wavepacket centers, i.e., they map the whole parameter space , of the model (24) onto a finite region. The momentum-dependent terms are given by

(27a) | |||

(27b) | |||

(27c) | |||

(27d) |

The comparison of the present model calculation with the Zajc model (13) is not straightforward. As mentioned in the sequel of Eq. (19), the wavepackets used in both models can be compared by setting the wavepacket centers . However, the integral over in (12b) performs an average over the positions while in the model (24) we do not average over the position but over the centers of wavepackets , . Due to these different starting points, the Zajc model is not a simple limiting case of (24). Nevertheless, main features of the Zajc model (13) can be reproduced qualitatively in the present model. The leading contribution of the momentum-independent Pratt terms shows again a power law behaviour. Also, the -dependence of the terms and of the Zajc model is recovered in certain limiting cases,

(28a) | |||||

(28b) |

In general, however, the -dependence of the terms and is much weaker in the present model. Especially, there is no limit in which both and . These differences between both models may provide a first idea about the extent to which the choice of the model distribution affects our conclusions.

### iv.1 Weighting multiparticle contributions

The normalization is not a direct physical observable, but it determines the weights with which multiparticle correlations contribute to the particle spectra. To see this, we consider the one-particle spectrum,

(29) |

The -th order contributions are normalized to one, and the weights add up to unity,

(30a) | |||||

(30b) |

The lowest order contribution , by which the one-particle spectrum is typically approximated, contributes a fraction only, the value characterizes the importance of higher order contributions. For a quantitative analysis we now determine the dependence of the normalization and the weights on the event multiplicity and the phase space density of the emission region.

We consider the terms , the building blocks of . For the present model, these are given in (26a). The factor in this equation ranges between and , and can be written as

(31a) | |||||

(31b) |

The correction factor appears only linearly in the expressions for , rather than as an -th power, and it is of order (it is exactly for a choice of parameters and such that ). This allows for the approximation

(32a) | |||||

(32b) |

with which the normalization takes the simple form

(33) | |||||

Here, the combinatorical factors denote the number of permutations of elements which contain exactly cycles. They are commonly refered to as Stirling numbers of the first kind. We have used their generating function in terms of -functions [25].

We can now determine the weights of multiparticle contributions to the one- and two-particle spectra. For the one-particle spectrum, we find by inserting (33) into the ZP-algorithm

(34a) | |||||

(34b) | |||||

(34c) |

The approximation in the last line is valid for large multiplicities, when . Similarly, the weights for the different contributions to the two-particle momentum spectrum can be calculated. Using the power law behaviour , we find

(35a) | |||

(35b) |

where

(36) | |||||

Again, the approximation in the last line is valid for large multiplicities, when . To sum up: multiparticle correlations account for a fraction of the one-particle spectrum. For the two-particle spectrum, they are somewhat more important: the pure pairwise correlations receive only a weight .

For sufficiently large event multiplicities, the weights and of multiparticle contributions are not separate functions of and , but functions of the product only. The physics entering can be most easily illustrated in the large phase space volume limit, when

(37) |

We hence call a “phase space density of emission points”. This notion should not be taken too literally: the product of the volumes of three-dimensional spheres in position and momentum space is , and hence, is for large sources approximately a factor 10 larger than the particle number per unit phase space cell. Also, for realistic source sizes, the value of deviates significantly from the approximation (37), and a calculation of starting from (32b) is preferable.

One can ask whether in the large limit, the normalization becomes a function of only. To clarify this, we recall that a product with has a -limit if and only if converges. For the normalization (33), we find

(38) |

There is no physical reason why should remain finite. It is not an observable. What matters in a quantitative study of multiparticle correlation effects is the phase space density of emission points and not the particle multiplicity, what matters are the weights and , and not the normalization .

We finally estimate realistic values for the phase space density in heavy ion collisions. For a choice of model parameters fm, fm, MeV say, we find . With multiplicities of like-sign pions in the hundreds, this leads to . The present static, spherically symmetric model is however unrealistic in so far that it does not take the strong longitudinal expansion into account which significantly increases the volume out of which particles are emitted. From these heuristic considerations, we expect realistic phase space densities to lie in the range

(39) |

Depending on the precise value of in this range, the importance of multiparticle contributions to the one- and two-particle spectrum varies significantly. For , higher order contributions start dominating, while they account for 10 % of the signal if .

### iv.2 The momentum dependence of multiparticle contributions

The two -dependent terms and in (27) control the dependence of on the relative pair momentum and the average pair momentum . They determine the momentum dependence of all particle spectra. In Figure 2, these factors are shown as a function of for a source with fm and fm. The first message of this plot is that even for very high multiparticle contributions (e.g. ), the momentum dependence of all building blocks can be calculated exactly. Secondly, the factors and show the interesting property that irrespective of the value chosen for (and hence for ), they rapidly converge to an -independent quantity. In contrast to the limiting case (28a), the -dependence of is much weaker for realistic model parameters , ,

(40) |

Analogously, for realistic model parameters , ,

(41) |

while for the limit (28b) of the parameter space, a strong -dependence remains.