Supersymmetry and Generic BSM Models in PYTHIA 8

We describe the implementation of supersymmetric models in PYTHIA 8, including production and decay of superparticles and allowing for violation of flavour, CP, and R-parity. We also present a framework for importing generic new-physics matrix elements into PYTHIA 8, in a way suitable for use with automated tools. We emphasize that this possibility should not be viewed as the only way to implement new-physics models in PYTHIA 8, but merely as an additional possibility on top of the already existing ones. Finally we address parton showers in exotic colour topologies, in particular ones involving colour epsilon tensors and colour sextets.


Introduction
The Standard Model (SM) of particle physics has been enormously successful in describing interactions between fundamental particles. The only experimentally unverified component of the SM is the Higgs boson which is thought to underlie the electroweak symmetry breaking (EWSB). However, experimental evidence for dark matter, neutrino masses and the theoretical requirement of naturalness appeals for a theory beyond standard model (BSM). Monte Carlo generators fulfil an important role, both in testing the SM to high precision and in testing the prediction of new theories, by providing a systematic procedure of comparing theoretical prediction to experimental observation.
Pythia 8 [1] is a general-purpose Monte Carlo event generator [2] for a full simulation of high-energy collision events. It includes a comprehensive library of hard-scattering processes, particle decays, initial-and final-state parton-shower models [3,4], hadronization through string fragmentation [5] and models of beam remnants and multiple interactions [6,7]. It contains a native implementation of a wide variety of SM and BSM processes and also provides a standard interface [8,9] to external programs which may be used by a standalone generator.
We describe here, the updates to the Pythia 8 event generator to include the popular BSM model of supersymmetry (SUSY), additions made to parton showers and hadronization algorithms to allow for exotic colour topologies and generic enhancements made to enable interfacing to parton-level BSM generators. In section 2, we describe the implementation of supersymmetric models in Pythia 8, including production and decay of superparticles. In section 3, we present a framework for importing generic new-physics matrix elements into Pythia 8, in a way suitable for use with automated tools. In section 4, we discuss the treatment of parton showers in exotic colour topologies. Section 5 contains a brief summary and conclusions. 1

Supersymmetry in Pythia 8
Supersymmetry (see [10] for a pedagogical introduction) is considered one of the best motivated extensions of the SM due to its ability to address many outstanding theoretical and experimental issues. In particular, the Minimal Supersymmetric extension of the Standard Model (MSSM) is currently a popular candidate for a BSM theory. The MSSM extends the SM by the addition of one pair of SUSY generators which implies the presence of one superpartner to each SM state. The MSSM particle spectrum therefore has squarks (q i ), sleptons (l i ) and gauginos (B,W i andg) as the supersymmetric counterparts of quarks, leptons and gauge bosons respectively. The requirement of self-consistency of the theory via anomaly cancellation also demands two Higgs doublet fields H u and H d . After electroweak symmetry breaking (EWSB), we are left with five Higgs degrees of freedom viz. the CP-even h 0 and H 0 , the CP-odd A 0 and two charged Higgs bosons H ± . The superpartners of the Higgsesthe fermionic "Higgsinos" -mix with the gauginos to form neutralinos and charginos. In particular, the neutral Higgsinos (H 1 andH 2 ) mix with the neutral U (1) and SU (2) gauginos (B andW 3 ) to form the mass eigenstates called the neutralinos (χ 0 i ; i = 1 − 4.) Similarly, the charged Higgsino mixes with the charged SU (2) gaugino to form charginos (χ ± i ; i = 1, 2.) The next-to-minimal supersymmetric extension of the SM (nMSSM) extends this scenario by adding one extra singlet Higgs field. This adds another member to the neutralinos and the neutralino mixing matrix is enlarged to 5 × 5. The current implementation of Pythia 8 includes the nMSSM extension and allows processes with CP, flavour or R-parity violation.
Pythia 8 uses the standard PDG codes for numbering the superpartners [11] and the particle spectrum is read in via an SLHA file [12,13]. We use the super-CKM basis (in the conventions of the SLHA2 [13]) for describing the squark sector which allows non-minimal flavour violation. The mass-eigenstates of the squarks are then related to the left-and righthanded squarks via a 6 × 6 complex mixing matrix. Our implementation can therefore be used to study both CP violation and flavour violation in the squark sector.

