Update of Global Two-Higgs-Doublet Model Fits

We perform global fits of Two-Higgs-Doublet models with a softly broken $\mathbb{Z}_2$ symmetry to recent results from the LHC detectors CMS and ATLAS, that is signal strengths and direct search limits obtained at $\sqrt{s} = 8$ TeV and $\sqrt{s}=13$ TeV. We combine all available ATLAS and CMS constraints with the other relevant theoretical and experimental bounds and present the latest limits on the model parameters. We obtain that deviations from the so-called alignment limit $\beta-\alpha=\pi/2$ cannot be larger than $0.03$ in type I and have to be smaller than $0.02$ in the remaining three types. For the latter, we also observe lower limits on the heavy Higgs masses in the global fit. The splittings between these masses cannot exceed $200$ GeV in the types I and X and $130$ GeV in the types II and Y. Finally, we find that the decay widths of the heavy Higgs particles cannot be larger than $7\%$ of their masses if they are lighter than $1.5$ TeV.


Introduction
The discovery of a new scalar resonance with a mass around 125 GeV [1,2] in the Run 1 phase of the Large Hadron Collider (LHC) has paved the way for new directions in high-energy particle physics. Analyzing the properties of this particle has suggested strong evidence that it is the Higgs boson of the Standard Model (SM), i.e. a scalar CP-even state which has SM-like couplings to the other particles. Currently the combined analysis based on the Run 1 (7 and 8 TeV) LHC data shows that its couplings with the vector bosons are found to be compatible with those expected from the SM within a ∼ 10% uncertainty, whereas the coupling to the third generation fermions (top, bottom quarks and the τ lepton) is compatible within an uncertainty of ∼ 15 − 20% [3]. Thus the current status of the Higgs properties still allows to explore new interpretations of the observation coming from new physics of different underlying structures.
The Two-Higgs-doublet model (2HDM) [4][5][6] is one of such extensions of the Standard Model. Alike other popular NP models, it gets more and more constrained by recent experimental progress, especially by the LHC data . As the name suggests the 2HDM has two Higgs doublets in contrast to the single Higgs doublet in the SM. This extension of the Higgs sector leads to the existence of five scalar bosons, namely a heavy and light CP-even Higgs boson, H and h, a CP-odd Higgs boson, A, and a pair of charged Higgs bosons, H ± . Whether the scalar boson observed in the Run 1 of LHC is a part of an extended Higgs sector is an outstanding question and is at the cynosure of attention of the current Run 2 (13 TeV) phase of the LHC.
The questions we ask is: Which parts of the 2HDM parameter space are favoured after imposing the latest experimental data from the LHC? Compared to the Run 1 phase of the LHC, in Run 2 the situation has changed in several respects: Alongside the latest Higgs signal strength data by the ATLAS and CMS collaborations also many of the recent results of searches for additional heavy Higgs bosons are more constraining than the Run 1 data. Both experiments have performed dedicated searches for new signatures in various possible final states at the LHC: besides the fermionic final states tt, bb, τ + τ − , tb and τ + ν they include gauge bosons (γγ, Zγ, ZZ, W + W − ) and Higgs particles (hh, hZ, HZ, AZ) as the decay products of a heavy resonance. So far, these searches for heavy resonances have remained elusive in the ATLAS and CMS data, and thus the measurements put modelindependent 95% C.L. upper limits on the production cross section times branching ratios for different production processes and decay modes. In the present work, we assess the status of all four types of softly broken Z 2 symmetric 2HDM when all the experimental constraints coming from the latest LHC data are taken into account. We confront these with the theoretical constraints on these models (positivity, stability and next-to-leading order unitarity). Furthermore, we perform global Bayesian fits to all relevant constraints on these models, which also include electroweak precision and flavour observables, and highlight the complementarity between them.
This paper is organized as follows: The 2HDM is defined in Section 2. In Section 3 we list all relevant constraints and explain the fitting set-up. The results are presented in the subsequent sections, first taking into account only the Higgs signal strengths in Section 4 and the direct searches in Section 5, before combining them with the other constraints in Section 6. We conclude in Section 7. In Appendix A we explain how we treat the prior dependence of the massive parameters.

Model
The Two-Higgs-Doublet model with a softly broken Z 2 symmetry is characterized by the following scalar potential: where Φ 1 and Φ 2 are the two Higgs doublets. While writing the potential we have assumed that the scalar potential is CP conserving. Instead of the eight potential parameters from Eq. (2.1) we will use the physical parameters in the rest of this article. They consist of the vacuum expectation value v, the CP-even Higgs masses m h and m H , the CP-odd Higgs mass m A , the mass of the charged Higgs, m H + , the two diagonalization angles α and β, and the soft Z 2 breaking parameter m 2 12 . Assuming the observed scalar of mass ∼ 125 GeV at the LHC to be the lighter CP-even Higgs h, the first two of these can be treated as fixed. The rest of the scalar masses could in general even be lighter than 125 GeV, they are not necessarily in the decoupling limit [5]. Keeping in mind the discovery potential of the HL-LHC, in the following we will consider them to be in the range between 130 GeV and 1.6 TeV, that is beyond the region where the 125 GeV scalar was found. Moreover, we trade the angles α and β with β − α and tan β, since these combinations can be directly related to physical observables. All SM parameters were fixed to their best-fit values [33,34].
Neglecting the first two generations of fermions, the Yukawa part of the 2HDM Lagrangian reads as follows: In the above Lagrangian, by convention the top quark only couples to Φ 2 ; its Yukawa cou- With an unbroken Z 2 symmetry in the Yukawa sector, there are only four possibilities through which the Higgs fields couple to the bottom quark and tau lepton at tree-level. They are called type I, type II, type X or "lepton specific" and type Y or "flipped". In Table 1 we categorize the corresponding Yukawa coupling assignments. Type I Type II Type X ("lepton specific") Type Y ("flipped")

