††institutetext: Institute for Theoretical Particle Physics, Karlsruhe
Institute of Technology (KIT),
Wolfgang-Gaede-Str. 1, D-76131 Karlsruhe, Germany

#  $B_{s}\to\mu^{+}\mu^{-}$ in a Two-Higgs-Doublet Model with flavour-changing
up-type Yukawa couplings

Martin S. Lang,  Ulrich Nierste<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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

##  1 Introduction

Two-Higgs-Doublet models (2HDMs) Branco:2011iw ; Gunion:1989we  are a popular
extension of the Standard Model (SM) due to their relative simplicity,
involving no additional particles apart from a second Higgs doublet. Moreover,
a strong motivation to study 2HDMs also comes from theories in which a second
Higgs doublet is required due to symmetry arguments, e.g. axion models in the
context of the strong CP puzzle Kim:1986ax ; Peccei:1977hh ; Peccei:1977ur  or
minimal supersymmetry Haber:1984rc . 2HDMs differ in the structure of Higgs-
fermion Yukawa couplings. The historically most favoured variant is the so-
called type-II 2HDM in which up- and down-type quarks couple to separate Higgs
doublets exclusively. In most of the 2HDMs, flavour-changing neutral current
processes such as the decay $B_{s}\to\mu^{+}\mu^{-}$ are loop-suppressed and
therefore small masses of the additional Higgs bosons are in principle
possible, an appealing feature during the early LHC searches. A general 2HDM
exhibits, however, a much richer Yukawa structure Crivellin:2012ye ;
Crivellin:2013wna , in which flavour-changing neutral Higgs-boson couplings
are possible, in case quarks couple to multiple Higgs doublets at the same
time.

In this paper, we will consider the rare decay $B_{s}\to\mu^{+}\mu^{-}$ in the
context of a 2HDM with more general Yukawa couplings. The model under
consideration will feature flavour-changing up-type Yukawa couplings of the
additional neutral Higgs bosons, most notably a charm-top transition which may
eliminate the Cabibbo-Kobayashi-Maskawa (CKM) suppression that is present in
both the SM and the 2HDM type-II models. We put special emphasis on the
leading corrections for large values of $\tan\beta$.

The outline of the paper is as follows: In section 2, we will introduce the
effective operators contributing to the low-energy weak
$B_{s}\to\mu^{+}\mu^{-}$ decay. The Yukawa and Higgs sectors of the 2HDM are
presented in section 3, while section 4 is dedicated to a short description of
the computational setup used for the evaluation Feynman diagrams. The limiting
case of the type-II 2HDM is introduced in section 5, with the additional
contributions from flavour-changing neutral-Higgs Yukawa couplings being
presented in section 6. Finally, we will discuss the phenomenology of such
models in section 7, before summarizing in section 8.

##  2 The decay $\mathbf{B_{s}\to\mu^{+}\mu^{-}}$

The typical momentum scale for $B_{s}$ decays is of order $M_{B_{s}}$ or
smaller, so that weak $B_{s}$ decays can be described by an effective theory
in which the heavy $W,Z$ bosons, the top quark, and the Higgs bosons of the
2HDM are integrated out. The resulting $|\Delta B|=1$ Hamiltonian $H_{\rm
eff}$ describes the interactions mediated by these heavy particles in terms of
dimension-6 operators changing the beauty quantum number $B$ by one unit. The
piece of $H_{\rm eff}$ relevant for $B_{s}\to\mu^{+}\mu^{-}$ reads

$H_{\mathrm{eff}}=N\sum_{i=A,S,P}\left(C_{i}Q_{i}+C_{i}^{\prime}Q_{i}^{\prime}\right)\,.$
(1)

The operators in Eq. (1) are

$\displaystyle Q_{A}=$
$\displaystyle\left(\bar{b}\gamma_{\mu}P_{L}s\right)\left(\bar{\mu}\gamma^{\mu}\gamma_{5}\mu\right)\,,$
$\displaystyle\qquad Q_{A}^{\prime}=\,$
$\displaystyle\left(\bar{b}\gamma_{\mu}P_{R}s\right)\left(\bar{\mu}\gamma^{\mu}\gamma_{5}\mu\right)\,,$
(2) $\displaystyle Q_{S}=$
$\displaystyle\left(\bar{b}P_{L}s\right)\left(\bar{\mu}\mu\right)\,,$
$\displaystyle\qquad Q_{S}^{\prime}=\,$
$\displaystyle\left(\bar{b}P_{R}s\right)\left(\bar{\mu}\mu\right)\,,$
$\displaystyle Q_{P}=$
$\displaystyle\left(\bar{b}P_{L}s\right)\left(\bar{\mu}\gamma_{5}\mu\right)\,,$
$\displaystyle\qquad Q_{P}^{\prime}=\,$
$\displaystyle\left(\bar{b}P_{R}s\right)\left(\bar{\mu}\gamma_{5}\mu\right)\,,$

and are multiplied with their respective _Wilson coefficients_
$C_{A},\ldots,C_{P}^{\prime}$ which contain the dependence on the heavy
masses. The normalisation factor in Eq. (1) is

$\displaystyle N$
$\displaystyle=\;\frac{G_{F}^{2}M_{W}^{2}}{\pi^{2}}V_{ts}V_{tb}^{*}\,=\,\frac{G_{F}\alpha_{\mathrm{em}}\left(\mu\right)}{\sqrt{2}\pi\sin^{2}\theta_{w}}V_{ts}V_{tb}^{*},$
(3)

which complies with the conventions of Ref. Hermann:2013kca . The second “$=$”
sign only holds to lowest order in the electroweak interaction, while in
higher orders the relation between the Fermi constant $G_{F}$ and the
electromagnetic coupling $\alpha_{\mathrm{em}}=e^{2}/(4\pi)$, the weak mixing
angle $\theta_{w}$ and the $W$ mass $M_{W}$ is modified. Electroweak
corrections have been calculated in Ref. Bobeth:2013tba  and e.g. remove the
ambiguities related to the choice of the renormalisation scheme for these
parameters; in the second version for $N$ in Eq. (3) this issue also includes
the choice of the scale in the running $\alpha_{\mathrm{em}}$. We choose the
first definition $N\propto G_{F}^{2}$ in this paper, for which the electroweak
corrections to the SM contribution to $C_{A}$ are as small as $-2.4\%$
Bobeth:2013tba .

We introduce the perturbative expansion of the $C_{i}$ as

$C_{i}=C_{i}^{(0)}+\left(\frac{\alpha_{s}}{4\pi}\right)C_{i}^{(1)}+\ldots\,,$
(4)

where $C_{i}^{(0)}$ denotes the leading order (LO), arising in the SM from
one-loop electroweak diagrams. $C_{i}^{(1)}$ comprises the next-to-leading
order (NLO) QCD corrections. In the SM only $C_{A}$ is relevant,
$C_{A}^{\prime}$ is suppressed w.r.t. $C_{A}$ by the ratio
$m_{b}m_{s}/M_{W}^{2}$ involving strange and bottom masses $m_{\\{s,b\\}}$ and
$C_{S,P}^{(\prime)}$ receive additional suppression factors of
$M_{B_{s}}^{2}/M_{W}^{2}$. The leading contributions $C_{i}^{(1)}$ arise from
one-loop electroweak diagrams at order $G_{F}^{2}$ in the SM, hence the
branching ratio is rather small.

The average time-integrated branching ratio is given by DeBruyn:2012wj ;
DeBruyn:2012wk

$\overline{\mathcal{B}}\left(B_{s}\to\mu^{+}\mu^{-}\right)=|N|^{2}\frac{M_{B_{s}}^{3}f_{B_{s}}^{2}}{8\pi\Gamma_{H}^{s}}\beta\left[\absolutevalue{r\left(C_{A}-C_{A}^{\prime}\right)-u\left(C_{P}-C_{P}^{\prime}\right)}^{2}F_{P}+\absolutevalue{u\beta\left(C_{S}-C_{S}^{\prime}\right)}^{2}F_{S}\right]\,,$
(5)

with dimensionless quantities

$r=\frac{2m_{\mu}}{M_{B_{s}}}\,,\qquad\beta=\sqrt{1-r^{2}}\,,\qquad
u=\frac{M_{B_{s}}}{m_{b}+m_{s}}\,.$ (6)

Here, $\Gamma_{H}^{s}$ ($\Gamma_{L}^{s}$) denotes the decay width of the
heavier (lighter) $B_{s}$ mass eigenstate, and the factors $F_{P}$ and $F_{S}$
account for the mixing of the $B_{s}$-$\bar{B}_{s}$ system, given by

$\displaystyle F_{P}=\,$ $\displaystyle
1-\frac{\Gamma_{L}^{s}-\Gamma_{H}^{s}}{\Gamma_{L}^{s}}\sin^{2}\left[\arg\left(r\left(C_{A}^{\prime}-C_{A}\right)-u\left(C_{P}^{\prime}-C_{P}\right)\right)\right]\,,$
$\displaystyle F_{S}=\,$ $\displaystyle
1-\frac{\Gamma_{L}^{s}-\Gamma_{H}^{s}}{\Gamma_{L}^{s}}\cos^{2}\left[\arg\left(C_{S}^{\prime}-C_{S}\right)\right]$
(7)

in the absence of significant CP-violating New Physics contributions to the
$B_{s}-\bar{B}_{s}$ mixing amplitude. In writing Eq. (5), we have used the
pseudoscalar decay constant $f_{B_{s}}$ to rewrite the operator matrix
elements as

$\displaystyle\matrixelement{0}{\bar{b}\gamma_{\mu}\gamma_{5}s}{B_{s}(p)}=\,$
$\displaystyle\mathrm{i}p_{\mu}f_{B_{s}}\,,$
$\displaystyle\matrixelement{0}{\bar{b}\gamma_{5}s}{B_{s}(p)}=\,$
$\displaystyle-\mathrm{i}f_{B_{s}}\frac{M_{B_{s}}^{2}}{m_{b}+m_{s}}\,,$ (8)

where the second equation follows from the first one by use of the Dirac
equation.

##  3 2HDM with suppressed down-type FCNC couplings

###  3.1 General Yukawa sector

The Yukawa Lagrangian of the general (“type III”) 2HDM reads

$\displaystyle L_{Y}$
$\displaystyle=-\overline{Q}_{f}^{\prime}\Big{[}\overline{Y}_{fi}^{d}H_{d}+\overline{\epsilon}_{fi}^{d}H_{u}\Big{]}d_{iR}^{\prime}\;-\;\overline{Q}_{f}^{\prime}\Big{[}\overline{Y}_{fi}^{u}\epsilon
H_{u}^{*}+\overline{\epsilon}_{fi}^{u}\epsilon
H_{d}^{*}\Big{]}u_{iR}^{\prime}+h.c.$
$\displaystyle\equiv-\overline{Q}^{\prime}\Big{[}\overline{Y}^{d}H_{d}+\overline{\epsilon}^{d}H_{u}\Big{]}d_{R}^{\prime}\;-\;\overline{Q}^{\prime}\Big{[}\overline{Y}^{u}\epsilon
H_{u}^{*}+\overline{\epsilon}^{u}\epsilon H_{d}^{*}\Big{]}u_{R}^{\prime}+h.c.$
(9)

with

$\displaystyle H_{u,d}$ $\displaystyle=\begin{pmatrix}H_{u,d}^{+}\\\
H_{u,d}^{0}\end{pmatrix}$ $\displaystyle\epsilon H_{u,d}^{*}$
$\displaystyle=\begin{pmatrix}H_{u,d}^{0*}\\\ -H_{u,d}^{-}\end{pmatrix}$ (10)
$\displaystyle\mbox{and~{}~{}~{}}Q_{f}^{\prime}$
$\displaystyle=\begin{pmatrix}u_{fL}^{\prime}\\\
d_{fL}^{\prime}\end{pmatrix}.$ (11)

The subscripts $f,i=1,2,3$ label the generations, e.g.
$u_{3L}^{\prime}=t_{L}^{\prime}$. The notation in Eq. (9) follows Ref.
Crivellin:2013wna , except that our $H_{d}$ corresponds to $-\epsilon
H_{d}^{*}$ of that paper.

The vacuum expectation values (vevs) and related quantities are

$\displaystyle\langle H_{u}^{0}\rangle$ $\displaystyle=v_{u}=v\sin\beta,$
$\displaystyle\langle H_{d}^{0}\rangle$ $\displaystyle=v_{d}=v\cos\beta,$
$\displaystyle\tan\beta$ $\displaystyle:=\frac{v_{u}}{v_{d}},$ $\displaystyle
v$ $\displaystyle:=\sqrt{v_{u}^{2}+v_{d}^{2}}=174\,\mbox{GeV}.$ (12)

The quark mass matrices are

$\displaystyle M^{d}$
$\displaystyle=\overline{Y}^{d}v\cos\beta+\overline{\epsilon}^{d}v\sin\beta,$
$\displaystyle M^{u}$
$\displaystyle=\overline{Y}^{u}v\sin\beta+\overline{\epsilon}^{u}v\cos\beta,$
(13)

which we diagonalise in the usual way:

