AMY Lorentz invariant parton cascade - the thermal equilibrium case

We introduce the parton cascade Alpaca , which evolves parton ensembles corresponding to single events according to the effective kinetic theory of QCD at high temperature formulated by Arnold, Moore and Yaffe by explicitly simulating elastic scattering, splitting and merging. By taking the ensemble average over many events the phase space density (as evolved by the Boltzmann equation) is recovered, but the parton cascade can go beyond the evolution of the mean because it can be turned into a complete event generator that produces fully exclusive final states including fluctuations and correlations. The parton cascade does not require the phase space density as input (except for the initial condition at the starting time). Rather, effective masses and temperature, which are functions of time and are defined as integrals over expressions involving the distribution function, are estimated in each event from just the parton ensemble of that event. We validate the framework by showing that ensembles sampled from a thermal distribution stay in thermal equilibrium even after running the simulation for a long time. This is a non-trivial result, because it requires all parts of the simulation to intertwine correctly.


Introduction
According to our current understanding of collisions of heavy nuclei at collider energies the high density in the final state of these collisions leads to strong final state re-scattering that rapidly drives the system towards local thermal equilibrium [1,2].At (proper) time O(1 fm/c), the systems enters an extended phase of hydrodynamic evolution, during which it undergoes continued strong longitudinal and moderate transverse expansion.During this phase, it is well described by viscous hydrodynamics [3][4][5] with a low shear viscosity to entropy density ratio η/s and equation of state of a deconfined quark-gluon medium [6].At the pseudo-critical temperature T c ≃ 160 MeV [7] the partonic medium freezes out into hadrons that continue to interact until the kinetic freeze-out around 150 MeV.Correlations due the collective expansion of the system are imprinted onto the hadronic final state and give rise to an azimuthal anisotropy (quantified by the flow coefficients v n ), a near-side long-range correlation in rapidity known as the ridge, and a hardening of transverse momentum spectra with ⟨p ⊥ ⟩ increasing with hadron mass.The comparatively high temperature at chemical freeze-out also leads to enhanced strangeness production [8][9][10][11][12][13].
The interpretation of these observations as signals of final state collectivity is, however, challenged by the observation of these signatures in proton-nucleus and high-multiplicity protonproton collisions.Naively, these systems are too small and too dilute to develop a sizeable thermalized medium.Still, the question has to be asked to what extent there is final state re-scattering also in these small systems and whether this could lead to a partial equilibration that is visible in hadron distributions (see [14,15] for reviews).
In the case of heavy-ion collisions, the rapid equilibration to an extent that viscous hydrodynamics becomes applicable can be understood both in terms of strong-coupling dynamics [16][17][18][19], as well as in terms of weak-coupling dynamics [20][21][22].In the latter case it can be studied using the effective kinetic theory of QCD by Arnold, Moore and Yaffe (AMY) [23][24][25], describing the dynamics of quarks of gluons.Kinetic theory is a natural candidate for the description of small collision systems [26][27][28][29], since it is valid in-and out-of-equilibrium.
Progress in understanding similarities and differences between small and large collision systems will depend critically on detailed apples-to-apples comparisons between theory and experimental data.This is particularly important for the small systems, where the interpretation of data is complicated by biases and auto-correlations.Monte Carlo (MC) event generators are tools that can bridge the gap between theory and experiments and are indispensable for a thorough understanding of collider data in all areas.Compared to the perturbative regime of proton-proton physics, where sophisticated event generators built on the Standard Model (and beyond) are available, event generators for soft physics and heavy-ion physics are far less advanced.It is our aim to fill in this gap by developing a MC event generator representation of the AMY effective kinetic theory that is applicable to small and large collision systems with a focus on soft particle production.Alpaca (AMY Lorentz invariant PArton CAscade) is implemented as a module in the Sherpa [30] event generator.The main difference to existing parton cascade codes [31][32][33][34][35][36][37][38][39][40] is that it reproduces exactly the AMY kernels for all elastic scattering, splitting and merging processes.Furthermore, it is implemented in a way that is Lorentz-invariant by construction.This avoids problems with residual causality violation [41] and altered correlations and fluctuations due to particle subdivision [42].In this paper we focus on the implementation of the effective kinetic theory in a parton cascade code and show that it reproduces AMY in thermal equilibrium.A major challenge in this is to extract the screening masses and effective temperature T * , which are non-local in phase-space, from the parton ensemble.For the validation of the thermal equilibrium we initialise the system in thermal equilibrium and convince ourselves that the distribution function does not change.This will only be the case when detailed balance is respected.In particular, splitting and merging rates have to be the same, which requires a non-trivial interplay between elastic and inelastic processes.The scattering rates depend on the screening masses and the splitting and merging rates on the effective temperature, detailed balance is therefore a critical test of our dynamical extraction of these quantities.Running the simulation for a long time in thermal equilibrium is an important validation exercise showing that the dynamics are correctly encoded.The next step is then to study how a system that is initialised out of equilibrium approaches equilibrium.Since this is a different question we will discuss it in a separate publication.
The paper is organised as follows: in Section 2 we briefly review the method used to make the parton cascade Lorentz invariant and Section 3 gives a summary of the AMY effective kinetic theory.Details on the implementation in Alpaca are given in Section 4, with an introduction to the framework in 4.1, a discussion of the set-up used for the thermal equilibrium simulations in 4.2, details on elastic scattering in 4.3, splitting and merging in 4.4 and lastly the results of a long run with all the different components included is shown in 4.5.We conclude the discussion in Section 5 and give further technical details in the appendices.

A Lorentz invariant parton cascade
The no-interaction theorem [43] states that N particles moving in a 6N dimensional phase space cannot interact if Lorentz invariance (instead of Galilei invariance) is required.This theorem can, however, be circumvented by regarding the particles as moving in a 8N dimensional phase space [44], i.e. formulating the theory in terms of four-positions and four-momenta of the particles.With a physical choice of the quasi-potential describing the interactions the particles are classically off-shell (i.e.p2 i ̸ = m 2 i ) when they are within the range of the quasi-potential and on-shell otherwise.The four-vectors can be parameterised by an evolution parameter τ , that is a Lorentz scalar.The equations of motion are then given by where H is the classical Hamiltonian.
We here follow the approach of [37] and use a simplified Hamiltonian, in which there are only binary interactions of particles taking place at discrete points in τ 1 .The interactions thus lead to instantaneous (in τ ) changes of the momenta and between interactions the particles travel on free trajectories.We thus require only the free Hamiltonian, which for a system of N particles is the sum of N single particle free Hamiltonians.In Alpaca we neglect the (parametrically small) effective masses and treat all partons as massless.The Hamiltonian is then given by where λ is an arbitrary Lagrange multiplier [46].The equations of motion are easily solved to yield x iµ (τ ) = 2λp iµ (τ − τ 0i ) + x iµ (τ 0i ) .
The parameter τ does not have a physical interpretation.One can, however, chose an arbitrary but fixed reference frame and calculate the particles' momentum and position in this frame.Using x iµ = (ct, r i ) and p iµ = (ω i /c, k i ) with |k i | = ω i /c 1 A similar way of formulating a Lorentz invariant parton cascade has recently been proposed in [45].
In order to find the {τ 0i } at which the initial conditions {x i (τ 0i )} are specified, the parameter τ has to be related to some observer time.This can be done by generalising the procedure introduced in [44] for a 2-particle system to N particles 2 .The position of the centre-of-momentum X µ in the centre-of-momentum frame satisfies X c.m.
In order to decide whether a given pair of particles interacts, one needs to measure their distance.Following [37], the Lorentz invariant distance squared between two particles is defined as where x = x i − x j is the relative four-distance between the particles, and p = p i + p j is the total four-momentum of the pair.In the centreof-momentum frame of the pair, this reduces to the squared three-distance.When the distance becomes smaller than σ scat /π, where σ scat is the scattering cross section, an interaction occurs3 .
The interactions are ordered in τ in a frameindependent way.This avoids the problem of causality violation in parton cascades that appears when ordering interactions in observer time, in which case the ordering is frame dependent.
In general, when a pair of particles interacts, the condition d ij < σ scat /π is satisfied during a finite τ interval and one has to make a choice at which point in τ the interaction takes place.We here choose the point of closest approach between the particles, i.e. the τ corresponding to min(d ij (τ )).The value of τ for which the closest approach is reached can be calculated from Eq. (10) and is given by The parton cascade proceeds by moving from one scattering to the next.The 1 → 2 splitting processes are easily included by translating the splitting rate in terms of an observer time into a splitting rate in τ using Eq. ( 7) (see Section 4.4.1).