Constraints and fitting set-up
Our statistical analysis of the 2HDM is a Bayesian fit, in which the following priors are used for the previously defined parameters: Extreme tan β values outside the chosen prior are expected to be excluded due to the absence of strong 2HDM effects in certain flavour observables (see e.g. reference [35]); the aforementioned interval is a very conservative estimate. The only implicit assumption we make is that the 125 GeV scalar is the light CP-even Higgs particle of the 2HDM and that the other scalars should be heavier, yet in LHC reach. The focus of this article is on LHC Higgs observables, that is h signal strengths and searches for H, A and H + . Most details of the implementation of the corresponding observables can be found in our last article [31]. The modifications to this will be explained in the following.  Table 2. h signal strengths from Table 8, Table 13 and Figure 27 of the official ATLAS and CMS combination for Run 1 [3], based on 25 fb −1 of integrated luminosity. We neglect correlations below 0.1. The colours in the first column indicate the decay category in Figures 1 and 2.
For the signal strengths, we define µ decay production , where "production" stands for the ggF, VBF, Vh, Zh, Wh, tth or pp production channels of the h, while "decay" denotes the subsequent h decay products γγ, ZZ, W W , τ τ , bb, µµ or Zγ. 1 For the last one, only upper limits are available; we assign to this signal strength a central value of 0 and adjust the Gaussian error such that the likelihood distribution has the 95% limit at the value provided by the experimental collaborations. All h couplings are calculated at leading order: While the fermionic decays and the bosonic decays to W W and ZZ are possible at  tree-level, we apply one-loop expressions for the decays into final states including massless bosons (that is gg, γγ and Zγ) [36]. A list of the available experimental signal strength values from LHC Run 1 and 2 can be found in the Tables 2 (ATLAS and CMS combination for Run 1), 3 (ATLAS numbers for Run 2) and 4 (CMS measurements for Run 2). For the Run 2 data, we also list the corresponding integrated luminosities L. The numbers for the correlations in Table 2 can be found in the mentioned document. For Run 2, ATLAS provides correlations only for the combination of the γγ and ZZ decays; observing very similar numbers in the corresponding Run 1 data, we assume identical correlations for the γγ and ZZ final states. The correlation between µ W W tth and µ τ τ tth was extracted from Figure 17 in [41]. (We assume that the V V final state therein is dominated by W W .) Also the CMS correlations in Table 4 were reconstructed from the signal strength contours (or cross section times branching ratio contours) in the plane of VBF vs. ggF production. In Section 4 we discuss the individual  Table 4. Run 2 h signal strengths measured by CMS. The colours in the first column correspond to the ones in Figures 1 and 2, the ones in the fifth column to the amount of underlying data like in Table 3. impact of the signal strengths on the 2HDM parameters, ordered by the decay products.
Concerning the direct searches for the heavy CP-even, the CP-odd and the charged Higgs, we have updated the number of used LHC analyses from 16 in [31] to 50 Run 1 and 2 measurements in the present article. We calculate the product of the production cross section [104][105][106][107][108][109][110][111][112][113] and the branching ratio [114,115] of a specific decay, σ · B. In order to compare it with the experimental bounds, we assign Gaussian likelihoods with a central value of 0 to the ratio of the theoretical value and the observed upper limit of σ · B. This method agrees with the treatment of the upper limit of the Zγ signal strength mentioned above and coincides with our approach in [31] under the assumption that the observed upper limit does not deviate from the expected one. With no evidence for such a deviation in any of the searches, this approximation seems to be justified. The experimental input    Table 7. List of the available charged heavy Higgs searches from LHC Run 2 relevant for the 2HDM. For an explanation, see the description below Table 5.
instance (γγ)(bb) -is included in the upper limit, we list it in the table. If the final state is not included in the σ · B limits, but its information is needed to distinguish it from other searches, we write it in square brackets. The secondary decay products of one particle are combined in parentheses. In the case in which two concurring searches are available that are partially based on the same set of data, we use the limit which is derived from the larger amount of data. For instance, the latest CMS results of the searches for an H decaying via ZZ into two leptons and two neutrinos are available only for m H > 600 GeV. Lighter m H scenarios will be constrained using an older publication based on an integrated luminosity of 2.3 fb −1 . Also for the upper limits on H → hh → (bb)(bb) by CMS and on H + → tb by ATLAS we apply different searches depending on the masses. For gg → X → γγ, CMS combined their 8 and 13 TeV data; the limits are given for the 13 TeV production cross section. A detailed discussion of how the different searches constrain the 2HDM can be found in Section 5, where we show the results ordered by the decay products. Apart from the discussed tree-level Higgs observables the 2HDM scalars can also contribute to the quantum corrections of other observables, the most important ones being the electroweak precision observables, the b → sγ branching ratio and the mass difference in the B s meson system. While the implementation into HEPfit was already explained in [31], we updated the experimental values [34,116,117]. Also for the treatment of theoretical constraints we refer to [31], with two exceptions: We do not apply any constraints arising from the renormalization group evolution and define our model at the electroweak scale. And for the next-to-leading order unitarity bounds we chose the most conservative approach that appeared reasonable to us, namely requiring that the real and imaginary parts of the S-matrix eigenvalues should be between −0.5 and 0.5 and between 0 and 1, respectively. Moreover we impose perturbativity by discarding scenarios for which the one-loop contribution to these eigenvalues exceeds the tree-level term in magnitude.
As numerical set-up we use the open-source package HEPfit [118], interfaced with the release candidate of the Bayesian Analysis Toolkit (BAT) [119]. 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.

