$B_s \to \mu^+ \mu^-$ in a Two-Higgs-Doublet Model with flavour-changing up-type Yukawa couplings

We present a Two-Higgs-Doublet Model in which the structure of the quark Yukawa matrices is governed by three spurions breaking the flavour symmetries of the quark Yukawa sector. The model naturally suppresses flavour-changing neutral current (FCNC) amplitudes in the down-type sector, but permits sizable FCNC couplings in the up sector. We calculate the branching ratio of $B_s \to \mu^+ \mu^-$ to leading and next-to-leading order of QCD for the case with FCNC couplings of the heavy neutral Higgs bosons to up-type quarks and verify that all counterterms follow the pattern dictated by the spurion expansion of the Yukawa matrices. We find correlations between $B_s \to \mu^+ \mu^-$, $b\to s\gamma$, and the Higgs masses. The $B_s - \bar B_s$ mixing amplitude is naturally suppressed in the model but can probe a portion of the parameter space with very heavy Higgs bosons.


Introduction
Two-Higgs-Doublet models (2HDMs) [1,2] are popular extensions of the Standard Model (SM) due to their relative simplicity, involving no additional fields apart from a second Higgs doublet.Moreover, a strong motivation to study 2HDMs also comes from theories in which a second Higgs doublet is required due to symmetry arguments, e.g.axion models in the context of the strong CP puzzle [3][4][5] or minimal supersymmetry [6].2HDMs differ in the structure of Higgs-fermion Yukawa couplings.The historically most favoured variants are the so-called type-I and type-II 2HDM in which both up-type and down-type quarks only couple to one of the two Higgs doublets.In these types of 2HDMs, flavour-changing neutral current processes (FCNC) such as the decay B s → µ + µ − are loop-suppressed and therefore small masses of the additional Higgs bosons are in principle possible, an appealing feature in the time of early LHC searches.The generic 2HDM (also dubbed "type-III"), with most general Yukawa matrices, exhibits, however, a much richer phenomenology [7][8][9] in which flavour-changing neutral Higgs-boson couplings are possible, in case up-type or down-type quarks couple to more than one Higgs doublet.The generic 2HDM suffers from a low level of predictivity caused by its large number of parameters and moreover requires that FCNC couplings in the down-type sector are tuned to small values in an ad-hoc way to comply with the many experimental constraints from s → d, b → s, and b → d FCNC processes.Yet the economic 2HDM of type I and II cannot address current anomalies in rare bottom quark decays, for instance data favour an excess of B(B → Dτ ν)/B(B → Dµν) over its SM prediction [10][11][12][13][14][15][16][17][18][19][20][21] and this cannot be explained in these models.Also to address the observed deficit in b → sµ + µ − for low values of the lepton-pair invariant mass one resorts to the generic 2HDM [9].
In this paper, we propose a three-spurion 2HDM, a 2HDM which is more general than the popular type-I and type-II models (which it contains as limiting cases) but more constrained than the generic 2HDM.A process on which our model has characteristic imprints is B s → µ + µ − , which we study in detail in this paper, including the calculation of two-loop QCD corrections.The model under consideration will feature flavour-changing uptype Yukawa couplings of the additional heavy neutral Higgs bosons, most notably a charmtop transition which may eliminate the Cabibbo-Kobayashi-Maskawa (CKM) suppression in b → s transitions that is present in both the SM and the 2HDM type-II models.
The outline of the paper is as follows: The Yukawa and Higgs sectors of the threespurion 2HDM are presented in Section 2, with particular emphasis on some peculiarities of this specific model.In Section 3, we will introduce the effective operators contributing to the low-energy weak B s → µ + µ − decay, while Section 4 is dedicated to a short description of the computational setup used for the evaluation of Feynman diagrams.The limiting case of the type-II 2HDM is discussed in Section 5, with the additional contributions from flavour-changing neutral-Higgs Yukawa couplings being presented in Section 6.Finally, we will discuss the phenomenology of such models in Section 7 and summarize in Section 8.

2HDM with suppressed down-type FCNC couplings
In this section we introduce the so-called three-spurion 2HDM, which allows for significant flavour-changing Yukawa couplings in the up-type quark sector, while the ones in the downtype quark sector are naturally suppressed.

General Yukawa sector
Our starting point, the quark Yukawa Lagrangian of the general ("type III") 2HDM, reads with four general complex 3 × 3 matrices Y u,d and ϵ u,d as well as and 3) The subscripts f, i = 1, 2, 3 label the generations, e.g.u ′ 3L = t ′ L .The notation of Eq. (2.1) follows Ref. [8], except that our H d corresponds to −ϵH * d of that paper.The vacuum expectation values (vevs) and related quantities are (2.4) The quark mass matrices are which we diagonalise in the usual way with the help of unitary matrices S u,d L,R , with the unprimed fields corresponding to quark mass eigenstates.The gauge sector of the 2HDM is invariant under independent unitary rotations of the fields Q ′ , d ′ R , and u ′ R in flavour space.We use Eq.(2.6) and choose and find L Y in the so-called down basis: with the appropriately transformed Yukawa matrices and the CKM matrix (2.11) The Yukawa Lagrangian L Y in Eq. (2.9) is manifestly SU(2) invariant, with the SU(2) doublet Eq. (2.9) is our starting point; the Yukawa matrices are related to the diagonal mass matrices as (2.12) Non-vanishing off-diagonal entries of Y The doublets ϕ SM , ϕ new of the Higgs basis [22][23][24] are defined by a rotation of the two doublets H u and H d by the angle β such that ϕ new has no vev: One has Next we use Eq.(2.13) to express L Y in terms of ϕ SM and ϕ new : By using Eq.(2.12) to eliminate Y u and Y d one can write the couplings to the physical charged and neutral Higgs bosons in Eq. (2.15) as with the matrices [8] (2.17) The non-diagonal matrices g d and g u characterise the deviations from the popular type-II 2HDM (for which g d = g u = 0) and can induce flavour-changing couplings of neutral Higgs bosons.Note that the type-I model is also included in the formalism and recovered by using Eq.(2.12) with Y u = 0 in the expression for g d .For our loop calculation and the phenomenological analysis it is advantageous to work with g d,u rather than ϵ d,u , especially for the definition of the renormalisation prescriptions.We write g d i d j ≡ g d ij , where d i is the i-th down-type quark flavor, i = 1, 2, 3, and similarly for g u .
We restrict ourselves to the CP-conserving Higgs potential, such that A 0 is a pseudoscalar boson, while ϕ 0 and ϕ 0′ are scalar particles.The two Higgs bosons ϕ 0 and ϕ 0′ are in general not mass eigenstates.The latter are given by h 0 and H 0 , with The angle α is a priori arbitrary, but data on decays of the 125 GeV Higgs boson constrain cos(β − α) to be close to zero.