Couplings
Pythia 8 reads particle masses and mixing matrices via the SUSY Les Houches Accord (SLHA2) framework [13]. (Read-in of SLHA1 spectra [12] is also supported, but mixing the two standards is strongly discouraged, as the internal translation from SLHA1 to SLHA2 has only been designed with the original SLHA1 in mind.) The raw data read in by the SusyLesHouches class is accessed by the CoupSUSY class which uses the information to construct all the SUSY couplings. The couplings are defined according to [14] for all cases except for couplings of superparticles to Higgs bosons which are defined according to [15].
The running of electroweak and strong couplings is carried over from the corresponding one-loop calculations in the Standard Model. The GAUGE block can be used to set the boundary values of all three SM couplings at the SUSY breaking scale. By default, the masses of W and Z are assumed to be the pole masses and are used to calculate the on-shell value of If externally provided in the SLHA file, the value of sin θ W can be set to the running value using the flag SUSY:sin2thetaWMode = 2 (see the Pythia 8 HTML user reference included with the code [1]). The ratio of the two Higgs vacuum expectation values (tan β) is read in from the low scale value provided by the MINPAR and EXTPAR blocks. The default value of the Higgs mixing angle (α H ) is set to the SM limit (β − π/2) which is then overwritten by the contents of the HMIX block.
Since the SLHA interface has been extended and can now be used to pass information on any new particles and decays [9], the presence of the MODSEL block is used as an indicator of SUSY models and Pythia 8 will initialize the CoupSUSY class only if this block is present. Skipping the MODSEL block is acceptable for Les Houches Event files (LHEF) as long as the user supplies an external decay table for all required cascade decays.

R-parity violation
The most general MSSM superpotential allows both lepton and baryon-number violating processes. This is generally avoided by demanding invariance under an R-parity defined as (−1) 3B−L+2S . From this definition, all SM particles are even whereas all superpartners are odd under R-parity. A well known consequence of this is that the Lightest SUSY particle (LSP) must be stable. A neutral, weakly interacting LSP can therefore be a good candidate for dark matter. However, the imposition of R-parity can be considered an aesthetic requirement rather than a consistency requirement and possible R-parity violating interactions, if present, can be probed by collider experiments. We therefore include R-parity violating production and decay processes in our implementation.
In SLHA conventions, the R-parity violating superpotential is given by The µ-type terms correspond to bi-linear R-parity violation which causes a mixing between the leptons and neutralinos/charginos. The λ and λ ′ -type terms lead to lepton number violation whereas λ ′′ -type terms lead to baryon-number violation. The current implementation does not include the effects of the bi-linear term. The R-parity violating couplings λ ijk are antisymmetric under i ↔ j. Therefore only couplings for i > j are read and the rest are set by the symmetry property. Similarly, λ ′′ ijk is antisymmetric under j ↔ k and hence only couplings with j > k need to be provided. This implementation includes in particular, the resonant production of a squark via λ ′′ -type couplings which can be probed at hadron collider experiments. The changes made to showering and hadronization to account for the non-standard colour structure from such terms will be explicitly described in section 4.
resonant production of squarks via baryon number violating couplings. All available SUSY processes can be turned on using SUSY:all = on. Individual subprocesses can then be selected based on the final state by setting SUSY:idA = PDGcode and SUSY:idB = PDGcode. If only idA is provided, all processes with that particle in the final state are turned on. Alternatively, one or more production processes can be turned on using the string SUSY:processname = on, again with SUSY:idA and SUSY:idB providing a further level of subprocess selection. The available subprocess classes are listed in Table 1.
The squark-antisquark and squark-squark production processes include contributions from EW diagrams and their interferences. To estimate the size of these contributions, and/or for purposes of comparison to other codes that do not include them, the cross sections can be restricted to include only the strong-interaction contributions, using the following flags: • qqbar2squarkantisquark:onlyQCD = true.
The baryon number violating coupling λ ′′ ijk if present, can induce resonant squark production via the process d j d k →ũ * i which produces a resonant up-type antisquark or via u i d j →d * k or u i d k →d * j which produce a down-type antisquark. The expression for an up-type squark production process is The expression for down-type squarks is similar, taking into account the symmetry property λ ′′ ijk = −λ ′′ ikj . We implement this production process as qq2antisquark and the charge conjugate process (q iqj →q k ) is included by default.
The supersymmetric Higgs sector is identical in many ways to the Two-Higgs Doublet Model. The Higgs production processes have already been implemented in Pythia 8 in the SigmaHiggs class. The production of the Higgs bosons can be accessed by including the switch HiggsBSM:all=on. For specific Higgs processes, please refer to the HTML user reference included with the code [1].