h signal strengths
In this section we show the impact of the h signal strengths on the 2HDM parameters. The fits were done with the most up-to-date experimental inputs; for a comparison with the status before EPS-HEP 2017, see [120]. The differences in the Z 2 symmetry assignment to the fermions result in a type dependent treatment of their couplings to the light Higgs boson. The signal strength of the process with a given initial state i producing an h which decays to the final state f can be written as where r x is the ratio of the 2HDM and the SM partial width of an h decaying into x and B SM (h → x) is the corresponding SM branching ratio. From this equation one can see that every signal strength depends on the 2HDM h couplings of all decay products. In Figure 1 we show the individual impact of the signal strengths with a specific final state on the β − α vs. tan β plane as well as their combination in all four types of Z 2 symmetry. We have tried to adopt the colouring scheme from Figure 14 of the Run 1 combination [3]. All contours delimit the regions allowed with a probability of 95.4%. The upper limit on µ Zγ pp is included in the combination, but not shown separately as its effect is minimal.
In type I all fermions have the same relative coupling to h: r tt = r bb = r τ τ = r µµ = cos(β − α)/ tan β + sin(β − α). This can only deviate significantly from 1 if tan β is smaller than 1 and β − α is not close to the alignment limit π/2. In these regions, the di-photon signal strengths are the most constraining ones, see the upper left panel of Figure  1. For larger tan β values, the ZZ and W W signal strengths become the most important constraints. In the combined fit to all signal strengths, the largest possible deviation of β − α from π/2 is 0.26 at 95.4% if we marginalize over all other parameters.
The upper right panel of Figure 1 shows fit results with the same inputs for type II. Here, the relative down-type fermion and lepton couplings to h are different from the top coupling, and thus the fermionic signal strengths yield more powerful constraints. But also the signal strengths with a bosonic final state become stronger because of the modifications of the loop coupling r gg and the fermionic couplings in the denominator of Eq. (4.1). Especially the W W and ZZ signal strengths constrain β − α to be very close to π/2; the largest deviation from the alignment limit in the one-dimensional fit to all signal strengths is 0.055 at 95.4%. The so-called "wrong-sign" solution for the fermionic couplings, which is represented by the lower "branches" of the individual W W , ZZ, bb and τ τ signal strength fits for tan β > 3 and β − α < 1.5, can (cannot) be excluded in the combined fit to all signal strengths with a probability of 95.4% (99.7%). These scenarios have also been shown to be incompatible with the assumption that the 2HDM of type II is stable under a renormalization group evolution up to O(1) TeV [31,121].
In type X, the h couplings of the down-type quarks agree with the ones of the top quark, but the leptonic couplings are like in type II. Consequently, the contour of the bb decays in the lower left panel of Figure 1 has a similar shape as the one of type I, while the τ τ and µµ decays behave more like in type II for large tan β for β − α < π/2. For tan β > 2 the latter two are the dominant signal strengths. For very large tan β, the wrong-sign solution of the fermion couplings is allowed at 95.4%. However, no larger deviations of β − α from π/2 than 0.069 are allowed at the 95.4% level if we combine all signal strength information and marginalize over all other parameters.
Finally, the type Y fit can be found in the lower right panel of Figure 1. Like in type II, β −α has to be very close to the alignment limit with the bosonic signal strengths being the strongest constraints. But like in type X, the wrong-sign coupling of the fermions cannot be completely excluded at 95.4% in the fit combining all signal strengths, although it is only possible for very large tan β. In this type's combined fit and marginalizing over the other parameters, β − α cannot be further away from π/2 than 0.056 with a probability of 95.4%.
As compared to the status before EPS-HEP 2017 [120], the W W , γγ and bb signal strengths have become more constraining; the latter changed drastically due to additionally released data. In type II, a small spot of the "wrong sign" branch around tan β = 3 and β − α = 1 was allowed at the 95.4% before summer 2017 and has disappeared now.
The two angles α and β define all tree-level couplings of fermions and bosons to the light Higgs h, but the loop couplings to gluons and photons are more complicated. In order to analyse their allowed ranges, we show the r gg vs. r γγ plane in Figure 2. Apart from the individual fits to the different final states and their combination like in Figure 1, we also add the contours from a fit to only Run 1 and only Run 2 data, respectively. In all types, the combined fit to all signal strengths is dominated by the bosonic decays. While the µ γγ mainly delimit r γγ , r gg is constrained also by the µ W W and µ ZZ measurements. The maximal deviation of r γγ (r gg ) from its SM value is roughly 30% (20%). The wrong-sign solution for the fermion couplings can be seen in type II, X and Y: The regions for r gg > 1 in Figure 2 contain the lower branches of Figure 1. For type II, it has been shown that the wrong-sign couplings feature increased r gg and reduced r γγ [12]. In the lower right panel of Figure 2, this "second solution" is visible between 1.1 and 1.2 for r gg as spikes in the W W , ZZ and τ τ contours for large r γγ as well as in the combined signal strength fit in both r γγ directions. Comparing all Run 1 signal strengths with all Run 2 signal strengths, one can see that generally the Run 2 data is more constraining in all types. However, the Run 1 signal strengths prefer a smaller r gg and thus determine the upper limit of the gluon coupling ratio in the combined fit to all signal strengths. This can be seen in the types II and Y.
The run time for fits with 120 million iterations was 230 ± 80 CPU hours for the di-photon final state and 70 ± 20 CPU hours for the other single channels.

