Revisiting $B_s \rightarrow \mu ^+\mu ^-$ in the two-Higgs doublet models with $Z_2$ symmetry

Abstract

We revisit the rare leptonic decay $B_s \rightarrow \mu ^+ \mu ^-$ in the two-Higgs doublet models (2HDMs) with a softly broken $Z_2$ symmetry, namely type-I, type-II, type-X and type-Y 2HDMs. We have derived the relevant full one-loop Wilson coefficients of the four 2HDMs from the recent calculation in the aligned 2HDM by Li, Lu and Pich, which could be mapped to all the four 2HDMs for both large and small $\tan \beta$. It is found that a new term associated with the soft $Z_2$ symmetry breaking parameter M can be enhanced by $\tan ^2\beta$ in the type-II 2HDM, which has not been considered in the literature. Imposing both theoretical and experimental constraints, we have renewed the bounds on the parameter spaces of the four 2HDMs. Different from our previous paper, however, we find that all the four 2HDMs give sizable and similar contributions to $\overline{\mathcal B}(B_s \rightarrow \mu ^+ \mu ^-)$ within the stringently restricted parameter spaces, but very tiny as regards the mass-eigenstate rate asymmetry $\mathcal A_{\Delta \Gamma }$; this makes it unfeasible to discriminate the four types of 2HDM with the correlations between the observables in $B_s \rightarrow \mu ^+ \mu ^-$ decay.

1 Introduction

The discovery [1, 2] of a new boson with a mass close to 125 GeV has been well anticipated as the standard model Higgs boson [3–6] and provided the first experimental evidence of the Higgs mechanism [7–9]. It is a great triumph, but not an end, of the giant campaign for Higgs hunting in the development of particle physics. Although the subsequent more precise measurements [10–14] at the LHC have shown the properties of the Higgs boson are well consistent with the predictions of the standard model (SM), the precision of the current experimental data still leave open the possibility of an extended Higgs sector [15, 16]. Among many new physics scenarios beyond the SM, the two Higgs doublet models (2HDM) [17–19] are the simplest extensions of the SM.

In the 2HDMs, an additional Higgs doublet is introduced to the SM Higgs sector, which could result in rich phenomena, in collider physics [20–29], flavor physics [30–37], neutrino physics [38], dark matter [39–41], and cosmology [42, 43]. However, unlike the SM, unwanted tree-level flavor-changing neutral current (FCNC) interactions in the 2HDM are not forbidden by the Glashow–Illiopoulos–Maiani (GIM) mechanism. Besides some other solutions [44–48], this issue is usually addressed by the natural flavor conservation (NFC) hypothesis through imposing a discrete $Z_2$ symmetry [49]. According to different $Z_2$ charge assignments, there are four types of the NFC 2HDM, referred to as the type-I, type-II, type-X and type-Y 2HDM, respectively. Of course, there are new parameters in the 2HDMs to be determined or excluded by the measurements of electroweak processes. To this end, B-meson decays are usually employed to constrain their parameter spaces.

Among the rare B-meson decays, the leptonic processes $B_q \rightarrow \mu ^+ \mu ^-$ ($q=d$, or s) are of special interest [50, 51]. They suffer from very few hadronic uncertainties and are induced by FCNC transitions, which make them sensitive probes to the effects of physics beyond the SM, especially models with a non-standard Higgs sector [52–56]. Recently, the next-to-leading order (NLO) electroweak corrections and the next-to-next-to-leading order (NNLO) QCD corrections [57–59] in the SM have been calculated. On the BSM side, a full one-loop calculation in the aligned 2HDM (A2HDM) has been performed in Ref. [60].

Motivated by this progress, in this paper we perform a detailed study of the $B_s \rightarrow \mu ^+ \mu ^-$ decay within the 2HDMs with $Z_2$ symmetry. At present, this process is calculated in the type-II 2HDM in large $\tan \beta$ limit only [61–63]. Using the Higgs base correspondence between the A2HDM and the 2HDMs, we will derive the relevant full one-loop Wilson coefficients of the four variant 2HDMs contributing to the $B_s \rightarrow \mu ^+ \mu ^-$ decay from the recent A2HDM results [60] without the large $\tan \beta$ approximation. We also investigate the possibility to discriminate the four different types of 2HDM in the light of the recent collider and flavor physics data, as an update of our previous work [64]. We combine the constraints from $B_{s,d}\rightarrow \mu ^+ \mu ^-$, $B_{s,d}$–$\bar{B}_{s,d}$ mixing, $B\rightarrow \tau \nu$ and $\bar{B} \rightarrow X_s \gamma$ [65, 66], with the experimental data from the direct search for Higgs bosons at LEP [67], Tevatron [68, 69] and LHC [70, 71], and the constraints from perturbativity, tree-level vacuum stability and perturbative unitary. For the $B_s \rightarrow \mu ^+ \mu ^-$ decay, the correlations between its branching ratio and the mass-eigenstate rate asymmetry $\mathcal A_{\Delta \Gamma }$ are also reevaluated with the constrained parameter space of the 2HDMs obtained in this paper. We have found that $\mathcal A_{\Delta \Gamma }$ can slightly deviate from the SM prediction in the type-II 2HDM only, and that the ratio of time-integrated $\overline{\mathcal B}(B_s \rightarrow \mu ^+ \mu ^-)$ gets similar contributions from the four 2HDMs; this makes it very hard to discriminate the four types of 2HDMs with the correlation between $\mathcal A_{\Delta \Gamma }$ and $\overline{\mathcal B}(B_s \rightarrow \mu ^+ \mu ^-)$ as suggested in our previous work [64].

The paper is organized as follows. In Sect. 2, we give a brief overview of the $B_s \rightarrow \mu ^+ \mu ^-$ decay. In Sect. 3, full one-loop contributions from the 2HDMs with $Z_2$ symmetry are derived explicitly. In Sect. 4, we give our detailed numerical results and discussions. We conclude in Sect. 5. The relevant theoretical formulas are recapitulated in the Appendix.

2 $B_s \rightarrow \mu ^+ \mu ^-$ in the SM

In the SM, the leptonic decays $B_q \rightarrow \mu ^+ \mu ^-$ ($q=d$ or s) arise from the W box and Z penguin diagrams. Generally, these decays can be described by the low-energy effective Hamiltonian

$$\begin{aligned} \mathcal H_\mathrm{eff}=-\frac{G_F}{\sqrt{2}}\frac{\alpha _e}{\pi s_W^2}V_{tb} V_{tq}^* ( C_{10} \mathcal O_{10} + C_S \mathcal O_S + C_P \mathcal O_P ), \end{aligned}$$

(2.1)

where $\alpha _e$ denotes the QED fine-structure constant and $V_{ij}$ the CKM matrix elements. The semi-leptonic operators are defined as

$$\begin{aligned}&\mathcal O_{10}= ( \bar{q} \gamma _\mu P_L b ) ( \bar{\mu }\gamma ^\mu \gamma _5 \mu ),\nonumber \\&\mathcal O_S=\frac{m_\mu m_b}{m_W^2} ( \bar{q} P_R b ) ( \bar{\mu }\mu ) ,\nonumber \\&\mathcal O_P=\frac{m_\mu m_b}{m_W^2} (\bar{q} P_R b ) ( \bar{\mu }\gamma _5 \mu ). \end{aligned}$$

(2.2)