$\displaystyle u_{L,R}^{\prime}=S_{L,R}^{u}u_{L,R},$ $\displaystyle
d_{L,R}^{\prime}=S_{L,R}^{d}d_{L,R},$ (14) $\displaystyle
S_{L}^{d\dagger}M^{d}S_{R}^{d}=\hat{M}^{d}=$
$\displaystyle\begin{pmatrix}m_{d}&0&0\\\ 0&m_{s}&0\\\
0&0&m_{b}\end{pmatrix},$ $\displaystyle
S_{L}^{u\dagger}M^{u}S_{R}^{u}=\hat{M}^{u}=$
$\displaystyle\begin{pmatrix}m_{u}&0&0\\\ 0&m_{c}&0\\\
0&0&m_{t}\end{pmatrix},$ (15)

with the unprimed fields corresponding to quark mass eigenstates. The gauge
sector of the 2HDM is invariant under independent unitary rotations of the
fields $Q^{\prime}$, $d_{R}^{\prime}$, and $u_{R}^{\prime}$ in flavour space.
We use Eq. (14) and choose

$\displaystyle Q^{\prime}$ $\displaystyle\equiv S_{L}^{d}Q$ (16)

and find $L_{Y}$ in the so-called _down basis_ :

$\displaystyle L_{Y}$
$\displaystyle\equiv-\overline{Q}\Big{[}Y^{d}H_{d}+\epsilon^{d}H_{u}\Big{]}d_{R}\;-\;\overline{Q}\,\Big{[}Y^{u}\epsilon
H_{u}^{*}+\epsilon^{u}\epsilon H_{d}^{*}\Big{]}u_{R}+h.c.$ (17)

with the appropriately transformed Yukawa matrices

$\displaystyle Y^{u,d}$
$\displaystyle=S_{L}^{d\dagger}\overline{Y}^{u,d}S_{R}^{u,d},$
$\displaystyle\epsilon^{u,d}$
$\displaystyle=S_{L}^{d\dagger}\overline{\epsilon}^{u,d}S_{R}^{u,d}$ (18)

and the CKM matrix

$\displaystyle V$ $\displaystyle=S_{L}^{u\dagger}S_{L}^{d}.$ (19)

The Yukawa Lagrangian $L_{Y}$ in Eq. (17) is manifestly SU(2) invariant, with
the SU(2) doublet

$\displaystyle Q$
$\displaystyle=S_{L}^{d\dagger}Q^{\prime}=\begin{pmatrix}V^{\dagger}\,u_{L}\\\
d_{L}\end{pmatrix}.$

Eq. (17) is our starting point; the Yukawa matrices are related to the
diagonal mass matrices as

$\displaystyle\frac{\hat{M}_{d}}{v}$
$\displaystyle=\,Y^{d}\cos\beta+\epsilon^{d}\sin\beta,$
$\displaystyle\frac{\hat{M}_{u}}{v}$
$\displaystyle=\,V\left(Y^{u}\sin\beta+\epsilon^{u}\cos\beta\right).$ (20)

Non-vanishing off-diagonal entries of $Y^{d,u}$, $\epsilon^{d,u}$ give rise to
FCNC couplings of the neutral components of the Higgs doublets.

In a general 2HDM the quantity $\tan\beta$ has no physical meaning: One can
arbitrarily rotate $\begin{pmatrix}H_{u},H_{d}\end{pmatrix}$ in Eq. (9)
leading to a Lagrangian of the same form (yet with different Yukawa matrices)
and the rotation angle will add to $\beta$. The situation is different in
variants of the 2HDM in which $H^{u}$ and $H^{d}$ are distinguished by quantum
numbers which forbid such rotations of
$\begin{pmatrix}H_{u},H_{d}\end{pmatrix}$. Prominent examples are the 2HDM of
type I and II, in which two out of the four Yukawa matrices in Eqs. (9) and
(17) are absent. The type-I model corresponds to $\epsilon^{d}=Y^{u}=0$. The
type-II model is instead found for $\epsilon^{d}=\epsilon^{u}=0$.

The doublets $\phi_{\rm SM}$, $\phi_{\rm new}$ of the _Higgs basis_
Haber:2006ue ; Eriksson:2009ws ; Branco:1999fs  are defined by a rotation of
the two doublets $H_{u}$ and $H_{d}$ by the angle $\beta$ such that $\phi_{\rm
new}$ has no vev:

$\displaystyle\begin{pmatrix}\phi_{\rm new}\\\ \phi_{\rm SM}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\cos\beta&-\sin\beta\\\
\sin\beta&\cos\beta\end{pmatrix}\,\begin{pmatrix}H_{u}\\\ H_{d}\end{pmatrix}$
(21)

One has

$\displaystyle\phi_{\rm SM}$ $\displaystyle=\begin{pmatrix}\displaystyle
G^{+}\\\\[2.84526pt] \displaystyle
v+\frac{\phi^{0}+iG^{0}}{\sqrt{2}}\end{pmatrix},$ $\displaystyle\phi_{\rm
new}$ $\displaystyle=\begin{pmatrix}\displaystyle H^{+}\\\\[2.84526pt]
\displaystyle\frac{\phi^{0\prime}+iA^{0}}{\sqrt{2}}\end{pmatrix}\,.$ (22)

Next we use Eq. (21) to express $L_{Y}$ in terms of $\phi_{\rm SM}$ and
$\phi_{\rm new}$:

$\displaystyle L_{Y}=$
$\displaystyle\quad\overline{Q}\Big{[}-Y^{d}\cos\beta-\epsilon^{d}\sin\beta\Big{]}\phi_{\rm
SM}d_{R}\quad\,+\,\overline{Q}\Big{[}Y^{d}\sin\beta-\epsilon^{d}\cos\beta\Big{]}\phi_{\rm
new}d_{R}$
$\displaystyle+\overline{Q}V^{\dagger}\Big{[}-Y^{u}\sin\beta-\epsilon^{u}\cos\beta\Big{]}\epsilon\phi_{\rm
SM}^{*}u_{R}\,+\,\overline{Q}V^{\dagger}\Big{[}-Y^{u}\cos\beta+\epsilon^{u}\sin\beta\Big{]}\epsilon\phi_{\rm
new}^{*}u_{R}+h.c.$ (23)

By using Eq. (13) to eliminate $Y^{u}$ and $Y^{d}$ one can write the couplings
to physical charged and neutral Higgs bosons in Eq. (3.1) as

$\displaystyle\mathcal{L}_{Y}^{\rm phys}=\,$
$\displaystyle-\bar{u}_{L}\left[\frac{\hat{M}^{u}}{v}\cot\beta+g^{u}\right]u_{R}\frac{\phi^{0\prime}-\mathrm{i}A^{0}}{\sqrt{2}}+\bar{d}_{L}\left[\frac{\hat{M}^{d}}{v}\tan\beta+g^{d}\right]d_{R}\frac{\phi^{0\prime}+\mathrm{i}A^{0}}{\sqrt{2}}$
$\displaystyle+\bar{u}_{L}V\left[\frac{\hat{M}^{d}}{v}\tan\beta+g^{d}\right]d_{R}H^{+}+\bar{d}_{L}V^{\dagger}\left[\frac{\hat{M}^{u}}{v}\cot\beta+g^{u}\right]u_{R}H^{-}$
$\displaystyle-\overline{d}_{L}\frac{\hat{M}^{d}}{v}d_{R}\,\left(v+\frac{\phi^{0}}{\sqrt{2}}\right)\;-\;\overline{u}_{L}\frac{\hat{M}^{u}}{v}u_{R}\,\left(v+\frac{\phi^{0}}{\sqrt{2}}\right)\;+\;h.c.$
(24)

with the matrices Crivellin:2013wna

$\displaystyle g^{d}$
$\displaystyle=-\epsilon^{d}\sin\beta\left(\tan\beta+\cot\beta\right),$
$\displaystyle g^{u}$
$\displaystyle=-\epsilon^{u}\cos\beta\left(\tan\beta+\cot\beta\right).$ (25)

The non-diagonal matrices $g^{d}$ and $g^{u}$ characterise the deviations from
the popular type-II 2HDM (for which $g^{d}=g^{u}=0$) and can induce flavour-
changing couplings of neutral Higgs bosons. Note that the type-I model is also
included in the formalism and recovered by using Eq. (13) with $Y^{u}=0$ in
the expression for $g^{d}$. For our loop calculation and the phenomenological
analysis it is advantageous to work with $g^{d,u}$ rather than
$\epsilon^{d,u}$, especially for the definition of the renormalisation
prescriptions. We write $g^{d}_{ij}=g_{d_{i}d_{j}}$, where $d_{i}$ is the
$i$-th down-type quark flavor, $i=1,2,3$, and similarly for $g^{u}$.

We are interested in 2HDM scenarios in which current and forthcoming
measurements by LHCb and CMS can find deviations in
$B(B_{s}\to\mu^{+}\mu^{-})$ from the SM prediction. This consideration guides
us to lepton Yukawa couplings of the type-II form,

$\mathcal{L}_{Y}^{l}=+\bar{l}_{L}\left[\frac{\hat{M}^{l}}{v}\tan\beta\right]l_{R}\frac{\phi^{0\prime}+\mathrm{i}A^{0}}{\sqrt{2}}+\bar{\nu}_{L}\left[\frac{\hat{M}^{l}}{v}\tan\beta\right]l_{R}H^{+}+h.c.$
(26)

and large values of $\tan\beta$ to enhance the muon Yukawa coupling
$Y_{\mu}=m_{\mu}/v\tan\beta$ over its SM value. $\hat{M}^{l}$ in Eq. (26)
denotes the diagonal mass matrix of the charged leptons.

We restrict ourselves to the CP-conserving Higgs potential, in which diagrams
involving $A^{0}$ contribute to the Wilson coefficient $C_{P}$, while diagrams
with the CP-even neutral Higgs states contribute to $C_{S}$. The two Higgs
bosons $\phi^{0}$ and $\phi^{0\prime}$ are in general not mass eigenstates.
The latter are given by

$\begin{pmatrix}h^{0}\\\ H^{0}\\\
A^{0}\end{pmatrix}=\begin{pmatrix}\sin\left(\beta-\alpha\right)&\cos\left(\beta-\alpha\right)&0\\\
\cos\left(\beta-\alpha\right)&-\sin\left(\beta-\alpha\right)&0\\\
0&0&1\end{pmatrix}\,\begin{pmatrix}\phi^{0}\\\ \phi^{0\prime}\\\
A^{0}\end{pmatrix}\,.$ (27)

The angle $\alpha$ is an a priori arbitrary parameter of the theory that
serves to diagonalise the mass matrix of neutral Higgs bosons.

###  3.2 Spurion expansion

The 2HDM of type I and II (and their variants with modified lepton couplings)
invoke (softly broken) $\mathbb{Z}_{2}$ symmetries to forbid FCNC couplings of
the neutral Higgs bosons. This was motivated by the wish to find non-standard
Higgs bosons at modern colliders, because for generic values of the Yukawa
matrices in Eq. (17) constraints from FCNC processes like meson-antimeson
mixing push these masses to values outside the reach of of LEP, Tevatron, and
LHC. With the absence of discoveries of non-standard Higgs bosons this line of
arguments loses its appeal and the consideration of more general Yukawa
sectors is in order.

The type-II model is the most studied variant of the 2HDM for two reasons:
Firstly, it constitutes the tree-level Higgs sector of the Minimal
Supersymmetric Standard Model (MSSM), in which the holomorphy of the
superpotential enforces $\epsilon^{u,d}=0$. Secondly, the type-II model is
phenomenologically especially interesting, because in this model FCNC
processes are sensitive to loop effects of the charged Higgs boson. A
prominent example of the latter feature is the branching ratio
$\mathcal{B}(B\to s\gamma)$, which sets a stringent bound on the charged-Higgs
mass Misiak:2017bgg . The type-II 2HDM further permits the possibility of
large down-type Yukawa couplings, a scenario motivated by the possibility of
bottom-top Yukawa unification. Such large-$\tan\beta$ scenarios are
efficiently constrained by $\mathcal{B}(B_{s}\to\mu^{+}\mu^{-})$ Logan:2000iv
; Bobeth:2013uxa ; Hermann:2013kca  and we will come back to this topic in
section 5. Concerning the first above-mentioned motivation, the Higgs sector
of the MSSM is really richer than that of the type-II 2HDM: In the limit of
infinitely heavy superpartners one encounters non-decoupling loop-induced
Yukawa matrices $\epsilon^{u,d}$, an effect caused by the supersymmetry-
breaking terms Banks:1987iu ; Hall:1993gn ; Carena:1994bv . Despite the loop
suppression large effects are possible in FCNC processes with down-type quarks
which involve the product $\epsilon^{d}\tan\beta$ Babu:1999hn ; Isidori:2001fv
; Buras:2002vd ; Buras:2002wq ; Hofer:2009xb ; Gorbahn:2009pp  with huge
impact on $\mathcal{B}(B_{s}\to\mu^{+}\mu^{-})$ Babu:1999hn ; Dedes:2001fv ;
Isidori:2001fv ; Buras:2002wq ; Buras:2002vd .