Heavy Higgs searches
In the following, we will scrutinize the impact of the searches for heavy Higgs particles in all four types of a 2HDM with a softly broken Z 2 symmetry, ordered by their decay products. First we will address the fermionic decays to tt, bb and τ τ and the loop induced decays with γγ and Zγ in the final state. The searches for signals in these channels apply to both, H and A bosons. After that, we will turn towards the H specific decays into two massive vector bosons or two h bosons and the A specific channel with an h and a Z in the final state. Finally, the decays of a charged Higgs to τ ν and tb will be discussed. The limits C AZ 8 and C HZ 8 on the decays H → AZ and A → HZ were only used in the combination of all heavy searches; their impact on the m A vs. m H plane is found to be weak if we marginalize over all other parameters. The narrow width approximation will be applied throughout this section; we will comment on its validity at the end of the next section.
The grey shaded regions in the plots in this section depict the prediction of the 13 TeV σ · B for the corresponding channel without applying any theoretical or experimental constraints on the model; in other words they correspond to our priors. The black dashed lines delimit the available ranges of the σ · B for the corresponding channel when only the theoretical constraints defined in the Sec. 3 have been used in the fit. The areas within the various coloured solid lines depict the 95.4% posterior ranges of σ · B after imposing the experimental constraints from the LHC for a particular measurement. The legend of each plot refers to the channels described in Tables 5, 6 and 7. The horizontal coloured lines on the top of the panels mark the mass ranges analyzed at the LHC for each of the searches denoted in the legend. In the following plots the posterior prediction of σ · B after considering a particular direct search replicates the prior behaviour unless it deviates from rest of the posteriors for the same channel. In other words, the direct searches only have a visible impact on the 2HDM parameters if their contour is lower than the other coloured contours. Also the lower lines of the searches mostly represent the priors and do not pose any significant lower limit on σ · B.
The fits to the single experiments in this section took 60 ± 7 CPU hours with 120 million iterations.

H and A decays to tt
For H and A masses heavier than two top quarks the decay to tt is the dominant in the 2HDM, at least for moderate values of tan β. Unfortunately, a possible signal strongly interferes with the tree-level background process gg → tt. The only available experimental limits of an analysis which takes into account this interference is [122] for the 2HDM of type II. Its limits, however, constrain only a small region for tan β ≈ 1 and m H/A ≈ 500 GeV, which has been excluded by indirect constraints. Therefore, we do not take into account this direct measurement.

H and A decays to bb
The direct search for H → bb decay does not put any constraint on σ · B in all four 2HDM types considered in the analysis. In Figure 3 one can see that theoretical constraints provide a suppression of σ · B by roughly an order of magnitude compared to the fit without any constraint in type I and X. Similar to the previous case, pseudoscalar decaying to bb searches do not provide any stronger constraint on σ · B than the fit with theory constraint alone in all four types. For this search, theory constraints alone restrict σ · B by at least one order in magnitude with respect to the fit without any constraints in the parameter space analyzed for all the types. This suppression is more dominant in type II and Y compared to the other two cases. In type II and Y, in the regime m A 600 GeV the A → bb search from Run 1 suppresses σ · B compared to the fit without any constraints but it remains sub-dominant or at most of similar strength to the fit with theory constraint alone.