In the SM, the contributions from the scalar operators $\mathcal O_S$ and $\mathcal O_P$ are highly suppressed (the corresponding Wilson coefficients are given in Eq. (A.1)), but $C_{10}$ will play the dominant role. Its explicit expressions up to the NLO QCD corrections can be found in Refs. [72–74]. Recently, calculations of the NLO EW [58] and NNLO QCD [59] corrections have also been completed [57]. This progress will be incorporated into our calculations.

With the effective Hamiltonian Eq. (2.1), the branching ratio of $B_q \rightarrow \mu ^+ \mu ^-$ reads

$$\begin{aligned} \mathcal B(B_q\rightarrow \mu ^+\mu ^-)&=\frac{\tau _{B_q}G_F^4 m_W^4}{8\pi ^5}|V_{tb}V_{tq}^*|^2f_{B_q}^2m_{B_q}m_\mu ^2 \nonumber \\&\quad \;\times \sqrt{1-\frac{4 m_\mu ^2}{m_{B_q}^2}}(|P|^2+|S|^2), \end{aligned}$$

(2.3)

where $m_{B_q}$, $\tau _{B_q}$ and $f_{B_q}$ denote the mass, mean lifetime, and decay constant of $B_q$ meson, respectively. The short-distance contributions S and P are defined as

$$\begin{aligned} P&=C_{10}+\frac{m_{B_q}^2}{2m_W^2}\left( \frac{m_b}{m_b+m_q}\right) C_P,\nonumber \\ S&= \sqrt{1-\frac{4m_\mu ^2}{m_{B_q}^2}}\frac{m_{B_q}^2}{2m_W^2}\left( \frac{m_b}{m_b+m_q}\right) C_S. \end{aligned}$$

(2.4)

As discussed in the following section, there is no BSM phase in the 2HDMs with $Z_2$ symmetry. Therefore, we only consider the case that both S and P are real in this paper.

As pointed out in Ref. [75], the measured branching ratio of $B_q \rightarrow \mu ^+ \mu ^-$ should be the time-integrated one, denoted by $\overline{\mathcal B}(B_q \rightarrow \mu ^+ \mu ^-)$. In order to compare with the experimental measurements, the sizable effect of $B_s$–$\bar{B}_s$ oscillations should be taken into account [75, 76], and one has

$$\begin{aligned} \overline{\mathcal B}(B_s\rightarrow \mu ^+\mu ^-)&=\left( \frac{1+\mathcal A_{\Delta \Gamma }y_s}{1-y_s^2}\right) \mathcal B(B_s\rightarrow \mu ^+\mu ^-),\nonumber \\ \overline{\mathcal B}(B_d\rightarrow \mu ^+\mu ^-)&\approx \mathcal B(B_d\rightarrow \mu ^+\mu ^-), \end{aligned}$$

(2.5)

where the mass-eigenstate rate asymmetry $\mathcal A_{\Delta \Gamma }$ can be expressed as

$$\begin{aligned} \mathcal A_{\Delta \Gamma }=\frac{|P|^2-|S|^2}{|P|^2+|S|^2}. \end{aligned}$$

(2.6)

The observable $\mathcal A_{\Delta \Gamma }$ is independent of the branching ratio of $B_s \rightarrow \mu ^+ \mu ^-$ and provides complementary information on the short-distance structure of this decay. In the SM, $\mathcal A_{\Delta \Gamma }=+1$.

Following Ref. [75], it is convenient to introduce the ratio

$$\begin{aligned} R&\equiv \frac{\overline{\mathcal B}(B_s\rightarrow \mu ^+ \mu ^-)}{\mathcal B(B_s\rightarrow \mu ^+ \mu ^-)_\mathrm{SM}}\nonumber \\ {}&=\left( \frac{|P|^2}{1-y_s}+\frac{|S|^2}{1+y_s}\right) \frac{1}{|S_\mathrm{SM}|^2+|P_\mathrm{SM}|^2}, \end{aligned}$$

(2.7)

where both hadronic uncertainties and CKM matrix elements are canceled out.

3 $B_s \rightarrow \mu ^+ \mu ^-$ in the 2HDMs with $Z_2$ symmetry

In the 2HDMs with $Z_2$ symmetry, $b \rightarrow s \mu ^+ \mu ^-$ processes receive contributions from box diagrams with charged Higgs and penguin diagrams with Z boson and neutral Higgs bosons. The Wilson coefficient $C_{10}$ has been calculated in the type-II 2HDM [54]. For $C_S$ and $C_P$, only the leading contributions in the large $\tan \beta$ limit have been computed in the type-II model [61–63]. However, the remaining contributions could be important for some specific $\tan \beta$ values in the other types of 2HDMs. In this section, we first of all give a brief introduction to the 2HDMs with $Z_2$ symmetry, and then show that the Wilson coefficients could be derived explicitly from the recent full one-loop results of the A2HDM [60].

3.1 2HDMs with $Z_2$ symmetry

The 2HDM extends the SM Higgs sector with an additional scalar doublet. With the two Higgs doublets $\Phi _1$ and $\Phi _2$, the CP-conversing 2HDM potential with a softly broken $Z_2$ symmetry reads [19]

$$\begin{aligned} V&=+m_1^2\Phi _1^\dagger \Phi _1 +m_2^2\Phi _2^\dagger \Phi _2-m_3^2(\Phi _1^\dagger \Phi _2+\Phi _2^\dagger \Phi _1)\nonumber \\&\quad +\frac{\lambda _1}{2}(\Phi _1^\dagger \Phi _1)^2+\frac{\lambda _2}{2}(\Phi _2^\dagger \Phi _2)^2+\lambda _3(\Phi _1^\dagger \Phi _1)(\Phi _2^\dagger \Phi _2)\nonumber \\&\quad +\lambda _4(\Phi _1^\dagger \Phi _2)(\Phi _2^\dagger \Phi _1)\nonumber \\&\quad +\frac{\lambda _5}{2}[(\Phi _1^\dagger \Phi _2)^2+(\Phi _2^\dagger \Phi _1)^2], \end{aligned}$$

(3.1)

where $m_3^2(\Phi _1^\dagger \Phi _2 + \Phi _2^\dagger \Phi _1)$ is a soft $Z_2$ symmetry breaking term and the parameters $m_{1-3}$ and $\lambda _{1-5}$ are real. The two Higgs doublets $\Phi _1$ and $\Phi _2$ can be generally parameterized as

$$\begin{aligned} \Phi _i= \left( \begin{array}{c} \omega _i^+\\ \frac{1}{\sqrt{2}}(v_i+h_i-iz_i) \end{array} \right) , \end{aligned}$$

(3.2)

where the two vacuum expectation values (VEVs) $v_1$ and $v_2$ are real and positive. From the vacuum condition [77]

$$\begin{aligned} \begin{aligned} m_3^2 v_2 -m_1^2 v_1 -\frac{1}{2} \lambda _1 v_1^3 -\frac{1}{2}\lambda _{345} v_1 v_2^2&=0, \\ m_3^2 v_1 -m_2^2 v_2 -\frac{1}{2} \lambda _2 v_2^3 -\frac{1}{2}\lambda _{345} v_1^2 v_2&=0, \end{aligned} \end{aligned}$$

(3.3)