Sparticle Decays
SUSY Particle decays are handled by the class SUSYResonanceWidths. The user can choose to read in decay tables via SLHA or use the decay widths calculated by Pythia. As a default, Pythia does not calculate the decay width if a table is externally supplied. Note, however, that while Pythia's internal treatment can include sophistications such as matrix-elementbased phase-space weighting and running widths, channels read in from an SLHA decay table will be decayed purely according to phase space, with no matrix-element weighting. The internal treatment should therefore be preferable, in most cases, and an option for overriding the automatic read-in of decay tables is provided, by setting the flag SLHA:useDecayTable = false, see sec. 3.1.
The decay of a particular particle may be turned off manually using the standard Pythia 8 structure PDGcode:mayDecay = false or by setting its width to zero in the SLHA decay table. In the former case, the particle will still be distributed according to a Breit-Wigner distribution with non-zero width, whereas it will always be assigned its pole mass in the latter.
Individual decay modes may be switched on/off using the standard Pythia 8 methods, documented in the section on "The Particle Data Scheme" in the program's HTML documentation [1]. We discuss ways to switch modes on/off using SLHA decays tables in section 3.1.
The internal treatment of 2-body decays is so far restricted to on-shell particles. A mechanism for effectively generating 3-body decays via sequences of 1 → 2 decays involving off-shell particles is foreseen as an update in the near future (and will be announced in the Pythia 8 update notes). An equivalent mechanism is already implemented in Pythia 8, e.g., for h → ZZ decays for light Higgs bosons.
Currently the following R-parity conserving two-body decays are implemented: Besides these, we also include two-body R-parity violating decays of squarks via λ ′ (q → lq ′ ) and λ ′′ -type couplings (q → q ′ q ′′ ). We also include the three-body decays of neutralinos through λ ′′ -type couplings via an intermediate squark [16]. For certain final states in three body decays, partial decay via sequential two-body decays may also be kinematically allowed. In this case, we demand that only the off-shell components of the matrix element-squared are allowed to contribute to the three-body decay width. Any interferences between the off-shell and on-shell components are also turned off. The two-body sequential decays then proceed as normal.
The Higgs boson running widths are calculated in the associated classes ResonanceH for CP even (h 0 , H 0 ) and the CP odd (A 0 ) Higgses, and ResonanceHchg for charged Higgses(H ± ). By default, the Higgs decay tables are not overwritten even if they are read via SLHA because Pythia 8 performs a more accurate phase space calculation than the flat weighting that is performed for decay widths read in via SLHA. The decays of Higgses into SUSY particles will be included in a future update.

