Unitarity constraints on general scalar couplings with SARAH

We present an update of the Mathematica package SARAH to calculate unitarity constraints in BSM models. The new functions can perform an analytical and numerical calculation of the two-particle scattering matrix of (uncoloured) scalars. We do not make use of the simplifying assumption of a very large scattering energy, but include all contributions which could become important at small energies above the weak scale. This allows us to constrain trilinear scalar couplings. However, it can also modify (weakening or strengthening) the constraints on quartic couplings, which we show via the example of a singlet extended Standard Model.


Introduction
In a classic paper, Lee, Quigg and Thacker showed that the Higgs mass in the Standard Model (SM) must be below 1 TeV in order to maintain perturbative unitarity [1]. From the measurement of the Higgs mass at the LHC [2,3], we have learned that the quartic coupling in the SM is even well below 1, i.e. the scalar sector of the SM has very weak self-couplings. However, this is not necessarily true if one adds more fundamental scalars to the theory. The scalar potential of BSM models often involve many new parameters which are experimentally barely constrained. Therefore, theoretical conditions like the stability of the potential or the conservation of unitarity are very important to find physical viable parameter regions in these models. The constraints from tree-level perturbativity are often applied in well studied models such as ones with additional singlets [4][5][6], doublets [7][8][9][10][11][12][13], or triplets [14][15][16]. However, a e-mail: florian.staub@kit.edu very often the constraints are derived under the assumption that the scattering energy is much larger than the involved masses. In this limit, only point interactions are important, and all diagrams with propagators are neglected. As a consequence, cubic couplings do not enter the widely used constraints at all. However, it was already pointed out long ago that large cubic couplings of light scalars are dangerous [17], and also limits for the triple Higgs coupling in the SM were derived [18,19]. In this work, we shall present a general calculation of unitarity constraints without the assumption of a very large scattering energy, with the following motivation: 1. We want to place bounds on genuine trilinear couplings. 2. For theories where additional scalars couple to the Higgs, even if there are no trilinear couplings before electroweak symmetry breaking, they are generated after the Higgs takes a vev, and unitarity of scattering at finite s gives new constraints on these quartic couplings. 3. For theories defined with a low cutoff, scattering may never be in the regime where the energies are sufficiently large to neglect the s, t, u-channel processes. 4. Even for theories with a high cutoff, the infinite energy approximation is rarely justified since the couplings must run: if we take the energy sufficiently high to be able to neglect particle masses, the resummed couplings will typically have completely different values.
For this purpose we have extended the Mathematica package SARAH [20][21][22][23][24] by routines for the analytical and numerical study of the full tree-level unitarity constraints. While the analytical routines are helpful to obtain expressions for 2 → 2 scattering elements, a symbolic calculation of the full scalar scattering matrix could become slow and less illuminating. Therefore, for practical application the Fortran output for SPheno [25,26] has been also extended to obtain a numerical prediction for the maximal eigenvalue of the full scattering matrix.
We discuss in sect. 2 the underlying calculations to obtain unitarity constraints in generic BSM models, and the assumptions/restrictions that we shall apply. The importance of the full calculation is demonstrated in sect. 3 via the example of singlet extensions of the SM. In sect. 4 we show how the new routines are used. A brief summary is given in sect. 5.
2 Generic calculation of unitarity constraints 2.1 2 → 2 Scattering processes of scalars at finite momentum The derivation of unitarity constraints is elementary, but the derivation for finite momentum is rarely found in the literature -and there are many common misunderstandings -so we present a clear exposition in Appendix A. The result is that the partial wave constraint becomes where a J is a normal matrix related to the partial wave decomposition of 2 → 2 scattering matrix elements M ba from a scattering of a pair of particles The factor δ 12 (δ 34 ) is 1 if particles {1, 2}({3, 4}) are identical, and zero otherwise. P J are the Legendre polynomials, p i is the centre of mass three-momentum for particle i, and s = ( p 1 + p 2 ) 2 is the standard Mandelstam variable.
From the fundamental equation (1) different constraints can be derived; we shall only consider the zeroth partial wave, and denoting a i 0 as the eigenvalues of a 0 we shall apply The diagrams which contribute to 2 → 2 scalar scattering processes are shown in Fig. 1. For a general field theory consisting of real scalars φ i and couplings The integration over cos θ is trivial for the contact and schannel processes, and always straightforward for the others using where E i are the energies of the particles in the centre of mass frame, and p 1 , p 3 are the three-momenta. We shall express the results in terms of the function so that allowing us to define In terms of these, the (modified) zeroth partial waves are