H and A decays to τ τ
In the upper panel of Figure 4 we show that the H → τ τ searches suppress the σ · B limit by at least one order of magnitude compared to the fit with theory constraints alone in the regime where the heavy Higgs mass is below 250 GeV. In this regime the strongest constraint comes from Run 1 data and the suppression of σ · B is more pronounced in the types I and Y. Theory constraints alone restrict σ · B by roughly an order of magnitude compared to the fit without any constraints for type I and Y, whereas in type II and X, theory constraints raise the lower limit of σ · B for m H > 500 GeV. This can be understood as a sensitivity of the fit to fine-tuned scenarios with extreme tan β values, which are disfavoured in a fit to the theoretical bounds.
From the fit without any constraints in the lower panel of Figure 4 we see that the predicted ranges of σ · B for the pseudoscalar decaying to τ τ are quite narrow for type II and X compared to the other two 2HDM types. The theory constraints yield a suppression by at least one order of magnitude for σ · B compared to the fit without any constraints in the types I and Y whereas for type II and X the theory constraints push up the lower limit on σ · B by one order of magnitude compared to the fit without any constraints for m A > 350 GeV. The direct search limits for A → τ τ suppress σ · B by roughly one to two orders of magnitude compared to the fit with theory constraints alone for all the types except for type Y as long as m A 300 GeV. In type Y, the experimental upper limits on σ · B are stronger than the theoretical ones for pseudoscalar masses between 200 and 400 GeV, in type II even up to 1 TeV. Scenarios with very light A and σ · B values between 10 −5 σ · B 10 pb seem also to be excluded at 95% by the prior. However, since we have no experimental data on this region, this prior dependence is not an issue here.

H and A decays to γγ
The theory constraints on H/A → γγ suppress σ · B by one to three orders of magnitude compared to the fit without any constraints in all four 2HDM types, see Figure 5. Direct searches for a heavy CP-even Higgs decaying to two photons constrain σ · B by roughly one order of magnitude compared the fit with theory constraints for m H 250 GeV in all types. The searches in the di-photon decay channel of a pseudoscalar Higgs yield a suppression of σ · B by one to three orders of magnitude compared to the fit with theory constraints for m A 600 GeV for all four types considered. In the types II and Y, we observe again that certain intermediate σ · B regions for low m A are disfavoured by the prior.

H and A decays to Zγ
We see from the top panel in Figure 6 that for H → Zγ, theory constraints suppress σ · B by one to four orders of magnitude compared to the fit without any constraints for all types. In all four types, direct search for this channel does not provide any constraint on σ · B except for a very small window below m H 250 GeV, but it remains sub-dominant compared to the fit with theory bounds.
Similar to the heavy CP-even Higgs case, theory constraints yield a suppression of σ · B by one to two orders of magnitude for the decay A → Zγ compared to the fit without any constraints in all four types. Direct searches for this channel provide a suppression of σ · B by an order of magnitude compared to the fit with theory constraints in the mass window 250 m A 350 GeV in all four types. Again, some parts of the σ · B region for light A masses are disfavoured by the prior in the types II and Y.

H decays to ZZ or W W
The heavy Higgs decays to massive gauge bosons can be divided into searches for ZZ and W W , but the relative coupling of H to two vector bosons V V = ZZ, W W is universal and type independent. However, the production of the H differs between the types. We show the H → ZZ channels in the upper panel of Figure 7 and the searches for H → W W as well as the combined searches for H → V V in its lower panel. The σ · B are constrained by the theoretical bounds in the decoupling limit, where m H > 600 GeV. The direct LHC searches for this channel yield a strong suppression of σ · B by one to three orders of magnitude compared to the fit with theory constraint in the mass regime 150 m H 800 GeV (150 m H 700) for the ZZ (W W ) channel. For the ZZ searches the m H 250 GeV region is constrained by Run 1 data whereas Run 2 data determine the dominant limits for the rest of the region. For the W W searches, Run 1 data dictate the limit until m H 600 GeV and the high mass range is dominated by Run 2 data. Additionally the C V V 8 search for H → V V severely constrains σ · B in the 200 m H 250 GeV region for type II and Y.

H decays to hh
In the upper panel of Figure 8 we show that for the H decaying to two h bosons, theory constraints already yield a strong suppression of σ · B compared to the fit without any   constraints in all four types. Direct searches suppress σ·B by at most one order of magnitude compared to the fit with the theory constraints in the mass range 200 m H 600 (500) GeV in type I and X (type II and Y). The main constraints stem from Run 1 data and the Run 2 searches for hh resonances decaying to two photons and two bottom quarks. Although different searches at Run 1 and Run 2 continue to constrain the σ · B for this channel up to m H 1200 GeV, they remain sub-dominant to the limit from the fit to the theory constraints.

A decays to hZ
The searches for A decaying into hZ are shown in the σ · B vs. m A planes in the lower panel of Figure 8; more precisely, they are projected onto the σ · B of the decay to (bb)Z. As compared to the fit without any constraints, theory constraints effectively are important in the decoupling limit in all types as well as for m A 300 GeV in type I and X. Direct searches for this channel yield a strong suppression (one to three orders of magnitude) of σ ·B compared to the fit with theory constraints for all four types as long as m A 800 GeV. In all types, the strongest bounds come from searches in the A → hZ → (bb)Z channel; in type X we additionally observe that the A → hZ → (τ τ )Z limits yield the most important constraints for m A between 200 GeV and 300 GeV.

H + decays
At the present, charged Higgs decaying into τ ν does not provide any constraint on σ · B in the 2HDM's under consideration, see the upper panel of Figure 9. The theory constraints yield a suppression of σ · B by roughly one order of magnitude as compared to the fit without any constraint for all types but type X as well as for all types if the charged Higgs decays into tb (lower panel of Figure 9). In the latter case, the direct searches from Run 1 (Run 2) provide even stronger limits on σ · B between 180 GeV and 600 GeV (1 TeV), with only small differences between the 2HDM types.