Interfacing Generic BSM Models
The simplest way of implementing a new model may often be to just inherit from SM or BSM processes that have already been implemented in Pythia 8, modifying and generalizing them as appropriate, as described in the section on "Semi-internal Processes" in the main Pythia 8 documentation [1].
Alternatively, Pythia 8 can read in parton-level events generated by external matrixelement event generators [17][18][19][20][21], using the Les Houches Event File (LHEF) format [8,22]. If the events contain new particles, so-called QNUMBERS blocks [9,23,24], described in sec. 3.1 below, can be used to add information on the quantum numbers of new particles, and SLHA decay tables [12] may also be provided. (For SUSY models, in addition, complete SLHA spectra can be given, as discussed in sec. 2.) The encoding of colour flow is then particularly important, for the events to be showered and hadronized correctly. Some pedagogical examples, with illustrations, are given in the original LHA paper [22], and further explicit examples with colour-epsilon and colour-sextet structures are given in section 4 below. The LHEF paper [8] describes how to encode this in an LHE file, with examples of correct LHE files available, e.g., in Pythia's examples/ directory. When reading events from LHE files, the BSM/SLHA information may either be enclosed within the LHE file (preferred), or provided in a separate file. In the former case, the BSM/SLHA information should be included in the <header> part of the LHE file [8], inside an <slha> tag [9]. In the latter case, a separate BSM/SLHA file may be specified using the Pythia 8 command SLHA:file = fileName. The mode SLHA:readFrom gives the user some additional control over whether and from where BSM/SLHA information is read in. It should normally be left at its default setting, but can optionally be used either to switch off SLHA read-in entirely, or to force read-in from a specific file: SLHA:readFrom = 0 # do not read BSM/SLHA information at all = 1 # (default) read in from the <slha>...</slha> block of a LHEF, # if such a file is read during initialization, and else from # SLHA:file = 2 # read in from SLHA:file The framework described in section 3.2 represents a third option which combines features from both of the two possibilities above. It allows parameters and matrix-element code to be imported directly from external packages, to generate semi-internal processes in Pythia 8 (i.e., without an intermediate LHE file) in a fully automated and generic way, as long as the final-state parton multiplicity does not exceed the limitations of Pythia's internal hard-process phase-space generator [1]. A working interface between Pythia 8 and MadGraph 5 [20] has been constructed along these lines, for 2 → 2 processes, and will be reported on separately. Here, we focus on the Pythia 8 side of the interface.
The interface basically consists of two parts: 1) information about particles and couplings using a generalized SLHA format (section 3.1), and 2) accessing that information from within a semi-internal Pythia 8 process (section 3.2).

Information about New Particles
Information about particle quantum numbers, masses, couplings, and decays, can be given in an ASCII file, using a generalization of the SLHA [12] and BSM-LHEF [9] formats, whose name is provided to Pythia by setting the word SLHA:file = fileName.

QNUMBERS
The SLHA file should contain a QNUMBERS block [9] for each state not already associated with an ID code (a.k.a. PDG code, see [1,11] for a list) in Pythia 8. For a hypothetical electrically neutral colour-octet self-conjugate fermion (a.k.a. a gluino) that we wish to assign the code 7654321 and the name "balleron", the structure of this block should be Note that the name(s) given after the # mark in the block definition are optional and entirely up to the user. If present, they will be used, e.g., when printing out event records with Pythia's event.list() method.
The SM quantum numbers given in the QNUMBERS blocks are required by Pythia 8 for QED and QCD showering, and for colour-flow tracing. (Currently, Pythia does not make use of the spin information.) As a rule, we advise to avoid clashes with existing ID codes, to the extent possible in the implementation. A useful rule of thumb is to only assign codes higher than 3 million to new states, though one should be careful not to choose numbers larger than a 32-bit computer integer can contain, which puts a cap at ∼ 2 billion.

MASS
The file should also contain the SLHA block MASS, which must, as a minimum, contain one entry for each new state, in the form In principle, the block can also contain entries for SM particles. Here, some caution and common sense must be applied, however. Allowing SLHA spectra to change hadron and/or light-quark masses in Pythia 8 is strongly discouraged, as these parameters are used by the parton-shower and hadronization models. Changing the b-quark mass, for instance, should ideally be accompanied by a retuning of the b fragmentation parameters. Since this is not the sort of question a BSM phenomenology study would normally address, by default, therefore, Pythia 8 tries to protect against unintentional overwriting of the SM sector via the flag SLHA:keepSM, which is on by default. To be more specific, this flag causes particle data (including decay tables, see below) for ID codes in the ranges 1 -24 and 81 -999,999 to be ignored. Notably this includes Z 0 (ID 23), W ± (ID 24), and t (ID 6). The SM Higgs (ID 25), however, may still be modified by the SLHA input, as may other particles with ID codes in the range 25 -80 and beyond 1,000,000. If you switch off this flag then also SM particles are modified by SLHA input.
Alternatively, the parameter SLHA:minMassSM, with default value 100.0 GeV, can be specified to allow any particle with ID code below 1,000,000 to be modified, if its default mass in Pythialies below some threshold value, given by this parameter. The default value of 100.0 allows SLHA input to modify the top quark, but not, e.g., the Z 0 and W ± bosons.