The AMY Effective Kinetic Theory
In [23], P. Arnold, G. Moore, and L. Yaffe derived an effective kinetic theory of QCD at high temperatures.At sufficiently high T , where the QCD coupling g(T ) is small, the temperature T sets a hard scale, while gT is a softer scale.It is a starting assumption of the effective-theory setup that these parametrically separated scales are also quantitatively well-separated.Under this assumption the QCD plasma consists of quarks and gluons as welldefined quasi-particles with typical momenta of order T .These hard partons propagate as nearly free and nearly massless particles, since thermal masses are O(gT ) and the small-angle scattering rate is O(g 2 T ).The rate of large-angle scattering is O(g4 T ).
In addition to elastic 2 ↔ 2 scattering, 1 ↔ 2 splitting and merging processes have to be taken into account.In vacuum, these are kinematically allowed only when all three momenta are exactly collinear.In plasma, this restriction disappears the scattering probability, but did not observe significant differences.
when a soft scattering is involved.The process is then nearly collinear with a rate O(g 4 T ), i.e., comparable to the large-angle elastic-scattering rate.The process is not instantaneous, as it involves an off-shell propagator.The formation time of splitting/merging processes turns out to be O(1/g 2 T ), i.e., the same as the small-anglescattering mean free path.Consequently, additional soft collisions occur during the formation time making it an N + 1 ↔ N + 2 process.The additional soft collisions are not resolved individually by the radiation and cannot be treated as classically independent.Instead, they act coherently giving rise to the QCD analogue of the Landau-Pomeranchuk-Migdal (LPM) effect [47,48].The splitting/merging processes involving N coherent soft scatterings are denoted as "1 ↔ 2" processes.
The effective kinetic theory is formulated as a set of Boltzmann equations ) for the phase space densities f s (x, p, t) of partons of species s at position x with momentum p at time t.The collision kernels C 2↔2 s and C "1↔2" s encode the effect of elastic scattering and splitting/merging, respectively.
The elastic scattering kernel includes all 2 ↔ 2 processes, because all of them have parametrically the same hard scattering rate.It has the usual form where the dependence of the phase space densities f s on position and time has not been written out.
Capital letters denote four-vectors and the particles are taken as massless.The factor ν s is the number of spin and colour states (i.e.ν q = 6 for quarks and anti-quarks and ν g = 16 for gluons).
The effective in-medium matrix elements squared |M ab cd | 2 are taken as summed over the final state and averaged over the initial state 4 .It should be noted that due to screening in the plasma the matrix elements depend on the phase space densities.The short hand notation for the Lorentz-invariant momentum integration has been used, and the factor in front of the integral contains a factor of 1/2 for symmetry.The factors (1 ± f ), where the upper sign is for gluons and the lower is for (anti)quarks, encode Bose enhancement or Pauli blocking, depending on the species.The first term in Eq. ( 13) is a gain term accounting for partons arriving at momentum p due to elastic scattering, while the second is the loss term accounting for partons leaving momentum p.
The kernels for the effective splitting/merging processes can be written as where the (small) transverse momenta associated with the splitting or merging have been integrated out neglecting the deviations from exact collinearity in the phase space densities.The functions γ a bc are the differential splitting/merging rates for the a ↔ b + c processes and contain the effect of multiple soft scattering during the formation time.As in the elastic scattering case, they depend implicitly on the phase space densities.The rates γ a bc should be understood as summed (not averaged) over all external partons.
The effective kinetic theory is leading-order accurate in the coupling λ = g 2 N c = 4πα s N c for hard partons.
In order for any kinetic theory to be applicable the duration of scattering events must be small compared to the mean free path.In the case of the AMY effective kinetic theory in equilibrium this in ensured by the assumption that the temperature is sufficiently high for the theory to be weakly coupled.Out of equilibrium, the question of applicability of the theory requires a more careful discussion, which is provided in [23].The most important requirements are that the phase space densities support a separation of scales as outlined at the beginning of this section.This means that all relevant excitations should have momenta that are large compared to the screening masses.The effective masses, in turn, have to be much larger than all other mass scales (Λ QCD and quark masses).Furthermore, the phase space densities have to be sufficiently slowly varying, in particular they should not vary significantly over the length/time of the formation time of near the collinear splitting/merging processes.Lastly, the formulation of the AMY kinetic theory is technically leading-order accurate only for isotropic distributions due to certain plasma instabilities.However, numerical evidence points that the effect of the unstable modes may be quantitatively small [49].4 The Alpaca framework

General remarks
Alpaca (AMY Lorentz invariant PArton CAscade) is a parton cascade which solves the Boltzmann equation of the AMY effective kinetic theory indirectly by, in a Lorentz invariant way, evolving a parton ensemble according to the AMY collision kernels.To fully capture the dynamics of these kernels, three types of processes are considered: elastic scattering, quasi-collinear merging and quasi-collinear splitting.
In addition to being able to evolve finitesized systems, Alpaca also contains an implementation of infinitely sized ensembles in thermal equilibrium, used to verify simulations in a simpler setting with analytically known results.The infinite-sized system is implemented through a box with periodic boundary conditions, something that introduces further complications to the Lorentz invariant framework, which is discussed in detail in Section 4.2.
The evolution of the particle ensemble takes place in τ as defined in Section 2. Each event starts at τ = 0 and is set to end at some τ max .Since τ does not have a physical meaning it is convenient to also define an observer time t max that the system should reach.To ensure that all particles are evolved to t max the simulation will keep running until all particles have t ≥ t max even when this requires going to τ > τ max .There is no need for t max and τ max to be the same, but we here choose t max = τ max .After the initialization of an event, the next possible processes of the three described above are found and ordered in τ .Between any processes taking place partons are treated classically as point particles and are considered free streaming, according to Eq. ( 5).
The possible interactions of particles implemented in Alpaca to fully capture the dynamics given by AMY are elastic scattering and inelastic splitting/merging.However, due to the LPM effect present in the inelastic processes as well as the restriction of collinearity, the inelastic and elastic processes are treated differently.For the elastic case, an effective cross section extracted from the elastic collision kernel is used.In the current implementation we use a black disk approximation for the cross section, i.e., an interaction between a pair of particles occurs if their invariant distance at closest approach satisfies d ij < σ/π.If several processes are allowed, σ is the sum of the cross sections of the individual processes.It is then decided in a second step which processes is going to take place according to the relative contributions to the total cross section.
In the Lorentz invariant framework a problem arises in a box with periodic boundary conditions.As discussed in Section 4.2 the algorithm for finding the closest approaches between pairs of particles requires the cross section as input.The cross section however, depends on the local quantities m g/q (x, p) and f (x, p) which are not known a priori.The solution is to find an overestimate of the cross sections that is independent of local quantities and then reject interactions with a probability that is given by the true cross section (evaluated at the point of the potential interaction) divided by the overestimate 5  For the inelastic processes, the effect of multiple coherent scattering during the formation time (LPM effect) is already included in the splitting/merging rates γ a bc .We therefore don't explicitly simulate them, but treat the inelastic processes as effective 1 ↔ 2 processes.However, since we are dealing with an ensemble of discrete particles, a pair of particles will never be exactly collinear.We therefore have to allow for a small relative transverse momentum k ⊥ for the merging processes, and consequently also for splitting (the two types of processes have to cover the same phase space in order to preserve detailed balance).The k ⊥ is absorbed by a nearby parton, which we call the 'recoil parton'.This recoil parton thus effectively acts as the (possibly multiple) parton(s) with which the soft coherent scattering occurs during the formation time of the splitting/merging.Since merging processes initially consists of a pair they can be treated similarly to elastic scattering, where we extract an effective cross section from the inelastic collision kernels and then evaluate this using a black disk approach.
For collinear splittings, an implementation of the Veto algorithm is utilized, see Section 4.4.1 for details.Since the splitting and merging rate is infra-red divergent we have to introduce a small cut-off in x.We do this by introducing a minimum momentum for all particles, p min , so that no particle is initialized with p < p min or allowed to scatter/split/merge into a momentum lower than this.This does not compromise the theory's accuracy.
Currently, Alpaca is only utilizing the general event handling framework in Sherpa, due to the focus on testing the model in thermal equilibrium.We plan to integrate Alpaca in the event generation in Sherpa in the future to allow for hadronisation and inclusion of hard processes like jet production.

