On Multiparticle Production in Very High Energy Scattering

A phenomenological model of particle production and hadronisation in high energy collisions is formulated using Dirac fields in Yukawa-like interaction and the resulting stochastic equation is solved numerically. Different initial conditions are used to compare particleparticle ( ψ ψ ) and particle-antiparticle ( ψ* ψ ) interactions. It is shown that in this simplified view, there is a clear difference between the final multiplicity distributions resulting from the two initial conditions. To model the restricted phase space (limited pseudorapidity) measurements in experiment, a “loss” function is also proposed to account for the undetected particles close to the beam line.


Introduction
Particle production is a high energy phenomena that stems from the Einstein's relativistic physics. It is highly relevant in the domain of quantum systems interacting at energies that are significantly larger than the mass scales of known matter. Experiments done at the LHC and all previous collaborations aim to study matter at smaller and smaller length scales with each generation to probe the physics of the early universe. Particle production and hadronisation is an important aspect of what happened after the big bang and our understanding of it is incomplete. Standard model physics is perturbative by design and we do not have a complete map of what goes on in the large coupling sector of Quantum Chromodynamics (QCD), which is the relevant theory for particle production.
Our understanding of QCD in terms of perturbation theory is limited due to the running nature of the coupling constant. Studies of branching models circumvent this problem in an interesting way. Branching models derived from perturbation theories result in stochastic equations of evolution that are independent of the coupling. It is therefore a valid assumption that these equations are still applicable when the coupling is strong [1]. Over the years, such branching models have proved to be a helpful guide to studying branching processes in high energy experiments † . This is one of the non-perturbative approaches employed to model physics of the soft sector. The study of hadronic multiplicity at high energies is a highly contested topic as there is no first principles explanation from QCD to the observed data. This further strengthens the case for stochastic phenomenological studies that take into consideration the relevant physics while making certain simplifying assumptions to deal with the problem at hand.
The most popular parton branching equation was first discussed by Giovannini [2]. His branching model has been extensively studied in the final decades of the last century [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. Such branching models work under the following assumptions: particle production occurs in three distinct steps. First, partons showers are created through the branching process birthing from a bunch of highly virtual partons. Next, these partons hadronise through an unknown process. Finally, these hadrons are detected in our experiment subject to detector design and efficiency. This final step is best handled through combinatorics. Unfortunately, the process of hadronisation is still not very well understood.
In this work, we propose a stochastic branching model based on these principles with the interaction term λ ψ* φ ψ. This is similar to the effective nuclear force described by Yukawa interaction and an exchange of pions. Although Yukawa processes have been studied in quite some detail in the literature, as a branching process they are mostly studied as a φ 3 process [18]. This is one of the more successful models that has the Furry distribution as its solution [19,20]. However, since this is a scalar theory, there is no distinction between particles and antiparticles. It is therefore impossible to probe the difference, if any, between particle-particle and particle-antiparticle interactions [21]. Giovannini's parton branching equation also fails at this count. It does not make any distinction between quark and anti-quark branching probabilities. We shall address these shortcomings in this work and showcase how this distinction between particle-particle and particle-antiparticle interactions manifests itself in our model. Our setup aims to model the difference between pp and p̄p collisions [21] which Giovannini's parton branching equation [2] does not address. Multiplicity distributions from a quantum field theory can be derived using the Alterelli-Parisi equation [22] (also known as the DGLAP equation)

Stochastic Branching Models
where the integral runs from x to 1 and the derivative on the left is with respect to t. D(x,t) is the 1-particle inclusive distribution and P(x) is the splitting function that dictates the outgoing distribution of incoming momentum at any given vertex. The variables (z, x) are the momentum fraction of the state in question whereas t is the evolution parameter ‡ Figure is borrowed from reference [18]. measuring the virtuality of the initial state. 1-particle inclusive distribution D(x) describes the probability of producing 1 particle with momentum fraction x.
The physical idea is that we start out with a state that has very high virtuality and, as it branches into more and more particles, the final state approaches the mass shell (Fig. 1). There is a renewed interest in the DGLAP equation in recent years as it is used to study inclusive production and branching in conformal field theories [23].
An interesting note here is that the evolution equation for the multiplicity distribution can also be derived as a stochastic branching process. The evolution parameter t can in this sense be treated as a "time" parameter. As discussed in [1] and [18], the φ 3 branching process simply follows the Furry branching process of pure birth. The Furry distribution has been used rather successfully to model particle production through geometric models [19,20]. In fact, it is one of the limiting cases of Giovannini's branching model. However, it still has limitations described previously pertaining to limited physical applicability at higher energies. It nevertheless proves the strength of stochastic methods to modelling particle production. It is this approach that we shall follow in our current work. A more detailed analysis using the DGLAP equation is left to a future study.

The Evolution Equation
There are three fundamental processes possible in our model dictated by the same fundamental vertex: ψ Bremsstrahlung, ψ* Bremsstrahlung and ψ* ψ pair production. The equation for stochastic evolution of this system can be easily derived. We assign a probability to each of the three kinds of branching processes: A for ψ Bremsstrahlung, A* for ψ* Bremsstrahlung and B for ψ* ψ pair production. A state with n particles, m antiparticles and o scalars has probability given by P mno and evolves in a small step dt as Therefore, P' mno (t) = nA (P mn,o-1 -P mno ) + mA* ( P mn,o-1 -P mno ) where we have suppressed the t dependence on the right-hand side for brevity. This equation does not have a closed form analytic solution. We therefore use numerical methods to solve it. There is a small caveat to this distribution. In a high energy experiment, what we measure in the end state is the number of charged particles. Since we wish to model this charged particle multiplicity distribution P n , we need to sum over all neutral scalars and count all charged particles. P n is given by Another important note here is that for particle-particle and particle-antiparticle interactions, our initial condition always contains 2 initial charged particles. Our model therefore produces non-zero values only for even multiplicities since all branching processes increase the particle number by 2. This is expected as these numbers are for full phase space, i.e., -∞ < η < ∞ § . Since most pp and p̄p data is restricted phase space due to detector design, we introduce a "loss" function to account for particles lost down the beam line.

The Loss Function
To account of particles lost outside the detector region and down the beam line, we can perform a convolution of our multiplicity distribution with a probability distribution that describes the loss. The final multiplicity distribution P̄n should be calculated as follows: where p i is the probability of losing i particles. Following Occam's razor, we assume the simplest form of p i : where Binomial (i, n+i, p) is the probability mass function of the binomial distribution. p i , therefore, gives us the probability of losing i particles given that n+i particles are produced assuming that p is the probability of losing the particle. Employing such a loss function gives us non-zero values for all multiplicities, as is observed in restricted phase space data for pp and p̄p collisions. Ideally, what we expect is that this probability of loss will increase as we constrain our pseudorapidity window since more and more particles go undetected.

Method for Numerical Integration
As mentioned earlier, the evolution equation does not have any closed form solution. Numerical methods are therefore required to further study our model and make contact with data. Numerical integration was done by performing 4th order Runge-Kutta integration (RK-4) for 4 different values of the step-size (h = 0.1, 0.05, 0.01, 0.005) and then using polynomial interpolation to find the values for h = 0. This was done mainly to improve normalisation as RK-4 regularly under-predicts the solution. This error is dependent on the step size h and therefore, treating the total probability of the final state found by the integrator as a function of h should give better results through interpolation. Indeed, this § The UA5 data used for the fit is for a restricted pseudorapidity window. However, the full phase space data reported by UA5 is all even multiplicities as well, as anticipated by our distribution. resulted in significant improvements in normalisation while saving on computational complexity. Most of the numerical work was done in Python. However, C++ code was wrapped using Simplified Wrapper and Interface Generator (SWIG) [24] for performance and efficiency gains whenever required. For example, integrating from t = 0 to t = 4.0 took ~ 43 minutes in Python. The same calculation with wrapped C++ code in Python took ~ 13 seconds.

Results and Discussion
The following figure shows the effect of having different values of A and A* on the multiplicity distribution for ψ ψ interactions and ψ* ψ interactions. This is one of the main highlights of this work: from these plots we can conclude that in our model, matter-antimatter asymmetry shows up as a difference in the multiplicity distributions for the particle-particle and particle-antiparticle collisions. Some preliminary fitting using p̄p data from the UA5 collaboration [25] for pseudorapidity window -3 < η < 3 gives us very different values for A and A* as is tabulated below. For this fit, we get a very promising value for χ 2 / d.o.f. of 53.27/75. Our model, therefore, successfully captures the feature we were probing for. A detailed analysis of this result should include a fit to pp data from the LHC for the same pseudorapidity window and same energy. This is one of the approaches currently under consideration.

Conclusion
A phenomenological branching model for particle production that captures matterantimatter asymmetry is presented. The stochastic equation is solved numerically and preliminary fits to data from the UA5 collaboration seem promising. Future work on the model would involve performing Monte Carlo Markov Chain fitting to do Bayesian parameter estimation and fitting to other datasets. It is expected that the parameter p in the "loss" function should depend only on the pseudorapidity window over which the multiplicity is measured. Further, fitting to pp data from the LHC for the same energy and pseudorapidity window should prove to be a good test for the model. Once the model is rigorously tested against datasets, there are all kinds of phenomenological studies that can be performed. An interesting direction of recent study has been towards modified combinants of multiplicity distribution [26][27][28]. Some preliminary calculations show promising results from our distribution towards these interesting quantities. These directions are currently being explored.