Spurion expansion
The 2HDM of type I and II (and their variants with modified lepton couplings) invoke (softly broken) Z 2 symmetries to forbid FCNC couplings of the neutral Higgs bosons.This was motivated by the wish to find non-standard Higgs bosons at modern colliders, because for generic values of the Yukawa matrices in Eq. (2.9) constraints from FCNC processes like meson-antimeson mixing push these masses to values outside the reach of LEP, Tevatron, and LHC.With the absence of discoveries of non-standard Higgs bosons this line of arguments loses its appeal and the consideration of more general Yukawa sectors is in order.
The type-II model is the most studied variant of the 2HDM for two reasons: Firstly, it constitutes the tree-level Higgs sector of the Minimal Supersymmetric Standard Model (MSSM), in which the holomorphy of the superpotential enforces ϵ u,d = 0. Secondly, the type-II model is phenomenologically especially interesting, because in this model FCNC processes are sensitive to loop effects of the charged Higgs boson.A prominent example of the latter feature is the branching ratio B(B → sγ), which sets a stringent bound on the charged-Higgs mass [25].The type-II 2HDM further permits the possibility of large down-type Yukawa couplings, a scenario motivated by the possibility of bottomtop Yukawa unification.Such large-tan β scenarios are efficiently constrained by B(B s → µ + µ − ) [26][27][28] and we will come back to this topic in Section 5. Concerning the first abovementioned motivation, the Higgs sector of the MSSM is really richer than that of the type-II 2HDM: In the limit of infinitely heavy superpartners one encounters non-decoupling loopinduced Yukawa matrices ϵ u,d , an effect caused by the supersymmetry-breaking terms [29][30][31].Despite the loop suppression large effects are possible in FCNC processes with downtype quarks which involve the product ϵ d tan β [32][33][34][35][36][37] with huge impact on B(B s → µ + µ − ) [32][33][34][35]38].
The phenomenological constraints from meson-antimeson mixing and rare (semi-)leptonic decays place severe bounds on the off-diagonal elements of ϵ d , while those of ϵ u are essentially unconstrained except for ϵ u 12 and ϵ u 21 .This situation calls for a variant of the general 2HDM in which Y u,d and ϵ u are arbitrary, while ϵ d is suppressed.A strong motivation for such a model is the possibility of spontaneous CP violation, implemented through a Higgs potential developing complex vevs and real Yukawa matrices.Spontaneous CP violation is not possible with a 2HDM of type I or II, but requires at least three out of the four matrices in Eq. (2.9) to be non-zero.Avoiding fine-tuning implies that the dominant piece of the needed effect stems from ϵ u , while ϵ d can be neglected [39].However, in such a model, the mixing of the neutral Higgs fields is different from Eq. (2.18) and instead involves all three fields.Yet for the CP-conserving observables considered in this paper this feature is of minor relevance.The 2HDM scenario with sizable Y u,d and ϵ u has the appealing feature that it simultaneously permits both measurable effects in FCNC processes and sufficiently light masses of the non-standard Higgs particles enabling their discovery at the LHC.
Setting ϵ d naively to zero leads to a non-renormalisable model, because there are UVdivergent loops involving up-type quarks with Y u,d and ϵ u couplings, which require counterterms proportional to elements of ϵ d .Whenever one seeks to constrain the elements g u jk of Eq. (2.16) from FCNC transitions of down-type quarks, one must foresee such a counterterm to find a meaningful prediction.For example, B s → µ + µ − is a b → s transition constraining g ct and the corresponding loop contribution requires a counterterm for g sb .The minimal renormalisable theory is found by invoking flavour symmetries and systematically expanding in terms of the spurions breaking these symmetries [40,41].The 2HDM gauge sector is invariant under unitary rotations 3 and the Yukawa sector is formally invariant under this flavour symmetry if one transforms the matrices in Eq. (2.9) as We propose to categorise the classes of renormalisable 2HDM in terms of the spurions breaking the SU (3) 3 flavour symmetry of the quark sector. 1 The minimal choice are two spurions, with just two physically distinct possibilities.Both comply with the definition of minimal flavour violation (MFV) as defined in Ref. [41].The first possibility is to take Y u,d as spurions and express the other two Yukawa matrices as , where P u and P d are polynomials.This variant is discussed in Ref. [41] and amounts to a generalisation of the 2HDM of type II.It also constitutes a renormalisable extension of the aligned 2HDM of Ref. [42][43][44], in which P u,d are constants.If the 2HDM is the low-energy limit of a more fundamental theory obeying the considered two-spurion symmetry-breaking pattern, the latter will naturally induce ϵ u,d in the described way as well.An example for such a UV theory is the MSSM with a flavour-blind supersymmetry breaking mechanism (such as gauge mediation).The second possibility is to choose Y d , ϵ u as spurions, leading to a generalisation of the type-I 2HDM.
There are two possibilities for a 2HDM with three spurions, which can be taken as Y u,d plus either ϵ u or ϵ d .The first possibility is the phenomenologically interesting one and studied in this paper.The expansion up to third order reads with complex coefficients c, . . ., b 22 .Concerning Eq. (2.20) several comments are in order: • The spurion expansion is only meaningful, if the contributions to the off-diagonal elements of ϵ d from higher electroweak orders (with five or more Yukawa matrices) are small, so that they can be neglected.A sufficient condition for this is realised in scenarios in which c 11 , . . ., b 22 are induced by one-loop contributions in either the UV completion or the 2HDM, while terms with (2n + 1) spurions are only generated at n-loop order and beyond.We consider this scenario throughout this paper. 2Additional QCD corrections (e.g. an extra loop with a gluon) comply with the pattern in Eq. (2.20), i.e.QCD renormalises the coefficients, but does not induce new ones.
• By rotating (H u , H d ) in Eq. (2.9) one can eliminate c in Eq. (2.20).But in general radiative corrections bring this term back and a counterterm to c is needed, unless one corrects the rotation in each order of perturbation theory.It is therefore advisable to keep c in Eq. (2.20); we treat it as a perturbative quantity with c = 0 at tree level.
• The decay B s → µ + µ − , which is the focus of the phenomenological analysis in this paper, involves the FCNC vertex • With Eq. (2.12) we can trade Y u,d in Eq. (2.20) for the quark masses and CKM elements.Compared to the SM we find 14 additional complex parameters, the 9 entries of ϵ u and c 11 , . . .b 22 .Yet it is much more convenient to express observables in terms of g u jk of Eq. (2.16) instead of ϵ u jk and then use Eqs.(2.17) and (2.20) (with Eq. (2.12)) to calculate the g d jk in terms of the coefficients of the spurion expansion.While this procedure is needed in a global analysis of all available data -which is beyond the scope of this paper-, the study of b → s transitions alone will simply involve g ct and, with CKM suppression, g ut and g tt .
In a practical calculation it is cumbersome to implement the renormalisation procedure in the described way, by providing counterterms to ϵ u jk and c, . . ., b 22 .Instead, it is much easier ϵ u are almost aligned, so that Eq. (2.12) means small off-diagonal matrix elements of these matrices in the chosen basis.Since Y u 33 , ϵ u 33 = O(1) and further Y d 33 = O(1) is possible, some terms with five spurions must be added to Eq. (2.20), just as in the MFV case of Ref. [41].
to renormalise the g u,d jk .If one renormalises all g u,d jk in the MS scheme, one automatically complies with SU(2)×U(1) gauge symmetry.Therefore it is sufficient to choose the g d jk in such a way that Eq. (2.20) is obeyed at tree level.In our calculation we will only need a counterterm to g sb (in addition to the usual QCD counterterms for the SM parameters), if m s is set to zero.If m s is kept non-zero, an additional counterterm to g bs is required.
In the three-spurion 2HDM the parameter tan β is well-defined, because unitary rotations of (H u , H d ) would lead to ϵ d ̸ = 0 at tree level and spoil the spurion expansion.We are interested in phenomenologically interesting scenarios, in which the rare decays B(B s → µ + µ − ) or b → sγ deviate from the SM predictions at a level probed in current and forthcoming measurements.The corresponding decay amplitudes involve a helicity flip and come with the Yukawa matrix Y d , which grows with tan β, so that we will consider the case that tan β is large.An interesting feature of the three-spurion 2HDM is that the above-mentioned amplitudes scale differently with tan β than in the type-II model, due to new loop contributions with g ct .
Next we briefly discuss the leptonic Yukawa Lagrangian.The extremely stringent experimental bounds on FCNC transitions like ℓ j → ℓ k γ suggest that only one spurion is present in the lepton sector, meaning that the lepton Yukawa couplings of the two Higgs doublets are automatically diagonal in the mass eigenstate basis. 3For simplicity, we choose the leptonic Lagrangian of type-II form (i.e.we set g l = 0) , with the familiar tan β-enhanced non-standard-Higgs couplings to charged leptons: The diagonal mass matrix of the charged leptons is denoted by M l and the couplings of ϕ SM are omitted in Eq. (2.21).
3 The decay B s → µ + µ − The typical momentum scale for B s decays is of order M Bs or smaller, so that weak B s decays can be described by an effective theory in which the heavy W, Z bosons, the top quark, and the Higgs bosons of the 2HDM are integrated out.The resulting |∆B| = 1 Hamiltonian H eff describes the interactions mediated by these heavy particles in terms of dimension-6 operators changing the beauty quantum number B by one unit.The piece of The operators in Eq. (3.1) are and are multiplied with their respective Wilson coefficients C A , . . ., C ′ P which contain the dependence on the large masses.The normalisation factor in Eq. (3.1) is which complies with the conventions of Ref. [28].The second "=" sign only holds to lowest order in the electroweak interaction, while in higher orders the relation between the Fermi constant G F and the electromagnetic coupling α em = e 2 /(4π), the weak mixing angle θ w and the W mass M W is modified.Electroweak corrections have been calculated in Ref. [45] and e.g.remove the ambiguities related to the choice of the renormalisation scheme for these parameters; in the second version for N in Eq. (3.3) this issue also includes the choice of the scale in the running α em .In the following, we choose the first definition N ∝ G 2 F in this paper, for which the electroweak corrections to the SM contribution to C A are as small as −2.4% [45].
We introduce the perturbative expansion of the C i as where C (0) i denotes the leading order (LO), arising in the SM from one-loop electroweak diagrams.C F in the SM, hence the branching ratio is rather small.The average time-integrated branching ratio is given by [46,47] with dimensionless quantities Here, Γ s H (Γ s L ) denotes the decay width of the heavier (lighter) B s mass eigenstate, and the factors F P and F S account for the mixing of the B s − Bs system, given by in the absence of significant CP-violating New Physics contributions to the B s − Bs mixing amplitude.In writing Eq. (3.5), we have used the pseudoscalar decay constant f Bs to rewrite the operator matrix elements as where the second equation follows from the first one by use of the equations of motion.The Higgs-mediated contributions in the SM can be neglected due to the tiny Yukawa couplings of external particles.However, in the 2HDM with large Yukawa couplings of Higgs bosons to right-handed down-type quarks and leptons, Feynman diagrams with Higgs bosons are known to contribute significantly to C S and C P , as will be discussed in the following.