Simulating a system in equilibrium
In numerical studies of systems of particles in thermal equilibrium, one usually considers the system in a box with periodic boundary conditions.The position of a particle leaving the box at one end is then shifted such that the particle re-enters the box at the other end (the momentum remains unchanged).The justification for this procedure is that a system in thermal equilibrium is homogeneous and therefore, if one imagines the space to be divided into equally sized cells, all cells look the same on average.The periodic boundary condition makes all cells look exactly the same.
To simulate the evolution of a system of particles, one has to consider the interactions of each particle with all other particles.Since all cells are identical, it is enough to consider all interactions inside one cell.One therefore has to identify all pairs of particles that interact inside the box within a pre-defined time interval.In a parton cascade considering only binary interactions, one thus has to compute the time and distance of closest approach for each pair of particles.The distance of closest approach decides whether the pair will scatter, and the time at which this happens decides whether the scattering will take place inside the box (or whether one or both particles will have left the box before the encounter).To correctly account for all possible interaction inside the box, one has to consider scattering of a particle inside the box with particles coming from neighbouring cells.Both the time and distance of closest approach depend only on the relative distance between the two particles.One thus captures all possible encounters of particles 1 and 2 and all their copies in other cells by considering the copy of particle 1 located in the box and all copies of particle 2 on other cells.
When the system with periodic boundary conditions evolves in an observer time, e.g. the time in the rest frame of the box, a particle in the box can only collide inside the box with a particle coming from inside the box or one of the neighbouring 26 cells6 , which are copies of the original cell, surrounding the box.If the second particle comes from farther away, the closest approach of the pair will always happen outside the box.Such an evolution is, however, not Lorentz invariant, but frame dependent.As discussed in Section 2 the evolution of a multi-particle system can be simulated in a frame independent way by ordering scattering events not in an observer time t, but in a generalised time τ , which is a Lorentz scalar.In this way the evolution of the system becomes Lorentz invariant, but now a particle in the box can scatter inside the box with a copy of another particle coming from any cell.One therefore has to consider a countably infinite number of encounters for this pair.Of these we are only interested in those that bring the particles close enough for a scattering.And of those that lead to a scattering we need to find the first one, i.e. the one with the smallest τ starting from the current τ of the evolution.
We have developed an algorithm that finds the earliest closest approach for a given pair of particles that is close enough for an interaction to occur.It systematically checks closest approaches between the first particle in the box with copies of the second particle coming from other cells ordered in increasing invariant time τ .This algorithm needs the scattering cross section as input, because when a closest approach is not close enough for a scattering to occur it has to continue checking other copies of the second particles.The algorithm is specified in full detail in Appendix A.

General considerations
As mentioned in Section 4.1, to determine if a parton pair scatters elastically, and if so which type of elastic scattering occurs, multichannel accept/reject Monte Carlo sampling is used.The scattering matrix elements are related to the elastic cross section through where is the phase-space measure of the outgoing particles.Note that we pick up a factor of 1/2 in the first equality due to either gluons being indistinguishable or from double counting of fermionic processes.The factor within the bracket can be expressed using the relative angle θ ∆pk between p and k as s/2|p||k| = 1 − cos(θ ∆pk ) and enforces a vanishing contribution for particles moving collinearly.In an isotropic system this expression integrates to unity.For hard scattering screening effects are negligible and the standard vacuum matrix elements can be used.For soft scattering (small t or u), however, the matrix elements receive large corrections from in-medium physics.Following the strategy of [50,51] we replace the infra-red divergent 1/t 2 terms in the matrix elements by where m s is the effective mass of the parton in the propagator and the prefactors ζ g = e 5/3 /4 and ζ q = e 2 /4 contain the conversion factor between (asymptotic) effective mass and screening mass and a scaling factor that ensures that the modified matrix elements reproduce to leading order the HTL results for drag and momentum broadening [50,52].Strictly speaking, the factors ζ s should be re-calculated for our choice of making the replacement, which differs from [50,51].We leave this for the future.We list the matrix elements used in Alpaca in Appendix C.1.
In addition to the scattering matrix elements the collision kernels contain factors of (1 ± f ) for the outgoing partons, which encode Bose enhancement and Pauli blocking.Thus the integral over the outgoing momenta of the matrix elements does not give the standard cross sections, but effective cross sections containing the Bose enhancement/Pauli blocking factors.The effective cross sections can be used for simulating the scatterings in the standard way.However, also the Bose enhancement/Pauli blocking factors depend on the position of the scattering and are not known beforehand.We therefore include them in the rejection step by including an overestimate of these factors when integrating the effective cross sections and later rejecting the scattering with the ratio of the true Bose enhancement/Pauli blocking factor and the overestimate.
In practice, if the closest approach occurs at some τ and the distance is smaller than σ/π where σ is the overestimated effective cross section, m 2 s (x, p) and f s (x, p) are extracted as described in Section 4.3.2.In combination with the matrix element, as given in Appendix C, the full cross section is found and the process is accepted or rejected based on the ratio σ ab 2→2 /σ ab 2→2 .If the process is accepted, t is sampled from the differential cross section (using accept/reject Monte Carlo with the overestimates of the matrix elements given in Appendix C) and the particles' momenta are updated.

The effective masses and phase space densities
The screening mass m s for a particle of species s, represents the medium screening effects and contributes in a highly non-trivial way to the evolution of the parton ensemble.It is responsible for regulating the otherwise divergent matrix elements in a local way and is in turn essential to capture interesting dynamics in the non-equilibrium case.The (position averaged) definition of the effective mass for a gluon is ) and for a fermion of species s where the average is over a spatial volume V .Since d 3 p/|p| and d 3 x/V are Lorentz invariant measures, and f (p, x) is a Lorentz invariant scalar [53,54], it follows that the screening masses are Lorentz invariant quantities.The screening masses have to be calculated for each scattering and at fixed observer time by integrating/summing over nearby partons.In any reference frame the two colliding partons will in general experience the scattering at different observer times.We choose to calculate the screening mass at the average of these times, which we denote by t D .A more serious complication is that when evolving the system in τ (rather than observer time) the partons' positions and momenta are known at fixed τ , which corresponds to a different observer time for every parton.Partons that are known at observer times later than t D for some scattering can be propagated back in time, but for partons with times earlier than t D this is not possible.We therefore choose to let the earlier partons free-stream to t D .This violates causality, as we are ignoring potential scatterings during that time which would change the partons' positions and momenta.This can be mitigated by solving the time evolution iteratively, as discussed in Section 4.5.
To extract m 2 s in a local way we have the following prescription.For any particle pair with their closest approach at some τ , we define t D to be the average observer time of the particle pair at τ .We then identify the N inc closest (in terms of Euclidean spatial distance in the local rest frame) particles, w.r.t. to the midpoint of the original particle pair at t D .This gives a volume V = 4πr 3 max /3, where r max is the Euclidian distance to the particle furthest away.Due to the semi-classical approach of Alpaca the partons are treated as pointlike and hence the density of particles can be expressed as which gives an expression for f (x, p).For a gluon propagator, following Eq.( 19), each of the N inc included particles then adds a factor of to m 2 g depending on the species s of the included particle.Similarly, following Eq.( 20) for a quark propagator, each surrounding particle adds a factor to m 2 q/q depending on whether the included particle is a gluon or quark (of the same species as the virtual quark).
In order to evaluate the Bose enhancement and Pauli blocking factors, 1 ± f a (x, p), that appear in the AMY collision kernels, we also need explicit approximations of f a (x, p).Consequently, we approximate the phase space density by once again treating the particles as pointlike and looking at the phase space volume V = V x V p where we have a spherical spatial volume around x, the point of interest, and a spherical shell around the origin containing the (absolute) momentum p in momentum space, i.e.
(24) Here, r is the spatial distance between the included particle and x and is set to some fixed value, while ∆p is the allowed differences in absolute momentum p between the original and included particles, also set to some fixed value.We then simply count the N inc particles that exist within V and the extracted phase space density follows through the relation given in Eq. ( 21), We also set an upper bound for the extracted phase space density, f a (p) = min(f a (p), T * /p), since the underlying theory is not valid for f (p) > T * /p.Note that the implementation above is set up for an isotropic distribution, which our test case of thermal equilibrium described later will be.For anisotropic distributions the method is easily generalized by restricting the volume in momentum space to a sphere around p.
The phase space density also has to be extracted at fixed observer time and the same freestreaming procedure as for the effective masses is used (with the same caveat concerning causality violation applying, although less severely since f is local in momentum space and therefore the spread in times is much smaller, cf.Eq. ( 9)).