they can be expressed as other parameters in the Higgs potential, where $\lambda _{345}=\lambda _3+\lambda _4+\lambda _5$ is defined. By introducing the VEV v ($v=v_\mathrm{SM}=246\,\mathrm{GeV}$), the mixing angle $\beta$ and the soft $Z_2$ symmetry breaking parameter M as $v_1=v\cos \beta$, $v_2=v\sin \beta$ and $M^2=m_3^2/s_\beta c_\beta$, we can use $(v,\beta ,M,\lambda _{1-5})$ as independent 2HDM potential parameters.

Physical Higgs states are obtained by the following rotations:

$$\begin{aligned}&\left( \begin{array}{ll} h_1\\ h_2 \end{array} \right) =R(\alpha ) \left( \begin{array}{ll} H\\ h \end{array} \right) ,\quad \left( \begin{array}{ll} z_1\\ z_2 \end{array} \right) =R(\beta ) \left( \begin{array}{ll} G^0\\ A \end{array} \right) ,\quad \nonumber \\&\left( \begin{array}{l@{\quad }l} \omega _1^+\\ \omega _2^+ \end{array} \right) =R(\beta ) \left( \begin{array}{ll} G^+\\ H^+ \end{array} \right) , \end{aligned}$$

(3.4)

where the rotation matrix is given by

$$\begin{aligned} R(\theta )= \left( \begin{array}{l@{\quad }l} \cos \theta &{}-\sin \theta \\ \sin \theta &{}\cos \theta \end{array} \right) . \end{aligned}$$

(3.5)

The mixing angle $\alpha$ is determined by the Higgs potential of Eq. (3.1) [77],

$$\begin{aligned} \tan 2\alpha&=\frac{(M^2-\lambda _{345}v^2)s_{2\beta }}{(M^2-\lambda _1 v^2)c_\beta ^2 - (M^2-\lambda _2 v^2)s_\beta ^2}. \end{aligned}$$

(3.6)

In the 2HDM with $Z_2$ symmetry, the physical Higgs spectrum consists of five degrees of freedom: two charged scalars $H^\pm$, two CP-even neutral scalars h and H, and one CP-odd neutral scalar A. The quartic couplings $\lambda _i$ in the Higgs potential can be expressed in terms of their masses as [77]

$$\begin{aligned} \lambda _1&=\frac{1}{v^2 c_\beta ^2} (-s_\beta ^2 M^2 + s_\alpha ^2 m_h^2 + c_\alpha ^2 m_H^2 ) , \nonumber \\ \lambda _2&=\frac{1}{v^2 s_\beta ^2} (-c_\beta ^2 M^2 + c_\alpha ^2 m_h^2 + s_\alpha ^2 m_H^2 ) , \nonumber \\ \lambda _3&=-\frac{M^2}{v^2} +2\frac{m_{H^\pm }^2}{v^2} +\frac{1}{v^2}\frac{s_{2\alpha }}{s_{2\beta }} (m_H^2-m_h^2),\nonumber \\ \lambda _4&=\frac{1}{v^2} (M^2 + m_A^2 -2 m_{H^\pm }^2), \nonumber \\ \lambda _5&=\frac{1}{v^2} (M^2-m_A^2). \end{aligned}$$

(3.7)

Therefore, the eight parameters in the Higgs potential $m_{1-3}$ and $\lambda _{1-5}$ can be rewritten equivalently by the four physical Higgs masses $m_h$, $m_H$, $m_A$, $m_{H^\pm }$, the two mixing angles $\alpha$ and $\beta$, the VEV $v=v_\mathrm{SM}$, and the $Z_2$ symmetry breaking parameter M. In the case of $\lambda _1=\lambda _2$, which is considered in Refs. [62, 63], M can be eliminated and the 2HDM potential parameters can be expressed by seven parameters $(\alpha ,\beta ,v,m_h, m_H, m_A,m_{H^\pm })$ as

$$\begin{aligned} \lambda _1=\lambda _2&=\frac{1}{2v^2}(m_h^2+m_H^2)-\frac{1}{2v^2} \frac{c_{2\alpha }}{c_{2\beta }} (m_h^2-m_H^2),\nonumber \\ \lambda _3&=-\frac{1}{2v^2}(m_h^2+m_H^2-4m_{H^\pm }^2)\nonumber \\&\quad -\frac{1}{2v^2}(m_h^2-m_H^2)\left( \frac{c_{2\alpha }}{c_{2\beta }}+2\frac{s_{2\alpha }}{s_{2\beta }}\right) ,\nonumber \\ \lambda _4&=\frac{1}{v^2}(m_A^2-2m_{H^\pm }^2)+\frac{1}{2v^2}(m_h^2+m_H^2)\nonumber \\&\quad +\frac{1}{2v^2}\frac{c_{2\alpha }}{c_{2\beta }}(m_h^2-m_H^2),\nonumber \\ \lambda _5&=-\frac{m_A^2}{v^2}+\frac{1}{2v^2}(m_h^2+m_H^2) +\frac{1}{2v^2}\frac{c_{2\alpha }}{c_{2\beta }}(m_h^2-m_H^2),\nonumber \\ M^2&=\frac{1}{2}(m_h^2+m_H^2) + \frac{1}{2}\frac{c_{2\alpha }}{c_{2\beta }}(m_h^2-m_H^2). \end{aligned}$$

(3.8)

In the interaction basis, the general Yukawa Lagrangian of the 2HDM can be written as

$$\begin{aligned} - \mathcal L_Y&= \bar{Q}_L( Y_1^d\Phi _1+Y_2^d\Phi _2) d_R +\bar{Q}_L(Y_1^u \tilde{\Phi }_1+Y_2^u \tilde{\Phi }_2)u_R \nonumber \\&\quad +\bar{L}_L(Y_1^\ell \Phi _1+Y_2^\ell \Phi _2) e_R +\text {H.c.}, \end{aligned}$$

(3.9)

where $\tilde{\Phi }_i=i\sigma _2 \Phi _i^*$, $Q_L$ and $L_L$ denote the SM quark and lepton doublets, and $u_R$, $d_R$, and $e_R$ are the right-handed up-type quark, down-type quark, and lepton singlet, respectively. The Yukawa coupling matrices $Y_i^{u,d,\ell }$ are $3\times 3$ complex matrices in flavor space.

Table 1 Charge assignments of the $Z_2$ symmetry in the four types of 2HDM

In order to avoid tree-level FCNC, a discrete $Z_2$ symmetry is introduced [49]. All the possible nontrivial $Z_2$ charge assignments are listed in Table 1, which define the four well-known types of 2HDM, i.e. type-I, type-II, type-X, and type-Y. In the mass-eigenstate basis, the Yukawa interactions can be written in the form

$$\begin{aligned} -{\mathcal L}_Y&=+\sum _{f=u,d,\ell } \left[ \phantom {\left. +\left( \frac{m_f}{v}\xi _h^f \bar{f} fh+\frac{m_f}{v}\xi _H^f \bar{f} fH-i\frac{m_f}{v}\xi _A^f \bar{f} \gamma _5fA \right) \right] }m_f \bar{f} f\right. \nonumber \\&\quad \left. +\left( \frac{m_f}{v}\xi _h^f \bar{f} fh+\frac{m_f}{v}\xi _H^f \bar{f} fH-i\frac{m_f}{v}\xi _A^f \bar{f} \gamma _5fA \right) \right] \nonumber \\&\quad +\frac{\sqrt{2}}{v}\bar{u} (m_u V \xi _A^u P_L+ V m_d\xi _A^d P_R )d H^+\nonumber \\&\quad +\frac{\sqrt{2}m_\ell \xi _A^\ell }{v}\bar{\nu }_L \ell _R H^+ +\text {H.c.}, \end{aligned}$$