Computational setup
In this section, we describe the chain of programs used to generate and evaluate the corresponding Feynman diagrams at leading and next-to-leading order.We use FeynRules with the Universal FeynRules Output (UFO) [48][49][50] to obtain Feynman rules for the model.The output is processed by tapir [51] into a Lagrangian file, which is used with qgraf [52] in order to generate all Feynman diagrams.
We compute the one-loop diagrams for general electroweak gauge parameters for the W and Z bosons and check that they drop out in the final result.This is a welcome check for the conversion from UFO to qgraf and tapir.If diagrams with all different quark flavours are explicitly calculated, we have at one-loop level < O(100) Feynman diagrams in the SM, 4 and an additional O(300) diagrams from contributions with at least one non-SM Higgs boson, most of which are Higgs-penguin diagrams.At two-loop order, we perform the calculation in the Feynman gauge for the W and Z bosons, but we keep the gluon gauge parameter general and verify that it drops out of the final results.
The diagrams are then converted with the help of tapir into FORM [53] code using Feynman rule definitions that were also produced in the conversion from UFO to the qgraf Lagrangian file.The individual expressions for the diagrams are mapped onto integral families using exp [54,55] and a custom FORM setup is used to perform the remaining computational steps.Since the Wilson coefficients are independent of the momenta of the external particles, we set the latter to zero, so that only vacuum integrals need to be computed.(An exception are diagrams with an FCNC self-energy in an external leg, which are calculated with on-shell b quark and m b ̸ = 0 before the subsequent limit m b → 0 is taken.)We use a FORM implementation of the algorithm presented in Ref. [56], see Ref. [57].
5 The decay B s → µ + µ − in the Two-Higgs-Doublet Model of type-II In the SM the leading-order result was obtained in Ref. [58] and next-to-leading order QCD corrections have been presented in Refs.[59][60][61][62].Higher-order QCD and electroweak b s Sample diagrams contributing to C A at leading and next-to-leading order in the SM.
corrections have been computed in Refs.[27,28,45,63,64].Next-to-leading corrections in the type-II 2HDM have been calculated in Refs.[26,[65][66][67].We have reproduced these results for the present paper, also as a check of the automated setup.The result can be expressed in terms of dimensionless mass ratios of the top quark, W boson and charged Higgs boson masses,  and next-to-leading order they are given by where the Z-penguin contributions include the flavour-changing quark self-energy diagrams required to make the penguin diagrams finite.Moreover, since charm and up quark masses can be neglected compared to M W and m t , after summation over internal up-type quarks all contributions are proportional to the CKM factor V ts V * tb (contained in the normalisation factor N defined in Eq. (3.3)) due to the unitarity of the CKM matrix.
In the 2HDM of type-II the leading new effects stem from O tan 2 β contributions to S and C (′) P , arising from penguin diagrams with a neutral Higgs boson, see e.g.Fig. 2, and the W + -H + box diagram [26].The Higgs-penguin contributions to C S are given in Feynman gauge for the W -boson by where (in C A ) terms arising from Z-penguin diagrams with a charged Higgs boson.Box diagrams with a single charged Higgs boson also contribute to C S,P (C ′ S,P ), with Wilson coefficients proportional to m b (m s ).We do not explicitly list these contributions here; they can be found in Refs.[65][66][67].With the exception of the doubly muon-mass suppressed H + -H − box contributions to C A , we include all of these additional terms in our analysis.The feature that C S and C P are proportional to m b while their primed counterparts are proportional to m s holds beyond the type-II 2HDM in our more general 2HDM with three spurions because of A rather remarkable feature of the type-II 2HDM is the fact that the leading terms in tan β depend only on tan β and the charged-Higgs boson mass M H + , that is they are independent of the parameters of the neutral Higgs sector [26].In the type-II 2HDM, the leading tan 2 β contributions to the pseudoscalar and scalar Wilson coefficients satisfy the rather simple relation (5.7) In Fig. 3 we show the branching ratio B (B s → µ + µ − ) in the type-II 2HDM as a function of the charged Higgs boson mass.The horizontal green, grey, and violet bands correspond to the experimental measurements of LHCb [85] and CMS [70] as well as the theory prediction [27], respectively, including 2σ uncertainties for the experimental values (cf.Table.1).The theory prediction of Ref. [27] uses |V cb | = 0.0424 ± 0.0009, which is close to today's value inferred from inclusive b → cℓν decays.If one uses |V cb | = 0.03936 ± 0.00068 from exclusive B decays [86] the central value of the theory prediction for 10 9 B (B s → µ + µ − ) drops from 3.65 to 3.15.
The coloured lines are predictions from the 2HDM for different values of tan β.It is interesting to note that for tan β ≲ 25 low values of M H ± are required to reproduce the central value of the LHCb measurement.Of course, in the limit M H ± → ∞ all 2HDM curves approach the SM prediction.Note that in the type-II model there is a tan βindependent 95% C.L. lower bound on M H ± in the range 570-800 GeV from B → X s γ [25].This bound can be easily weakened with our model's additional couplings discussed in the following section.
Recently it has been pointed out that LHC data still permit M H ± ≤ 400 GeV with couplings compatible with solutions of the b → cτ ν flavour anomalies [87,88].Charged-Higgs explanations of the latter have been found viable in Refs.[89,90] and are invigorated by recent LHCb data on B → D ( * ) τ ν [91], see Refs.[92,93].It should be clearly stated that this solution to the b → cτ ν puzzle is not realised in the type-II model, in which for instance B(b → cτ ν) is suppressed rather than enhanced over B(b → cµν) as preferred by data.Also while charged Higgs searches at the LHC are compatible with the quoted H ± masses [94], the lower bound on M H ± inferred from the data on gg/b b → A 0 [→ τ + τ − ] searches is larger than 1 TeV [95] in perturbed versions of the type-II 2HDM.Also in our model we cannot substantially weaken this bound; to this end one must modify the A 0 coupling to τ 's by deviating from the type-II form in Eq. (2.21) and permitting the third lepton generation to couple to the other Higgs doublet.

Additional contributions in a model with flavour-changing neutral Yukawa couplings
In a model with a Yukawa Lagrangian given by Eq. (2.16) there are additional tree-level contributions bs → h 0 (H 0 , A 0 ) → µ + µ − of order tan β.A sample Feynman diagram is shown in Fig. 4. At loop-level there are O(tan 3 β) terms due to diagrams in which the neutral Higgs boson couples to the b line.At LO these are diagrams with a FCNC selfenergy and we also refer to them as self-energy diagrams at NLO, even if a gluon connects b s µ + µ − h 0 , H 0 , A 0 the FCNC loop with the s quark as in the diagram on the right in Fig. 2.These O tan 3 β terms occur because a tan β-enhanced coupling in the self-energy diagrams is not cancelled by a factor cot β in the second charged-Higgs coupling, which is a distinguishing feature compared to the type-II model.If these self-energy diagrams involve a helicity flip of the internal fermion line, see Fig. 5, they come with a factor of tan β and are linear in the flavour-changing Yukawa matrix g u , which enters the result through the charged-Higgs coupling in the fourth term of Eq. (2.16).In fact, the dominant dependence on g ct stems from this source and not from diagrams in which a neutral Higgs boson couples to charm and top in a vertex diagram, which have a factor of tan β less.Note that contributions from diagrams without helicity flip are either already included in the type-II model or quadratic in g u (and without the factor tan β), and we will consequently neglect the latter in what follows.

Pseudoscalar Wilson coefficient C P
At tree level the Wilson coefficient C P originating from the diagram in Fig. 4 with a virtual A 0 boson is given by Since there is no loop-suppression, even small values of g bs and g * sb can have significant impact on the branching ratio, but the spurion expansion in Eq. (2.20) naturally leads to (parametrically suppressed) small couplings.For convenience, let us define the linear combinations g ± bs ≡ g bs ± g * sb .At one-loop order only the self-energy diagrams contain an enhancement factor tan 3 β whereas the vertex contributions contain at most a quadratic term.The one-loop selfenergy contributions to C P and C ′ P are ultra-violet divergent.The corresponding counterterm is generated from Eq. ( 6.1) and is of the form g +,0 bs = g + bs + δ A 0 ,bs + . . . .(6.2) In the MS renormalisation scheme the one-loop contribution is given by δ Defining analogously A 0 ,bs + . . ., (6.4) the counterterms δ A 0 ,bs and δ A 0 ,bs are found from δ A 0 ,bs and δ A 0 ,bs in Eqs.(6.3) and ( 6 with normalisation factor the weak coupling constant g, and Our results in Eqs.(6.1-6.7)confirm the result in Eq. (3.16) of Ref. [9].At NLO in QCD (i.e.two-loop order), the counterterm required to obtain a finite result is given by Note that in addition to the top quark mass m t also the flavour-changing Yukawa couplings are minimally renormalised.The two-loop contribution to the Wilson coefficients reads with C(1) P involves no large logarithm for the choice µ ∼ M H + .To verify this expand Eq. (6.7) as This feature is generic for all log µ 2 terms stemming from loops with Yukawa couplings, because heavy particles like H + do not contribute to the renormalisation group functions (i.e.β functions and anomalous dimensions) for µ < M H + since heavy particles are integrated out at scales of the order of their masses.The same loop diagrams yielding C(0) P also determine the piece of the β functions of g sb proportional to the Yukawa couplings of top and bottom quarks.The running of g sb (µ) from this source is compensated by the explicit logarithm in Eq. (6.11) and the remaining Yukawa-µ dependence is a tiny two-loop effect.
The situation is different with the µ dependence stemming from gluon loops and related to the familiar QCD running of couplings and quark masses.The QCD running of g sb is trivial, since the combination g sb (µ)/m b (µ), which enters C P •⟨Q P ⟩, is independent of µ (see Eq. (3.8)).To study the µ dependence of our two-loop result C (1) P , we keep the log µ 2 term stemming from the Yukawa interaction (i.e. the analogue of the logarithm in Eq. (6.11)) fixed and vary µ otherwise.C (0) P depends on µ implicitly through the µ-dependence of m t and g it and this dependence should be compensated by the explicit log µ 2 terms in C (1) P , reducing the µ-dependence to the three-loop (one-loop Yukawa correction and NNLO QCD) level. 5In Fig. 6 P (M H + ) by 24%, while the corresponding value for C (1) Wilson coefficient C P for a particular choice of M H + and g ct ≡ g ct ( mt ).All parameters that are not running in QCD have been fixed in this figure and we have neglected the contributions from g ut and g tt , which have the same running as g ct , such that the only running parameters are in m t g * ct C(0,1) P .The figure shows a significant improvement of the QCD scale dependence with the inclusion of next-to-leading order QCD corrections.The LO result does not permit a reliable prediction of C P , while C NLO P merely changes by ±5% around its central value when µ is varied between m t and M H + .
From Fig. 6 it is clear that our calculated QCD corrections are needed for a reliable prediction.Next we discuss uncalculated higher-order corrections involving Yukawa couplings, obtained by dressing the LO diagrams with an additional Higgs boson.The large O(1) couplings are the coupling of the SM-like Higgs boson h 0 to the top quark and the A 0 , H 0 couplings to the bottom quark. 6The former contributions are already present in the SM, contained in the electroweak corrections of Ref. [45].Since they are very small in the SM, they will constitute an even smaller correction to the extra diagrams of the 2HDM.The dominant contribution from an extra loop with A 0 or H 0 is expected from diagrams in which both ends of the additional Higgs line are connected to a b line.Contrary to the QCD case, all A 0 , H 0 diagrams involve momenta which are far off-shell, because M A 0 ,H 0 is much larger than m b .Thus the additional loop will give a correction in the few-% region to the leading 2HDM term, which is numerically constrained to the range between SM prediction and experimental value and thereby constitute a small correction to a small LO 2HDM term.

Scalar Wilson coefficient C S
The scalar Wilson coefficients receives contributions from both neutral CP-even Higgs mass eigenstates h 0 and H 0 .At tree level, the diagrams with h 0 and H 0 give rise to the Wilson coefficients where contains the dependence on the neutral Higgs-boson masses.The counterterms required to cure the divergences at one and two loops can be obtained from Eqs. (6.3) and (6.8) analogously, i.e. the combination g + bs − g − bs renormalises C S , while the combination g + bs + g − bs renormalises C ′ S .The renormalised Wilson coefficients C S and C ′ S are related to the pseudoscalar ones by Note that there are no QCD corrections to C (′),tree S,P since they cancel in the matching calculation.

Phenomenology
In this section we discuss the possible size of B (B s → µ + µ − ) in our 2HDM under the constraint that other b → s processes comply with the data.We include the tree-level contributions from diagrams with flavour-changing down-type couplings (cf.Fig. 4), as well as the SM results and the leading quadratic tan β contributions in the type-II 2HDM to which the additional diagrams of order tan 3 β discussed in the previous section add as corrections.
In a generic 2HDM the tree-level couplings g ± bs will drastically increase the branching ratio for B s → µ + µ − due to the missing loop suppression.In our model the spurion expansion suppresses g ± bs in a controlled way, but still permits large enough contributions to get phenomenologically interesting effects in B s → µ + µ − : Even rather small up-type couplings g ct significantly modify the Wilson coefficients C (′) P and C (′) S , as they feature a CKM factor V cs instead of V ts .In the following, we will restrict ourselves to the experimentally favoured scenario [96,97] of aligned Higgs doublets, in which sin (β − α) ≈ 1.For the numerical analysis, we will set sin (β − α) = 1, and therefore have R M = M −2 H 0 .In this case, the Higgs mass eigenstates h 0 and H 0 coincide with ϕ 0 and −ϕ 0′ , respectively, and only the latter possesses non-SM-like couplings to fermions.This reduces the number of relevant non-Yukawa-type parameters of the extended Higgs sector to four, namely tan β, M H + ,

Constraints from b → sγ decays
An important constraint on the magnitude of flavour-changing Yukawa couplings in the up-type quark sector arises from the inclusive rare decays B → X s γ.This process is mediated at the quark level by b → sγ through a top quark loop with a charged W boson in the Standard Model and receives additional contributions through charged-Higgs boson diagrams in the 2HDM, see Fig. 7.In order to eliminate the dependence on V cb ≃ −V ts , one traditionally works with the quantity Note that in the SM and the type-II 2HDM the branching ratio of b → dγ is much smaller than that of b → sγ.
In the type-II 2HDM, the contribution of the charged Higgs diagrams to the decay rate is positive, moving the theoretical prediction for the branching ratio away from the experimental value averaged in [25] to R γ,exp.= (3.22 ± 0.15) where a lower cutoff of E 0 ≥ 1.6 GeV has been put on the photon energy.The practical independence of R γ of tan β in the largest part of the parameter space in the type-II 2HDM allowed for the extraction of a lower limit on the mass of the charged Higgs boson in [25], see Section 5.In the general 2HDM, where, contrary to the type-II case, the factor tan β from the btH − vertex is not cancelled by a factor cot β from the stH − vertex, the quantity R γ depends on tan β. 7 The combination of the tL b R H + and sL t R H − Yukawa couplings arising in the process is then given by where we have defined the short-hand notation For the case of the type-II model (or ramifications of it) these searches already exclude a significant portion of the M A 0 − tan β plane at present.In Fig. 9, we show a recent collection of exclusion limits obtained through various different searches by the ATLAS experiment [95] in a modification of the type-II model.These limits imply that large values of tan β can only be realised if at the same time the additional Higgs bosons are quite heavy, M A 0 ≳ 1.5 TeV.
It should be noted that for heavy M A 0 these bounds also approximately apply for M H + and M H 0 , since the masses become degenerate in the heavy-Higgs limit (for a thorough analysis of the allowed mass splittings see [71]; see also the vertical dashed lines in Fig. 3).Note that the bounds from neutral Higgs searches are more constraining than those for H + searches [94].The most stringent constraints are from final states with τ 's and also apply to our lepton Yukawa Lagrangian in Eq. (2.21) which we have chosen of type-II.(However, the choice of the τ Yukawa couplings is of no relevance for the phenomenology of the b → s FCNC processes discussed in this paper.)The presence of the additional couplings g u jk will increase some branching ratios at the expense of those of decays into τ 's, so that these bounds can be somewhat weakened, but the general trend remains valid.Our scenarios discussed below comply with the ATLAS bounds.

B s − Bs mixing
The effective sL b R A 0 and sL b R H 0 vertices are the sum of the tree-level couplings ∝ g sb and the loop contributions involving vertex and self-energy diagrams ts .The coefficients C P and C S are proportional to this effective vertex, which is therefore the relevant quantity constrained by B(B s → µ + µ − ).Now the B s − Bs mixing amplitude depends quadratically on this effective vertex, while both quantities decrease quadratically with M A 0 and M H 0 , so that the quantities give complementary information on the parameter space.B s − Bs mixing is mediated in the Standard Model by the |∆B| = 2 operator Q V LL = bL γ µ s L bL γ µ s L and in the 2HDM by three additional effective operators The Wilson coefficient of Q SLL only involves the effective vertex governing B s → µ + µ − and the tree-level propagators of A 0 and H 0 , while the coefficients of the other operators also involve the corresponding chirality-flipped effective sR b L A 0 and sR b L H 0 vertices.The coefficient C SLL (C SRR ) is proportional to m 2 b (m 2 s ), while the coefficient C SLR is proportional to m b m s ; we therefore drop the operator Q SRR in our analysis.However, considering the ATLAS constraints implying large masses H 0 is heavily suppressed due to the small splitting of Higgs masses in the large-mass limit, and C SLR becomes relevant.This feature results from the fact that Q SLL violates hypercharge by two units of the Higgs hypercharge and the coefficient C SLL is therefore suppressed by a factor of v 2 /M 2 A 0 compared to C SLR multiplying Q SLR which conserves hypercharge and SU(2).As a consequence, C SLR is more important than C SLL .This feature has been widely studied in the context of the effective 2HDM emerging from integrating out superpartners in the MSSM [33][34][35]37].
In the following we correlate B(B s → µ + µ − ) with B(B → X s γ) and the B s − Bs oscillation frequency ∆M Bs , which is proportional to the magnitude of the B s − Bs mixing amplitude.Only the effective sL b R A 0 and sL b R H 0 vertices are physical, by e.g.changing the renormalisation condition for g sb we can shift pieces between g sb and the renormlised loop.For simplicity, we set the tree-value of g sb to zero, i.e. consider the case that this coupling is only generated radiatively.The only non-trivial Yukawa structure entering the considered observables is then g eff st .In Fig. 10 and Fig. 11 we illustrate the dependence of the branching ratio B (B s → µ + µ − ) on the flavour-changing Yukawa couplings for some exemplary numerical values.We choose to show our results separately for the LHCb [68,69] and CMS [70] measurements, which lie on different sides of the SM prediction calculated with the value |V cb | taken from inclusive decays.The interval allowed on the horizontal axis (the real part Re g eff st ) is for each choice of Higgs masses determined by the range allowed by b → sγ.Recall that g eff st carries a factor of tan β, which needs to be taken into account if the bounds on g eff st were to be converted into direct bounds on g ct .A priori generic ≤ O(0.1) values of g ct could make g eff st in Eq. (7.4) as large as O(50).For the considered Higgs masses such large values of |g eff st | are forbidden by b → sγ and the range for Im (g eff,bsγ ) shown on the y axis in plots (a) to (d) only serves the purpose to show the constraint from B (B s → µ + µ − ).Low values for Higgs masses enforce small values for tan β from the LHC searches, but for these scenarios b → sγ forbids any measurable effect in B (B s → µ + µ − ).This situation changes if one considers large values for tan β and the Higgs masses, because B (B s → µ + µ − ) grows faster with tan β than B(B → X s γ), see plots (e) and (f).We find that in the considered scenarios the bounds from B(B → X s γ) are stronger than those from ∆M Bs , which we always find in the band allowed by the current theoretical uncertainty.This situation changes if we go to even larger Higgs masses, permitting larger values of |g eff st | in B(B → X s γ).The           Table 1.Numerical input used for the phenomenological analysis.For the QCD running of the top-quark mass and the renormalisation scale, we have used RunDec [108,109].The numerical values of the CKM matrix elements have been taken from the updates provided on the CKMfitter [110] web page.The numerical values of the B s decay widths have been taken from the online updates provided at the HFLAV web page [86].
corresponding diagram (FCNC self-energy with the neutral Higgs attached to the b line) is enhanced by a factor of tan β compared to the vertex diagram, resulting in a O tan 3 β contribution to the B s → µ + µ − amplitude which is absent in the type-II model.This feature makes B s → µ + µ − a sensitive probe of the model even for Higgs masses well above the lower bounds found by the LHC experiments.For small Higgs masses, however, b → sγ precludes large effects in B s → µ + µ − .In our model the dominant contribution to B s − Bs mixing is naturally small due to a suppression factor of m s m b /v 2 ; nevertheless B s − Bs mixing sets constraints on the parameter space for very large Higgs masses.
3).The dominant one-loop vertex diagram involves an internal H d line and the product Y d Y d † Y d or ϵ u ϵ u † Y d stemming from the three H d Yukawa couplings.The UV divergences can be cancelled by counterterms to c 11 and b 11 .
comprises the next-to-leading order (NLO) QCD corrections.In the SM onlyC A is relevant, C ′ A is suppressed w.r.t.C A by the ratio m b m s /M 2 W involving the strange and bottom quark masses m s,b and C (′) S,P receive additional suppression factors of M 2 Bs /M 2 W .The leading contributions C (0) i arise from one-loop electroweak diagrams at order G 2

