Next-to-leading order unitarity fits in Two-Higgs-Doublet models with soft $\mathbb{Z}_2$ breaking

We fit the next-to-leading order unitarity conditions to the Two-Higgs-Doublet model with a softly broken $\mathbb{Z}_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 $\mathbb{Z}_2$ symmetry with a probability of $95\%$. All fits are performed using the open-source code HEPfit.


I. 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 symmetry breaking, e.g. the Higgs sector of the minimal supersymmetric Standard Model (MSSM) [7-9], twin Higgs models [10][11][12][13][14], composite Higgs models [15,16]. Also, there is no fundamental reason for the minimality of the SM scalar sector, and multiple generations are known to exist in the fermion sector. Furthermore, the uncertainties in the SM Higgs coupling measurement [17] do not exclude the presence of additional scalars.
After the SM, the simplest and most straightforward extension of the SM is the addition of another Higgs doublet, the so-called Two-Higgs-Doublet model (2HDM) [18][19][20], which has been analysed in great detail in the literature, see for instance . It is interesting to study the unitarity bounds in the 2HDM because the scale at which new particles are expected to appear based on naturalness arguments is the same scale as the Lee, Quigg, Thacker upper limit [45,46] on the Higgs mass in the SM, which is of order of 1 TeV. In fact, there exists a large number of works studying the tree-level unitarity bounds on the quartic couplings, λ i , and Higgs masses of the 2HDM, see e.g. [47][48][49][50][51][52][53]. Unlike in the SM, extracting the bounds on the masses of the Higgs boson from the bounds on the quartic couplings is not straightforward because in the 2HDM the quartic couplings are in general functions of more parameters than just the masses of the Higgs bosons and their corresponding vacuum expectation values (VEVs). Recently, the perturbative unitarity bounds in the CP-conserving softly-broken Z 2 symmetric 2HDM were analyzed at the oneloop level [54]. This calculation settled a particular issue regarding how to estimate higherorder effects on available upper limits on the quartic couplings: In the SM, the unitarity bounds had been determined beyond the leading order (LO), and the typical result was that the bounds on the (RG-improved) quartic coupling of the SM were improved by a factor of a few with respect to the tree-level analysis [55][56][57][58]. While it was known that the treelevel unitarity bounds in NP models were likely to be overly conservative, there was no well defined way to implement stricter bounds using only tree-level results. Inspired by the SM results, it was advocated to re-scale the tree-level conditions by a factor of 1/4 to estimate higher order contributions [32]. A renormalization group analysis at next-to-leading order (NLO) confirmed this prescription if one wants a stable Higgs potential beyond 10 TeV [41].
However, now that an explicit NLO computation is at hand, this uncertainty on how to treat higher order corrections to the partial-wave amplitudes has been removed.
In this article, we improve on the results currently available in the literature in two main ways: regarding the unitary constraints, we go beyond the leading order precision by employing one-loop corrections which are enhanced, O(λ i λ j /16π 2 ), in the limit s |λ i |v 2 M 2 W , s m 2 12 to all the 2 → 2 longitudinal vector boson and Higgs boson scattering amplitudes. Secondly, we perform global parameter fits including the most up-to-date Run I ATLAS and CMS results, rather than only using a handful of benchmark scenarios, which might not cover the whole spectrum of interesting features.
The structure of this article is as follows: We give a short introduction to the model and its constraints in sections II and III, respectively. The statistical framework is presented in section IV. Section V contains the results of the fits. We conclude in section VI. Supplementary figures can be found in Appendix A, while we list the formulae for the NLO unitarity criteria and the fit inputs in Appendix B and Appendix C, respectively.

II. 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 flavourchanging 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 tree-level; the literature refers to this as the alignment limit [19,[59][60][61].
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

III. 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)

A. Theoretical constraints
On the theory side, constraints on the 2HDM come from the following requirements: • the Higgs potential must be bounded from below [62] between M Z and 750 GeV, • the minimum of the Higgs potential at 246 GeV should be the global minimum [63], • 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 [54].
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) [41].
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 [54].
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 an inequality, The 2HDM one-loop corrections necessary to use this inequality were recently computed in Ref. [54]. 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 [47][48][49][50][51][52][53]. Comparing with the discussion of higher order corrections in the SM, stronger bounds were estimated and used for the 2HDM [32,41], 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. [54] 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. [54] 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 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 [54], which can also be found in Appendix B. In order to satisfy unitarity, the a Z 2 Y σ± have to individually fulfill the condition (5). Note that at LO, the eigenvalues are related to the ones defined in [51] Another constraint is the requirement that higher order corrections to the partial wave amplitudes are suppressed. In particular, following [54] 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. [57]. 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 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 (5) 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.