The phenomenological constraints from meson-antimeson mixing and rare
(semi-)leptonic decays place severe bounds on the off-diagonal elements of
$\epsilon^{d}$, while those of $\epsilon^{u}$ are essentially unconstrained
except for $\epsilon^{u}_{12}$ and $\epsilon^{u}_{21}$. This situation calls
for a variant of the 2HDM in which $Y^{u,d}$ and $\epsilon^{u}$ are arbitrary,
while $\epsilon^{d}$ is suppressed. A strong motivation for such a model is
the possibility of spontaneous CP violation, implemented through a Higgs
potential developing complex vevs and real Yukawa matrices. Spontaneous CP
violation is not possible with a 2HDM of type I or II, but requires at least
three out of the four matrices in Eq. (17) to be non-zero. Avoiding fine-
tuning implies that the better part of the needed effect stems from
$\epsilon^{u}$, while $\epsilon^{d}$ can be neglected Nierste:2019fbx .
However, in such a model, the mixing of the neutral Higgs fields is different
from Eq. (27) and instead involves all three fields. Yet for the CP-conserving
observables considered in this paper this feature is of minor relevance. The
2HDM scenario with sizable $Y^{u,d}$ and $\epsilon^{u}$ has the appealing
feature that it simultaneously permits both measurable effects in FCNC
processes _and_ sufficiently light masses of the non-standard Higgs particles
enabling their discovery at the LHC.

Setting $\epsilon^{d}$ naively to zero leads to a non-renormalisable model,
because there are UV-divergent loops involving up-type quarks with $Y^{u,d}$
and $\epsilon^{u}$ couplings, which require counterterms proportional to
elements of $\epsilon^{d}$. Whenever one seeks to constrain the elements
$g^{u}_{jk}$ of Eq. (3.1) from FCNC transitions of down-type quarks, one must
foresee such a counterterm to find a meaningful prediction. For example,
$B_{s}\to\mu^{+}\mu^{-}$ is a $b\to s$ transition constraining $g_{ct}$ and
the corresponding loop contribution requires a counterterm for $g_{sb}$. The
minimal renormalisable theory is found by invoking flavour symmetries and
systematically expanding in terms of the spurions breaking these symmetries
Chivukula:1987py ; DAmbrosio:2002vsn . The 2HDM gauge sector is invariant
under unitary rotations $Q\to U_{Q}Q$, $d_{R}\to U_{d}d_{R}$, $u_{R}\to
U_{u}u_{R}$ in flavour space with $(U_{Q},U_{d},U_{u})\in SU(3)^{3}$ and the
Yukawa sector is formally invariant under this flavour symmetry if one
transforms the matrices in Eq. (17) as

$\displaystyle Y^{u,d}\to
U_{Q}Y^{u,d}U_{u,d}^{\dagger},\qquad\qquad\epsilon^{u,d}\to
U_{Q}\epsilon^{u,d}U_{u,d}^{\dagger}.$ (28)

We propose to categorise the classes of renormalisable 2HDM in terms of the
spurions breaking the $SU(3)^{3}$ flavour symmetry of the quark sector.111The
generalisation to the lepton sector is straightforward, but not relevant for
the calculations in this paper. The minimal choice are two spurions, with just
two physically distinct possibilities. Both comply with the definition of
minimal flavour violation (MFV) as defined in Ref. DAmbrosio:2002vsn . The
first possibility is to take $Y^{u,d}$ as spurions and express the other two
Yukawa matrices as
$\epsilon^{u,d}=P_{u,d}(Y^{u}Y^{u\dagger},Y^{d}Y^{d\dagger})\,Y^{u,d}$, where
$P_{u}$ and $P_{d}$ are polynomials. This variant is discussed in Ref.
DAmbrosio:2002vsn  and amounts to a generalisation of the 2HDM of type II. It
also constitutes a renormalisable extention of the aligned 2HDM of Ref.
Tuzon:2010vt ; Eberhardt:2020dat , in which $P_{u,d}$ are constants. If the
2HDM is the low-energy limit of a more fundamental theory obeying the
considered two-spurion symmetry-breaking pattern, the latter will naturally
induce $\epsilon^{u,d}$ in the described way as well. An example for such a UV
theory is the MSSM with a flavour-blind supersymmetry breaking mechanism (such
as gauge mediation). The second possibility is to choose $Y^{d}$,
$\epsilon^{u}$ as spurions, leading to a generalisation of the type-I 2HDM.

There are two possibilities for a 2HDM with three spurions, which can be taken
as $Y^{u,d}$ plus either $\epsilon^{u}$ or $\epsilon^{d}$. The first
possibility is the phenomenologically interesting one and studied in this
paper. The expansion up to third order reads

$\displaystyle\epsilon^{d}=$
$\displaystyle\;cY^{d}\,+\,c_{11}Y^{d}Y^{d\dagger}Y^{d}$
$\displaystyle\phantom{\;cY^{d}}\,+\,b_{11}\epsilon^{u}\epsilon^{u\dagger}Y^{d}\,+\,b_{12}\epsilon^{u}Y^{u\dagger}Y^{d}\,+\,b_{21}Y^{u}\epsilon^{u\dagger}Y^{d}\,+\,b_{22}Y^{u}Y^{u\dagger}Y^{d}$
(29)

with complex coefficients $c,\ldots,b_{22}$.

Concerning Eq. (29) several comments are in order:

  * •

The spurion expansion is only meaningful, if the contributions to the off-
diagonal elements of $\epsilon^{d}$ from higher electroweak orders (with five
or more Yukawa matrices) are small, so that they can be neglected. A
sufficient condition for this is realised in scenarios in which
$c_{11},\ldots,b_{22}$ are induced by one-loop contributions in either the UV
completion or the 2HDM, while terms with $(2n+1)$ spurions are only generated
at $n$-loop order and beyond. We consider this scenario throughout this
paper.222A different application of Eq. (29), in which
$c_{11},\ldots,b_{22}={\cal O}(1)$ is allowed, is the case that $Y^{u}$ and
$\epsilon^{u}$ are almost aligned, so that Eq. (20) means small off-diagonal
matrix elements of these matrices in the chosen basis. Since
$Y^{u}_{33},\epsilon^{u}_{33}={\cal O}(1)$ and further $Y^{d}_{33}={\cal
O}(1)$ is possible, some terms with five spurions must be added to Eq. (29),
just as in the MFV case of Ref. DAmbrosio:2002vsn . Additional QCD corrections
(e.g. an extra loop with a gluon) comply with the pattern in Eq. (29), i.e.
QCD renormalises the coefficients, but does not induce new ones.

  * •

By rotating $\begin{pmatrix}H_{u},H_{d}\end{pmatrix}$ in Eq. (17) one can
eliminate $c$ in Eq. (29). But in general radiative corrections bring this
term back and a counterterm to $c$ is needed, unless one corrects the rotation
in each order of perturbation theory. It is therefore advisable to keep $c$ in
Eq. (29); we treat it as a perturbative quantity with $c=0$ at tree level.

  * •

The decay $B_{s}\to\mu^{+}\mu^{-}$ involves the FCNC vertex
$\overline{Q}_{j}$-$d_{Rk}$-$H_{d}$ with $(j,k)=(2,3)$. The dominant one-loop
vertex diagram involves an internal $H_{d}$ line and the product
$Y^{d}Y^{d\dagger}Y^{d}$ or $\epsilon^{u}\epsilon^{u\dagger}Y^{d}$ stemming
from the three $H_{d}$ Yukawa couplings. The UV divergences can be cancelled
by counterterms to $c_{11}$ and $b_{11}$.

  * •

With Eq. (20) we can trade $Y^{u,d}$ in Eq. (29) for the quark masses and CKM
elements. Compared to the SM we find 14 additional complex parameters, the 9
entries of $\epsilon^{u}$ and $c_{11},\ldots b_{22}$. Yet it is much more
convenient to express observables in terms of $g_{jk}^{u}$ of Eq. (3.1)
instead of $\epsilon^{u}_{jk}$ and then use Eqs. (25) and (29) (with Eq. (20))
to calculate the $g_{jk}^{d}$ in terms of the coefficients of the spurion
expansion. While this procedure is needed in a global analysis of all
available data —which is beyond the scope of this paper—, the study of $b\to
s$ transitions alone will simply involve $g_{ct}$ and $g_{tt}$.

In a practical calculation it is cumbersome to implement the renormalisation
procedure in the described way, by providing counterterms to
$\epsilon^{u}_{jk}$ and $c,\ldots,b_{22}$. Instead, it is much easier to
renormalise the $g_{jk}^{u,d}$. If one renormalises all $g_{jk}^{u,d}$ in the
$\overline{\rm MS}$ scheme, one automatically complies with SU(2)$\times$U(1)
gauge symmetry. Therefore it is sufficient to choose the $g_{jk}^{d}$ in such
a way that Eq. (29) is obeyed at tree level. In our calculation we will only
need a counterterm to $g_{sb}$ (in addition to the usual QCD counterterms for
the SM parameters), if $m_{s}$ is set to zero. If $m_{s}$ is kept, one finds
primed Wilson coefficients which require a counterterm to $g_{bs}$.

##  4 Computational setup

In this section, we describe the chain of programs used to generate and
evaluate the corresponding Feynman diagrams at leading and next-to-leading
order. We use FeynRules with the Universal FeynRules Output (UFO)
Christensen:2009jx ; Alloul:2013bka ; Degrande:2011ua  to obtain Feynman rules
for the model. The output is processed by tapir Gerlach:2022qnc  into a
Lagrangian file, which is used with qgraf Nogueira:1991ex  in order to
generate all Feynman diagrams.

We compute the one-loop diagrams for general electroweak gauge parameters for
the $W$ and $Z$ bosons and check that they drop out in the final result. This
is a welcome check for the conversion from UFO to qgraf and tapir. If diagrams
with all different quark flavours are explicitly calculated, we have at one-
loop level $<\mathcal{O}(100)$ Feynman diagrams in the SM,333From the
technical perspective, it is convenient to split the electroweak gauge bosons
into two different “particles”, with different propagator denominators
$\left(k^{2}-M^{2}\right)^{-1}$ and $\left(k^{2}-\xi M^{2}\right)^{-1}$,
respectively, and treat them as different diagrams; hence the large number of
diagrams in the 1-loop calculation. and an additional $\mathcal{O}(300)$
diagrams from contributions with at least one non-SM Higgs boson, most of
which are Higgs-penguin diagrams.

At two loop-order, we perform the calculation in the Feynman gauge for the $W$
and $Z$ bosons. The diagrams are then converted with the help of tapir into
FORM Ruijl:2017dtg  code using Feynman rule definitions that were also
produced in the conversion from UFO to the qgraf Lagrangian file. The
individual expressions for the diagrams are mapped onto integral families
using exp Harlander:1998cmq ; Seidensticker:1999bb  and a custom FORM setup is
used to perform the remaining computational steps. Since the Wilson
coefficients are independent of the momenta of the external particles, we set
the latter to zero, so that only vacuum integrals need to be computed. (An
exception are diagrams with an FCNC self-energy in an external leg, which are
calculated with on-shell $b$ quark and $m_{b}\neq 0$ before the subsequent
limit $m_{b}\to 0$ is taken.) We use a FORM implementation of the algorithm
presented in Ref. Davydychev:1992mt , see Ref. Salomon:2012iha .

##  5 The decay $\mathbf{B_{s}\to\mu^{+}\mu^{-}}$ in the Two-Higgs-Doublet
Model of type-II

Figure 1: Sample diagrams contributing to $C_{A}$ at leading and next-to-
leading order in the SM.

In the SM next-to-leading order corrections are available from Refs.
10.1143/PTP.65.297 ; BUCHALLA1993285 ; BUCHALLA1993225 ; Buchalla:1998ba ;
Misiak:1999yg  and higher order QCD and electroweak corrections have been
computed in Refs. Hermann:2013kca ; Bobeth:2013tba ; Bobeth:2013uxa ;
Beneke:2017vpq ; Beneke:2019slt . Next-to-leading corrections in the type-II
2HDM have been computed in Refs. Logan:2000iv ; Bobeth:2001sq ; Bobeth:2001jm
; Hermann:2014lga , which we have reproduced for the present paper, also as a
check of the automated setup. The result can be expressed in terms of
dimensionless mass ratios of the top quark, $W$ boson and charged Higgs boson
masses,

$\displaystyle x_{t}$
$\displaystyle=\,\frac{m_{t}^{2}}{M_{W}^{2}}\,,\qquad\qquad
r_{H}=\,\frac{m_{t}^{2}}{M_{H^{+}}^{2}}\,.$ (30)

In the SM, only the Wilson coefficient $C_{A}$ receives significant
contributions from diagrams such as the ones depicted in Fig. 1. The Wilson
coefficients $C_{A}^{\prime},C_{S},\ldots,C_{P}^{\prime}$ are suppressed by
powers of ratios of light (external) masses and $M_{W}$. The leading SM
contribution is, moreover, independent of the Yukawa couplings of all external
particles, that is we can set $m_{b}=m_{s}=m_{\mu}=0$ for the contributions
from $W$-box and $Z$-penguin diagrams. At leading and next-to-leading order
they are given by