. 1 )
In the SM, only the Wilson coefficient C A receives significant contributions from diagrams such as the ones depicted in Fig. 1.The Wilson coefficients C ′ A , C S , . . ., C ′ P are suppressed by powers of ratios of light (external) masses and M W .The leading SM contribution is, moreover, independent of the Yukawa couplings of all external particles, that is we can set m b = m s = m µ = 0 for the contributions from W -box and Z-penguin diagrams.At leading b s µ

Figure 2 .P
Figure 2. Sample two-loop Feynman diagrams contributing to the Wilson coefficients C (′) S (H 1,2 ) and C (′)P (A 0 ).In the left diagram, only flavour-diagonal transitions i = j are possible in the type-II 2HDM, while the more general 2HDM allows also for transitions with i ̸ = j, e.g.transitions from a charm quark into a top quark.

. 6 ) 4 W and m 2 t m 2 µ /M 4 W 2 W
The Wilson coefficients for the right-handed operators C ′h,(i) S can be obtained by the replacement r b → r s in Eqs.(5.4) and (5.5).Further contributions include terms of order m b m s m 2 µ tan 4 β/M entering the Wilson coefficients C ′ A and C A , respectively, as well as O m b m s tan 2 β/M (in C ′ A ) and O m 2 t cot 2 β/M 2 W