The thermal equilibrium case
In order to verify that Alpaca reproduces the dynamics of AMY we put our particles in a box with periodic boundary conditions and sample our initial distributions in two different scenarios.First, the Boltzmann distribution which has the benefit of providing us with a simple setup that allows for many quantities to be found analytically and compared to, e.g. the scattering rate.This distribution does not give us the full picture though, since it does not factor in quantum effects like Bose enhancement and Pauli blocking.For some cases these effects are not of importance in our verification of Alpaca, and we can just exclude them in the simulations.For splitting and merging however, this is not possible 7 .Hence, we will also study the scenario of Bose-Einstein and Fermi-Dirac distributions, and To tackle the divergence for vanishing momentum in f BE (and later also divergent splitting/merging rates for low energy x) we introduce a cutoff in momentum, p min , so that no particle is initialized with momentum lower than this, and no particle is allowed to scatter/split/merge into a state with lower absolute momentum than this.
Both of the infinite thermalized systems above provide simplified scenarios where one can find analytical (and simpler numerical) results to compare against.To validate our implementation we show the following: (i) The elastic scattering rate dN coll./dt in the case of thermal equilibrium with no Bose enhancement and Pauli blocking, and with a fixed σ 2↔2 , is consistent with the expected results, both for the classical and Bose-Einstein and Fermi-Dirac case.(ii) The elastic scattering rate dN coll./dt with a non-constant σ 2↔2 (implemented as described in Section 4.3.1)produces the correct rate w.r.t. the elastic collision kernel in AMY EKT.This is also checked both for the classical and Bose-Einstein and Fermi-Dirac case.
(iii) The screening masses m 2 s and phase space densities f s are extracted correctly, i.e., they match the analytical values, both for the classical and Bose-Einstein and Fermi-Dirac case.
(i) The number of particles of species s, dN s , contained in an infinitesimal phase-space volume d 6 ξ = d 3 xd 3 p around x and p can be related to the phase-space density f s (x, p) through The collision kernels C s [f ] are related to the rate of change of the number of particles dN s in d 6 ξ between t and t + dt, i.e.
with C s [f ] + and C s [f ] − corresponding to the gain and loss terms in the collision kernel.The gain (or loss) term of the collision kernel can now be related to the rate of total number of scatterings as where we pick up a factor of 1/2 due to either integrating over both incoming momentum of indistinguishable particles when s is a gluon, or double counting processes in the sum when s is a quark or antiquark.Here we also define the momentum integrated gain term of the collision kernel for a species s as I ± s (with + and − for the gain and loss part of the collision kernel respectively), which are quantities of interest when evaluating (ii).Considering the case of a purely gluonic system without any Bose-enhancement we can extract a simple scaling behaviour for the rate of total number of collisions, where we have used the definition of the total cross section given in Eq. ( 16).The expression within the brackets integrates to unity in isotropic systems, as explained in Section 4.3.1, and so for the thermal case with a constant cross section we have the scaling behaviour with N g being the total number of gluons in the system, and n g the gluon density.From the equations above we also note that The scaling behaviour shown above has been examined in Alpaca through scanning the parameter space, and exhibits the correct scaling behaviour, as shown in Fig. 1.
(ii) Including quarks and/or Bose enhancement/Pauli blocking is straight-forward but the collision rate does then not exhibit the same simple scaling behaviour.However, due to the (approximate) time independence of f (x, p) in thermal equilibrium we can integrate ] numerically, using Monte Carlo methods with stratified sampling, and compare to the results of Alpaca through the relation found in Eq. ( 29).This numerical integration has been crosschecked against the code used in e.g.[28].The results of the comparison of the full numerical integration are shown in Tables 1 and 2 for Boltzmann and Bose-Einstein/Fermi-Dirac distributions, respectively.A comparison of the rates of each sub-process (for the loss-terms) are shown in Tables 3 and 4, again for the different distributions.The explicit expression of the integrated collision kernels are shown in Appendix C.2.As can be seen in the tables, the rates coincides well for f C , while for f BE/FD they are within ∼ 10% of the numerical results.This can be attributed to the cutoff p min which effectively changes the value of m 2 g/q that we extract in Alpaca, while this shift in screening mass is not included in the numerical results.
(iii) Lastly, we compare the implementation of m 2 s (x, p), described in Section 4.3.2 to the analytical expectations for thermal equilibrium, which for a system with both quarks and gluons are and These comparisons are shown in Fig. 2 and 3.The numerical values are extracted directly after initialisation when the partons all have the same observer time and the complication due to causality violation mentioned above is absent.As shown in Fig. 3, the mean value converges to the analytical expectation for both quarks and gluons, however, it requires a rather large number of particles to be included in the calculation.Approximating partons using Gaussian distributions (see Section 4.4.3)instead of Dirac delta functions might improve on this convergence rate.In Fig. 4 the scaling behaviour of the extracted m 2 s against different temperatures can be seen.
Finally, a comparison of f (x, p) can be seen in Fig. 5, which is shown to converge well to the analytical expectation.
All runs in this section are initialized in a box with volume V x = 7.6 3 GeV −3 and α s = 0.04.
Table 1 The full integrals of the elastic collision kernels 123)-( 126), are compared between their values produced by Alpaca (through the collision rate, see Eq. ( 32)) and by numerically integrating the collision kernels using Monte Carlo methods.The comparison is preformed in the thermal case with T = 0.4 GeV, αs = 0.3, τmax = 1 GeV −1 and box volume Vx = 10.1 3 GeV −3 .The errors of the numerical results are negligible compared to the results from Alpaca, and hence not presented.Table 2 The full integrals of the elastic collision kernels

Integral [GeV
q/q [f BE/FD ] (I ± q ), see Eqs. ( 123)-(126), are compared between their values produced by Alpaca (through the collision rate, see Eq. ( 32)) and by numerically integrating the collision kernels using Monte Carlo methods.The comparison is preformed in the thermal case, including quantum corrections in the initial distribution, with T = 0.4 GeV, p min = 0.1 GeV, αs = 0.04, τmax = 1 GeV −1 and box volume Vx = 7.6 3 GeV −3 .The errors of the numerical results are negligible compared to the results from Alpaca, and hence not presented.

Quasi-collinear splitting and merging processes
The effective "1 ↔ 2" splitting/merging rates have been obtained to leading order in [55] in terms of solutions to a linear integral equation.
An explicit analytical solution has been found in [56] for very energetic partons, where the LPM effect dominates the dynamics.In the incoherent Bethe-Heitler limit it is also possible to find an analytical solution.Between these two limiting cases no analytical solution is known.We here Table 3 The full integrals of the different subprocesses in the loss terms of the elastic collision kernels 123) and ( 125), are compared between their values produced by Alpaca (through the collision rate of each subprocess, see Eq. ( 32)) and by numerically integrating the collision kernels using Monte Carlo methods.The comparison is preformed in the thermal case with T = 0.4 GeV, αs = 0.3, τmax = 1 GeV −1 and box volume Vx = 10.1 3 GeV −3 .The errors of the numerical results are negligible compared to the results from Alpaca, and hence not presented.4.6 • 10 −5 0.9(1) I −,q q↔gg q 0.08(4) Table 4 The full integrals of the different subprocesses in the loss terms of the elastic collision kernels q/q [f BE/FD ] (I − q ), see Eqs. ( 123) and (125), are compared between their values produced by Alpaca (through the collision rate of each subprocess, see Eq. ( 32)) and by numerically integrating the collision kernels using Monte Carlo methods.The comparison is preformed in the thermal, including quantum corrections in the initial distribution, with T = 0.4 GeV, p min = 0.1 GeV, αs = 0.04, τmax = 1 GeV −1 and box volume Vx = 7.6 3 GeV −3 .The errors of the numerical results are negligible compared to the results from Alpaca, and hence not presented.
use an interpolation that reproduces the analytical results in the Bethe-Heitler and LPM regimes and agrees well with numerical solutions between them.Details are given in Appendix B.
As mentioned in Section 3, the splitting/merging rates describe the splitting (merging) of a (two) parton(s) in conjunction with an arbitrary number of coherent soft scatterings.They thus depend on the density of partons generating soft fluctuations of the background gauge field.It turns out that the phase space densities enter via a single parameter T * , which is the temperature of a system in thermal equilibrium that has the same infra-red behaviour as the system under consideration.In thermal equilibrium it thus reduces to the temperature.The extraction of T * in Alpaca is discussed in Section 4.4.3.

