How heavy can dark matter be? Constraining colourful unitarity with SARAH

We describe the automation of the calculation of perturbative unitarity constraints including scalars that have colour charges, and its release in SARAH 4.14.4. We apply this, along with vacuum stability constraints, to a simple dark matter model with colourful mediators and interesting decays, and show how it leads to a bound on a thermal relic dark matter mass well below the classic Griest-Kamionkowski limit.

To illustrate the relationship between the two, consider 2 → 2 scattering processes from states a ≡ (i, j) to b ≡ (k, l) with matrix elements M ba and centre-of-mass momenta p a , p b . We decompose them into partial waves with where δ a (δ b ) is 1 for identical i = j(k = l) and 0 otherwise; and z the cosine of the angle between the three-momenta p a , p b . Then using unitarity of the corresponding S-matrix S ∼ 1 + iM, we find 1 2i In limiting the dark matter mass, the factor of 2 δa is compensated for non-identical particles by having two different species. These bounds should be contrasted with the typical "perturbative" ones; for example, consider a toy model dark matter candidate S with a Z 2 symmetry that annihilates to a charged scalar X via a quartic interaction: (1.7) If we consider high-energy scattering as s → ∞ then we obtain a ba 0 = − λtoy 16π √ 2 and we find the bound λ toy < 8π √ 2 at tree level. This leads to the bound (1.8) Consider now non-relativistic annihilation of the singlet S into relativistic X, so |p a | ≈ m S v, |p b | ≈ m S , s ≈ 4m 2 S , then we have the perturbative bound compared to the "absolute" bound of (1.10) Clearly even for this trivial case, for v 1 the perturbative bound is stronger and will lead to a lower limit on the DM mass, since we have taken the bound on λ toy at s → ∞ and applied it for small s. Crucially, though, this bound is really a measure of the perturbativity of the theory, since we only derived it with tree-level information, so it is entirely possible that a theory would saturate the "absolute" bound in the non-perturbative regime.
In our toy example, we included for simplicity only a quartic coupling and took s → ∞. This is rather typical in the literature among calculations of unitarity constraints. These ignore the contributions from, in particular, scalar trilinear couplings -which have enormous implications for dark matter phenomenology, since they are responsible for all s/t/u channel interactions. However, a framework within the package SARAH [53,54] for automatically calculating the constraints on scalar trilinears was introduced in ref. [38], which can automatically scan over scattering momentum to find the best limit on the couplings of the theory. This has since been applied in e.g. ref. [10-12, 14, 15, 39-41, 46]. As we saw above even in a trivial example, this will lead to generally stronger bounds on the dark matter mass than in ref. [1]. However, the calculation in ref. [38] was until now limited to colour neutral scalars. In this paper we shall describe the extension in SARAH v4.14.4 to colourful scalars, where all group theory factors are automatically calculated, and use this to place constraints on scalar trilinear couplings that are relevant for a simple dark matter model with colourful mediators.
Unitarity, however, is not the only constraint on trilinear couplings: they can also lead to alternative vacua, which in the case of charged fields mean charge-or colour-breaking minima of the potential. These are offset by having larger quartic couplings to stabilise the vacuum at the origin in field space. The typical approach to constraining a new model with such scalars, therefore, would be to use vacuum stability to constrain the size of cubic couplings, which in turn push the theory to large quartic couplings; large scattering-momentum unitarity to give an upper bound on the quartic couplings; and the dark matter annihilation cross-section is then limited by the values of both (since it can proceed via both quartic and s/t/u-channel interactions).
This reasoning is reinforced, as discussed for example in ref. [38], by the fact that for a single neutral scalar field with both cubic and quartic couplings, the full bounds from unitarity on the cubic coupling are generally less constraining that those from vacuum stability plus the upper limit on the quartic from unitarity. On the other hand, this naive picture does not necessarily hold for models with colourful states, or more scalars, but up until now there was no simple way of deriving the unitarity constraints for such theories. To our knowledge, such bounds had only been applied in a model with a colour octet in ref. [44,55,56] 2 (in the large scattering momentum limit only); and in the (N)MSSM in ref. [25,26] (with a scan over scattering momentum as discussed here) and [57] (using an earlier version of the code described in this paper). In the latter reference, a comparison of unitarity and vacuum stability bounds was performed for the Higgs-squark sector where the conclusion was that the unitarity constraints on the trilinear and quartic couplings between scalars were irrelevant in the MSSM (where the quartic couplings are given only by gauge and Yukawa couplings) but were complementary to the vacuum stability constraint in the NMSSM. However, in those models the colourful scalar sectors interact only with the Higgs scalars, which cannot provide a dark matter candidate. We also point to ref. [58], which makes use of the routines described here to constrain models of radiative fermion mass generation.
In this paper we shall investigate in detail the (genuine) complementarity of the requirements of (full) unitarity including finite momentum scattering, vacuum stability and relic density to place an upper bound on a scalar dark matter model with colourful mediators for the first time, which will allow us to put an upper bound on the dark matter mass well below the Griest-Kamionkowski limit. In section 2 we describe our model and how we have calculated vacuum stability bounds for it; in sec. 3 we describe the automatisation of the group theory calculations as we have implemented in SARAH v4.14.4; in sec. 4 we describe the procedure that we used to investigate the parameter space of our model and show the results, giving an upper bound on the mass of the dark matter particle.

