Bose-Einstein effects

A crude option for the simulation of Bose-Einstein effects is included since long, but is turned off by default. In view of its shortcomings, alternative descriptions have been introduced that try to overcome at least some of them [Lön95].

The detailed BE physics is not that well understood, see e.g. [Lör89]. What is offered is an algorithm, more than just a parameterization (since very specific assumptions and choices have been made), and yet less than a true model (since the underlying physics picture is rather fuzzy). In this scheme, the fragmentation is allowed to proceed as usual, and so is the decay of short-lived particles like $\rho$. Then pairs of identical particles, $\pi^+$ say, are considered one by one. The $Q_{ij}$ value of a pair $i$ and $j$ is evaluated,

$$Q_{ij} = \sqrt{ (p_i + p_j)^2 - 4m^2} ~,$$ (260)

where $m$ is the common particle mass. A shifted (smaller) $Q'_{ij}$ is then to be found such that the (infinite statistics) ratio $f_2(Q)$ of shifted to unshifted $Q$ distributions is given by the requested parameterization. The shape may be chosen either exponential or Gaussian,

$$f_2(Q) = 1 + \lambda \exp \left( - (Q/d)^r \right), ~~~~r = 1~\mathrm{or}~2 ~.$$ (261)

(In fact, the distribution has to dip slightly below unity at $Q$ values outside the Bose enhancement region, from conservation of total multiplicity.) If the inclusive distribution of $Q_{ij}$ values is assumed given just by phase space, at least at small relative momentum then, with $\d ^3 p / E \propto Q^2 \, \d Q / \sqrt{Q^2 + 4m^2}$, $Q'_{ij}$ is found as the solution to the equation