Splitting
Soft scatterings within the medium can bring partons slightly off-shell which kinematically allows for them to split nearly collinearly into two new partons.The average formation time of this process coincides with the mean free time for soft elastic collisions within the medium and so all soft scatterings occurring during the formation time of the splitting have to be considered coherently.Destructive interference leads to a suppression of the splitting rate compared to the incoherent case, which is known as the LPM effect in QCD.Probability distribution f BE/FD : m 2 g/q in thermal equilibrium, T = 0.4 GeV, N inc = 30 ALPACA: mean(m 2 g ) = 0.0604(9) GeV 2 Analytical: m 2 g = 0.0603 GeV 2 ALPACA: mean(m 2 q ) = 0.0268(3) GeV 2 Analytical: m 2 q = 0.0268 GeV 2 Fig. 2 Probability distribution of extracted m 2 g/q in thermal equilibrium, for T = 0.4 GeV, and 30 included particles.The errors presented for the mean values correspond to the error of the mean of the distribution (in contrast, the width of the distribution is caused by including a finite number of particles in the extraction and does not translate into an uncertainty on the mean).Left: f C .Right: f BE/FD .
We here assume that subsequent splittings are independent and the time evolution of a parton undergoing repeated splitting can be formulated as a Markov chain 8 .This is a well studied problem in event generators, and a common approach to sample the occurrence of particle branchings is the Veto Algorithm.This algorithm can be employed when the splitting probability is of the form for a particle of species a, where x (and 1 − x) correspond to the fractions of energy going to the two new particles.Here, P a (t, x) is the probability for the particle to split at time t with energy fraction x, N a (t) is the probability that a particle of species f BE/FD : m 2 g/q vs T, p min = 0.01 GeV, N inc = 30 Analytical: m 2 g = 3παsT 2 Analytical: m 2 q = 4παsT 2 /3 ALPACA: mean(m 2 g ) ALPACA: mean(m 2 q ) Fig. 4 Scaling behaviour of the mean value of extracted m 2 g/q in thermal systems, w.r.t different temperatures of the system.Number of included particles in each calculation is fixed to 30.Left: f C .Right: f BE/FD .a has not split yet at t, also known as the Sudakov factor, and h a (t, x) is the instantaneous splitting probability (density) at t and x.We solve the splittings in Alpaca by using an overestimate version of the Veto Algorithm.This means the following: We find an overestimate ha (t, x) ≥ h a (t, x) of the instantaneous splitting probability and then use ha (t) = ha (t, x)dx to sample t split .The energy fraction x is then sampled uniformly between x min and x max .The draw is accepted or rejected based on the ratio h a (t, x)/ ha (t, x).If the splitting is rejected we do not reset in time, but set the previously sampled t split to be our new initial time, and redo the sampling of splitting time and energy fraction until we find a splitting that is accepted.When multiple decay channels are available, t split is sampled from the sum of the corresponding splitting probabilities, and the correct decay channel is selected based on the ratio of the splitting probabilities, once the splitting has been accepted.To implement this algorithm, two quantities are needed: h a (t, x) and an overestimate ha (t, x) ≥ h a (t, x).We have previously defined the number of particles of species a in a phase space volume d 6 ξ at time dt in Eq. ( 29), and the number of particles of species a leaving d 6 ξ during dt in Eq. (31).Hence, the (total) instantaneous splitting probability for a particle in d 6 ξ of species a during time dt is The term in the splitting/merging kernel, Eq. ( 15), which corresponds to a particle of species a with momentum p splitting to two particles with momentum p ′ and k ′ in terms of the quasicollinear branching rate γ a bc is and so The instantaneous splitting probability as a function of x is then extracted from Eq. ( 41) as In Appendix B a detailed account of how to calculate γ a bc is given, and an overestimate of the instantaneous splitting probability given above can be found in Appendix C.3.Note that the overestimates of h a (x) depend on T * , which we assume does not change significantly over our relevant timescales.Hence, for the overestimates we sample T * at current τ (again taking the system to constant observer time using the same procedure as for the effective masses) and multiply with some fixed numerical factor c > 1 such that our overestimate effective temperature is cT * .The effective temperature is then re-calculated at τ split when evaluating h a (x).
The quantities h a (x) and ha (x) are used in the Veto Algorithm to sample x and t split , as described in the beginning of this section, and this gives us a corresponding τ split .Once a τ split is found, a formation time t form is calculated according to the corresponding expression given in [23], A ∆τ interval of length τ form corresponding to t form is placed around τ split such that τ split is distributed uniformly in the interval.During this formation time the particle cannot scatter elastically or quasi-collinearly merge because elastic scattering during the formation time is re-summed in the splitting kernel and hard elastic scattering and merging are parametrically rare so that neglecting them does not affect the theory's accuracy.When the parton later splits at τ split two new partons are created with kinematics following the massless case of final-state emitter and final-state spectator of Sherpa's Catani-Seymour shower described in [58], and they are assigned the remaining formation time of the original particle.The splitting probabilities h a (x) does not factor in any relative transverse momentum k ⊥ for the new parton pair, since they are derived from the AMY collision kernel which assumes a nearly exactly collinear process.Deviating slightly from this assumption, we allow for a small k ⊥ for the new parton pair sampled from the soft gluon spectrum dk 2 ⊥ /(k 2 ⊥ + k 2 ⊥,reg ) for some fixed k 2 ⊥,reg ≪ 1.The recoil is absorbed by an additional parton such that all partons remain on-shell.This recoil parton is chosen as the (spatially) closest located parton to the parton splitting (which in addition allows the outgoing partons to all have a momentum larger than our cutoff p min ).

