Two Higgs doublet models and $b\rightarrow s$ exclusive decays

Abstract

We computed the leading order Wilson coefficients relevant to all the exclusive $b\rightarrow s\ell ^+\ell ^-$ decays in the framework of the two Higgs doublet model (2HDM) with a softly broken $\mathbb {Z}_2$ symmetry by including the $\mathcal {O}(m_b)$ corrections. We elucidate the issue of appropriate matching between the full and the effective theory when dealing with the (pseudo-)scalar operators for which keeping the external momenta different from zero is necessary. We then make a phenomenological analysis by using the measured ${\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)$ and ${\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)_{\mathrm {high}-q^2}$, for which the hadronic uncertainties are well controlled, and we discuss their impact on various types of 2HDM. A brief discussion of the decays with $\tau$-leptons in the final state is provided too.

1 Introduction

Physical processes driven by the flavor changing neutral currents (FCNC) are forbidden in the Standard Model (SM) at tree level. Since they occur through loops, their measurements offer a low-energy window to the particle content in the loops. In other words, they do not only represent a fine test of validity of the SM, but they also offer an opportunity to look for the effects of physics (particles) beyond the SM (BSM) at low energies. The main obstacle to the accurate comparison between the SM theory and the experimental data lies in the fact that the non-perturbative QCD effects are not under full theoretical control. While the solution to non-perturbative QCD is lacking, in some situations the hadronization effects can be solved by means of numerical simulations of QCD on the lattice (LQCD). Over the past couple of decades we witnessed huge progress in reducing the uncertainties in the LQCD results. Nowadays, an excellent theoretical control of the neutral meson mixing processes promoted those FCNC processes to viability tests of the New Physics (NP) model candidates. Besides the oscillation frequencies of the neutral meson systems, the processes based on $b\rightarrow s$ transitions received a great deal of attention in the particle physics community. While the inclusive and exclusive processes based on the penguin-induced $b\rightarrow s\gamma$ decay have been, and still are, a very significant constraint when building a NP model, the processes based on $b\rightarrow s\ell ^+\ell ^-$ received huge attention because they allow one to access another types of penguin and box diagrams. With the advent of the Large Hadron Collider (LHC) the measurement of ${\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)$ became possible [1] and the result appeared to be somewhat lower than predicted [2]. The spectrum of $d{\mathcal {B}}(B\rightarrow K\mu ^+\mu ^-)/dq^2$ has been measured [3] too and in the range of large $q^2$’s it appears to be larger than predicted [4,5,6,7,8,9]. A full angular analysis of ${\mathcal {B}}(B\rightarrow K^*\mu ^+\mu ^-)$ [3, 10] and ${\mathcal {B}}(B_s\rightarrow \phi \mu ^+\mu ^-)$ [11] revealed discrepancies in several observables with respect to their SM predictions [12,13,14,15,16,17,18,19]. Moreover, the recent measurements of $R_K= {\mathcal {B}}^\prime (B\rightarrow K \mu ^+\mu ^-)/{\mathcal {B}}^\prime (B\rightarrow K e^+e^-)$ [20] and $R_{K^*}= {\mathcal {B}}^\prime (B\rightarrow {K^*} \mu ^+\mu ^-)/{\mathcal {B}}^\prime (B\rightarrow {K^*} e^+e^-)$ [21] were shown to be significantly lower than predicted [22].¹ Those new experimental data helped discarding several NP models and are currently used as constraints in building a NP model.