DECAY
The file may also include one or more SLHA decay tables [12]. New BSM particles without decay tables will be treated as stable by Pythia 8. For coloured states, this may result in errors at the hadronization stage, and/or in the possibly unintentional production of so-called R-hadrons [25], with a reasonably generic model for the latter available in Pythia 8 [26]. On the other hand, a redefinition of Pythia's treatment of the decays of SM particles, like Z 0 and W ± may be undesirable, since Pythia's internal treatment is normally more sophisticated (discussed briefly in sec. 2.4). Thus, again, caution and common sense is advised when processing (B)SM particles through Pythia, with the protection parameters SLHA:keepSM and SLHA:minMassSM described above also active for decay tables. An option for overriding the automatic read-in of decay tables is also provided, by setting the flag SLHA:useDecayTable = false.
The format for decay tables is [12] # ID WIDTH in GeV DECAY 7654321 2.034369169E+00 # balleron decays # BR NDA ID1 ID2 ID3 9.900000000E-01 Note that the branching ratios (BRs) must sum up to unity, hence zeroing individual BRs is not a good way of switching modes off. Instead, Pythia 8 is equipped to interpret a negative BR as a mode which is desired switched off for the present run, but which should be treated as having the corresponding positive BR for purposes of normalization.
Finally, a note of warning on double counting. This may occur if a particle can decay via an intermediate on-shell resonance. An example is H 0 → q 1q2 q 3q4 which may proceed via H 0 → W W followed by W → qq ′ . If branching ratios for both H 0 → W W and H 0 → q 1q2 q 3q4 are included, each with their full partial width, a double counting of the on-shell H 0 → W W contribution would result. (This would also show up as branching ratios summing to a value greater than unity.) Such cases should be dealt with consistently, e.g., by subtracting off the on-shell contributions from the H 0 → q 1q2 q 3q4 partial width.

Accessing the Information from a Semi-Internal Process
Already the original SLHA1 [12] allowed for the possibility to create user-defined blocks, beyond those defined by the accord itself. The only requirement is obviously that the block names already defined in the accord(s) should not be usurped. The SLHA interface in Pythia 8 will store the contents of all blocks, both standard and user-defined ones, as internal vectors of strings.
By default, Pythia's internal BSM implementation only extracts numerical content from those blocks it recognizes (i.e., the standard SLHA 1&2 blocks and QNUMBERS), and uses those to initialize its couplings and particle data arrays. However, generic methods are also provided, that can be used access to the contents of any block, whether standard or user-defined, from inside any class inheriting from Pythia's SigmaProcess class (i.e., in particular, from any semi-internal process written by a user), through its SLHA pointer, slhaPtr, by using the following methods: bool slhaPtr->getEntry(string blockName, double& val); bool slhaPtr->getEntry(string blockName, int indx, double& val); bool slhaPtr->getEntry(string blockName, int indx, int jndx, double& val); bool slhaPtr->getEntry(string blockName, int indx, int jndx, int kndx, double& val); This particular example assumes that the user wants to read the entries (without index, indexed, matrix-indexed, or 3-tensor-indexed, respectively) in the user-defined block blockName, and that the entry value, val, should be interpreted as a double. In fact, the last argument is templated, and hence if anything other than a double is desired to be read, the user has only to give the last argument a different type. Since the user presumably knows what type of content his/her own user-defined blocks contain, this solution allows the content to be accessed in the correct format, without Pythia needing to know what that format is beforehand. If anything goes wrong (i.e., the block does not exist, or it does not have an entry with that index, or that entry cannot be read as a double), the method returns false; true otherwise. This effectively allows input of completely arbitrary parameters using the SLHA machinery, with the user having full control over names and conventions. Of course, it is then also the user's responsibility to ensure complete consistency between the names and conventions used in the SLHA input, and those assumed in any user-written semi-internal process code.
Note also that the special SLHA block SMINPUTS (containing SM parameters [12]) will always be accessible through the methods above, regardless of whether a corresponding SLHA block has been read in or not. The SMINPUTS block is initialized starting from PYTHIA's own internal default values, with subsequent modifications as dictated by updates to PYTHIA's particle and parameter databases before initialization and/or by SLHA read-in. This functionality is intended to give a generic BSM implementation access to the SM parameters contained in SMINPUTS in a universal way.
To give a specific example, the interface to Madgraph 5 was structured in the following way. Among the possible output formats available for matrix elements in Madgraph 5, one is a mode called pythia8. When invoked, this mode writes out the corresponding matrix element(s) in exactly the format required by Pythia 8's semi-internal process machinery. The resulting code can therefore be imported directly into Pythia 8, and Madgraph even provides explicit instructions and a Makefile for doing precisely that. In general, however, such matrix elements may contain parameters that refer, e.g., to couplings in a model unknown to Pythia. A central question was therefore how to provide information on such parameters at runtime, in a sufficiently generic manner. The solution is that Madgraph writes out the relevant parameters as custom-made SLHA-like blocks in a BSM/SLHA file included together with the matrix-element code. It then also inserts appropriate calls to slhaPtr->getEntry() in the cross-section expressions, so that each parameter can be retrieved when needed, without any user intervention required.
Note that this entirely circumvents a particularly troublesome issue that before was thought to make any truly universal "BSM Accord" impractical, the problem of agreeing on a common standard for names and parameters for completely arbitrary models. In the Madgraph-Pythia interface, it is sufficient that Madgraph itself assigns some unique names and contents to each block. It has complete freedom in choosing which conventions to use, as long as it correctly inserts the corresponding readEntry() calls in its matrix-element output. This effectively generalizes the SLHA data structure to apply to completely general BSM models.
The interface has been tested by the authors (in collaboration with our Madgraph colleagues) to work for importing both a few trivial examples of models, such as a 4-generation model, to more exotic ones, such as a model with colour-sextet diquarks and one with a baryon-number-violating vertex. Showers in such topologies are the topic of section 4.