A model of colourful mediators
To illustrate the new capabilities in SARAH and test the idea of a maximum dark matter mass, we shall take a model with colourful scalar mediators, but where the dark matter candidate is the usual scalar singlet S with a Z 2 symmetry. The scalar mediator fields Q E and Q O both have quantum numbers (3, 1) −1/3 under (SU (3), SU (2)) Y ; the difference between them is that Q E is even under the Z 2 , and Q O is odd. Then the most general lagrangian where the hidden sector 2 We thank Junjie Cao for bringing the first of these to our attention after the first version of this paper. respects CP symmetry is Here q i are the (3, 2) 1/6 Weyl fermions representing left-handed SM quarks. This model has several interesting features. The first, which is the main point of considering it, is the trilinear coupling κ: this entirely controls the s/t/u-channel processes for dark-matter annihiliation and is crucial for the unitarity and vacuum stability analysis. The next is the baryonic coupling Y ij Q : the mediators carry baryon number, which is respected by the model (perturbatively). It also means that the state Q E decays to pairs of quarks; we shall take it to predominantly couple to the third generation, i.e. decays to a tb pair. Therefore it is somewhat hard to search for at the LHC, being constrained mainly by ttbb searches for which no BSM reanalysis is yet possible, so we expect its mass to be only bounded to be larger than 1 TeV (rather than 2 TeV and above for other colourful scalars that decay to the first two generations of quarks). This choice also makes the model somewhat safe from direct detection constraints (provided that the Higgs portal coupling λ HS is small). In this work, we shall be considering in any case much larger masses, so collider and direct searches are not relevant.
Another interesting feature is that the state Q O can only decay to the singlet plus Q E , requiring it to be heavier than the singlet. In addition, there are three operators containing two pairs of Q O , Q E , namely the λ 7 , λ 8 and λ C terms. It is now possible within SARAH to specify all of these and for them to be properly taken into account in the unitarity constraints; however, for our analysis we shall only consider λ 7 and take λ C , λ 8 to be zero. This is mildly relevant for unitarity and vacuum stability constraints -but not at all for the dark matter density.
Since we are considering heavy dark matter that has little interaction via the Higgs portal, the relevant part of the scalar potential for this model involves the fields S, Q E and Q O . These can develop expectation values and a colour-breaking minimum if κ is large enough; however, finding the minimum of the potential involves solving coupled cubic equations and is not analytically tractable except for the the point where the masses and couplings are equal. To find possible true minima we wrote a small Python code which we briefly describe in appendix A.1. This uses HOM4PS2 [59] to quickly find all minima of the set of coupled minimisation conditions for our chosen field directions. We found this simpler than installing the no-longer-supported Vevacious [60], especially since there is a potentially large separation of scales between our dark matter sector and the Higgs sector; also note that we are only interested in the tree-level minima because we are explicitly searching for points which have large trilinear couplings where perturbativity may break down.