Simultaneously with the research of FCNC processes, the LHC experiments allowed observing the missing ingredient of the SM, the Higgs boson, the mass of which has been found to be $m_h = 125.09(24)$ GeV [23]. While this was a milestone of the LHC, the pending question of hierarchy of scales remains open and a quest for physics BSM continues. One of the minimalistic approaches to building a model of physics BSM is to extend the Higgs sector by introducing an extra Higgs doublet. Phenomenology in the scenarios with two Higgs doublets appears to be very rich and the associated models are generically called the two Higgs doublet models (2HDM), cf. e.g. [24,25,26,27]. Nowadays the experimental search of the additional Higgs bosons is one of the main goals at LHC, in particular that of the charged Higgs boson [28]. Like in the SM, introducing fermions to the 2HDM context results in a plethora of new parameters. To restrain the number of those parameters and to prevent the appearance the FCNC at tree level it is common to assume a peculiar pattern of Yukawa couplings. To test those assumptions one needs to compare many measured observables with theoretical expressions derived in SM with the extended Higgs sector. In this paper we elaborate a few lessons one can learn from the measured $b\rightarrow s\mu ^+\mu ^-$ processes about 2HDM with a softly broken $\mathbb {Z}_2$ symmetry. In doing so we will use two observables, namely ${\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)$ and ${\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)_{\mathrm {high}-q^2}$, which are very well measured experimentally and for which the theoretical control of the corresponding hadronic uncertainties is established by the LQCD computations [29]. For other observables the theoretical uncertainties are not as well assessed and one might run a risk of interpreting the unknown hadronic uncertainties as signals of physics BSM. We should also emphasize that 2HDM alone cannot explain $R_{K^{(^*)}}^\mathrm {SM}>R_{K^{(^*)}}^\mathrm {exp}$ because only the scalar and pseudoscalar operators can generate lepton flavor universality violation in this framework, which is not enough to accommodate the observed deviation from the SM. A study along the line we are pursuing here has been initiated in Ref. [30] in which the authors computed the Wilson coefficients in the Aligned 2HDM (A2HDM), for the operators relevant to the $B_s\rightarrow \mu ^+\mu ^-$ decay. In this paper we revisit their computation and extend it to include the operators that are needed for the phenomenological analysis of $B\rightarrow K^{(*)}\ell ^+\ell ^-$ and other similar decays. While we broadly agree with the results of Ref. [30], there are a couple of points in which we disagree. We will examine those points, compute the remaining Wilson coefficients and use our results to discuss the phenomenological consequences on the 2HDM scenarios by comparing ${\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {2HDM}$ and ${\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)_{\mathrm {high}-q^2}^\mathrm {2HDM}$ with their experimental values. We will then discuss the consequences on the similar decays with $\tau$-leptons in the final state.

The remainder of this paper is organized as follows: In Sect. 2 we remind the reader of the main general constraints on the spectrum of scalars in 2HDM and perform a scan of parameters by assuming the lowest CP-even Higgs state to be the one measured at LHC. In Sect. 3 we write the low-energy effective theory and present our results for all the Wilson coefficients in Sect. 4. We compare our results with the existing ones (in the limits in which the comparison can be made) in Sect. 5 and elucidate the subtleties related to the matching procedure in the between the full (2HDM) and effective theories in Sect. 6. Phenomenological discussion is made in Sects. 7 and 8. We briefly conclude in Sect. 9.

2 General constraints on 2HDM

In this section we remind the reader of the basic ingredients of 2HDM, enumerate the parameters of the model and list the main general constraints on the spectrum of scalars which are then used to perform a scan of allowed parameters to obtain the allowed ranges of the Higgs masses and couplings.

2.1 2HDM

We consider a general CP-conserving 2HDM with a softly broken $\mathbb {Z}_2$ symmetry. The most general potential can then be written as

$$\begin{aligned} V(\Phi _1,\Phi _2)&= m_{11}^2\Phi _1^\dagger \Phi _1 +m_{22}^2\Phi _2^\dagger \Phi _2 + m_{12}^2(\Phi _1^\dagger \Phi _2 +\Phi _2^\dagger \Phi _1)\nonumber \\&\quad +\dfrac{\lambda _1}{2}(\Phi _1^\dagger \Phi _1)^2+\dfrac{\lambda _2}{2}(\Phi _2^\dagger \Phi _2)^2\nonumber \\&\quad +\lambda _3 \Phi _1^\dagger \Phi _1 \Phi _2^\dagger \Phi _2+\lambda _4 \Phi _1^\dagger \Phi _2 \Phi _2^\dagger \Phi _1\nonumber \\&\quad +\dfrac{\lambda _5}{2}[ (\Phi _1^\dagger \Phi _2)^2+(\Phi _2^\dagger \Phi _1)^2], \end{aligned}$$

(1)

where the term proportional to $m_{12}^2$ accounts for the soft breaking of $\mathbb {Z}_2$.² The scalar doublets $\Phi _a$ ($a=1,2$) can be parameterized as

$$\begin{aligned} \Phi _a(x) = \begin{pmatrix} \phi _a^+(x) \\ \frac{1}{\sqrt{2}}\left[ v_a+\rho _a(x)+i \eta _a(x)\right] \end{pmatrix}, \end{aligned}$$

(2)

with $v_{1,2}\ge 0$ being the vacuum expectation values, satisfying $v^\mathrm {SM}=\sqrt{v_1^2+v_2^2}$; it is already known from experiment that $v^\mathrm {SM}=246.22$ GeV [32]. In the following, for notational simplicity, we will drop the argument of the Higgs fields. Two of the six fields are Goldstone bosons, while the remaining ones are four massive scalars: two CP-even states (h, H), one CP-odd state (A), and one charged Higgs ($H^\pm$). These fields are defined as

$$\begin{aligned}&\begin{pmatrix} \phi _1^+\\ \phi _2^+ \end{pmatrix} =\begin{pmatrix} \cos \beta &{} - \sin \beta \\ \sin \beta &{} \cos \beta \end{pmatrix} \begin{pmatrix} G^+\\ H^+ \end{pmatrix},\nonumber \\&\begin{pmatrix} \eta _1\\ \eta _2 \end{pmatrix} =\begin{pmatrix} \cos \beta &{} - \sin \beta \\ \sin \beta &{} \cos \beta \end{pmatrix} \begin{pmatrix} G^0\\ A \end{pmatrix}, \end{aligned}$$

(3)

and

$$\begin{aligned} \begin{pmatrix} \rho _1\\ \rho _2 \end{pmatrix}&=\begin{pmatrix} \cos \alpha &{} - \sin \alpha \\ \sin \alpha &{} \cos \alpha \end{pmatrix} \begin{pmatrix} H\\ h \end{pmatrix}. \end{aligned}$$

(4)

The mixing angles $\alpha$ and $\beta$ satisfy

$$\begin{aligned}&\tan \beta = \frac{v_2}{v_1},\qquad \nonumber \\&\tan 2\alpha = \dfrac{2(-m_{12}^2+\lambda _{345}v_1 v_2)}{m_{12}^2(v_2/v_1-v_1/v_2)+\lambda _1 v_1^2-\lambda _2 v_2^2},\nonumber \\ \end{aligned}$$

(5)

with $\lambda _{345}\equiv \lambda _3+\lambda _4+\lambda _5$. The masses of the physical scalars can be written in terms of parameters which appear in the potential as

$$\begin{aligned}&m_H^2 = M^2 \sin ^2(\alpha -\beta )+\biggl (\lambda _1 \cos ^2\alpha \cos ^2\beta \nonumber \\&\quad \qquad +\lambda _2 \sin ^2\alpha \sin ^2\beta +\frac{\lambda _{345}}{2}\sin 2\alpha \sin 2\beta \biggr )v^2, \end{aligned}$$

(6)

$$\begin{aligned}&m_h^2 = M^2 \cos ^2(\alpha -\beta )+\biggl (\lambda _1 \sin ^2\alpha \cos ^2\beta \nonumber \\&\quad \qquad +\lambda _2 \cos ^2\alpha \sin ^2\beta -\frac{\lambda _{345}}{2}\sin 2\alpha \sin 2\beta \biggr )v^2, \end{aligned}$$

(7)

$$\begin{aligned}&m_{A}^2 = M^2-\lambda _5 v^2, \end{aligned}$$

(8)

$$\begin{aligned}&m_{H^\pm }^2 = M^2-\frac{\lambda _{4}+\lambda _5}{2} v^2, \end{aligned}$$

(9)

where the $\mathbb {Z}_2$ breaking term is now parameterized via $M^2\equiv \dfrac{m_{12}^2}{\sin \beta \cos \beta }$.

In the Yukawa sector, the $\mathbb {Z}_2$ symmetry becomes particularly important as it prevents the flavor changing processes to appear at tree level [33]. Furthermore it enforces each type of the right-handed fermion to couple to a single Higgs doublet. Four choices are then possible and they are called Type I, II, X (lepton specific) and Z (flipped) 2HDM [26, 27].³ To be more specific, we first write the Yukawa Lagrangian as

$$\begin{aligned} \mathcal {L}_Y&= - \dfrac{\sqrt{2}}{v} H^+ \lbrace \bar{u}~[\zeta _d \, V m_d P_R-\zeta _u \, m_u V P_L]~d +\zeta _\ell \, \bar{\nu } m_\ell P_R \ell \rbrace \nonumber \\&\quad -\,\dfrac{1}{v}\sum _{f,\varphi _i^0\in \{h,H,A\}} \xi _f^{\varphi _i^0} \varphi _i^0 \, [\bar{f} m_f P_R f ]+\mathrm {h.c.,} \end{aligned}$$

(10)

where u and d stand for the up- and down-type quark, $\ell$ is a lepton flavor, f stands for a generic fermion, V for the Cabibbo–Kobayashi–Maskawa (CKM) matrix, and $P_{L,R}= (1\mp \gamma _5)/2$. A specific choice of parameters $\zeta _f$ corresponds to the above-mentioned types of 2HDM, which we also summarize in Table 1. Notice that the couplings $\xi _f^{\varphi _i^0}$ appearing in the neutral Lagrangian part can be mapped onto the charged ones via

$$\begin{aligned} \xi ^h_f&= \sin (\beta -\alpha ) +\cos (\beta -\alpha ) \zeta _f, \nonumber \\ \xi ^H_f&= \cos (\beta -\alpha ) -\sin (\beta -\alpha ) \zeta _f, \nonumber \\ \xi ^A_{u}&=-i\zeta _u,\qquad \xi ^A_{d,\ell }=i \zeta _{d,\ell }. \end{aligned}$$

(11)

Table 1 Couplings $\zeta _f$ in various types of 2HDM

2.2 General constraints and scan of parameters

To perform a thorough scan of scalars in a general 2HDM we use the general constraints summarized below.

• Stability: To ensure that the scalar potential is bounded from below, the quartic couplings should satisfy the relations [25]

  $$\begin{aligned}&\lambda _{1,2}>0,\qquad \lambda _3>-(\lambda _1 \lambda _2)^{1/2},\qquad \mathrm {and} \qquad \nonumber \\&\lambda _3+\lambda _4- \vert \lambda _5 \vert >-(\lambda _1 \lambda _2)^{1/2}. \end{aligned}$$

  (12)

  Furthermore, the stability of the electroweak vacuum implies that

  $$\begin{aligned} m_{11}^2+\dfrac{\lambda _1 v_1^2}{2}+\dfrac{\lambda _3 v_2^2}{2}&= \frac{v_2}{v_1} \left[ m_{12}^2 - (\lambda _4+\lambda _5)\dfrac{v_1 v_2}{2}\right] , \nonumber \\\end{aligned}$$

  (13)
  $$\begin{aligned} m_{22}^2+\dfrac{\lambda _2 v_2^2}{2}+\dfrac{\lambda _3 v_1^2}{2}&= \frac{v_1}{v_2} \left[ m_{12}^2 - (\lambda _4+\lambda _5)\dfrac{v_1 v_2}{2}\right] ,\nonumber \\ \end{aligned}$$

  (14)

  which then allows us to express $m_{11}^2$ and $m_{22}^2$ in terms of the soft $\mathbb {Z}_2$ breaking term $m_{12}^2$ and the quartic couplings $\lambda _{1-5}$. These constraints should be combined with the necessary and sufficient condition that the minimum developed at $(v_1,v_2)$ is global [34]:

  $$\begin{aligned} m_{12}^2 \left( m_{11}^2-m_{22}^2 \sqrt{\lambda _1/\lambda _2} \right) \left( \tan \beta - \root 4 \of {\lambda _1/\lambda _2}\right) >0. \end{aligned}$$

  (15)

• Perturbative unitarity: An important constraint on the spectrum of scalars within 2HDM stems from the unitarity requirement of the S-wave component of the scalar scattering amplitudes. That condition implies the following inequalities [35, 36]:

  $$\begin{aligned} |a_\pm |, |b_\pm |, |c_\pm |, |f_\pm |, |e_{1,2}|, |f_1|, |p_1| < 8 \pi , \end{aligned}$$

  (16)

  where

  $$\begin{aligned} a_\pm&= \dfrac{3}{2}(\lambda _1+\lambda _2)\pm \sqrt{\dfrac{9}{4}(\lambda _1-\lambda _2)^2+(2\lambda _3+\lambda _4)^2},\nonumber \\ b_\pm&= \dfrac{1}{2}(\lambda _1+\lambda _2)\pm \dfrac{1}{2} \sqrt{(\lambda _1-\lambda _2)^2+4\lambda _4^2},\nonumber \\ c_\pm&= \dfrac{1}{2}(\lambda _1+\lambda _2)\pm \dfrac{1}{2} \sqrt{(\lambda _1-\lambda _2)^2+4\lambda _5^2},\nonumber \\ e_1&= \lambda _3 + 2 \lambda _4 -3\lambda _5,\qquad e_2 = \lambda _3-\lambda _5,\nonumber \\ f_+&= \lambda _3+2 \lambda _4+3\lambda _5, \qquad f_- =\lambda _3+\lambda _5,\nonumber \\ f_1&= \lambda _3+\lambda _4, \qquad p_1 = \lambda _3-\lambda _4. \end{aligned}$$

  (17)

• Electroweak precision tests: Finally, the additional scalars contribute to the gauge boson vacuum polarization. As a result, the electroweak precision data provide important constraint. In particular the T parameter bounds the mass splitting between $m_H$ and $m_{H^\pm }$ in the scenario in which h is identified with the SM-like Higgs, cf. Ref. [37] for example. The general expressions for the parameters S, T and U in 2HDM can be found in Ref. [38]. To derive the bounds on the scalar spectrum we consider the following values and the corresponding correlation matrix [39]:

  $$\begin{aligned}&\begin{aligned}&\Delta S^\mathrm{SM} = 0.05\pm 0.11, \\&\Delta T^\mathrm{SM} = 0.09\pm 0.13, \\&\Delta U^\mathrm{SM} = 0.01\pm 0.11, \\ \end{aligned} \qquad \qquad \nonumber \\&\quad \mathrm {corr} = \left( \begin{array}{c@{\quad }c@{\quad }c} 1 &{} 0.90 &{} -0.59 \\ 0.90 &{} 1 &{} -0.83 \\ -0.59 &{} -0.83 &{} 1 \end{array}\right) . \end{aligned}$$

  (18)

  The $\chi ^2$ function is then expressed as

  $$\begin{aligned} \chi ^2= \sum _{i,j}(X_i - X_i^\mathrm{SM})(\sigma ^2)_{ij}^{-1}(X_j - X_j^\mathrm{SM}), \end{aligned}$$

  (19)

  where the vectors of central values and uncertainties are denoted $X=(\Delta S, \Delta T, \Delta U)$ and $\sigma =(0.11,0.13,0.11)$, while the elements of the covariance matrix are obtained via $\sigma _{ij}^2\equiv \sigma _i \mathrm {corr}_{ij} \sigma _j$.

Fig. 1

Results of the scan described in the text

As mentioned above, we identify the lightest CP-even state h with the SM-like scalar observed at the LHC with mass $m_h=125.09(24)$ GeV [32]. To forbid the dangerous decays $h\rightarrow A A$ which could over-saturate the total width of h ($\simeq \Gamma _h^\mathrm {SM}$), we assume that $m_A> m_h/2$. Moreover, we impose the alignment condition $|\cos (\beta -\alpha )|\le 0.3$, in order to ensure that the couplings of h to $V=W,Z$ remain consistent with the values measured so far, which appear to be in good agreement with the SM predictions [40]. The above-mentioned constraints are then imposed onto a set of randomly generated points in the intervals:

$$\begin{aligned}&\tan \beta \in (0.2,50),\qquad \alpha \in \left( -\frac{\pi }{2},\frac{\pi }{2}\right) , \qquad \left| M^2\right| \le (1.2~\mathrm {TeV})^2,\nonumber \\&m_{H^\pm }\in (m_W, 1.2~\mathrm {TeV}),\qquad m_{H}\in (m_h, 1.2~\mathrm {TeV}),\nonumber \\&\qquad m_{A}\in \left( m_h/2, 1.2~\mathrm {TeV}\right) . \end{aligned}$$

(20)

A scan of parameters consistent with the constraints listed above favors the moderate and small values of $\tan \beta \in (0.2,15]$. To see that the larger values of $\tan \beta$ cannot be discarded it is sufficient to examine Eq. (6) in the alignment limit. For that reason, and in addition to the free scan, we perform a second scan with $m_H\approx |M|$, which helps us probing higher values of $\tan \beta$, and we then combine results of both scans. The combined results are shown in Fig. 1 in two planes, $(\tan \beta , m_{H^\pm })$ and $(m_A, m_{H^\pm })$. From the right panel of Fig. 1 we observe that the additional scalars become mass degenerate in the decoupling region ($M^2 \gg v^2$), as it can be easily deduced from Eqs. (6)–(9). We should also emphasize that the results of our scans agree with what has been previously reported in the literature, cf. [41,42,43,44].

In Sect. 8 we will confront the points allowed by our scan with the experimental measurements of exclusive $b\rightarrow s$ decays.

3 Effective Hamiltonian

The most general effective Hamiltonian describing the $b\rightarrow s \ell \ell$ transitions, made of dimension six operators, is given by [45]

$$\begin{aligned}&\mathcal {H}_\mathrm {eff} = - \frac{4 G_F}{\sqrt{2}}V_{tb}V_{ts}^*\nonumber \\&\qquad \quad \times \sum _{i}(C_i(\mu )\mathcal {O}_i(\mu )+C_i^\prime (\mu )\mathcal {O}_i^\prime (\mu )) + \mathrm {h.c.}, \end{aligned}$$

(21)

where

$$\begin{aligned} \mathcal {O}_9&= \frac{e^2}{(4\pi )^2}(\bar{s}\gamma _\mu P_L b)(\bar{\ell }\gamma ^\mu \ell ), \qquad \mathcal {O}_S = \frac{e^2}{(4\pi )^2}(\bar{s} P_R b)(\bar{\ell } \ell ), \end{aligned}$$

(22)

$$\begin{aligned} \mathcal {O}_{10}&= \frac{e^2}{(4\pi )^2}(\bar{s}\gamma _\mu P_L b)(\bar{\ell }\gamma ^\mu \gamma _5 \ell ) ,\qquad \mathcal {O}_P = \frac{e^2}{(4\pi )^2}(\bar{s} P_R b)(\bar{\ell }\gamma _5 \ell ), \end{aligned}$$

(23)

$$\begin{aligned} \mathcal {O}_T&= \frac{e^2}{(4\pi )^2}(\bar{s}\sigma _{\mu \nu } b)(\bar{\ell }\sigma ^{\mu \nu } \ell ), \qquad \mathcal {O}_{T5} = \frac{e^2}{(4\pi )^2}(\bar{s}\sigma _{\mu \nu } b)(\bar{\ell }\sigma ^{\mu \nu } \gamma _5\ell ), \end{aligned}$$

(24)

and $\mathcal {O}_7=e/(4\pi )^2 m_b (\bar{s} \sigma _{\mu \nu }P_R b)F^{\mu \nu }$ is the electromagnetic penguin operator. The operators with a flipped chirality, $\mathcal {O}^\prime _{7,9,10,S,P}$, are obtained from $\mathcal {O}_{7,9,10,S,P}$ by replacing $P_L \leftrightarrow P_R$ in the quark current.

Fig. 2

Photon penguin diagrams generated by the charged Higgs bosons

The dimension six operators appearing in Eq. (21) are sufficient to match the one-loop amplitude when the external fermion momenta are neglected. This, however, is not true if the computation is made with external momenta different from zero which is, in general, necessary when dealing with (pseudo-)scalar operators. For example, in order to get a correct expression for the Wilson coefficient $C_P$ one needs to consider the external momenta, which then leads to $C_P\propto m_\ell m_b/m_W^2$, c.f. [30]. We therefore need to select all situations in which one can obtain the terms of the form $(m_\ell m_b/m_W^2) \mathcal {O}_P$, such as

(25)

which can obviously be reduced to $(m_\ell m_b/m_W^2) \mathcal {O}_P$. As an example,

(26)

We should stress at this point that in the course of our computation the contributions proportional to the chirality flipped operator $\mathcal {O}_{S,P}^\prime$ are suppressed by a factor of $m_s$ which is why they are readily neglected. A complication arises when encountering the operators with insertion of in the leptonic current, with the convention $b(p_b)\rightarrow s(p_s) \ell ^-(p_-)\ell ^+(p_+)$, where we also use $q=p_b-p_s=p_++p_-$. A way to deal with that, adopted in Ref. [30], consists in setting $p_s=0$, so that , and in this way one can again, like in the previous example, use the equations of motion. That way to deal with the problem in hands, however, leads to an incomplete expression for $C_P$, for example. If, instead, one keeps all the momenta nonzero, we get a complete result. At this point we just emphasize that the matching should be performed by keeping all the external momenta different from zero and the contributions stemming from dimension-seven operators can be neglected at the very end of computation. We further elucidate this problem in Sect. 6 where we also propose a general framework for the appropriate matching between the full and effective theories in a case in which the (pseudo-)scalar bosons are explicitly taken into account. Before closing this section we could also argue that a term $(m_\ell m_b/m_W^2) \mathcal {O}_P$ could be obtained from the dimension-eight operators, such as

(27)

Such terms, however, never appear in our calculations and we will limit our discussion to the dimension-seven operators only.

4 Wilson coefficients

After unambiguously matching the full with the effective theories we obtain the one-loop expressions for the Wilson coefficients generated by the additional scalar particles. We summarize our results in this section. For clarity we will write them as

$$\begin{aligned} C_{7}&= C_7^{\mathrm {NP}\,,\gamma } , \end{aligned}$$

(28)

$$\begin{aligned} C_{9}&= C_9^{\mathrm {NP}\,,\gamma }+C_9^{\mathrm {NP},\,Z}, \end{aligned}$$

(29)

$$\begin{aligned} C_{10}&= C_{10}^{\mathrm {NP},\,Z}, \end{aligned}$$

(30)

$$\begin{aligned} C_{P}&= C_{P}^{\mathrm {NP},\,\mathrm {box}}+C_{P}^{\mathrm {NP},\,Z}+C_{P}^{\mathrm {NP},\,A} \end{aligned}$$

(31)

$$\begin{aligned} C_{S}&= C_{S}^{\mathrm {NP},\,\mathrm {box}}+C_S^{\mathrm {NP},\,h}+C_S^{\mathrm {NP},\,H} \end{aligned}$$

(32)

where the superscripts denote the types of diagrams that contributes to a given Wilson coefficient, namely, the box diagrams, the $\gamma ,Z$-penguins and the (pseudo-)scalar penguins. These coefficients should be added to the (effective) ones obtained in the SM: $C_{7}=-0.304$, $C_9=4.211$, $C_{10}=-4.103$, and $C_{S,P}\simeq 0$ [46].⁴

Henceforth, we neglect the s-quark mass and give all our results in the unitary gauge. To check the consistency of our formulas, we also performed the computation in the Feynman gauge. In the remainder of this section we present our resulting expressions for each of the coefficients appearing in Eqs. (30)–(32). We use the standard notation,

$$\begin{aligned} x_q = \dfrac{m_q^2}{m_W^2},\qquad \quad x_{H^\pm } = \dfrac{m_{H^\pm }^2}{m_W^2},\qquad \quad x_{\varphi _i^0}=\dfrac{m_{\varphi _i^0}^2}{m_W^2}, \end{aligned}$$

(33)

where $q\in \{ b,t\}$, and $\varphi _i^0 \in \lbrace h,H,A \rbrace$.

Fig. 3

Z penguin diagrams generated by the additional scalars

Fig. 4

Box diagrams generated by the additional scalars

4.1 $\gamma$-penguins in 2HDM

The $\gamma$-penguin diagrams induced by the charged Higgs are shown in Fig. 2. The off-shell and on-shell contributions can be matched onto the Wilson coefficients $C_7$ and $C_9$, respectively, we obtain

$$\begin{aligned} \begin{aligned} C_7^{\mathrm {NP}, \gamma } =&-|\zeta _u|^2 \frac{x_t}{72}\Bigg [\frac{7 x_{H^\pm }^2-5 x_{H^\pm }x_t-8 x_t^2}{ (x_{H^\pm }-x_{t})^3}\\&\quad +\dfrac{6 x_{H^\pm }x_t(3x_t-2 x_{H^\pm })}{(x_{H^\pm }-x_{t})^4}\log \left( \frac{x_{H^\pm }}{x_{t}} \right) \Bigg ]\\&\quad -\zeta _u^*\zeta _d \frac{x_t}{12} \Bigg [ \frac{3 x_{H^\pm }-5 x_t}{(x_{t}-x_{H^\pm })^2}+\dfrac{2 x_{H^\pm }(3x_{t}-2x_{H^\pm })}{(x_{t}-x_{H^\pm })^3}\\&\qquad \quad \times \log \left( \frac{x_{t}}{x_{H^\pm }} \right) \Bigg ] \end{aligned} \end{aligned}$$

(34)

and

$$\begin{aligned} \begin{aligned} C_9^{\mathrm {NP}, \gamma }&= |\zeta _u|^2 \frac{x_t}{108 }\left[ \frac{38 x_{H^\pm }^2-79 x_{H^\pm }x_t+47 x_t^2}{(x_{H^\pm }-x_{t})^3}\right. \\&\quad \left. -\dfrac{6(4 x_{H^\pm }^3-6 x_{H^\pm }^2 x_t+3 x_t^3)}{(x_{H^\pm }-x_{t})^4}\log \left( \frac{x_{H^\pm }}{x_{t}}\right) \right] \\&\quad +\zeta _u^*\zeta _d\dfrac{x_t x_b}{108}\left[ \frac{-37 x_{H^\pm }^2+8x_{H^\pm }x_t+53 x_t^2}{(x_{H^\pm }-x_t)^4}\right. \\&\quad \left. +\frac{6(2 x_{H^\pm }^3+6 x_{H^\pm }^2 x_t-9 x_{H^\pm }x_t^2-3 x_t^3)}{(x_{H^\pm }-x_t)^5}\log \left( \frac{x_{H^\pm }}{x_t}\right) \right] . \end{aligned} \end{aligned}$$

(35)

The dominant terms in both $C_7^{\mathrm {NP}, \gamma }$ and $C_9^{\mathrm {NP}, \gamma }$ come from the top quark contribution and are proportional to $|\zeta _u|^2$. The term proportional to $\zeta _u^*\zeta _d$ in $C_9^{\mathrm {NP}, \gamma }$ is suppressed by $m_b^2$, thus it indeed is subdominant.

4.2 Z-penguins in 2HDM

The Z-penguin diagrams contribute significantly to the Wilson coefficients $C_P$, $C_9$ and $C_{10}$ through the diagrams shown in Fig. 3. The leading order expressions for $C_9$ and $C_{10}$ read

$$\begin{aligned} C_{9}^{\mathrm {NP},Z}&= C_{10}^{\mathrm {NP},Z}(-1+4\sin ^2 \theta _{W}), \end{aligned}$$

(36)

$$\begin{aligned} C_{10}^{\mathrm {NP},Z}&= |\zeta _u|^2 \frac{x_t^2}{8 \sin ^2\theta _W }\Bigg [\frac{1}{x_{H^\pm }-x_{t}}-\frac{x_{H^\pm }}{(x_{H^\pm }-x_{t})^2}\nonumber \\&\qquad \qquad \qquad \qquad \qquad \times \log \left( \frac{x_{H^\pm }}{x_{t}}\right) \Bigg ]\nonumber \\&+\zeta _u^*\zeta _d \frac{x_t x_b}{16 \sin ^2\theta _W}\Bigg [\frac{x_{H^\pm }+x_{t}}{(x_{H^\pm }-x_{t})^2}-\frac{2 x_t x_{H^\pm }}{(x_{H^\pm }-x_{t})^3}\nonumber \\&\qquad \qquad \qquad \qquad \qquad \times \log \left( \frac{x_{H^\pm }}{x_{t}}\right) \Bigg ]. \end{aligned}$$

(37)

Similarly, for $C_P$ we obtain

$$\begin{aligned} C_P^{\mathrm {NP},Z}&= \zeta _u^*\zeta _d \dfrac{\sqrt{x_b x_\ell } \,x_t}{16 \sin ^2\theta _W}\Bigg [\frac{x_t-3 x_{H^\pm }}{(x_{H^\pm }-x_t)^2} +\frac{2 x_{H^\pm }^2}{(x_{H^\pm }-x_t)^3}\log \left( \frac{x_{H^\pm }}{x_t}\right) \Bigg ] \nonumber \\&\quad +|\zeta _u|^2\frac{\sqrt{x_b x_\ell } \,x_t}{216}\nonumber \\&\quad \Biggl \lbrace \frac{38 x_{H^\pm }^2+54 x_{H^\pm }^2 x_t-79 x_{H^\pm }x_t-108 x_{H+} x_t^2+47 x_t^2+54 x_t^3}{(x_{H^\pm }-x_t)^3 } \nonumber \\&\quad -\dfrac{6(4 x_{H^\pm }^3+9 x_{H^\pm }^3 x_t-6 x_{H^\pm }^2 x_t-18 x_{H^\pm }^2x_t^2+9 x_{H^\pm }x_t^3+3 x_t^3)}{(x_{H^\pm }-x_t)^4}\nonumber \\&\quad \times \log \left( \frac{x_{H^\pm }}{x_t}\right) -\dfrac{3}{2\sin ^2{\theta _W}}\nonumber \\&\quad \times \Biggl [\dfrac{2 x_{H^\pm }^2 +36 x_{H^\pm }^2 x_t- 7 x_{H^\pm }x_t-72 x_{H^\pm } x_t^2+11 x_t^2+36 x_t^3}{(x_{H^\pm }-x_t)^3}\nonumber \\&\quad -\dfrac{6 x_t(6 x_{H^\pm }^3-12 x_{H^\pm }^2x_t+6 x_{H^\pm }x_t^2+x_t^2}{(x_{H^\pm }-x_t)^4} \log \left( \frac{x_{H^\pm }}{x_t}\right) \Biggr ]\Biggr \rbrace . \end{aligned}$$

(38)

4.3 Charged Higgs boxes in 2HDM

The box diagrams, peculiar for 2HDM, are drawn in Fig. 4. At low energy they contribute to the Wilson coefficients $C_{S}$ and $C_{P}$; we have

$$\begin{aligned} \begin{aligned} C_S^\mathrm {NP,\,box}&= \dfrac{\sqrt{x_\ell x_b}\,x_t}{8(x_{H^\pm }-x_t)\sin ^2\theta _W}\\&\qquad \times \Biggl \lbrace \zeta _\ell \zeta _u^*\left( \frac{x_t}{x_t-1}\log x_t-\frac{x_{H^\pm }}{x_{H^\pm }-1}\log x_{H^\pm }\right) \\&\quad +\zeta _u \zeta _\ell ^*\left[ 1-\frac{x_{H^\pm }-x_t^2}{(x_{H^\pm }-x_t)(x_t-1)}\log x_t\right. \\&\quad \left. - \frac{x_{H^\pm }(x_t-1)}{(x_{H^\pm }-x_t)(x_{H^\pm }-1)}\log x_{H^\pm } \right] \\&\quad +2\zeta _d\zeta _\ell ^*\log \left( \dfrac{x_t}{x_{H^\pm }}\right) \Biggr \rbrace \end{aligned} \end{aligned}$$

(39)

and

$$\begin{aligned} C_P^\mathrm {NP,\,box}&= \dfrac{\sqrt{x_\ell x_b}\,x_t }{8(x_{H^\pm }-x_t)\sin ^2\theta _W}\nonumber \\&\quad \times \Biggl \lbrace \zeta _\ell \zeta _u^*\left( \frac{x_t}{x_t-1}\log x_t-\frac{x_{H^\pm }}{x_{H^\pm }-1}\log x_{H^\pm }\right) \nonumber \\&\quad -\zeta _u \zeta _\ell ^*\left[ 1-\frac{x_{H^\pm }-x_t^2}{(x_{H^\pm }-x_t)(x_t-1)}\log x_t \right. \nonumber \\&\quad \left. - \frac{x_{H^\pm }(x_t-1)}{(x_{H^\pm }-x_t)(x_{H^\pm }-1)}\log x_{H^\pm } \right] \nonumber \\&\quad -2\zeta _d\zeta _\ell ^*\log \left( \frac{x_t}{x_{H^\pm }}\right) \Biggr \rbrace . \end{aligned}$$

(40)

In addition to $C_{S,P}^\mathrm {NP,\,box}$, the tensor and (axial-)vector operators receive contributions but suppressed by the lepton mass, i.e. by $x_\ell =m_\ell ^2/m_W^2$. These coefficients are negligible even for decays with $\tau$’s in the final state as it can be verified by using the expressions we provide in Appendix C.2.

4.4 Scalar penguins in 2HDM

We now turn to the effective coefficients $C_P^{\mathrm {NP},\, A}$, $C_S^{\mathrm {NP},\, h}$ and $C_S^{\mathrm {NP},\, H}$, generated by the scalar penguin diagrams shown in Fig. 5. We recall that the total ultraviolet divergence coming from these diagrams is proportional to the factor $(1+\zeta _u\zeta _d)(\zeta _u-\zeta _d)$, which vanishes due to the $\mathbb {Z}_2$ symmetry (cf. Table 1).⁵

The penguins with the CP-odd Higgs give rise to

$$\begin{aligned} C_P^{\mathrm {NP},\, A}&= -\dfrac{\sqrt{x_\ell x_b}}{\sin ^2\theta _W}\dfrac{\zeta _\ell x_t}{2 x_A }\nonumber \\&\quad \times \Bigg \lbrace \dfrac{\zeta _u^3 x_t}{2} \Bigg [ \dfrac{1}{x_{H^\pm }-x_t}-\dfrac{x_{H^\pm }}{(x_{H^\pm }-x_t)^2}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ]\nonumber \\&\quad +\frac{\zeta _u}{4}\Bigg [ -\dfrac{3 x_{H^\pm }x_t-6 x_{H^\pm }-2 x_t^2+5x_t}{(x_t-1)(x_{H^\pm }-x_t)}\nonumber \\&\quad +\dfrac{x_{H^\pm }(x_{H^\pm }^2-7 x_{H^\pm }+6 x_t)}{(x_{H^\pm }-x_t)^2(x_{H^\pm }-1)}\log x_{H^\pm }\nonumber \\&\quad -\dfrac{x_{H^\pm }^2(x_t^2-2 x_t+4)+3 x_t^2(2x_t-2 x_{H^\pm }-1)}{(x_{H^\pm }-x_t)^2(x_t-1)^2} \log x_t\Bigg ]\Bigg \rbrace , \end{aligned}$$

(41)

where we used $\zeta _f \in \mathbb {R}$; $(1+\zeta _u\zeta _d)(\zeta _u-\zeta _d)=0$. Similarly, the penguins with the CP-even Higgs lead to

$$\begin{aligned} C_S^{\mathrm {NP},\, h}&= \dfrac{\sqrt{x_\ell x_b}}{\sin ^2 \theta _W}\dfrac{x_t}{2 x_h }\left[ \sin (\beta -\alpha )+\cos (\beta -\alpha )\zeta _\ell \right] \nonumber \\&\quad \times \Bigg [g_1\sin (\beta -\alpha )+g_2\cos (\beta -\alpha )\nonumber \\&\quad -g_0 \frac{2v^2}{m_W^2}\lambda _{H^+ H^-}^h \Bigg ],\nonumber \\ C_S^{\mathrm {NP},\, H}&= \dfrac{\sqrt{x_\ell x_b}}{\sin ^2 \theta _W}\dfrac{x_t}{2 x_H }\left[ \cos (\beta -\alpha )-\sin (\beta -\alpha )\zeta _\ell \right] \nonumber \\&\quad \times \Bigg [g_1\cos (\beta -\alpha )-g_2\sin (\beta -\alpha )\nonumber \\&\quad -g_0\frac{2v^2}{m_W^2}\lambda _{H^+ H^-}^H\Bigg ], \end{aligned}$$

(42)

where $\lambda _{H^+H^-}^{\varphi _i^0}$ are the trilinear couplings defined in Appendix B. The functions $g_{0,1,2}$ are given in Appendix C along with the amplitudes generated by each of the diagrams shown in Fig. 5.

Fig. 5

Higgs penguin diagrams generated by the additional scalars

5 Comparison with other computations

In this section we compare our Wilson coefficients with the results obtained in previous studies. Before doing so we should emphasize the novelties of the present work:

• (i) The result for $C_9$ in a general 2HDM with a $\mathbb {Z}_2$ symmetry is new.

• (ii) The subleading terms $\mathcal {O}(m_{b})$ to $C_{9,10}$ have been neglected in the previous computations, and they are included here.

• (iii) We provided an independent computation of the coefficients $C_S$ and $C_P$, and elucidate inconsistencies present in Ref. [30], cf. Sect. 6, where we propose a general prescription for matching procedure when the external momenta are not neglected.

The effective coefficients $C_S$ and $C_P$, in the context of Type II 2HDM, were first computed in Refs. [47,48,49,50,51,52]. In these papers $\tan \beta$ was assumed to be very large, which considerably simplifies the computation because in that case only the box diagrams give significant contributions. We agree with these results if we keep only the leading terms in $\tan \beta$ in our expressions, namely

$$\begin{aligned} C_P=-C_S&\simeq \tan ^2\beta \dfrac{\sqrt{x_\ell x_b}}{4\sin ^2\theta _W}\dfrac{x_t}{x_{H^\pm }-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) . \end{aligned}$$

(43)

Along the same lines, the leading order QCD corrections to the same coefficients were included in Ref. [53]. Recently, the computation of $C_S$ and $C_P$ was extended to the context of a general A2HDM, which comprises all four types of 2HDM with $\mathbb {Z}_2$ symmetry discussed here but without the usual assumption of large $\tan \beta$ [30]. We agree with their general results, except for the expression for $C_P^{\mathrm {NP},\, Z}$ which differs from the one reported in the present paper. The disagreement comes from the fact that the authors of Ref. [30] worked with the assumption $p_s=0$, which appears not to be fully appropriate.⁶ By keeping $p_s\ne 0$ one realizes that the computation of Z-penguin leads to two independent terms, one proportional to $p_H=p_b+p_s$ and the other to $q=p_b-p_s$. By using the equations of motion, $C_{P,S}$ correctly receive contributions from the terms proportional to q, but not from those proportional to $p_H$. With $p_s=0$ only one invariant appears, because $p_H\equiv q$, and thus the resulting $C_{P,S}$ also receive spurious contributions from $p_H$.

Regarding the other Wilson coefficients, the first computations of $C_7$ for a general 2HDM have been performed in Ref. [54], then in Refs. [55,56,57] where the leading and subleading QCD corrections were included too. Our results are consistent with those, as well as with the expression for $C_{10}$ presented in Ref. [51] and more recently in Ref. [30]. The only difference with respect to those results is that we include the subleading terms in $m_b$.

6 Matching procedure

In this section we discuss in more detail the matching of the one-loop amplitudes when the nonzero external momenta are considered. We stress once again that keeping external momenta nonzero is necessary to obtain the correct values for the Wilson coefficients $C_{S,P}$. As we mentioned in Sect. 3 the insertion of external momenta result in dimension-seven operators which can be simplified by using equations of motion, except in the cases when the lepton momenta are to be contracted with the quark current and/or the quark momenta to be contracted with the lepton current. The amplitudes which need a special treatment, which give rise to the terms $\propto m_\ell m_b/m_W^2$, are

(44)

where $i,j=L,R$ and $s,b,\ell$ are the fermion spinors. Note again that our convention is $b(p_b)\rightarrow s(p_s) \ell ^-(p_-)\ell ^+(p_+)$, and $q=p_b-p_s=p_++p_-$. In our calculation, specifically, the amplitude $\mathcal {A}_{ij}^q$ appeared in the computation of the Z-penguin diagrams, while $\mathcal {A}^\ell _{ij}$ and $\mathcal {A}_{ij}^{Vq}$ are encountered when computing the box diagrams.

In order to keep our discussion general, we first extend the Hamiltonian (21) and include the following operators:

$$\begin{aligned} \mathcal {H}_\mathrm {eff}^{\prime }&= - \frac{4 G_F}{\sqrt{2}}V_{tb}V_{ts}^*\sum _{i,j=L,R}(C_{ij}^{\mathcal {T}\ell }(\mu )\mathcal {O}_{ij}^{\mathcal {T}\ell }(\mu )\nonumber \\&\quad +C_{ij}^{\mathcal {T}q}(\mu )\mathcal {O}_{ij}^{\mathcal {T}q}(\mu ))+ \mathrm {h.c.}, \end{aligned}$$

(45)

where

$$\begin{aligned} \begin{aligned} \mathcal {O}_{ij}^{\mathcal {T}\ell }&=\frac{e^2}{(4\pi )^2}\frac{1}{m_W}(\bar{s}\gamma ^\mu P_i b)\partial ^\nu (\bar{\ell }\sigma _{\mu \nu } P_j\ell ),\\ \mathcal {O}_{ij}^{\mathcal {T}q}&=-\frac{e^2}{(4\pi )^2}\frac{1}{m_W}{\partial ^\nu }(\bar{s}\sigma _{\mu \nu } P_i b)(\bar{\ell }\gamma ^\mu P_j \ell ), \end{aligned} \end{aligned}$$

(46)

with $i,j=L,R$.⁷ We reiterate that even though these operators are suppressed by $1/m_W$, they are necessary to unambiguously match the loop induced amplitudes with the effective field theory. The above choice of the basis of dimension-seven operators is convenient since they do not contribute to $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)$, while for the other decays their hadronic matrix elements are easy to calculate.

By using the Fierz rearrangement and by applying the field equations, the amplitudes (44) are reduced to

$$\begin{aligned} \mathcal {A}^\ell _{LL}&\leftrightarrow - \mathcal {O}_{LL}^{\mathcal {T}\ell }+\mathcal {O}_9\frac{m_\ell }{m_W}, \end{aligned}$$

(47)

$$\begin{aligned} \mathcal {A}^\ell _{LR}&\leftrightarrow - \mathcal {O}_{LR}^{\mathcal {T}\ell }+\mathcal {O}_9\frac{m_\ell }{m_W}, \end{aligned}$$

(48)

$$\begin{aligned} \mathcal {A}^{V\ell }_{LL}&\leftrightarrow -\mathcal {O}_{LL}^{\mathcal {T}q}+ \left( \mathcal {O}_S^\prime -\dfrac{\mathcal {O}_T-\mathcal {O}_{T5}}{4}\right) \dfrac{m_\ell }{m_W}, \end{aligned}$$

(49)

$$\begin{aligned} \mathcal {A}^{V\ell }_{LR}&\leftrightarrow \mathcal {O}_{LR}^{\mathcal {T}q}+ \left( \mathcal {O}_S^\prime +\dfrac{\mathcal {O}_T-\mathcal {O}_{T5}}{4}\right) \dfrac{m_\ell }{m_W}, \end{aligned}$$

(50)

$$\begin{aligned} \mathcal {A}^q_{LL}&\leftrightarrow \mathcal {O}_{LL}^{\mathcal {T}q}+\dfrac{\mathcal {O}_9^\prime -\mathcal {O}_{10}^\prime }{2}\dfrac{m_b}{m_W}+\dfrac{\mathcal {O}_9-\mathcal {O}_{10}}{2}\dfrac{m_s}{m_W}, \end{aligned}$$

(51)

$$\begin{aligned} \mathcal {A}^q_{LR}&\leftrightarrow \mathcal {O}_{LR}^{\mathcal {T}q}+\dfrac{\mathcal {O}_9^\prime +\mathcal {O}_{10}^\prime }{2}\dfrac{m_b}{m_W}+\dfrac{\mathcal {O}_9+\mathcal {O}_{10}}{2}\dfrac{m_s}{m_W}, \end{aligned}$$

(52)

$$\begin{aligned} \mathcal {A}^{Vq}_{LL}&\leftrightarrow \mathcal {O}_{LL}^{\mathcal {T}\ell }+\dfrac{\mathcal {O}_S-\mathcal {O}_P}{2}\dfrac{m_b}{m_W}\nonumber \\&\quad +\left( \mathcal {O}_S^\prime -\mathcal {O}_P^\prime - \dfrac{\mathcal {O}_T-\mathcal {O}_{T5}}{2}\right) \dfrac{m_s}{2 m_W}, \end{aligned}$$

(53)

$$\begin{aligned} \mathcal {A}^{Vq}_{LR}&\leftrightarrow -\mathcal {O}_{LR}^{\mathcal {T}\ell }+\dfrac{\mathcal {O}_S^\prime +\mathcal {O}_P^\prime }{2}\dfrac{m_s}{m_W}\nonumber \\&\quad +\left( \mathcal {O}_S+\mathcal {O}_P + \dfrac{\mathcal {O}_T+\mathcal {O}_{T5}}{2}\right) \dfrac{m_b}{2 m_W}. \end{aligned}$$

(54)

To remain completely general, in the above equations we also kept the lepton mass and the mass of s-quark different from zero. As an example we show the validity of Eq. (49). Using $p_--p_+=2p_--q$, and by the multiple use of the field equations, we can write

(55)

By applying the Fierz identity once again, we arrive at

$$\begin{aligned} \mathcal {A}^{V\ell }_{LL}&{\mathop {\rightarrow }\limits ^{\text {Fierz}}} \frac{m_\ell }{m_W} \left( \mathcal {O}_S^\prime - { \mathcal {O}_T - \mathcal {O}_{T5}\over 4}\right) - \mathcal {O}_{LL}^{\mathcal {T}q}. \end{aligned}$$

(56)

Clearly, for the appropriate matching of these amplitudes to the effective theory, the operators appearing in Eq. (21) are not enough and the extended basis given in Eq. (45) is necessary. Once the matching is performed, the operators from Eq. (21) could be neglected since they are $1/m_W$ suppressed with respect to the dominant (dimension six) ones.

This delicate point can then be verified explicitly by computing the Wilson coefficients $C_{RL}^{\mathcal {T}q}$ and $C_{RR}^{\mathcal {T}q}$ which come from the Z-penguin diagrams and the coefficients $C^{\mathcal {T}\ell }_{LL}=(C^{\mathcal {T}\ell }_{LR})^*$ generated by the box diagrams. Their explicit expression is given in Appendix C.1.

We can now easily understand the source of our disagreement with Ref. [30]. If one sets $p_s=0$ in $\mathcal {A}_{RR}^q$ of Eq. (44), then just like in Ref. [30] one could write which, by means of the equations of motion, yields

$$\begin{aligned} \mathcal {A}_{RR}^q = \frac{m_\ell }{m_W} \dfrac{\alpha }{4\pi } (\bar{s}P_R b) \left( \bar{\ell } (P_R-P_L)\ell \right) = \sqrt{x_\ell }\, \mathcal {O}_P \ , \end{aligned}$$

(57)

which then in the actual computation gives a contribution to $C_P$. With our procedure, we understand that this contribution does not come from $C_P$ but actually from $\sqrt{x_\ell } C_{RL}^{\mathcal {T}q}$. In other words, by using our definition of operators and of the effective Hamiltonian, we find

(58)

Had we set $p_s=0$ we would have missed the contribution of the dimension-seven operator. We emphasize, once again, that $\mathcal {A}_{RR}^q$ is a non-trivial operator with derivative which cannot be straightforwardly simplified by means of equations of motion.

Finally, after a comparison between ours and the result for $C_P$ presented in Ref. [30] we find⁸

$$\begin{aligned} C_P^\mathrm {Ref.[17]} = \left[ C_P +\frac{\sqrt{x_\ell }}{2 \sin ^2\theta _W} C_{RR}^{\mathcal {T}q}\right] ^\mathrm {(this\,work)}\!\!\!. \end{aligned}$$

(59)

In other words, the Wilson coefficient $C_P$ of Ref. [30] contains the Wilson coefficient of the operator $\mathcal {O}_{RR}^{\mathcal {T}q}$, the matrix element of which is not equal to the matrix element of the operator $\mathcal {O}_P$ but is, instead, suppressed by $m_W$ as we explicitly check in the next section. For that reason the Wilson coefficient of Ref. [30] is not well defined unless the basis of dimension-seven operators is explicitly specified.

7 $B_s \rightarrow \mu ^+ \mu ^-$ and $B \rightarrow K \mu ^+\mu ^-$ in 2HDM

In this section we give the expressions for $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)$ and $\mathcal {B}(B\rightarrow K \mu ^+\mu ^-)$ to which we also include the contributions of the operators given in Eq. (46). Those additional operators were necessary for the appropriate matching procedure between the full and the effective theories. However, since they are suppressed by $1/m_W$ they are expected to be negligible with respect to the dominant operators entering the effective Hamiltonian (21). The purpose of this exercise is to check whether or not the size of the matrix elements of the operators (46) is indeed numerically insignificant for phenomenology.

