Baryon Anticorrelations in PYTHIA

We present a computational investigation of a problem of hadron collisions from recent years, that of baryon anticorrelations. This is an experimental dearth of baryons near other baryons in phase space, not seen upon examining numerical Monte Carlo simulations. We have addressed one of the best known Monte Carlo codes, PYTHIA, to see what baryon (anti)correlations it produces, where they are originated at the string-fragmentation level in the underlying Lund model, and what simple modifications could lead to better agreement with data. We propose two ad-hoc alterations of the fragmentation code, a “one-baryon” and an “always-baryon” policies that qualitatively reproduce the data behaviour, i.e anticorrelation, and suggest that lacking Pauli-principle induced corrections at the quark level could be the culprit behind the current disagreement between computations and experiment.


I. INTRODUCTION A. Baryon anticorrelations
The baryon anticorrelation problem was exposed in a series of works by the ALICE collaboration [1,2], investigations derived therefrom [3], and further corroborated by the STAR collaboration [4].Basically it is a severe disagreement about the probability of two baryons exiting a high-energy collision closely in phase-space (similar momenta).Computational simulations using event generators fail to see an anticorrelation clearly present in the data.To understand it, let us first define the appropriate two-particle correlations that are measured.
Let us begin with the definition of the single and pair densities of species α and β as where N α 1 represents the number of particles of species α, N αβ 2 represents the number of pairs of particles of species α and β, and η i and ϕ i represent the corresponding pseudorrapidity, η, and azimuthal angle, ϕ, of the involved particle.Clearly the magnitude measures the degree of correlation between particles of species α in phase space bin (η 1 , ϕ 1 ) and particles of species β in phase space bin (η 2 , ϕ 2 ).Two-particle correlations are usually reported, normalizing the Eq. ( 3) magnitude, as normalized second order cumulants, and in relative separation in pseudorapidity, ∆η = η 1 − η 2 , and azimuth, ∆ϕ = Experimentally the expression in the denominator is usually obtained directly in ∆η, ∆ϕ by using the mixed events technique which builds the ρ 2 distribution considering one pair component from one event and the second pair component from a different event.By using this technique, the expression (3) is identically zero, because particle from different events are not correlated, and is obtained.Then, what is usually reported as measured two-particle correlation function is where S is the ratio of the distribution density in relative rapidity ∆η and azimutal angle ∆ϕ of presumably correlated pairs (because they are taken from the same event), to the total number of pairs, S(∆η, ∆ϕ) = 1 and B is the equivalent distribution of uncorrelated pairs (stemming from separate events, or "mixed" in the field's jargon) The mixed event distribution in Eq. ( 8) corrects for the effects of limited acceptance in pseudorapidity on S but, apart from a factor scale, it does not affect its azimuthal shape.For this reason, S may be studied alone to investigate the behaviour in ∆ϕ of the correlation function C. Now let us describe the baryon anticorrelation problem.Contrasting experimental data with simulations from three tunes of the PYTHIA Monte Carlo event generator (two Perugia and the Monash one), the outcome was rather disappointing [1].The disagreement is not particularly due to PYTHIA because its substitution by PHOJET, an alternative Monte Carlo event generator, did not bring improvement.Figure 1, taken from [2], shows the twoparticle correlation functions.Correlation functions are integrated over ∆η (and normalized), which is a more convenient to compare several results.FIG.1: Left: Typical meson-meson correlation function from ALICE data [2] as a function of ∆ϕ, showing a large forward peak (two Kaons are more likely to be detected with ∆ϕ = ϕ 2 − ϕ 1 0) and a sizeable backward peak at ∆ϕ = π.Simulation and experimental data are in reasonable qualitative agreement.Center and right: Experimental baryon-baryon+anti-baryon-anti-baryion correlation functions become negative (data, grey squares) when ∆ϕ 0, showing a forward anticorrelation.Monte Carlo simulations seem unable to reproduce that behaviour.Note that this happens not only for identical baryons such as pp (center) but also for distinguishable ones such as pΛ (right).Reproduced from [2] by the ALICE collaboration under the terms of the Creative Commons License 4.0 (http://creativecommons.org/licenses/by/4.0/);no changes have been effected.(The red and green lines correspond to PYTHIA 6 with two Perugia tunes, and play no further role in our discussion).
The OPAL collaboration discussed in [5] what was known on baryon correlations at the time of wrapping up the LEP experimental program.Their conclusion was that, due to limited statistics and insensitivity to dynamics by lacking rapidity in the analysis, the correlations due to discrete quantum number (flavor, charge, baryon number conservation were clearly visible, but not much beyond that: quark/fragmentation-level correlations were not clearly established.The situation has now changed due to the impressive statistics and tracking capabilities of more modern devices, and dynamical effects will need to be addressed.

B. Common use of the PYTHIA software to predict correlations
PYTHIA is a Monte Carlo event-generator software program, widely used to simulate possible outcomes of experiments (like collisions between two objects, be it fundamental particles or nuclei) and help understand detector responses.Through various physical model mechanisms, PYTHIA evolves the initial input state into a number of output particles.These are listed by their flavor, spin and momentum in its original format described in [6].
Such simulation packages are popular because of their convenience to explore possible experimental outcomes when designing a measurement or during its later analysis.Moreover, PYTHIA allows one to quickly test theoretical ideas, for example, adding some Beyond the Standard Model feature, to start assessing its impact on the data, even to guess physics at energies that current colliders cannot reach.Finally, finding a way to simulate a process can help understand it, so the shortcomings of PYTHIA can sometimes reveal shortcomings in current knowledge.
PYTHIA addresses those goals by quickly generating a great number of simulated events.The main requirements on the software are speed of execution and fidelity of the data produced at the end of the simulation compared to the data from experiment.On the contrary, fidelity to physical behaviour in the intermediates steps is secondary: the description of physical mechanisms used by PYTHIA can be highly phenomenological, notwithstanding its inspiration in the ideas of flux tubes forming among color sources and breaking in Quantum Chromodynamics [7].It is therefore more descriptive than predictive, and can disagree with experiment when confronted with a new type of data.The phenomenological parametrizations deployed tend to use many parameters not set by any theory, but in order to fit the data.There are several possible ensembles of their values, called "tunes".
Returning to the left plot of Fig. 1 for meson-meson correlations, we see that the results of the PYTHIA simulation are not far from experiment.The usual features are a strong correlation at small phase-space separation (∆ϕ 0 ∆η) for particles that exit the collision together, perhaps from the same jet; a positive correlation for back-to-back particles with opposite azimuthal angles (a peak at ∆ϕ = π from two-jet events) and often (not visible in this flat depiction of the ∆ϕ dependence alone) a ridge-shaped correlation extending in ∆η often attributed to string-like flux tubes stretched by the separating nuclei.
For baryon and anti-baryon correlations, shown in the center and right panels of Fig. 1, however, the Monte Carlo simulation is unable to obtain a satisfactory description.The worst agreement is for baryon-baryon correlations, where PYTHIA does not even produce the right qualitative behaviour.At around ∆ϕ = 0, all the tunes considered predict a peak instead of the salient depression (negative correlation, that is, anticorrelation) shown by the experiment.We would like to understand what modification of PYTHIA would suppress that peak in baryon-baryon correlations and produce an anticorrelation instead.
Changing the "tune" (particular parameter set) of the simulation does not really improve the predicted correlations.Tunes can balance the influence of the different mechanisms inside PYTHIA, but cannot change them drastically.Whatever reason for the failure of PYTHIA in reproducing the data, it has to be deeper.Therefore we will limit ourselves to using the Monash tune in this work.

C. Recent attempts: the afterburners method
To correct the results of the simulation, a postprocessing procedure of the data has been proposed [3].In the following, it will be named the "afterburners" procedure (in analogy to motor technology).This correction to the PYTHIA output produces encouraging results for generic kinematics, but in this work we will only evaluate the very specific problem of proton-proton anticorrelation at ∆ϕ = 0, only one of several physical situations [3] but a problematic one: the correction to the pp correlation (see Fig. 2 in that article) is insufficient to account for the anticorrelation.
The idea behind this procedure is to write in addition a piece of Fortran code that from the output of PYTHIA associates to each pair of particles a weight.The same-event correlation function S is computed and binned as an histogram, which allows to weight each pair.In the simplest situation (no coupled channels, non identical particles), the weight is equal to the square modulus of the following wave function with The variables k * and r * are the momentum and position vectors in the pair rest frame, γ is Euler's constant and f 0 , d 0 and a are parameters given for each pair.Further, F 0 and G 0 are the regular and singular s-wave Coulomb functions.For further explanations, see [3] and [8].
The implementation of this final-state wavefunction allows to correct the correlation function C to account for the antisymmetry upon exchanging the final-state baryons, unlike in the raw PYTHIA output.That computation of the weights is thus grounded both on quantum statistics and on final-state interactions due to Coulomb and nuclear forces.The exact formula giving the weight of a pair from the experiment, from the momentum and position of each particle (classically treated by PYTHIA, as explained below in section II) is quite complex.
We have attempted to redeploy that proposed method in our computations: the resulting modified correlation functions are shown in Fig. 2. One can see at first glance that the data behaviour is not fully reproduced: the ridge at ∆ϕ = 0 present in the panel labeled as (a) is not totally suppressed in (b); and further, the correlation now strongly peaks at (∆η, ∆ϕ) = (0, 0).Though this method does not seem to alleviate the baryon-anticorrelation problem in our pedestrian simulation, the existence of this new peak is not completely problematic.The middle plot of Fig. 1 shows a tiny peak in the middle of the anticorrelation valley (probably not a statistical blip, it is also clearly visible in three-dimensional renderings where it stands out for ∆η = 0.It appears that the afterburners correction might be predicting such peak.But (at least in our implementation) it does not produce the wide depression surrounding it.
Anyway, the afterburners procedure corrects the simulation output for all possible correlation functions: both pp and ππ are affected (in this case, the wavefunction is of course symmetrized).
Given this quick calculation, and the more extensive one [3] reported by the Warsaw Technical University investigators, we think an additional way to correct the PYTHIA output needs to be found.

II. FRAGMENTATION IN PYTHIA
To motivate to the reader the modification that we are going to effect, we briefly expose our understanding of the PYTHIA algorithms to produce final state hadrons.Then in the next section III we will show the correlations that PYTHIA produces before modification.

A. Simplified working method
A wise initial focus to minimize intervention in a time-tested code is to try to solve the problem at one single step of the algorithm.Since the afterburner postprocessing in subsection I C proved insufficient, this step would probably be at a previous level inside PYTHIA (The earlier a modification is applied, the larger its effect upon cascading over the rest of the simulation.) On the one hand, unstable-hadron decays cannot be the source of the problem, for they are simple processes in the code mimicking experimentally measured decays, so a disagreement with experiment is hard to expect.
On the other hand, computed meson-meson correlations broadly agree with experiment hinting that we should modify only the baryon-production algorithm.To effect a difference between baryons and mesons supposes that hadrons have emerged from the colliding system at this step.
Additionally, QCD amplitudes are also well understood so the perturbative matrix elements affecting the initial stages of the process are unlikely to be the culprit.
So we adopt the point of view that if there is one single mechanism to be corrected in order to obtain the experimental behaviour of the correlation function C for proton pairs in Fig. 1, it must occur at the step of the formation of hadrons.
Because of the reported pΛ anticorrelation (right panel of Fig. 1), the wanted mechanism has to be generic, not tied to the nucleon.
We then proceed in the next subsection to summarize the way [9] PYTHIA obtains prompt hadrons given initial partons.

B. Lund string model
The PYTHIA level just before the formation of hadrons handles quarks and gluons organized into color singlets.Any interactions involving different such singlet clusters had been taken into account earlier during the parton level computations and can probably be ignored for hadronization.At the stage of hadron formation, the code takes as input a set of partons and returns hadrons, modelling low-energy QCD interactions.Due to its complexity and for speed of algorithm execution, PYTHIA uses a phenomenological model called the Lund string model.Being far from perfect, this is convenient to treat hadron formation with sufficient richness of detail while maintaining acceptable turnout.
Summarizing the next few paragraphs, PYTHIA views the input color singlet made of quarks and gluons as a string that fragments into pieces (the hadrons).The simplest possible string is made of one quark, one antiquark and no gluons.For brevity, we will describe the process of fragmentation only in that case.For a complete description, we refer to the literature [9].
In the internal machinery of PYTHIA, each particle is modelled semiclassically, with exactly defined position and momentum, at odds with Heisenberg's uncertainty.The Lund model thus makes sense on average; this is convenient if a great number of fragmentations are performed, as is the regular usage of PYTHIA, at the cost of microscopic interpretability.
A model string stretches between a quark and an antiquark of opposite color, called "ends"; their flavor is not significant for the QCD interaction.Due to our understanding of QCD, the binding energy is taken as proportional to the distance that separates the ends.
If a quark-antiquark pair of the right color appeared at some point of the string (a tunnel-like effect), each parton of the pair could bind with an end to form its color singlet.Thus, the system would become two strings.If the total length of these two new strings was lower than the length of the original one, the situation would be energetically favorable.
The momentum of these partons emerging from the vacuum is chosen such that one of the daughter strings turns into a physical hadron of a certain flavour and invariant mass.This constitutes "an iteration".During the next iteration, the remaining daughter string becomes the mother one.This loop stops when all of the string's energy has been taken away by hadrons.
Instead of creating a qq pair, PYTHIA can also opt to produce a diquark-antidiquark one.PYTHIA views the diquark made of two quarks as a single particle.Such diquark can form a color singlet with a quark, just like an antiquark would.The only difference is that the hadron produced is a baryon instead of a meson 1 .
The hadrons that have been produced by this process can be of any flavour, for example containing charm or other types of spatial or spin excitations.Most of them are going to decay in the next step of PYTHIA, to finally give particles stable against strong-force mediated decays, such as pions, nucleons, kaons, Lambdas etc.

III. CORRELATIONS AMONG STRINGS AND AMONG BARYONS IN PYTHIA FRAGMENTATION
In section II we described the passage from a string to prompt hadrons in terms of the color flow.Flavour is of lesser interest because at the end, most heavy-flavoured quarks later disappear following weak decays and are often not detected.
More determining is the momentum of the objects implied in the fragmentation.It so happens that the theoretical study of the hadrons' momentum distribution function is burdensome (so that how to apply some level of Pauli suppression depending on the momentum distributions is not obvious to us).Earlier work [9] addresses momentum distributions at the level of one iteration, but extracting the momentum distributions from the full fragmentation loop is not trivial.
Because theoretical computation of the momentum distributions (both string input and hadron output to the fragmentation process) is unassailable, plotting the output data from simulations is the only way to understand their global behaviour.

Hadron distributions
To qualitatively understand the proton-proton two-particle correlation function (6) in ∆ϕ, we put aside the denominator B, to focus on the numerator S, which is basically a two-particle distribution affected by the limited longitudinal, η, acceptance which introduces the characteristic triangular shape in ∆η.We expect the protons present in the final state to have been produced during the fragmentation of some string or by the decay of a heavier baryon, itself due to the fragmentation.Roughly, half of the baryons are supposed to decay into protons, half into neutrons and some into Λ hyperons.Thus, the distribution of protons in the final step should be comparable to the distribution of baryons just after the fragmentation step (respectively antiprotons to anti-baryons).
What we want to extract is the distribution S of the pairs of baryons coming from the same event.Such sameevent distribution will receive contributions from baryons ejected from the same string as well as baryons coming each from a separate string within the same event.S is then equal to the sum of the same-string baryon correlation plus the different-string baryon correlations.Because we seek to understand whether/why there are not suppressed correlations of baryons near to each other in phase space we will focus on baryons stemming from the same string, denoting the corresponding correlation with S frag .For later comparability with data, the study is carried out in the lab frame.
The resulting numerical data is shown in Fig. 3.Among the plots there, (a), (c), (d) and (f) will contribute to the proton-proton or antiproton-antiproton S (if they were followed to the very end of the fragmentation).
Independently of the baryon number combination, all pairs are positively correlated.The only structure that seems to break the isotropy in ∆ϕ is the peak at around (∆η, ∆ϕ) = (0, 0).Such peak accepts an easy interpretation: when a string fragments, baryons are preferentially produced along the same direction of the string, due to momentum conservation.Apart from that peak, (anti-)baryons from the pair are emitted with random ϕ in the lab frame.
This strong forward peak at zero relative angle and pseudorapidity is precisely what seems to be suppressed in experimental data for BB and B B correlation functions, and we thus see the unwanted feature is already present at the level of individual string fragmentation.
The baryon string (bottom row panels) also produces two curious ridges most visible in panel (d) around ∆η = ±4 for 13 TeV, which are outside ALICE's acceptance.

Internal string distribution
To ascertain whether the "wrong" baryon correlations in PYTHIA are produced already at the string level, before the actual production of final-state particles, we have computed string-string correlations and report them in this section.As advanced in Fig. 3, we distinguish two types of string, defined by their baryon number.If the initial two ends of a string are a quark and an anti-quark, the global baryon number is equal to 0, and we call it a meson string.If the ends hold a quark and a diquark, the baryon number is 1, and we call it a baryon string.Other cases happen to be so rare in PYTHIA that they do not affect statistical distributions and are not further considered.
Each string is viewed as a particle whose momentum (and thus its pseudorapidity η and azimuth ϕ) is the sum over the string's constituents.The angular dependence of the two-string correlation (6) (unlike the absolute normalization) can be inferred from the same-event correlation S str , which we plot in Fig. 4. In it we find that in the azimuthal dimension the distribution is mostly isotropic, except around ∆ϕ = π, where we find a sizable peak.
The relative importance of the isotropic and the peaking contributions varies with the participating strings.For meson-meson strings (a), the peak is tiny and the isotropic background offset is large.For baryon-baryon strings (c), it is the structure at ∆ϕ = π (now separated in two peaks at large relative pseudorapidity) that dominates the graph, and the azimuthal isotropic offset is negligible in comparison.In (b), for meson-baryon strings, the situation is intermediate.The pair of twin peaks at large relative pseudorapidity of order 10 or more in (c) is somewhat surprising.It must be that the simulated collision between two protons produces, not only isotropically emitted baryons, but also two opposite beams.Whether a program artifact or an effect of the parent colliding protons remains to be seen, but in any case it is not currently comparable to the experimental data due to the limited rapidity coverage of the ALICE detector.

Interpretation
To summarize, the picture in ∆ϕ is the following: strings are quite uncorrelated except that they are likely to appear with opposite momenta (subsection III 2).
Baryons and anti-baryons from the same string (subsection III 1) are likewise often randomly distributed, though also with a colinear component.A part of the explanation is that between every baryon pair an antibaryon needs to be produced (to balance quark number) and this erases azimuthal baryon-baryon correlations, for example.Baryons from the same string tend to be emitted in separate directions in the string's reference frame, but the Lorentz boost given to the string upon passing to the laboratory frame appears to overcome this effect, so that baryons stemming from the same string are actually strongly collimated and yield a forward peak in the final correlation.
The part of the hadrons that is uniformly distributed along ∆ϕ either from the same or from different strings are not supposed to provide any particular correlation in the final output.They only contribute to the isotropic offset, that influences the height of each peak or depth of each valley.As they do not change the qualitative behaviour, we do not discuss them any further in the following.
The fraction of baryons and anti-baryons from the same string, colinearly emitted, will cause a peak around ∆ϕ = 0 in the final picture.The baryons and anti-baryons from different strings will contribute to the isotropic offset if their strings are uncorrelated and to a peak at around ∆ϕ = π if their strings are anti-colinear.
Given that the number of uncorrelated strings is much larger than that of opposite-momentum strings, providing a large combinatorial background, it appears that the balance of the peaks at around ∆ϕ = 0, both for protonproton and for proton-antiproton, are made of the correlations between the (anti-)baryons emitted by the same string 2 .This suggests to us that a modification to bring the Monte Carlo in better agreement with experimental data would be most efficient if applied to the production of baryons from the same string (that are naturally more affected by Fermi statistics, even starting at the quark level).

IV. SUGGESTED MODIFICATIONS TO PYTHIA THAT SEEM TO REPRODUCE BARYON ANTICORRELATIONS AS SEEN IN DATA
A. First modification: one-baryon policy Following the explanation in subsection III 3, the existence of the peak at around ∆ϕ = 0 in pp + pp correlation functions might be linked to strings giving birth to more than one hadron.
To test this idea, we completely forbid it in the program with a Von Neumann rejection step: anytime that a string fragments into more than one baryon (idem for anti-baryons), the fragmentation is recalculated.

Physical justification
One can see this procedure as a way to impose the Pauli principle to baryons (which are fermions).The algorithm is not influenced by the nature of the baryons nor their spin.Instead, it could be influenced by the Fermi statistics of the quarks themselves.The similarity of the center (pp) and right (pΛ) plots in Fig. 1 hints that no difference should be made among species of baryons.At the baryon level this may sound strange since p and Λ are distinguishable, but at the quark level, since both contain u and d valence quarks may be more understandable.
In the current PYTHIA 8.3 version, Fermi-Dirac statistics does not seem to be implemented in any way.It would have been surprising that the absence of such an important principle was not falsified at some point by contrasting the generated data against experiment: it so has happened that multiple-baryon final states have not been thoroughly analyzed until this last decade.
The appropriate theoretical procedure would be to anti-symmetrize an appropriate multi-fermion wavefunction or quantum amplitude, so we write down some formalism that could serve as illustration to fix ideas, although our computer-algorithm modifications will be much cruder.A possible output of a fragmentation [11] is associated with such an amplitude, derived from the wavefunctions of the involved entities through Were A, B and C are strings or hadrons, which are equivalent in this description.γ a bc (r, w) is a function that can be computed with methods discussed in the literature [11].
It is easy to write a (perhaps long) expression for the semi-inclusive two-baryon production amplitude of the form That allows to anti-symmetrize the baryon part of the wavefunction as in an atomic Slater determinant Taking the perhaps drastic but much simplifying approximation | M | |M | → 0 leads us to forbid the production of two baryons from the same string.This is the content of the first "one-baryon only" correction.
If one wanted to apply less dramatic simplification, one would have to redevelop parts of the Lund string model formalism to make it compatible with, for example, the formalism of [11].

Results
That simplest procedure of imposing a one-baryon per string policy can be very time consuming.Combining it with the afterburners correction worsens this and additionally leads to technical difficulties that makes it unpractical to gather sufficient Monte Carlo statistics in a limited university departmental cluster.pp one.
In the afterburners procedure, the lengthy code, together with the number of pairs increasing with the square of the number of events, can be a drawback.It brings about a rise of the whole time of execution, increases the frequency of software crashes and produces random pairs that occasionally give nonsense results (like negative weights) and must be discarded.If the afterburners method is used alone this is not dramatic and can be handled.But the difficulties compound if this addition to PYTHIA is combined with another.
As it is, we renounce to deploy the two methods together.It would be interesting to find a time-economic implementation, but we have not managed to design one in a reasonably small time.Our one-baryon policy correction gives results that are sufficiently encouraging to confirm the relevance of our method and suggest that further work will be interesting.
Figure 5 (right top plot marked as b) shows that the first demand on the output is fulfilled: the peak at ∆ϕ = 0 has been suppressed by imposing to each string the one-baryon policy.Indeed, that peak seems to be produced by baryons coming from the same string.Still, this policy has a major side effect, which is to strongly reduce the number of protons in the final state (and also of Λ hyperons and generically all baryon species).This was predictable because the method flatly reduces the average number of baryons produced by each string.B. Second, balancing correction: no baryonless fragmentation or "always baryons" policy.
To repurpose PYTHIA to correct for the dearth of protons in the final state, following the one-baryon policy, a second correction is tried.
Recognizing that the overall baryon number has been reduced because each fragmentation is now constrained to produce one baryon at most, a simple way to rebalance it is to force each meson string to fragment into exactly one baryon and one anti-baryon (two baryons and one anti-baryon for baryon strings).
This second "always baryons" correction is a purely practical expedient that no physical explanation seems to justify within the string model.And yet, the results are satisfactory.Figure 5 (lower-left plot marked c) shows that the relative features (ridges, valleys, absence of peaks) are equivalent to those of the one-baryon policy (upper right plot marked b), and it is only the overall number of baryons and anti-baryons which is increased.Any small remaining differences can probably be solved by a change of PYTHIA tune: indeed, each change to PYTHIA must induce a new calibration of its tune, but we leave this for future investigators and comment on it only briefly.
Among the parameters to be tuned, one should be carefully considered: P Q→QQ , the probability to engender diquark/anti-diquark pairs emerging from the vacuum during the string fragmentation, instead of quark/anti-quark ones.Changing this relative ratio is not independent of restraining the number of baryons emitted during the fragmentation.With our second always-baryons modification, this number of baryons is imposed so this parameter has no impact on the final output.
Nevertheless, an unwise choice of this P Q→QQ parameter has an unwelcome consequence.If one fragmentation iteration is unlikely to give exactly one baryon, it will be relaunched multiple times until it branches into a solution that satisfies the one-baryon policy.This may affect dramatically the length of execution of each fragmentation program loop.
But if only the first "one-baryon policy" correction is imposed, during fragmentation of meson strings the number of baryon-anti-baryon pairs produced can be either one or none.and that control parameter affects the relative probability of the two outcomes, becoming relevant to the output and not only to the execution time.
Choosing a high value for this P Q→QQ parameter can be viewed as getting closer to the second correction "always-baryon policy" even without imposing it.But increasing it causes the mean number of baryons in a stringfragmentation iteration to rise, and with it the possibility of discarding the fragmentation and having to rerun it due to exceeding the one-baryon policy.The effect on the global execution time could then be dramatic.Hence, the "always-baryon policy" can be viewed as an effective way of achieving the same effect without having to increase this parameter nor pay the extra execution time.
On the downside, this always-baryon policy is far from perfect from a computing implementation point of view.Compared to the one-baryon policy alone, it makes the code slower and more likely to crash.As it is, this implementation makes a combination with the afterburners method impractical.For completeness Fig. 6 shows, top row, the correlation function C from PYTHIA8 pp collisions, after applying the always-baryon policy, for pp, ΛΛ, Λp, and Λp particle pairs.Bottom row shows the correlation function azimuthal projections compared to ALICE data [2].Observe the excellent agreement between data and simulation for baryonbaryon correlations where the anticorrelation observed in data is perfectly reproduced by the simulation.The strong discrepancy in the baryon-antibaryon correlation function compared to data shown in Fig. 6 is not produced by the application of the always-baryon policy.The discrepancy is already observed from PYTHIA8 without introducing any modification as was also realized by the ALICE collaboration [2].

V. CONCLUSION AND PERSPECTIVES
We have examined the problem of baryon anticorrelations, a severe current disagreement between experimental data from the ALICE collaboration (lately also from the STAR collaboration) and theoretical Monte Carlo simulations of hadron collisions.Such anticorrelations are typical of lower-energy nuclear physics, for example in two-proton radiactivity (see [12], right panel of Fig. 18, that shows a decrease at cos θ Y = −1 when the protons are emitted back to back, respect to cos θ Y > −1, an effect ascribed to Coulomb repulsion and the Pauli principle) but are now also routinely obtained in high-energy collisions, due to the improved particle identification; they should, and are not, reproduced by current event generators.
We believe that the "fault" lies on the theory side and that we have made a contribution to the discussion by exposing that modification of the Monte Carlo event generator package PYTHIA 8 at the fragmentation level can qualitatively account for anticorrelation.In our not very refined approach, this has been achieved by means of a "one-baryon policy" applied to each string.Because, with standard PYTHIA parameters, this depresses the baryon to meson ratio (the total number of baryons produced is naturally smaller), a temporary correcting "all-baryon policy" has also been adopted, so that in our implementation given in the appendix, each string fragments into exactly one baryon (plus an unspecified number of mesons).
If this problem was to be reassessed in future work, one could try to instead of suppress baryons from the same string, give them a natal kick to separate them (or more selectively suppress them depending on their produced momentum).Another interesting theoretical extension would entail an extended probability distribution along the string, in order to mimic baryon or quark wave functions therealong.
Implementing a mechanism to take into account Fermi statistics in some way in PYTHIA is really difficult if we want to keep its computer efficiency, reasonably fit the data and give a convincing physical explanation.This is largely due to the fact that the Lund string model is mostly classical with point like particles while Fermi statistics is a quantum phenomenon based on wavefunctions.These two description are both essential, the first one for speed of execution and the second one for physical meaning.Thus, for the problem to be solved in a cleaner way, a translation between the two formalisms is necessary.For this, one must perhaps wait for a deeper understanding of the soft QCD processes that string fragmentation approximates.Conversely, finding a convenient procedure explaining major data features can give hints towards understanding of those deeper processes.
The LEP baryon correlation data could be reasonably fit by PYTHIA as is, given that the color string did form linking a back-to-back primary quark-antiquark pair; this means that baryons from the same string did not form positive correlations near ∆η 0 ∆ϕ in OPAL data, as they were somewhat randomized, with the string frame not too far from the laboratory frame.
At the LHC strings are however formed at various rapidities and azimuths, with a natal Lorentz boost.Because of that string boost, two baryons formed from the same string will create that positive correlation in the laboratory frame.Therefore, to avoid it and bring about the anticorrelation seen in the data, two-baryon production from the same string should be suppressed: our way of achieving it is the very rough pair of policies (one-baryon and all-baryon) that certainly need to be improved in future work.s t r i n g V e r t i c e s .c l e a r ( ) ; l e g M i n V e r t i c e s .c l e a r ( ) ; l e g M i d V e r t i c e s .c l e a r ( ) ; posEnd .hadSoFar = 0 ; negEnd .hadSoFar = 0 ; // For c l o s e d g l u o n s t r i n g : p i c k f i r s t b r e a k u p r e g i o n .system .setUp ( i P a r t o n , e v e n t ) ; stopMassNow = stopMass ; int n E x t r a J o i n = 0 ; // F a l l b a c k l o o p , when j o i n i n g i n t h e m i d d l e f a i l s .B a i l o u t i f s t u c k .f o r ( int iTry = 0 ; ; ++iTry ) { i f ( iTry > NTRYJOIN) { i n f o P t r −>errorMsg ( " E r r o r i n S t r i n g F r a g m e n t a t i o n : : fragment : " " s t u c k i n j o i n i n g " ) ; i f ( h a s J u n c t i o n ) ++n E x t r a J o i n ; i f ( n E x t r a J o i n > 0 ) e v e n t .popBack ( n E x t r a J o i n ) ; return f a l s e ; } // A f t e r s e v e r a l f a i l e d t r i e s j o i n some ( e x t r a ) n e a rb y p a r t o n s .i f ( iTry == NTRYJOIN / 3 ) n E x t r a J o i n = e x t r a J o i n ( 2 ., e v e n t ) ; i f ( iTry == 2 * NTRYJOIN / 3 ) n E x t r a J o i n += e x t r a J o i n ( 4 ., e v e n t ) ; // A f t e r s e v e r a l f a i l e d t r i e s g r a d u a l l y a l l o w l a r g e r s t o p mass .// Check how many n e ar b y s t r i n g p i e c e s t h e r e a r e f o r t h e n e x t hadron .double nNSP = ( c l o s e P a c k i n g ) ? n e a r S t r i n g P i e c e s ( nowEnd , r a p P a i r s ) : 0 .; // The FlavourRope t r e a t m e n t c h a n g e s t h e f r a g m e n t a t i o n p a r a m e t e r s .i f ( f l a v R o p e P t r ) { i f ( ! fl avR ope Ptr −>doChangeFragPar ( f l a v S e l P t r , z S e l P t r , pTSelPtr , ( fromPos ?hadMomPos .m2Calc ( ) : hadMomNeg .m2Calc ( ) ) , i P a r t o n , ( fromPos ?idPos : idNeg ) ) ) i n f o P t r −>errorMsg ( " E r r o r i n S t r i n g F r a g m e n t a t i o n : : fragment : " " FlavourRope f a i l e d t o change f r a g m e n t a t i o n p a r a m e t e r s ." ) ; } // P o s s i b i l i t y f o r a u s e r t o change t h e f r a g m e n t a t i o n p a r a m e t e r s .i f ( ( userHooksPtr != 0 ) && userHooksPtr−>canChangeFragPar ( ) ) { i f ( ! userHooksPtr−>doChangeFragPar ( f l a v S e l P t r , z S e l P t r , pTSelPtr , ( fromPos ?idPos : idNeg ) , ( fromPos ?hadMomPos .m2Calc ( ) : hadMomNeg .m2Calc ( ) ) , i P a r t o n , &nowEnd ) ) i n f o P t r −>errorMsg ( " E r r o r i n S t r i n g F r a g m e n t a t i o n : : fragment : " " f a i l e d t o change h a d r o n i s a t i o n p a r a m e t e r s ." ) ; } // Check w h e t h e r t o use s p e c i a l hard d i q u a r k h a n d l i n g i n beam remnant .

FIG. 2 :
FIG. 2: Changes to the prediction of the proton-proton R 2 correlation function induced by the "afterburners" correction.(a) Proton-proton correlation predicted directly from the PYTHIA Monte Carlo program, without correction.(b) Predicted proton-proton correlation corrected by the afterburners method.Note that the positive correlation at the surroundings of ∆ϕ 0 has not been erased.
FIG.3: Various two-particle baryon and/or antibaryon correlations extracted at the level of strings in PYTHIA.By meson string we mean one with zero baryon number, while a baryon string contains baryon number one (and thus, one net baryon is produced in the fragmentation).By S frag we mean that the two particles in each term of the correlation comes from the same string within the same event.

( a )
FIG.4: String-string S str correlations in pp collisions from PYTHIA showing clear ∆ϕ = π peaks that would eventually evolve to back-to-back jets.No kinematic cuts have been applied.As the strings carry a lot of energy, the pseudorapidity scale is not comparable to measured hadrons.Because high-η hadrons escape ALICE's acceptance, the contribution of the beams remnant from the original collided nuclei are erased in experimental data, but are clearly visible here in the uncut simulation, as a pair of ridges at |∆η| > 10.
(a) pp C correlation with unmodified PYTHIA (b) pp C correlation, one-baryon per string policy (c) pp C correlation, always-baryon policy Baryons produced (10 6 ) p pp pp Uncorrected PYTHIA 5.44 1.52 2.61 One-baryon policy 1.55 0.14 0.44 Always-baryon policy 5.20 1.32 2.28 3rd and 4th columns : number of same-event pairs.Less pairs entail a noisier figure.
FIG. 5: pp correlation function C in pp collisions from standard PYTHIA 8 (top left) shows clear forward and backward ridges for ∆ϕ = 0, π and all relative pseudorapidities.The one-baryon policy modification (top right) forbids two baryons from the same string, erasing the forward correlation but depressing the total number of baryons.The additional modification of the always-baryon policy (lower left) solves this problem by increasing the overall number of baryons.The backward correlation is still there; the forward correlation has now been turned into an anticorrelation; and the number of baryons produced is more adequate.
i s C l o s e d = c o l C o n f i g [ iSub ] .i s C l o s e d ; i f ( i s C l o s e d ) i P a r t o n = f i n d F i r s t R e g i o n ( iSub , c o l C o n f i g , e v e n t ) ; // For j u n c t i o n t o p o l o g y : f r a g m e n t o f f two o f t h e s t r i n g l e g s .// Then i P a r t o n o v e r w r i t t e n t o re mainin g l e g + l e f t o v e r d i q u a r k .pJunctionHadrons = 0 .; h a s J u n c t i o n = c o l C o n f i g [ iSub ] .h a s J u n c t i o n ; i f ( h a s J u n c t i o n && !fragmentToJunction ( e v e n t ) ) return f a l s e ; int j u n c t i o n H a d r o n s = hadrons .s i z e ( ) ; i f ( h a s J u n c t i o n ) { idPos = e v e n t [ i P a r t o n [ 0 ] ] .i d ( ) ; idNeg = e v e n t .back ( ) .i d ( ) ; pSum −= pJunctionHadrons ; } // S e t up k i n e m a t i c s o f s t r i n g e v o l u t i o n ( = motion ) .
bool forbidPopcornNow = f a l s e ; i f ( f o r b i d P o p c o r n && !h a s J u n c t i o n && ( nowEnd .hadSoFar == 0 ) ) { int iNow = ( fromPos ) ? i P o s : iNeg ; i f ( e v e n t [ iNow ] .i s D i q u a r k ( ) ) { bool motherInBeam = ( e v e n t [ iNow ] .s t a t u s A b s ( ) == 6 3 ) ; bool grannyInBeam = ( e v e n t [ e v e n t [ iNow ] .mother1 ( ) ] .s t a t u s A b s ( ) == 6 3 ) ; i f ( motherInBeam | | grannyInBeam ) forbidPopcornNow = true ; } } // C o n s t r u c t t r i a l hadron and c h e c k t h a t e n e r g y remains .nowEnd .newHadron (nNSP , forbidPopcornNow ) ; i f ( energyUsedUp ( fromPos ) ) break ; // O p t i o n a l l y a l l o w a hard baryon f r a g m e n t a t i o n i n beam remnant .bool useInputZ = f a l s e ; double zUse = 0 .5 ; i f ( forbidPopcornNow && hardRemn ) { useInputZ = true ; zUse = z S e l P t r −>zLund ( aRemn , bRemn ) ; } // C o n s t r u c t k i n e m a t i c s o f t h e new hadron and s t o r e i t .Vec4 pHad = nowEnd .k i n e m a t i c s H a d r o n ( system , s t r i n g V e r t i c e s , useInputZ , zUse ) ; int statusHad = ( fromPos ) ? 83 : 8 4 ; nowEnd .hadSoFar += 1 ; // Change s t a t u s code i f hadron from j u n c t i o n .
i f ( abs ( nowEnd .idHad ) > 1000 && abs ( nowEnd .idHad ) < 1 0 0 0 0 ) { i f ( fromPos && e v e n t [ i P o s ] .s t a t u s A b s ( ) == 74 && !usedPosJun ) { statusHad = 8 7 ; usedPosJun = true ; } i f ( ! fromPos && e v e n t [ iNeg ] .s t a t u s A b s ( ) == 74 && !usedNegJun ) { statusHad = 8 8 ; usedNegJun = true ; } i f ( ! fromPos && h a s J u n c t i o n && !usedNegJun ) { statusHad = 8 8 ; usedNegJun = true ; } } // P o s s i b i l i t y f o r a u s e r t o v e t o t h e hadron p r o d u c t i o n .i f ( ( userHooksPtr != 0 ) && userHooksPtr−>canChangeFragPar ( ) ) { // P r o v i d e f u l l p a r t i c l e i n f o f o r v e t o d e c i s i o n .i f ( userHooksPtr−>doVetoFragmentation ( P a r t i c l e ( nowEnd .idHad , statusHad , iPos , iNeg , 0 , 0 , 0 , 0 , pHad , nowEnd .mHad) , &nowEnd ) ) continue ; } // B oo kk ee pi ng o f momentum t a k e n away .i f ( fromPos ) hadMomPos += pHad ; e l s e hadMomNeg += pHad ; // Append produced hadron .int colHadOld = nowEnd .c o l O l d ; int colHadNew = nowEnd .colNew ; i f ( ! nowEnd .fromPos ) swap ( colHadOld , colHadNew ) ; hadrons .append ( nowEnd .idHad , statusHad , iPos , iNeg , 0 , 0 , colHadOld , colHadNew , pHad , nowEnd .mHad ) ; i f ( pHad .e ( ) < 0 .) break ; // Update s t r i n g end and remai ning momentum .nowEnd .update ( ) ; pRem −= pHad ;// End o f f r a g m e n t a t i o n l o o p .} // Check how many n ea r b y s t r i n g p i e c e s t h e r e a r e f o r t h e l a s t hadron .double nNSP = ( c l o s e P a c k i n g ) ? n e a r S t r i n g P i e c e s ( ( ( rndmPtr−>f l a t ( ) < 0 .5 ) ? posEnd : negEnd ) , r a p P a i r s ) : 0 .; // When done , j o i n i n t h e m i d d l e .I f t h i s works , t h e n r e a l l y done .i f ( f i n a l T w o ( fromPos , event , usedPosJun , usedNegJun , nNSP) ) { i f ( ! r e d o ) break ; // redo i f t h e number o f b a r y o n s i s not r i g h t i f ( iTry > 10 && t o l e r a n c e ) { i n f o P t r −>errorMsg ( " E r r o r i n S t r i n g F r a g m e n t a t i o n : : fragment : " " u s e t o l e r a n c e " ) ; break ; } e l s e { int nBAB = 0 ; // number o f b a r y o n s + a n t i b a r y o n s int idHad ; f o r ( int kount = 0 ; kount < hadrons .s i z e ( ) ; ++kount ) { idHad = hadrons [ kount ] .i d ( ) ; i f ( idHad > 1000 && idHad < 1 0 0 0 0 ) ++nBAB;