Colourful unitarity bounds
Unitarity bounds on colourful scattering amplitudes for the MSSM were considered in [26] where a derivation of the colour factors was given case by case for the different representations and amplitudes present. Here we shall give a description of the general procedure that we use, that applies to the scattering of any states.
Let us suppose that our initial (or final) states can be labelled A i , B j and transform non-trivially under a non-Abelian group, let us say with dimensions d A , d B . This means that we multiply the number of rows that it takes up in the scattering matrix by d A × d B . Clearly, however, we can break this into irreducible representations: where n is the total number of irreducible representations. Obviously the scattering matrix will only be non-zero when the incoming and outgoing pairs are in the same irrep, so then we need to apply a unitary transformation on the d A d B states to split them into n blocks; these are given by (generalised) Clebsch-Gordan coefficients. These can be built from invariant tensors, that is a mapping of A ⊗ B ⊗ C * → 1; we can denote this as (t C ) ij a so that A i B j C a (t C ) ij a is invariant under group tranformations. By considering infinitesimal transformations it is easy to see that contracting different invariant tensors together make another invariant tensor, and since the only invariant with just one representation and its conjugate is a Kronecker delta, then we must have 3 However, there could be more than one copy of any given representation in the decomposition above -the most relevant example here being for a product of two octet representations, for which where the relevant bit is the appearance of two 8 reps; this is more familiarly understood as the existence of two invariants, d abc and f abc , which contract the symmetric and antisymmetric combinations. Hence if we have two or more copies of a given representation, we can label them C and D and have Now we are free to diagonalise the basis of invariants and normalise them appropriately. Since the scattering matrix is an isomorphism of the initial to final colour rep, by Schur's lemma it is proportional to the identity. Then each matrix will just be d C copies of this along the diagonal. So then we need to do a unitary transformation R ij,i on the scattering matrix to split it into blocks. For it to be unitary, we need Note that the second line involves the sum over all representations present. From the above, it is clear that we can construct these matrices from our diagonalised basis of invariants, and the first condition means that we must take g CD = δ CD and R ij a = ⊕ C (t C ) ij a .
Translating this to amplitudes, for i, j → k, l we have a scattering matrix M kl ij or equivalently (a 0 ) kl ij upon which we are free to make unitary tranformations of the states to get since outgoing states are equivalent to conjugated incoming ones. So, once we have constructed the invariants, we contract them with our scattering matrices to obtain a block-diagonal form.
We now have a choice to extract a (C) 0 : we can take the trace over the remaining indices a, b, pick one example, or construct a 0 a † 0 on colour space and take the square root of the diagonal entries. In SARAH we take the simplest choice and put a = b = 1 as constraints in the evaluation of the amplitudes as it is by far the least computationally expensive. However, it should be noted that, if some of the couplings/invariants are specified by the user in a different basis, then there could in principle be a rotation between the incoming and outgoing states which would then yield incorrect results here.