7.1 $B_s \rightarrow \mu ^+ \mu ^-$

On the basis of Lorentz invariance and invariance of the strong interaction with respect to parity, one can easily verify that $B_s\rightarrow \mu ^+\mu ^-$ is not affected by the operators $\mathcal {O}_{i,j}^{\mathcal {T}q}$ and $\mathcal {O}_{i,j}^{\mathcal {T}\ell }$, with $i,j=L,R$. The expression for the decay rate of this process remains the standard one

$$\begin{aligned} \mathcal {B}(B_s\rightarrow \ell ^+\ell ^-)^\mathrm {th}&= \tau _{B_s}\dfrac{\alpha ^2 G_F^2 m_{B_s}\beta _\ell }{16 \pi ^3} \left| V_{tb}V_{ts}^*\right| ^2 \nonumber \\&\quad \times f_{B_s}^2 m_\ell ^2 \left[ \left| C_{10}-C_{10}^\prime +\dfrac{m_{B_s}^2 (C_P-C_{P}^\prime )}{2 m_\ell (m_b+m_s)} \right| ^2\right. \nonumber \\&\quad \left. + \left| C_S-C_S^\prime \right| ^2\dfrac{m_{B_s}^2(m_{B_s}^2-4 m_\ell ^2)}{4 m_\ell ^2(m_b+m_s)^2}\right] , \end{aligned}$$