All heavy Higgs searches
In Figure 10 we show the available parameter space for 2HDM masses and angles from the fit where the heavy Higgs searches are taken into account. The region inside the various coloured patches are disfavoured by the corresponding search category denoted in the legend. The central areas inside the solid grey line mark the 95.4% allowed regions when all heavy Higgs searches are considered in the fit, including also C AZ 8 and C HZ 8 . In the panels in the first row the black dashed lines mark the limit from the fit to theory constraints only. The combination of all H/A/H + searches is represented by the orange/blue/green dashed contours. The channels described in the previous sub-sections which are not constraining or very weakly constraining in this mass vs. angles plane are not shown in the figure.
From the first row of Figure 10 we can see that the region around β − α = π/2 remains unconstrained in all four types of 2HDMs when all the heavy Higgs searches are taken in account. From moderate to high masses, the di-Higgs channels dominate the excluded regions in all four types whereas H → τ τ , H → γγ and H → V V are the most important constraints below the hh threshold. In type I and X, the final exclusion region is mainly   dominated by the heavy Higgs to two light Higgs channel, while in type II and Y this is only true if β − α > π/2; for β − α < π/2 the final constraint is weak except for m H ≈ 250 GeV. Although the di-photon searches alone disfavour m H 600 GeV for regions near β − α π/2 this region is allowed when considering all the heavy Higgs searches in the fit. The H → V V decays only are susceptible to the tan β 1 region and for m H 800 GeV in all four types; this again is an effect of fine-tuning. In the displayed m H vs. tan β ranges, all the H searches become ineffective for m H 1000 GeV in type II and Y, whereas in type I and X the tan β 1 regions remain inaccessible even if the heavy Higgs mass is as large as 1500 GeV. This is an effect of the combination of all heavy Higgs searches, in which scenarios with small tan β are sampled less by the fitter. This feature is not worrisome, because these regions are also suppressed by the flavour observabels as we will show in the next section.
In the pseudoscalar mass vs. tan β planes we see that the di-photon channel constrains low tan β and m A up to 600 GeV for all four types. The A → hZ channel can exclude tan β values up to 10 and m A almost as heavy as 1000 GeV in all four types. The exclusion also applies for the large tan β regions in the types II and Y. The next most important channel for the pseudoscalar searches in type II is the A → τ τ channel, which efficiently excludes high tan β regions for m A as large as 1000 GeV. In type X this channel is susceptible to tan β 20 and m A 400 GeV. The exclusion from this channel is weaker for type I and Y. The contour for all pseudoscalar searches is mainly dominated by the hZ channel in type I and Y, and a combination of hZ and τ τ in type II and X. Combining all the pseudoscalar Higgs searches the present data constrain certain regions where m A 1000 GeV. For type X, the m A < 400 GeV region remains available only if tan β 10. In the same mass regions a very narrow range of intermediate tan β remains accessible for type II when all constrains are taken into account. The bound on large pseudoscalar masses is similar to the heavy Higgs case when all constraints are taken into account.
As described in the previous sub-section, the main channel which constrains the charged Higgs mass is H + → tb. The exclusion region of this channel is shown in red in the last row of Figure 10. From the figures we see that the present searches for the charged Higgs mass can only constrain the regions with tan β 1 and m H + 1000 GeV in all four types. The inclusion of all the searches in the fit yields stronger constraints in the m H + vs. tan β plane than the fit to all charged Higgs searches only. In type II and Y this is due to the H → hh → 4b searches, which are particularly sensitive to large tan β and disfavour certain regions featuring tan β 15. For type I and X it is more difficult to pinpoint one particular channel for the seeming exclusion of tan β values up to 2 for charged Higgs masses above 1 TeV, but in the next section we see that these bounds can be relaxed once we take into account also other constraints.