(3.10)

where $P_{L,R}=(1\mp \gamma _5)/2$. The Yukawa couplings $\xi _{h,H,A}^f$ in the four types of 2HDM are listed in Table 2. In addition, the couplings of the light CP-even Higgs boson h to gauge bosons $W^+W^-$ or ZZ can be written as $g_{hVV}=\sin (\beta -\alpha )g_{hVV}^\mathrm{SM}$, which is normalized to the corresponding couplings of the SM Higgs boson $g_{hVV}^\mathrm{SM}$ [19].

Recently, the LHC Run I data confirm the SM Higgs-like nature of the $125\,\mathrm{GeV}$ boson discovered at the LHC [3–6]. If the light CP-even Higgs h in the 2HDM is identified with the observed $125\,\mathrm{GeV}$ boson, global fits to the LHC Higgs data suggest that all four types of 2HDM should lie close to the so-called alignment limit [78–84]

$$\begin{aligned} \sin (\beta -\alpha )=1, \end{aligned}$$

(3.11)

where both the Yukawa and the gauge couplings of h are identical to the values of the SM Higgs boson. From Eqs. (3.3) and (3.6), the alignment limit can be achieved when the quartic couplings in the Higgs potential satisfy [85–87]

$$\begin{aligned} \tan ^2\beta =\frac{\lambda _1-\lambda _{345}}{\lambda _2-\lambda _{345}} , \quad \mathrm{or}\quad \lambda _1=\lambda _2=\lambda _{345}. \end{aligned}$$

(3.12)

For recent studies on the alignment limit in the 2HDM, we refer to Refs. [85, 86].

Table 2 Yukawa couplings in the four types of 2HDM

Since the 2HDMs with $Z_2$ symmetry are particular cases of the A2HDM [48], there exists a one-to-one correspondence for Yukawa couplings between these two models. However, the correspondence is not so straightforward for Higgs cubic couplings. Unlike the 2HDMs with $Z_2$ symmetry, the A2HDM potential is usually defined in the so-called “Higgs basis” [88], in which only one Higgs doublet gets a nonzero VEV. Therefore, the parameter $\tan \beta$ defined in the NFC 2HDMs is not a physical parameter in the A2HDM [89].

3.2 $B_s \rightarrow \mu ^+ \mu ^-$ in the 2HDMs with $Z_2$ symmetry

In both the A2HDM and the NFC 2HDMs, $B_s \rightarrow \mu ^+ \mu ^-$ decay is induced by gauge boson Z, Goldstone boson $G^0$, and Higgs bosons $\varphi \equiv \lbrace h, H, A \rbrace$ penguin diagrams, as well as box diagrams mediated with $W^\pm$, $H^\pm$, and $G^\pm$. At one-loop level, their contributions to the Wilson coefficients are divided into the following different parts:

$$\begin{aligned} C_{10}= & {} ( C_{10}^{Z,\,\mathrm{SM}} + C_{10}^\mathrm{box,\,\mathrm{SM}}) + ( C_{10}^{Z,\,\mathrm{2HDM}}),\nonumber \\ C_S= & {} (C_S^\mathrm{box,\,SM} + C_S^\mathrm{box,\, 2HDM} + C_S^{\varphi ,\,\mathrm{2HDM}} ),\nonumber \\ C_P= & {} (C_P^\mathrm{box,\, SM}+C_P^{Z,\, \mathrm{SM}} + C_P^{G,\,\mathrm{SM}})\\&+\, (C_P^{Z,\,\mathrm{2HDM}} + C_P^{G,\,\mathrm{2HDM}} )\nonumber \\&+\, (C_P^\mathrm{box,\, 2HDM} + C_P^{\varphi ,\,\mathrm{2HDM}} ),\nonumber \end{aligned}$$

(3.13)

where each part in the parentheses is gauge invariant. This gauge invariance is validated by the actual calculation in both the Feynman and the unitary gauges in the A2HDM [60]. The Wilson coefficients labeled with “SM” denote the contributions from the diagrams involved with only the SM fields (with the Goldstone bosons but not the Higgs boson), whose expressions are given in Appendix A. Those with “2HDM” contain the Higgs contributions. For simplicity, their explicit expressions are given in the unitarity gauge in the following, where the Goldstone boson contributions are absent.

The Higgs bosons affect the box and Z penguin diagrams with Yukawa interactions. Their contributions to Wilson coefficients in the NFC 2HDMs can easily be obtained from the A2HDM results with replacement of the Yukawa couplings,

$$\begin{aligned}&C_{S,P,\,\mathrm Unitary}^\mathrm{box,\, 2HDM}=C_{S,P,\, \mathrm Unitary}^\mathrm{box,\, A2HDM}\Bigr |_{(\varsigma _u,\varsigma _d,\varsigma _\ell )\rightarrow (-\xi _A^u,\xi _A^d,\xi _A^\ell )},\\&C_{10,P,\,\mathrm Unitary}^{Z,\,\mathrm{2HDM}}=C_{10,P,\, \mathrm Unitary}^{Z\,\mathrm penguin, \, A2HDM} \Bigr |_{(\varsigma _u,\varsigma _d,\varsigma _\ell )\rightarrow (-\xi _A^u,\xi _A^d,\xi _A^\ell )}.\nonumber \end{aligned}$$

(3.14)

For making this article self-contained, we present the Wilson coefficients after the correspondences made in Appendix A.

The Higgs penguin diagrams involve Yukawa couplings as well as Higgs–gauge couplings and Higgs cubic couplings. Therefore, their Wilson coefficients cannot be derived from the A2HDM results so straightforwardly as in the box and Z penguin diagrams, as discussed in previous section. Since the A2HDM Wilson coefficients are given for individual Higgs penguin diagrams in Ref. [60], we use the following approach. For every Higgs penguin diagram in the NFC 2HDMs, its contribution is derived from the A2HDM results with the replacement of the Higgs–gauge vertex and the triple Higgs vertex. Then the total contributions to the Wilson coefficients are obtained,

$$\begin{aligned} C_{S,\,\mathrm Unitary}^{\varphi ,\,\mathrm{2HDM}}&=+\frac{x_t\xi _h^\ell }{2x_h}\left( -s_{\alpha -\beta }g_1^{(a)}+c_{\alpha -\beta }g_2^{(a)}+\frac{2v^2}{m_W^2}\lambda _{H^+H^-}^hg_0\right) \nonumber \\&\quad +\frac{x_t\xi _H^\ell }{2x_H}\left( +c_{\alpha -\beta }g_1^{(a)}+s_{\alpha -\beta }g_2^{(a)}+\frac{2v^2}{m_W^2}\lambda _{H^+H^-}^Hg_0\right) , \nonumber \\ C_{P,\,\mathrm Unitary}^{\varphi ,\,\mathrm{2HDM}}&=-\frac{x_t\xi _A^\ell }{2x_A}g_3^{(a)}, \end{aligned}$$

(3.15)

where $x_t = m_t^2/m_W^2$, $x_{h,H,A}=m_{h,H,A}^2/m_W^2$, the functions $g_{0-3}^{(a)}\equiv g_{0-3}^{(a)}\bigl (x_t,x_{H^\pm },-\xi _A^u,\xi _A^d \bigr )$ defined in Eq. (A.4), and the Higgs cubic couplings are defined as