(60)

where $\beta _\ell =\sqrt{1-4 m_\ell ^2/m_{B_s}^2}$. To compare Eq. (60) with the available experimental value, one needs to take into account the effects of $B_s- \overline{B}_s$ oscillations which, to a good approximation, amounts to [58]

$$\begin{aligned} \mathcal {B}(B_s\rightarrow \ell ^+\ell ^-)^\mathrm {exp} \approx \dfrac{1}{1-y_s}\mathcal {B}(B_s\rightarrow \ell ^+\ell ^-)^\mathrm {th}, \end{aligned}$$

(61)

where $y_s=\Delta \Gamma _{B_s}/(2 \Gamma _{B_s})=0.061(9)$, experimentally established by the LHCb Collaboration [59]. As we mentioned before, the dimension-seven operators (46) were chosen in such a way that they do not contribute the $B_s\rightarrow \ell ^+\ell ^-$ decay amplitude.

7.2 $B \rightarrow K \mu ^+\mu ^-$

In contrast to $B_s\rightarrow \ell ^+\ell ^-$, the decay $B\rightarrow K \ell ^+\ell ^-$ receives contributions from the operators of the extended basis (46). To write the decay amplitude in a compact form, it is convenient to use the formalism of helicity amplitudes (HAs). In the absence of the (pseudo-)scalar operators, the total amplitude can be schematically written as