Combination of all constraints
After discussing the individual effects of the h signal strengths and the searches for H, A and H + on the 2HDM, we want to confront these constraints with the other bounds on the parameters. In Figure 11 we copy the information about all heavy Higgs searches from Figure 10 in the mass vs. angle planes, and add the bounds from signal strengths and theoretical constraints to the m H vs. β −α planes and the impact of the flavour observables to the planes with tan β. Finally, the global fit to all constraints is represented by the grey regions.
The signal strength bounds do not depend on the masses of the heavy Higgs particles; their limits on β − α differ for each type, see Section 4. The theory conditions force the 2HDM's into the alignment limit for m H > 600 GeV, decoupling the heavy Higgs particle from physics around the electroweak scale. The type dependence of this effect is negligible as it only enters via sub-leading Yukawa terms in the beta function parts of the NLO unitarity conditions. Besides the obvious consequences for certain constellations of m H and β − α, the combination of the signal strengths with the theory constraints also disfavours large values of tan β. This is because extreme values for the latter result in a destabilization of the unitarity conditions [31]. Specific combinations of the 2HDM angles can still fulfil the theoretical constraints, but these solutions are highly fine-tuned and thus have a low posterior probability. The flavour constraints have been discussed many times in the literature; the summary is that in all types the B s mass difference sets lower limits on tan β (at least for masses within the reach of the LHC), while the branching ratio of b → sγ processes enforces m H + 580 GeV at 95% C.L. in type II and Y [123]. In combination with theory and electroweak precision bounds, these limits on the charged Higgs mass can be translated to lower limits on the neutral masses. (The individual impact of ST U will be explained below as it is not visible in these two-dimensional projections of the parameter space.) The combination of all constraints is more intricate than the naive superposition of all individual bounds. First of all, we should mention that it also depends on the prior we choose for the masses: while the direct experimental observables depend on the masses of the heavy Higgs bosons, the theoretical bounds and the loop-induced effects are only sensitive to the mass squares. Since we want to combine both, we need to decide whether we want to use a flat prior for the masses or the mass squares. A detailed discussion can be found in Appendix A. The contours we show here are a superposition of a fit with flat mass priors and a fit with flat mass square priors in order to be as conservative as possible. In the m H vs. β − α planes of type I and X, the combination of signal strengths and theory only leave a very small strip around the alignment limit of β − α = π/2. Heavy Higgs searches additionally exclude m H < 380 GeV in type X. The reason for this are not only the H → τ τ searches as observed in Figure 10, as they are similar to the bounds in type I, but mainly it is an interplay of the strong bounds of the A searches for low tan β and the above-mentioned exclusion of large tan β due to signal strengths and theory, which together disallow m A < 400 GeV. This bound translates to a limit on m H using the unitarity and electroweak precision constraints, since these bounds delimit the mass splittings, see below. Also in the type II and Y planes we see a lower limit of m H > 550 GeV, which in this case derives from the lower bound on m H + from the b → sγ measurements and the fact that the mass difference m H − m H + cannot be too large. The absolute maximal deviation of β − α from π/2 is 0.03 in type I and 0.02 in the types II, X and Y. (This corresponds to 1 − sin(β − α) < 5 · 10 −4 and < 2 · 10 −4 , respectively.) In the m H vs. tan β planes one can see that the fine-tuning for large tan β scenarios disfavours these regions and pushes the allowed contours towards smaller values of tan β. Only in type I and for m H < 350 GeV, the heavy Higgs searches have a visible impact on this plane, excluding tan β 2.5. More or less the same holds for m A vs. tan β, where in type I tan β < 3 is excluded by direct A searches if m A < 350 GeV. We already mentioned above that in type X, the interplay between fine-tuning and A searches sets a lower limit of 400 GeV on m A . Having a look at the m H + vs. tan β planes, one can see the lower bounds on the charged Higgs mass in type II and Y from b → sγ, which we quantify to be 600 GeV in our fit. Also here, we observe that large tan β values are disfavoured and the posterior regions are shifted towards small tan β.
The fit only to electroweak precision data does not exclude any region in the twodimensional mass vs. angle projections in Figure 11. What it does constrain are the mass differences between H, A and H + . That is why in Figure 12 we show the difference between the pseudoscalar and charged Higgs mass, once depending on the H mass (left column) and once against the m H − m H + difference. In the m A − m H + vs. m H planes the dominant constraints come from the theory bounds, at least in the decoupling limit m H > 600 GeV. The ST U pseudo-observables are stronger if m H < 600 GeV and m H + > m A . In type I, they yield a lower bound on m A − m H + in the global fit for m H < 250. If we look at the m A − m H + vs. m H − m H + planes, we observe that here the oblique parameters are the strongest constraint on m A −m H + if the charged Higgs mass is larger than m H . Combining all constraints and marginalizing over all other parameters, we obtain the following ranges for the mass differences allowed with a probability of 95%: For all heavy Higgs search limits, we implicitly assumed the narrow width approximation. In a simultaneous fit to all constraints except for these direct search limits we find that with a probability of 95% the decay widths of H, A and H + never exceed 5.5% of the mass of the particle in the types II and Y. For masses below 1 TeV the maximal decay widths are less than 3.5% in these two types. In type I and X and for φ = H, H + the fits yield Γ φ /m φ < 3.5% (< 5%) if m φ < 1 TeV (< 1.5 TeV). Only the ratio Γ A /m A can reach 7% for m A ≈ 550 GeV, but in the decoupling limit m A > 600 GeV similar bounds apply as for H and H + . We would like to stress here that all these 95% limits are maximally allowed values and that in a typical 2HDM scenario the widths are significantly smaller. Therefore we conclude that the narrow width approximation is a reasonable choice for 2HDM scenarios. Finally, addressing the last variable of our chosen parametrization, we also observe limits on the soft Z 2 breaking parameter m 2 12 . The upper limits strongly depend on the maximally allowed physical Higgs masses and are around (1 TeV) 2 in all types. Due to the lower mass limits on the physical Higgs particles in the types II, X and Y, we also observe that m 2 12 is limited from below, the respective minimal values being (280 GeV) 2 , (170 GeV) 2 and (240 GeV) 2 . Only in type I an unbroken Z 2 symmetry is still compatible with all constraints.
The run time for the global fits to all constraints with 240 million iterations was 550 ± 130 CPU hours.