Examples
The general technique that we use here is different from the approach in ref. [26], and so it is instructive to give some simple examples. We did cross-check all of the colour factors produced by the SARAH in the (N)MSSM with the results there. However, since the colour representations available in those models are not different from ours, we instead give examples directly in the model here and in appendix B.2.
Consider first our dark matter annihilation channel S, S → (Q E ) i , (Q E ) j . We can decompose the final state into a singlet and an octet, but here we can only find the singlet representation. To find the projectors we can consider the SU (N ) identity the projectors for the singlet and octet are 1 √ 3 δ i j and √ 2(T a ) i j in order for the above equation to become equation (3.6). In our model we only have t/u-channel annihilation via Q O exchange, so the diagram is proportional to κ 2 1 δ j i and so Consider now the interaction with coupling λ 5 with scattering of Q E , Q E pairs to each other. The vertex in this case is −2λ 5 (δ j i δ l k + δ l i δ j k ). So for this diagram via the singlet and octet channels we have Hence in the s → ∞ limit we have the strongest limit from the singlet representation, and a limit of λ 5 ≤ π; the same limit applies for λ 6 .
If we consider Q E , Q E scattering then we can use the same vertex, but now we decompose the representations into 3 + 6. The projector for the antisymmetric combination can be taken to be 1 √ 2 (δ 1i δ 2k − δ 2k δ 1i ) for incoming states (and (i ↔ j, k ↔ l) for outgoing) and for the symmetric one we can just take δ 1i δ 1k or equivalently 1 Hence again these give weaker bounds than the singlet representation.