$\displaystyle C_{A}^{W+Z,(0)}=\,$
$\displaystyle\frac{x_{t}\left(x_{t}-4\right)}{8\left(x_{t}-1\right)}+\frac{3x_{t}^{2}\log
x_{t}}{8\left(x_{t}-1\right)^{2}}\,,$ (31) $\displaystyle C_{A}^{W+Z,(1)}=\,$
$\displaystyle\frac{2x_{t}\left(2x_{t}^{2}+5x_{t}+5\right)}{3\left(x_{t}-1\right)^{2}}-\frac{x_{t}\left(x_{t}^{2}+5x_{t}+2\right)\log
x_{t}}{\left(x_{t}-1\right)^{3}}-\frac{x_{t}\left(x_{t}^{2}+2\right)\mathrm{Li}_{2}\left(1-\frac{1}{x_{t}}\right)}{\left(x_{t}-1\right)^{2}}$
$\displaystyle+\left(\frac{x_{t}\left(x_{t}^{2}+x_{t}+4\right)}{\left(x_{t}-1\right)^{2}}-\frac{6x_{t}^{2}\log
x_{t}}{\left(x_{t}-1\right)^{3}}\right)\log\frac{\mu^{2}}{m_{t}^{2}}\,,$ (32)

where the $Z$-penguin contributions include the flavour-changing quark self-
energy diagrams required to make the penguin diagrams finite. Moreover, since
charm and up quark masses can be neglected compared to $M_{W}$ and $m_{t}$,
after summation over internal up-type quarks all contributions are
proportional to the CKM factor $V_{ts}V_{tb}^{*}$ (contained in the
normalisation factor $N$ defined in Eq. (3)) due to the unitarity of the CKM
matrix.

Figure 2: Sample two-loop Feynman diagrams contributing to the Wilson
coefficients $C_{S}^{(\prime)}$ ($H_{1,2}$) and $C_{P}^{(\prime)}$ ($A^{0}$).
In the left diagram, only flavour-diagonal transitions $i=j$ are possible in
the type-II 2HDM, while the more general 2HDM allows also for transitions with
$i\neq j$, e.g. transitions from a charm quark into a top quark.

In the type-II Two-Higgs-Doublet model, the leading terms are given by
$\mathcal{O}\left(\tan^{2}\beta\right)$ contributions to $C_{S}^{(\prime)}$
and $C_{P}^{(\prime)}$, arising from penguin diagrams with a neutral Higgs
boson, see e.g. Fig. 2, and the $W^{+}$-$H^{+}$ box diagram. The Higgs-penguin
contributions to $C_{S}$ are given by

$\displaystyle C_{S}^{h,(0)}=\,$ $\displaystyle-
r_{b}\left[\frac{r_{H}x_{t}\log
r_{H}}{4\left(r_{H}-1\right)\left(r_{H}-x_{t}\right)}-\frac{r_{H}x_{t}\log
x_{t}}{4\left(x_{t}-1\right)\left(r_{H}-x_{t}\right)}\right]\,,$ (33)
$\displaystyle C_{S}^{h,(1)}=\,$ $\displaystyle-
r_{b}\left[-\frac{8r_{H}x_{t}}{3\left(r_{H}-1\right)\left(x_{t}-1\right)}+\frac{2r_{H}x_{t}\left(3r_{H}-7\right)\log
r_{H}}{3\left(r_{H}-1\right)^{2}\left(r_{H}-x_{t}\right)}-\frac{2r_{H}x_{t}\left(3x_{t}-7\right)\log
x_{t}}{3\left(x_{t}-1\right)^{2}\left(r_{H}-x_{t}\right)}\right.$
$\displaystyle-2\left.\log\left(\frac{\mu^{2}}{m_{t}^{2}}\right)\left(\frac{r_{H}x_{t}}{\left(r_{H}-1\right)\left(x_{t}-1\right)}+\frac{r_{H}x_{t}\log
r_{H}}{\left(r_{H}-1\right)^{2}\left(r_{H}-x_{t}\right)}-\frac{r_{H}x_{t}\log
x_{t}}{\left(x_{t}-1\right)^{2}\left(r_{H}-x_{t}\right)}\right)\right.$
$\displaystyle+\left.\frac{2r_{H}x_{t}\mathrm{Li}_{2}\left(1-\frac{1}{r_{H}}\right)}{r_{H}-x_{t}}-\frac{2r_{H}x_{t}\mathrm{Li}_{2}\left(1-\frac{1}{x_{t}}\right)}{r_{H}-x_{t}}\right]\,,$
(34)

where

$r_{q}=\frac{m_{\mu}m_{q}\,\tan^{2}\beta}{M_{W}^{2}}\,.$ (35)

The Wilson coefficients for the right-handed operators $C_{S}^{\prime h,(i)}$
can be obtained by the replacement $r_{b}\to r_{s}$ in Eqs. (33) and (34).
Further contributions include terms of order
$m_{b}m_{s}m_{\mu}^{2}\tan^{4}\beta/M_{W}^{4}$ and
$m_{t}^{2}m_{\mu}^{2}/M_{W}^{4}$ entering the Wilson coefficients
$C_{A}^{\prime}$ and $C_{A}$, respectively, as well as
$\mathcal{O}\left(m_{b}m_{s}\tan^{2}\beta/M_{W}^{2}\right)$ (in
$C_{A}^{\prime}$) and
$\mathcal{O}\left(m_{t}^{2}\cot^{2}\beta/M_{W}^{2}\right)$ (in $C_{A}$) terms
arising from $Z$-penguin diagrams with a charged Higgs boson. Box diagrams
with a single charged Higgs boson also contribute to $C_{S,P}$
($C_{S,P}^{\prime}$), with Wilson coefficients proportional to $m_{b}$
($m_{s}$). We do not explicitly list these contributions here; they can be
found in Refs. Bobeth:2001sq ; Bobeth:2001jm ; Hermann:2014lga . With the
exception of the doubly muon-mass suppressed $H^{+}$–$H^{-}$ box contributions
to $C_{A}$, we include all of these additional terms in our analysis. The
feature that $C_{S}$ and $C_{P}$ are proportional to $m_{b}$ while their
primed counterparts are proportional to $m_{s}$ holds beyond the type-II 2HDM
in our more general 2HDM with three spurions because of
$\epsilon^{d}=P_{d}\left(\epsilon^{u},Y^{u},Y^{d}\right)Y^{d}$, entailing
factors of $m_{d_{j}}$ in Yukawa couplings of $d_{jR}$. A rather remarkable
feature of the type-II 2HDM is the fact that the leading terms in $\tan\beta$
depend only on $\tan\beta$ and the charged-Higgs boson mass $M_{H^{+}}$, that
is they are independent of the parameters of the neutral Higgs sector
Logan:2000iv . In the type-II 2HDM, the leading $\tan^{2}\beta$ contributions
to the pseudoscalar and scalar Wilson coefficients satisfy the rather simple
relation

$C_{S}=C_{P}\,,\qquad C_{S}^{\prime}=-C_{P}^{\prime}\,.$ (36)

Figure 3: Branching ratio $\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)$ in
the type-II Two-Higgs-Doublet model at next-to-leading order in QCD, for
different values of $\tan\beta$. All running parameters are evaluated at the
scale $\mu=\overline{m}_{t}\left(\overline{m}_{t}\right)$. The black dashed
and dotted lines denote the central value and the experimental $2\sigma$ error
of the LHCb measurement LHCb:2021vsc ; LHCb:2021awg , respectively, the blue
band shows the corresponding information for the CMS result CMS:2022dbz . The
purple dashed and dotted bands indicate the SM prediction with uncertainties,
presented in Ref. Bobeth:2013uxa . Note that the SM predictions in the purple
band were obtained including also $\mathcal{O}\left(\alpha_{s}^{2}\right)$ and
$\mathcal{O}\left(\alpha_{\mathrm{em}}\right)$ contributions.

In Fig. 3 we show the branching ratio
$\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)$ in the type-II 2HDM as a
function of the charged Higgs boson mass. The horizontal gray, blue, and
violet bands correspond to the experimental measurements of LHCb Zyla:2020zbs
and CMS CMS:2022dbz  as well as the theory prediction Bobeth:2013uxa ,
respectively, including 2$\sigma$ uncertainties for the experimental values
(cf. Table. 1). The theory prediction of Ref. Bobeth:2013uxa  uses
$|V_{cb}|=0.0424\pm 0.0009$, which is close to today’s value inferred from
inclusive $b\to c\ell\nu$ decays. If one uses $|V_{cb}|=0.03936\pm 0.00068$
from exclusive $B$ decays HFLAV:2019otj  the central value of the theory
prediction for $10^{9}\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)$ drops
from $3.65$ to $3.15$.

The coloured lines are predictions from the 2HDM for different values of
$\tan\beta$. It is interesting to note that for $\tan\beta=25$ low values of
$M_{H^{\pm}}$ are required to reproduce the central value of the LHCb
measurement. Recently it has been pointed out that LHC data still permit
$M_{H^{\pm}}\leq 400\,$GeV with couplings compatible with solutions of the
$b\to c\tau\nu$ flavour anomalies Blanke:2022kpi . Charged-Higgs explanations
of the latter have been found viable in Ref. Blanke:2018yud ; Blanke:2019qrx
and are invigorated by recent LHCb data on $B\to D^{(*)}\tau\nu$ LHCbSem , see
Refs. Iguro:2022yzr ; Fedele:2022iib . Of course, in the limit
$M_{H^{\pm}}\to\infty$ all 2HDM curves approach the SM prediction. Note that
in the type-II model there is a $\tan\beta$-independent 95% C.L. lower bound
on $M_{H^{\pm}}$ in the range 570-800 GeV from $B\to X_{s}\gamma$
Misiak:2017bgg . This bound can be easily weakened with our model’s additional
couplings discussed in the following section. While charged Higgs searches at
the LHC are compatible with the quoted $H^{\pm}$ masses ATLAS:2018gfm , the
lower bound on $M_{H^{\pm}}$ inferred from the data on $gg/b\bar{b}\to
A^{0}[\to\tau^{+}\tau^{-}]$ searches is larger than 1 TeV ATLAS:2022yvl  in
perturbed versions of the type-II 2HDM, but can be weakened if e.g. the
$A^{0}$ coupling to $\tau$’s is modified w.r.t. to the type-II form in Eq.
(26).

##  6 Additional contributions in a model with flavour-changing neutral Yukawa
couplings

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

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

Figure 5: The loop-induced change of flavours $s\to b$ via a quark self-energy
diagram. The helicity flip denoted by the cross ensures a factor of
$\tan\beta$ and linearity in $\hat{g}^{u}$.

###  6.1 Pseudoscalar Wilson coefficient $C_{P}$

At tree level the Wilson coefficient $C_{P}$ originating from the diagram in
Fig. 4 with a virtual $A^{0}$ boson is given by

$C_{P}^{(0),\mathrm{tree}}=-\frac{\pi^{2}}{G_{F}^{2}M_{W}^{2}V_{tb}^{*}V_{ts}}\frac{m_{\mu}\,\tan\beta}{2vM_{A^{0}}^{2}}g_{sb}^{*}=-\frac{g_{sb}^{*}}{g_{bs}}C_{P}^{\prime(0)}\,.$
(37)

Since there is no loop-suppression, even small values of $g_{bs}$ and
$g_{sb}^{*}$ can have significant impact on the branching ratio, but the
spurion expansion in Eq. (29) naturally leads to (parametrically suppressed)
small couplings. For convenience, let us define the linear combinations
$g_{bs}^{\pm}\equiv g_{bs}\pm g_{sb}^{*}$.

At one-loop order only the self-energy diagrams contain an enhancement factor
$\tan^{3}\beta$ whereas the vertex contributions contain at most a quadratic
term. The one-loop self-energy contributions to $C_{P}$ and $C_{P}^{\prime}$
are ultra-violet divergent. The corresponding counterterm is generated from
Eq. (37) and is of the form

$\displaystyle g_{bs}^{+,0}$ $\displaystyle=$ $\displaystyle
g_{bs}^{+}+\delta_{A^{0},bs}^{(0)}+\left(\frac{\alpha_{s}}{4\pi}\right)\delta_{A^{0},bs}^{(1)}+\ldots\,.$
(38)

In the $\overline{\rm MS}$ renormalisation scheme the one-loop contribution is
given by

$\delta_{A^{0},bs}^{(0)}=-\frac{1}{\epsilon}\frac{\sqrt{2}G_{F}m_{t}\tan^{2}\beta}{8\pi^{2}}\left[m_{s}g_{ct}V_{ts}V_{cb}^{*}+\left(m_{s}g_{tt}V_{ts}+m_{b}g_{ct}^{*}V_{cs}+m_{b}g_{tt}^{*}V_{ts}\right)V_{tb}^{*}\right]\,.$
(39)

Here, the terms with $m_{b}$ renormalise the Wilson coefficient $C_{P}$, while
the ones with $m_{s}$ renormalise the Wilson coefficient $C_{P}^{\prime}$. The
renormalised finite pseudoscalar Wilson coefficients read