Handling of poles
In the neighbourhood of poles, the tree-level amplitude diverges which signals that we need to take higher-order corrections into account (which will effectively modify the divergent propagator to include the width, cutting off the divergence). Moreover, since we are only calculating unitarity constraints at tree-level, in the presence of large couplings large quantum corrections to masses may mean that the physical location of the poles is a long way away from the tree-level mass parameters. Both of these issues imply that we should not trust our results in such cases. We therefore apply the following conditions: 1. s-Channel poles Obviously, s-channel poles are present if any propagator mass is close to √ s. In order to cut out this region, we set the entire irreducible scattering matrix to zero if the condition is violated. By default, a value of 0.25 is used for the cut variable C S . C S can also be changed by the user. This condition is simpler than the one proposed in Ref. [27] which calculates the width Γ of a particle, and puts the condition | √ s − m| > amΓ (Q = bm) using suitable coefficients a, b 1. However, the impact on the scattering amplitude is expected to be small. Only for models with scalars which have a very large it might be necessary to adjust C S . 2. t-/u-Channel poles Particles in the t and u channels can become on-shell. For a t-channel diagram, this can happen if holds. Similar conditions exist, for 1 ↔ 2 and 3 ↔ 4. The conditions for u-channels are obtained by exchanging 3 ↔ 4. These conditions (used in [27]) are only necessary to have a pole, but not sufficient -they are too conservative. In fact, the presence of such a pole also demands that the scattering energy s is smaller than a given value. The general conditions for the minimal scattering energy s min to avoid poles are From this we find that for an often appearing, kinematic configuration with m 3 = m 1 and m 4 = m 2 , a t-channel pole only shows up for We will include three different treatments of such poles, which the user can select depending on taste: (a) Only the matrix element for which such a poles appears is set to zero, but all other entries of the scattering matrix are kept. This gives the most aggressive limits. (b) A partial diagonalisation of the scattering matrix is performed as proposed in Ref. [27]: assume that X is the set of all kinematically accesible states at a given energy √ s and Y is the subset which involve states that suffer from a t or u channel pole. In that case, one can diagonalise the scattering matrix for the set X \Y with a unitarity matrix U to obtain the eigenvalues a i,X \Y 0 . The condition equation (3) is changed to (c) The entire irreducible scattering matrix is set to zero. This gives the weakest limits.

The role of the Goldstone boson equivalence theorem
To apply unitarity constraints, in principle we should consider all coupled channels for all particles. However, in practice there are a large number available, most of which will not contribute in a meaningful way to constraints, and so in the interest of computational speed it is necessary to impose some simplifying assumptions. These are: 1. We can neglect all contributions proportional to gauge couplings, and scatter at energies well above the mass of any gauge bosons; clearly for smaller scattering energies this would mean we would be in the neighbourhood of an abundance of poles. Furthermore, light bosons mediate infra-red unsafe scattering, so our above formalism would require modification, and it is therefore reasonable to eliminate them. 2. We neglect all fermionic contributions. The above assumption partly justifies this, as the contributions to scattering from Standard Model fermions should be small at energies well above their masses. 3. To avoid an abundance of group structures of the scattering pairs, we do not consider any particles that transform under any unbroken symmetries except for the electric charge. In particular, this excludes any strongly coupled particles (such as top partners).
Assumption (1) is the most reasonable, and also most powerful: since amplitudes involving transverse gauge bosons are always proportional to gauge couplings, we can neglect them. For longitudinal gauge bosons of mass m V , whose polarisation vectors can be taken to be the scattering amplitudes contain factors of 1/m V and hence inverse powers of the gauge couplings, so that they can have a finite amplitude as the gauge couplings are taken to zero. On the other hand, since we scatter at energies well above their masses, the Goldstone Boson equivalence theorem allows us to instead replace all external longitudinal gauge bosons with the Goldstone boson, with the important proviso that it has a physical mass equal to the gauge boson mass (so not equal to ξ m V ). It then turns out that, since we neglect contributions proportional to the gauge couplings, we can also neglect gauge boson propagators -but only if we work in Feynman gauge; we discuss in appendix B why this is so and what happens in other gauges. Taken together, then, the above assumptions, and working in Feynman gauge, enable us to consider scattering amplitudes where all states are scalars. While it would be an interesting if time-consuming task to relax some of these assumptions (which we leave to future work), they are already very powerful and allow us to study a wide range of theories.