Limiting the dark matter mass
Now that we have assembled the relevant machinery, in this section we will finally search for an upper bound on the dark matter mass in our model. However, to find the maximal dark matter mass with these constraints in our model with three heavy scalars, a cubic coupling and several quartic couplings involves a search on a multidimensional parameter space. We are interested in the mass hierarchy m O > m S > m E and in exploring ranges of m S up to O(300) TeV. Moreover, the quartic couplings should naively be bounded by λ S ≤ 2π 3 , λ 5,6 ≤ π. However, as seen in [38], cancellations between the contributions from quartics and cubic couplings, and the effect of a finite momentum cutoff could in principle allow somewhat larger values. Therefore, in a series of Markov Chain Monte Carlo (MCMC) scans we explored larger values with the final, finer scan of one million points having an upper limit of λ 5,6 ≤ 3.5. We chose this rather generous upper limit (instead of, say, λ 5,6 ≤ 3.2) to make sure that there are no unexpected phenomena in a theoretically excluded range. These are the most important quartic couplings since they control the overall stability of the potential. As described in appendix A.1, for our vacuum stability determination, we consider field directions along S ≡ x, Q E ≡ 1 √ 2 y, Q O ≡ 1 √ 2 z where x, y, z are real. We see that taking λ 5 = λ 6 = 4λ S renders the potential symmetric in x, y, z at large field values (when other couplings vanish) so for simplicity we impose this condition in our search which leaves us with a scan over and we simply take the other quartic couplings to zero except for λ 7 which we, quite arbitrarily, set to 0.1 although this has no impact on the search, except perhaps a very slight influence on vacuum stability. In other words, we are allowing self-couplings of the mediators and the singlet, respectively, and some coupling between the mediators. On the other hand, we are ignoring quartic couplings among the singlet and the mediators, and those where a Higgs boson is involved, since we are interested in the model with a t/u-channel mediator and not in the quartic quartic coupling channel -or as a Higgs-portal model which have been extensively studied in the literature and has larger direct detection prospects. Moreover, for simplicity we take Y 33 Q = 1 and zero for other Baryonic couplings, so that our Q E field only decays to a top and bottom quark pair. This leaves us with five parameters, four of which are dimensionful. In principle λ 2 would also have an important impact on the annihilation of singlets to mediators, while changing the relationships of the quartic couplings may have some impact on the stability results. In future it would be interesting to perform a more sophisticated scan to allow for a more high-dimensional parameter space.
To explore our parameter space, we performed a series of scans, starting from a uniform grid, then implementing several parallel Markov Chain Monte Carlo scans via the Metropolis-Hastings algorithm distributed across multiple cores on a cluster. Since we are interested only in the upper bound on the singlet mass, we construct a likelihood function L as a product: where the first three likelihoods are sigmoid functions that cut off smoothly above the upper limit: and δ stability is 1 for a stable vacuum and 0 otherwise. This amounts to fixed large bias for stable over unstable points 4 without categorically excluding unstable points. The second term of the combined likelihood corresponds to the unitarity constraint. The last term is a bias on the dark matter mass, forcing the scan to probe heavier singlets: The value of m S differs depending on the scan. After completing the MCMC scans, we select the points of the sample that strictly satisfy our constraints, which are therefore imposed as "hard cuts". Employing MCMC scans bears the advantage that a valid parameter space can be proposed more efficiently than a grid or random scan because the latter focus on regions that are allowed and avoid wasting computational resources on regions that are clearly excluded. In all MCMC scans, we select the largest partial wave amplitude to get a "good" point.  Figure 1a shows the distribution of the singlet mass after our scan, including only those points which passed all cuts. In table 2, we list the amount of points that pass after each combination of cuts. Hereby, the cut on the mass hierarchy ensures that m S ≤ m O , and that λ S ≥ 0.5. There is a clear cutoff at around m S 47 TeV, after which we found no more valid points. This implies a considerable amount of this mass range could be covered with a 100-TeV-collider. This is also the central result of this paper.  Table 2: Points left over after each cut. The raw sample contained one million points. D refers to the cut on the dark matter density, U to that on unitarity, and V to that on vacuum stability. Details see text.
We show in figures 2a-2c the effect of the separate cuts on the remaining points on the parameters m S , κ and Λ. We see that Λ is bound by the naive unitarity constraint of π.
As was expected when setting up the model, we find a pretty clear relation between the strength of the coupling κ and the masses of the involved particles. This can be seen in figure 1b. There is a clear peak around 3.5 for κ/m S , although there are some outliers towards higher values. If instead One can see that the cutoff at Λ ∼ π is due to unitarity. The y-axis shows, in all three plots, how many of one million scan points made it through each cut, respectively. In contrast to figure 1, these plots do not contain any information about which points make it through two or more cuts.
we take κ in relation to the largest mass of each datapoint, i.e., one chooses the largest out of m S and m O , the outliers disappear ( figure 1c). Instead, we find a peak at a ratio of about around 2.5. Figure 3a shows the valid points in the κ − m S -plane after each individual cut. One can see a clear correlation between the two, and the peak of κ/m S at 3.5 (figure 1b) is manifest. The outliers with a higher κ/m S ratio tend to be concentrated around the lower end of the distribution, where κ is around 50 TeV and m S is below 10 TeV. One can see that the vacuum stability constrains the allowed area from the bottom, i.e., the valid points are situated above the diagonal passing through (κ, m S ) = (100TeV, 10TeV) and (150TeV, 30TeV). Likewise, the dark matter criterion constrains the allowed area from the top, i.e., the valid points are below the diagonal passing through (50TeV, 25TeV) and (150TeV, 45TeV).
In figure 3b, one can see the distribution of valid points in the Λ − m S -plane after every cut.Vacuum stability eliminates points with low values of Λ or m S . The dark matter cut, by itself does not have much impact on the shape of the distribution. As expected, the cutoff Λ ≤ π is ensured by unitarity (third panel of figure 3b). After all cuts, the points in the lower m S range are excluded, as expected, but also those above the diagonal passing through (Λ, m S ) = (1.5, 25TeV) and (3.0, 50TeV). The latter is a compound effect from the cuts on dark matter and unitarity, which shows that the dark matter cut does play a role after all.
Finally, figure 3c shows the distribution of valid points in the κ − m O -plane. After each of the individual cuts, the resulting shape is bordered by three diagonals: one almost vertical one on the low-κ end, and two more or less parallel ones going from the bottom left to the top right of the respective panel. The distribution of valid points after all cuts can be deduced almost directly from the overlap of the distributions after the three individual cuts.
Finally, figure 4 shows the distribution of κ/m max as a function of Λ after various cuts. One can see that vacuum stability imposes κ/m max Λ + 1. Unitarity cuts away at some of the higher values of Λ and κ/m max , and cuts off at Λ ≤ π.

