Updated global fit of the aligned two-Higgs-doublet model with heavy scalars

An updated global fit on the parameter-space of the Aligned Two-Higgs-Doublet model has been performed with the help of the open-source package HEPfit, assuming the Standard-Model Higgs to be the lightest scalar. No new sources of CP violation, other than the phase in the CKM matrix of the Standard Model, have been considered. A similar global fit was previously performed in Ref. [JHEP 05 (2021) 005] with a slightly different set of parameters. Our updated fit incorporates improved analyses of the theoretical constraints required for perturbative unitarity and boundedness of the scalar potential from below, additional flavour observables and updated data on direct searches of heavy scalars at the LHC, Higgs signal strengths and electroweak precision observables. Although not included in the main fit, the implications of the CDF measurement of the $W^\pm$ mass are also discussed.


Introduction
With the discovery of the Higgs boson [2,3], the Standard Model (SM) has become a well-established theory describing the interactions among elementary particles in a very elegant way.Nonetheless, there are ample evidences indicating that the SM cannot explain all aspects of nature, and hence the existence of physics beyond the Standard Model (BSM) is indispensable.In order to explain such BSM phenomena, the SM is usually enlarged with some additional particles or gauge groups that respect all its fundamental principles.While the augmentation of the SM with extra generations of quarks and leptons receives severe experimental constraints from unitarity-triangle data [4] or Zboson branching fractions [5], extensions with additional SU (2) L scalar doublets do not suffer from such stringent bounds. 1mong such scalar extensions, the simplest one is the two-Higgs-doublet model (THDM), where one additional scalar doublet with the same quantum numbers than the SM Higgs doublet is appended [7][8][9].In addition to the three needed Goldstone bosons, this model includes three neutral and one pair of charged scalars.Such rich scalar spectrum opens various interesting possibilities such as new sources of CP violation [10][11][12][13][14], axion-like phenomenology [15][16][17], dark matter aspects [18][19][20], neutrino mass generation [21,22], electroweak baryogenesis [23][24][25], stability of the scalar potential till the Planck scale [26][27][28], etc.Moreover, the THDM can also be thought of as a low-energy effective field theory framework for several models with larger symmetry groups (like supersymmetry).
One major shortcoming of the most general THDM is the emergence of tree-level flavour-changing neutral currents (FCNC), which are observed to be tightly constrained experimentally.This problem is usually avoided by imposing discrete Z 2 symmetries so that each type of right-handed fermions couples to one scalar doublet only [29,30].However, the absence of tree-level FCNC can be guaranteed with a much weaker requirement: the alignment of Yukawa couplings in the flavour space, so that the interactions of the two scalar doublets acquire the same flavour structure [31][32][33][34].This provides a more generic theoretical framework, known as the 'Aligned Two-Higgs Doublet model' (ATHDM), where flavour violation is minimal [35,36] and (highly-suppressed) FCNC only appear at higher perturbative orders [31][32][33][37][38][39][40][41][42][43][44].In this model, all the fermion-scalar interactions become proportional to the masses of the corresponding fermions, leading to a quite compelling phenomenology in high-energy colliders [44][45][46] as well as in low-energy flavour experiments [38,47,48].The ATHDM constitutes a generic platform for THDMs; all Z 2 -symmetric THDM scenarios can be explored as special cases of the ATHDM [31].In addition, the ATHDM can accommodate new sources of CP violation beyond the CKM matrix in both the scalar and Yukawa sectors [31].
The parameter space of the THDM has been extensively scrutinized in the literature, considering LHC data, LEP data, flavour bounds and theoretical constraints [1,38,[44][45][46].A global fit of the ATHDM with no new sources of CP violation (and a slightly different set of parameters) was performed in Ref. [1], taking into account various theoretical and experimental bounds.In this paper we reanalyse the global fit of this scenario in great detail, assuming that the SM Higgs boson is the lightest scalar.The fits are performed with the help of the open-source code HEPfit [77].Compared with the previous study, we have updated the Higgs signal strengths and added direct searches for all the scalars from the most recent LHC data.Additionally, we have improved the analysis of theoretical constraints by including necessary and sufficient conditions for boundedness of the scalar potential from below (i.e. the potential never tends to negative infinity), and incorporated the branching fractions for semileptonic decays of B, D and K mesons in the flavour sector too.We have also analysed the implications of the CDF measurement of the W ± mass [78], although we have not included it in our global fits in view of its currently unresolved discrepancy with other measurements.A similar conservative attitude has been taken for the muon g −2 anomaly [79], in view of the current controversy with lattice [80][81][82] and τ -decay [83,84] data; although not included in the global fit, we have also studied the ensuing constraints on the parameter space of the ATHDM.
The paper is organized in the following way.Next section (Section 2) briefly describes the ATHDM scenario.The set up for the global fit is discussed in Section 3.While in Section 4, we illustrate different theoretical as well as experimental constraints that are considered for this study, in the subsequent section (Section 5) we present the numerical results from the fits.Finally we conclude in Section 6.All the data and various other details used in the global fit are incorporated in the Appendix section.

Scalar Sector
We extend the SM with a second complex scalar doublet having the same hypercharge as the SM scalar doublet, i.e.Y = 1/2.The neutral components of both scalar doublets may acquire non-zero complex vacuum expectation values (vev); nonetheless, with a suitable SU (2) L ⊗U (1) Y global transformation, one can always rotate the basis of the scalar space such that only the neutral component of the first doublet acquires a non-zero (real) vev.Working in this so-called Higgs-basis, we can write down the two doublet scalar fields as: where Φ 1 gets the vev v = 246 GeV.The components G ± and G 0 act as Goldstone bosons, providing the masses to the W ± and Z bosons.Thus, we are left with one pair of charged scalars H ± , two neutral scalars S 1,2 and one neutral pseudoscalar S 3 .The three neutral particles S j mix with each other through an orthogonal transformation to produce the mass eigenstates φ 0 i .The explicit form of this orthogonal transformation depends on the scalar potential.In general, if the scalar potential is not invariant under the CP symmetry, all the three S j fields mix together giving no definite CP to the mass eigenstates.
The most general scalar potential, allowed by the SM gauge symmetry, takes the form: where µ 3 , λ 5 , λ 6 and λ 7 are complex parameters and the rest are real.Minimizing the potential with respect to the neutral fields at their corresponding minima, one obtains: implying that µ 1 and µ 3 are not independent quantities.Redefining the phase of Φ 2 , one can make any one of the three parameters λ 5 , λ 6 or λ 7 real.Therefore, the scalar potential involves eleven independent real parameters: µ 2 , v, λ 1,2,3,4 , |λ 5,6,7 | and the two relative phases between λ 5,6,7 .
The quadratic terms of the potential, which give rise to the scalar masses, can be written as: with In order to obtain the physical mass eigenstates, we need to diagonalise the neutral mass matrix (2.5).Imposing CP conservation in the scalar sector, all the parameters in our potential become real, reducing to nine the number of independent inputs.Moreover, in this case the mass eigenstates have definite CP: we get two CP-even (h, H) and one CP-odd (A) fields.From the mass matrix, we explicitly see that in this case S 3 = A does not mix with the other scalars and we only need to diagonalise a 2 × 2 matrix.The squared masses are then given by: Note also that, since the trace of the mass matrix must be invariant, the scalar masses must satisfy the following relation: which can also be easily confirmed by substituting the expressions of the masses.The mixing between the two CP-even neutral scalars, is given by: (2.10) Using the above equations, we can trade five of the nine parameters of the potential (µ 2 , v and λ 1−7 ) by the four scalar masses (M H ± , M h , M H , M A ) and the mixing angle α, which are more physical.Clearly λ 2 must be kept since it does not relate with any of the new parameters.Since µ 2 and λ 3 always appear in the same combination (M H ± ), we must keep one of them which we choose to be λ 3 .From Eq. (2.10) we can relate λ 1 and λ 6 with α, M 2 h , M 2 H and v. Then we can use Eq.(2.8) to obtain λ 4 and Σ to obtain λ 5 .With this procedure we trade µ 2 , λ 1 , λ 4 , λ 5 and λ 6 by the four masses and the mixing angle.Our choice for the nine input parameters is then: v, M H ± , M h , M H , M A , α, λ 2 , λ 3 and λ 7 .The remaining parameters of the scalar potential can be easily obtained in terms of these inputs: While considering interactions of neutral scalars with gauge bosons, S 1 plays the role of the SM Higgs boson.This implies that the couplings of the scalar mass eigenstates with a pair of gauge bosons are given by: where V V ≡ W + W − , ZZ.

Yukawa Sector
The interactions of the fermion mass eigenstates with the scalar fields read: where the generation indices are suppressed and the subscripts L, R denote the usual left and right chiral fields.Here, M f (f ≡ u, d, l) are diagonal mass matrices for the up-type quark, down-type quarks and charged leptons, respectively, originated through the Yukawa interactions with the doublet Φ 1 .Y f are the Yukawa matrices parametrizing the fermionic couplings with the second doublet Φ 2 , and V is the usual CKM matrix required for the quark mixing.In general, Y f could be arbitrary 3 × 3 complex matrices, which leads to unwanted FCNCs at tree level.This can be easily avoided imposing the alignment condition of M f and Y f in flavour space [31,32], i.e.
where the Yukawa alignment parameters ς f could be arbitrary complex numbers.In terms of the scalar and fermion mass eigenstates, the Yukawa Lagrangian takes then the form: (2.15) where P L,R are chirality projection operators and φ 0 i are the scalar mass eigenstates.In the following, we will assume a CP-conserving potential and that there are no additional sources of CP violation beyond the CKM phase, i.e. we will only consider real alignment parameters.The Yukawa couplings of the neutral scalars are then given by: ) (2.17) The usual THDM scenarios based on Z 2 symmetries are just particular cases of the more general ATHDM framework; they can be retrieved by imposing µ 3 = λ 6 = λ 7 = 0 along with the following conditions: The alignment requirement (2.14) remains stable under renormalisation when it is protected by Z 2 symmetries [37].Otherwise, higher-order quantum corrections create a misalignment of M f and Y f that generates loop-suppressed FCNC effects.However, the special Yukawa structure of the ATHDM strongly constrain the possible FCNC interactions, making those effects numerically small [31,32,38].Assuming exact alignment at some high-energy scale (even at the Planck mass), the small misalignment generated by running down to low energies remains well below the current experimental limits [33,39,40,43].