$$\begin{aligned}&\left[ \begin{array}{l} \lambda ^h_{H^+H^-}\\ \lambda ^H_{H^+H^-}\\ \lambda ^A_{H^+H^-} \end{array} \right] \nonumber \\&\quad =\frac{1}{2v^2s_{2\beta }} \left[ \begin{array}{ccc} (m_h^2-2m_{H^\pm }^2)c_{\alpha -3\beta }+(-4M^2+3m_h^2+2m_{H^\pm }^2)c_{\alpha +\beta }\\ (m_H^2-2m_{H^\pm }^2)s_{\alpha -3\beta }+(-4M^2+3m_H^2+2m_{H^\pm }^2)s_{\alpha +\beta }\\ 0 \end{array}\right] , \end{aligned}$$

(3.16)

where the soft $Z_2$ symmetry breaking parameter M has been defined in Sect. 3.1.

In the literature [61–63], it is found that the Wilson coefficients can receive large $\tan \beta$ enhancement only in the type-II 2HDM and the branching ratio with large $\tan \beta$ depends only on the Higgs masses $m_{H^\pm }$, $m_H$, $m_h$ and the mixing angle $\alpha$. However, as shown by Eqs. (3.15) and (3.16), a term proportional to $M^2/m_H^2$ in our full one-loop Wilson coefficient $C_S$ is also enhanced by $\tan ^2\beta$, which comes from the heavy Higgs H penguin diagrams mediated by charged Higgs bosons. Using the parameter $m_3$ in the Higgs potential of Eq. (3.1) directly, this term is proportional to $m_3^2/m_H^2$ and enhanced by $\tan ^3\beta$. This M dependent term has not been considered yet in the previous studies in the literature. Therefore, its effects are worthy of a detailed investigation.

The soft $Z_2$ symmetry breaking parameter M is associated with the spontaneous CP breaking [17, 90–92] and characterizes the masses of all the Higgs bosons [77]. This parameter enters the $B_s \rightarrow \mu ^+ \mu ^-$ decays through the Higgs penguin diagrams. However, it is found that the M term cannot make more significant contributions than other terms of the Wilson coefficient $C_S$. Here, we would choose h as the Higgs boson discovered by ATLAS [1] and CMS [2] and take the alignment limit $\beta -\alpha =\pi /2$, which is favored by the current 2HDM fits [78–84]. Then the cubic couplings in Eq. (3.16) read

$$\begin{aligned} \left[ \begin{array}{l} \lambda ^h_{H^+H^-}\\ \lambda ^H_{H^+H^-}\\ \lambda ^A_{H^+H^-} \end{array}\right] \doteq \frac{1}{v^2} \left[ \begin{array}{cc} -2M^2+2m_{H^\pm }^2+m_h^2\\ \cot {2\beta } (2M^2-2m_H^2)\\ 0 \end{array}\right] . \end{aligned}$$

(3.17)

Focusing on the coupling $\lambda _{H^+ H^-}^h$, it can be seen from Eqs. (2.4) and (3.17) that large contributions from this coupling would require $|M^2-m_{H^\pm }^2 |/v^2\gg m_W^2/m_B^2$. However, we know $|M^2-m_{H^\pm }^2 |/ v^2=|\lambda _4+\lambda _5 |/2 < 4\pi$ from the 2HDM vacuum condition [77] and perturbativity [93]. It is also noted that the Higgs penguin diagrams can be enhanced by very large $\tan \beta$ or $\cot \beta$. In all the four types of 2HDMs, the $\lambda _{H^+H^-}^h$ contributions could be enhanced by large $\cot ^2\beta$. In practice, $\cot \beta \gtrsim 3$ has been excluded by the perturbativity [93]. Similarly, the coupling $\lambda _{H^+H^-}^H$ can make a large contribution if $M^2/m_H^2 \gg m_W^2/m_B^2$. Among the four models, this contribution is enhanced by $\tan ^2\beta$ only in type-II 2HDM. However, the ratio $M^2/m_H^2$ still suffers from the theoretical constraints, which will be discussed with numerical results in the following section.

Although the effects from the operators $\mathcal O_S$ and $\mathcal O_P$ are suppressed by $m_B^2/m_W^2$, these two scalar operators can make significant contributions in the two parameter regions: (i) in the type-II 2HDM, since both $C_S$ and $C_P$ contain $\tan \beta$ enhanced terms, the effects of the scalar operators are enhanced in the parameter space with large $\tan \beta$. (ii) The contributions from the CP-odd Higgs penguin diagrams are inversely proportional to the mass of the CP-odd Higgs boson A. Thus, the Wilson coefficient $C_P$ becomes much more significant in the region with small values of $m_A$ ¹ in all the four 2HDMs.

In the particular case of the type-II 2HDM, our result of $C_{10}$ agrees with the one calculated in Ref. [54]. For the Wilson coefficients $C_S$ and $C_P$ in the 2HDM, the calculations have been performed by various groups [52, 53, 61–63, 94–97]. The latest results are presented in these three papers [61–63], where the 2HDM contributions are computed in the type-II model in some specific cases. In Ref. [61], the Wilson coefficients are calculated in large $\tan \beta$ limit, i.e., only $\tan ^2\beta$ enhanced terms are kept. However, the Higgs penguin diagrams with trilinear $hH^+ H^-$ and $HH^+ H^-$ couplings are not considered. In Refs. [62, 63],² after including these penguin diagrams, the calculations are performed again in the large $\tan \beta$ limit but with the assumption $\lambda _1 = \lambda _2$ for the couplings in the Higgs potential.³ Considering only terms proportional to $\tan ^2\beta$, our result agrees with the one of Ref. [61] in the case of $\lambda _{H^+H^-}^h=\lambda _{H^+H^-}^H=0$, and those of Refs. [62, 63] in the case of $\lambda _1=\lambda _2$. Generally, the 2HDM contains eight free parameters, i.e., $m_{1-3}$ and $\lambda _{1-5}$ in the Higgs potential of Eq. (3.1). They can be rewritten equivalently in terms of the Higgs masses $m_h$, $m_H$, $m_A$, $m_{H^\pm }$, the mixing angles $\alpha$ and $\beta$, the parameter M, and the VEV $v=v_\mathrm{SM}$. If the condition $\lambda _1=\lambda _2$ is assumed, M can be expressed by the other parameters, as shown in Eq. (3.8). It is the reason why terms depending on the $Z_2$ symmetry breaking parameter M were absent in the previous calculations [61–63] but are present in this paper.

4 Numerical analysis

Searches for $B_{s,d} \rightarrow \mu ^+ \mu ^-$ decays have been performed at the BaBar, Belle, and Tevatron (for a review, see Ref. [98]). At the LHC, measurements by CMS [100] and LHCb [101] collaborations with the full data of LHC Run I have resulted in the averaged values for the time-integrated branching ratios [102]

$$\begin{aligned} \overline{\mathcal B}(B_s \rightarrow \mu ^+ \mu ^-)&= 2.8_{-0.6}^{+0.7} \times 10^{-9},\\ \overline{\mathcal B}(B_d \rightarrow \mu ^+ \mu ^-)&= 3.9_{-1.4}^{+1.6} \times 10^{-10}, \end{aligned}$$