$$\begin{aligned} \mathcal {M}=\mathcal {M}_\mu ^L \,\bar{\ell } \gamma ^\mu P_L \ell +\mathcal {M}_{\mu \nu }^L \,\bar{\ell } \sigma ^{\mu \nu } P_L \ell + (L \leftrightarrow R). \end{aligned}$$

(62)

By describing the decay mode as $B\rightarrow K V^*\rightarrow K \ell ^+\ell ^-$, where $V^*$ is a virtual vector boson, one can decompose the total decay amplitude in terms of HAs,

$$\begin{aligned} A_{m}^{L(R)}&= \mathcal {M}_{\mu }^{L(R)}\varepsilon _V^{\mu *}(m), \qquad \text {and} \qquad A_{mn}^{L(R)}\nonumber \\&=\mathcal {M}_{\mu \nu }^{L(R)}\varepsilon _V^{\mu *}(m)\varepsilon _V^{\nu *}(n), \end{aligned}$$

(63)

where $\varepsilon _V^{\mu }(m)$ (with $m,n=0,t,\pm$) are the $V^*$-boson polarization vectors, explicitly defined in Appendix A. We repeat that the above decomposition is valid as long as the scalar and the pseudoscalar operators are not present. To incorporate those contributions unambiguously one can assume the lepton masses to be unequal ($m_{\ell _1}\ne m_{\ell _2}$) and then apply the Ward identities,

$$\begin{aligned} \bar{\ell }_1 \gamma _5 \ell _2 =\dfrac{q^\mu }{m_{\ell _1}+m_{\ell _2}}\bar{\ell }_1 \gamma _\mu \gamma _5 \ell _2,\;\, \bar{\ell }_1 \ell _2 =\dfrac{q^\mu }{m_{\ell _1}-m_{\ell _2}}\bar{\ell }_1 \gamma _\mu \ell _2 , \end{aligned}$$

(64)

to absorb the (pseudo-)scalar terms in the time-like coefficients ${A}_{t}^{L(R)}$. By taking the limit $m_{\ell _1}=m_{\ell _2}$ in the final expression one ends up with the desired HAs and the total decay amplitude, from which is then easy to compute the decay rate [60]. Notice that the contributions from $C_{S,P}^{(\prime )}$ enter the amplitudes $A_S$ and $A_t$ defined as

$$\begin{aligned} A_t&= \lim _{m_{\ell _1}\rightarrow m_{\ell _2}}(A_t^L - A_t^R ), \end{aligned}$$

(65)

$$\begin{aligned} A_S&= \lim _{m_{\ell _1}\rightarrow m_{\ell _2}} \Bigg [ \dfrac{m_{\ell _1}-m_{\ell _2}}{\sqrt{q^2}} (A_t^L + A_t^R ) \Bigg ]. \end{aligned}$$

(66)

More details regarding this point can be found in Ref. [60]. We also need to stress that all the helicity amplitudes are the $q^2$-dependent functions, $A_i\equiv A_i(q^2)$. By applying the method briefly sketched above we obtain

$$\begin{aligned}&\dfrac{\mathrm {d}}{\mathrm {d}q^2}\mathcal {B}(B \rightarrow K \ell ^+ \ell ^-)^\mathrm {th}\nonumber \\&\quad = \dfrac{2(q^2-m_\ell ^2)}{3}[|A_0^L|^2+|A_0^R|^2]\nonumber \\&\qquad +2 m_\ell ^2 \left| A_t\right| ^2+\dfrac{q^2-4m_\ell ^2}{2}|A_S|^2\nonumber \\&\qquad +\dfrac{q^2+2 m_\ell ^2}{3}[|A_{t0}^L-A_{0t}^L|^2+|A_{t0}^R-A_{0t}^R|^2]\nonumber \\&\qquad +4m_\ell ^2 \mathrm {Re}[A_0^{L*} A_0^R]\nonumber \\&\qquad +\dfrac{8 (q^2-4m_\ell ^2)}{3}\left| A_{T5}\right| ^2+\dfrac{4(q^2-4 m_\ell ^2)}{3}\nonumber \\&\qquad \quad \times \mathrm {Re}[A_{T5}^*(A_{t0}^L-A_{0t}^L)-(L\leftrightarrow R)] +4 m_\ell ^2 \nonumber \\&\qquad \quad \times \mathrm {Re}[A_{0t}^{L*}(A_{0t}^R-A_{t0}^R)-A_{t0}^{L*}(A_{0t}^R-A_{t0}^R) ]\nonumber \\&\qquad -2 m_\ell \sqrt{q^2}\,\mathrm {Im} [(A_0^L+A_0^R)^*(A_{t0}^L-A_{0t}^L+(L \leftrightarrow R)) ] , \end{aligned}$$

(67)

where the explicit expressions for the helicity amplitudes are

$$\begin{aligned} A_0^{L(R)}(q^2)&= \mathcal {N}_K\dfrac{\lambda _B^{1/2}}{2\sqrt{q^2}} \left[ f_+(q^2)[(C_9+C_9^\prime )\mp (C_{10}+C_{10}^\prime )] \right. \nonumber \\&\quad +f_T(q^2) \dfrac{2 m_b }{m_B+m_K}(C_7+C_7^\prime )\nonumber \\&\quad \left. -f_T(q^2) \dfrac{q^2}{m_W (m_B+m_K)} [C_{L,L(R)}^{\mathcal {T}q}+C_{R,L(R)}^{\mathcal {T}q}] \right] , \end{aligned}$$

(68)

$$\begin{aligned} A_t (q^2)&= - \mathcal {N}_K f_0(q^2) \dfrac{m_B^2-m_K^2}{\sqrt{q^2}} \nonumber \\&\qquad \times \Bigg [C_{10}+C_{10}^\prime +\dfrac{q^2 \left( C_P+C_P^\prime \right) }{2 m_\ell (m_b-m_s)} \Bigg ], \end{aligned}$$

(69)

$$\begin{aligned} A_S(q^2)&= \mathcal {N}_K f_0(q^2) \dfrac{m_B^2-m_K^2}{m_b-m_s}\left( C_S+C_S^\prime \right) , \end{aligned}$$

(70)

$$\begin{aligned} A_{0t}^{L(R)} (q^2)&= i \mathcal {N}_K\lambda ^{1/2}_B\Bigg [f_T(q^2)\dfrac{C_T}{m_B+m_K}\nonumber \\&\qquad +f_+(q^2)\dfrac{C_{L,L(R)}^{\mathcal {T}\ell }+ C_{R,L(R)}^{\mathcal {T}\ell }}{2 m_W}\Bigg ], \end{aligned}$$

(71)

$$\begin{aligned} A_{t0}^{L(R)} (q^2)&= -i \mathcal {N}_Kf_T(q^2) \dfrac{C_T \lambda _B^{1/2}}{m_B+m_K}, \end{aligned}$$

(72)

$$\begin{aligned} A_{T5} (q^2)&\equiv A^{L(R)}_{+-} = i \mathcal {N}_K f_T(q^2) \dfrac{C_{T5} \lambda _B^{1/2}}{m_B+m_K}, \end{aligned}$$

(73)

where the normalization factor also accounts for the remaining phase space, namely

$$\begin{aligned} |\mathcal {N}_K(q^2)|^2 = \tau _{B_d}\dfrac{\alpha _\mathrm {em}^2 G_F^2 | V_{tb}V_{ts}^*|^2}{512 \pi ^5 m_B^3}\dfrac{\lambda _q^{1/2}}{q^2}\lambda _B^{1/2}. \end{aligned}$$

(74)

For brevity, in the above formulas, we used $\lambda _q=\lambda (\sqrt{q^2},m_\ell ,m_\ell )$ and $\lambda _B=\lambda (m_B,m_K,\sqrt{q^2})$, where $\lambda (a,b,c)\equiv [a^2-(b-c)^2][a^2-(b+c)^2]$. The kinematic conventions and the form factor definitions are collected in Appendix A. In the limit in which the derivative operators vanish we retrieve the usual expression for differential branching fraction [60]. The choice of dimension-seven operators (46) is convenient also because their matrix elements are proportional to the original hadronic matrix elements multiplied by $i q^\mu$. As it can be seen from the above expressions the coefficients $C_{i,j}^{\mathcal {T}\ell }$ and $C_{i,j}^{\mathcal {T}q}$ enters the above formulas with the explicit $1/m_W$-suppression factor. In other words, with the above formulas and by using the Wilson coefficients presented in the previous sections, we see that the derivative operators (46) are indeed irrelevant for phenomenology. Their presence is therefore essential for the unambiguous matching procedure in the computation of Wilson coefficients but they do not alter the phenomenological analysis even at the sub-percent level.

Fig. 6

Results of the scan given in Fig. 1 after imposing the constraints coming from $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp}$ and $\mathcal {B}(B\rightarrow K \mu ^+ \mu ^-)_{\mathrm {high}\,q^2}^\mathrm {exp}$ to $3\sigma$ accuracy. Blue points are allowed by all observables, while gray points are excluded by $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)$, and the red ones are excluded by $\mathcal {B}(B\rightarrow K \mu ^+ \mu ^-)_{\mathrm {high}\,q^2}$

8 Phenomenology and discussion