Merging
In the case of evaluating if two particles merge our initial state is a pair, and so we cannot utilize the same methods as we have for particles splitting.Instead, we derive an effective cross section from the merging rate and follow the same approach that we use for elastic scatterings, i.e. to evaluate if two particles merge we consider their Lorentz invariant distance d 2 ij at their closest approach, and if it falls within the effective cross section of merging.In Appendix C.4 it is shown that the quasi-collinear splitting rate, γ a bc , is related to the effective matrix element |M a bc | 2 through It follows, by definition, that for two particles of momentum p and k merging into a particle p ′ the differential cross section can be expressed as with ) being the phase space measure for the outgoing particle.The resulting cross section is then where k 2 ⊥ is the relative transverse momentum between p and k.The Dirac delta present in the cross section enforces the exact collinearity of the merging processes, which would be the only kinematically allowed merging in vacuum.As with the splitting however, we circumvent this condition by introducing a recoil parton to absorb the transverse momentum.Hence, our merging process will allow partons with a non-zero k ⊥ to merge.We relate the Dirac-delta to the k 2 ⊥ -distribution used for sampling splittings as δ(k is the k 2 ⊥ -distribution normalised to unity and we pick up a factor of 1/2 since δ(k 2 ⊥ ) is not defined for k 2 ⊥ < 0 and so only integrates to 1/2.This gives us our final expression of the merging cross section, As in the case of elastic scattering the merging rates γ depend on the local quantity T * and the kernels also contain Bose enhancement/Pauli blocking factors, which we deal with in the same way as before: when integrating the effective cross section we overestimate the merging rate and the Bose enhancement/Pauli blocking factors and later reject mergings with the ratio of the true to and the overestimated effective cross section.
To calculate the cross section at the closest approach for a parton pair the energy fraction x is also needed, for which we use the Lorentzinvariant definition where P µ and K µ are the 4-momenta of the incoming partons, and R µ is the 4-momentum of the recoiler, as in [58].
As mentioned earlier, since the splitting rate is infra-red divergent we have to introduce a small cut-off in x.For consistency we apply the same cut-off to the merging rate.This does not compromise the theory's accuracy. 9nce a closest approach which falls within the cross section of the merging is found, a formation time is assigned to the particle pair in the same manner as for a splitting parton, see Section 4.4.1.If the merging is accepted in the rejection step to correct to the true cross section, the momenta are updated following (inversely) the massless case of final-state emitter and final-state spectator in [58], where the final-state spectator, or reocil parton, is chosen to be the parton with the closest spatial distance to the merging pair (that also allow for kinematics where no parton ends up with p < p min ).The remainder of the formation time of the original pair is also assigned to the new parton.

The effective temperature T *
In the quasi-collinear splitting/merging rate γ the quantity T * appears as a local variable (see Appendix B), this is the effective temperature at which an equilibrium system would have the same soft elastic scattering rate as the system under consideration.The definition of T * averaged over a spatial volume V is given by where we note that the denominator is equal to m 2 g , and the integrals again have to be calculated at fixed observed time and in the local rest frame.Our prescription to extract T * locally follows along the lines of the method used to find m 2 g/q , see Section 4.3.2,though with a key difference.In our treatment of m 2 g/q we considered all particles to be point-like, an assumption that for T * introduces problems with squared Dirac delta terms of the same variables.To remedy this we instead consider partons as Gaussian distributions with width σ p and σ x in the phase space density, i.e.
where the sum is over all partons i of species a.Assuming σ i,x = σ x and σ i,p = σ p for all partons i, we have for m 2 g T * that with σ x σ p ≥ 1/2 due to the uncertainty principle.With T * expressed as above we can proceed to extract it locally in the same manner as in Section 4.3.2.We include the N inc closest particles, where each particle i contributes with a term and each particle pair i, j of the same species s adds a term

The thermal equilibrium case
To verify the implementation of the splitting/merging kernel in Alpaca we once again utilize the simplicity of an infinite thermalized system.Since the splitting and merging rates are controlled by T * , which only equals T for f BE/F D , we will mainly focus our analysis to this initial distribution in the rest of this section, with the exception of looking at fixed γ.We will verify the following: (i) The total splitting and merging rates dN split /dt and dN merge /dt are equal, and correspond to independently extracted numerical values.(ii) The effective temperature T * is extracted so that it matches the analytical expected value.The effects of formation time are not of interest for these verifications, since we do not run inelastic and elastic processes simultaneously here, and so formation time for both splitting and merging have been set to zero.
(i) It follows from Eq. ( 31) that for an isotropic system the splitting rate can be related to the collision kernel as For the case of a purely gluonic system with constant γ g gg and no Bose enhancement it follows that the splitting and merging rates will not be the same for f BE/FD .Hence, we instead consider the same scenario for f = f C , where it follows that The thermal equilibrium splitting and merging rates have been examined in Alpaca for fixed γ in a gluonic system with no Bose-terms through a scan of the parameter space and the results are shown in Fig. 6.As can be seen, the rates deviate from the analytical value given above.This follows from the fact that the system relaxes to a different distribution when no Bose enhancement is included.The new distribution that the system relaxes into is extracted at the end of the run and used in Eq. ( 56) to find the corresponding expected splitting/merging rate, which our results fall within errors of.Fig. 6 Scaling behaviour of the splitting and merging rates in a gluonic system with fixed γ and no Bose enhancement factors are included.The system relaxes to a different phase space distribution, which is extracted at the end of the run.The resulting splitting/merging rate using the extracted distribution in Eq. ( 56) is shown as a purple line, with error bands in blue.The run parameters vary T between 0.3 − 0.4 GeV and γ between 0.001 − 0.005 GeV 2 .Fixed run parameters are p min = 0.01 GeV, ∆t = 250 GeV −1 and box volume Vx = 10.1 3 GeV −3 .
We have also looked at the case of a dynamic γ, extracted as described in Appendix B for a purely gluonic system with f = f BE .The result of this is found in Fig. 7 and shown to correspond well with the (independently) numerically extracted rate.(ii) The thermal equilibrium provides us with a simple analytical solution to T * for f BE/FD , which is The convergence of the extracted T * towards the analytical expectation w.r.t. the number of included particles is shown in Fig. 8.Note that it converges slower than m 2 s .This is due to the different implementation for extracting the numerator of T * compared to the denominator m 2 g .The probability distribution of the extracted T * , including a fixed number of particles, is shown in Fig. 9.The width of the distribution is due to the discrete sampling of f BE/FD , and will not decrease by increasing number of events.Lastly, a scan over different temperatures has been done to ensure that the effective temperature exhibits the correct scaling behaviour, as it is shown to do in Fig. 10.All values for T * shown in this section are extracted at initial τ in systems with box volume V x = 7.6 3 GeV −3 .

Putting everything together
Having verified all relevant parts of our implementation separately, we now proceed to the final verification, which is to combine all the pieces and let it run over longer timescales to ensure the dynamics of the system are correct, allowing it to remain in a thermal equilibrium.There are two caveats to this long run though.
we have introduced a global cutoff in momentum, p min , which no particle is initialized below, and each process which produces outgoing momenta k < p min is resampled until k ≥ p min .This means that the we cannot expect our system to relax exactly into f BE/FD .Hence, we initialize our system using a shifted initial distribution, This causes the analytical values of m 2 g/q and T * to shift compared to the regular Bose-Einstein and Fermi-Dirac case as well, and we extract them separately numerically to compare to our results in Alpaca.
Secondly, as mentioned in Section 4.3.2,during the evolution the state of the system is known at fixed τ , which means that all particles have different times.But the effective masses m 2 g/q and temperature T * have to be calculated at fixed time.The evolution of particles with later times than the required time can be backtraced.Since the future of particles with times earlier than the required time is still unknown, we let these particles free stream to the required time.This does, however, introduce a bias since splittings tend to take place at earlier τ than merging 10 .
The problem outlined above can be solved in an iterative fashion, by first running an event with m 2 g/q , T * and f (p) extracted with the currently available particles at any given τ .When the times t i for all particles i have passed some t max , the run ends, and the history of position, momenta and flavour for each particle over all t is saved.The system is then restarted and initialized with the same initial distribution that was sampled during the first event.During the course of the second event, m 2 g/q , T * and f (p) at any given t is then extracted at fixed time t using the history of the particles from the previous run.The increase/decrease of the dynamic quantities in the first run will affect the dynamics of the evolution, and so the values extracted in the second run will not be exactly correct, though it will iteratively converge toward the correct values by repeating the same procedure, always extracting the dynamic values from the previous run.However, though this procedure should converge in an equilibrium setting there is no guarantee that it will for a system initialized out of equilibrium.A possible solution to remedy this for non-equilibrium systems could instead be to utilize the Lorentz invariant parton cascade setup presented in [45], which is a modified version of what is described in Section 2. The study of non-equilibrium distributions is however left for future publications.
In this chapter we present the results side by side from a run without iteration and a run with one iteration step.We define these cases as N iter = 0 and N iter = 1 respectively.
For the full run presented in this chapter the following is set/included.
• Initialized with both gluons and quarks using f BE/FD,shift , T = 0.4 GeV, p min = 0.1 GeV, α s = 0.04, box volume V x = 7.6 3 GeV −3 , t min = 0 GeV −1 and t max = 1000 GeV −1 , 10 Since high energy particles are more likely to split while low energy particles are more likely to merge, we will at any fixed tau on average have had more splittings than mergings following Eq.( 9).For any fixed t, the number of particles is the same on average, but at fixed τ there are on average more particles in the system.This is not fully corrected by letting earlier particles free-stream and hence, when calculating screening masses and effective temperature in this way we have slightly more particles in the system than we should have.This increases the value of m 2 g/q which in turn decreases T * .
for 100 separate runs with different initial configurations.• All channels in 2 → 2 elastic scattering are allowed and evaluated including Bose/Pauli factors.• All channels in 1 → 2 inelastic splitting/merging are allowed and evaluated including Bose/Pauli factors.• The quantities m 2 g/q , T * and f (p) are extracted dynamically from the parton ensemble in each run.In the run where N iter = 1, the quantities are extracted using the parton position and momentum history from a previous run.
• Formation times are included for splitting and merging.The resulting distributions that the system relaxes into for gluons and quarks is shown in Fig. 11.As can be seen there, for gluons there is no significant difference in distribution at the end compared to the start, apart from a small drop in amplitude which can be attributed to the number of particles oscillating around the initial value, see Fig. 12.For quarks the distribution is shifted slightly to larger energies, as can also be observed Fig. 13, but is still within error bars of the original distribution.
In Table 5 see the mean splitting/merging rate over time.The total rates for inelastic splitting and merging are consistent with each other, for both the non-iterative and iterative runs, with a slight decrease in the rate for the latter (which can be attributed to a shift in T * , see below).We also see that detailed balance is preserved as the average splitting and merging rates for each channel overlap within errors.The total scattering rates for the two different runs can be found in Table 6, and we see an increase for N iter = 1 as expected, since too large m 2 g/q (see below) will suppress the scattering rate.
The rates are also visualized in Fig. 14 as number of scatterings, splittings and mergings as functions of time.No significant variation over time is observed, with the exception of a moderate increase of the scattering rate at early times, which could be due to the system settling into its actual equilibrium configuration (something similar can be seen in T * , cf.Fig. 16).The corresponding mean free time per particle for elastic scattering is t free,scatter = 38(3) GeV −1 and t free,split/merge = 53(4) GeV −1 for the noniterative case, and t free,scatter = 28(2) GeV −1 and t free,split/merge = 57(4) GeV −1 for the case of one iteration.
Fig. 12 shows that the total number of particles stays reasonably constant around the initialized values.In Fig. 13 the average energy per particle is shown to increase slightly for both the iterative and non-iterative runs, though seemingly remaining more constant in the latter case once the new equilibrium has been found.
Lastly, as discussed in the beginning of this section, we observe a drastic change in m 2 g/q and T * within the first ∆t = 100 GeV −1 in Figs. 15  and 16.As the particles spread out more in t the system settles into extracting too large values for m 2 g/q and consequently too low values T * for the non-iterative case.For the case of just one iteration, we see how the correct value of both m 2 g/q and T * is extracted, within oscillations.Since the effective masses and temperature are consistent with the true values after just one iteration, it can be expected that further iterations will not lead to significant changes in the evolution.When the true values are not known, one has to find out how many iterations are needed by comparing the results from subsequent iterations.
To summarize, when running our system initialized in an infinite thermal equilibrium, with all relevant quantities needed to faithfully reproduce the AMY Kinetic Theory extracted dynamically during the run, the system remains in this equilibrium (up to small oscillations in particle number and minor shifts in average energy per particle). ] ]  N iter Type dN 1↔2,tot /dt GeV dN g↔gg /dt GeV dN q↔gq /dt GeV dN g↔q q /dt GeV 0 Split Merge Table 5 Comparison of the splitting and merging rates for the long run with N iter ∈ {0, 1}.
N iter dN 2→2,tot /dt GeV 0 3.9(3) 1 5.3(4) Table 6 Comparison of the total elastic scattering rate for the long run with N iter ∈ {0, 1}.Nq N q,analytical,pmin=0 =92.4 Ng: Standard deviation between events Nq: Standard deviation between events Fig. 12 Mean number of particles in the system as a function of t.The error bands correspond to the standard deviation between events.Left: N iter = 0. Right: N iter = 1.[GeV] N iter = 1, f BE/FD,shift : m 2 g/q vs t, T = 0.4 GeV, p min = 0.1 GeV Analytical: m 2 g = 0.047 GeV 2 Analytical: m 2 q = 0.021 GeV 2 ALPACA: mean(m 2 g ) ALPACA: mean(m 2 q ) Fig. 15 Mean extracted m 2 g/q as a function of time.The error bars correspond to the standard deviation between events.Left: N iter = 0. Right: N iter = 1.