where the errors are dominated by the statistical uncertainties and expected to be significantly reduced in the near future. Both of them are in good agreement with the latest updated SM predictions [57], $\overline{\mathcal B}(B_s \rightarrow \mu ^+ \mu ^-)=(3.65\pm 0.23)\times 10^{-9}$ and $\overline{\mathcal B}(B_d \rightarrow \mu ^+ \mu ^-)=(1.06\pm 0.09)\times 10^{-10}$, in which the NLO EW [58] and the NNLO QCD [59] corrections have been included. Thus, strong constraints on the 2HDM parameters are expected.

Fig. 1

a The M dependence of the branching ratio of $B_s \rightarrow \mu ^+ \mu ^-$ in the type-II 2HDM for $\tan \beta =20$ (solid) and $\tan \beta =40$ (dashed). The SM prediction (dotted) and $2\sigma$ experimental range (dot–dashed) are also shown. b Allowed regions of the parameter space $(\tan \beta ,m_{H^\pm })$ from $\overline{\mathcal B} (B_s \rightarrow \mu ^+ \mu ^-)$ for the four types of 2HDM

In the NFC 2HDMs, the relevant parameters are the two mixing angles $\alpha$ and $\beta$, four Higgs mass parameters $m_{H^\pm }$, $m_h$, $m_H$, and $m_A$. In the $B_{s,d} \rightarrow \mu ^+ \mu ^-$ decays, the $Z_2$ symmetry breaking parameter M also enters into the decay amplitude and is independent from these parameters. As discussed in Refs. [103, 104], we choose the light neutral Higgs h in the 2HDM as the SM Higgs observed at the LHC and adopt the alignment limit $\sin (\beta -\alpha )=1$. Then the model parameters are reduced to $(m_H,m_A,m_{H^\pm },M,\tan \beta )$. As discussed in Ref. [64], we shall restrict these parameters in the following ranges:

$$\begin{aligned}&m_H \in [m_h,1000]\,\mathrm{GeV}, \quad m_{H^\pm },m_A,M\in [1,1000]\,\mathrm{GeV},\quad \nonumber \\&\tan \beta \in [0.1, 100]. \end{aligned}$$

(4.1)

Starting from these parameter spaces, we will start our numerical scan.

In the numerical analysis, we impose experimental constraints in the same way as in Ref. [64]. To constrain the 2HDM parameters, we have taken into account (i) flavor processes: $B_{s,d}$–$\bar{B}_{s,d}$ mixing, $\bar{B} \rightarrow X_s \gamma$, $B \rightarrow \tau \nu$ and $B_{s,d}\rightarrow \mu ^+ \mu ^-$ decays, (ii) direct searches for Higgs bosons at LEP [67], Tevatron [68, 69] and LHC [70, 71], both of which have been discussed in detail in our previous work [64]. Additionally, we also consider the oblique parameter $\Delta \rho$ in the EW precision measurement [105–110] and require the couplings $\lambda _{1-5}$ to satisfy (iii) theoretical constraints: perturbativity [93], tree-level vacuum stability [90, 111, 112] and perturbative unitarity [19, 113] (see Ref. [114] for the expressions).

Fig. 2

a Combined constraints on the parameter space of the four types of 2HDM, plotted in the $(\tan \beta ,m_{H^\pm })$ plane. b Correlations between R and $\mathcal A_{\Delta \Gamma }$ in the four types of 2HDM

For $B_s \rightarrow \mu ^+ \mu ^-$ decay, both the NNLO QCD and the NLO EW corrections in the SM and the full one-loop contributions in the 2HDM are included. As discussed in Sect. 3, the effects of the soft $Z_2$ symmetry breaking parameter M can be enhanced by large $\tan \beta$ in the type-II 2HDM. The M dependence of the branching ratio $\overline{\mathcal B}(B_s \rightarrow \mu ^+ \mu ^-)$ is shown in Fig. 1a in the type-II 2HDM for various $\tan \beta$ and $m_H$ values. As expected, the effects of M become significant when the two ratios $M^2/m_H^2$ and $\tan \beta$ are large. However, it is found that the theoretical constraints from perturbativity, vacuum stability, and perturbative unitarity have put the bound $M^2/m_H^2 \lesssim 1$ (and $M \lesssim 1\,\mathrm{TeV}$) in the parameter space of Eq. (4.1). Therefore, the soft $Z_2$ symmetry breaking parameter M cannot make more significant effects than the other $\tan \beta$ enhanced terms in $C_S$ and $C_P$.

After considering the current experimental data, the allowed parameter spaces of all the four 2HDMs are obtained. Since the constraints from $B_d \rightarrow \mu ^+ \mu ^-$ appear to be more or less weaker than those from $B_s \rightarrow \mu ^+ \mu ^-$, we only show the results from the latter one, which are plotted in the $(\tan \beta ,m_{H^\pm })$ plane in Fig. 1b. Compared to our previous results [64], the parameter space with small $\tan \beta$ is excluded for all the four types of 2HDMs. This change is caused by the contributions with small $\tan \beta$ neglected in the previous calculations [61–63] but included in the present full one-loop computation as discussed in Sect. 3. For the large $\tan \beta$ region, only the type-II model is bounded, which is in agreement with our previous result but still weaker than the one from $\mathcal B (B \rightarrow \tau \nu )$.

Combining all the constraints aforementioned, we obtain the survived parameter space of all the four types 2HDMs, as an update of our previous results [64], which is shown in the $(\tan \beta , m_{H^\pm })$ plane in Fig. 2a. It is found that the small $\tan \beta$ region is restricted for all the four models by $B_s$–$\bar{B}_s$ mixing and $B \rightarrow X_s \gamma$, while the large $\tan \beta$ region is constrained only in the type-II 2HDM by $B\rightarrow \tau \nu$ and $B_s \rightarrow \mu ^+ \mu ^-$ decays. Compared to our previous results, the current constraints on the large $\tan \beta$ region in the type-II 2HDM are more stringent. This is mainly because the theoretical constraints are included in the current analysis.

In these constrained parameter spaces of the four 2HDMs, the correlations between the observables $\mathcal A_{\Delta \Gamma }$ and R defined in Eqs. (2.6) and (2.7) are reevaluated, which are presented in Fig. 2b. Unlike our previous results [64], the correlations in the four different types of 2HDMs are almost indistinguishable. The allowed ranges of R are the same for all the four models, while $\mathcal A_{\Delta \Gamma }$ can deviate slightly from the SM prediction only in the type-II 2HDM. It is found that the difference from our previous results is mainly caused by the theoretical constraints and the new full one-loop Wilson coefficients considered in the current analysis. In the type-II 2HDM, the bounds on $\tan \beta$ are more stringent compared to our previous results as discussed above. Thus, the allowed range of $C_S$ is restricted more stringently in the current analysis. In this case, $\mathcal A_{\Delta \Gamma }$ can deviate from the SM prediction very tiny, which can be seen from Eq. (2.6). As discussed in Sect. 3, our results of the Wilson coefficients can also be applied to the small $\tan \beta$ region in all the four models, while some terms are not included in $C_P$ used in our previous analysis. In the case of small $m_A$, $C_P$ is enhanced and these terms make the allowed regions of R in the type-I and type-Y 2HDMs as large as the one in the type-X model. Meanwhile, the value of R is almost independent of $\mathcal A_{\Delta \Gamma }$ in the type-II 2HDM.

5 Conclusion