In this section we use our results for Wilson coefficients and compare the experimental data for the exclusive $b\rightarrow s\ell ^+\ell ^-$ modes with various types of 2HDM. We decided to focus on $\mathcal {B}(B_s\rightarrow \mu ^+ \mu ^-)^\mathrm {exp}=(2.8^{+0.7}_{-0.6})\times 10^{-9}$ [1], and $\mathcal {B}(B\rightarrow K \mu ^+ \mu ^-)_{\mathrm {high}\,q^2}^\mathrm {exp}=(8.5\pm 0.3 \pm 0.4)\times 10^{-8}$ [3], where “high $q^2$” means that the decay rate has been integrated over the interval $q^2\in [15,22]~\mathrm {GeV}^2$. The reason for opting for these decay modes is that the relevant hadronic uncertainties are under good theoretical control. The hadronic quantity entering the $B_s\rightarrow \mu ^+ \mu ^-$ decay amplitude is the decay constant, $f_{B_s}$. It has been abundantly computed by means of numerical simulations of QCD on the lattice (LQCD) and its value is nowadays one of the most accurately computed hadronic quantities as far as $B_{(s)}$-mesons are concerned [29]. The hadronic form factors entering the $B\rightarrow K \mu ^+ \mu ^-$ decay amplitude have been directly computed in LQCD only in the region of large $q^2$’s [62, 63], which explains why we use $\mathcal {B}(B\rightarrow K \mu ^+ \mu ^-)_{\mathrm {high}\,q^2}^\mathrm {exp}$ to do phenomenology. Furthermore, since the bin corresponding to $q^2\in [15,22]~\mathrm {GeV}^2$ is rather wide and away from the very narrow charmonium resonances, the assumption of quark–hadron duality is likely to be valid [61]. By using the recent LQCD results for the form factors provided by the HPQCD [62] and MILC Collaborations [63], the SM results are

$$\begin{aligned}&\mathcal {B}(B\rightarrow K \mu ^+\mu ^-)_{\mathrm {high}\,q^2}\nonumber \\&\quad =\left\{ \biggl . (10.0\pm 0.5)\times 10^{-8}\biggr |_\mathrm {HPQCD} ,\biggl .(10.7\pm 0.5)\times 10^{-8}\biggr |_\mathrm {MILC} \right\} , \end{aligned}$$

(75)

both being about $2 \sigma$ larger than the experimental value measured at LHCb.⁹ Since the current disagreement between theory and experiment needs to be corroborated by more data, we decided to impose all the constraints to $3\sigma$ accuracy. We will then discuss the impact of $\mathcal {B}(B\rightarrow K \mu ^+ \mu ^-)_{\mathrm {high}\,q^2}^\mathrm {exp}$ on 2HDM if the current discrepancy remains, i.e. by requiring the 2HDM to compensate the disagreement between theory (SM) and experiment at the level of $2\sigma$ and more. Notice also that the measured $\mathcal {B}(B_s\rightarrow \mu ^+ \mu ^-)^\mathrm {exp}$ is slightly smaller than predicted, $\mathcal {B}(B_s\rightarrow \mu ^+ \mu ^-)^\mathrm {SM} =(3.65\pm 0.23)\times 10^{-9}$[2].

We now use the results of our scan from Sect. 2.2, and we require the $3\sigma$ agreement between experiment and theory, which means that we add the generic 2HDM Wilson coefficients derived in the previous section to the SM values. The result, in the $(\tan \beta ,m_{H^\pm })$ plane, is shown in Fig. 6 for each type of 2HDM discussed in Sect. 2. We learn that both $\mathcal {B}(B_s\rightarrow \mu ^+ \mu ^-)$ and $\mathcal {B}(B\rightarrow K \mu ^+ \mu ^-)_{\mathrm {high}\,q^2}$ exclude the low $\tan \beta \lesssim 1$ region regardless of the type of 2HDM considered. The limit of exclusion of low $\tan \beta$ coming from $\mathcal {B}(B\rightarrow K \mu ^+ \mu ^-)_{\mathrm {high}\,q^2}$ is slightly larger than the one arising from $\mathcal {B}(B_s\rightarrow \mu ^+ \mu ^-)$. The limit on low $\tan \beta$ obtained in this way for each of our four models is given in Table 2.

Table 2 Allowed values of low $\tan \beta$ (at $99\%$ CL) for the different 2HDMs. See text for details

Fig. 7

Results of the scan in Fig. 1 after imposing the $b\rightarrow s$ constraints to $2\sigma$ accuracy. The hatched area is excluded by $\mathcal {B}(B\rightarrow X_s \gamma )$ at 95% [64]. See Fig. 6 for the color code

Fig. 8

Same as in Fig. 7 but in the $(m_A,m_{H^\pm })$ plane

Besides excluding $\tan \beta \lesssim 1$, it may appear as a surprise that the large $\tan \beta$ are not excluded by these data. The reason for that is the fact that the (pseudo-)scalar Wilson coefficient, with respect to the dominant (axial-)vector one, comes with a term proportional to $(m_{B_s}/m_W)^2$ which suppresses the large $\tan \beta$ values. This feature can be easily verified in the Type II model for which the coefficients $C_{S,P}$, in the large $\tan \beta$ limit, are given in Eq. (43). This is why only a small number of points have been eliminated from our scan of Type II model at large $\tan \beta$ but relatively light $m_{H^\pm }$.

Since the SM value is in slight tension with $\mathcal {B}(B\rightarrow K \mu ^+ \mu ^-)_{\mathrm {high}\,q^2}^\mathrm {exp}$ at the $2.1\sigma$ level, we can now check which of the models discussed in this paper can be made consistent with the experimental data if any disagreement beyond $2\sigma$ between theory (SM) and experiment is to be attributed to 2HDM. It turns out that two such models are Type II and Type Z 2HDM, which we illustrate in Fig. 7. For the other two scenarios (Type I and Type X) the NP contributions are either too small or already in conflict with $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp}$. From Figs. 7 and 8 we see that in order to explain the discrepancy one needs a relatively light charged scalar: (i) $m_{H^\pm } \lesssim 735~\mathrm {GeV}$ and $\tan \beta > 2.3$ in the Type II scenario, and (ii) $m_{H^\pm } \lesssim 380~\mathrm {GeV}$ and $\tan \beta > 3.5$ for the Type Z scenario. Since the masses of the additional scalars are correlated, we see that $m_H$ and $m_A$ become bounded as well, cf. Fig. 8. In the case of Type II and Type Z 2HDM an additional bound on the charged Higgs has been recently derived from the inclusive mode $\mathcal {B}(B\rightarrow X_s \gamma )$. After comparing the experimental spectra with theoretical expressions in which the higher order QCD corrections have been included, the lower bound $m_{H^\pm } >570$ GeV (95% CL) was obtained in Ref. [64] (c.f. also Ref. [65]). This bound is superposed on our results in Figs. 7 and 8, which then also eliminates Type Z 2HDM. Furthermore, we can say that the requirement of agreement between theory and experiment to $2\sigma$, for the quantities discussed in this section, reduces the available space of parameters for Type II 2HDM to $m_{H^\pm } \in (570, 735)~\mathrm {GeV}$, and $\tan \beta \in (16, 35)$, while the available range of values for the mass of the CP-odd Higgs becomes $m_A\in (145, 865)~\mathrm {GeV}$.

Fig. 9

We show the branching fractions of the decay to $\tau$-leptons with respect to their SM predictions, as obtained in the Type II 2HDM, consistent with experimental results for the decays to muons in the final state

In what follows we will assume that the $2\sigma$ disagreement of the measure ${\mathcal {B}}(B\rightarrow K\mu ^+\mu ^-)_{\mathrm {high}\,q^2}^\mathrm {exp}$ with respect to the SM prediction indeed remains as such in the future and discuss the consequences on the decays ${\mathcal {B}}(B_s\rightarrow \tau ^+\tau ^-)$ and ${\mathcal {B}}(B\rightarrow K\tau ^+\tau ^-)_{\mathrm {high}\,q^2}$ if the Type II 2HDM is used to explain the disagreement. From Eq. (60) we can see that

$$\begin{aligned} {{\mathcal {B}}(B_s\rightarrow \tau ^+\tau ^-)\over {\mathcal {B}}(B_s\rightarrow \tau ^+\tau ^-)^{\mathrm {SM}} }= & {} {{\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)\over {\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)^{\mathrm {SM}} }\nonumber \\&- { |C_S^{\tau \tau }|^2\over |C_{10}^\mathrm {SM}|^2 } {m_{B_s}^2\over (m_b+m_s)^2}, \end{aligned}$$

(76)

where the only remaining $m_\ell$ dependence comes from the last numerator in the factor multiplying $|C_S-C_S^\prime |^2$ in Eq. (60). In Fig. 9 we illustrate the validity of the above equality. Notice that a tiny departure from equality comes from the large $\tan \beta$ values which enhance the $C_S$ contribution. In other words, the current experimental result ${\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp}$, which is slightly lower than the one predicted in the SM, is expected to lead to ${\mathcal {B}}(B_s\rightarrow \tau ^+\tau ^-)^\mathrm {exp}$ compatible or slightly lower than predicted in the SM. The cancellation of the lepton mass in ${\mathcal {B}}(B_s\rightarrow \ell ^+\ell ^-)^{\mathrm {2HDM}}$, discussed above, does not occur in ${\mathcal {B}}(B\rightarrow K \ell ^+\ell ^-)^{\mathrm {2HDM}}_{\mathrm {high}-q^2}$. As a result we obtain

$$\begin{aligned} {{\mathcal {B}}(B\rightarrow K \tau ^+\tau ^-)^{\mathrm {Type\,II}}\over {\mathcal {B}}(B\rightarrow K \tau ^+\tau ^-)^{\mathrm {SM}} } \lesssim {{\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)^{\mathrm {Type\,II}}\over {\mathcal {B}}(B\rightarrow K\mu ^+\mu ^-)^{\mathrm {SM}} }\,, \end{aligned}$$

(77)

where we omitted the subscript “high-$q^2$” to avoid too heavy a notation. Illustration is provided in Fig. 9. We can rephrase this observation with an equivalent statement:

$$\begin{aligned} {{\mathcal {B}}(B\rightarrow K \tau ^+\tau ^-)^{\mathrm {Type~II}}\over {\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)^{\mathrm {Type~II}}} < {{\mathcal {B}}(B\rightarrow K \tau ^+\tau ^-)^{\mathrm {SM}}\over {\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)^{\mathrm {SM}}}. \end{aligned}$$

(78)

To be fully explicit, we obtain

$$\begin{aligned} \left. {{\mathcal {B}}(B\rightarrow K \tau ^+\tau ^-)\over {\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)}\right| _{\mathrm {high}-q^2}\!\! \in (1.12,1.14)_\mathrm {SM}, (1.0,1.1)_{\mathrm {Type~II}}.\nonumber \\ \end{aligned}$$

(79)

9 Conclusion

In this paper we computed the leading order Wilson coefficients relevant to the exclusive $b\rightarrow s\ell ^+\ell ^-$ decays in the framework of 2HDM with a softly broken $\mathbb {Z}_2$ symmetry. Most of these Wilson coefficients have been computed previously but in the limit of large $\tan \beta$, which we extend here to a generic setup. We also included $\mathcal {O}(m_b)$ corrections, which were neglected in the previous computations. Regarding the (pseudo-)scalar Wilson coefficients, we elucidated the issue of unambiguous matching of the one-loop amplitudes between the full and effective theories which requires extending the basis of operators in the effective theory by including two types of operators suppressed by $1/m_W$ (altogether, eight new operators). We pointed out that for the appropriate identification of the Z-penguin contribution to the Wilson coefficient $C_P$ it is necessary to keep all external momenta different from zero.

After having computed the full set of Wilson coefficients we were able to make a phenomenological analysis by focusing on ${\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)$ and ${\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)_{\mathrm {high}-q^2}$, the quantities which are measured at LHC and for which the hadronic uncertainties are under good theoretical control (computed in LQCD). After carefully scanning the parameter space of 2HDM with a softly broken $\mathbb {Z}_2$ symmetry, we tested various types of 2HDM against the experimental data for ${\mathcal {B}}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp}$ and ${\mathcal {B}}(B\rightarrow K \mu ^+\mu ^-)_{\mathrm {high}-q^2}^\mathrm {exp}$, and found that to $3\sigma$ the values of low $\tan \beta \lesssim 1$ are excluded for all types of 2HDM considered here.

If, instead, we require the $2\sigma$ agreement with experiment, then only Type II and Type Z models provide a viable description of the data. After combining ours with the bound on the charged Higgs deduced from the inclusive $b\rightarrow s\gamma$ decay, we find that the Type Z model can be discarded and

$$\begin{aligned}&\mathrm {Type~II}\ : \quad m_{H^\pm } \in (570,735)~\mathrm{GeV}, \qquad m_A\in (145,865)~\mathrm{GeV}, \nonumber \\&\qquad \tan \beta \in (16,35). \end{aligned}$$

(80)

We also discussed the repercussions of the current results on the decays ${\mathcal {B}}(B_s\rightarrow \tau ^+\tau ^-)$ and ${\mathcal {B}}(B\rightarrow K \tau ^+\tau ^-)_{\mathrm {high}-q^2}$.

Appendices

A Conventions and kinematics

1.1 Angular conventions

We adopt the same angular conventions for $B(p)\rightarrow K (k) \ell ^+(p_+)\ell ^-(p_-)$ as those used in Refs. [60, 66]. In the B-meson rest frame, the leptonic and hadronic four-momenta are chosen as $q^\mu =p_+^\mu +p_-^\mu =(q_0,0,0,q_z)$ and $k^\mu =(k_0,0,0,-q_z)$, where