Conclusions
The AMY effective kinetic theory of QCD at high temperatures has been successfully applied to various aspects of heavy ion collisions and is a natural candidate for the description of small collisions systems.We therefore believe it is time to turn it into a phenomenology tool by constructing a Monte Carlo event generator that solves the AMY Boltzmann equations by explicit simulation.Monte Carlo event generators have the advantage that in principle any observable can be calculated from an event sample (instead of having to do a new calculation for each observable) and that results can be directly compared to experimental data.
We have introduced the parton cascade Alpaca as a part of the multi-purpose event generator Sherpa.Alpaca encodes the AMY collision kernels in a Lorentz invariant framework for simulating the dynamics of multi-particle systems.The ensemble average is thus a solution of the Boltzmann equation.Frame independence is achieved by ordering collision and splitting events not in a frame dependent observer time but in a frame independent generalised time [44].
The partons can undergo elastic scattering and effective quasi-collinear splitting or merging, where in the latter coherent multiple soft scattering during the formation has to be taken into account.Both types of processes, elastic and inelastic, depend on the phase space densities of partons through the effective masses m g and m q related to screening effects in the medium and the effective temperature T * .These quantities are local in position space but require integration over the momentum.We show that an estimate of the effective masses and temperature can be obtained from a single event without the need to provide further information about the phase space densities.In thermal equilibrium the analytical values are recovered with increasing precision as more and more partons are included in the estimate.However, in order to get an accurate value a large number of up to 30 particles has to be included.It remains to be seen whether this is good enough for out-of-equilibrium systems.In a similar way estimates of the local phase space density needed for Bose enhancement and Pauli blocking factors can be obtained.
To further validate the Alpaca framework we show that in thermal equilibrium • the elastic scattering rate for fixed cross section exhibits the expected scaling with temperature, volume and cross section • the effective masses have the correct value and the expected dependence on temperature • the elastic scattering rates with dynamical matrix elements reproduce the expected results within a few percent • for fixed γ the splitting/merging rate exhibits the expected scaling with γ, temperature and volume • the effective temperature scales as expected with temperature • with dynamical γ the splitting and merging rates are the same and agree well with numerical results • when putting everything together and running over longer time scales the system remains in thermal equilibrium This gives us confidence that Alpaca indeed faithfully reproduces the AMY dynamics.Going from thermal equilibrium to out-of-equilibrium systems requires only small extensions.The next step will thus be to apply Alpaca to out-ofequilibrium systems to study equilibration and to the modeling of heavy ion collisions.

A Box with periodic boundary conditions in Lorentz invariant framework
The Lorentz invariant distance (squared) between to particles i and j is given by where the positions of the particles are functions of τ x µ (τ ) = 2λp µ (τ − τ 0 ) + x 0µ (61) with x 0µ = x µ (τ 0 ) = (t 0 ; x 0 ) and the arbitrary but fixed Lagrange multiplier λ.
The closest approach (i.e. the minimum of d 2 ij ) is reached for τ = τij , which is given by The distance of closest approach is then simply given by d 2 ij (τ ij ), which reduces to From now on we drop the index µ from fourvectors and introduce the notation [ab] for the scalar product of two four-vectors a and b.Threevectors are denoted by boldface symbols.
We consider here a cubic box 11 with side length L. A copy of a particle j in another cell is generated by shifting the initial position x 0 along the coordinate axes by multiples of the side length L x 0j −→ x 0j + (0; k, l, m)L = x 0j + s (64) with integer k, l, and m.
The "invariant time" of closest approach τij with shifted particle j is then given by where τij,0 is the invariant time of closest approach without shift (i.e.k = l = m = 0).The invariant squared distance becomes where x i/j (τ ) refers to the particles' positions without shift and d 2 ij,0 (τ ) denotes the squared invariant distance without shift.Similarly, the invariant distance squared at closest approach becomes Our task is to find the triple (k, l, m) that produces the first scattering, i.e. the first encounter with d 2 ij (τ ij ) < σ/π.To achieve this, it is helpful to first consider k, l, and m as continuous variables that span a three-dimensional space (R 3 ).In this space τij is constant on planes, whose position and orientation are determined by Eq. ( 65).Equating this to some τ 0 and solving for m yields the value of m for which, given k and l, the invariant time of closest approach equals τ 0 .
It is also worth noting that the direction of fastest increase of τij is given by p i − p j (the gradient of τij w.r.t.(k, l, m)).Furthermore, in the (k, l, m) space, the invariant distance of closest approach d ij (τ ij ) vanishes on a line.It can be found from the condition that the gradient (again w.r.t.(k, l, m)) vanishes.The partial derivatives of Eq. ( 67) are given by where Inserting this into Eq.( 68) one has all three coordinates of the point of vanishing d 2 ij (τ ij ) for a given τij = τ 0 .
The strategy for finding the next scattering of a pair of particles is now to start at this point, choosing τ 0 to be the current invariant time τ of the evolution.From this point one then moves along the line of vanishing d 2 ij (τ ij ) in the direction of increasing τij .For the points with integer k, l and m lying closest to the line it is checked whether they lead to a scattering, i.e. whether d ij (τ ij ) < σ/π.The first such point is the sought after next scattering 12 .In practice, the checking of nearby points is done by imagining the (k, l, m) space to be divided into cubes the corners of which are the points with integer k, l and m.One can then move from cube to cube by calculating through which sides the line enters and leaves the cube using Eqs.( 72) -(77).The corners of the cubes intersected by the line are the points that have to be checked.