Figure 3 .
Figure3.Branching ratio B (B s → µ + µ − ) in the type-II Two-Higgs-Doublet model at next-toleading order in QCD, for different values of tan β.All running parameters are evaluated at the scale µ = m t (m t ).The green dashed and dotted lines denote the central value and the experimental 2σ error of the LHCb measurement[68,69], respectively, the grey band shows the corresponding information for the CMS result[70].The purple dashed and dotted bands indicate the SM prediction with uncertainties, presented in Ref.[27].Note that the SM prediction in the purple band was obtained including also O α 2 s and O (α em ) contributions and corresponds to |V cb | = 0.0424±0.0009.The perturbativity of the scalar 2HDM potential constrains M 2 A 0 − M 2 H ± ≤ |λ 4 − λ 5 |v 2 ≲ 8v 2 [71-84], which can be converted into a lower bound on M H ± by use of the experimentally excluded area of the (M A 0 , tan β) plane such as the ones shown in Fig. 9 used by us.The vertical dashed lines indicate these lower limits on M H ± for each value of tan β, whereas the vertical dotted lines show the same limits obtained with the less stringent constraint |λ 4 − λ 5 | ≤ 8π.The plot shows that B s → µ + µ − will only give contraints on the type-II 2HDM competitive with the Higgs searches, once the experimental precision on B (B s → µ + µ − ) substantially improves.

