Improved unitarity constraints in Two-Higgs-Doublet-Models

Two-Higgs-Doublet-Models (THDMs) are among the simplest extensions of the standard model and are intensively studied in the literature. Using on-shell parameters such as the masses of the additional scalars as input, corresponds often to large quartic couplings in the underlying Lagrangian. Therefore, it is important to check if these couplings are for instance in agreement with perturbative unitarity. The common approach for doing this check is to consider the two-particle scattering matrix of scalars in the large centre-of-mass energy limit where only point interactions contribute. We show that this is not always a valid approximation: the full calculation including all tree-level contributions at finite energy can lead to much more stringent constraints. We show how the allowed regions in the parameter space are affected. In particular, the light Higgs window with a second Higgs below 125 GeV completely closes for large values of the $Z_2$ breaking parameter $|M_{12}|$. We also compare against the loop corrected constraints, which use also the large $\sqrt{s}$ approximation, and find that (effective) cubic couplings are often more important than radiative corrections.


I. INTRODUCTION
The discovery of a scalar boson at the Large Hadron Collider with a mass of around 125 GeV was a milestone for particle physics [1,2]. This state has all expected properties of the long searched-for Higgs boson, and all particles predicted by the standard model of particle physics (SM) have finally been found. Even if no additional, fundamental scalar has been observed so far at the LHC, it is much too early to give up the possibility that more Higgs-bosons exist which are involved in electroweak symmetry breaking (EWSB). There are several possibilities what the origin and the properties of such states could be. A very attractive and well studied scenario is that a second Higgs doublet exists. After EWSB, the two Higgs doublets yield one particle which has all the properties of the discovered state, but they also predict the presence of one charged and two neutral additional bosons. There exist several constraints on this kind of models: the LHC measurements must be reproduced, the absence of any other signal must be explained, including modifications to rare decay processes. From the theoretical point of view, these models are usually confronted with two conditions: (i) the electroweak vacuum must be stable or at least sufficiently long-lived [3][4][5][6][7][8][9][10][11][12], (ii) unitarity should not be violated [13][14][15][16][17][18][19]. In order to probe unitarity in BSM models, the standard procedure in the literature is to calculate the scattering matrix for 2 → 2 processes involving scalars. Usually, only point interactions are included, which do not vanish for very large scattering energies √ s. For extensions of the Standard Model, the contributions from scalar trilinear couplings * goodsell@lpthe.jussieu.fr † florian.staub@kit.edu have only been considered for singlet extensions and the minimal supersymmetric standard model [20][21][22]. Therefore, it is time to check if the large √ s approximation in THDMs is valid or under which circumstances it might give misleading results. This letter is organised as follows: we show our conventions for THDMs in sec. II, before we briefly summarise our approach to calculate the tree-level unitarity constraints in sec. III. The impact on the parameter space is discussed in sec. IV. In sec. V, we compare against previously derived one-loop results; and rederive the constraints for different unitarity conditions. We conclude in sec. VI.

II. MODEL
The scalar potential of a CP conserving THDM with softly broken Z 2 symmetry reads After EWSB, the neutral components of the two Higgs states receive vacuum expectation values (VEVs) of with v 2 1 + v 2 2 = v 246 GeV and tan β = v2 v1 . The mass spectrum consists of superposition of these gauge eigenstates, i.e.
Here, G and G + are the Goldstone modes of the Z and W boson. The mixing in these sectors is fixed by tan β, while in the CP-even sector arXiv:1805.07310v1 [hep-ph] 18 May 2018 a rotation angle α defines the transition from gauge to mass eigenstates. In practical applications, one can trade the physical masses m h , m H , m A and m H + as well as tan β and tan α for the quartic couplings. The necessary relations are with t β = tan β and t α = tan α. This has the advantage that physical observables instead of Lagrangian parameters can be chosen as input. However, one needs to be careful since a randomly chosen set of masses could easily correspond to a problematic set of quartic couplings: for very large couplings perturbativity will be spoilt and also unitarity can be violated.