Fit set up
Our numerical analyses have been performed with the open-source HEPfit package [77].This code has been widely used, due to its efficiency and versatility that allows making global fits, both within the SM [85] or in general BSM extensions such as the SM effective field theory [86,87] or particular models of new physics [1,88], as it is the case of this paper.HEPfit works within the Bayesian statistics framework and, therefore, we need to choose carefully the priors of our variables.We have a total of ten new degrees of freedom with respect to the SM: the physical masses of the additional scalars (M H ± , M H and M A ), the mixing angle of the CP-even neutral scalars (α), three quartic couplings of the potential (λ 2 , λ 3 and λ 7 ) and the three Yukawa alignment parameters (ς u , ς d and ς l ).Since our fits include all the available information, we have used uniform distributions as priors.
In general, the priors are chosen to cover the region of the parameters that is physically relevant.For the mixing angle we have chosen the prior in such a way that at least the 5σ region of the posterior probability is contained within the selected range.The quartic couplings are mainly constrained from theory assumptions, which impose a hard cut in the values that these parameters can take.In this case we have chosen a prior wide enough to include all points allowed by the theory constraints.The priors of the Yukawa couplings were also set within the limits allowed by perturbative constraints.
There is a huge freedom when choosing ranges for the scalar masses.In this analysis we are assuming that the SM Higgs is the lightest scalar, i.e. that M H , M A and M H ± are larger than 125 GeV.The complementary scenario where the SM Higgs is not the lightest scalar is phenomenologically very intriguing and capable of displaying very distinct collider signatures at the LHC [58][59][60].However, these signatures and the overall phenomenology are strongly dependent on the assumed hierarchy of scalar masses (which BSM particles are considered to be lighter than the SM Higgs).Moreover, a detailed analysis of the different possible scenarios requires the inclusion of additional experimental data from LEP and flavour factories, and a proper theoretical treatment of hadronic resonances in the mass region below 10 GeV.These intricacies for the light BSM particles are beyond the scope of this paper and therefore, here we focus only on the heavy BSM scenario. 2igher scalar masses are obviously preferred by the data because no clear deviations from the SM have been observed so far.Taking into account that many direct searches have been studying mass ranges up to 1 TeV, this number seems a reasonable higher cut-off for our global analysis.However, we will also provide some results in which the highest value of the scalar masses have been set to 1.5 TeV, so that we can get a feeling on how much our results depend on the priors.Finally, we would also like to comment that, instead of taking the masses as fundamental parameters, using the masses squared could also be a well-justified choice.However, as can be seen in Ref. [1], when using a uniform distribution for the masses squared the results depend much more on the ranges of masses analysed.Therefore, we have decided to impose our priors linearly on the scalar mass parameters.The chosen priors can be found in Table 1. Priors Table 1.Priors chosen for the new-physics parameters.