Examples: singlet extentions of the Standard Model
We want to demonstrate the importance of the unitarity constraints beyond the large s approximation by a brief example.

Pure singlet model
First we shall consider the simplest possible BSM model: the SM extended by a real singlet S. To illustrate point (1) in the introduction, if we just consider the singlet and assume that its couplings to the Higgs sector are small relative to its self-couplings, we can take the Lagrangian: There are two additional minima away from the origin if but the origin remains the true minimum if As a probe of genuine trilinear couplings, taking the minimum at the origin is most interesting, because once the singlet obtains an expectation value, stability constraints the trilinears to be rather small compared to the physical mass. It is simple to derive a 0 for this case: We show the constraints on this for √ s = 4000 GeV, m S = 500 and 1000 GeV (or other arbitrary units) in Fig. 2. Clearly the constraint from s → ∞ would give λ S < 4π behaviour of |a 0 | as √ s is varied above threshold on the right-hand plot: there is always a rapid increase followed by logarithmic behaviour.
The finite s constraints clearly consist of two different regimes that intersect, and come from the fact that the t/u channel contribution has opposite sign to the s-channel and quartic term. As a function of s, a 0 grows sharply from 0 at s = 4m 2 S , before decreasing again and tending to the large s value. So the constraints come both from the maximum allowed s, and around s = 6m 2 S ; for λ S = 0 it occurs at s 5.6m 2 S . From the maximum s, we obtain the curved regions in the plot that have a minimum value for λ S as it passes through κ = 0. These therefore show a difference when we change m S . On the other hand, the overlapping curves that pass into the unstable region |κ/m S | > 3 √ λ S come from taking s near 6m 2 S . Note that the value s = 5.6m 2 S is not near any pole value, and the corrections to the singlet mass are well under control; at one loop they are so for κ = 5m S , we have δm 2 S ∼ 0.3m 2 S . Moreover, the scattering energy is sufficiently large that the produced particles are relativistic, so we are not in a regime where e.g. Sommerfeld enhancements would play a significant role. Hence the enhancement to the partial wave amplitude is an effect that we can use to constrain the couplings of the theory. In particular, it gives an upper bound on κ for which trustworthy results are calculable, independent of vacuum stability considerations -especially for larger values of λ S . We expect this to be a general feature: from an inspection of the righthand plot of Fig. 2 we see that the strongest limits (away from poles) to a given model will either come from near-threshold production or at large s.
This model also allows us to simply illustrate our points (3) and (4) in the introduction. What we are interested in constraining are the values of κ/m S and λ S at low energies. However, the partial waves receive quantum corrections which can be very significant for large couplings, and if the scattering energy is large, we should certainly resum the logarithms and place constraints only on couplings evaluated at a renormalisation scale of √ s (see e.g. [28]). In this model, the one-loop β-function for the quartic coupling gives which can be solved exactly, and gives a Landau pole at For λ S (500 GeV) = 4, this is at 1500 GeV! Hence we cannot apply the infinite-energy scattering limit to this coupling. Put another way, since we must understand the limits in Fig. 2 to be evaluated at μ = √ s, if λ max S (4000 GeV) = 4, then λ max S (500 GeV) = 1.4.
3.2 Singlet extended SM with conserved Z 2 While the above model is trivial, it contains most of the ingredients that we find in more complicated models, in particular the partial cancellation between the channels. Now we will turn to a more physical example: a singlet that couples to the Higgs, but with a Z 2 symmetry which stabilises it and prevents mixing with the Higgs. The potential reads This theory contains no trilinear scalar couplings before electroweak symmetry breaking. As such, it is a useful prototype of popular extensions of the SM such as the Two Higgs Doublet Model, NMSSM, etc, as well as being phenomenologically interesting in its own right (for example, it provides a dark matter candidate). However, once the Higgs obains an expectation value so that we can write the neutral Higgs is generated. Thus we will have s, t, u-channel scattering processes in the scalar sector which will modify the unitarity constraints! In the large s limit we have Max |λ H S | , |λ H | , We want to compare this with the full calculation. Results for the scattering processes are already given in literature [4,5], but we disagree with both references in different channels. Therefore, we list all matrix elements in appendix C. In the following, analytical discussion, we concentrate only on the parts involving CP-even states. The scattering matrix involving only the the Higgs and the singlet is If we assume for the moment λ H S λ, λ S , the dominant contribution is the hS → hS scattering.
The result reads 16πa 0 (hS → hS) In order to simplify this expression we consider the limit of small m S and large v 2 λ 2 Thus, for s ∼ m 2 h , this scales as which can be significantly larger than the limit from point interactions only. This is also confirmed by our numerical calculation with SPheno. In Fig. 3 we compare the limits from including point interaction only with the full calculation in the (λ S , λ H S ) plane for different singlet masses. We see that the unitarity limits on λ H S become much stronger for m S < m h and small λ S . However, even for larger masses a pronounced effect is visible. Even for m S = 500 GeV, the limits are stronger by a factor of two. This brief example demonstrates the importance of going beyond the large s approximation when considering unitarity constraints in BSM models. Detailed discussions of these effects in other, phenomenologically more interesting models will be given elsewhere [29,30].