Showers and Hadronization in Exotic Colour Topologies
In this section, we describe Pythia's treatment of QCD radiation in topologies containing colour-epsilon tensors (section 4.1) and ones involving particles with colour-sextet quantum numbers (section 4.2). This applies regardless of whether the event is generated as an internal, semi-internal, or LHEF process. We also comment briefly on hadronization aspects, pointing out relevant sources of further information.
As an aid to implementations using LHEF, a few examples of how to arrange Les Houches colour tags in colour-epsilon and colour-sextet cases are given in Figs. 1 and 1c, respectively. For completeness, we also show the status and ID codes, and the mother information, for each particle. These are explained in more detail in [22], where also more illustrations (including both standard and baryon-number violating ones) can be found.
For completeness, Fig. 1d shows a situation which Pythiais not yet capable of handling. The illustration shows a complicated baryon-number-violating cascade decay of a hypothetical fourth-generation top quark (assigned ID code 8) involving both supersymmetric, fourthgeneration, and SM particles, to produce a situation with a total of three colour-connected baryon-number-violating vertices. At the moment, Pythia's junction fragmentation model [27] is at most capable of handling up to two connected colour junctions (specifically, single junctions and junction-antijunction systems), hence only if a g → qq branching in the shower happens to break up the triple-junction system into smaller colour-singlet subsystems would Pythia's fragmentation model be able to deal with it.
A somewhat less pathological case in which multi-junction topologies may result is if a single baryon-number violating vertex becomes colour-connected to both of the junctions in the (baryon) beam remnants. This may happen some small fraction of the time through multiple parton interactions. In such cases, the following error message will be printed and the generation of the event restarted, Error in ColConfig::insert: junction topology too complicated; too many junction legs