Fit constraints
We aim to use as much information as we can in order to constrain the parameter space of the ATHDM.We have combined in this work the whole set of theoretical constraints with all relevant experimental limits coming from LHC direct searches, Higgs data, electroweak precision data, and flavour observables.

Theoretical considerations
Regarding theoretical constraints, the two conditions that are usually demanded are: a scalar potential bounded from below and perturbative unitarity of the S matrix.The requirement that the scalar potential is bounded from below indicates that it should not go to large negative values, which would make it unstable, for any configuration of the fields.For this purpose, it is useful to recast the scalar potential V in the following Minkowskian form [89]: with and Diagonalisation of the mixed-symmetric matrix3 Λ µ ν produces one "timelike" (Λ 0 ) and three "spacelike" (Λ 1,2,3 ) eigenvalues. 4The scalar potential remains bounded from below when the following two conditions are satisfied [89,90]: Several necessary conditions for the scalar potential to be bounded from below have been listed in Ref. [91].However, they are comprehended in the two above-mentioned necessary and sufficient conditions.
One additional constraint could be imposed on the quartic couplings by demanding that the vacuum of the scalar potential is a stable neutral minimum. 5For this purpose one defines the discriminant of the matrix ξ I 4 − Λ µ ν as: where the Lagrange multiplier ξ is determined by minimizing the scalar potential V with the constraint6 r µ r µ = 0.The condition for a global minimum is found to be [90]: D > 0, or D < 0 with ξ > Λ 0 .In our analysis we have included this condition, requiring a stable neutral minimum.However, removing this constraint does not change our results significantly.
The unitarity of the S matrix ensures that, in order to conserve the total probability, scattering amplitudes do not grow monotonically with energy.Since unitarity emerges from the basic formulation of QFT, it is bound to hold in a complete renormalisable theory while dealing with the full S matrix.However, the unitarity bound does not need to be satisfied by the perturbative calculation of the S matrix to any particular order.The stronger requirement of perturbative unitarity enforces the unitarity of the S matrix to be obeyed at every order of perturbation theory.It indirectly indicates that the couplings, in terms of which the perturbative expansion is performed, are not very large and one can safely neglect higher-order contributions.Therefore, imposition of perturbative unitarity on all the 2 → 2 scattering amplitudes involving the scalars (both the massive ones and the Goldstone bosons, which account for the W ± L and Z L ) will restrict the quartic couplings λ i , ensuring that the perturbative expansion of the S matrix does not diverge at high energies.Tree-level unitarity is enforced on the 2 → 2 scattering matrix of scalars by demanding that: where a 0 is the matrix of tree-level partial-wave amplitudes and a 0 j the corresponding eigenvalues in the j th partial wave.Nevertheless, while considering the scattering of scalars at very high energy, only the S-wave amplitude with j = 0 becomes most significant at tree level.Constructing two-scalar scattering states with definite hypercharge and weak isospin (Y, I) and grouping the ones with the same set of quantum numbers, the 25-dimensional7 matrix a 0 can be expressed in a block-diagonal form: where the superscript indicates the total charge of the initial or final states.Thus the scattering matrix a 0 is comprised of the following four submatrices [91,94]: The perturbative-unitarity condition (4.5) can be traded in terms of the eigenvalues (e i ) of the above four submatrices by demanding: On the other hand, considering the coupling of fermions to the charged scalars, we vary the value of ς f in the perturbative range satisfying A more detailed study on the perturbativity of the alignment parameters could be performed based on the procedure mentioned in Ref. [95], which would extend the range of ς f a bit more.However, for our analysis it is sufficient to consider the above inequality.