$\displaystyle C_{P}^{(0)}=$ $\displaystyle
C_{P}^{(0),\mathrm{tree}}+\,\tilde{N}\,\frac{m_{b}}{M_{A^{0}}^{2}}\,\left(g_{ct}^{*}\frac{V_{cs}}{V_{ts}}+g_{tt}^{*}\right)\,\tilde{C}_{P}^{(0)}\,,$
(40) $\displaystyle C_{P}^{\prime(0)}=$ $\displaystyle
C_{P}^{\prime(0),\mathrm{tree}}\,-\tilde{N}\,\frac{m_{s}}{M_{A^{0}}^{2}}\left(g_{ct}\frac{V_{cb}^{*}}{V_{tb}^{*}}+g_{tt}\right)\,\tilde{C}_{P}^{(0)}\,,$

with normalisation factor

$\tilde{N}=\frac{m_{t}m_{\mu}\tan^{3}\beta}{G_{F}\sqrt{2}M_{W}^{2}v}\,=\,\frac{gm_{t}m_{\mu}\tan^{3}\beta}{2G_{F}M_{W}^{3}},$
(41)

the weak coupling constant $g$, and

$\tilde{C}_{P}^{(0)}=-\frac{1}{8}\left[1+\log\left(\frac{\mu^{2}}{m_{t}^{2}}\right)-\frac{\log
r_{H}}{r_{H}-1}\right]\,.$ (42)

At NLO in QCD (i.e. two-loop order), the counterterm required to obtain a
finite result is given by

$\displaystyle\delta_{A^{0},sb}^{(1)}=$
$\displaystyle-\frac{G_{F}m_{t}\tan^{2}\beta}{\sqrt{2}\pi^{2}}\left[m_{s}g_{ct}V_{ts}V_{cb}^{*}+\left(m_{s}g_{tt}V_{ts}+m_{b}g_{ct}^{*}V_{cs}+m_{b}g_{tt}^{*}V_{ts}\right)V_{tb}^{*}\right]$
$\displaystyle\times\left[\frac{1}{\epsilon^{2}}+\frac{4}{3\epsilon}\right]\,.$
(43)

Note that in addition to the top quark mass $m_{t}$ also the flavour-changing
Yukawa couplings are minimally renormalised. The two-loop contribution to the
Wilson coefficients reads

$\displaystyle C_{P}^{(1)}=$
$\displaystyle\tilde{N}\,\frac{m_{b}}{M_{A^{0}}^{2}}\,\left(g_{ct}^{*}\frac{V_{cs}}{V_{ts}}+g_{tt}^{*}\right)\,\tilde{C}_{P}^{(1)}\,,$
(44) $\displaystyle C_{P}^{\prime(1)}=$
$\displaystyle-\tilde{N}\,\frac{m_{s}}{M_{A^{0}}^{2}}\left(g_{ct}\frac{V_{cb}^{*}}{V_{tb}^{*}}+g_{tt}\right)\,\tilde{C}_{P}^{(1)}\,,$

with

$\displaystyle\tilde{C}_{P}^{(1)}=$
$\displaystyle-\frac{4(r_{H}-2)}{3(r_{H}-1)}+\frac{(3r_{H}-7)\log(r_{H})}{3(r_{H}-1)^{2}}+\left(\frac{7-4r_{H}}{3(r_{H}-1)}-\frac{\log(r_{H})}{(r_{H}-1)^{2}}\right)\log\left(\frac{\mu^{2}}{m_{t}^{2}}\right)$
$\displaystyle-\frac{1}{2}\log^{2}\left(\frac{\mu^{2}}{m_{t}^{2}}\right)+\text{Li}_{2}\left(1-\frac{1}{r_{H}}\right)\,.$
(45)

###  6.2 Scalar Wilson coefficient $C_{S}$

The scalar Wilson coefficients receives contributions from both neutral CP-
even Higgs mass eigenstates $h^{0}$ and $H^{0}$. At tree level, the diagrams
with $h^{0}$ and $H^{0}$ give rise to the Wilson coefficients

$C_{S}^{(0),\mathrm{tree}}=-\frac{\pi^{2}}{G_{F}^{2}M_{W}^{2}V_{tb}^{*}V_{ts}}\frac{m_{\mu}\,\tan\beta}{2v}R_{M}g_{sb}^{*}=\frac{g_{sb}^{*}}{g_{bs}}C_{S}^{\prime(0)}\,,$
(46)

where

$R_{M}=\frac{M_{h^{0}}^{2}\sin^{2}\left(\beta-\alpha\right)+M_{H^{0}}^{2}\cos^{2}\left(\beta-\alpha\right)}{M_{h^{0}}^{2}M_{H^{0}}^{2}}$
(47)

contains the dependence on the neutral Higgs-boson masses. The counterterms
required to cure the divergences at one and two loops can be obtained from
Eqs. (39) and (6.1) by the simple replacements $g_{ct}\to-g_{ct}$ and
$g_{tt}\to-g_{tt}$. The renormalised Wilson coefficients $C_{S}$ and
$C_{S}^{\prime}$ are related to the pseudoscalar ones by

$C_{S}^{(i)}=\frac{R_{M}}{M_{A^{0}}^{2}}C_{P}^{(i)}\,,\qquad
C_{S}^{\prime(i)}=-\frac{R_{M}}{M_{A^{0}}^{2}}C_{P}^{\prime(i)}\,.$ (48)

Note that there are no QCD corrections to $C_{S,P}^{(\prime),\mathrm{tree}}$
since they cancel in the matching calculation.

##  7 Phenomenology

In this section we discuss the possible size of
$\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)$ in our 2HDM under the
constraint that other $b\to s$ processes comply with the data. We include the
tree-level contributions from diagrams with flavour-changing down-type
couplings (cf. Fig. 4), as well as the SM results and the leading quadratic
$\tan\beta$ contributions in the type-II 2HDM to which the additional diagrams
of order $\tan^{3}\beta$ discussed in the previous section add as corrections.

In a generic 2HDM the tree-level couplings $g_{bs}^{\pm}$ will drastically
increase the branching ratio for $B_{s}\to\mu^{+}\mu^{-}$ due to the missing
loop suppression. In our model the spurion expansion suppresses $g_{bs}^{\pm}$
in a controlled way, but still permits large enough contributions to get
phenomenologically interesting effects in $B_{s}\to\mu^{+}\mu^{-}$: Even
rather small up-type couplings $g_{ct}$ significantly modify the Wilson
coefficients $C_{P}^{(\prime)}$ and $C_{S}^{(\prime)}$, as they feature a CKM
factor $V_{cs}$ instead of $V_{ts}$. In the following, we will restrict
ourselves to the experimentally favoured scenario ATLAS:2019nkf ; CMS:2018uag
of aligned Higgs doublets, in which $\sin\left(\beta-\alpha\right)\approx 1$.
For the numerical analysis, we will set $\sin\left(\beta-\alpha\right)=1$, and
therefore have $R_{M}=M_{H^{0}}^{-2}$. In this case, the Higgs mass
eigenstates $h^{0}$ and $H^{0}$ coincide with $\phi^{0}$ and
$-\phi^{0\prime}$, respectively, and only the latter possesses non-SM-like
couplings to fermions. This reduces the number of relevant non-Yukawa-type
parameters of the extended Higgs sector to four, namely $\tan\beta$,
$M_{H^{+}}$, $M_{A^{0}}$ and $M_{\phi^{0\prime}}\equiv M_{H^{0}}$.

###  7.1 Constraints from $\mathbf{b\to s\gamma}$ decays

Figure 6: Sample Feynman diagram for the rare decay $b\to s\gamma$ at one loop
in the 2HDM. The cross indicates a chirality flip $t_{L}\to t_{R}$.

An important constraint on the magnitude of flavour-changing Yukawa couplings
in the up-type quark sector arises from the inclusive rare decays $B\to
X_{s}\gamma$. This process is mediated at the quark level by $b\to s\gamma$
through a top quark loop with a charged $W$ boson in the Standard Model and
receives additional contributions through charged-Higgs boson diagrams in the
2HDM, see Fig. 6. In order to eliminate the dependence on $V_{cb}\simeq-
V_{ts}$, one traditionally works with the quantity

$R_{\gamma}\equiv\frac{\mathcal{B}\left(b\to
s\gamma\right)+\mathcal{B}\left(b\to d\gamma\right)}{\mathcal{B}\left(b\to
cl\nu\right)}.$ (49)

Note that in the SM and the type-II 2HDM the branching ratio of $b\to d\gamma$
is much smaller than that of $b\to s\gamma$.

In the type-II 2HDM, the contribution of the charged Higgs diagrams to the
decay rate is positive, moving the theoretical prediction for the branching
ratio away from the experimental value averaged in Misiak:2017bgg  to

$R_{\gamma,\mathrm{exp.}}=\left(3.22\pm 0.15\right)\cdot 10^{-3}\,,$ (50)

where a lower cutoff of $E_{0}\geq$1.6\text{\,}\mathrm{GeV}$$ has been put on
the photon energy. The practical independence of $R_{\gamma}$ of $\tan\beta$
in the largest part of the parameter space in the type-II 2HDM allowed for the
extraction of a lower limit on the mass of the charged Higgs boson in
Misiak:2017bgg , see section 5. In the general 2HDM, where the factor
$\tan\beta$ from the $btH^{-}$ vertex need not be cancelled by a factor
$\cot\beta$ from the $stH^{-}$ vertex, the quantity $R_{\gamma}$ depends on
$\tan\beta$. The combination of the $\bar{t}_{L}b_{R}H^{+}$ and
$\bar{s}_{L}t_{R}H^{-}$ Yukawa couplings arising in the process is then given
by

$\left(\frac{m_{b}\tan\beta\,V_{tb}}{v}\right)\left(\frac{m_{t}\,V_{ts}^{*}}{v\tan\beta}+V_{ts}^{*}\,g_{tt}+V_{cs}^{*}\,g_{ct}\right)=\frac{m_{b}m_{t}\,V_{tb}V_{ts}^{*}}{v^{2}}\left(1+g_{st}^{\mathrm{eff}}\right)\,,$
(51)

where we have defined the short-hand notation

$g_{st}^{\mathrm{eff}}\equiv\frac{v\,\tan\beta}{m_{t}}\left(g_{ct}\frac{V_{cs}^{*}}{V_{ts}^{*}}+g_{tt}\right)\,.$
(52)

Note that $g_{st}^{\mathrm{eff}}$ carries a factor of $\tan\beta$. We stress
that the limit $g_{st}^{\mathrm{eff}}=0$ corresponds to the type-II 2HDM,
while the case $g_{st}^{\mathrm{eff}}=-1$ is the Standard Model limit. In
order to constrain our new flavour-changing couplings, we use the results from
Hermann:2012fc  and Misiak:2017bgg  with this trivial change of the
$\bar{s}_{L}t_{R}H^{-}$ Yukawa coupling.

For large enough values of $M_{H^{\pm}}$, the central value of $R_{\gamma}$ is
approximately given (at $\mu=\overline{m}_{t}\left(\overline{m}_{t}\right)$)
by

$\displaystyle R_{\gamma}\approx
10^{-4}\cdot\Bigg{\\{}\left.33.10\right|_{\mathrm{SM}}+$
$\displaystyle\left(1+\,\mathrm{Re}\,g_{st}^{\mathrm{eff}}\right)\left[r_{H}\left(-48.93-47.60\log
r_{H}-0.99\left(\log r_{H}\right)^{2}\right.\right.$
$\displaystyle\left.\left.\left.-0.15\left(\log
r_{H}\right)^{3}+4.71\,\mathrm{Li}_{2}\left(1-\frac{1}{r_{H}}\right)\right)\right.\right.$
$\displaystyle\left.\left.+r_{H}^{2}\left(-53.82-98.18\log
r_{H}+4.79\,\mathrm{Li}_{2}\left(1-\frac{1}{r_{H}}\right)\right)\right.\right.$
$\displaystyle\left.+r_{H}^{3}\left(-56.04-150.43\log
r_{H}+3.17\,\mathrm{Li}_{2}\left(1-\frac{1}{r_{H}}\right)\right)\right]\Bigg{\\}}\,,$
(53)