III. UNITARITY CONSTRAINTS
Perturbative unitarity constraints come from applying the unitarity of the S-matrix for 2 → 2 scalar field scattering amplitudes. We calculate a matrix a ba 0 given by which is derived proporional to the zeroth partial wave of scattering pairs of scalars a to pairs b having matrix element M(cos θ), where θ is the angle between the incoming and outgoing three-momenta (p a , p b respectively) in the centre-of-mass frame. The factor δ 12 (δ 34 ) is 1 if particles {1, 2}({3, 4}) are identical, and zero otherwise. We then find the eigenvalues of this matrix, which we denote a i 0 , and insist that they must satisfy Classic unitarity constraints for the THDM have been calculated in the limit of large scattering energies, in which case only the quartic couplings contribute to scattering and the momentum dependence of the prefactor of the integrand in (8) disappears; moreover all diagrams with propagators are suppressed by the collision energy squared and can be neglected, so the final result appears superficially independent of the scattering energy. This has been applied at tree [] and one-loop [] level. The limits on the quartic couplings at tree level in this approximation are Max |λ 3 ± λ 4 | , λ 1 + λ 2 ± (λ 1 − λ 2 ) 2 + λ 2 4 , |λ 3 ± λ 5 | , However, it has not been tested if the large s approximation is valid in all BSM models in which it is applied. It could be that large contributions are present at smaller s which then rule out given parameter regions in the considered model. The theory could develop a Landau pole before s is sufficiently large to neglect the masses, or could be defined with a low cutoff. And at large values of the couplings, their running is usually sufficiently fast so that the values of the couplings at an energy scale √ s are vastly different from those at lower energies. So in order to be able to test unitarity at finite s, the Mathematica package SARAH has now been extended. The salient features are: (i) all tree-level diagrams with internal and external scalars are included to calculate the full scattering matrix; (ii) We neglect all gauge couplings, and treat Goldstone bosons as physical particles with mass equal to the gauge boson; (iii) the calculation is done in terms of mass eigenstates, i.e. the full VEV-dependence is kept; (iv) the numerical evaluation is done with the Fortran code SPheno [23,24]; (v) large enhancements close to poles are cut in order not to overestimate the limits. This is demonstrated at oneexample in Fig. 1. More details and derivations of our full procedure are given in the accompanying paper [25].

IV. RESULTS
In this section we shall study the impact of the improved unitarity constraints on the two Higgs doublet model at tree level. We have chosen for our discussion type-I, but the results hold also for other models, becuase our we omit fermions from our scattering processes. Hence there is only an indirect difference between the constraints for type-I and type-II: the limits from flavour observables are stronger for light charged Higgs masses for type-II. Hence, the m H + must be larger in general for type-II [26]. On the other hand, we include the constraints from Higgs searches via HiggsBounds [27][28][29], which can vary to a lesser extent between type I and II models.
Our numerical analysis is based on the SPheno [23, 24] interface of SARAH [30][31][32][33][34]. By default, SPheno calculates the mass spectrum at the full one-loop level and includes all important two-loop corrections to the neutral scalar masses [35][36][37]. However, we shall not make use of these routines in the following but work at tree-level, or equivalently under the assumption that an OS calculation works in principle (with all the caveats discussed in Ref. [38]). This is because we cannot (yet) calculate quantum corrections to unitarity at finite s, and when the couplings are large in almost all cases the quantum corrections to masses/couplings become very large: this gives further motivation for including only constraints at finite s!