B Splitting/merging rates
The splitting/merging rates γ a bc are products of the spin averaged Altarelli-Parisi splitting functions and a function µ 2 .For the possible splitting/merging processes one gets where η = x(1 − x)λT * p/m 2 g with the coupling constant given by λ = g 2 N c = 4πα s N c .The collinear splitting/joining rates are dictated by 2-dimensional quantum mechanics in the transverse plane of the splitting particle.The rate splitting rate is related to µ 2 (η, x; , s 1 , s 2 , s 3 ) = √ 2(4π)m 2 g Im(f (0)) where f is the solution to the radially symmetric two-dimensional Schrödinger equation where and the boundary conditions are given by f The effective mass appearing in the equation is the mass difference 13 The collision kernel C appearing in the Schrödinger equation is given by with where m D is the Debye mass.
The quantity µ can be easily evaluated in two limits, that of small η/| M 2 | corresponding to the Bethe-Heitler limit, and the large η limit corresponding to the deep LPM limit [56].The 13 For the interesting special case of d A = 8, C A = 3, N f = 3, C F = 4/3, i.e.QCD with three flavours of quarks, the ration of the effective masses of quarks and gluons in independent of the phase space density and given by m 2 g /m 2 q = 9/4.
expression in the Bethe-Heitler limit read 14 with where 17.Expanding the expression further for small x gives a particularly simple form that can easily be used in analytical estimates In the limit | M 2 | → 0, to the next to leading logarithmic order 15 , (105) 14 The Green's function of the free equation is G Using the Green's function the small η solution can be written as G1(b, 0). 15 The expression here differs from that of [56] by log(x) → log(x + 1).The finite shift of the argument of the log does not affect the Next to Leading Logarithmic order but makes the expression well behaved at small values of η.Fig. 17 Function Q BH (r) appearing in the solution of the splitting/merging rates in the Bethe Heitler limit.
For values of η and | M | for which neither of the limits gives a good approximation, the equation can be solved numerically usingWe find that the numerical result is well approximated by the following interpolating function, which is indeed what we use in the simulation: (106) Fig. 18 shows the numerical values of µ along with the two limits and the interpolating function for the representative cases of µ 2 (η, x = 0.5; A, A, A), µ 2 (η, x = 0.5; F, A, F ) and µ 2 (η, x = 0.5; A, F, F ), which corresponds to all different splitting processes at x = 0.5.The interpolating function gives a good approximation of the numerical solution.

C.1 Elastic scattering matrix elements
Purely in terms of Mandelstam variables the total cross section becomes (107) and we overestimate this as Below the different matrix elements used, their overestimates and integrated expressions are presented.The prefactors ζ g = e 5/3 /4 and ζ q = e 2 /4 are discussed in Section 4.3.1.It is assumed that ν q = ν q .Note that the equations for q q → gg are not symmetric, to get the corresponding expressions for gg → q q one has to exchange the degeneracy factors in front as ν 2 q → ν 2 g . .At small values of η the scale is well described by the Bethe-Heitler approximation (green line) whereas at large values the curve approaches (logarithmically) the deep-LPM approximation of [56] (purple line).The interpolation function constructed from the two limits describes the numerical data adequately at all values of η and x for N f = 3 QCD.Black lines are extracted from Alpaca, while other lines are analytical values, with an almost complete overlap between the two.Bottom row: same as top for µ 2 (η, x; F, A, F ) (left) and µ 2 (η, x; A, F, F ) (right).
× − 2ζ 4 g m 8 g + 3ζ , (114) , (115)  t(5ζ q m 2 q + 13s) s + ζ q m 2 q + 2t 2 + 4(s 2 + ζ 2 q m 4 q ) log(s + t + ζ q m 2 q ) s + ζ q m 2 q , (120) In order to verify the elastic scattering rates of Alpaca in the thermal case with a non-constant cross section, we have numerically integrated the collision kernels given in AMY and compared the results to the scattering rates of Alpaca through the relation found in Eq. (32).We have looked separately at C[f ] + and C[f ] − , which corresponds to the gain and loss terms in the collision kernel.The notation below for the coefficients is as follows, 2 i is due to the particles in the initial state being indistinguishable, 2 d is due to double-counting of the incoming particles, 2 b is due to quark and antiquark both contributing and 3 f (or 2 f ) is due to the three flavours currently in Alpaca.Hence, for particles scattering out of p we have for gluons that the loss is and the gain is × δ (4) The factor of 2 b in front of gq ↔ gq is due to that full row being equivalent to g q ↔ g q, assuming f q = f q = f g .For a quark of species s we have that the loss is = I −,qiqj ↔qiqj q + I −,qiqi↔qiqi q + I −,qi qi↔qi qi q + I −,qi qi↔qj qj q + I −,qg↔qg q + I −,q q↔gg q (125) and the gain is where the factor 2 b on the first matrix element is due to the fact that q i q j ↔ q i q j is equivalent to q i qj ↔ q i qj .Lastly, C −,2↔2 q [f, p] = C −,2↔2 q [f, p] and C +,2↔2 q [f, p] = C +,2↔2 q [f, p].The integrals above have been solved numerically, separately from Alpaca, using Monte Carlo methods with stratified sampling, and compared to the elastic scattering rate in Alpaca through the relation shown in Eq. (32).The results of the full integration are shown in Section 4.3.3.

C.3 Splitting probability overestimate
To find an overestimate ha (x) ≥ h a (x), with h a (x) as defined in Eq. ( 42), we note the inequality relation (128) and so Note also that for all splitting channels it holds that P a bc (x) ≤ P a bc (x min ), with the condition that x max = 1 − x min .Combining all of the inequalities presented above we end up with the three following overestimates, h g→gg (x) ≤ hg→gg = (2π) h g→q q (x) ≤ hg→qq = (2π) 3 2|p|ν q P g q q (x min )μ (137)

fFig. 5
Fig. 5 Extracted phase space density, f (x, p), for thermal systems with temperature T = 0.4 GeV.The values are extracted at initial t.Left: f C .Right: f BE/FD .

Fig. 7
Fig.7The splitting and mering rate for a gluonic thermal system with γ calculated as described in Appendix B. The figure shows results from 125 runs with different initial configurations sampled from a thermal distribution with τmax = 2500 GeV −1 , αs = 0.04 and box volume Vx = 10.1 3 GeV −3 .

fFig. 8 fFig. 9
Fig.8Scaling behaviour of T * and m 2 g T * (see Eq. (51)) as a function of number of particles included, N inc .

fFig. 10
Fig. 10 Scaling behaviour of T * as a function of different initialization temperatures T for f BE/FD .

Fig. 11
Fig. 11 Distribution of the system extracted after the run is finished, averaged over the last ∆t = 200 GeV −1 .The blue error bands correspond to mean of the fluctuations within each event over time, while the green error bands correspond to the standard deviation between events.Left: N iter = 0. Right: N iter = 1.Top: fg(p).Bottom: fq(p).

Fig. 13 NNFig. 14 2 ]N
Fig. 13 Mean between events of average energy per particle as a function of t.Left: N iter = 0. Right: N iter = 1.

NNFig. 16
Fig. 16 Mean of T * as a of time.The error bars correspond to the standard deviation between events.Left: N iter = 0. Right: N iter = 1.
) Setting Eqs.(78) and (79) equal to zero and solving for k and l yields k

Fig. 18
Fig.18Top: Scale µ 2 (η, x; A, A, A) entering the splitting rate of gluons for QCD.The dots represent numerical evaluation of µ 2 .At small values of η the scale is well described by the Bethe-Heitler approximation (green line) whereas at large values the curve approaches (logarithmically) the deep-LPM approximation of[56] (purple line).The interpolation function constructed from the two limits describes the numerical data adequately at all values of η and x for N f = 3 QCD.Black lines are extracted from Alpaca, while other lines are analytical values, with an almost complete overlap between the two.Bottom row: same as top for µ 2 (η, x; F, A, F ) (left) and µ 2 (η, x; A, F, F ) (right).