$$\begin{aligned} q_0&= \dfrac{m_B^2+q^2-m_K^2}{2 m_B},\qquad k_0=\dfrac{m_B^2+m_{K}^2-q^2}{2m_B},\quad \mathrm {and}\\ q_z&=\dfrac{\lambda ^{1/2}(m_B,m_K,\sqrt{q^2})}{2 m_B}. \end{aligned}$$

In the dilepton rest frame the components of the leptonic four-vectors are given by

$$\begin{aligned} p_-^\mu&= (E_\ell ,|p_\ell |\sin \theta _\ell ,0,|p_\ell |\cos \theta _\ell ),\\ p_+^\mu&= (E_\ell ,-|p_\ell |\sin \theta _\ell ,0,-|p_\ell |\cos \theta _\ell ), \end{aligned}$$

where $E_\ell =\sqrt{q^2}/2$, and $\theta _\ell$ is the angle between $\ell ^-$ (in the dilepton rest frame) and the line of flight of the two leptons (in the B-meson rest frame). The momentum $p_\ell$ is given by

$$\begin{aligned} |p_\ell |=\dfrac{\lambda ^{1/2}(\sqrt{q^2},m_\ell ,m_\ell )}{2 m_B}. \end{aligned}$$

(81)

1.2 Polarization vectors

In the B-meson rest frame we take the polarization vectors of the virtual vector boson $V^*$ to be

$$\begin{aligned} \varepsilon ^\mu _V(\pm )&= \dfrac{1}{\sqrt{2}}(0,\pm 1,i,0),\quad \varepsilon ^\mu _V(0) = \dfrac{1}{\sqrt{q^2}}(q_z,0,0,q_0),\nonumber \\ \varepsilon ^\mu _V(t)&= \dfrac{1}{\sqrt{q^2}}(q_0,0,0,q_z). \end{aligned}$$

(82)

These vectors are orthonormal and satisfy the completeness relation,

$$\begin{aligned} \sum _{n,n^\prime } \varepsilon _V^{\mu *}(n)\varepsilon _V^\nu (n^\prime )g_{nn^\prime }=g^{\mu \nu }, \end{aligned}$$

(83)

where $n,n^\prime \in \lbrace t,0,\pm \rbrace$, and $g_{n n^\prime }=\mathrm {diag}(1,-1,-1,-1)$.

1.3 Hadronic matrix elements

For completeness we also give the definitions of the decay constant ($f_{B_s}$) and of the form factors [$f_{+,0,T}(q^2)$], quantities which parametrize the hadronic matrix elements relevant to the processes discussed in this paper:

$$\begin{aligned} \langle 0 | \bar{b} \gamma _\mu \gamma _5 s | B_s(p) \rangle&= i p_\mu f_{B_s},\nonumber \\ \langle \bar{K}(k)|\bar{s}\gamma _\mu b|\bar{B}(p)\rangle&= \Big [(p+k)_\mu - \frac{m_B^2-m_K^2}{q^2}q_\mu \Big ]f_+(q^2)\nonumber \\&\quad +\frac{m_B^2-m_K^2}{q^2} q_\mu f_0(q^2),\nonumber \\ \langle \bar{K}(k)|\bar{s} b|\bar{B}(p)\rangle&= {1 \over m_b-m_s}q^\mu \langle \bar{K}(k)|\bar{s} \gamma _\mu b|\bar{B}(p)\rangle \nonumber \\&= {m_B^2-m_K^2 \over m_b-m_s}f_0(q^2) ,\nonumber \\ \langle \bar{K}(k)|\bar{s}\sigma _{\mu \nu } b|\bar{B}(p)\rangle&= -i (p_\mu k_\nu -p_\nu k_\mu )\frac{2 f_T(q^2,\mu )}{m_B+m_K}, \end{aligned}$$

(84)

where for $B\rightarrow K\ell ^+\ell ^-$ the kinematically accessible $q^2$ values lie in the interval $4m_\ell ^2\le q^2\le (m_B-m_K)^2$. Notice that we do not write explicitly the scale dependence of the quark masses, nor of the scalar and tensor densities and of the form factor $f_T(q^2)$. In the actual computations the ${\overline{\mathrm{MS}}}$ values of these quantities are taken at $\mu =m_b$.

B Feynman rules

In this appendix we collect the Feynman rules used in our computation. For the couplings between the gauge bosons and the scalars we have

(85)

where $\lambda _{W^+W^-}^h=\sin (\beta -\alpha )$, $\lambda _{W^+W^-}^H=\cos (\beta -\alpha )$ and $\lambda _{W^+W^-}^A=0$. Similarly, we also have

(86)

(87)

where $\lambda _{H^\pm W^\mp }^h=\cos (\beta -\alpha )$, $\lambda _{H^\pm W^\mp }^H=-\sin (\beta -\alpha )$, and $\lambda _{H^\pm W^\mp }^A=\mp i$, depending on the charges of the initial particles. For the trilinear scalar interactions, we have

(88)

where the trilinear couplings read

$$\begin{aligned} \lambda _{H^+ H^-}^h&=- \dfrac{m_h^2[3 \cos (\alpha +\beta )+\cos (\alpha -3\beta )]+4\sin (2\beta ) \sin (\beta -\alpha )m_{H^\pm }^2-4M^2\cos (\alpha +\beta )}{2 v^2 \sin (2\beta )},\nonumber \\ \lambda _{H^+ H^-}^H&= - \dfrac{m_H^2[3 \sin (\alpha +\beta )+\sin (\alpha -3\beta )]+4\sin (2\beta ) \cos (\beta -\alpha )m_{H^\pm }^2-4M^2\sin (\alpha +\beta )}{2 v^2 \sin (2\beta )},\nonumber \\ \lambda _{H^+ H^-}^A&= 0. \end{aligned}$$

(89)

These results agree with the ones given in Refs. [30, 67] after the appropriate change of basis and/or conventions.¹⁰

C Scalar penguins and auxiliary functions

In this appendix we give the expressions for the Wilson coefficients generated by each diagram shown in Fig. 5. We also give the expressions for the auxiliary functions ($f_i$ and $g_i$) used in this paper.

The penguins arising from coupling to $\varphi _i^0\in \lbrace h,H,A\rbrace$ contribute to the effective coefficient $C_{S,P}$ and can be generically written as

$$\begin{aligned} C_{S}^{\mathrm {NP},\varphi _i^0}&= \dfrac{\sqrt{x_b x_\ell }}{\sin ^2\theta _W} \sum _{k=1}^{18} \dfrac{m_t^2}{m_{\varphi _i^0}} \mathrm {Re}\left( \xi _\ell ^{\varphi _i^0}\right) \widehat{C}^{k,\varphi _i^0}, \end{aligned}$$

(90)

$$\begin{aligned} C_{P}^{\mathrm {NP},\varphi _i^0}&= \dfrac{\sqrt{x_b x_\ell }}{\sin ^2\theta _W}\sum _{k=1}^{18} \dfrac{m_t^2}{m_{\varphi _i^0}} i \,\mathrm {Im}\left( \xi _\ell ^{\varphi _i^0}\right) \widehat{C}^{k,\varphi _i^0}, \end{aligned}$$

(91)

where $\widehat{C}^{k,\varphi _i^0}$ is a common coefficient generated by the diagram k, with $k=1,\dots ,18$. Since, in our framework, $\zeta _{\ell }^h,\zeta _{\ell }^H\in \mathbb {R}$ and $\zeta _{\ell }^A\in i\times \mathbb {R}$, it is clear that the CP-even scalars h and H contribute only to $C_S$, whereas the CP-odd Higgs A contributes only to $C_P$, as expected from the assumption of CP conservation. We obtain in the unitary gauge

$$\begin{aligned} \widehat{C}^{1,\varphi _i^0}&= \dfrac{\xi _u^{\varphi _i^0}}{4}\Bigg \lbrace \zeta _d \zeta _u^*\dfrac{x_t}{x_{H^\pm }-x_t}\Bigg [1- \dfrac{x_{H^\pm }}{x_{H^\pm }-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ]\nonumber \\&\quad +|\zeta _u|^2 \dfrac{x_t}{2(x_{H^\pm }-x_t)^2}\Bigg [\dfrac{3 x_t-x_{H^\pm }}{2}+\dfrac{x_{H^\pm }(x_{H^\pm }-2x_t)}{x_{H^\pm }-x_t}\nonumber \\&\qquad \quad \times \log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ]\Bigg \rbrace \nonumber \\&\quad + \dfrac{\xi _u^{\varphi _i^0*}}{4}\Bigg \lbrace \zeta _d\zeta _u^*\Bigg [\Lambda -\dfrac{x_t}{x_{H^\pm }-x_t}-\dfrac{x_{H^\pm }^2}{(x_{H^\pm }-x_t)^2}\log x_{H^\pm }\nonumber \\&\quad +\dfrac{x_t(2 x_{H^\pm }-x_t)}{(x_{H^\pm }-x_t)^2}\log x_t\Bigg ]\nonumber \\&\quad +|\zeta _u|^2\dfrac{x_t}{2(x_{H^\pm }-x_t)^2}\Bigg [\dfrac{3 x_{H^\pm }-x_t}{2}-\dfrac{x_{H^\pm }^2}{x_{H^\pm }-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ]\Bigg \rbrace , \end{aligned}$$

(92)

$$\begin{aligned} \widehat{C}^{2,\varphi _i^0}&= -\varepsilon _{\varphi _i^0} \dfrac{\sin ^2\theta _W\lambda _{H^+H^-}^{\varphi _i^0}}{4\pi \alpha (x_{H^\pm }-x_t)}\Bigg \lbrace \zeta _d\zeta _u^*\Bigg [\dfrac{x_t}{x_{H^\pm }-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) -1\Bigg ] \nonumber \\&\quad +|\zeta _u|^2\Bigg [\dfrac{x_t^2}{2(x_{H^\pm }-x_t)^2}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) +\dfrac{x_{H^\pm }-3 x_t}{4(x_{H^\pm }-x_t)}\Bigg ]\Bigg \rbrace , \end{aligned}$$

(93)

$$\begin{aligned} \widehat{C}^{3,\varphi _i^0}&= \dfrac{\xi _d^{\varphi _i^0}}{4}\zeta _d \zeta _u^*\Bigg [-\Lambda +\dfrac{x_{H^\pm }}{x_{H^\pm }-x_t}\log x_{H^\pm }-\dfrac{x_t}{x_{H^\pm }-x_t}\log x_t \Bigg ], \end{aligned}$$

(94)

$$\begin{aligned} \widehat{C}^{4,\varphi _i^0}&= 0, \end{aligned}$$

(95)

$$\begin{aligned} \widehat{C}^{5,\varphi _i^0}&= \dfrac{1}{4}\Bigg \lbrace \xi _u^{\varphi ^0_{i}*}\Bigg [\Lambda -\dfrac{5x_t^2-13 x_t-2}{4(x_t-1)^2}-\dfrac{2x_t^3-6x_t^2+9x_t-2}{2(x_t-1)^3}\log x_t\Bigg ]\nonumber \\&\quad +\xi _u^{\varphi ^0_{i}}\Bigg [\dfrac{\Lambda }{2}-\dfrac{2x_t^2-x_t-7}{4(x_t-1)^2}-\dfrac{x_t^3-3x_t^2+3x_t+2}{2(x_t-1)^2}\log x_t\Bigg ]\Bigg \rbrace , \end{aligned}$$

(96)

$$\begin{aligned} \widehat{C}^{6,\varphi _i^0}&= \varepsilon _{\varphi _i^0}\dfrac{\lambda _{W^+W^-}^{\varphi _i^0}}{8}\Bigg [-3\Lambda +\dfrac{x_t^2-2x_t-11}{2(x_t-1)^2}\nonumber \\&\quad +\dfrac{3 x_t(x_t^2-3x_t+4)}{(x_t-1)^3}\log x_t\Bigg ], \end{aligned}$$

(97)

$$\begin{aligned} \widehat{C}^{7,\varphi _i^0}&= \widehat{C}^{8,\varphi _i^0}=0, \end{aligned}$$

(98)

$$\begin{aligned} \widehat{C}^{9,\varphi _i^0}&= \dfrac{\lambda _{H^+W^-}^{\varphi _i^0}}{8}\zeta _u^*\Bigg [\dfrac{1}{2}-\Lambda + \dfrac{x_{H^\pm }(x_{H^\pm }+2)\log x_{H^\pm }}{(x_{H^\pm }-1)(x_{H^\pm }-x_t)}\nonumber \\&\quad -\dfrac{x_t(x_t+2)\log x_t}{(x_t-1)(x_{H^\pm }-x_t)}\Bigg ], \end{aligned}$$

(99)

$$\begin{aligned} \widehat{C}^{10,\varphi _i^0}&= \dfrac{\lambda _{H^+W^-}^{\varphi _i^0\,*}}{4}\Bigg \lbrace -\dfrac{\zeta _u}{2}\Bigg [\dfrac{x_t(x_{H^\pm }x_t-4x_{H^\pm }+3x_t)}{(x_t-1)(x_{H^\pm }-x_t)^2}\log x_t\nonumber \\&\quad -\dfrac{x_{H^\pm }(x_{H^\pm }x_t-3 x_{H^\pm }+2 x_t)}{(x_t-1)(x_{H^\pm }-x_t)^2}\log x_{H^\pm }\nonumber \\&\quad +\dfrac{x_{H^\pm }}{x_{H^\pm }-x_t}\Bigg ]+\zeta _d\Bigg [-\Lambda +\dfrac{x_{H^\pm }\log x_{H^\pm }}{x_{H^\pm }-x_t}-\dfrac{x_t\log x_t}{x_{H^\pm }-x_t}\Bigg ]\Bigg \rbrace , \end{aligned}$$

(100)

where the couplings $\lambda _{W^+W^-}^{\varphi _i^0}$ and $\lambda _{H^\pm W^\mp }^{\varphi _i^0}$ are defined below Eqs. (85) and (87), respectively. The coefficient $\varepsilon _{\varphi _i^0}=-1$ for $\varphi _i^0 =A$, and $+1$ otherwise. Moreover, $\Lambda =- \frac{2 \mu ^{D-4}}{D-4}-\gamma _E+\log 4\pi -\log \left( \frac{m_W^2}{\mu ^2}\right) +1$ contains an ultraviolet divergence which cancels out after summing up all the diagrams. The diagrams (9.11)–(9.18) do not contribute in our computation, owing to the fact that we work in the unitary gauge. To make sure that our resulting (total) expressions are gauge independent we performed the computation in the Feynman gauge too. In comparison with Ref. [30], we only disagree with one of the signs in the expression for $\widehat{C}^{5,\varphi _i^0}$, which we believe is a typo.

The auxiliary functions $g_{0,1,2}$ used in Eq. (42) are defined by

$$\begin{aligned} g_0&= \dfrac{1}{4(x_{H^\pm }-x_t)}\Bigg \lbrace \zeta _d \zeta _u^*\Bigg [ \dfrac{x_t}{x_{H^\pm }-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) -1\Bigg ] \nonumber \\&\quad + |\zeta _u|^2 \Bigg [ \dfrac{x_t^2}{2(x_{H^\pm }-x_t)^2}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) +\dfrac{x_{H^\pm }-3 x_t}{4(x_{H^\pm }-x_t)} \Bigg ] \Bigg \rbrace \end{aligned}$$