In this paper, we have performed an updated analysis of the rare leptonic decay $B_s \rightarrow \mu ^+ \mu ^-$ in the 2HDM with a softly broken $Z_2$ symmetry. We have derived the full one-loop Wilson coefficients $C_{10}$, $C_S$, and $C_P$ from the recent A2HDM results [60], which can be applied to the contributions of all the four types of 2HDMs for both large and small $\tan \beta$ value. Our main conclusions are summarized as follows:

• Compared to $C_{10}$, the Wilson coefficients $C_S$ and $C_P$ are negligible in the entire 2HDM parameter space, except for large $\tan \beta$ in the type-II 2HDM or small CP-odd Higgs mass $m_A$ in the four models. In addition, only the Wilson coefficients $C_S$ and $C_P$ in the type-II 2HDM can be enhanced by large $\tan \beta$.

• The soft $Z_2$ symmetry breaking parameter M enters into the Higgs penguin diagrams and affects the Wilson coefficient $C_S$. The dominant contributions are proportional to $M^2/m_H^2$ and enhanced by $\tan ^2\beta$ in the type-II 2HDM, which have not been considered in the literature [61–63]. However, after combing the theoretical constraints from perturbativity, vacuum stability and perturbative unitarity, we have found that the parameter M cannot make more significant contributions than other terms in the Wilson coefficients.

• After imposing the experimental constraints, regions with small $\tan \beta$ are excluded for all the four types of 2HDM, which are quite different from our previous results [64]. As expected, large $\tan \beta$ region is only excluded in the type-II 2HDM.

As an update of our previous analysis [64], we have also investigated the possibility to distinguish the four types of 2HDM in light of the recent updated flavor physics data, the collider data from the direct searches for Higgs bosons and the theoretical progresses. The combined bounds on the 2HDM parameters have been derived for the four models. In the survived parameter regions, the correlations between $\mathcal A_{\Delta \Gamma }$ and R in all the four of the 2HDMs are almost indistinguishable from each other. In the 2HDMs with $Z_2$ symmetry, $\mathcal A_{\Delta \Gamma }$ can only have a very tiny deviation from the SM prediction, while R could deviate from the SM one sizably. This could be tested by the much more precise measurement of $B_s \rightarrow \mu ^+ \mu ^-$ at the LHC in the coming years.

The Wilson coefficients in the SM and the 2HDMs

In this appendix, we recapitulate the relevant expressions of the Wilson coefficients in the SM and the four types of the 2HDMs for completeness, which are obtained from Ref. [60].

In the SM, the one-loop Wilson coefficients of the scalar operators can be written as

$$\begin{aligned} \begin{aligned} C_S^\mathrm{SM}&=C_S^\mathrm{box,\,SM}+C_S^{h,\,\mathrm{SM}},\\ C_P^\mathrm{SM}&=C_P^\mathrm{box,\,\mathrm{SM}}+C_P^{Z,\,\mathrm{SM}}+C_P^{G,\,\mathrm{SM}}. \end{aligned} \end{aligned}$$

(A.1)

In the unitary gauge, their expressions read

$$\begin{aligned} C_{S,\,\mathrm Unitary}^\mathrm{box,\,SM}&=-\frac{x_t(x_t+1)}{48(x_t-1)^2} - \frac{(x_t-2)(3x_t^2-3x_t+1)}{24(x_t-1)^3}\, \ln x_t, \nonumber \\ C_{S,\,\mathrm Unitary}^{h,\,\mathrm{SM}}&=-\frac{3 x_t}{8 x_h}, \\ C_{P,\,\mathrm Unitary}^\mathrm{box,\,\mathrm{SM}}&=+\frac{x_t(71x_t^2-172 x_t-19)}{144(x_t-1)^3}\nonumber \\&\quad + \frac{x_t^4-12x_t^3+34 x_t^2-x_t-2}{24(x_t-1)^4} \ln x_t , \nonumber \\ C_{P,\,\mathrm Unitary}^{Z,\,\mathrm{SM}}&=+\frac{1}{12}\left[ \frac{x_t(18 x_t^3-137 x_t^2+262 x_t-95)}{6 (x_t-1)^3} \right. \nonumber \\&\quad \left. + \frac{8 x_t^4-11 x_t^3-15 x_t^2+12 x_t-2}{(x_t-1)^4}\ln x_t \right] \nonumber \\&\quad - \frac{s_W^2}{36}\left[ \frac{x_t(18 x_t^3-139 x_t^2+274 x_t-129)}{2 (x_t-1)^3} \right. \nonumber \\&\quad \left. + \frac{24 x_t^4-33 x_t^3-45 x_t^2+50 x_t-8}{(x_t-1)^4} \ln x_t\right] , \nonumber \end{aligned}$$

(A.2)

where $C_{S,\,\mathrm Unitary}^{h,\,\mathrm{SM}}$ denotes the contributions from the SM Higgs penguin diagrams. The other Wilson coefficients $C_{S,P,\,\mathrm Unitary}^\mathrm{box,\,\mathrm{SM}}$ and $C_{P,\,\mathrm Unitary}^{Z,\,\mathrm{SM}}$ are same in the SM and the 2HDMs.

In the four types of the 2HDMs, the various contributions in the Wilson coefficients of Eq. (3.13) are obtained by the replacement of the Yukawa couplings in Eq. (3.14), which are given in the unitary gauge,

$$\begin{aligned} C^\mathrm{box,\,\mathrm{2HDM}}_{S,\,\mathrm Unitary}=&-\frac{\xi _A^u\xi _A^\ell x_t}{8(x_{H^\pm }-x_t)}\left[ 1-\frac{x_{H^\pm }\ln (x_{H^\pm }/x_t)}{(x_{H^\pm }-x_t)}\right] \nonumber \\&-\xi _A^d\xi _A^\ell \frac{x_t\ln (x_{H^\pm }/x_t)}{4(x_{H^\pm }-x_t)}, \\ C_{P,\,\mathrm Unitary}^\mathrm{box,\,\mathrm{2HDM}}&= +\frac{\xi _A^u\xi _A^\ell x_t}{8(x_{H^\pm }-x_t)}\nonumber \\&\quad \times \left[ 1 +\frac{2x_t^2-x_{H^\pm }x_t-x_{H^\pm }}{(x_t-1)(x_{H^\pm }-x_t)}\ln x_t\right. \nonumber \\&\quad \left. + \frac{x_{H^\pm }(1-2x_t+x_{H^\pm })}{(x_{H^\pm }-1)(x_{H^\pm }-x_t)}\ln x_{H^\pm }\right] \nonumber \\&\quad +\xi _A^d\xi _A^\ell \frac{x_t\ln (x_{H^\pm }/x_t)}{4(x_{H^\pm }-x_t)}, \nonumber \\ C^{Z,\, \mathrm{2HDM}}_{10,\,\mathrm Unitary}&= + (\xi _A^u)^2\,\frac{x_t^2}{8}\left[ \frac{1}{x_{H^\pm }-x_t} - \frac{x_{H^\pm }\ln (x_{H^\pm }/x_t)}{(x_{H^\pm }-x_t)^2}\right] , \nonumber \\ C_{P,\,\mathrm Unitary}^{Z,\,\mathrm{2HDM}}&= +\frac{x_t}{4(x_{H^\pm }-x_t)^2}\,\nonumber \\&\quad \times \biggl \lbrace -\xi _A^d\xi _A^u\biggl [- \frac{x_t +x_{H^\pm }}{2} + \frac{x_t x_{H^\pm }}{x_{H^\pm }-x_t}\ln \frac{x_{H^\pm }}{x_t}\biggr ] \nonumber \\&\quad + (\xi _A^u)^2\biggl [\frac{x_{H^\pm }^2-8x_{H^\pm }x_t-17x_t^2}{36(x_{H^\pm }-x_t)}-x_t(x_{H^\pm }-x_t) \nonumber \\&\quad +\left( \frac{x_t^2(3x_{H^\pm }+x_t)}{6(x_{H^\pm }-x_t)^2}+x_tx_{H^\pm }\right) \ln \frac{x_{H^\pm }}{x_t}\biggr ]\biggr \rbrace \nonumber \\&\quad + \frac{s_W^2 x_t}{6(x_{H^\pm }-x_t)^2}\,\nonumber \\&\quad \times \biggl \lbrace -\xi _A^d\xi _A^u\,\biggl [\frac{5 x_t-3 x_{H^\pm }}{2}\nonumber \\&\quad + \frac{x_{H^\pm }(2x_{H^\pm }-3x_t)}{x_{H^\pm }-x_t}\ln \frac{x_{H^\pm }}{x_t}\biggr ] \nonumber \\&\quad - (\xi _A^u)^2\biggl [\Bigg (\frac{4 x_{H^\pm }^3-12 x_{H^\pm }^2 x_t + 9 x_{H^\pm }x_t^2 + 3 x_t^3}{6(x_{H^\pm }- x_t)^2}\nonumber \\&\quad +\frac{3}{2}x_tx_{H^\pm }\Bigg )\ln \frac{x_{H^\pm }}{x_t} \nonumber \\&\quad - \frac{17 x_{H^\pm }^2-64 x_{H^\pm }x_t +71 x_t^2}{36(x_{H^\pm }-x_t)} -\frac{3}{2}x_t(x_{H^\pm }-x_t)\biggr ]\biggr \rbrace .\nonumber \end{aligned}$$

