Global fits in the Aligned Two-Higgs-Doublet model

We present the results of a global fit to the Aligned Two-Higgs Doublet Model, assuming that there are no new sources of CP violation beyond the quark mixing matrix. We use the most constraining flavour observables, electroweak precision measurements and the available data on Higgs signal strengths and collider searches for heavy scalars, together with the theoretical requirements of perturbativity and positivity of the scalar potential. The combination of all these constraints restricts the values of the scalar masses, the couplings of the scalar potential and the flavour-alignment parameters. The numerical fits have been performed using the open-source HEPfit package.


Introduction
The particle content of the Standard Model (SM) was confirmed in 2012 with the discovery of the Higgs boson [1,2]. However there are still some aspects of Nature that can not be explained by the SM alone and leave room for new physics (NP). In particular, scalars transforming as doublets under the SU(2) L group, and therefore satisfying the successful mass relation M W = M Z cos θ, are appropriate candidates for building extended electroweak models.
The simplest of these extensions is the two-Higgs doublet model (2HDM), containing a second Higgs doublet with the same quantum numbers as the SM one. In order to avoid dangerous flavour-changing neutral-current (FCNC) transitions, the usual implementations of the model [3][4][5] assume specific discrete Z 2 symmetries to constrain the Yukawa sector [6], so that only one scalar doublet can couple to a given right-handed fermion.

The Aligned Two-Higgs-Doublet model
Let us consider the SM extended with a second complex scalar doublet of hypercharge Y = 1 2 . In general, the neutral components of both doublets can acquire vacuum expectation values. However, making a global SU(2) transformation in the scalar space spanned by the two doublets, it is always possible to work in the so-called Higgs basis, , , (2.1) where only one doublet has non-zero vacuum expectation value, with v = ( √ 2 G F ) −1/2 ≈ 246 GeV. The field Φ 1 plays the role of the SM Higgs doublet with G 0 and G ± the electroweak Goldstone bosons. The scalar spectrum contains five degrees of freedom: the charged scalars H ± and three neutral fields S i .
The quadratic terms in the potential determine the physical scalar masses [17]: 3) where and (2.6) A = S 3 is a CP-odd neutral scalar, while the CP-even neutral mass eigenstates are linear combinations of S 1 and S 2 , (2.8) The couplings of a single neutral scalar with a pair of gauge bosons are identical to the SM ones, with the field S 1 taking the role of the SM Higgs. Therefore (V V = W + W − , ZZ), The complete list of gauge couplings and scalar interactions can be found in ref. [17].
In terms of fermion mass eigenstates the Yukawa Lagrangian reads: where all fermionic fields are written as 3-dimensional flavour vectors, M f (f = d, u, ) are the diagonal mass matrices and V CKM is the usual Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix. In general, the Yukawa matrices Y f of the second doublet are not related to the fermion mass matrices and their elements can take arbitrary values, yielding FCNCs which are tightly constrained phenomenologically [52]. The dangerous FCNC transitions can be easily avoided at tree level, imposing that only a single flavour structure is present for each right-handed fermion, i.e., that the Yukawa matrices are aligned in the flavour space [8],