Colour-Epsilon Topologies
For colour topologies involving the epsilon tensor in colour space (i.e., colour topologies with non-zero baryon number) we first consider the example oft →qq in the RPV-SUSY model.
The Lagrangian for the UDD-type interaction terms is . Bottom Left: production of a colour-sextet particle, assigned the fictitious ID code 6000001; the negative anti-colour tag (drawn below the sextet propagator) is interpreted as an additional (positive) colour tag. Bottom Right: A complicated baryon-number-violating cascade decay (of a hypothetical fourth-generation fermion) producing a total of three colour-connected baryon-number-violating vertices; such topologies (with three or more interconnected colour junctions) cannot yet be handled by Pythia's string fragmentation model [27]. To extract the behaviour of the radiation function, we look at the ratio of exact matrix element fort R (p 1 ) →d(p 2 )s(p 3 ) + g(q) via λ ′′ 312 to the matrix element fort R (p 1 ) →d(p 2 )s(p 3 ) and retaining only the parts that are soft-or collinear-singular (i.e., which diverge for one or more q · p i → 0). Since momentum is explicitly conserved in the shower branching process, the pre-and post-emission momenta must be related by with p 2 1 = m 2 1 =ŝ the invariant mass of the decaying squark. The Born-level matrix element squared is given by: Three diagrams (shown in Fig 2) contribute to the process where one gluon is emitted from this configuration. The matrix element corresponding to this process i.e.t R (p 1 ) → d(p 2 )s(p 3 )g(q) is denoted by M 1 and, for massless decay products (p 2 2 = p 2 3 = 0), is given by where The antenna pattern represented by equation (7) can be characterized as follows: the terms on the second line represent three soft-eikonal dipole factors (see, e.g., [2]), one for each of the three possible two-particle combinations. The factor 1/(N c − 1) in front of the dipole factors implies that the normalization of each of these eikonals is half as large as that of the eikonal term in an ordinary qq antenna, see, e.g., [28][29][30]. The two first terms on the last line of equation (7) correspond to additional purely collinear singularities for each of the quarks. The factor 1/(N c −1) is here absent; the collinear singularities have the same strength as those of an ordinary qq antenna. The X and Y terms contain subleading-color and a quasi-collinear term for the decayingt, respectively. The eikonal terms (including the leading part of the Y term, ∝ m 2 /(p 1 · q) 2 ) agree with the expression in [16,31], which is used to generate radiation for this type of colour topology in Herwig [32]. (Note that, in the Herwig implementation, the pattern is generated using ordinary full-strength radiation functions, by selecting randomly between each two-particle combination, thereby reproducing the full pattern when summing over events [16,31].) For the implementation in Pythia, we have chosen a different strategy, as follows. First, using momentum conservation, we may rewrite the antenna pattern above to only contain the final-state particle momenta, This reduces the eikonal part of expression to a single antenna between the two final-state quarks, plus subleading leftover terms. The eikonal and the collinear terms then correspond exactly to the standard radiation pattern from a qq dipole with an extra term of O( 1 Nc ). For the present work, we therefore take the radiation pattern of a standard-strength dipole spanned between the two final-state quarks as our starting point. Using s ij = 2p i · p j , this radiation function is given by [28] The Pythia showers are not based directly on equation (10), but rather on Altarelli-Parisi (AP) splitting kernels, which partition the radiation pattern onto two terms, each of which is governed by the q → qg splitting function, with z the energy fraction retained by the quark after emitting the gluon. The energy fractions of the final-state quarks, 2 and 3, are defined as in [33], The expression actually used in the Pythia showering is the sum of the AP contributions, The full matrix-element ratio, |M 1 | 2 /|M 0 | 2 , as well as the various approximate forms discussed here, are illustrated in Fig. 3, with the mass of the decayingt arbitrarily set to m 1 = 300 GeV. On the left-hand pane, we show the size of the radiation function (without the overall factor 8πα s C F ) as a function of the opening angle between the final-state gluon and one of the quarks, for a fixed (∼soft) gluon energy E g = 10 GeV. On the right-hand pane, we show the dependence on energy, for a fixed (∼collinear) opening angle θ qg = 20 • . The bottom row shows the ratio of each approximation to the matrix-element result.
The thick solid (blue) line represents the fullt * → qq + g matrix element, equation (7). For comparison, the thin solid (red) line shows the pattern obtained for a standard dipole, equation (10). The dashed (brown) curve shows the Pythia approximation to the dipole pattern, given by the sum of the AP splitting kernels in equation (13). Finally, the light dot-dashed (green) curve shows the eikonal approximation to the matrix element, used by Herwig.
In the soft limit (left-hand pane of Fig. 3), all the expressions agree in the two extremal points, in which the gluon is both soft and collinear. For wide-angle soft emissions, e.g. at 90 • opening angle, the standard dipole pattern (as well as its DGLAP variant) overestimate the full matrix element by up to a factor ∼ 1.5. That is, the Pythia shower will generate slightly too many soft wide-angle gluons. By contrast, as would be expected in the soft limit, the eikonal approximation works well for all opening angles.
In the collinear limit (right-hand pane of Fig. 3), the x axis is now the gluon energy, with the opening angle held fixed. All the expressions again agree for small gluon energies, in the double soft-and collinear limit. For intermediate gluon energies, the standard dipole pattern (as well as its DGLAP variant) again slightly overestimate the full matrix element, while they again agree with the matrix element in the hard collinear limit, on the right-hand edge of the plot. The eikonal, however, does not include the collinear-singular terms on the last line of equation (7) and hence does not reproduce the rise of the other curves in the hard-collinear limit.
In summary, our shower model will slightly overestimate the total amount of radiation, in particular at large angles, while the Herwig model underestimates hard-collinear radiation. We therefore regard the two as complementary. We note that the neglected terms could still subsequently be incorporated into Pythia 8 as a matrix-element correction [34,35], presumably mostly relevant if B-violating processes should indeed be observed in nature.
For the caseχ 0 1 → qqq, the corresponding expression is similar to equation (7) with three half-strength eikonals between the quark from the neutralino-quark-squark vertex and the two quarks from the RPV vertex [16]. Besides these, each quark has the corresponding full collinear singularity. Herwig treats this situation by randomly connecting each quark in the final state to either of the other two quarks. We have chosen instead to implement it in Pythia 8 as three genuinely half-strength dipoles spanned between the three final-state quarks.
For the case of three-body RPV gluino decay,g → qqq, only the resonant parts,g → qq * → qqq, have so far been implemented in Pythia 8, cf. section 2.4. For future off-shell contributions, the emission structure of the non-resonant parts will be obtained from the relative strengths of the intermediate off-shellg → qq * contributions.
In all cases, the subsequent hadronization phase makes use of Pythia's ability to handle string topologies including colour junctions [27], and hence issues such as baryon-number flow should be treated at least semi-realistically, allowing studies at the individual-particle level.