(101)

$$\begin{aligned} g_1&= -\dfrac{3}{4}+\zeta _d \zeta _u^*\dfrac{x_t}{x_{H^\pm }-x_t}\Bigg [ 1-\dfrac{x_{H^\pm }}{x_{H^\pm }-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ]\nonumber \\&\quad +|\zeta _u|^2 \dfrac{x_t}{2(x_{H^\pm }-x_t)^2}\Bigg [ \dfrac{x_{H^\pm }+x_t}{2}-\dfrac{x_{H^\pm }x_t}{x_{H^\pm }-x_t} \log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ], \end{aligned}$$

(102)

$$\begin{aligned} g_2&= \zeta _d(\zeta _d\zeta _u^*+1) f_1 (x_t,x_{H^\pm })+\zeta _d\left( \zeta _u^*\right) ^2 f_2 (x_t,x_{H^\pm })\nonumber \\&\quad +\zeta _d \left| \zeta _u\right| ^2 f_3 (x_t,x_{H^\pm })\nonumber \\&\quad +\zeta _u \left| \zeta _u\right| ^2 f_4 (x_t,x_{H^\pm })- \zeta _u^*\left| \zeta _u\right| ^2 f_5 (x_t,x_{H^\pm })\nonumber \\&\quad +\zeta _u f_6 (x_t,x_{H^\pm })-\zeta _u^*f_7 (x_t,x_{H^\pm }), \end{aligned}$$

(103)

with

$$\begin{aligned} f_1(x_t,x_{H^\pm })&=\frac{1}{2(x_{H^\pm }-x_t)}[-x_{H^\pm }+x_t\nonumber \\&\quad +x_{H^\pm }\log x_{H^\pm }-x_t \log x_t], \end{aligned}$$

(104)

$$\begin{aligned} f_2(x_t,x_{H^\pm })&=\frac{1}{2(x_{H^\pm }-x_t)}\Bigg [x_t-\frac{x_{H^\pm }x_t}{x_{H^\pm }-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ], \end{aligned}$$

(105)

$$\begin{aligned} f_3(x_t,x_{H^\pm })&=\frac{1}{2(x_{H^\pm }-x_t)}\Bigg [x_{H^\pm }-\frac{x_{H^\pm }^2 \log x_{H^\pm }}{x_{H^\pm }-x_t}\nonumber \\&\quad + \frac{x_t(2 x_{H^\pm }-x_t)\log x_t}{x_{H^\pm }-x_t}\Bigg ], \end{aligned}$$

(106)

$$\begin{aligned} f_4(x_t,x_{H^\pm })&= \frac{1}{4(x_{H^\pm }-x_t)^2} \Bigg [ \frac{x_t (3 x_{H^\pm }-x_t)}{2}\nonumber \\&\quad -\frac{x_{H^\pm }^2 x_t}{x_{H^\pm }-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ] , \end{aligned}$$

(107)

$$\begin{aligned} f_5(x_t,x_{H^\pm })&= \frac{1}{4(x_{H^\pm }-x_t)^2} \Bigg [ \frac{x_t (x_{H^\pm }-3 x_t)}{2}\nonumber \\&\quad - \frac{x_{H^\pm } x_t(x_{H^\pm }-2x_t)}{x_{H+}-x_t}\log \left( \dfrac{x_{H^\pm }}{x_t}\right) \Bigg ], \end{aligned}$$

(108)

$$\begin{aligned} f_6(x_t,x_{H^\pm })&= \frac{1}{2(x_{H^\pm }-x_t)}\Bigg [\frac{x_t(x_t^2-3 x_{H^\pm } x_t+9 x_{H^\pm }-5 x_t-2)}{4(x_t-1)^2} +\frac{x_{H^\pm }(x_{H^\pm }x_t-3 x_{H^\pm }+2x_t)\log x_{H^\pm }}{2 (x_{H^\pm }-1)(x_{H^\pm }-x_t)}\nonumber \\&+\frac{x_{H^\pm }^2(-2 x_t^3+6 x_t^2-9 x_t+2)+3 x_{H^\pm }x_t^2(x_t^2-2 x_t+3)-x_t^2(2 x_t^3-3 x_t^2+3x_t+1)}{2 (x_t-1)^3(x_{H^\pm }-x_t)}\log x_t\Bigg ], \end{aligned}$$

(109)

$$\begin{aligned} f_7(x_t,x_{H^\pm })&= \frac{1}{2(x_{H^\pm }-x_t)}\Bigg [\frac{(x_t^2+x_t-8)(x_{H^\pm }-x_t)}{4(x_t-1)^2}-\frac{x_{H^\pm }(x_{H^\pm }+2)}{2(x_{H^\pm }-1)}\log x_{H^\pm } \nonumber \\&+\frac{x_{H^\pm }(x_t^3-3x_t^2+3x_t+2)+3(x_t-2)x_t^2}{2(x_t-1)^3}\log x_t\Bigg ]. \end{aligned}$$

(110)

Notice that in the above expressions we assumed the couplings $\zeta _f \in \mathbb {C}$ in order to keep our formulas as general as possible, although in the body of the paper we consistently used $\zeta _f\in \mathbb {R}$.

1.1 C.1 Wilson coefficients for the derivative operators

In this subsection we present the explicit expressions for the Wilson coefficients relevant to the derivative operators given in Eq. (46). From the Z-penguins we obtain

$$\begin{aligned} C_{RR}^{\mathcal {T}q}= & {} |\zeta _u|^2\dfrac{\sqrt{x_b}x_t}{72}\Bigg \lbrace \frac{3(x_{H^\pm }^2-5 x_{H^\pm }x_t-2x_t^2)}{(x_{H^\pm }-x_t)^3}\nonumber \\&+\frac{18 x_{H^\pm }x_t^2}{(x_{H^\pm }-x_t)^4}\log \left( \frac{x_{H^\pm }}{x_t}\right) \nonumber \\&-2 \sin ^2\theta _W \Bigg [\dfrac{7 x_{H^\pm }^2-5 x_{H^\pm }x_t-8 x_t^2}{(x_{H^\pm }-x_t)^3}\nonumber \\&-\dfrac{6 x_{H^\pm } x_t(2 x_{H^\pm }-3x_t)}{(x_{H^\pm }-x_t)^4}\log \left( \frac{x_{H^\pm }}{x_t}\right) \Bigg ] \Bigg \rbrace \nonumber \\&+\zeta _u^*\zeta _d\dfrac{\sqrt{x_b}x_t}{24}\Bigg \lbrace \dfrac{3(x_{H^\pm }-3 x_t)}{(x_{H^\pm }-x_t)^2}-\dfrac{6 x_{H^\pm }(x_{H^\pm }- 2 x_t)}{ x_{H^\pm }-x_t}\nonumber \\&\times \log \left( \frac{x_{H^\pm }}{x_t}\right) \nonumber \\&+4\sin ^2\theta _W\left[ \dfrac{5 x_t-3 x_{H^\pm }}{(x_{H^\pm }-x_t)^2}+\dfrac{2 x_{H^\pm }(2x_{H^\pm }-3x_t)}{(x_{H^\pm }-x_t)^3}\right. \nonumber \\&\left. \times \log \left( \frac{x_{H^\pm }}{x_t}\right) \right] \Bigg \rbrace , \end{aligned}$$

(111)

and $C_{RL}^{\mathcal {T}q} = C_{RR}^{\mathcal {T}q} \left( 1-\dfrac{1}{2 \sin ^2\theta _W}\right)$. Similarly, from the box diagrams we get

$$\begin{aligned} C^{\mathcal {T}\ell }_{LL}= & {} -\zeta _u\zeta _\ell ^*\frac{\sqrt{x_\ell }x_t}{4(x_{H^\pm }-x_t)\sin ^2\theta _W}\nonumber \\&\Bigg [-\frac{1}{(x_{H^\pm }-1)}+\dfrac{x_{H^\pm }(1-x_{H^\pm })\log x_t}{(x_{H^\pm }-x_t)(x_t-1)(x_{H^\pm }-1)}\nonumber \\&-\dfrac{x_{H^\pm }(x_t+1-2 x_{H^\pm })\log x_{H^\pm }}{(x_{H^\pm }-x_t)(x_{H^\pm }-1)^2}\Bigg ], \end{aligned}$$

(112)

and $C^{\mathcal {T}\ell }_{LL}=(C^{\mathcal {T}\ell }_{LR})^*$.

1.2 C.2 Wilson coefficients suppressed by $m_\ell$

In addition to the Wilson coefficients given in Sect. 2, in the computation of the box diagrams one gets contributions suppressed by the lepton mass. For completeness, we give these contributions here. We obtain

$$\begin{aligned} C_{T(5)}^\mathrm {NP,\,box}&=\zeta _u^*\zeta _\ell \; \dfrac{\sqrt{x_b x_\ell }x_t}{32(x_{H^\pm }-x_t)\sin ^2\theta _W} \nonumber \\&\times \Bigg [ \frac{x_t \log x_t}{(x_t-1)(x_{H^\pm }-x_t)}\nonumber \\&\quad -\frac{x_{H^\pm } \log x_{H^\pm }}{(x_{H^\pm }-1)(x_{H^\pm }-x_t)} +\frac{x_t-\log x_t-1}{(x_t-1)^2}\Bigg ] \end{aligned}$$

(113)

and

$$\begin{aligned} C_9^\mathrm {NP,\,box}&= \frac{x_\ell x_t}{16\sin ^2\theta _W}\Bigg \lbrace |\zeta _u|^2 |\zeta _\ell |^2\Bigg [-\frac{1}{x_{H^\pm }-x_t}\nonumber \\&\quad +\frac{x_t}{(x_{H^\pm }-x_t)^2}\log \left( \frac{x_{H^\pm }}{x_t}\right) \Bigg ]\nonumber \\&\quad +2 \mathrm {Re}[\zeta _u\zeta _\ell ^*] \Bigg [\dfrac{(x_{t}+2)\log x_t}{(x_{H^\pm }-x_t)(x_t-1)} \nonumber \\&\quad -\dfrac{(x_{H^\pm }+2)\log x_{H^\pm }}{(x_{H^\pm }-x_t)(x_{H^\pm }-1)} \Bigg ]\Bigg \rbrace +2\sqrt{x_\ell } \,\mathrm {Re}\left( C_{LL}^{\mathcal {T}\ell }\right) , \end{aligned}$$

(114)

$$\begin{aligned} C_{10}^\mathrm {NP,\,box}&= \frac{x_\ell x_t}{16\sin ^2\theta _W}\Bigg \lbrace |\zeta _u|^2 |\zeta _\ell |^2\Bigg [-\frac{1}{x_{H^\pm }-x_t}\nonumber \\&\quad +\frac{x_t}{(x_{H^\pm }-x_t)^2}\log \left( \frac{x_{H^\pm }}{x_t}\right) \Bigg ]\nonumber \\&+2 \mathrm {Re}[\zeta _u\zeta _\ell ^*] \Bigg [\dfrac{(x_{t}-2)\log x_t}{(x_{H^\pm }-x_t)(x_t-1)} \nonumber \\&\quad -\dfrac{(x_{H^\pm }-2)\log x_{H^\pm }}{(x_{H^\pm }-x_t)(x_{H^\pm }-1)} \Bigg ]\Bigg \rbrace . \end{aligned}$$

(115)