Implementation in SARAH
We have extended the Mathematica package by the results and procedures summarised in the previous sections; in particular, the restrictions/assumptions that we apply are described in section 2.3. The user has two possibilities to use the new functionality: (i) during the Mathematica session analytical expressions for specific scattering processes or the full scattering matrix are available; (ii) the necessary routines for a numerical calculation of the unitarity constraints

Commands
In order to obtain analytical expressions for 2 → 2 scattering processes or the entire scattering matrix new commands are available in SARAH with version 4.13.0. G e t S c a t t e r i n g D i a g r a m s Here, incoming1, incoming2 are the incoming particles and outgoing1, outgoing2 the outgoing ones. One needs to use for these variables the names of fields in SARAH. Optionally, a generation index can also be given.

§ ¤
In [1] G e t S c a t t e r i n g D i a g r a m s [{ hh , hh } -> { Ah , Ah }] ¦ ¥ This returns the scattering element for hh → A A. If the scalar h and pseudo-scalar Ah appears in several generations in the given model, the indices in1, in2, out1, out2 are used. -Explicit generation indices, e.g.

§ ¤
In [1] G e t S c a t t e r i n g D i a g r a m s [{ hh [1] , hh [1]} -> { Ah [2] , This sets the generation indices of the incoming fields to 1 and of the outgoing fields to 2.
The result of GetScatteringDiagrams is a function of the couplings and masses in the model. In addition, keywords are introduced to make it possible to trace back the origin of the different terms: s, t and u-channel diagrams as well as point interactions are multiplied with a variable: -sChan for a s-channel diagram -tChan for a t-channel diagram -uChan for a u-channel diagram -qChan for quartic interactions Thus, one can easily remove specific diagrams in order to check their impact by setting the corresponding variables to zero. 3. Scattering matrix: the full scattering matrix is return by running § ¤ In [1] B u i l d S c a t t e r i n g M a t r i x ¦ ¥ All generation indices of the external fields are explicitly inserted.

Example
We show via the example of the SM how the new commands are used in practice. First, SARAH needs to be loaded and the SM be initialised: § ¤ In [1] << SARAH -4.13.0/ SARAH . m .
In [2] Start [ " SM " ]; ¦ ¥ Afterwards, one can start to play with the unitarity constraints. Here, we want to replace the quartic coupling λ, which is usually used in the vertices, by the Higgs mass. That's done during the initialisation process of the unitarity constraints.

§ ¤
In [1] InitUnitarity Here, we make here the same assumptions as above. Moreover, we set all masses but the one of the CP even Higgs to zero, i.e. we take the limit m Z = m H + = 0. The outcome is: § ¤ Out [1] {{( mh2 ((2 mh2 -3  where we only have shown the first two rows. In order to see the basis in which the matrix is given, one can check § ¤ In [1] s c a t t e r i n g P a i r s ¦ ¥ which reads in our case § ¤

Including the unitarity constraints in the SPheno Output
While it might be helpful to obtain analytical expressions for some specific channels, in practice a numerical calculation is often more useful. However, Mathematica is not the preferred environment for exhaustive, numerical calculations.
Therefore, it was natural to extend the existing SPheno output of SARAH by the new function. So far, SARAH is already producing Fortran source code which can be compiled with SPheno. This provides the possibility to calculate many things for a given model very quickly, e.g. two-loop RGEs, one-and two-loop masses [31][32][33], flavour and precision constraints [34], two-and three-body decays at tree-level, loop corrections to two-body decays [35], and so on.

Generating the Fortran code
The properties of the new spectrum generator based on SPheno are defined within SARAH by using the input file SPheno.m. SPheno.m contains for instance the information about the free input parameters expected from the user, the boundary conditions at different scales, choices for involved scales and several other settings. With SARAH 4.13.0 the following settings are supported: AddTreeLevelUnitarityLimits=True;

¦ ¥
This enables the output of all routines to calculate the treelevel unitarity constraints. By default, this generates the full scattering matrix involving all scalar fields in the model which are colourless. In the case that some particles should not be included, they can be explicitly removed via SPheno.m § ¤ 1 RemoveParticlesFromScattering={Se,Sv};

¦ ¥
Here, we have for instance decided not to include charged and neutral sleptons in the case of a supersymmetric model. Once SPheno.m for a given model has been edited, one can proceed as usual to obtain the source code and compile:

Run § ¤
In [1] MakeSPheno [] ¦ ¥ to obtain the source code 2. Copy the code to a new SPheno sub-directory For the last step, a Les Houches input file must be provided which includes the numerical values for the input parameters as well as settings for SPheno.

Configuring the unitarity calculations
If the unitarity constraints are turned on the in the SPheno output, several new settings in the Les Houches input file are available: 2 : what is the value for √ s at which the scattering is maximised 11-13 : this repeats the input for √ s min , √ s max and the number of steps.

Summary
We have presented an extension of the Mathematica package SARAH to calculate unitarity constraints in BSM models. It is now possible to obtain predictions for the maximal element of the scattering matrix in a wide range of models without making use of the large s approximation. We have provided generic expressions for the calculations, along with pedagogical derivations, and clarified some technical issues concerning additional gauge bosons and the choice of gauge. We have briefly shown the importance of these improved constraints in the example of the real singlet extended SM. More detailed discussions of the effects of the new constraints in doublet and triplet extensions will be given elsewhere [29,30]. 02, ANR-10-LABX-63). We would like to thank Sophie Williamson and Manuel Krauss for helpful discussions and collaboration on related topics.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

A: Derivation of the partial wave unitarity constraint
In this appendix we will present an elementary derivation of unitarity constraints, retaining finite momentum factors that are less widely known (and absent from e.g. [5]).
First, we define the S-matrix in terms of the interaction matrix T as S = 1 + i T . Then in terms of matrix elements of scattering from (multiparticle) states a with a set of momentum { p} to a set of states b with a set of momenta {k} we have and so Now S must be a unitary matrix, and so the the constraints from unitarity come from Then we insert a complete set of states to evaluate T † T : {k, b}|T † T |{ p, a} Now specialising to the case of 2 → 2 scattering, we can rewrite the equation as Here δ c = 0 if the particles in c are not identical, and 1 if they are identical, to allow us to keep the same phase space region of integration (otherwise we double count), and p c is the three-momentum in the centre of mass frame for the pair c. Now for the partial wave analysis, we define the threevectors vectors p a , k b , p c to lie along the unit vectorŝ k a = (1, 0, 0) We decompose the matrices into partial waves: where P J are the Legendre polynomials, satisfying Next we require the identity where the spherical harmonics satisfy Y J m ∝ e imφ P m J (cos θ), Y J 0 = 2J + 1 4π P J (cos θ).

B: Scattering amplitudes and partial waves from gauge boson propagators
In this appendix we will clarify the fate of scattering amplitudes among scalars where there is a gauge boson propagator. We neglect all contributions to the final amplitude that are proportional to gauge couplings, but if we work away from the Feynman gauge (it may be desirable to define a theory in that way) then the vector propagators have factors of 1/m 2 V where m V is the vector boson mass -which is proportional to the gauge couplings. In which case we would necessarily need to included gauge boson amplitudes as well as the Goldstone bosons. We write the couplings of a massive gauge boson to real scalars φ i as To these, we should add the contributions from the Goldstone bosons. In [35] it was shown that, for scalars coupling to the corresponding golstone boson with the same index V that the couplings are related by Hence when we add the contribution from the Goldstone bosons, we just obtain M ba (Vector propagators + Goldstones)

(B.21)
This result is manifestly gauge invariant. Setting the gauge couplings g V 12 to zero, we have a remaining piece which is just equal to the contribution from the Goldstone bosons in Feynman gauge! However, we should note that in other gauges it is necessary to include the gauge boson propagators; for example, if we work in unitary gauge then there are no Goldstone boson propagators! As an aside, if we want to include the contributions from heavy gauge bosons (i.e. not neglect their couplings), it is simple to perform the angular integrations for these contributions, using: .