Colour-Sextet Particles
Within the leading-N c dipole approach to radiation adopted in Pythia, we represent a coloursextet charge as the sum of two colour-triplet charges, in much the same way as octet charges (e.g., gluons) are represented as the sum of a triplet and an antitriplet charge. Each triplet charge is independently colour-connected to an antitriplet charge. Hence a sextet may be colour-connected either by a "double bond" to an anti-sextet (in an overall singlet 66 configuration), or by two "single bonds" to two independent antitriplet charges, depending on the colour flow in the event. Each such "bond" is interpreted as an ordinary QCD dipole, with the sextet end treated as a massive quark.
At the technical level, we note that the Les Houches colour-tag standard was not originally designed to deal with sextet colour configurations. This is easy to remedy, however. Since a sextet never carries an anticolour, its anticolour tag is effectively available for use. To distinguish an additional colour (i.e., a sextet) from the ordinary anticolour (octet) case, we adopt the convention that a negative anticolour tag is interpreted as an additional colour, and vice versa for anti-sextets, as was illustrated in Fig. 1c. This appears to violate no present use of the standard (negative colour tags were so far never used in practise, as far as we are aware).
We note that a more complete treatment of the radiation and phenomenology of sextet diquarks was published while this manuscript was in preparation, see [36]. SLHA interface. This interface can therefore be used for implementation of generic BSM models without requiring a previous agreement on standardization of names and parameters.
Finally, we have commented on how Pythia 8 handles showering in non-standard colour topologies, such as the epsilon topologies encountered in R-parity violating models and in sextet di-quark ones.

Appendix B: Validated cross sections
We present validated cross sections for point SPS1a. All sparticle decays are turned off. The non-default parameters used were chosen mostly for simplicity, and to enable direct comparison with both the Pythia 6 [38,39] and Xsusy [14]