B. Experimental constraints
The experimental constraints included in our analysis are: 1 Ref. [54] used the differential operator D GMU = 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µ).
• the Peskin-Takeuchi parameters S, T , and U [64], • the h signal strengths, • the non-observation of H and A at the LHC and • the B s meson mass difference ∆m Bs [65,66] and the branching ratio B(B → X s γ) [67].
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 [64]. 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 [68], see Table II 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 [69]; the SM branching ratios were calculated with HDecay 6.10 [70]. In order to express the 2HDM cross sections and branching ratios in terms of the SM ones, we make use of the formulae of [71] for the loop induced decays of the neutral Higgs bosons. Central values, errors (Gaussian approximation) and correlations for the signal strengths in (11) and (12) were obtained from Fig. 13 and Table 14 of [17] and can be found in Table III 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 .
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 IV in Appendix C. Most SM Higgs production cross sections are taken from the LHC Higgs Cross Section Working Group [72]; the remaining ones are calculated with HIGLU 4.34 [73], Sushi 1.5 [74], and Madgraph5 2.2.2 [75]. The branching ratios were calculated with HDecay 6.10 [70] while the decay widths for both Higgs-to-Higgs decays and Higgs decays into a Higgs boson and a gauge boson are taken from [76].
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 [65,66] at LO. For the inclusive measurement of B(B → X s γ), NNLO corrections are important [67]. 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 V in Appendix C.

IV. HEPFIT
As numerical set-up we use the open-source code HEPfit [77], interfaced with the release version of the Bayesian Analysis Toolkit (BAT) [78]. 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), and the posterior distribution of the parameters in the Higgs potential (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.

V. 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 Sec. III A, 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 Fig. 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 (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 Sec. III A is shown in the left panel of Fig. 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. [41], 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 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 Fig. 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 Fig. 3. Like in Fig. 1, we observe a hierarchy between LO unitarity, NLO unitarity and NLO unitarity with the perturbativity requirement R 1 < 1, with the first set of constraints being the weakest bound and the last being the strongest bound. While LO unitarity allows for maximal |m H − m A |, |m H − m H + | and |m A − m H + | of 500 GeV, the perturbative NLO unitarity conditions sets upper limits on the mass splittings of around 360 GeV. m H + > m H and m A > m H are already almost excluded by LO unitarity for λ 3 < 0 and λ 5 > 0, respectively; we can see that after the inclusion of NLO unitarity with the R 1 condition also other constellations like λ 3 = 0 feature significantly smaller possible mass differences. Fig. 4 contains the same quartic coupling planes as Fig. 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 III A and is identical with the blue contours of the previous figures. The unfilled contours have been obtained using only one of the following  Fig. 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 Z 2 ,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. 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"). In both types, the first set is most constraining for negative λ 3 or 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 Fig. 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 Fig. 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 0.15 and 0.04 in type I and II, respectively. (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 Fig.s 8 and 9 in Appendix A.
In Fig. 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%), 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 re- maining 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 Fig.s 10 and 11 in Appendix VI. 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 Fig. 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.

VI. 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. [54], 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 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.
As a next step, we have added all the 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 0 ≤ λ 1 < 4.2, 0 ≤ λ 2 < 1.2, −1.3 < λ 3 < 3.1, −2.5 < λ 4 < 2.9, −2.7 < λ 5 < 0.5 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  Fig. 4 with the difference that also the grey type I contour was filled here.
latter case, also an unbroken Z 2 symmetry can be ruled out at 95 Fig.s 8 to 11, and the quartic couplings of the so-called Higgs basis in Fig.s 12 and 13.
In Fig. 8, we show the effect of the h signal strengths on the tan β vs. β − α plane for 2HDM of type I. In the top left panel of Fig. 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 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 Fig. 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. Fig. 9 displays the same panels as Fig. 8 but for type II.  [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.  interplay: also in the m H/A vs. tan β planes, small tan β values are excluded in the fit with all constraints, because the masses of the neutral Higgs bosons cannot be very different from the H + mass due to unitarity (see Fig. 6).
Instead of the general parametrisation of the potential in (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 [79,80], and its potential can be written as Only five of the seven quartic couplings Z i are linearly independent. One can see from Fig. 12 that they get constrained by the different unitarity conditions in a similar way than the λ i in Fig. 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 R 1 ) sets upper limits of approximately 8 (5) on the absolute values of the Z i . Analogous to Fig. 4, we also show the impact of the experimental constraints on the Z i vs. Z j planes in Fig. 13. Especially Z 1 , Z 6 and Z 7 suffer strong additional restrictions from the experiments.

APPENDIX B
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. [54]. The complete expressions for a 0 can be found in Appendices B and C of Ref. [54].
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.  Fig. 13 and Table 14 of [17].

APPENDIX C
In the Tables II to V we list all used experimental inputs for our fits with their corresponding references.