which agrees with the exact result within $1\,\%$ in the complete subdomain of
$\left(M_{H^{\pm}},\mathrm{Re}\,g_{st}^{\mathrm{eff}}\right)\in\left[$500\text{\,}\mathrm{GeV}$,\infty\right)\times\left[-5,5\right]$
in which $R_{\gamma}$ lies within the band allowed by the experimental and
theoretical uncertainties (see below). The approximate formula in Eq. (7.1)
only includes the interference of the new-physics contribution with the SM
result and neglects the squared new-physics contribution. In deriving Eq.
(7.1) we have used Eq. (10) of Ref. Misiak:2015xwa . We adopt the estimate of
the theoretical uncertainty of about $6.73\,\%$ given in Misiak:2017bgg ,
consisting of individual uncertainties of $5\,\%$ (non-perturbative),
$1.5\,\%$ (parametric), $3\,\%$ (higher-order), and $3\,\%$ (interpolation in
the charm quark mass $m_{c}$). In Fig. 7, we illustrate the ratio $R_{\gamma}$
in the $\left(M_{H^{\pm}},\mathrm{Re}\,g_{st}^{\mathrm{eff}}\right)$ plane.
The thick dashed line corresponds to the central experimental value in Eq.
(50) and $R_{\gamma,\mathrm{exp}}\pm\Delta_{\mathrm{exp}+\mathrm{th}}$, where
$\Delta_{\mathrm{exp}+\mathrm{th}}=2\sigma_{\mathrm{exp}}+\delta_{\mathrm{th}}$
is the sum of the experimental $2\sigma$ uncertainty intervals and the
theoretical uncertainty of $6.73\,\%$. For small values of $M_{H^{\pm}}$, the
allowed range of $g_{st}^{\mathrm{eff}}$ is tightly constrained around
$g_{st}^{\mathrm{eff}}=-1$, due to the closeness of the Standard Model
prediction and the experimental central value. At larger $M_{H^{\pm}}$, a
significantly wider range of flavour-changing Yukawa couplings is allowed.
From $R_{\gamma}$ we derive $M_{H^{\pm}}$-dependent upper and lower bounds on
$\mathrm{Re}\left(g_{st}^{\mathrm{eff}}\right)$ which we will use in order to
constrain the real part of the flavour-changing up-type couplings appearing in
$B_{s}\to\mu^{+}\mu^{-}$.

Figure 7: The ratio $R_{\gamma}$ as a function of the charged Higgs boson mass
$M_{H^{\pm}}$ and $g_{st}^{\mathrm{eff}}$. The dashed lines correspond to the
experimental central value and the $\pm 2\sigma$ intervals, to which we have
also added the theoretical uncertainties.

Figure 8: Exclusion limits in the $M_{A^{0}}-\tan\beta$ plane from different
searches performed by the ATLAS experiment. Taken from ATLAS:2022yvl . The
shown limits are for a different variant of the type-II 2HDM than ours, but
qualitatively apply as well to our model.

###  7.2 Higgs searches

Searches for heavy Higgs bosons at the LHC put powerful lower limits on the
masses of new Higgs particles, but these depend on the Yukawa sector of the
considered 2HDM. For the case of the type-II model (or ramifications of it)
these searches already exclude a significant portion of the
$M_{A^{0}}-\tan\beta$ plane at present. In Fig. 8, we show a recent collection
of exclusion limits obtained through various different searches by the ATLAS
experiment ATLAS:2022yvl  in a modification of the type-II model. These limits
imply that large values of $\tan\beta$ can only be realised if at the same
time the additional Higgs bosons are quite heavy,
$M_{A^{0}}\gtrsim$1.5\text{\,}\mathrm{TeV}$$.

It should be noted that for heavy $M_{A^{0}}$ these bounds also approximately
apply for $M_{H^{+}}$ and $M_{H^{0}}$, since the masses become degenerate in
the heavy-Higgs limit (for a thorough analysis of the allowed mass splittings
see Cacchio:2016qyh ). Note that the bounds from neutral Higgs searches are
more constraining than those for $H^{+}$ searches ATLAS:2018gfm . The most
stringent constraints are from final states with $\tau$’s and also apply to
our lepton Yukawa Lagrangian in Eq. (26) which we have chosen of type-II.
(However, the choice of the $\tau$ Yukawa couplings is of no relevance for the
phenomenology of the $b\to s$ FCNC processes discussed in this paper.) The
presence of the additional couplings $g_{jk}^{u}$ will increase some branching
ratios at the expense of those of decays into $\tau$’s, so that these bounds
can be somewhat weakened, but the general trend remains valid. Our scenarios
discussed below comply with the ATLAS bounds.

###  7.3 $\mathbf{B_{s}-\bar{B}_{s}}$ mixing

The effective $\bar{s}_{L}b_{R}A^{0}$ and $\bar{s}_{L}b_{R}H^{0}$ vertices are
the sum of the tree-level couplings $\propto g_{sb}$ and the loop
contributions involving vertex and self-energy diagrams $\propto
g_{ct}V_{tb}V_{cs}^{*}+g_{tt}V_{tb}V_{ts}^{*}$. $C_{P}$ and $C_{S}$ are
proportional to this effective vertex, which is therefore the relevant
quantity constrained by $\mathcal{B}(B_{s}\to\mu^{+}\mu^{-})$. Now the
$B_{s}-\bar{B}_{s}$ mixing amplitude depends quadratically on this effective
vertex, while both quantities drop quadratically with $M_{A^{0}}$ and
$M_{H^{0}}$, so that the quantities give complementary information on the
parameter space. $B_{s}-\bar{B}_{s}$ mixing is mediated in the Standard Model
by the $\Delta B=2$ operator
$Q_{VLL}=\left(\bar{b}_{L}\gamma_{\mu}s_{L}\right)\left(\bar{b}_{L}\gamma^{\mu}s_{L}\right)$
and in the 2HDM by three additional effective operators

$\displaystyle Q_{SLL}=$
$\displaystyle\left(\bar{b}_{R}s_{L}\right)\left(\bar{b}_{R}s_{L}\right)\,,$
(54) $\displaystyle Q_{SLR}=$
$\displaystyle\left(\bar{b}_{R}s_{L}\right)\left(\bar{b}_{L}s_{R}\right)\,,$
$\displaystyle Q_{SRR}=$
$\displaystyle\left(\bar{b}_{L}s_{R}\right)\left(\bar{b}_{L}s_{R}\right)\,,$

The Wilson coefficient of $Q_{SLL}$only involves the effective vertex
governing $B_{s}\to\mu^{+}\mu^{-}$ and the tree-level propagators of $A^{0}$
and $H^{0}$, while the coefficients of the other operators also involve the
corresponding chirality-flipped effective $\bar{s}_{R}b_{L}A^{0}$ and
$\bar{s}_{R}b_{L}H^{0}$ vertices. The coefficient $C_{SLL}$ ($C_{SRR}$) is
proportional to $m_{b}^{2}$ ($m_{s}^{2}$), while the coefficient $C_{SLR}$ is
proportional to $m_{b}m_{s}$; we therefore drop the operator $Q_{SRR}$ in our
analysis. However, considering the ATLAS constraints implying large masses
$M_{A^{0},H^{0}}$, $C_{SLL}\propto 1/M_{A^{0}}^{2}-1/M_{H^{0}}^{2}$ is heavily
suppressed due to the small splitting of Higgs masses in the large-mass limit,
and $C_{SLR}$ becomes relevant. This feature results from the fact that
$Q_{SLL}$ violates hypercharge by two units of the Higgs hypercharge and the
coefficient $C_{SLL}$ is therefore suppressed by a factor of
$v^{2}/M_{A^{0}}^{2}$ compared to $C_{SLR}$ multiplying $Q_{SLR}$ which
conserves hypercharge and SU(2). As a consequence, $C_{SLR}$ is more important
than $C_{SLL}$. This feature has been widely studied in the context of the
effective 2HDM emerging from integrating out superpartners in the MSSM
Isidori:2001fv ; Buras:2002wq ; Buras:2002vd ; Gorbahn:2009pp .

In the following we correlate $\mathcal{B}(B_{s}\to\mu^{+}\mu^{-})$ with
$\mathcal{B}(B\to X_{s}\gamma)$ and the $B_{s}-\bar{B}_{s}$ oscillation
frequency $\Delta M_{B_{s}}$, which is proportional to the magnitude of the
$B_{s}-\bar{B}_{s}$ mixing amplitude. Only the effective
$\bar{s}_{L}b_{R}A^{0}$ and $\bar{s}_{L}b_{R}H^{0}$ vertices are physical, by
e.g. changing the renormalisation condition for $g_{sb}$ we can shift pieces
between $g_{sb}$ and the renormlised loop. For simplicity, we set the tree-
value of $g_{sb}$ to zero, i.e. consider the case that this coupling is only
generated radiatively. The only non-trivial Yukawa structure entering the
considered observables is then $g_{st}^{\mathrm{eff}}$.

In Fig. 9 and Fig. 10 we illustrate the dependence of the branching ratio
$\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)$ on the flavour-changing
Yukawa couplings for some exemplary numerical values. We choose to show our
results separately for the LHCb LHCb:2021vsc ; LHCb:2021awg  and CMS
CMS:2022dbz  measurements, which lie on different sides of the SM prediction
calculated with the $|V_{cb}|$ from inclusive decays. The interval allowed on
the horizontal axis (the real part
$\mathrm{Re}\left(g_{st}^{\mathrm{eff}}\right)$) is for each choice of Higgs
masses determined by the range allowed by $b\to s\gamma$. Recall that
$g_{st}^{\mathrm{eff}}$ carries a factor of $\tan\beta$, which needs to be
taken into account if the bounds on $g_{st}^{\mathrm{eff}}$ were to be
converted into direct bounds on $g_{ct}$ and $g_{tt}$. A priori generic
$\leq{\cal O}(0.1)$ values of $g_{ct}$ could make $g_{st}^{\mathrm{eff}}$ in
Eq. (52) as large as ${\cal O}(50)$. For the considered Higgs masses such
large values of $|g_{st}^{\mathrm{eff}}|$ are forbidden by $b\to s\gamma$ and
the range for $\mathrm{Im}\left(g_{\mathrm{eff},bs\gamma}\right)$ shown on the
$y$ axis in plots (a) to (d) only serves the purpose to show the constraint
from $\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)$. Low values for Higgs
masses enforce small values for $\tan\beta$ from the LHC searches, but for
these scenarios $b\to s\gamma$ forbids any measurable effect in
$\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)$. This situation changes if
one considers large values for $\tan\beta$ and the Higgs masses, because
$\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)$ grows faster with $\tan\beta$
than $\mathcal{B}(B\to X_{s}\gamma)$, see plots (e) and (f). We find that in
the considered scenarios the bounds from $\mathcal{B}(B\to X_{s}\gamma)$ are
stronger than those from $\Delta M_{B_{s}}$, which we always find in the band
allowed by the current theoretical uncertainty. This situation changes if we
go to even larger Higgs masses, permitting larger values of
$|g_{st}^{\mathrm{eff}}|$ in $\mathcal{B}(B\to X_{s}\gamma)$. The quadratic
dependence of $\Delta M_{B_{s}}$ on $g_{st}^{\mathrm{eff}}$ makes $\Delta
M_{B_{s}}$ a good probe of the parameter region in which both the Higgs masses
and $g_{st}^{\mathrm{eff}}$ are large.

We do not study the case of sizable imaginary parts of the new FCNC couplings
here. These imaginary parts impact CP asymmetries such as $A_{CP}^{\rm
mix}(B_{s}\to J/\psi\phi)$, which will be investigated in a follow-up paper.
All numerical SM input parameters are given in Tab. 1. Note that we do not
take into account uncertainties in the input parameters as the dependence of
the branching ratio is much weaker than the dependence due to the variation of
flavour-changing Yukawa couplings.

We stress here that the presented calculation equally applies to
$B_{d}\to\mu^{+}\mu^{-}$ with the change $V_{ts}\to V_{td}$, but this is not
true anymore when considering the constraint from $B_{d}-\bar{B}_{d}$ mixing
because of the additional dependence on the light quark mass. Since $m_{d}$ is
negligibly small, the $B_{d}-\bar{B}_{d}$ mixing amplitude becomes insensitive
to the effective $\bar{b}dA^{0}$ and $\bar{b}dH^{0}$ couplings.

(a) $M_{H^{\pm}}=M_{H^{0}}=M_{A^{0}}=$600\text{\,}\mathrm{GeV}$$,
$\tan\beta=7$.

(b) $M_{H^{\pm}}=M_{A^{0}}=$800\text{\,}\mathrm{GeV}$$,
$M_{H^{0}}=$600\text{\,}\mathrm{GeV}$$, $\tan\beta=7$.

(c) $M_{H^{\pm}}=M_{H^{0}}=M_{A^{0}}=$1200\text{\,}\mathrm{GeV}$$,
$\tan\beta=12$.

(d) $M_{H^{\pm}}=$1200\text{\,}\mathrm{GeV}$$,
$M_{H^{0}}=M_{A^{0}}=$1100\text{\,}\mathrm{GeV}$$, $\tan\beta=12$.

(e) $M_{H^{\pm}}=M_{H^{0}}=M_{A^{0}}=$2000\text{\,}\mathrm{GeV}$$,
$\tan\beta=60$.

(f) $M_{H^{\pm}}=$1950\text{\,}\mathrm{GeV}$$,
$M_{H^{0}}=$2050\text{\,}\mathrm{GeV}$$,
$M_{A^{0}}=$2000\text{\,}\mathrm{GeV}$$, $\tan\beta=60$.

Figure 9: Branching ratio of $B_{s}\to\mu^{+}\mu^{-}$ for different values of
Higgs masses and $\tan\beta$. All quantities are evaluated at
$\mu=\overline{m}_{t}\left(\overline{m}_{t}\right)$. The red dashed and dotted
lines indicate the experimental central value and the $2\sigma$ uncertainties
of LHCb branching ratio measurement LHCb:2021vsc ; LHCb:2021awg ,
respectively. The constraint of $\mathcal{B}(B\to X_{s}\gamma)$ on
$|\mathrm{Im}\,g_{st}^{\mathrm{eff}}|$ is not shown.