Fit setup and constraints
For our analysis we consider the CKM matrix as the only source of CP violation. Thus, we are assuming that the couplings of the scalar potential in eq. (2.2) and the alignment parameters in eq. (2.11) are real. The parameter space of the A2HDM is then characterized by twelve real quantities: the three alignment parameters and nine degrees of freedom in the scalar potential which we choose to be v, the four scalar masses, the CP-even mixing angleα and the quartic couplings λ 5,6,7 . Two inputs are already empirically determined: the vacuum expectation value and the Higgs mass m h = 125.10 ± 0.14 GeV [53]. 1 The numerical values of the relevant SM parameters entering the fits are compiled in table 1.
Our fits have been performed with the open-source HEPfit package [55,56], 2 which uses a Markov-Chain Monte-Carlo implementation based on the Bayesian Analysis Toolkit [57]. We assume the following priors for the fitted parameters: The priors of the remaining CP-even scalar mass depend on the scenario studied. Light (heavy) scenario refers to the case in which the observed Higgs with a mass around 125 GeV (h) is the lightest (heaviest) CP-even scalar of the model. In this paper we will focus on the light scenario, selected with the boolean flag SMHiggs set to true, 3 and adopt as mass priors for the non-SM Higgs (H) Nevertheless, in section 5 we will also discuss the implications of our fitted data set on the heavy scenario, selected with the boolean flag SMHiggs set to false, taking as mass priors A more detailed analysis of the heavy scenario, including additional data from light scalar searches is deferred to a future work. The scalar masses are chosen in a range such that they are relevant for the future LHC searches. The selected priors for the scalar potential parameters λ i are conservative, since larger values are excluded by theoretical constraints. The mixing angleα is varied in its full domain, and the alignment parameters ς f (f = u, d, ) are varied within their perturbative ranges, i.e., √ 2 ς f m f /v ≤ 1. Bayesian statistics does not provide an unambiguous way to determine the prior distributions. A rule of thumb would be considering as flat priors the ones appearing linearly in our observables. However, for the mass parameters this does not give a unique choice: while direct searches depend linearly on the heavy scalar masses, loop-induced processes appearing in flavour observables and in the Higgs signal strengths depend on the masses squared. To avoid a possible bias in the choice of these priors, we have performed fits with two different mass parametrizations. These two choices of mass priors are selected with the boolean flag use_sq_mass. If it is set to true (false) squared (linear) mass priors are used. The effect of the choice of mass priors will be commented in the cases of interest. When the choice of the mass priors is irrelevant, squared mass priors will be used.
The global fit includes the theoretical and experimental constraints discussed below. Theoretical constraints, electroweak precision observables and some of the Higgs constraints and direct searches were already included in the THDM model class and have been adapted

JHEP05(2021)005
to the more general A2DHM case. The relevant flavour observables and the more recent experimental information on direct searches and Higgs signal strengths have been implemented in the HEPfit package for this analysis.

Theoretical constraints
To assure that the scalar potential is bounded from below, one must impose the following positivity constraints on the quartic couplings λ i [58,59]: These necessary conditions restrict the allowed pattern of scalar masses. By imposing perturbative unitarity of the S-matrix we avoid that a given combination of parameters results in a too large scattering amplitude that violates the unitarity limit at a given perturbative order. Thus, we are actually requiring that the perturbative series does not break down. Here, unitarity is enforced for two-to-two scattering of scalar particles at leading order (LO), using [60] where a (0) j are the tree-level contributions to the j partial wave amplitude. For the highenergy scattering of scalars, only the S-wave amplitude (j = 0) is relevant at LO. The corresponding matrix of partial wave amplitudes is given by (3.6) and the a (0) j are the eigenvalues of a 0 . Again, these conditions are relevant to constrain the scalar potential parameters λ i .

Electroweak constraints
The electroweak precision observables (EWPOs) measured at LEP and SLC are also included in the analysis. Since the choice of nuisance parameters does not affect significantly the results, we employ best-fit fixed values for the SM inputs M Z , m t , α s and ∆α (5) The study of the oblique parameters S, T and U [61][62][63], which are very sensitive to the scalar mass splittings, is not enough to disentangle the A2HDM contributions because of the presence of additional Z-vertex corrections [64,65]. The most relevant ones are the quantum corrections to Γ(Z → bb), which are enhanced by the large value of the top-quark mass [66][67][68]. We take this into account by making first a combined fit of EWPOs, excluding the ratio R b ≡ Γ(Z → bb)/Γ(Z → hadrons) [68,69]. The updated results of this fit can be seen in table 2, which updates the analysis of ref. [70]. These allowed ranges for the oblique parameters are then used, together with the measured value of R b , to constrain the A2HDM.