Direct searches
Many searches of additional scalars have been developed at the LHC.In order to use those results we have compared the theoretical prediction of the cross-section times branching ratio, σ • B, for several processes with the exclusion limits obtained by the CMS and ATLAS experimental collaborations.
The experimental data are provided in the form of tables, which compile the values of the 95% upper limits on σ • B, as a function of the resonance mass.In these tables, HEPfit performs a linear interpolation if needed.In order to compare the theoretical results with the experimental data, we define the ratio of the theoretical prediction over the experimental upper limit.To this ratio we assign a Gaussian distribution (restricted to positive values) centered at 0 such that the value 1 is excluded with a 95% probability.
All the channels included in our fit can be found in Appendix A. In Tab. 4 we show the channels included for the charged scalars, while those of the neutral scalars are shown in Tabs.5, 6 and 7. Specifically, in Tab. 5 are shown the decays into neutral scalars (including the SM Higgs boson); in Tab.6 the decays into fermions, γγ and Zγ; and in Tab.7 the decays into weak gauge bosons.When the limits on σ • B provided in the experimental papers consider also a subsequent decay of the SM particles, we show this final decay with parenthesis.When the limits on σ • B are provided directly on the decay width of the NP particles to SM particles, but using particular decays of the produced SM particles, we show such SM decays with square brackets.

Higgs data
The presence of additional scalars generates relevant effects on the production and decay of the SM Higgs.The one-loop h → 2γ amplitude receives additional contributions from the charged scalar.Furthermore, due to the mixing among the CP-even scalars, the coupling of the SM Higgs with the weak bosons is modified as shown in Eq. (2.12), which changes the Higgs production through vector boson fusion (VBF) and the associated production with vector bosons (Vh).The decays of the Higgs boson to fermions are also sensitive to the scalar mixing, as shown in Eqs.(2.16) and (2.17), which also modifies the Higgs production through gluon fusion (ggF) and its associated production with t t pairs (tth).
The production and the subsequent decay of the SM Higgs boson have been measured (or bounds have been set) at the LHC for the most relevant production modes (ggF, VBF, Vh and tth) and decay channels (cc, b b, γγ, µ + µ − , τ + τ − , W W , Zγ and ZZ).These data are parametrised in terms of the Higgs signal strengths, which are defined as the measured cross section times branching ratio for a given production and decay Higgs channel, in units of the SM prediction.Table 3, in Appendix A, compiles the experimental papers from which we have taken the values of the different Higgs signal strengths relevant to our analysis.

Electroweak precision observables
The presence of additional scalars generates contributions to the oblique parameters (also known as Peskin-Takeuchi parameters [96,97]) S, T and U [98].The experimental values of these parameters are obtained from global fits of electroweak precision data, using observables directly measured (mainly) at LEP and SLC.Among those observables, we highlight the ratio R b ≡ Γ(Z → b b)/Γ(Z → hadrons) [99,100], which is also affected by the additional scalars.Using the values of the oblique parameters from the PDG would be inconsistent in our study because those values do not take into account the NP contamination on R b .Moreover, since we include R b directly in our global fit of the ATHDM, the information from this observable cannot be employed to determine the oblique parameters.
In order to obtain non-contaminated values for the oblique parameters, we have repeated the electroweak fit removing R b from the list of fitted measurements, obtaining the values of S, T and U that will be included as inputs to our analysis.Those values are summarised in Tab. 9 of Appendix B.2.
The values of the oblique parameters are also highly dependent on the value of the W -boson mass (M W ). For our baseline results on these parameters we have used the value of M W quoted by the PDG 2022 [101].However, in April of 2022 the CDF collaboration announced a controversial new measurement of M W [78], which is incompatible with the SM prediction by 7σ.Since there is not yet consensus in the community, we have not included this measurement in our main global fit but, instead, we will show how much the results change with the values for the oblique parameters obtained incorporating this new measurement in a global electroweak fit performed with HEPfit [102,103]. 8inally, we have also noticed that using the three oblique parameters S, T and U or fitting only S and T (assuming U to be negligible) gives very similar results.Therefore, we use only S and T , since it is well known that the THDM contributions to U are highly suppressed [98].

Flavour observables
The parameter space of the ATHDM can also be constrained from flavour observables, since the new scalars generate relevant contributions to many of them.However, in order to be able to determine these contributions, we first need to know the numerical values of the CKM parameters.For this work we have adopted the Wolfenstein parametrisation [104], so we need to provide the corresponding values as inputs in our analysis.The world averages quoted in the PDG 2022 [101] originate in SM fits from the CKMfitter [105] and UTfit [4,106] collaborations, which make use of several flavour transitions that could be affected by the additional scalars.The UTfit collaboration also provides the values of the Wolfenstein parameters removing the loop processes from the fit [107], which gives the correct values for models in which the tree-level effects of NP are negligible.Nevertheless, for some (small) regions of the parameter space in the ATHDM there could still be some contamination to some of the tree-level processes used to determine the Wolfenstein parameters.For this reason we have decided to repeat this CKM fit, using only processes that are not contaminated by the additional scalars.Details on this fit, as well as the resulting numerical values for the CKM parameters, can be found in Appendix B.1.Note also that we use as observables for our ATHDM global fit some processes that are usually taken into account in the CKM fit, like the pseudoscalar-meson leptonic decays, so it is clear that we should remove these processes from our CKM determination.
We have not included in our main fit the anomalous magnetic moment of the muon, (g − 2) µ [117], because there is no full consensus in the community on the SM prediction of this observable.The most recent lattice computations of the hadronic vacuum polarisation [80][81][82] and the estimates from τ decay [83,84] seem to find a SM result much closer to the experimental value [118,119] than the dispersive e + e − prediction [79], and additional hints of a possible dispersive underestimate are suggested by a recent QCD analysis of the Adler function [120] and the CMD3 data [121].Nevertheless, this observable is included in our code and we will show the individual constraints obtained with it, comparing those results with the ones obtained using all the other flavour observables that we have mentioned here.