(A.3)

The Higgs penguin contributions $C_{S,P,\,\mathrm Unitary}^{\varphi ,\,\mathrm{2HDM}}$ have been given in Eq. (3.15), where the functions $g_{0-3}^{(a)}$ are defined as

$$\begin{aligned}&g_0(x_t,x_{H^\pm },-\xi _A^u,\xi _A^d) =-\frac{1}{4x_{H^\pm }}\nonumber \\&\qquad \times \left[ -\xi _A^u\xi _A^d(f_1+f_2+f_3+1)+(\xi _A^u)^2\right. \nonumber \\&\qquad \times \left. \left( f_4-f_5-\frac{1}{4}\right) \right] , \nonumber \\&g_1^{(a)}(x_t,x_{H^\pm },-\xi _A^u,\xi _A^d) = -\frac{\,3\,}{4}\nonumber \\&\qquad -\xi _A^u\xi _A^d (f_1+f_2+f_3) + (\xi _A^u)^2 (f_4-f_5), \nonumber \\&g_2^{(a)}(x_t,x_{H^\pm },-\xi _A^u,\xi _A^d) = -(\xi _A^d)^2\xi _A^u f_1 +\xi _A^d(\xi _A^u)^2 ( f_3 + f_2)\nonumber \\&\qquad +(\xi _A^u)^3 (f_5 - f_4) +\xi _A^u (f_7 - f_6) +\xi _A^d f_1, \nonumber \\&g_3^{(a)}(x_t,x_{H^\pm },-\xi _A^u,\xi _A^d) = -(\xi _A^d)^2 \xi _A^u f_1 +\xi _A^d(\xi _A^u)^2 (f_3 - f_2)\nonumber \\&\qquad -(\xi _A^u)^3 ( f_5 + f_4) -\xi _A^u ( f_7 + f_6 ) +\xi _A^d f_1. \end{aligned}$$

(A.4)

Here the one-loop functions $f_i$ are abbreviated as $f_i\equiv f_i(x_t,x_{H^\pm })$ with the definitions

$$\begin{aligned} f_1(x,y )&= \frac{1}{2(y- x)^{\phantom {2}}}\, [x-y +y \ln y-x \ln x], \nonumber \\ f_2(x,y)&= \frac{1}{2(y-x)^{\phantom {2}}} \left[ x - \frac{y x}{y-x}(\ln y-\ln x)\right] , \nonumber \\ f_3(x,y)&=\frac{1}{2(y-x)^{\phantom {2}}} \left[ y-\frac{y^2 \ln y}{y-x}+\frac{x (2y- x) \ln x}{y-x}\right] , \nonumber \\ f_4(x,y)&=\frac{1}{4(y-x)^2} \left[ \frac{x \left( 3 y-x\right) }{2}-\frac{y^2 x}{y- x}(\ln y-\ln x)\right] , \\ f_5(x,y)&=\frac{1}{4(y-x)^2} \left[ \frac{x (y-3 x)}{2}-\frac{y x (y-2 x)}{y- x}(\ln y-\ln x)\right] , \nonumber \\ f_6(x,y)&=\frac{x\left( x^2-3yx+9y-5 x-2\right) }{8(x-1)^2(y-x)} +\frac{y \left( y x-3y+2 x\right) }{4 (y-1)(y- x)^2}\ln y\nonumber \\&\quad +\frac{y^2 \left( -2 x^3+6 x^2-9 x+2\right) +3 y x^2 (x^2-2 x+3)-x^2 \left( 2 x^3-3 x^2+3 x+1\right) }{4 (x-1)^3 (y- x)^2}\ln x, \nonumber \\ f_7(x,y)&=\frac{\left( x^2+x-8\right) }{8 (x-1)^2}-\frac{y (y+2) }{4 (y-1)(y-x)}\ln y\nonumber \\&\quad +\frac{y \left( x^3-3 x^2+3 x+2\right) +3 \left( x-2\right) x^2}{4 (x-1)^3(y-x)}\ln x.\nonumber \end{aligned}$$

(A.5)

It is noted that the divergence in the Higgs penguin diagrams at one-loop level is canceled by a FCNC local operator in the A2HDM [60]. In the 2HDMs with $Z_2$ symmetry, we find that the divergence automatically vanishes after adding all the Higgs penguin contributions.

For the four types of 2HDM, the values of the relevant Yukawa couplings are listed in Table 2. When deriving the expressions of the Higgs penguin diagrams in Eq. (3.15), the following identities have been used:

$$\begin{aligned} \begin{aligned} \xi _h^u&=-\cos (\alpha -\beta ) \xi _A^u - \sin (\alpha -\beta ),\quad \\ \xi _H^u&=-\sin (\alpha -\beta ) \xi _A^u + \cos (\alpha -\beta ),\\ \xi _h^d&=+\cos (\alpha -\beta ) \xi _A^d - \sin (\alpha -\beta ),\quad \\ \xi _H^d&=+\sin (\alpha -\beta ) \xi _A^d + \cos (\alpha -\beta ), \end{aligned} \end{aligned}$$

(A.6)

and

$$\begin{aligned} (\xi _A^u+\xi _A^d)(\xi _A^u\xi _A^d-1)=0, \end{aligned}$$

(A.7)

which can be obtained from Table 2. It should be noted that there is a freedom in the definitions of the functions $f_i$, since adding the LHS of Eq. (A.7) to Eq. (A.4) does not change $g_{0-3}^{(a)}$.