Highest singlet mass
The point that we found where the singlet was heaviest, and the main result of this section, has  Figure 3: Left column: selected two-dimensional planes of the parameter space after cutting for vacuum stability. Second column: same, but cut for dark matter. Third column: same, but cut for unitarity. Right column: same, but after all three cuts. Top: distribution of κ against m S . Middle: distribution of Λ against m S . There is a clear cutoff at Λ π due to unitarity. Bottom: distribution of κ against m O , after each cut. As in 2c, one sees that the cutoff at about π is due to unitarity. The coloured regions indicate the regions where there were valid points after each cut, normalised for each plot (so a direct comparison of the colours between plots is not possible). earlier discussion), evaluated at √ s = 141 TeV, well away from any poles. This point is on the cusp of being ruled out by the unitarity calculation, which is dominated by the coupling κ; we find that decreasing the coupling Λ changes a 0 very little at this point but leads to an unstable vacuum already at Λ = 3, while increasing κ to 180 TeV leads to a 0 > 0.5 (and also an unstable vacuum).

Trilinears excluded by unitarity alone
Finally we wish to highlight that, although most points conformed to the naive expectation that we could apply the limit Λ < π from unitarity and constrain κ just from vacuum stability, there are exceptions that underline the complementarity of the unitarity calculation. For example, This point has Ωh 2 = 0.12 and maximum a 0 = 0.51 (again from the singlet submatrix) and the vacuum stability equations have no other solutions than the origin. In fact, this point is typical of a whole branch of points where m S ∼ m O m E for which this is true -these points are excluded by unitarity because of the size of κ, but we would not have seen this either from classic unitarity bounds where s → ∞ or from the vacuum stability constraints.