Fit results
Here we provide the results for all the fits we have performed in this analysis.We discuss first in several subsections the limits obtained using only some subsets of observables, in order to show the relevance that each observable has in the global fit.Afterwards, we present the global fit, which provides the main results of this work.We show also how the results of our global fit would change if the new CDF measurement of M W [78] is included.

Theoretical bounds
The theoretical considerations generate constraints on the parameters of the scalar potential.As shown in Section 2, we have a total of nine degrees of freedom in the CP-conserving potential.The scalar vacuum expectation value and M h are fixed by the measurement of the Fermi constant in muon decay [101] and the Higgs mass measurement at the LHC [101].We then have a total of seven degrees of freedom which could be possibly constrained with the theory assumptions.As mentioned before, we trade some of these scalar potential parameters with more physical inputs.Our chosen parameters are the masses of the new scalars, the mixing angle among the CP-even scalars, and the three quartic couplings λ 2 , λ 3 and λ 7 .
The theoretical constraints on the parameters of the potential can be translated into limits on the scalar mass splittings.Fig. 1 displays the resulting constraints on the correlation among the different masses and the correlation of the charged scalar mass with the mixing angle.Note that in this case we are showing the 100% probability region, since all the points that satisfy these constraints are equally valid.In general, we observe that for scalar mass values below 700 GeV, the constraints on the mass splittings and the mixing angle are rather weak.This can be easily understood, since for low masses the mass differences cannot be large enough to reach the theoretical bounds.However, above this threshold, especially beyond 1 TeV, the constraints become considerably more stringent.Since differences of masses squared are proportional to combinations of quartic couplings times the vev squared, the naive estimate √ 4πv ≈ 870 GeV gives indeed a good idea of the scale where these limits become relevant.The constraints on the quartic couplings are shown in Fig. 2, where we can see that negative values of λ 2 are forbidden while slightly negative values of λ 3 could be allowed, provided they are very close to zero.Finally, λ 7 is clearly constrained to be smaller than 3 in modulus.

Experimental bounds
Fig. 3 shows that the measured Higgs signal strengths generate tight constraints in the α − ς f planes.At 99.7% probability some allowed regions, at the corners of the central and right plots of Fig. 3, appear really far away from the SM solution.These regions correspond to down-type and/or lepton Yukawas of opposite sign to the SM Higgs couplings.Since these observables are sensitive only to the modulus of the Yukawa couplings, this could be a possible solution, although following a bayesian approach the prior probability of this kind of solutions would be quite small.No such regions appear for up-type quarks because the relative sign between the top Yukawa and the Higgs coupling to the gauge bosons gets constrained by the h → 2γ decay width and the range of allowed |ς u | values is much smaller.
The oblique parameters (S and T ) generate significant constraints on the mass splitting of the additional scalars.These constraints can be found in Fig. 4, where we compare the results obtained when S and T are fitted with the PDG 2022 value of M W (PDG-M W ) with those using the modified S and T values emerging from the world average [103] after including the CDF M W (CDF-M W ) measurement [78].The allowed regions using the PDG value of M W are compatible with small or even zero mass splittings.However, once the CDF value is included in the fit, a non-zero mass difference between the new scalars is definitely required to explain the deviation from the SM prediction of M W .This scalar mass difference could be accommodated within this model, since we are able to generate the needed contribution to M W , even when we impose agreement with all other observables in combination, as we will discuss in Section 5.4.
Tree-level combined 95% All Flavour 95% All Flavour 95% g 2 95% ς d , ς l ) and the mass of the charged Higgs M H ± .Fig. 5 displays the correlations among ς u , ς d and ς l , while in Fig. 6 we show the correlation of those parameters with M H ± .The dominant NP contribution to the B → X s γ amplitude is proportional to the product ς u ς d , which explains the strong correlation between these two alignment parameters observed in the left panel of Fig. 5.The first two panels in Fig. 6 exhibit also the presence of two additional bands emerging from the central B → X s γ allowed region; they correspond to solutions with a NP contribution equal to minus two times the SM amplitude.The strongest constraint on the ς u − M H ± plane comes from meson mixing, although R b , B → X s γ and B s → µ + µ − also help in constraining this plane.In the ς d − M H ± plane the tree-level leptonic and semileptonic decays of pseudoscalar mesons (and tau decays into pseudoscalars) generate the strongest limits for low values of the scalar masses.However, the treelevel NP contributions rapidly drop when increasing the mass of the charged scalar, so for larger mass values the loop processes dominate.B → X s γ provides the strongest limits for intermediate values of M H ± around 500 GeV, while for heavier masses B s → µ + µ − produces the most relevant constraints because it also gets contributions from neutral scalar exchanges that are sizeable at large |ς d |.In the ς l − M H ± plane the constraints are very poor, except at very low values of M H ± where useful limits are obtained from the tree-level processes.As shown in the first two panels of Fig. 6, small values of M H ± are allowed when ς u,d → 0, but when marginalising over these two alignment parameters their probability becomes very small, which explains the lower bound on the charged scalar mass observed in the third panel.Notice that the maximum and minimum values allowed for the different parameters are slightly different in Figs. 5 and 6; this is just a consequence of having non-Gaussian distributions.Finally, in Fig. 7 we compare the constraints on ς l obtained from the combination of all flavour observables with the results needed to explain the anomalous magnetic moment of the muon.Since we are only considering scalar masses above 125 GeV, quite large values of the leptonic alignment parameter are required, but there is, indeed, room to satisfy all the flavour constraints in combination with (g − 2) µ .Nevertheless, we decided to exclude this observable from our global fit because its SM prediction is currently under debate, as discussed in Section 4.5.In this section we provide the results emerging from the global fit, using all the experimental observables and theoretical considerations at the same time.As mentioned in Section 3, we have considered masses up to 1 and 1.5 TeV, in such a way that we can also explicitly see the dependence on the adopted priors.The marginalised probabilities, given in Tab. 2, exhibit indeed some dependence on the priors, specially regarding the mass limits.This behaviour is reasonable since the likelihood is maximised for higher values of the masses; smaller masses would always be disfavoured, moving the 95% region further away when we allow for higher mass values.The allowed ranges of the Yukawa alignment parameters are also wider when heavier masses are allowed because the scalar contributions to the flavour observables decrease with increasing mass values and, therefore, higher values of ς u,d,l become possible, as can be observed in Fig. 6.In contrast, the preferred region of the mixing angle shrinks when we scan over heavier masses because the theoretical constraints, shown in Fig. 1, become more severe for higher values of the scalar masses.