Figure 4 .Figure 5 .
Figure 4. Tree-level diagrams with flavour-changing neutral Higgs couplings.The b − s − H coupling is denoted by a dot on the vertex.

. 8 )
below by the replacement m b → −m b .Thus the terms with m b renormalise the Wilson coefficient C P ∝ g + bs − g − bs , while the ones with m s renormalise the Wilson coefficient C ′ P ∝ − g + bs + g − bs .The renormalised finite pseudoscalar Wilson coefficients read C

Figure 7 .
Figure 7. Sample Feynman diagram for the rare decay b → sγ at one loop in the 2HDM.The cross indicates a chirality flip t L → t R .

Figure 8 .
Figure 8.The ratio R γ as a function of the charged Higgs boson mass M H ± and g eff st .The dashed lines correspond to the experimental central value and the ±2σ intervals, to which we have also added the theoretical uncertainties.

Figure 10 .
Figure 10.Branching ratio of B s → µ + µ − for different values of Higgs masses and tan β.All quantities are evaluated at µ = m t (m t ).The red dashed and dotted lines indicate the experimental central value and the 2σ uncertainties of LHCb branching ratio measurement[68,69], respectively.The constraint of B(B → X s γ) on |Im g eff st | is not shown.
d,u , ϵ d,u give rise to FCNC couplings of the neutral components of the Higgs doublets.In a general 2HDM the quantity tan β has no physical meaning: One can arbitrarily rotate H u , H d in Eq. (2.1) leading to a Lagrangian of the same form (yet with different Yukawa matrices) and the rotation angle will add to β.The situation is different in variants of the 2HDM in which H u and H d are distinguished by quantum numbers which forbid such rotations of H u , H d .Prominent examples are the 2HDM of type I and II, in which two out of the four Yukawa matrices in Eqs.(2.1) and (2.9) are absent.The type-I model corresponds to ϵ d .10)Next we discuss the dependence of our result on the renormalisation scale µ at which the 2HDM result is matched to the effective |∆B| = 1 Hamiltonian and the size of higher-order corrections.In the case m t ∼ M H + the choice µ = O(m t , M H + ) leads to the absence of large logarithms and the variation of µ between m t and M H + does not constitute a relevant source of theoretical uncertainty.Thus we are left to the phenomenologically interesting case M H + ≫ m t .The log µ 2 terms in the Wilson coefficients have two different origins, divergent contributions involving Yukawa couplings or the QCD coupling, respectively.C P stems from a loop with Yukawa couplings and C we illustrate the QCD scale dependence of the O α 0 s and O (α s ) Figure6.QCD scale dependence of the Wilson coefficient C P at leading and next-to-leading order in α s .In this plot we have used M H + = 1.5 TeV and g ct (m t ) = 1 and have further set V ts /V cs → 0. All µ-independent prefactors have been fixed.C P (µ) depends logarithmically on the masses m t and M H + , so that every choice of µ in the interval [m t , M H + ] seems justified.The plot shows, however, that choosing µ ∼ M H + in C P (µ) would badly underestimate the NLO result, while µ ∼ m t would sizably overestimate it.C [70] → µ + µ − ) SM 3.65(23) × 10 −9 [27] B (B s → µ + µ − ) LHCb 3.09 +0.48 −0.44 × 10 −9 [68, 69] B (B s → µ + µ − ) CMS3.95 +0.52 −0.47 × 10 −9[70]