Higgs constraints
The Higgs signal strengths are defined as the ratio of the production cross section σ i times the branching ratio B f , over the SM prediction, for a given production channel (i = ggF, VBF, VH, ttH) and decay mode (f =bb, γγ, µ + µ − , τ + τ − , W W, Zγ, ZZ), where r i,f are the ratios of the production cross section σ i and decay width Γ f , respectively, with respect to their SM predictions. The signal strengths are calculated in the narrow-width approximation and depend on the alignment parameters, the mixing angleα and the scalar potential parameters. The input used contains LHC data (Run I and II) from the ATLAS and CMS collaborations, and data collected by D0 and CDF at the Tevatron. The data entering our fit are detailed in appendix A (table 4).
Information about heavy Higgs searches of ATLAS and CMS, both at Run I and II, is summarized in tables 5, 6, 7 and 8, also in appendix A. The analyses provided are quoted as 95% upper limits, for different production and decay channels, on either σ · B or (σ · B) / (σ · B) SM , as functions of the resonance masses in the narrow width approximation.
Since low-energy constraints are not considered in this work, direct searches from LEP are not included in the fits. Upcoming direct searches from LHC can be easily added and the fits shown below can be updated.

Flavour constraints
Since most of the standard CKM fits assume the SM and this would not be consistent with the study of NP, the choice of the CKM parameters is subtle. To avoid inconsistencies, a fit to the CKM entries is performed. V ud is extracted from superallowed (0 + → 0 + ) nuclear β decays [71]. Given the very small value of V ub , this fixes V us ≈ λ through CKM unitarity. 4 |V ub | and V cb ≈ Aλ 2 are obtained by combining exclusive and inclusive measurements of b → uν and b → cν transitions [75]. Finally, the apex (ρ,η) of 4 Owing to a recent recalculation of the nucleus-independent radiative corrections to superallowed nuclear β decays [72,73], the PDG 2020 [53] value of V ud is about 2σ smaller than the one quoted in the PDG 2018 compilation [74], which implies a large (> 3σ) violation of unitarity in the first row of the CKM matrix. If confirmed, this violation could not be accommodated within the A2HDM where the unitarity of the CKM matrix is exact. Improved estimates of radiative corrections are needed to resolve this issue. Meanwhile, we have adopted the PDG 2018 value of V ud that fits better with the kaon determination of Vus. the unitarity triangle is determined with the additional information of the ratio |V td /V ts |, extracted from ∆M Bs /∆M B d [75] that is not sensitive to charged scalar contributions [13]. The CKM inputs obtained in this way and later used in our global fits are summarized in table 3 and figure 1. 5 Charged-scalar exchanges contribute to neutral meson mixing through one-loop box diagrams [13,77,78]. The corrections induced by virtual top quarks are quite sizeable, specially for ∆M B s,d and ε K , and provide strong constraints on |ς u | (also R b ) as function of M H ± . The weak radiative decay B → X s γ [13,[79][80][81][82][83][84] gives also important correlated constraints on ς u and ς d , specially for large values of |ς u ς d |. The region ς u ς d < 0 is actually excluded, except for very small values of the alignment parameters [13]. NNLO corrections [85][86][87] are quite relevant for this observable and should be taken into account.

JHEP05(2021)005
A complete one-loop calculation within the A2HDM of the decay B s → µ + µ − was performed in [88,89]. This observable involves the b → sµ + µ − four-fermion operators O 10 , O S and O P . The decay amplitude receives contributions from both charged and neutral scalars, and provides complementary information on the alignment parameters ς u,d, and the scalar masses. It also includes small contributions from higher-order FCNC local interactions, needed to reabsorb UV divergences, which are assumed to be negligible here. A study of these effects can be found in [10]. Our fits include the constraints from B s → µ + µ − , B → X s γ and ∆M B s,d (and R b , which is discussed together with electroweak precision observables).

Figure 2.
Left panel: allowed regions for the scalar mass splittings coming from theoretical constraints at 100% probability (blue), from EWPOs at 95.5% probability (in orange, squared mass priors and in light purple, linear priors), and combining all constraints at 95.5% probability (linear mass priors in purple, squared mass priors in red and squared mass priors with lower masses in brown). The "All constraints" contains only the right-sign branch discussed in section 4.5. Right panel: two-dimensional bounds on the λ 5 , λ 6 and λ 7 parameters of the potential resulting from imposing theoretical constraints (blue, 100% probability), and considering all constraints (in dark red 95.5% probability, in red 68% probability).
Finally, the muon anomalous magnetic moment, calculated within the A2HDM in refs. [90,91], is of interest because it shows a deviation with respect to the SM that, if confirmed, would strongly constrain the leptonic alignment parameter ς . Its implications will be discussed in section 4.

Results: light scenario
In this section, we present our main results, obtained in the light scenario which assumes that the observed SM Higgs is the lightest CP-even scalar of the model. The complementary possibility (the observed scalar is the heaviest) will be briefly discussed in section 5. We analyse first the separate implications of the different types of constraints, before combining all of them into a final global fit to the data.

Theoretical constraints
Perturbative unitarity and positivity of the scalar potential set strong limits on the scalar masses and the quartic parameters λ i . The mass differences among H, A and H ± are strongly constrained, as shown by the allowed blue regions in the left panel of figure 2 (electroweak and combined constraints, also present in the plots, will be commented later in sections 4.2 and 4.5):  There is a clear correlation between the masses of any two scalar particles: a large mass for one scalar implies that the other scalar mass is also restricted to be large. This effect is stronger for higher values of the scalar masses. The allowed mass splittings decrease as the average mass scale increases. This is easily understood from the scalar mass relations in eqs. (2.3) and (2.4), since M H ± , M H and M A become degenerate in the limit µ 2 λ i v 2 . The impact of the assumed mass priors is further studied in figure 3. Both linear (light orange) and squared (dark orange) options give rise to allowed regions with similar shapes, although they are larger for the squared priors.
The theoretical constraints also restrict the allowed ranges of the scalar quartic couplings. Two-dimensional plots in the space (λ 5 , λ 6 , λ 7 ) are shown in the right panel of figure 2, which displays the correlations among these three parameters of the scalar potential. The positivity relations imply bounds on |λ 5 | and |λ 6 + λ 7 | because these two quantities induce negative contributions to the two last conditions in eq. (3.4). This implies an anti-correlation between λ 6 and λ 7 , which is clearly manifest by the allowed blue area in the λ 6 −λ 7 plane. The blue regions in figure 2 satisfy the bounds derived in previous works [4].

Electroweak constraints
The EWPOs restrict the individual masses of the scalar particles in the low-mass range, and are very useful to constrain their mass splittings. The oblique parameters are very sensitive to the scalar mass differences, which results in strong limits for the masses. This can be clearly observed in figures 4 and 7. The information from EWPOs complements in a very useful way the theoretical constraints discussed before.
The allowed regions obtained from EWPOs present a strong dependence on the mass priors, as can be seen in figure 4, which displays the limits resulting from different choices of mass priors. Independently of the priors, large values for the masses and small splittings are favoured. The light and dark blue regions show a very strong dependence on the assumed ranges for mass-squared priors. If the scalar masses are varied until 1500 GeV, masses below approximately 750 GeV are not allowed at a 68% probability, while if a lower range below 1000 GeV is adopted, scalar masses of 500 GeV are allowed at the same probability. The same tendency is observed if lower mass regions are chosen. If the fit is repeated with linear mass priors (purple regions), masses as low as 10 GeV are allowed for the charged and CP-odd neutral scalars, at a 68% probability. In this case, the dependence on the input mass ranges is also weaker. The allowed regions for linear mass priors up to 1000 GeV (1500 GeV) are indicated in dark (light) purple colour.

JHEP05(2021)005
The orange and light purple regions in the left panel of figure 2 display the constraints from EWPOs on the scalar mass splittings with squared and linear priors, respectively. These allowed regions have been obtained varying the mass priors in their full range up to 1500 GeV.

Higgs signal strengths
Since the measured Higgs signal strengths are consistent with the SM, within their current uncertainties, the gauge and Yukawa couplings of the SM-like Higgs boson should be close to the SM limit. In particular, the measured data on the W W * and ZZ * decay modes imply that cosα cannot deviate much from one and, therefore,α should be small. A similar comment applies to the Hff interactions. However, most Higgs observables are not sensitive to the signs of the Yukawa couplings and, therefore, the LHC data only require the modulus of |y h f | − 1 to be smaller than about 0.1-0.2. This gives two different types of solutions for the Yukawa couplings: there will be a broad range of allowed values of ς f withα ≈ 0, corresponding to y h f ≈ 1, and another region with somewhat larger values of the mixing angle corresponding to y h f ≈ − 1.  amplitude of the Higgs that involves one-loop contributions from virtual W ± , t and H ± . Assuming that the charged-scalar correction is small, the measured H → γγ signal strength determines the relative sign between y h u and g hW W to be positive. Therefore, only the regioñ α ς u 1 is allowed in this case. In the following we will distinguish among the two different possibilities: the "rightsign" solution, corresponding to y h d, ≈ 1 and the "wrong-sign" one corresponding to y h d, ≈ − 1. The former was previously analysed in the A2HDM [23] and, more recently, in the particular case of Z 2 symmetric models [47]. For the "right-sign" solution we find that the value ofα is strongly constrained (radian units): |α| ≤ 0.023 (95.5% probability).

Direct searches
The negative results from direct searches restrict the masses of the scalar particles. In order to access to the information that these observables provide, we first calculate the theoretical production cross section times branching ratio σ · B in the A2HDM. We consider then the ratio R ≡ (σ·B) theo /(σ·B) obs between the theoretical value and the observed limit, to which we assign a Gaussian likelihood with zero central value, which is in agreement with the null results obtained so far in the searches of heavy scalars. The corresponding standard deviation of the likelihood is adjusted in a way that the value R = 1 can be excluded with a probability of the 95%. The production cross sections and branching ratios for the other scalar particles are calculated in a similar way to the SM Higgs, taking into account the kinematically allowed region and the CP quantum number of the particle. In general, the data from direct searches favour heavier scalars and help us to restrict lower masses. However, since there are less experimental searches in the low-mass range, one gets less restrictive constraints for masses below 100 GeV. The constraints available so far seem to indicate that low masses are still allowed, so information from direct searches in that region would be crucial to understand the phenomenology at low masses.

Flavour constraints
Flavour observables are useful to constrain the Yukawa alignment parameters ς f . Figure 6 displays the allowed (95% probability) two-dimensional regions in the three-dimensional ς f space from independent analyses of the most relevant flavour measurements: B 0 mass mixing (dark pink), B s → µ + µ − (red), B → X s γ (blue) and (g − 2) µ (light pink). For clarity the observables that do not give relevant constraints in a given plane are omitted from the corresponding plot.
As already known from previous works [90][91][92], the (g − 2) µ anomaly requires sizeable NP contributions, which translates into non-zero values for ς that are rather large, while it is insensitive to ς d . This can be clearly seen in the central and right panels of figure 6 where the light-pink regions exclude values of |ς | below 10-20. The precise size of the discrepancy with the SM expectation relies, however, in a phenomenological evaluation of the hadronic contribution to the photon vacuum polarization (and a smaller light-by-light hadronic correction), involving a very subtle combination of different experimental data sets, which not always are in good agreement [93]. Thus, the statistical relevance of the (g − 2) µ anomaly could be magnified by underestimated uncertainties. The most recent lattice calculations seem in fact to suggest that the SM prediction could be much closer to the current (g − 2) µ measurement [94]. While waiting for a possible confirmation (or not) of this intriguing anomaly, and in order not to bias the results, we will not include (g − 2) µ in the global fit. We will discuss later whether the large values of |ς | currently required to accommodate (g − 2) µ are compatible with the parameter ranges emerging from a global fit to the other observables.
For similar reasons, the recent flavour anomalies observed in b → cτ ν and b → sµ + µ − transitions [95] will not be included either in our global fit. A recent model-independent analysis of the b → cτ ν anomaly has been already given in [96,97], where references to previous works can be found. If confirmed, these anomalies would provide clear evidence of NP with non-universal lepton couplings. The remaining flavour observables do not show significant deviations from the SM and, therefore, lead to allowed regions in figure 6 with smaller values of the alignment parameters, including the null SM solution. B s → µ + µ − is the only observable that constrains the leptonic couplings, excluding large values for |ς u,d ς | as indicated by the magenta areas in the figure. A similar restriction on the product |ς u ς d | can be appreciated in the left panel.

JHEP05(2021)005
The mass differences between the neutral B 0 eigenstates, ∆M B d and ∆M Bs , are dominated by virtual top-quark contributions. This results in a strong upper limit on |ς u |, which corresponds to the vertical pink bands in the left and central panels of figure 6. Similar limits emerge from ε K and R b [13]. However, it can be seen that values of |ς u | excluded at a 95.5% probability from ∆M B are no-longer excluded at the same probability when all flavour observables are combined. A stronger constraint on the ς u − ς d plane can be derived from the radiative decay b → sγ. The blue region in the left panel of the figure shows the strong correlation between the two quark alignment parameters, which prevents them to be large simultaneously.

Global fit
After discussing the separate effect of each type of observables, let us analyse the limits emerging from the global fit to all experimental and theoretical inputs.
The combined constraints on the scalar masses and mass differences are shown in the left panel of figure 2 and in figure 7. From these plots, it can be seen that theoretical and EWPOs constraints are complementary and by combining them with the remaining observables, light values for the masses are disfavoured. As for the electroweak constraints, there is a clear dependence on the mass priors. The global fits adopting squared mass priors with JHEP05(2021)005 However, the strong dependence on the mass priors indicates that these mass constraints should be taken with some care. The right panel in figure 2 shows the resulting allowed regions (68% probability in red and 95.5% probability in dark red) for the scalar potential parameters λ i when all constraints are included in the global fit. The addition of the Higgs signal strengths and the direct searches restricts significantly the parameter space obtained before from theoretical observables (blue areas). This effect is specially strong for λ 7 .
Combining the information from the Higgs signal strengths with the other observables turns out to be a bit subtle because the fine-tuned "wrong-sign" solutions discussed in section 4.3 lead to a very slow numerical convergence of the fit algorithm. To solve that, we have performed the fits shown in this section with the condition y h d, ≈ 1. The negative branch solution will be discussed separately in section 4.6.
For the positive branch, once we add the rest of observables to the Higgs signal strengths, the constraints of figure 5 get modified into the ones of figure 8. The global fit gives stronger limits for the alignment parameters, as expected, but leaves a somewhat wider allowed range for the mixing angle (radian units):   Figure 9 compares the constraints on the alignment parameters resulting from the combination of flavour observables (not including (g − 2) µ ) with the regions allowed by the global fit at 68% and 95.5% probability. The limits on the down-quark and lepton couplings become stronger, once all constraints are considered. This is mainly due to the combined effect of the Higgs and flavour observables. The strong correlation between the up-quark and down-quark alignment parameters observed before remains also in the global fit, with tighter upper bounds on |ς d |: larger values of the up coupling require smaller values of the down coupling and vice versa. A similar but weaker effect can be observed in the ς u − ς and ς d − ς planes.

"Wrong-sign" solution
The regions allowed by the global fit with the "wrong-sign" solution for the Higgs signal strengths are displayed in figure 10, at different probabilities. Since these constraints have been obtained imposing the "wrong-sign" solution, whose probability is smaller than 100%, the final probability would be (probability of the "wrong-sign" solution)×(probability of figure 10). The fine-tuned conditionας d, ∼ −2, emerging from y h d, ∼ −1, can be only satisfied in a small portion of the parameter space. The null value for the mixing angle is not reached, since it corresponds to y h f = 1 and, therefore, belongs to the normal "right-sign" branch discussed in the previous subsection.

Results: heavy scenario
In the previous section we have described the situation in which the observed Higgs corresponds to the lightest CP-even scalar of the model. In this section we will analyse the complementary situation, i.e. the heaviest CP-even scalar is the SM Higgs and there is an additional neutral scalar with mass below 125 GeV.
The theoretical constraints show the same tendency as for the light scenario. Since the mass of the CP-even scalar is now bounded to be light, the remaining two scalar masses cannot be heavier than 700 GeV. This can be seen in figure 11 for squared mass priors. Linear mass priors give very similar theoretical constraints, so they are omitted from the plot.  . Two-dimensional constraints on the scalar masses in the heavy scenario. The different allowed regions correspond to theoretical constraints (blue, squared priors, 100% probability), EWPOs (95.5% probability; orange, squared priors; violet, linear priors) and the combined global fit (red, squared priors, 95.5% probability).

JHEP05(2021)005
The constraints from EWPOs are also similar to the ones of the light scenario. Lower masses and large mass splittings are excluded when squared priors are adopted. Now M H is bounded to be smaller than 125 GeV, so the allowed regions become narrower. The worrisome difference between the results obtained with linear and squared mass priors is also present in this heavy scenario. Squared priors give very strong constraints for light masses that are no-longer found when linear priors are used. These constraints are displayed in figure 11, which shows that rather low values for all scalar masses are indeed allowed by the linear-priors fit. In the plane M A − M ± H the allowed regions for linear (violet) and squared (orange) priors overlap, so it not easy to distinguish them.
The analytical expressions of the gauge and Yukawa couplings of the light and heavy CP-even neutral scalars in eqs. (2.9) and (2.13) can be shifted with the change of variablẽ α =β − π 2 (notice that this bringsβ outside our previous convention forα), up to a global JHEP05(2021)005 Figure 12. Constraints from Higgs signal strengths in the heavy scenario. One can distinguish the "right-branch" withα ≈ − π 2 , corresponding to y h d, ≈ 1, and the "wrong branch" with (α + π 2 ) ς d, ≈ −2 for y h d, ≈ −1.
minus sign in the so-far unmeasured couplings of the additional neutral scalar. Therefore, the constraints onβ from the Higgs signal strengths would be similar to the ones obtained forα in the light scenario. Adopting the convention that g hV V should be positive, the hW W * and hZZ * measurements imply now thatα should be close to − π 2 . The "rightbranch" where y h d, have the same sign as y h u corresponds also to the region withα ≈ − π 2 , while the "wrong-branch" where y h d, have the opposite sign satisfies (α + π 2 ) ς f ≈ −2. As for the light scenario, the up Yukawa has always the same sign as g hV V . These allowed regions are displayed in figure 12. The sharp cut close toα = − π 2 in the negative branch is a consequence of the correlation betweenα and the down coupling ς d . Lower values of the mixing angle would require ς d < −50, which is not allowed by our priors.
Finally, most flavour constraints are independent of the neutral scalar masses, so they are identical in the light and heavy scenarios. The only relevant dependence appears in B s → µ + µ − for large values of ς d, [88]. The allowed regions in figure 6 remain then also valid in the heavy scenario, except the magenta areas that get slightly distorted.
The final constraints on the scalar masses from the global fit to this heavy scenario are displayed in figure 11, assuming squared mass priors (red regions). Large masses for the neutral CP-odd and charged scalars are not allowed at a 95.5% probability. Since the CP-even scalar is forced to have a small mass, electroweak constraints restrict the mass spitting between the two other scalars to be small. This is clearly seen in the M A − M H ± plane of figure 11.
Global fit results for theα − ς f planes in the positive branch are similar to the ones in the light scenario (see figure 8) shifting the mixing angleα →α − π 2 .

Summary
Several fits to the currently available data have been performed within the A2HDM, using the HEPfit tool that is based on Bayesian statistics. To reduce the number of fitted parameters, we have considered a CP-conserving scalar potential and real alignment couplings.

JHEP05(2021)005
We have included in the fit EWPOs, Higgs signal strengths, collider searches for additional scalars and flavour constraints, together with the theoretical requirements of perturbativity and vacuum stability. After analysing the separate implications of each type of constraints, we have combined all of them into a global fit to the model parameters.
Our main fits, discussed in section 4, assume that the observed Higgs at 125 GeV is the lightest CP-even scalar of the model. The theoretical requirements strongly constrain the mass differences among the three additional scalars, H, A and H ± , to be below 600 GeV. This upper bound becomes tighter as the masses increase, as shown in figure 2. The EWPOs further restrict the masses and mass differences, favouring small mass splittings and large masses. However, the precise bounds from EWPOs turn out to be very sensitive to the adopted priors. Taking squared mass priors with masses varied until 1500 GeV, one finds a lower bound around 750 GeV for the three masses, which gets reduced to 500 GeV if a lower mass range up to 1000 GeV is adopted as prior. Taking instead linear priors in the same mass range, masses as low as 10 GeV become allowed.
The agreement of the measured Higgs signal strengths with the SM expectations translates into a very strong constraint on the scalar mixing angle,α ≤ 0.023 rad (95.5% probability), and tight bounds on the alignment parameters, given in figure 5, which are further reinforced by the flavour constraints shown in figure 6. Combining all constraints into a global fit, one gets finally the allowed regions shown in figures 8 and 9. The up-quark alignment parameter must satisfy |ς u | < 1.5 (95.5% probability), while |ς d, | can take larger values provided the products |ς u ς d, | and |ς d ς | remain small. These figures assume that the down-quark and lepton Yukawa couplings do not deviate much from their SM values. However, since the current data on Higgs signal strengths cannot determine the signs of these two couplings, there is in addition a fine-tuned solution with "wrong-sign" Yukawas, shown in figure 10.
The global fit to all data does not solve the prior dependence of the fitted mass spectrum. As shown in figures 2 and 7, squared mass priors put stronger lower bounds on the scalar masses that depend on the assumed prior range, while linear priors still allow for quite low values of the masses. Clearly, the current negative results from collider searches are not yet stringent enough to discard the presence of new scalar states with masses near the electroweak scale.
We have also attempted a first study of the opposite scenario, where the 125 GeV Higgs is assumed to be the heaviest CP-even scalar. Using the same data set, we have found the constraints shown in figures 11 and 12 for the masses and alignment parameters, respectively. A very strong correlation between M A and M H ± is observed at high masses, as expected, because the mass splittings cannot become large. However, no useful lower bounds on the scalar masses can be extracted because they are again too sensitive to the adopted priors. A more detailed analysis of this scenario, including LEP searches and lowenergy data, could provide additional constraints that we plan to study in future works.
Our results are currently the most general global fit to the A2HDM. While previous phenomenological analyses of two-Higgs doublet models focused on particular cases based on discrete Z 2 symmetries or used only small subsets of observables, we have worked with a more generic theoretical framework with the only assumptions of flavour alignment and JHEP05(2021)005 real scalar and alignment parameters. A more general analysis, including the new sources of CP violation provided by the A2HDM, will be attempted in future publications.
Concerning the current flavour anomalies, it is worth to compare the fitted constraints on the alignment parameters in figure 9 with the parameter region able to accommodate (g − 2) µ , shown in figure 6. The global fit does not exclude the large values of |ς | which would be needed. However, in order to fit (g − 2) µ one also needs a quite light pseudoscalar with M A 50 GeV [90][91][92]. Specific searches for light scalar and pseudoscalar particles could be very relevant to investigate this possibility.
An explanation of the b → cτ ν and b → sµ + µ − anomalies within the context of two-Higgs doublet models would require non-universal lepton couplings. The generalised A2HDM [10], with family-dependent alignment parameters, provides a viable theoretical framework to address this type of phenomena. A scalar interpretation of the b → cτ ν data [98,99] seems still possible [96,97], but it would imply higher values of Br(B c → τ ν) than usually assumed. A detailed analysis of the b → sµ + µ − anomalies within the generalised A2HDM would provide very useful complementary information.

A Data compilation
The following tables detail the collider data sources employed in our global fit (data obtained at 2, 8 and 13 TeV are marked in green, purple and yellow, respectively). Table 4 compiles the LHC and Tevatron data sources on Higgs signal strengths. The information on heavy scalar searches at the LHC is collected in tables 5, 6, 7 and 8. These searches are applied either to the charged Higgs boson H ± or to the neutral scalars ϕ 0 i = H, A. Direct searches related to the charged Higgs boson are displayed in table 5.    [196] [0.05;1] 19.8 [197] [0.13;0.8] 36.1 [197] [0.13;0.8] 36.1