Masses up to 1 TeV
The correlations among the scalar masses are shown in Fig. 8, which exhibits diagonal bands enforced by the theoretical constraints and the oblique parameters.The bounds on the masses get significantly relaxed, getting close to disappear, in the limit of zero mass splittings between the charged and the neutral scalars.This is due to the reduction on the constraints from the oblique parameters, which makes these points much more favoured.The correlations among the mass splittings are shown in Fig. 9.In the middle plot one observes that the two mass differences between the charged scalar and the neutral CP-even and CP-odd bosons cannot be simultaneously large; when one increases the other must decrease.From the left and right plots we can also see that when the mass splitting among the neutral scalars gets large the mass splitting among the charged scalar and one of the neutral scalars must go to zero.Here again, this behaviour is generated by the oblique parameters, whose NP contribution tends to zero when the mass splitting between any neutral scalar and the charged scalar vanishes.The last interesting plots regarding the masses are their correlations with the Yukawa alignment parameters and the mixing angle, shown in Fig. 10.Here we can see how the higher the masses get, the higher the Yukawa alignment parameters can be.Nevertheless, for masses above 700 GeV, smaller values of the mixing angle are preferred due to the theoretical constraints, as explained before.Finally, for very small values of the Yukawa alignment parameters or the mixing angle the limits on the masses also start to disappear, as expected.
The correlations of the Yukawa alignment parameters among themselves and with the mixing angle are shown in Figs.11 and 12, respectively.In the first one we observe that the regions get wider when higher masses are allowed, since all experimental constraints become softer, as also shown in Fig. 10.Moreover, once one of the Yukawa alignment parameters tends to high values, the others get highly suppressed.Indeed, for masses up to 1 TeV, in order to obtain values of ς d or ς l of order 10, ς u must be of order 0.1, or smaller.We observe a similar behaviour looking at the correlations of the Yukawa alignment parameters with the mixing angle.The higher the mixing angle, the smaller the Yukawa alignment parameters must be; and the higher the allowed range of masses, the higher the Yukawa alignment parameters can be.However, as can also be seen in Fig. 10, the mixing angle does not grow with the masses, and slightly smaller values of α are preferred when heavier masses are allowed in the fit.Nevertheless, performing a fit with a wider mass region does not significantly change the shape of the correlations between the Yukawa alignment parameters and the mixing angle with respect to the masses of the scalars.Integrating over the masses, like in Figs.11 and 12, the shape stays basically unchanged and just rescales when a different mass range is covered.Indeed, performing the fits for these two different mass ranges allows us to get an idea of the shape of the multidimensional correlation of the masses against the other parameters which, unfortunately, cannot be printed in a two dimensional plot.
In Fig. 11 we also show the restrictions (Eq.(2.18)) that each Z 2 THDM would impose on the Yukawa alignment couplings.In these models all alignment parameters are related through an angle β, whose tangent measures the ratio of the vevs of the two scalar fields in the basis where the discrete Z 2 symmetry is imposed.Therefore, each Z 2 -symmetric model corresponds to specific curves (either straight lines or hyperbolas) in the ς f 1 − ς f 2 planes, each point on which represents different values of tan β.The curves in green, cyan, orange and black depict the Type-I, Type-II, Type-X and Type-Y THDMs, respectively, whereas the inert THDM corresponds to an isolated point at the origin.The first panel in Fig. 11 indicates that the type-II and type-Y models have some tension with the (mainly flavour) data.Note, however, that this is just a qualitative comparison.In order to analyse the parameter space of any of these particular models, we would need to change the priors of our analysis by considering tan β as the free parameter instead of the alignment parameters along with the constraint µ 3 = λ 6 = λ 7 = 0.A comprehensive investigation of the parameter space for each of these models, though valuable, falls outside the scope of this study.We refer the interested reader to the more recent specific analyses in Refs.[50,51].Finally, in Fig. 13 we plot the correlations among the quartic couplings of the scalar potential.These couplings are mainly constrained by the theoretical bounds (perturbative unitarity, boundness from below and vacuum stability) and therefore their constraints do not depend on the mass range studied.Indeed, these limits are approximately the same as those of Fig. 2, in which only the theoretical considerations have been taken into account.

CDF W mass measurement
As mentioned before, last year the CDF collaboration released a new measurement of the mass of the W boson [78] with a 7σ tension with the SM prediction.Since there is still a controversy in the community, we decided not to include this last measurement in our main analysis but we provide here the results of a global fit including that information.In particular, we take the values of the oblique parameters (S and T ) from Ref. [103].Using that data as an input, we obtain a posterior value for M W = 80.4178 ± 0.0091 GeV which is compatible with the CDF measurement (M CDF W = 80.433 ± 0.009) within a 95% probability.The additional scalars of the ATHDM would then be able to explain this result, if confirmed. 9As can be seen in Figs. 14 and 15, in order to generate this additional contribution to M W we need a mass difference among the charged scalar and the neutral scalars of, at least, few tens of GeVs.Besides this mass difference, the results for the other parameters are very similar to the ones obtained in the baseline fit, as expected.