Conclusions
We have described the calculation and implementation of constraints from unitarity of scattering for 2 → 2 processes involving scalars of any representation under the strong gauge group, and finite scattering momentum. Since these unitarity constraints automatically constrain all the scalar couplings of a theory, they are now very straightforward to include for a whole new class of models.
We also illustrated the utility of these routines and the complementarity of the information that they provide for studying dark matter models compared to vacuum stability and both naive infinite momentum perturbative unitarity constraints, and the "absolute" bound of ref. [1]. We showed that there are points for which vacuum stability and "naive" unitarity are insufficient, i.e. the full perturbative unitarity calculation is indispensable. We introduced a toy model with a baryonic coupling and colourful mediators that decay in an interesting way to a top-bottom quark pair, that is a very simple example of the sort of models that can now be explored with these constraints. It would be very interesting to explore models with more complicated gauge representations.
The work also paves the way for several further extensions in future work: additional unbroken gauge groups; fermions and/or vectors in the scattering matrix; loop corrections. Moreover, our dark matter model had a maximum mass of 47 TeV, and coupled to colourful states, so much of the allowed parameter space would be accessible to a future 100 TeV. It would therefore be interesting to consider dark matter-collider complementarity in terms of both its signatures at such a collider; but also the low-mass bounds at the LHC, since it could be searched for in the ttbb channel.
ing the colour factors in the unitarity routines in SARAH. MDG also thanks Michael Baker for correspondence about those routines, and helping to identify bugs in the beta-version.
A SARAH implementation of our model: SM-SQQ Since we implemented our model in SARAH we list here the relevant parts of the model file (which is now also made public with version v4.14. where the last line is the Z 2 symmetry charge. So QP, QM correspond to the fields Q E , Q O in the body of the paper respectively. The Lagrangian is then given by the terms We did not explicitly give the colour structure for Lambda7 and hence we do not also include Lambda8 which has the same fields but a different contraction of the indices.

A.1 Vacuum stability calculation
Here we describe our routines for computing the vacuum stability constraints. We begin with the potential for the fields S, z for x, y, z real and the other components of Q E , Q O zero (since this is the most unstable direction in field space). This yields a potential and then take the derivatives. These give three equations for which all the solutions can be found with HOM4PS2. However, we first rescale all of the dimensionful terms by We then calculate the numerical value of the potential at each of the solutions that we find and, if the minimum value is not at the origin of field space, we note that the vacuum is not stable.

B New routines in SARAH
With the release of version 4.14.4, SARAH contains updated routines to calculate unitarity constraints for scalars including colourful states (for now, no other unbroken non-abelian groups are considered for the unitarity routines). The algorithm used is as described in section 3, with the group invariants hard-coded for certain common representations (in particular for octet representations the f abc and d abc matrices) but otherwise calculated by the included routines from Susyno [66]. From the point of view of the user, the calculation functions exactly as in ref. [38] except that colourful scalars are automatically included, unless they are explicitly removed from the scattering (as is done by default for many models). However, some new features have been added to aid performance/use/inspection of the results which will be described here.

B.1 A new option for the cut-level of poles
In [38], several settings for the unitarity routines were outlined. For completeness, and to correct a misprint there, we give here the correct and updated complete options: The relative proximity to s-channel poles that is allowed, C s 5 .
With the new version, setting SPhenoInput 446 to 4 will cause the unitarity constraints to be disregarded for the value of s whenever a pole is found in an s, t or u channel of any diagram in any scattering submatrix (the program will continue to scan over the range of values of s in the hope of finding a valid constraint). This is the choice made in e.g. [28] and is the most conservative condition that can be placed, especially if coupled with a large √ s min .

B.2 Storage of symbolic form for each diagram
As part of the upgrade to the unitarity routines, couplings having more than one colour structure are properly taken into account and stored in new routines within SPhenoCouplings.f90. While not all functionality of SARAH will handle such couplings correctly yet (notably the loop decays) the unitarity routines and spectrum generation will give correct results.
To enable cross-checks and reuse of the new routines, SARAH writes a file Unitarity.m which contains symbolic information about each scattering diagram computed. The format is {{ s1 , s2 , s3 , s4 , prop , Type , c o l o u r r e p , { dyn1 , dyn2 , dyn3 , dyn4 } } , { c o u p l i n g s }} prop is the field appearing in the propagator (or just 1 for a quartic coupling); Type is one of Q, S, T,U meaning quartic, s/t/u-channel; dyn1, dyn2, dyn3, dyn4 are the dynkin indices of the fields s1,s2,s3,s4. The fields are given as incoming states, and the ordering is such that for s-channel and quartic interactions s1,s2 → s3,s4 while for the t-channel it is s1,s3 → s2,s4 and for u-channel it is s1,s3 → s4,s2. Some typical lines for the model described in this paper in the SARAH notation of appendix A would be The third and fourth lines show a quartic coupling, so prop is given as 1 (since there is no propagator). However, the colour representations are 3 and 6 respectively, and the two possible colour structures of the quartic coupling are involved. Recall that the lagrangian term is −λ 5 |Q E | 4 so the two structures are δ ij δ kl and δ il δ kj , both with coupling −2λ 5 in this case. Hence we the net result for the antisymmetric term is 0, and for the symmetric one it is −4λ 5 , exactly as we found in section 3.1.
One final note about the routines, to avoid confusion, is that, in SARAH we actually calculate the matrix for −a 0 and then take the absolute values of the eigenvalues.

B.3 Splitting into CP eigenstates
With the inclusion of more fields, the program must find the eigenvalues of a larger scattering matrix and it is desirable to find simplifications where possible. In much the same way that we use the representations under charge and the strong force to decompose into scattering blocks, we can also use CP to reduce the rank of our matrices. If the user places the line UNITARITYCP=True ; in the file SPheno.m for the model, SARAH will attempt to assign CP charges for the states and decompose the scattering matrices accordingly. The user should find that the result is entirely unchanged, but for more complicated models some performance improvement may be found.