$$\int_0^{Q_{ij}} \frac{Q^2 \, \d Q}{\sqrt{Q^2 + 4 m^2}} = \int_0^{Q'_{ij}} f_2(Q) \, \frac{Q^2 \, \d Q}{\sqrt{Q^2 + 4 m^2}} ~.$$ (262)

The change of $Q_{ij}$ can be translated into an effective shift of the three-momenta of the two particles, if one uses as extra constraint that the total three-momentum of each pair be conserved in the c.m. frame of the event. Only after all pairwise momentum shifts have been evaluated, with respect to the original momenta, are these momenta actually shifted, for each particle by the (three-momentum) sum of evaluated shifts. The total energy of the event is slightly reduced in the process, which is compensated by an overall rescaling of all c.m.-frame momentum vectors. It can be discussed which are the particles to involve in this rescaling. Currently the only exceptions to using everything are leptons and neutrinos coming from resonance decays (such as $\mathrm{W}$'s) and photons radiated by leptons (also in initial-state radiation). Finally, the decay chain is resumed with more long-lived particles like $\pi^0$.

Two comments can be made. The Bose-Einstein effect is here interpreted almost as a classical force acting on the `final state', rather than as a quantum mechanical phenomenon on the production amplitude. This is not a credo, but just an ansatz to make things manageable. Also, since only pairwise interactions are considered, the effects associated with three or more nearby particles tend to get overestimated. (More exact, but also more time-consuming methods may be found in [Zaj87].) Thus the input $\lambda$ may have to be chosen smaller than what one wants to get out. (On the other hand, many of the pairs of an event contains at least one particle produced in some secondary vertex, like a $\mathrm{D}$ decay. This reduces the fraction of pairs which may contribute to the Bose-Einstein effects, and thus reduces the potential signal.) This option should therefore be used with caution, and only as a first approximation to what Bose-Einstein effects can mean.

Probably the largest weakness of the above approach is the issue how to conserve the total four-momentum. It preserves three-momentum locally, but at the expense of not conserving energy. The subsequent rescaling of all momenta by a common factor (in the rest frame of the event) to restore energy conservation is purely ad hoc. For studies of a single $\mathrm{Z}^0$ decay, it can plausibly be argued that such a rescaling does minimal harm. The same need not hold for a pair of resonances. Indeed, studies [Lön95] show that this global rescaling scheme, which we will denote $\mathrm{BE}_0$, introduces an artificial negative shift in the reconstructed $\mathrm{W}$ mass, making it difficult (although doable) to study the true BE effects in this case. This is one reason to consider alternatives.

The global rescaling is also running counter to our philosophy that BE effects should be local. To be more specific, we assume that the energy density of the string is a fixed quantity. To the extent that a pair of particles have their four-momenta slightly shifted, the string should act as a `commuting vessel', providing the difference to other particles produced in the same local region of the string. What this means in reality is still not completely specified, so further assumptions are necessary. In the following we discuss four possible algorithms, whereof the last two are based strictly on the local conservation aspect above, while the first two are attempting a slightly different twist to the locality concept. All are based on calculating an additional shift $\delta\mathbf{r}_k^l$ for some pairs of particles, where particles $k$ and $l$ need not be identical bosons. In the end each particle momentum will then be shifted to $\mathbf{p}_i' = \mathbf{p}_i + \sum_{j \neq i} \delta \mathbf{p}_i^j + \alpha\sum_{k \neq i} \delta \mathbf{r}_i^k$, with the parameter $\alpha$ adjusted separately for each event so that the total energy is conserved.

In the first approach we emulate the criticism of the global event weight methods with weights always above unity, as being intrinsically unstable. It appears more plausible that weights fluctuate above and below unity. For instance, the simple pair symmetrization weight is $1 + \cos (\Delta x \cdot \Delta p)$, with the $1 + \lambda \exp(-Q^2R^2)$ form only obtained after integration over a Gaussian source. Non-Gaussian sources give oscillatory behaviours.

If weights above unity correspond to a shift of pairs towards smaller relative $Q$ values, the below-unity weights instead give a shift towards larger $Q$. One therefore is lead to a picture where very nearby identical particles are shifted closer, those somewhat further are shifted apart, those even further yet again shifted closer, and so on. Probably the oscillations dampen out rather quickly, as indicated both by data and by the global model studies. We therefore simplify by simulating only the first peak and dip. Furthermore, to include the desired damping and to make contact with our normal generation algorithm (for simplicity), we retain the Gaussian form, but the standard $f_2(Q) = 1 + \lambda \exp(-Q^2R^2)$ is multiplied by a further factor $1 + \alpha \lambda \exp(-Q^2R^2/9)$. The factor $1/9$ in the exponential, i.e. a factor 3 difference in the $Q$ variable, is consistent with data and also with what one might expect from a dampened $\cos$ form, but should be viewed more as a simple ansatz than having any deep meaning.

In the algorithm, which we denote $\mathrm{BE}_3$, $\delta\mathbf{r}_i^j$ is then non-zero only for pairs of identical bosons, and is calculated in the same way as $\delta\mathbf{p}_i^j$, with the additional factor $1/9$ in the exponential. As explained above, the $\delta\mathbf{r}_i^j$ shifts are then scaled by a common factor $\alpha$ that ensures total energy conservation. It turns out that the average $\alpha$ needed is $\approx -0.2$. The negative sign is exactly what we want to ensure that $\delta\mathbf{r}_i^j$ corresponds to shifting a pair apart, while the order of $\alpha$ is consistent with the expected increase in the number of affected pairs when a smaller effective radius $R/3$ is used. One shortcoming of the method, as implemented here, is that the input $f_2(0)$ is not quite 2 for $\lambda = 1$ but rather $(1 + \lambda) (1 + \alpha \lambda) \approx 1.6$. This could be solved by starting off with an input $\lambda$ somewhat above unity.

The second algorithm, denoted $\mathrm{BE}_{32}$, is a modification of the $\mathrm{BE}_3$ form intended to give $f_2(0) = 1 + \lambda$. The ansatz is

$$f_2(Q) = \left\{ 1 + \lambda \exp(-Q^2R^2) \right\} \left\... ...xp(-Q^2R^2/9) \left( 1 - \exp(-Q^2R^2/4) \right) \right\} ~,$$ (263)

applied to identical pairs. The combination $\exp(-Q^2R^2/9) \left( 1 - \exp(-Q^2R^2/4) \right)$ can be viewed as a Gaussian, smeared-out representation of the first dip of the $\cos$ function. As a technical trick, the $\delta\mathbf{r}_i^j$ are found as in the $\mathrm{BE}_3$ algorithm and thereafter scaled down by the $1 - \exp(-Q^2R^2/4)$ factor. (This procedure does not quite reproduce the formalism of eq. (), but comes sufficiently close for our purpose, given that the ansatz form in itself is somewhat arbitrary.) One should note that, even with the above improvement relative to the $\mathrm{BE}_3$ scheme, the observable two-particle correlation is lower at small $Q$ than in the $\mathrm{BE}_0$ algorithm, so some further tuning of $\lambda$ could be required. In this scheme, $\langle \alpha \rangle \approx -0.25$.

In the other two schemes, the original form of $f_2(Q)$ is retained, and the energy is instead conserved by picking another pair of particles that are shifted apart appropriately. That is, for each pair of identical particles $i$ and $j$, a pair of non-identical particles, $k$ and $l$, neither identical to $i$ or $j$, is found in the neighbourhood of $i$ and $j$. For each shift $\delta\mathbf{p}_i^j$, a corresponding $\delta\mathbf{r}_k^l$ is found so that the total energy and momentum in the $i,j,k,l$ system is conserved. However, the actual momentum shift of a particle is formed as the vector sum of many contributions, so the above pair compensation mechanism is not perfect. The mismatch is reflected in a non-unit value $\alpha$ used to rescale the $\delta\mathbf{r}_k^l$ terms.

The $k,l$ pair should be the particles `closest' to the pair affected by the BE shift, in the spirit of local energy conservation. One option would here have been to `look behind the scenes' and use information on the order of production along the string. However, once decays of short-lived particles are included, such an approach would still need arbitrary further rules. We therefore stay with the simplifying principle of only using the produced particles.

Looking at $\mathrm{W}^+ \mathrm{W}^-$ (or $\mathrm{Z}^0 \mathrm{Z}^0$) events and a pair $i, j$ with both particles from the same $\mathrm{W}$, it is not obvious whether the pair $k,l$ should also be selected only from this $\mathrm{W}$ or if all possible pairs should be considered. Below we have chosen the latter as default behaviour, but the former alternative is also studied below.

One obvious measure of closeness is small invariant mass. A first choice would then be to pick the combination that minimizes the invariant mass $m_{ijkl}$ of all four particles. However, such a procedure does not reproduce the input $f_2(Q)$ shape very well: both the peak height and peak width are significantly reduced, compared with what happens in the $\mathrm{BE}_0$ algorithm. The main reason is that either of $k$ or $l$ may have particles identical to itself in its local neighbourhood. The momentum compensation shift of $k$ is at random, more or less, and therefore tends to smear the BE signal that could be introduced relative to $k$'s identical partner. Note that, if $k$ and its partner are very close in $Q$ to start with, the relative change $\delta Q$ required to produce a significant BE effect is very small, approximately $\delta Q \propto Q$. The momentum compensation shift on $k$ can therefore easily become larger than the BE shift proper.

It is therefore necessary to disfavour momentum compensation shifts that break up close identical pairs. One alternative would have been to share the momentum conservation shifts suitably inside such pairs. We have taken a simpler course, by introducing a suppression factor $1 - \exp(-Q_k^2 R^2)$ for particle $k$, where $Q_k$ is the $Q$ value between $k$ and its nearest identical partner. The form is fixed such that a $Q_k = 0$ is forbidden and then the rise matches the shape of the BE distribution itself. Specifically, in the third algorithm, $\mathrm{BE}_m$, the pair $k,l$ is chosen so that the measure

$$W_{ijkl} = \frac{ (1 - \exp(-Q_k^2 R^2))(1 - \exp(-Q_l^2 R^2))}% {m^2_{ijkl}}$$ (264)

is maximized. The average $\alpha$ value required to rescale for the effect of multiple shifts is 0.73, i.e. somewhat below unity.

The $\mathrm{BE}_{\lambda}$ algorithm is inspired by the so-called $\lambda$ measure [And89] (not the be confused with the $\lambda$ parameter of $f_2(Q)$). It corresponds to a string length in the Lund string fragmentation framework. It can be shown that partons in a string are colour-connected in a way that tends to minimize this measure. The same is true for the ordering of the produced hadrons, although with large fluctuations. As above, having identical particles nearby to $k,l$ gives undesirable side effects. Therefore the selection is made so that

$$W_{ijkl} = \frac{ (1 - \exp(-Q_k^2 R^2))(1 - \exp(-Q_l^2 R^... ...mutations)}} (m_{ij}m_{jk}m_{kl},% m_{ij}m_{jl}m_{lk},\ldots)}$$ (265)

is maximized. The denominator is intended to correspond to $\exp(\lambda)$. For cases where particles $i$ and $j$ comes from the same string, this would favour compensating the energy using particles that are close by and in the same string. We find $\langle\alpha\rangle\approx 0.73$, as above.

The main switches and parameters affecting the Bose-Einstein algorithm are MSTJ(51) - MSTJ(57) and PARJ(91) - PARJ(96).