Next-to-leading order unitarity fits in Two-Higgs-Doublet models with soft ℤ2 breaking

We fit the next-to-leading order unitarity conditions to the Two-Higgs-Doublet model with a softly broken ℤ2 symmetry. In doing so, we alleviate the existing uncertainty on how to treat higher order corrections to quartic couplings of its Higgs potential. A simplified approach to implementing the next-to-leading order unitarity conditions is presented. These new bounds are then combined with all other relevant constraints, including the complete set of LHC Run I data. The upper 95% bounds we find are 4.2 on the absolute values of the quartic couplings, and 235 GeV (100 GeV) for the mass degeneracies between the heavy Higgs particles in the type I (type II) scenario. In type II, we exclude an unbroken ℤ2 symmetry with a probability of 95%. All fits are performed using the open-source code HEPfit.


Introduction
Run I of the Large Hadron Collider (LHC) has concluded with the discovery of the last missing piece of the Standard Model (SM) -the Higgs boson [1,2]. It has tested the validity of the SM in a previously unexplored regime of energy, and has not found any significant deviations from the SM, hinting at a gap in the mass spectrum from the SM to whatever lies beyond it. This picture is consistent with the absence of indications of New Physics (NP) coming from indirect searches (e.g. electroweak precision or flavour observables). Certainly Run I of the LHC did not address all the shortcomings of the SM -among which are the hierarchy problem, dark matter and an explanation for the flavour pattern. In order to find a solution for these problems the presence of NP is inevitable. One of the key questions for Run II of the LHC is then at what scale the NP appears. The measured value of the Higgs mass, m h ≈ 125 GeV [3], hints at an answer to this: it tells us that the Higgs potential of the SM is not stable up to very high energy scales [4,5]. Thus one has to introduce an additional mechanism or new degrees of freedom to stabilize the Higgs potential if one wants to exclude the possibility of vacuum metastability. Furthermore, to keep the SM Higgs mass natural, new degrees of freedom around the TeV scale are required [6].
One well motivated direction for discovering physics beyond the SM is to search for additional Higgs bosons. These particles often arise in natural theories of electroweak JHEP11(2016)026