Conclusion
In this work we have performed an extensive phenomenological analysis of the CP-conserving Aligned Two-Higgs-Doublet model, using the HEPfit package.These results update the previous work of Ref. [1],10 although we use here a slightly different set of parameters.The observed SM Higgs boson at a mass of 125 GeV has been assumed to be the lightest scalar particle and we have, therefore, focused our study in heavy new-physics scalars.More details on the region of the parameter space that we have scanned can be found in Section 3.
We have incorporated the theoretical requirements of perturbative unitarity, boundedness of the scalar potential from below and vacuum stability, and the most relevant experimental constraints, including direct searches at the LHC, Higgs signal strengths, electroweak precision observables and flavour data.The main results of the global fit, combining all theoretical and experimental constraints at the same time, are shown in Section 5.3.
We have a total of ten new-physics parameters to fit.Three of them are the Yukawa alignment parameters for up-type quarks (ς u ), down-type quarks (ς d ) and leptons (ς l ).The other seven come from the scalar potential, where we have chosen the three heavy scalar masses (M H , M A and M H ± ), the mixing angle among the neutral scalars ( α), and three quartic couplings (λ 2 , λ 3 and λ 7 ).The only relevant constraints on these quartic couplings come from theoretical considerations.The mixing angle is mainly constrained by the Higgs signal strengths and lower limits on the masses of the new scalars are obtained from direct searches and flavour observables (which constrain the mass of the charged scalar).The mass difference among the charged scalar and (at least one of) the neutral scalars is highly constrained by the oblique parameters.The theoretical requirements also provide relevant constraints on the mass splittings and on the mixing angle for high values of the charged scalar mass.Finally, the Yukawa alignment parameters are mainly constrained by the direct searches and flavour observables.
The marginalised probabilities for all the parameters can be found in Tab. 2. In this table we show results for a fit in which the scan on the masses is done up to 1 TeV and one in which the scan is done up to 1.5 TeV.These ranges of masses have been chosen according to the ranges testable at the LHC.Obviously, the results differ according to the range considered, since higher masses are always favoured.Note also that the lower bounds on the masses are highly correlated with the values of the Yukawa alignment parameters.If we set these parameters to sufficiently low values, the constraints on the masses vanish.For masses up to 1 TeV, values of |ς u | > 0.7 , |ς d | > 12, |ς l | > 30 and |α| > 0.06 lay outside the 95% probability region.
Besides performing this global fit with the state of the art data, we have also studied two additional measurements that, at the moment, deviate from the SM prediction: the muon (g − 2) [118,119] and the new measurement of M W from the CDF collaboration [78].In both cases we have not included them in the global fit because they are still controversial, although for different reasons.Nevertheless, we have shown that the ATHDM has enough flexibility to accommodate both measurements within its parameter space, while being still compatible with all other constraints.

A Data Included
The following tables compile all experimental references included in the fit for the direct searches and the Higgs signal strengths.

Mass range L
Table 7. Direct searches for neutral heavy scalars, φ 0 i = H, A, with vector-boson final states.V = W, Z, ℓ = e, µ.The parenthesis show the final decay of the SM particles produced from the NP particles.The square brackets are used when the values of σ • B are shown in terms of the primary decay (i.e. the NP particle decay) but a particular decay channel of the SM particles is used to obtain those values.See Section 4.2 for more details.

B Non-contaminated Standard Model inputs
In this work we have included as inputs the entries of the CKM matrix, in the Wolfenstein parametrisation [104], and the oblique parameters [96,97].The experimental values of these parameters quoted by the Particle Data Group [101] have been obtained through a SM fit of multiple observables, neglecting any new-physics contributions.Unfortunately, this assumption is not satisfied in our model for all the observables included in the fits, and we need to repeat those fits removing the problematic observables.

B.1 Non-contaminated CKM fit
We extract the Wolfenstein parameters from the measured values of the CKM matrix elements (or ratios among them).

B.1.1 Determination of |V ud |
For |V ud | we use directly the value quoted by the PDG [101], since this value comes from superallowed 0 + → 0 + nuclear β decays [237] with completely negligible contamination from our additional scalars.