(a) $M_{H^{\pm}}=M_{H^{0}}=M_{A^{0}}=$600\text{\,}\mathrm{GeV}$$,
$\tan\beta=7$.

(b) $M_{H^{\pm}}=M_{A^{0}}=$800\text{\,}\mathrm{GeV}$$,
$M_{H^{0}}=$600\text{\,}\mathrm{GeV}$$, $\tan\beta=7$.

(c) $M_{H^{\pm}}=M_{H^{0}}=M_{A^{0}}=$1200\text{\,}\mathrm{GeV}$$,
$\tan\beta=12$.

(d) $M_{H^{\pm}}=$1200\text{\,}\mathrm{GeV}$$,
$M_{H^{0}}=M_{A^{0}}=$1100\text{\,}\mathrm{GeV}$$, $\tan\beta=12$.

(e) $M_{H^{\pm}}=M_{H^{0}}=M_{A^{0}}=$2000\text{\,}\mathrm{GeV}$$,
$\tan\beta=60$.

(f) $M_{H^{\pm}}=$1950\text{\,}\mathrm{GeV}$$,
$M_{H^{0}}=$2050\text{\,}\mathrm{GeV}$$,
$M_{A^{0}}=$2000\text{\,}\mathrm{GeV}$$, $\tan\beta=60$.

Figure 10: Same as Fig. 9 for the CMS branching ratio measurement CMS:2022dbz
.

parameter | numerical value | reference
---|---|---
$M_{B_{s}}$ | $5.367\text{\,}\mathrm{GeV}$ | LHCb:2019rmd
$m_{b}$ | $4.18\text{\,}\mathrm{GeV}$ | Zyla:2020zbs
$m_{s}$ | $0.093\text{\,}\mathrm{GeV}$ | Zyla:2020zbs
$m_{\mu}$ | $0.106\text{\,}\mathrm{GeV}$ | Mohr:2015ccw
$\overline{m}_{t}\left(\overline{m}_{t}\right)$ | $162.622\text{\,}\mathrm{GeV}$ | Zyla:2020zbs
$\alpha_{s}\left(M_{Z}=$91.1876\text{\,}\mathrm{GeV}$\right)$ | $0.1179$ | Zyla:2020zbs
$G_{F}$ | $1.166\text{\times}{10}^{-5}\text{\,}\mathrm{G}\mathrm{eV}^{-2}$ | MuLan:2010shf
$M_{W}$ | $80.37\text{\,}\mathrm{GeV}$ | ATLAS:2017rzl
$f_{B_{s}}$ | $0.230\text{\,}\mathrm{GeV}$ | Aoki:2021kgd
$\Gamma_{H}^{s}$ | $\left($1.616\text{\,}\mathrm{ps}$\right)^{-1}$ | HFLAV
$\Gamma_{L}^{s}$ | $\left($1.427\text{\,}\mathrm{ps}$\right)^{-1}$ | HFLAV
$\absolutevalue{V_{ts}}$ | $0.041$ | CKMfitter
$\absolutevalue{V_{tb}}$ | $0.999$ | CKMfitter
$\absolutevalue{V_{cs}}$ | $0.974$ | CKMfitter
$\absolutevalue{V_{cb}}$ | $0.041$ | CKMfitter
$\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)_{\mathrm{SM}}$ | $3.65(23)\text{\times}{10}^{-9}\text{\,}$ | Bobeth:2013uxa
$\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)_{\mathrm{LHCb}}$ | $3.09_{-0.44}^{+0.48}\times 10^{-9}$ | LHCb:2021vsc ; LHCb:2021awg
$\mathcal{B}\left(B_{s}\to\mu^{+}\mu^{-}\right)_{\mathrm{CMS}}$ | $3.95_{-0.47}^{+0.52}\times 10^{-9}$ | CMS:2022dbz

Table 1: Numerical input used for the phenomenological analysis. For the QCD
running of the top-quark mass and the renormalisation scale, we have used
RunDec Chetyrkin:2000yt ; Herren:2017osy . The numerical values of the CKM
matrix elements have been taken from the updates provided on the CKMfitter
Charles:2004jd  web page. The numerical values of the $B_{s}$ decay widths
have been taken from the online updates provided at the HFLAV web page
HFLAV:2019otj .

##  8 Summary

We have presented a 2HDM with three flavour-breaking spurions in the quark
Yukawa sector. The model contains the established type-I and type-II models as
limiting cases and otherwise permits large FCNC couplings in the up-type
sector while naturally suppressing FCNC effects in down sector, as required by
data. Despite the large number of parameters the model makes characteristic
predictions, such as correlations between $b\to s\mu^{+}\mu^{-}$, $b\to
s\gamma$, and $B_{s}-\bar{B}_{s}$ mixing, all of which involve the same
combination $g_{st}^{\mathrm{eff}}$ of fundamental parameters. Also all
couplings to right-handed down-type quarks are proportional to the quark
masses as in the type-II model.

We have studied in detail the rare decay $B_{s}\to\mu^{+}\mu^{-}$, calculated
the Wilson coefficients of the effective operators $Q_{S}$ and $Q_{P}$, and
demonstrated the consistency of the model by showing that the UV counterterms
follow the pattern of the spurion expansion. Next we have calculated next-to-
leading order (two-loop) QCD corrections to this process to (i) verify that
higher-order QCD corrections can be correctly included (e.g. all UV
divergences could be renormalised in the usual way plus countertems for our
new couplings) and (ii) tame the sizable renormalisation-scale dependence of
the Yukawa couplings. Then we have studied the phenomenology of an FCNC
coupling $g_{ct}$ of the heavy neutral Higgs bosons to top and charm quarks.
We have found that —contrary to expectation— the dominant contribution to the
loop-induced $\bar{s}_{L}b_{R}A^{0}$ and $\bar{s}_{L}b_{R}H^{0}$ couplings do
not come from vertex diagrams with the neutral Higgs coupling to the internal
top-charm line, but from charged-Higgs couplings which inherit the dependence
on $g_{ct}$ through SU(2) symmetry. The corresponding diagram (FCNC self-
energy with the neutral Higgs attached to the $b$ line) is enhanced by a
factor of $\tan\beta$ compared to the vertex diagram, resulting in a
$\mathcal{O}\left(\tan^{3}\beta\right)$ contribution to the
$B_{s}\to\mu^{+}\mu^{-}$ amplitude which is absent in the type-II model. This
feature makes $B_{s}\to\mu^{+}\mu^{-}$ a sensitive probe of the model even for
Higgs masses well above the lower bounds found by the LHC experiments. For
small Higgs masses, however, $b\to s\gamma$ precludes large effects in
$B_{s}\to\mu^{+}\mu^{-}$. In our model the dominant contribution to
$B_{s}-\bar{B}_{s}$ mixing is naturally small due to a suppression factor of
$m_{s}m_{b}/v^{2}$; nevertheless $B_{s}-\bar{B}_{s}$ mixing sets constraints
on the parameter space for very large Higgs masses.

###### Acknowledgements.

M.L. would like to thank Florian Herren, Syuhei Iguro, Lucas Kunz and Mustafa
Tabet for helpful discussions. The authors are grateful to Matthias
Steinhauser for collaboration in an early stage of the project. This research
was supported the BMBF grants 05H18VKCC1 and 05H21VKKBA as well as Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) within the
Collaborative Research Center _Particle Physics Phenomenology after the Higgs
Discovery (P3H)_ (project no. 396021762 – TRR 257). All Feynman diagrams
within this paper have been produced using TikZ-FeynHand Dohse:2018vqo ;
Ellis:2016jkw .