Model
The two-Higgs-doublet model (2HDM) [18][19][20] is a simple and straightforward extension of the Standard Model (SM), obtained through the addition of a second Higgs doublet to the SM field content. A characteristic of general 2HDMs is the existence of flavour-changing neutral currents (FCNC) mediated by tree-level exchange of neutral Higgs bosons. A natural way to eliminate these potentially dangerous FCNC is to require that the Yukawa interactions respect a discrete Z 2 symmetry, which can be broken softly. The Z 2 symmetry can be chosen in four independent ways, depending on the Z 2 charge assignments for quarks and charged leptons; this lead to four different types of 2HDM which are referred to as type I, type II, type X (lepton specific) and type Y (flipped). The type II model is of particular interest because the Higgs sector of the MSSM is a 2HDM of type II. The models we focus on in this paper are the CP-conserving 2HDM of type I and II with a softly broken Z 2 symmetry. The most general Higgs potential in a CP-conserving 2HDM with a softly broken Z 2 symmetry reads where Φ 1 and Φ 2 are the two complex SU(2) L Higgs doublets with hypercharge Y = 1/2 and the eight scalar potential parameters are real to avoid explicit CP-violation, with m 2 12 being the Z 2 soft-breaking parameter. At the global minimum of the scalar potential V the neutral components of Φ 1 and Φ 2 acquire VEVs, v 1 / √ 2 and v 2 / √ 2, respectively, which are fixed by the minimization of the scalar potential and must satisfy v 2 1 +v 2 2 ≡ v 2 ≈ (246 GeV) 2 . The ratio of the two VEVs is defined as tan β = v 2 /v 1 , where 0 ≤ β ≤ π/2. Assuming no CP-violation in the Higgs sector, the physical scalar spectrum consists of two CP-even states h and H with m h < m H , the CP-odd state A and the charged state H ± . The masses of these scalar bosons are denoted as m φ with φ ∈ {h, H, A, H ± }. Throughout this paper we interpret the observed Higgs resonance as the light CP-even scalar h and thus treat m h as fixed by measurements, m h = 125.09 GeV [3]. We choose the independent physical parameters of the model to be where α is the mixing angle of the neutral CP-even 2HDM Higgs bosons. In this parametrization, the tree-level couplings of the Higgs bosons to vector bosons and fermions only depend on tan β and β − α. Moreover, for β − α = π/2 the couplings of h to SM fermions and vector bosons are SM-like and H does not couple to vector bosons at treelevel; the literature refers to this as the alignment limit [19,[62][63][64].
Considering only the third generation of fermions, the Yukawa Lagrangian under the above-mentioned Z 2 symmetry takes the following shape: where the top quark couples to Φ 2 by convention and the index k is 2 in type I and 1 in type II. The top Yukawa coupling is related to the SM value Y SM

Constraints
In this section we list and discuss the theoretical and experimental constraints we impose on the 2HDM parameter space. Since we want to combine them in a Bayesian fit, we list the priors on the parameters from (2.2):

Theoretical constraints
On the theory side, constraints on the 2HDM come from the following requirements: • the Higgs potential must be bounded from below [65] between M Z and 750 GeV, • the minimum of the Higgs potential at 246 GeV should be the global minimum [66], • the 2HDM quartic couplings λ i (i = 1, 2, 3, 4, 5) and the Yukawa couplings are assumed to be perturbative (i.e. smaller than 4π and √ 4π in magnitude, respectively) at least up to 750 GeV, • the S-matrix of 2 → 2 scattering processes for Higgs bosons and longitudinal vector bosons should be unitary up to NLO, and its NLO eigenvalues should not exceed the LO eigenvalues in magnitude [57].
Requiring positivity and perturbativity of the couplings to hold at least up to 750 GeV is motivated by the fact that this scale is well above the electroweak symmetry breaking scale and we can safely use the NLO unitarity conditions. For the renormalization group running we use NLO renormalization group equations (RGE) [44]. The first three bullet points have already been used in the literature and will be referred to as "stability up to 750 GeV" in the following. The fourth set of constraints has never been applied in a general 2HDM fit, which is why we want to explain the details of our approach in the following, referring to [57].
The unitarity of the S-matrix leads to constraints on the partial wave amplitudes of a theory, where a 2→n are the eigenvalues of the matrix of the -th, 2 → n partial wave amplitudes, a . Considering only 2 → 2 scattering (and dropping the superscript) this constraint becomes JHEP11(2016)026 an inequality, The 2HDM one-loop corrections necessary to use this inequality were recently computed in ref. [57]. Prior to this computation, the inequalities |Re(a 0 )| ≤ 1 2 or |a 0 | ≤ 1 were used to constrain the tree-level partial wave amplitudes [50][51][52][53][54][55][56]. Comparing with the discussion of higher order corrections in the SM, stronger bounds were estimated and used for the 2HDM [34,44], but this ansatz was controversial. Having at hand the calculated NLO unitarity conditions, we can determine the upper bound on the quartic couplings without any ambiguity of method.
The computation of ref. [57] was performed in the high energy limit, s |λ i |v 2 M 2 W , s m 2 12 , where the SU(2) L × U(1) Y symmetry is manifest. In this limit, a 0 is block diagonal at leading order, with blocks of definite weak isospin (σ) and hypercharge (Y ) (a >0 = 0 at leading order in this limit). Furthermore, the Z 2 -even and -odd states do not mix at tree-level, leading to smaller blocks. Due to the manifest symmetry at high energies, the calculation can be simplified by computing the amplitudes in the Z 2 basis using the non-physical Higgs fields, w ± j , n ( * ) The elements of a 0 are given by where, for example, In general, the block diagonal structure of a 0 does not hold beyond tree-level. However, it turns out that in the high energy limit, this structure is only broken by diagrams that correct the wavefunctions of the external legs, not by 1PI diagrams. Ref. [57] showed that the external wavefunction corrections are numerically subdominant with respect to the 1PI diagrams in some special cases. We confirm this and find it to be generalizable for all 2HDM scenarios with a softly broken Z 2 symmetry. Due to this relative numerical unimportance, we neglect the external wavefunction corrections throughout this work. In

JHEP11(2016)026
this approximation, the one-loop eigenvalues take the following form, 1 with the eigenvalues labeled as follows, a Z 2 Y σ± , and dropping the index = 0. B N is the block-diagonal element, (a 0 ) i,f , from eq. (B.N) in [57], which can also be found in appendix B. In order to satisfy unitarity, the a Z 2 Y σ± have to individually fulfill the condition (3.2). Note that at LO, the eigenvalues are related to the ones defined in [54] by a Z 2 Y σ± = −32π 2 Λ Z 2 Y σ∓ for λ 5 > 0. Another constraint is the requirement that higher order corrections to the partial wave amplitudes are suppressed. In particular, following [57] we define, where the (N)LO label denotes the pure (N)LO contribution. Similar criteria were used in the perturbative unitarity analysis of the SM in ref. [60]. Assuming that the power series is perturbatively stable, we want to require the NLO contribution to be smaller than the LO expression, hence R 1 < 1. However, we need to avoid the exclusion of accidentally small leading-order contributions. (For instance, a odd,LO 10 = (λ 4 − λ 3 )/(8π) is small if λ 3 ≈ λ 4 , while a odd,NLO 10 also depends on the other quartic couplings.) Therefore, we decided to use the R 1 criterion only if |a Z 2 ,LO Y σ± | > 0.02 ≈ 1/(16π). The 2HDM is unitary, so let us explain what we mean when we say unitarity constraints. Inequality (3.2) requires the couplings of a theory to be smaller than a certain value in magnitude, or else the theory will no longer appear to be unitary at the finite order of the perturbative expansion to which we are working. In this sense both the "perturbativity bound," R 1 , and the "unitarity bound" test the same thing, namely where the breakdown of perturbation theory occurs.

Experimental constraints
The experimental constraints included in our analysis are: • the Peskin-Takeuchi parameters S, T , and U [67], 1 Ref. [57] used the differential operator DGMU = 16π 2 µ 2 (d/dµ 2 ) in its definition of the beta functions.
In this work we use the traditional definition of the beta function, β λ i ≡ Dλi = µ(dλi/dµ).

JHEP11(2016)026
• the h signal strengths, • the non-observation of H and A at the LHC and • the B s meson mass difference ∆m Bs [68,69] and the branching ratio B(B → X s γ) [70].
As we saw in the previous section, the 2HDMs introduce new Higgs bosons which couple to the gauge bosons and which, thereby, can give contributions, through loops, to the gauge boson self-energies. Thus, the 2HDMs yield new contributions to S, T and U that generally move them away from their SM values. For the 2HDM predictions of the Peskin-Takeuchi parameters in the CP-conserving limit we make use of the formulae of [67]. As input values for the oblique parameters S, T and U and their correlation coefficients we take the most recent results obtained in a fit to electroweak precision data with HEPfit [71], see table 2 in appendix C. In order to confront the 2HDM with the latest ATLAS and CMS Run I data on Higgs signal strengths, we compute in the narrow-width approximation for each final state f ∈ {γγ, ZZ, W W, bb, τ τ } the signal strengths 2 having grouped the Higgs production modes in just two effective modes, ggF + tth and VBF + Vh, where "ggF", "tth", "VBF" and "Vh" stand for "gluon fusion", "tt associated production", "vector boson fusion" and "Higgstrahlung", respectively. The SM Higgs boson production cross sections are taken at a centre-of-mass energy of 8 TeV from [72]; the SM branching ratios were calculated with HDecay 6.10 [73]. In order to express the 2HDM cross sections and branching ratios in terms of the SM ones, we make use of the formulae of [74] for the loop induced decays of the neutral Higgs bosons. Central values, errors (Gaussian approximation) and correlations for the signal strengths in (3.8) and (3.9) were obtained from figure 13 and table 14 of [17] and can be found in table 3 in appendix C. Direct H and A searches are taken into account as follows: given the X → H/A → Y process, we define the ratio where σB = σ(X → H/A) · B(H/A → Y ) and the subscripts denote the theoretical 2HDM value of σB and its observed and expected exclusion limit at 95% CL by the experiments. With this definition, we can assume the R (X→H/A→Y ) Gauss ratios to be Gaussian with a standard deviation of 1. Note that these quantities depend on m H/A ; furthermore we neglect the error on σB| 95%,exp .

JHEP11(2016)026
The H and A search exclusion limits included in our analysis and their mass ranges, along with the exclusion plots from which they were digitalized, are listed in table 4 in appendix C. Most SM Higgs production cross sections are taken from the LHC Higgs Cross section Working Group [75]; the remaining ones are calculated with HIGLU 4.34 [76], Sushi 1.5 [77], and Madgraph5 2.2.2 [78]. The branching ratios were calculated with HDecay 6.10 [73] while the decay widths for both Higgs-to-Higgs decays and Higgs decays into a Higgs boson and a gauge boson are taken from [79].
From the plethora of flavour observables we only use the most relevant two for our 2HDM discussion: the mass difference in the B s meson system, ∆m Bs and the branching fraction B(B → X s γ). The former is calculated according to [68,69] at LO. For the inclusive measurement of B(B → X s γ), NNLO corrections are important [70]. As for fixed SM parameters this observable only depends on the two 2HDM parameters tan β and m H ± , we store the B(B → X s γ) values for various inputs of these two parameters in tables, and interpolate them linearly in the fits. A theoretical error of 7% is applied, which corresponds to the uncertainty in the SM parameters. The experimental inputs for the flavour observables can be found in table 5 in appendix C.

HEPfit
As numerical set-up we use the open-source code HEPfit [80], interfaced with the release version of the Bayesian Analysis Toolkit (BAT) [81]. The former calculates all mentioned 2HDM observables and feeds them into the parallelized BAT, which applies the Bayesian fit with Markov chain Monte Carlo simulations. The complete global fit with all theoretical bounds runs for approximately 60 hours with 12 parallel chains generating 2 · 10 7 iterations each. Adding the experimental observables as described above slows down the same fit to roughly 90 hours.
A fundamental difference between the Bayesian and the frequentist approach is the treatment of fine-tuning: if one changes the parametrization of a model, flat priors on the former parameters usually do not translate into flat distributions of the new basis in a Bayesian fit. Some values for a new parameter might only be obtained by a very specific constellation of the old parameters, which in that sense would mean that they require a certain amount of fine-tuning. A frequentist fit is not sensitive to this bias, but one could argue that it is also less natural. HEPfit makes use of the Bayesian approach assuming flat priors for the physical parameters (2.2), and the posterior distribution of the parameters in the Higgs potential (2.1) are "distorted" by the Jacobian of the change of parameters. However, the posterior intervals only have a well-defined meaning once experimental data is included, and that is when the dependence on the priors disappears. In the first part of the following section, when we only discuss theoretical constraints, the reader should bear in mind that our results depend on the priors (and thus on the parametrization). Also, we will present the 99.7% allowed regions for fits to only theoretical constraints, while after the inclusion of experimental data we show the 95.4% probability contours, which then have a statistical meaning.

JHEP11(2016)026 5 Results
In the following we will present the results of our fits of theoretical and experimental constraints to the 2HDM of type I and II. Before we address the physical 2HDM parameters we want to compare the effects of the unitarity constraints. As explained in section 3.1, we impose these bounds at a scale of 750 GeV; thus a stable Higgs potential at least up to this scale is implicitly assumed. Nonetheless, all quantities shown in the figures are to be understood at the scale M Z .
In figure 1 we show the 99.7% probability regions for all λ i vs. λ j planes with three different unitarity conditions: the green areas are allowed if we impose only LO unitarity, the red regions show the remaining parameter space if we use the NLO unitarity conditions, and the blue contours result from additionally requiring the NLO unitarity conditions to be perturbative (R 1 < 1). Generally one can see that perturbative NLO unitarity is a stronger constraint than NLO unitarity with arbitrary R 1 , which itself is always stronger than LO unitarity. Numerically the largest possible absolute values for any of the quartic couplings from (2.1) are 8.10, 7.21, or 5.75, if we apply LO, NLO or perturbative NLO unitarity, respectively. Especially in the λ 4 vs. λ 3 plane, but also in the λ 5 vs. λ 3 and λ 5 vs. λ 4 planes one can see that for particular constellations, NLO unitarity features sharp incisions towards the origin of the plane, whereas those indentations are absent if we use LO unitarity.
A closer look at different variations of our conditions explained in section 3.1 is shown in the left panel of figure 2 in the λ 4 vs. λ 3 plane. The green, red and blue solid lines correspond to the contours of the same colour in the previous figure; all lines are the 99.7% probability boundaries. As explained above, previous studies used 1/4 rather than 1/2 as upper limit for the real part of the LO unitarity eigenvalues. This choice is represented by the green dashed line and is almost always less stringent than R 1 -perturbative NLO unitarity. The red dashed line uses LO RGE instead of the NLO RGE which apply in all other cases. As already stated in ref. [44], the NLO RGE "stabilize" the potential with respect to the LO expressions in the sense that for the same starting point one runs into non-perturbative values for the quartic couplings at much lower scales with LO RGE. That is why larger values for the λ i are accessible at the electroweak scale if one uses the NLO RGE. What happens if one only requires that R 1 < 1 without imposing NLO unitarity can be seen at the pink contour. In other words, the blue line should be interpreted as the combination of the red (NLO unitarity) and the pink one. Only for λ 3 > 4, one can see that the combination of both sets of constraints is stronger than their individual impacts. Finally, the cyan contour is the result of using a ten times smaller threshold for the LO part of R 1 . However, compared to our typical limit of 0.02 this is not substantially different.
Again in the λ 4 vs. λ 3 plane, we also show the individual contributions of the relevant single NLO eigenvalues in the right panel of figure 2. The shaded areas are excluded at 99.7% probability by the eigenvalues indicated in the legend. It is worth noting that only "−" solutions seem to be important in this plane.
Going from the potential parametrisation to the physical parameters, one can see how the different constraints on the λ i couplings translate into restrictions on the mass differences between the heavy Higgs bosons H, A and H + in figure 3. Like in figure 1, we observe a hierarchy between LO unitarity, NLO unitarity and NLO unitarity with the JHEP11(2016)026   figure 1. The dashed green curve shows the effect of (arbitrarily) requiring the LO unitarity condition to be more restrictive. The pink curve demonstrates the impact of the perturbativity bounds without the unitarity bounds. The cyan curve requires |a Z2,LO Y σ± | > 0.002 in order for the perturbativity bounds to be enforced, rather than 0.02, leading to no significant change in the allowed parameter space. Right: breakdown of the single effects of the unitarity constraints. Only the most constraining eigenvalues are displayed. Figure 4 contains the same quartic coupling planes as figure 1, but additionally the experimental data has been taken into account. The blue region survives all theoretical constraints as mentioned at the beginning of section 3.1 and is identical with the blue contours of the previous figures. The unfilled contours have been obtained using only one of the following three sets of inputs in combination with the requirement that the scalar potential is stable up to 750 GeV: the oblique parameters (labelled "STU"); h signal strengths and H and A searches ("Higgs"); ∆m Bs and B(B → X s γ) ("Flavour"). For the latter two the coloured dashed lines represent type I fits, while the solid lines are the type II contours. In both types, the first set is most constraining for negative λ 3 or JHEP11(2016)026   positive λ 4 , the second set excludes λ 5 > 0.4, and the third set of inputs yields λ 2 < 3.5 and λ 3 > −2. Finally, the grey regions denote the combination of all theoretical and experimental constraints in type II and the grey dashed lines correspond to the type I fits. The allowed λ 1 and λ 4 intervals are similar to the ones obtained in figure 1, but the other three receive significant additional restrictions from the experimental bounds: with a probability of 95.4%, λ 2 cannot exceed 1.6 (1.2), λ 3 has to be within −1.6 and 3.0 (−1.3 and 3.1) and the allowed λ 5 interval is between −2.7 and 0.3 (−2.7 and 0.5), when marginalizing over all parameters in type I (type II).
Again turning towards the physical parameters, we show the allowed parameter space in the β − α vs. tan β plane for type I and II in the left and right panels of figure 5, respectively. The most important bounds in this plane comes from the h signal strengths, the heavy Higgs searches and ∆m Bs ; their 95.4% bounds are also depicted individually. The Higgs observables strongly constrain the difference between α and β. In the final fit with all constraints, the deviation from the alignment limit, |β − α − π/2|, cannot exceed JHEP11(2016)026 0.15 in type I and 0.04 in type II. (This corresponds to a maximal deviation of sin(β − α) from 1 by 0.01 and 7 · 10 −4 , respectively.) The mass difference in the B s system sets a type independent lower bound on β for the chosen mass priors. More details about the effects of the signal strengths can be found in figures 8 and 9 in appendix A.
In figure 6, we plot the allowed ranges for the heavy Higgs boson masses and their mass differences after imposing the theoretical and experimental constraints for type I and type II. The green, red and blue regions depict the 99.7% allowed parameter space for the various unitarity conditions discussed above. The orange region is the allowed by the ST U observables at 95.4%. Finally, the grey region is the available parameter space after all the theoretical and experimental constraints are taken into consideration. Even if the perturbative NLO unitarity contour represents a larger probability boundary than the oblique parameter contour (99.7% for the former, 95.4% for the latter), it is more stringent for masses above 400 GeV and thus the dominant constraint in the high mass regime. It allows for maximal mass splittings between m H , m A and m H + of around 100 GeV for masses above 600 GeV. After the inclusion of the LHC searches for heavy neutral Higgs bosons, we observe that the remaining parameter space is disconnected. The largest gap occurs around m H,A ≈ 550 GeV. The reason for this discontinuity is that our fits are incompatible with the observed ATLAS and CMS diphoton cross sections around this mass. For details, we refer to figures 10 and 11 in appendix A. With a probability of 95.4%, H and H + can be as light as 210 GeV, and m A cannot be smaller than approximately 400 GeV in type I. In type II, masses below 600 GeV are excluded at 95.4% after the inclusion of B(B → X s γ) to the fit.
Finally, we address the soft Z 2 symmetry breaking parameter m 2 12 . In figure 7 we show its dependence on tan β and the H Higgs mass in the two discussed types. While for the theoretical set of constraints a strong correlation between the heavy Higgs mass and m 2 12 is visible, this gets somewhat relaxed if one adds experimental data to the fit. This is due to the flavour constraints, which favour larger tan β and m H + values. The most important result here is that an unbroken Z 2 symmetry can be excluded with a probability of 95.4% in the combined fit to the type II; the single sets of constraints are individually compatible with an exact Z 2 symmetry. The lowest 95.4% allowed value for m 2 12 is (370 GeV) 2 , if we marginalize over all other parameters.

Conclusions
The determination of the NLO unitarity constraints to the 2HDM with a softly broken Z 2 symmetry mitigates the problem of how to tame higher order contributions involving large quartic couplings. The expressions have been derived in ref. [57], and in this article we perform the first general fits to them in the 2HDM of type I and II, making use of the publicly available package HEPfit. One important result is that wavefunction renormalization contributions can be safely neglected in these models. In our fits we also apply the suppression of non-perturbative higher order contributions with the R 1 condition, requiring that the NLO part cannot be larger in magnitude than the LO contribution if the latter is not accidentally small. We find that both steps, going from LO to NLO unitarity and comparing NLO unitarity with R 1 -perturbative NLO unitarity, individually put strong JHEP11(2016)026   figure 4 with the difference that also the grey type I contour was filled here.

JHEP11(2016)026
bounds on the 2HDM parameters. If we add all other relevant theoretical constraints, that is stability and positivity of the scalar potential up to a scale of 750 GeV, the quartic λ i couplings cannot exceed 5.8 in magnitude and the mass differences between m H , m A and m H + cannot be larger than approximately 360 GeV. (The latter even reduces to maximally 100 GeV for heavy Higgs masses above 800 GeV.) To our knowledge, this currently represents the strongest reliable bound on the mass splittings.
Moreover, we have added the most relevant experimental constraints to the fit: the electroweak precision data in form of the oblique parameters, the complete set of LHC Run I results and the most important flavour observables. These bounds constrain the quartic couplings even further: the allowed intervals for the quartic couplings are 0 ≤ λ 1 < 4.2, 0 ≤ λ 2 < 1.6, −1.6 < λ 3 < 3.0, −2.5 < λ 4 < 2.9, −2.7 < λ 5 < 0.3 in type I and in type II with a probability of 95.4%. For the physical parameters, we find that tan β cannot be smaller than 1 in both discussed types of the 2HDM. The deviation from the alignment limit |β − α − π/2| cannot exceed 0.15 (0.04) in type I (type II). In type I the global fit produces lower 95.4% bounds of 210 GeV for m H and m H + and 410 GeV for m A , while these limits are around 650 GeV for all three heavy Higgs masses in type II. In the latter case, also an unbroken Z 2 symmetry can be ruled out at 95.4%; the soft Z 2 breaking parameter m 2 12 has to be larger than (370 GeV) 2 .

Acknowledgments
We

A Additional figures
In this appendix, we present supplementary figures of the 2HDM parameter space: dedicated fits of the different signal strength measurements and H and A searches in type I and II in figures 8 to 11, and the quartic couplings of the so-called Higgs basis in figures 12 and 13.
In figure 8, we show the effect of the h signal strengths on the β − α vs. tan β plane for 2HDM of type I. In the top left panel of figure 8, the effect of considering all the five signal strengths in the "ggF+tth" production modes on this plane is represented by the orange shaded region, considering all the five signal strengths in the "VBF+VH" production JHEP11(2016)026  [17], which are also approximated by the black dotted ellipses. In the top left panel, we also mark the alignment limit β − α = π/2 by a grey dotted line.

JHEP11(2016)026
modes are shown in the pink region, and the light green shaded region depicts the allowed parameter space when all the ten signal strengths are taken into consideration. In order to compare the latter with the effect of each of the five decay modes individually, we separately plot the single decay modes at a time on the rest of the panels of figure 8 overlaid with the fit with all signal strengths. For each of these additional panels we also indicate the latest 8 TeV signal strength correlation contours at 68% CL taken from ref. [17]. In all the panels, the filled regions with solid (dashed) lines represent the 68.3% (95.4%) probability contours as obtained from the fits. Figure 9 displays the same panels as figure 8 but for type II. Figures 10 and 11 compare the most important constraints on the heavy Higgs masses vs. tan β planes (left column) and on the heavy Higgs masses vs. β −α planes (right column) in type I and type II, respectively. In the first four panels, the relevant H (top row) and A (second row) searches are represented by the shaded regions, which they exclude. For the attribution of the colours, we refer to the legends. The left panel of the bottom row of figures 10 and 11 shows in beige the constraint from the mass difference in the B s system on the charged Higgs mass vs. tan β plane, which excludes tan β < 1 for the chosen m H + range. Light m H and m A are excluded in type II (figure 11) because the masses of the neutral Higgs bosons cannot be very different from the H + mass due to unitarity (see figure 6).
Instead of the general parametrisation of the potential in (2.1), one is free to choose a basis in which only one of the two transformed doublets, H 1 and H 2 , obtain a VEV. This basis is called the Higgs basis [82,83], and its potential can be written as Only five of the seven quartic couplings Z i are linearly independent. One can see from figure 12 that they get constrained by the different unitarity conditions in a similar way than the λ i in figure 1, with the R 1 -perturbative NLO expressions being stronger than simple NLO unitarity, which itself is an improvement of LO unitarity. While the latter does not allow for |Z i | > 9, NLO unitarity (with   the experimental constraints on the Z i vs. Z j planes in figure 13. Especially Z 1 , Z 6 and Z 7 suffer strong additional restrictions from the experiments.

B Diagonal entries of the NLO S-matrix
For the reader's convenience we list the minimal set of elements of the matrix a 0 needed to write its eigenvalues, a 0 , at next-to-leading order accuracy in the limit that the wavefunction renormalization contribution is neglected. In what follows, each B N corresponds in this approximation to eq. (B.N) of appendix B of ref. [57]. The complete expressions for JHEP11(2016)026 a 0 can be found in appendices B and C of ref. [57].
It is worth mentioning that here only the LO expressions for the β functions should be used in order to be consistent with the order of perturbation theory. For the running in the fits we apply NLO RGE.

C Experimental inputs
In the tables 2 to 5 we list all used experimental inputs for our fits with their corresponding references.

JHEP11(2016)026
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.