A. The light Higgs window
We start with a discussion of the effects in the case that both CP even Higgs states have masses of 125 GeV or below. A comparison between the 'classic' -equation (10) -and new constraints is given in Fig. 2. Obviously, one finds much stronger constraints in two different cases once finite s is considered: (i) for smallish |M 12 | the wedge M A = M H + m H disappears; (ii) for larger |M 12 | the scattering amplitude in the overall (m A , m H + grows signifcantly. The responsible channels and best scattering energies causing these effects are quite different: • Small |M 12 |: consider the following simplified hierarchy, together with tan β = −1/(tan α) = 1. The dominant channels are those with heavy external states and a light Higgs exchange. For instance, the am-plitude AA → AA can be approximated as From that, we get that the ratio compared to the old constraints This ratio becomes maximal slightly above the kinematic threshold s Threshold = 4m 2 A and an enhancement of 2-3 is possible. Thus, the best scattering energy √ s is around 1-2 TeV.
• Large |M 12 |: in this case we can consider the following, simplified hierarchy Now, the dominant scattering processes are those with light external scalars only. The maximal eigenvalue of the full scattering matrix is roughly given by diagonalising the submatrix with CP-even states only By doing that, we find that the ratio between the old and new results scales as Thus, this ratio grows very quickly with increasing M 12 and one finds very strong unitarity constraints already at scattering energies √ s of a few hundred GeV.
B. Heavier Scalar a. Stronger Constraints We turn now to the case that all new scalars are heavier than the SM-like Higgs. We start with a short analytical estimate for parameter regions in which difference between the our calculation and previous results show up. This is, for instance, the case for the configuration Assuming again the tan β = −1/ tan α = 1 for the moment, the maximal eigenvalue for the scattering matrix in the large s limit is We want to compare this with the scattering HA → HA scattering process which includes diagrams with the SMlike Higgs in the propagator. We find Thus, for s = 5m 2 A close to the kinematic threshold we find an enhancement of roughly | 2 √ 5 log m h m A | compared to the large s approximation. For m A = 700 GeV this corresponds to nearly a factor of 2. We can confirm this by making use of the full numerical machinery. In Fig. 3 we show the impact on the maximal allowed value for m H in the (m A , m H + ) plane while scanning over all other parameters as indicated in the caption. We see that this value shrinks significantly and a large region of the plane which is allowed by the old constraints is no longer accessible.
b. Weaker Constraints If we consider the scattering up to a finite √ s, we can find that the scatter eigenvalues become smaller compared to the limit √ s → ∞ for several reasons: (i) the dominant channels can be kinematically forbidden; (ii) there can be a negative interference between the point interactions and the propagator diagrams; (iii) the dominant channels can be cut out because of possible resonances in order not to overestimate the unitarity constraints. Due to these effects, one needs to ask the question to which energy scale we have actually probed scattering processes of scalars at the LHC. Of course, the LHC is running with √ s = 14 TeV. However, it is unrealistic to assume that the full energy is available in the 2 → 2 scattering of scalars. Moreover, there are different options to handle the t-and u-channel poles, which can appear if internal states become onshell, depending on how aggressive or conservative the limits should be: if we remove these poles either completely or only by a partial diagonalisation of the scattering matrix, large contributions to the scattering can be dropped at small s. We demonstrate via one example in Fig. 4 where the maximal eigenvalue as a function of √ s is shown. If we completely ignore the t-and u poles we see a huge enhancement close to some kinematic thresholds. In contrast, if we work with a partial diagonalisation as proposed in Ref. [22] we see that we find the eigenvalue of the large s approximation only for √ s > 10 TeV. This might be rather surprising since all involved masses are below 1 TeV!

V. COMPARISON WITH LOOP CORRECTIONS
Since one of our motivations for considering finite s scattering is that the quantum corrections to masses and couplings become large as we increase the scattering energy, it is also important to examine the effect of loop corrections to unitarity. Moreover, the boundary of unitarity may also coincide with a loss of perturbativity. In general, loop corrections to unitarity have been very little studied in BSM models; however, in the context of the THDM, they were considered in Ref. [39] in the limit of √ s much larger than the masses in the theory. We can therefore make a direct comparison. In that paper, they presented general formulae for the loop corrections to a 0 in terms of the quartic couplings of the theory evaluated at the scale √ s, which are effectively independent of the particle masses. Results for two scenarios, one with an SO(3) symmetry and another with "MSSM-like" couplings, were presented.
We shall make our comparison with the "MSSM-like couplings"; in SARAH conventions this means With these restrictions the 'classic' tree-level constraints of equation (10) simplify to which describe a rhombus inlcuding the origin. Requiring stability of the potential requires which, when we combine the two, leaves a portion of the parameter space where λ 1 is at most 4 3 π, and λ 3 < 4π. In the previous sections, we applied the unitarity constraint |Re(a i 0 )| < 1/2, but in [39] they apply a different constraint, which we shall now examine. The starting point is the equation (for an elementary derivation see [25]). Naively this gives simply |a i 0 | ≤ 1, which is a constraint sometimes appllied, but with a little rearranging we have which gives the classic limit (9). This limit makes no assumption of perturbativity, and indeed when Re(a i 0 ) obtains its maximum value then Im(a i 0 ) = |Re(a i 0 )|. Since Im(a i 0 ) is only generated at first at one loop order, then saturating this bound would potentially require violating perturbativity. On the other hand, rearranging again, we can write the above as If we have complete ignorance of Im(a i 0 ) then we just recover the same constraint as above. However, if we have calculated a 0 at one loop and assume that perturbativity holds, then we can use our calculated values for the real and imaginary parts of a i 0 and use the above constraint. Focussing on one eigenvalue, let us write and expand eq. (26) then we find Now Ref. [39] then appeal to perturbation theory so that b I = (a 0 0 ) 2 + higher order terms (29) and then obtain This can then be a very strong constraint. In Fig. 5 we show the constraints from applying eq. (26) as done by Ref. [39], with the constraints from our trilinear couplings and the tree-level constraints for comparison. The tree-level quartic-only and one-loop constraints are independent of tan β and all of the mass scales (except that they should be interpreted as couplings evaluated at a renormalisation scale √ s), whereas for our scan we choose two values of tan β (marked on the plot) and fix the tree-level lightest Higgs mass to be 125 GeV -this is enough to determine all of the remaining free parameters once λ 1 and λ 3 are specified.
We see from Fig. 5 that even though the loop-level constraints seem extremely severe, our tree-level trilinear constraint still removes a significant chunk of the remaining parameter space.
However, these one-loop constraints have the curious feature of excluding couplings near the origin, which arises from the regions where one scattering eigenvalue vanishes at tree level. Indeed, from eq. (30) we see that if a 0 0 = 0 (as can happen for linear combinations of the couplings) then unitarity is apparently violated. In the FIG. 5. Tree-level and one-loop constraints on λ1 and λ3 in the "MSSM-like" THDM. Quartic-only tree-level constraints are shown as black dot-dashed lines, the vacuum stability constraint λ3 ≥ −2λ1 is the red dashed line; our tree-level contraints including trilinears are labelled with tan β = 2 and tan β = 30. The one-loop allowed region from [39] is the white region enclosed by the solid purple and orange curves. The second plot is a zoom into the first one.
notation of Ref. [39] the purple curve corresponds to the eigenvalue a 110odd 0 , which derives from the scattering of where τ 3 = 1 2 1 0 0 −1 and gives the scattering eigenvalue at tree level of a 0 0 = 2λ 1 + 2λ 3 . The orange curve corresponds to a 000even 0− from scattering which give the scattering eigenvalues at tree level of We therefore see that the one-loop constraints arise starting from the lines λ 1 + λ 3 = 0 and 4λ 1 + λ 3 = 0.
The reason for this is, however, assuming that the higherorder terms in eq. (29) are not important. Indeed, in the cases where a 0 0 = 0 for λ i = 0 we would apparently badly violate perturbation theory -but this is just because we have only computed up to one loop, and have a tuned cancellation at tree level. Since eq. (30) compares a tree-level and one-loop amplitude this seems particularly bad. Hence, if we examine the perturbation series more closely, specialising to the case of only quartic couplings for simplicity, and define λ to be a number of O(λ i ) as a perturbation series parameter, so that we see that a 2→n 0 is nonzero first for 2 → 4 processes at order λ 2 . Hence defining n>2 |a 2→n | 2 ≡ |X| 2 λ 4 we have, order by order in perturbation theory up to λ 4 : We see that the origin of eq. (30) depends on neglecting b 2,I , but if we include the information from eq. (34) then we would have obtained instead of eq. (30): b 2 R + higher order terms of indeterminate sign = 0.
Furthermore, when a 0 0 = 0 we simply recover b 3,I ≥ 0 and b 2,I = 0, which we could surmise from a 0 being of O(λ 2 ) and the standard unitarity relation. We do not obtain any new constraint beyond |Re(a 0 )| ≤ 1 2 . Hence in Fig. 6 we recompute the constraints at oneloop applying instead |Re(a 0 )| ≤ 1 2 for Re(a 0 ) = a 0 0 + b 1,R for the same scattering processes listed above. We use the expressions in the appendix from Ref. [39] to obtain the one-loop scattering amplitudes neglecting the wavefunction renormalisation contributions. These are mass-dependent and were found to be small in Ref. [39]. The reason is that, in the limit that √ s is much larger than all masses, only diagonal self-energies appear in the results which consist of expressions of the form Due to the presence of the trilinear couplings in these terms they appear at the same order as box and triangle diagrams. The one-loop constraint is then stronger than the "naive" tree-level one in some cases, and weaker in others; but we find that our tree-level constraints including the effect of trilinears are stronger than both in almost all cases.
For comparison, one could also check the one-loop allowed region for the sometimes used criterion |a 0 | ≤ 1. These are almost universally weaker than the tree-level constraints, implying that they are not sufficiently conservative, as can be expected.
FIG. 6. Tree-level and one-loop constraints on λ1 and λ3 in the "MSSM-like" THDM. Quartic-only tree-level constraints are shown as black dot-dashed lines, the vacuum stability constraint λ3 ≥ −2λ1 is the red dashed line; our tree-level contraints including trilinears are labelled with tan β = 2 and tan β = 30. The one-loop allowed region applying the constraint |Re(a0)| ≤ 1 2 is the white region enclosed by the solid orange curve.

VI. CONCLUSION
We have revised the tree-level perturbative constraints in THDMs by including the contributions from (effective) trilinear couplings, and provide an extension of the package SARAH which makes it possible to include these constraints in phenomenological studies in THDMs and many other BSM models. We found that the obtained limits can be significantly stronger than the ones usually applied in literature which are only correct in the limit of large scattering energies √ s. The importance of the improved constraints has been demonstrated by two chosen examples: (i) it was shown that the values of M 12 are highly constrained in the light Higgs windows; (ii) one finds a stronger upper limit for the CP-even Higgs mass in scenarios with |M 12 | < m A , m H + . On the other side, we have also discussed that the restriction to maximal scattering energies of a few TeV can revive points which violated unitarity only at much higher energies. We also made comparison with previous constraints derived at one-loop level in the large s approximation. Our results indicate that the tree-level constraints including trilinear couplings are the most important for this class of models, and are not superseded by the one-loop largemomentum constraints; instead, it would be a very interesting if rather complicated task to include the effects of the trilinear couplings at one-loop order, which could potentially strengthen constraints on these models further. In other BSM models similar -or even larger -differences between the full calculation and the large s approximation can be seen. This is discussed for example in Ref. [40] for several triplet extensions.