## References

  * (1) G.C. Branco, P.M. Ferreira, L. Lavoura, M.N. Rebelo, M. Sher and J.P. Silva, _Theory and phenomenology of two-Higgs-doublet models_ , _Phys. Rept._ 516 (2012) 1 [1106.0034].
  * (2) J.F. Gunion, H.E. Haber, G.L. Kane and S. Dawson, _The Higgs Hunter’s Guide_ , vol. 80 (2000).
  * (3) J.E. Kim, _Light Pseudoscalars, Particle Physics and Cosmology_ , _Phys. Rept._ 150 (1987) 1.
  * (4) R.D. Peccei and H.R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * (5) R.D. Peccei and H.R. Quinn, _Constraints Imposed by CP Conservation in the Presence of Instantons_ , _Phys. Rev. D_ 16 (1977) 1791.
  * (6) H.E. Haber and G.L. Kane, _The Search for Supersymmetry: Probing Physics Beyond the Standard Model_ , _Phys. Rept._ 117 (1985) 75.
  * (7) A. Crivellin, C. Greub and A. Kokulu, _Explaining $B\to D\tau\nu$, $B\to D^{*}\tau\nu$ and $B\to\tau\nu$ in a 2HDM of type III_, _Phys. Rev. D_ 86 (2012) 054014 [1206.2634].
  * (8) A. Crivellin, A. Kokulu and C. Greub, _Flavor-phenomenology of two-Higgs-doublet models with generic Yukawa structure_ , _Phys. Rev. D_ 87 (2013) 094031 [1303.5877].
  * (9) T. Hermann, M. Misiak and M. Steinhauser, _Three-loop QCD corrections to $B_{s}\to\mu^{+}\mu^{-}$_, _JHEP_ 12 (2013) 097 [1311.1347].
  * (10) C. Bobeth, M. Gorbahn and E. Stamou, _Electroweak Corrections to $B_{s,d}\to\ell^{+}\ell^{-}$_, _Phys. Rev. D_ 89 (2014) 034023 [1311.1348].
  * (11) K. De Bruyn, R. Fleischer, R. Knegjens, P. Koppenburg, M. Merk and N. Tuning, _Branching Ratio Measurements of $B_{s}$ Decays_, _Phys. Rev. D_ 86 (2012) 014027 [1204.1735].
  * (12) K. De Bruyn, R. Fleischer, R. Knegjens, P. Koppenburg, M. Merk, A. Pellegrino et al., _Probing New Physics via the $B^{0}_{s}\to\mu^{+}\mu^{-}$ Effective Lifetime_, _Phys. Rev. Lett._ 109 (2012) 041801 [1204.1737].
  * (13) H.E. Haber and D. O’Neil, _Basis-independent methods for the two-Higgs-doublet model. II. The Significance of tan $\beta$_, _Phys. Rev. D_ 74 (2006) 015018 [hep-ph/0602242].
  * (14) D. Eriksson, J. Rathsman and O. Stål, _2HDMC: Two-Higgs-Doublet Model Calculator Physics and Manual_ , _Comput. Phys. Commun._ 181 (2010) 189 [0902.0851].
  * (15) G.C. Branco, L. Lavoura and J.P. Silva, _CP Violation_ , vol. 103 (1999).
  * (16) M. Misiak and M. Steinhauser, _Weak radiative decays of the B meson and bounds on $M_{H^{\pm}}$ in the Two-Higgs-Doublet Model_, _Eur. Phys. J. C_ 77 (2017) 201 [1702.04571].
  * (17) H.E. Logan and U. Nierste, _$B_{s,d}\to\ell^{+}\ell^{-}$ in a two Higgs doublet model_, _Nucl. Phys. B_ 586 (2000) 39 [hep-ph/0004139].
  * (18) C. Bobeth, M. Gorbahn, T. Hermann, M. Misiak, E. Stamou and M. Steinhauser, _$B_{s,d}\to l^{+}l^{-}$ in the Standard Model with Reduced Theoretical Uncertainty_, _Phys. Rev. Lett._ 112 (2014) 101801 [1311.0903].
  * (19) T. Banks, _Supersymmetry and the Quark Mass Matrix_ , _Nucl. Phys. B_ 303 (1988) 172.
  * (20) L.J. Hall, R. Rattazzi and U. Sarid, _The Top quark mass in supersymmetric SO(10) unification_ , _Phys. Rev. D_ 50 (1994) 7048 [hep-ph/9306309].
  * (21) M. Carena, M. Olechowski, S. Pokorski and C.E.M. Wagner, _Electroweak symmetry breaking and bottom - top Yukawa unification_ , _Nucl. Phys. B_ 426 (1994) 269 [hep-ph/9402253].
  * (22) K.S. Babu and C.F. Kolda, _Higgs mediated $B^{0}\to\mu^{+}\mu^{-}$ in minimal supersymmetry_, _Phys. Rev. Lett._ 84 (2000) 228 [hep-ph/9909476].
  * (23) G. Isidori and A. Retico, _Scalar flavor changing neutral currents in the large tan beta limit_ , _JHEP_ 11 (2001) 001 [hep-ph/0110121].
  * (24) A.J. Buras, P.H. Chankowski, J. Rosiek and Ł. Sławianowska, _$\Delta M_{d,s},B^{0}{d,s}\to\mu^{+}\mu^{-}$ and $B\to X_{s}\gamma$ in supersymmetry at large $\tan\beta$_, _Nucl. Phys. B_ 659 (2003) 3 [hep-ph/0210145].
  * (25) A.J. Buras, P.H. Chankowski, J. Rosiek and Ł. Sławianowska, _Correlation between $\Delta M_{s}$ and $B^{0}_{s,d}\rightarrow\mu^{+}\mu^{-}$ in Supersymmetry at large $\tan\beta$_, _Phys. Lett._ B546 (2002) 96 [hep-ph/0207241].
  * (26) L. Hofer, U. Nierste and D. Scherer, _Resummation of tan-beta-enhanced supersymmetric loop corrections beyond the decoupling limit_ , _JHEP_ 10 (2009) 081 [0907.5408].
  * (27) M. Gorbahn, S. Jäger, U. Nierste and S. Trine, _The supersymmetric Higgs sector and $B-\bar{B}$ mixing for large tan $\beta$_, _Phys. Rev. D_ 84 (2011) 034030 [0901.2065].
  * (28) A. Dedes, H.K. Dreiner and U. Nierste, _Correlation of $B_{s}\to\mu^{+}\mu^{-}$ and (g-2) ($\mu$) in minimal supergravity_, _Phys. Rev. Lett._ 87 (2001) 251804 [hep-ph/0108037].
  * (29) U. Nierste, M. Tabet and R. Ziegler, _Cornering Spontaneous $CP$ Violation with Charged-Higgs-Boson Searches_, _Phys. Rev. Lett._ 125 (2020) 031801 [1912.11501].
  * (30) R.S. Chivukula and H. Georgi, _Composite Technicolor Standard Model_ , _Phys. Lett. B_ 188 (1987) 99.
  * (31) G. D’Ambrosio, G.F. Giudice, G. Isidori and A. Strumia, _Minimal flavor violation: An Effective field theory approach_ , _Nucl. Phys. B_ 645 (2002) 155 [hep-ph/0207036].
  * (32) P. Tuzon and A. Pich, _The Aligned two-Higgs Doublet model_ , _Acta Phys. Polon. Supp._ 3 (2010) 215 [1001.0293].
  * (33) O. Eberhardt, A.P.n. Martínez and A. Pich, _Global fits in the Aligned Two-Higgs-Doublet model_ , _JHEP_ 05 (2021) 005 [2012.09200].
  * (34) N.D. Christensen, P. de Aquino, C. Degrande, C. Duhr, B. Fuks, M. Herquet et al., _A Comprehensive approach to new physics simulations_ , _Eur. Phys. J. C_ 71 (2011) 1541 [0906.2474].
  * (35) A. Alloul, N.D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250 [1310.1921].
  * (36) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, _UFO - The Universal FeynRules Output_ , _Comput. Phys. Commun._ 183 (2012) 1201 [1108.2040].
  * (37) M. Gerlach, F. Herren and M. Lang, _tapir: A tool for topologies, amplitudes, partial fraction decomposition and input for reductions_ , _Comput. Phys. Commun._ 282 (2023) 108544 [2201.05618].
  * (38) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279.
  * (39) B. Ruijl, T. Ueda and J. Vermaseren, _FORM version 4.2_ , 1707.06453.
  * (40) R. Harlander, T. Seidensticker and M. Steinhauser, _Complete corrections of Order alpha alpha-s to the decay of the Z boson into bottom quarks_ , _Phys. Lett. B_ 426 (1998) 125 [hep-ph/9712228].
  * (41) T. Seidensticker, _Automatic application of successive asymptotic expansions of Feynman diagrams_ , in _6th International Workshop on New Computing Techniques in Physics Research: Software Engineering, Artificial Intelligence Neural Nets, Genetic Algorithms, Symbolic Algebra, Automatic Calculation_ , 5, 1999 [hep-ph/9905298].
  * (42) A.I. Davydychev and J.B. Tausk, _Two loop selfenergy diagrams with different masses and the momentum expansion_ , _Nucl. Phys. B_ 397 (1993) 123.
  * (43) J.E. Salomon, _Das Laufen und die Vereinigung der Eichkopplungen des Standardmodells zur Drei-Schleifen-Ordnung_ , Ph.D. thesis, KIT, Karlsruhe, 2012\.  10.5445/IR/1000030116.
  * (44) T. Inami and C.S. Lim, _Effects of Superheavy Quarks and Leptons in Low-Energy Weak Processes $K_{L}\to\mu\bar{\mu}$, $K^{+}\to\pi+\nu\bar{\nu}$ and $K^{0}\leftrightarrow\bar{K}^{0}$_, _Progress of Theoretical Physics_ 65 (1981) 297 [https://academic.oup.com/ptp/article-pdf/65/1/297/5252099/65-1-297.pdf].
  * (45) G. Buchalla and A.J. Buras, _QCD corrections to the sdZ vertex for arbitrary top quark mass_ , _Nuclear Physics B_ 398 (1993) 285.
  * (46) G. Buchalla and A.J. Buras, _QCD corrections to rare K- and B-decays for arbitrary top quark mass_ , _Nuclear Physics B_ 400 (1993) 225.
  * (47) G. Buchalla and A.J. Buras, _The rare decays $K\to\pi\nu\bar{\nu}$, $B\to X\nu\bar{\nu}$ and $B\to l^{+}l^{-}$: An Update_, _Nucl. Phys. B_ 548 (1999) 309 [hep-ph/9901288].
  * (48) M. Misiak and J. Urban, _QCD corrections to FCNC decays mediated by Z penguins and W boxes_ , _Phys. Lett. B_ 451 (1999) 161 [hep-ph/9901278].
  * (49) M. Beneke, C. Bobeth and R. Szafron, _Enhanced electromagnetic correction to the rare $B$-meson decay $B_{s,d}\to\mu^{+}\mu^{-}$_, _Phys. Rev. Lett._ 120 (2018) 011801 [1708.09152].
  * (50) M. Beneke, C. Bobeth and R. Szafron, _Power-enhanced leading-logarithmic QED corrections to $B_{q}\to\mu^{+}\mu^{-}$_, _JHEP_ 10 (2019) 232 [1908.07011].
  * (51) C. Bobeth, T. Ewerth, F. Krüger and J. Urban, _Analysis of neutral Higgs boson contributions to the decays $\bar{B}$( $s^{)}\to\ell^{+}\ell^{-}$ and $\bar{B}\to K\ell^{+}\ell^{-}$_, _Phys. Rev. D_ 64 (2001) 074014 [hep-ph/0104284].
  * (52) C. Bobeth, A.J. Buras, F. Krüger and J. Urban, _QCD corrections to $\bar{B}\to X_{d,s}\nu\bar{\nu}$, $\bar{B}_{d,s}\to\ell^{+}\ell^{-}$, $K\to\pi\nu\bar{\nu}$ and $K_{L}\to\mu^{+}\mu^{-}$ in the MSSM_, _Nucl. Phys. B_ 630 (2002) 87 [hep-ph/0112305].
  * (53) T. Hermann, _Seltene B-Meson-Zerfälle innerhalb und jenseits des Standardmodells_ , Ph.D. thesis, KIT, Karlsruhe, 2014.  10.5445/IR/1000041138.
  * (54) LHCb collaboration, _Analysis of Neutral B-Meson Decays into Two Muons_ , _Phys. Rev. Lett._ 128 (2022) 041801 [2108.09284].
  * (55) LHCb collaboration, _Measurement of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ decay properties and search for the $B^{0}\to\mu^{+}\mu^{-}$ and $B^{0}_{s}\to\mu^{+}\mu^{-}\gamma$ decays_, _Phys. Rev. D_ 105 (2022) 012010 [2108.09283].
  * (56) CMS collaboration, _Measurement of ${\rm B^{0}_{s}}\to\mu^{+}\mu^{-}$ decay properties and search for the ${\rm B}^{0}\to\mu\mu$ decay in proton-proton collisions at $\sqrt{s}=13~{}\rm{TeV}$, report CMS-PAS-BPH-21-006_, .
  * (57) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * (58) HFLAV collaboration, _Averages of b-hadron, c-hadron, and $\tau$-lepton properties as of 2018_, _Eur. Phys. J. C_ 81 (2021) 226 [1909.12524].
  * (59) M. Blanke and S. Iguro, _Collider probe of heavy additional Higgs bosons solving the muon $g-2$_, 2210.13508.
  * (60) M. Blanke, A. Crivellin, S. de Boer, T. Kitahara, M. Moscati, U. Nierste et al., _Impact of polarization observables and $B_{c}\to\tau\nu$ on new physics explanations of the $b\to c\tau\nu$ anomaly_, _Phys. Rev. D_ 99 (2019) 075006 [1811.09603].
  * (61) M. Blanke, A. Crivellin, T. Kitahara, M. Moscati, U. Nierste and I. Nišandžić, _Addendum to “Impact of polarization observables and $B_{c}\to\tau\nu$ on new physics explanations of the $b\to c\tau\nu$ anomaly”_, 1905.08253.
  * (62) LHCb collaboration, “First joint measurement of $r(d^{*})$ and $r(d^{0})$ at lhcb.” https://indico.cern.ch/event/1187939/.
  * (63) S. Iguro, T. Kitahara and R. Watanabe, _Global fit to $b\to c\tau\nu$ anomalies 2022 mid-autumn_, 2210.10751.
  * (64) M. Fedele, M. Blanke, A. Crivellin, S. Iguro, T. Kitahara, U. Nierste et al., _Impact of $\Lambda_{b}\to\Lambda_{c}\tau\nu$ measurement on New Physics in $b\to c\,l\nu$ transitions_, 2211.14172.
  * (65) ATLAS collaboration, _Search for charged Higgs bosons decaying via $H^{\pm}\to\tau^{\pm}\nu_{\tau}$ in the $\tau$+jets and $\tau$+lepton final states with 36 fb-1 of $pp$ collision data recorded at $\sqrt{s}=13$ TeV with the ATLAS experiment_, _JHEP_ 09 (2018) 139 [1807.07915].
  * (66) ATLAS collaboration, _Summary plots for beyond Standard Model Higgs boson benchmarks for direct and indirect searches, ATLAS UB note ATL-PHYS-PUB-2022-043._ , .
  * (67) ATLAS collaboration, _Combined measurements of Higgs boson production and decay using up to $80$ fb-1 of proton-proton collision data at $\sqrt{s}=$ 13 TeV collected with the ATLAS experiment_, _Phys. Rev. D_ 101 (2020) 012002 [1909.02845].
  * (68) CMS collaboration, _Combined measurements of Higgs boson couplings in proton–proton collisions at $\sqrt{s}=13\,\text{Te}\text{V}$_, _Eur. Phys. J. C_ 79 (2019) 421 [1809.10733].
  * (69) T. Hermann, M. Misiak and M. Steinhauser, _$\bar{B}\to X_{s}\gamma$ in the Two Higgs Doublet Model up to Next-to-Next-to-Leading Order in QCD_, _JHEP_ 11 (2012) 036 [1208.2788].
  * (70) M. Misiak et al., _Updated NNLO QCD predictions for the weak radiative B-meson decays_ , _Phys. Rev. Lett._ 114 (2015) 221801 [1503.01789].
  * (71) V. Cacchio, D. Chowdhury, O. Eberhardt and C.W. Murphy, _Next-to-leading order unitarity fits in Two-Higgs-Doublet models with soft $\mathbb{Z}_{2}$ breaking_, _JHEP_ 11 (2016) 026 [1609.01290].
  * (72) LHCb collaboration, _Observation of $B^{0}_{(s)}\to J/\psi p\overline{p}$ decays and precision measurements of the $B^{0}_{(s)}$ masses_, _Phys. Rev. Lett._ 122 (2019) 191804 [1902.05588].
  * (73) P.J. Mohr, D.B. Newell and B.N. Taylor, _CODATA Recommended Values of the Fundamental Physical Constants: 2014_ , _Rev. Mod. Phys._ 88 (2016) 035009 [1507.07956].
  * (74) MuLan collaboration, _Measurement of the Positive Muon Lifetime and Determination of the Fermi Constant to Part-per-Million Precision_ , _Phys. Rev. Lett._ 106 (2011) 041803 [1010.0991].
  * (75) ATLAS collaboration, _Measurement of the $W$-boson mass in pp collisions at $\sqrt{s}=7$ TeV with the ATLAS detector_, _Eur. Phys. J. C_ 78 (2018) 110 [1701.07240].
  * (76) Y. Aoki et al., _FLAG Review 2021_ , 2111.09849.
  * (77) K.G. Chetyrkin, J.H. Kühn and M. Steinhauser, _RunDec: A Mathematica package for running and decoupling of the strong coupling and quark masses_ , _Comput. Phys. Commun._ 133 (2000) 43 [hep-ph/0004189].
  * (78) F. Herren and M. Steinhauser, _Version 3 of RunDec and CRunDec_ , _Comput. Phys. Commun._ 224 (2018) 333 [1703.03751].
  * (79) CKMfitter Group collaboration, _CP violation and the CKM matrix: Assessing the impact of the asymmetric $B$ factories_, _Eur. Phys. J. C_ 41 (2005) 1 [hep-ph/0406184].
  * (80) M. Dohse, _TikZ-FeynHand: Basic User Guide_ , 1802.00689.
  * (81) J. Ellis, _TikZ-Feynman: Feynman diagrams with TikZ_ , _Comput. Phys. Commun._ 210 (2017) 103 [1601.05437].