Conclusions
In all four 2HDM types with a softly broken Z 2 symmetry we have presented global fits to the most recent data.
Focussing on the latest measurements from LHC, we have showed explicitly how the individual signal strengths affect the leading order h couplings at tree-level and at oneloop level. Combining all information about the signal strengths, we find that the quantity |β − α − π/2| cannot exceed 0.26, 0.055, 0.069 and 0.056 in the types I, II, X and Y. The one-loop couplings of the h to gluons and photons cannot differ by more than 20% and 30%, respectively, relative to their SM values.
In order to systematically discuss the searches for H, A and H + , we have categorized them according to their decay products and have compared the exclusion strength of the single available ATLAS and CMS analyses on the production cross section times branching ratio, depending on the masses. We have then combined all decay categories and have showed their impact on the 2HDM masses and mixing angles. For m H below 1 TeV we observe strong bounds on "extreme" values for the angles, that is if β − α is very different from the alignment limit π/2 or if tan β is smaller than 1 or larger than 10. The exact limits depend on the model type and m H . Also the LHC searches for pseudoscalars severely constrain the 2HDM parameters: For m A < 1 TeV, the lower limits on tan β reach values of around 10 in the types I and X. In the types II and Y, these limits are weaker, but there are also mass dependent upper limits. The bounds from charged Higgs searches are less constraining in comparison; nevertheless, they also start to be stronger than the indirect constraints in the regions with low m H + and low tan β.
Finally, we have confronted the LHC h signal strengths and heavy Higgs searches with all other relevant indirect constraints from theory and experiment. In detail, we have showed how stability and unitarity constraints and B physics observables set mass dependent limits on the 2HDM angles and on the differences between the heavy Higgs masses, while electroweak precision data only affect the latter. We have compared all different sets of constraints and have showed the results in the mass vs. angle planes as well as in the m A − m H + vs. m H (−m H + ) planes for all four types of Z 2 symmetric 2HDM's together with the simultaneous fit to all constraints. In this global fit we find the following 95% probability limits on the 2HDM parameters marginalizing over all other parameters: |β − α − π/2| cannot be larger than 0.03 in type I and 0.02 in the other types. In type II and Y, m H > 700 GeV, m A > 750 GeV, m H + > 740 GeV and m 2 12 > (240 GeV) 2 , while we observe lower mass limits of m H > 450 GeV, m A > 500 GeV, m H + > 460 GeV and m 2 12 > (170 GeV) 2 for type X. For the latter, it is the first time that a statistically significant lower limit on the massive parameters has been observed in a global fit for the analyzed mass ranges. Also, if we discard particularly fine-tuned scenarios, only the following ranges for tan β are allowed for masses below 1.6 TeV: [0.93; 10.5] in type I, [0.93; 5.0] in type II, [0.93; 5.2] in type X and [0.91; 5.6] in type Y. However, the upper limits are no strict bounds and have to be taken with a grain of salt. Moreover, we can put type dependent upper limits of order of 100 GeV on the differences between m H , m A and m H + , and thus kinematically exclude all decays of H or A into another heavy Higgs particle. As a consequence, the decay widths of H and H + cannot exceed 5.5% of their mass in all types, at least as long as we consider masses below 1.5 TeV. While in the types II and Y we see a similar limit for the A decay width, it can amount to up to 7% of m A in the types I and X.

Acknowledgments
We thank Enrico Franco, Ayan Paul, Maurizio Pierini and Luca Silvestrini for useful discussions. The fits were run on the Roma Tre Cluster. We are grateful for the availability of these resources and especially want to thank Antonio Budano for the support. We also thank Jose

A Prior dependence of the massive parameters
In this appendix we discuss the prior dependence of our analysis when all the constraints have been taken into consideration. In Figure 13 we compare the allowed parameter space in the mass versus angles planes and the mass difference versus mass (difference) planes from the fit with flat mass priors (blue solid and dashed curves) with the fit with flat mass square priors (red solid and dashed curves). At a first glance one can see that the fit with flat mass square priors prefers high mass regions and raises the lower limit on the masses by O(100) GeV compared to the fit with flat mass priors.
It is well known that Bayesian statistics do not provide a unique rule to determine the prior distribution and in general the posterior distribution is a prior dependent quantity. A thumb rule would be to choose a flat prior for the parameter on which the observables depend linearly. For example, if an observable quadratically depends on a particle mass, one would choose a flat mass square prior. Unfortunately, in the 2HDM the theoretical and indirect experimental constraints depend on the mass squares, whereas the direct experimental observables dependent on the masses. Not having sufficiently constraining Figure 13. The mass prior dependence of the fit with all constraints is shown in the m H/A/H + vs. tan β planes, in the m H vs. β − α planes (corresponding to Figure 11) and in the m A − m H + vs. m H and m A − m H + vs. m H − m H + planes (corresponding to Figure 12), from top to bottom. data makes assigning the mass priors in the fit with all constraints is a delicate task. This would not be problematic if the observables were measured with a high precision; in fact, we can see that in the types II and Y the difference between the two priors are considerably smaller as there are strong lower mass limits. In order to be as conservative and prior independent as possible, we decided to combine the 95.4% regions for both priors: The light grey contours for the fits with all constraints in Figures 11, 12 and 13 are obtained by superimposing a fit with flat mass priors and a fit with flat mass square priors. The corresponding numerical results mentioned in Section 6 are based on the more conservative fit; for instance, the limits for masses and mass differences were extracted from the fit with flat mass priors, while the upper limits on the decay widths were larger in the fits using flat mass square priors.   [94] ATLAS collaboration, Search for Higgs boson pair production in the bbγγ final state using