B.1.2 Determination of |V us |
The determination of |V us | in the PDG [101] is based on semileptonic decays of kaons, averaging the electronic and muonic channel, and, independently, on the pure leptonic decay of kaons and pions to muons which provide the ratio |V us /V ud |.The scalar contribution to the semileptonic decays is highly suppressed by the light lepton masses.Therefore, as long as the mass of the kaons is much higher than the mass of the decaying lepton we can safely neglect the NP contribution in the semileptonic decays.This assumption holds for the semileptonic decay to electrons but it is not fulfilled for the muonic channel.Therefore, for the semileptonic decays we only consider the decays to electrons, including the K L , K S and K + decays, which give the average value |V us f K + (0)| = 0.21626 ± 0.00040 [238].Taking the form factor average. f + (0) = 0.9698 ± 0.0017, from N f = 2 + 1 + 1 lattice QCD calculations [239], gives |V us | = 0.2230 ± 0.0006.
As mentioned above, the alternative determination of |V us | is based on the leptonic kaon decay.In this case, the SM contribution is helicity suppressed, making the scalar contribution relatively much higher.In general, the charged-scalar contribution to the leptonic decay of a pseudoscalar-meson takes the form [38] Γ(P with i, j the flavour indices corresponding to the valence quarks of the meson, the new-physics correction where f P is the meson decay constant and δ M ℓ2 em the electromagnetic radiative corrections.In particular, in order to determine |V us |, the ratio among the kaon and pion decay widths into muons is used, which gets a scalar contribution that is dominated by 2∆ us ≈ 2ς * l ς d m 2 K /M H ± .Since ς * l and ς d can, in general, reach quite high values, we cannot neglect the scalar contribution in this case and, therefore, we cannot utilize this for the |V us | determination.Note that the determination of |V us | obtained from the average of measurements involving leptonic decays (|V us | = 0.2252 ± 0.0005) is not compatible by more than two standard deviations with the determination derived from the average of measurements involving semileptonic decays (|V us | = 0.2231 ± 0.0006).When these two determinations are combined in the PDG average, the total uncertainty of the mean value is multiplied by a factor of two (|V us | = 0.2243 ± 0.0008) in order to account for the discrepancy among them.In our case, we could attribute the discrepancy to the new-physics effects.However, using only the determination from semileptonic decays generates a deviation from unitarity in the first row of the CKM matrix above the three-sigma level, when combined with |V ud | from nuclear β decays.This 'Cabibbo anomaly' is already present in the PDG average, although at a slightly lower level because the leptonic kaon decay pushes the central value of |V us | to a more favorable direction and the uncertainty is increased.In order to relax the deviation from unitarity below the three-sigma level, and also motivated by the PDG procedure, we follow a more conservative approach increasing also the uncertainty on |V us | by a factor of two (as well as the PDG does).The final value used in our fits is then Note, however, that the uncertainty is significantly higher than the one of |V us |, so the impact of this measurement on the CKM-fit is basically negligible.

B.1.4 Determination of |V cs |
Similarly to the previous case, the determination of |V cs | is obtained from measurements of semileptonic decays of D mesons and the leptonic decay of D s , provided that the form factors are obtained from lattice QCD computations.We have dismissed the determination from leptonic decays and we have used only the one coming from semileptonic decays: |V cs | = 0.972 ± 0.007 .(B.7) As happens for |V cd |, we have a much higher uncertainty compared to the light-quark data (|V ud |) so, again, this observable could be neglected from our fit leaving the results unchanged.meson mixing.Obviously, the additional scalars will contribute to this mixing but, once the ratio of the B d and B s transitions is taken, the new-physics contribution is highly suppressed since it is only present through SU(3)-breaking effects.We adopt the value for |V td /V ts | quoted by the PDG [101]:

B.2 Non-contaminated fit to the oblique parameters
For the determination of the oblique parameters a global fit of the EW observables is performed, removing the contribution from R b .We use the same inputs as described in Ref.  9. Results for the fit of the oblique parameters S, T and U , excluding the information from R b .

Figure 7 .
Figure 7.Comparison of the ς l region required to accommodate (g − 2) µ with the flavour constraints.

Figure 8 .
Figure 8. Correlations among the masses of the new scalars.

Figure 10 .
Figure 10.Correlations among the mixing angle and the Yukawa alignment parameters with the masses.

Figure 11 .
Figure 11.Correlation among the Yukawa alignment parameters.Also shown are the restrictions (2.18) imposed in the different Z 2 THDMs (the inert is omitted since in that case these couplings are simply set to zero).

Figure 12 .Figure 13 .
Figure 12.Correlation among the mixing angle and the Yukawa alignment parameters

Figure 14 .
Figure 14.Correlations among the masses of the new scalars, with and without including the CDF measurement of M W .

Figure 15 .
Figure 15.Correlations among the scalar mass splittings, with and without including the CDF measurement of M W .
ui + ς d m dj m ui + m dj , (B.3)and the SM contribution being related with the CKM matrix element byΓ SM (P + ij → l + ν l ) = G 2 F m 2 l f 2 P |V ij | 2

B. 1 . 3
|V us | = 0.2230 ± 0.0012 .(B.5) Determination of |V cd | The PDG [101] average is obtained from semileptonic decays of D mesons to light leptons, leptonic D decays to muons and taus, and from neutrino scattering data.The scalar contribution to the leptonic decay can be sizable and, therefore, we will only use the data coming from semileptonic decays (|V cd | = 0.2330 ± 0.0136) and the neutrino scattering data (|V cd | = 0.230 ± 0.011).Averaging these two values we obtain |V cd | = 0.231 ± 0.009 .(B.6)

Table 2 .
Marginalised individual results.The mass limits are at 95% probability while for the others we show the mean value and the square root of the variance.

Table 3 .
Higgs signal strengths measured by ATLAS and CMS.Note that not all the decays have been measured for all the production channels.

Table 4 .
Direct searches for charged scalars.

Table 5 .
Direct searches for neutral heavy scalars, φ 0 i = H, A, with final states including the SM Higgs boson or other neutral scalars.φ 3 denotes the heaviest scalar, V = W, Z, ℓ = e, µ.The parenthesis show the final decay of the SM particles produced from the NP particles.The square brackets are used when the values of σ • B are shown in terms of the primary decay (i.e. the NP particle decay) but a particular decay channel of the SM particles is used to obtain those values.See Section 4.2 for more details.

Table 6 .
Direct searches for neutral heavy scalars, φ 0 Zγ final states.The parenthesis show the final decay of the SM particles produced from the NP particles.The square brackets are used when the values of σ • B are shown in terms of the primary decay (i.e. the NP particle decay) but a particular decay channel of the SM particles is used to obtain those values.See Section 4.2 for more details.

Table 8 .
ts | = (0.207 ± 0.001 ± 0.003).Using as inputs all the measurements mentioned in the previous sections, we obtain the values of Tab. 8 for the Wolfenstein parameters.Wolfenstein parameters obtained from a fit to the CKM entries in Section B.1.