aainstitutetext: Theoretische Physik 1, Center for Particle Physics Siegen
(CPPS), Universität Siegen, Walter-Flex-Str. 3, 57068 Siegen,
Germanybbinstitutetext: INFN, Laboratori Nazionali di Frascati, C.P. 13,
100044 Frascati, Italyccinstitutetext: Université Paris-Saclay, CNRS/IN2P3,
IJCLab, 91405 Orsay, France

#  Enhancing $B_{s}\to e^{+}e^{-}$ to an Observable Level in the Two-Higgs-
Doublet Model

Matthew Black b Alexis D. Plascencia a,c Gilberto Tetlalmatzi-Xolocotzi
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>siegen.de

###### Abstract

As a result of the helicity suppression effect, within the Standard Model the
rare decay channel $B_{s}\to e^{+}e^{-}$ has a decay probability which is five
orders of magnitude below current experimental limits. Thus, any observation
of this channel within the current or forthcoming experiments will give
unambiguous evidence of Physics Beyond the Standard Model. In this work, we
present for the first time a New Physics scenario in which the branching
fraction $\bar{\mathcal{B}}r(B_{s}\to e^{+}e^{-})$ is enhanced up to values
which saturate the current experimental bounds. More concretely, we study the
general Two-Higgs-Doublet Model (2HDM) with a pseudoscalar coupling to
electrons unsuppressed by the electron mass. Furthermore, we demonstrate how
this scenario can arise from a UV-complete theory of quark-lepton unification
that can live at a low scale. This latter step allows us to establish
correlations between $B_{s}\to e^{+}e^{-}$ and the lepton-flavour-violating
decays $\tau^{-}\to\mu^{-}e^{+}e^{-}$ and $\tau\to\mu\gamma$.

SI-HEP-2022-23

P3H-22-092

##  1 Introduction

The rare decays $B_{s}\to\ell^{+}\ell^{-}$ for $\ell=e,\mu,\tau$ are
characterized by interesting properties which make them quite special and
suitable to test the Standard Model (SM) and to search for New Physics (NP).
For instance within the SM these transitions are only possible as loop-induced
processes. Moreover, they are extremely clean since only leptons are present
in the final state and all of the non-perturbative hadronic effects are
contained in the decay constant of the initial $B_{s}$ meson. As a matter of
fact the $B$ meson decay constants are currently known with a precision of
less than $1\%$ Aoki:2021kgd ; Bazavov:2017lyh ; ETM:2016nbo ; Dowdall:2013tga
; Hughes:2017spc .

One of the particular features of the $B_{s}$ meson rare processes is that
their decay probability in the SM is proportional to the mass of the final
state lepton; this is the so-called helicity suppression effect. For muons in
the final state, this leads to a SM branching fraction of
$\bar{\mathcal{B}}r(B_{s}\rightarrow\mu^{+}\mu^{-})=(3.55\pm 0.10)\times
10^{-9}$ which in spite of being already quite small has being measured by
different experimental collaborations leading to a combined result which is in
good agreement with the theoretical determination LHCb:2021awg ; LHCb:2021vsc
; ATLAS:2018cur ; CMS:2022dbz .

Due to the tiny mass of the electrons, for the channel $B_{s}\to e^{+}e^{-}$
the helicity suppression is maximal. Indeed, in the SM we have $\bar{\cal
B}r(B_{s}\to e^{+}e^{-})=(8.30\pm 0.36)\times 10^{-14}$ which is about four
orders of magnitude below the corresponding value for
$B_{s}\to\mu^{+}\mu^{-}$. Consequently, for a long time the experimental
search for this channel was not pursued. As a matter of fact, until 2020 the
only experimental result available was the upper bound reported by the CDF
collaboration CDF:2009ssr , which was then updated by the LHCb experiment
LHCb:2020pcv  with the result

$\bar{\cal B}r(B_{s}\to e^{+}e^{-})<9.4\times 10^{-9}.$ (1)

Due to the large difference between the most recent experimental bounds on
$\bar{\mathcal{B}}r(B_{s}\to e^{+}e^{-})$ and the corresponding SM prediction
we can conclude that any observation of this channel in the foreseeable future
can only be a manifestation of physics beyond the SM. Following a model-
independent approach, in Fleischer:2017ltw  it was shown how the presence of
NP pseudoscalar interactions could enhance the SM decay probability up to
values which can potentially saturate the known experimental bounds. One of
the main requirements to achieve this effect is that the NP couplings should
not be proportional to the mass of the electron $m_{e}$. This then excludes
models where the coupling between the NP pseudoscalars and the final state
electrons is determined at leading order by the mass $m_{e}$.

In this work, we present a minimal extension of the SM based on the type-III
Two-Higgs-Doublet Model (2HDM), in which a second Higgs is introduced with the
same quantum numbers as the SM Higgs and both scalar doublets are coupled to
quarks and leptons. We show how this scenario gives enough freedom to obtain
couplings between the electrons and the relevant scalar and pseudoscalar
particles that allow us to achieve large enhancements on
$\bar{\mathcal{B}}r(B_{s}\to e^{+}e^{-})$ while obeying all the relevant
phenomenological constraints. The literature on 2HDM is vast, for reviews on
the topic we refer the reader to Refs. Gunion:1989we ; Branco:2011iw .

Furthermore, we study how the type-III 2HDM scenario with the properties
outlined above can arise from a UV theory of quark-lepton unification. J. Pati
and A. Salam Pati:1974yy  postulated the idea of matter unification in which
the SM quarks and leptons belong to the same multiplet and this approach
remains as one of the best-motivated frameworks for physics beyond the SM.
However, since the top quark Yukawa coupling is predicted to be the same as
the Dirac neutrino coupling, then the seesaw mechanism Minkowski:1977sc ;
Yanagida:1979as ; GellMann:1980vs ; Mohapatra:1979ia  requires the energy
scale associated with the theory to be very high $\sim 10^{14}$ GeV, making it
hard to be phenomenologically tested. Consequently, here we consider the
theory proposed in Ref. Perez:2013osa , which can be regarded as a low energy
limit of the original Pati-Salam scenario, in which neutrinos acquire their
mass through the inverse seesaw mechanism Mohapatra:1986aw ; Mohapatra:1986bd
and the theory can be realized at a low energy scale.

This paper is structured as follows. In Section 2, we overview the
experimental and theoretical status regarding the $B_{s}$ meson rare decays.
In Section 3, we discuss the general 2HDM and study the Wilson operators
generated in this model. In Section 4, we present the corresponding
predictions for $\bar{\mathcal{B}}r(B_{s}\to e^{+}e^{-})$ and study the
phenomenological constraints on the parameter space considering different
observables. In Section 5, we present the theoretical motivation from a theory
of quark-lepton unification. Finally, our conclusions are presented in Section
6.

##  2 $B_{s}$ Meson Rare Decays

In order to address the decays $B_{s}\rightarrow\ell^{+}\ell^{-}$ we consider
the following effective Hamiltonian

$\mathcal{H}_{\rm
eff}=-\frac{G_{F}V_{tb}V_{ts}^{*}\alpha}{\sqrt{2}\pi}\left[C^{\ell\ell}_{10}O^{\ell\ell}_{10}+C_{S}^{\ell\ell}O^{\ell\ell}_{S}+C_{P}^{\ell\ell}O^{\ell\ell}_{P}++C_{10^{\prime}}^{\ell\ell}O^{\ell\ell}_{10^{\prime}}+C_{S^{\prime}}^{\ell\ell}O^{\ell\ell}_{S^{\prime}}+C_{P^{\prime}}^{\ell\ell}O^{\ell\ell}_{P^{\prime}}\right]+{\rm
h.c.}\,,$ (2)

where

$\displaystyle\mathcal{O}^{\ell\ell}_{10}$
$\displaystyle=\left(\bar{s}\gamma_{\mu}P_{L}b\right)\left(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell\right),$
$\displaystyle\mathcal{O}^{\ell\ell}_{10^{\prime}}$
$\displaystyle=\left(\bar{s}\gamma_{\mu}P_{R}b\right)\left(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell\right),$
$\displaystyle\mathcal{O}^{\ell\ell}_{S}$
$\displaystyle=m_{b}\left(\bar{s}P_{R}b\right)\left(\bar{\ell}\ell\right),$
$\displaystyle\mathcal{O}^{\ell\ell}_{S^{\prime}}$
$\displaystyle=m_{b}\left(\bar{s}P_{L}b\right)\left(\bar{\ell}\ell\right),$
$\displaystyle\mathcal{O}^{\ell\ell}_{P}$
$\displaystyle=m_{b}\left(\bar{s}P_{R}b\right)\left(\bar{\ell}\gamma^{5}\ell\right),$
$\displaystyle\mathcal{O}^{\ell\ell}_{P^{\prime}}$
$\displaystyle=m_{b}\left(\bar{s}P_{L}b\right)\left(\bar{\ell}\gamma^{5}\ell\right),$
(3)

for $\ell=e,\mu,\tau$.

The description of the $B_{s}$ meson rare decays
$B_{s}\rightarrow\ell^{+}\ell^{-}$ is given in terms of two measurable
quantities which offer complementary information. The first one is the time-
integrated branching fraction Fleischer:2017ltw ; DeBruyn:2012wk

$\displaystyle\bar{\cal
B}r(B_{s}\rightarrow\ell^{+}\ell^{-})=\frac{1}{2}\int^{\infty}_{0}\Braket{\Gamma(B_{s}(t)\rightarrow\ell^{+}\ell^{-})}dt,$
(4)

and the second one is the effective lifetime

$\displaystyle\tau_{\ell\ell}$ $\displaystyle\equiv$
$\displaystyle\frac{\int^{\infty}_{0}t\Braket{\Gamma(B_{s}(t)\rightarrow\ell^{+}\ell^{-})}dt}{\int^{\infty}_{0}\Braket{\Gamma(B_{s}(t)\rightarrow\ell^{+}\ell^{-})}dt},$
(5)

which is equivalent to the observable

$\displaystyle\mathcal{A}^{\ell\ell}_{\Delta\Gamma_{s}}$ $\displaystyle=$
$\displaystyle\frac{1}{y_{s}}\frac{(1-y^{2}_{s})\tau_{\ell\ell}-(1+y^{2}_{s})\tau_{B_{s}}}{2\tau_{B_{s}}-(1-y^{2}_{s})\tau_{\ell\ell}}.$
(6)

In Eq. (6), $\tau_{B_{s}}$ refers to the lifetime of the $B_{s}$ meson. In
addition, the neutral $B_{s}$ mixing effects are accounted for by

$\displaystyle y_{s}$ $\displaystyle\equiv$
$\displaystyle\frac{\Delta\Gamma_{s}}{2\Gamma_{s}},$ (7)

where $\Delta\Gamma_{s}$ is the decay width difference between the $B_{s}$ and
$\bar{B}_{s}$ mesons. Moreover, the untagged rate is defined as

$\displaystyle\Braket{\Gamma(B_{s}(t)\rightarrow\ell^{+}\ell^{-})}$
$\displaystyle\equiv$
$\displaystyle\Gamma(B^{0}_{s}(t)\rightarrow\ell^{+}\ell^{-})+\Gamma(\bar{B}^{0}_{s}(t)\rightarrow\ell^{+}\ell^{-})$
$\displaystyle=$ $\displaystyle\Gamma(B_{s}\rightarrow\ell^{+}\ell^{-})_{\rm
prompt}\times
e^{-t/\tau_{B_{s}}}\Bigl{(}\cosh(y_{s}t/\tau_{B_{s}})+\mathcal{A}^{\ell\ell}_{\Delta\Gamma_{s}}\sinh(y_{s}t/\tau_{B_{s}})\Bigl{)},$

where

$\displaystyle\Gamma(B_{s}\rightarrow\ell^{+}\ell^{-})_{\rm prompt}$
$\displaystyle=$
$\displaystyle\frac{G^{2}_{F}\alpha^{2}}{16\pi^{3}}|V_{ts}V^{*}_{tb}|^{2}f^{2}_{B_{s}}M_{B_{s}}m^{2}_{\ell}\sqrt{1-4\frac{m^{2}_{\ell}}{M^{2}_{B_{s}}}}|C^{\rm
SM}_{10}|^{2}\Bigl{(}|P_{\ell\ell}|^{2}+|S_{\ell\ell}|^{2}\Bigl{)}.$

The functions $P_{\ell\ell}$ and $S_{\ell\ell}$ are given by

$\displaystyle P_{\ell\ell}$ $\displaystyle\equiv$
$\displaystyle\frac{C^{\ell\ell}_{10}-C^{\ell\ell}_{10^{\prime}}}{C^{\rm
SM}_{10}}+\frac{M^{2}_{B_{s}}}{2m_{\ell}}\Bigl{(}\frac{m_{b}}{m_{b}+m_{s}}\Bigl{)}\Bigl{[}\frac{C^{\ell\ell}_{P}-C^{\ell\ell}_{P^{\prime}}}{C^{\rm
SM}_{10}}\Bigl{]},$ $\displaystyle S_{\ell\ell}$ $\displaystyle\equiv$
$\displaystyle\sqrt{1-4\frac{m^{2}_{\ell}}{M^{2}_{B_{s}}}}\frac{M^{2}_{B_{s}}}{2m_{\ell}}\Bigl{(}\frac{m_{b}}{m_{b}+m_{s}}\Bigl{)}\Bigl{[}\frac{C^{\ell\ell}_{S}-C^{\ell\ell}_{S^{\prime}}}{C^{\rm
SM}_{10}}\Bigl{]}.$ (10)

In the SM,
$C^{\ell\ell}_{P}=C^{\ell\ell}_{P^{\prime}}=C^{\ell\ell}_{S}=C^{\ell\ell}_{S^{\prime}}=0$,
leading to

$\displaystyle P^{\rm SM}_{\ell\ell}=1,$ $\displaystyle S^{\rm
SM}_{\ell\ell}=0,$ (11)

thus the branching fraction simplifies to

$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow\ell^{+}\ell^{-})_{\rm
SM}=\frac{1}{1-y_{s}}\frac{G^{2}_{F}\alpha^{2}}{16\pi^{3}}\tau_{B_{s}}|V_{ts}V^{*}_{tb}|^{2}f^{2}_{B_{s}}M_{B_{s}}m^{2}_{\ell}\sqrt{1-4\frac{m^{2}_{\ell}}{M^{2}_{B_{s}}}}|C^{\rm
SM}_{10}|^{2}.$ (12)

For real Wilson coefficients, the theoretical branching fraction for the
process $B_{s}\rightarrow\ell^{+}\ell^{-}$ is

$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow\ell^{+}\ell^{-})$
$\displaystyle=$
$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow\ell^{+}\ell^{-})_{\rm
SM}\times\Bigl{[}|P_{\ell\ell}|^{2}+\frac{1-y_{s}}{1+y_{s}}|S_{\ell\ell}|^{2}\Bigl{]}.$

An analogous expression in terms of $P_{\ell\ell}$ and $S_{\ell\ell}$ can also
be written for $\tau_{\ell\ell}$. However, due to the equivalence with
$\mathcal{A}^{\ell\ell}_{\Delta\Gamma_{s}}$ we only provide an explicit
expression for the latter:

$\displaystyle\mathcal{A}^{\ell\ell}_{\Delta\Gamma_{s}}$ $\displaystyle=$
$\displaystyle\frac{|P_{\ell\ell}|^{2}-|S_{\ell\ell}|^{2}}{|P_{\ell\ell}|^{2}+|S_{\ell\ell}|^{2}},$
(14)

and finally, $\tau_{\ell\ell}$ can be obtained by applying Eq. (6).

As can be seen in Eq. (12), in the SM, the decay probability $\bar{\cal
B}r(B_{s}\rightarrow\ell^{+}\ell^{-})_{\rm SM}$ is proportional to the square
of the mass of the final state lepton $m^{2}_{\ell}$. Since muons and
electrons are particularly light, for
$\bar{\mathcal{B}}r(B_{s}\rightarrow\mu^{+}\mu^{-})_{\rm SM}$ and
$\bar{\mathcal{B}}r(B_{s}\rightarrow e^{+}e^{-})_{\rm SM}$ the masses
$m_{\mu}$ and $m_{e}$ respectively act as suppression factors. Then the SM
predictions for the branching fractions for the different rare decays are:

$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow e^{+}e^{-})_{\rm SM}$
$\displaystyle=$ $\displaystyle(8.30\pm 0.22)\times 10^{-14},$ (15)
$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow\mu^{+}\mu^{-})_{\rm SM}$
$\displaystyle=$ $\displaystyle(3.55\pm 0.10)\times 10^{-9},$ (16)
$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow\tau^{+}\tau^{-})_{\rm SM}$
$\displaystyle=$ $\displaystyle(7.52\pm 0.20)\times 10^{-7}.$ (17)

For the experimental value of
$\bar{\mathcal{B}}r(B_{s}\rightarrow\mu^{+}\mu^{-})$ we update the result
presented in Altmannshofer:2021qrr  by performing a weighted average including
the measurements from LHCb, ATLAS and the latest value from CMS LHCb:2021awg ;
LHCb:2021vsc ; ATLAS:2018cur ; CMS:2022dbz :

$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow\mu^{+}\mu^{-})_{\rm
Exp}=(3.39\pm 0.29)\times 10^{-9}.$ (18)

In addition, LHCb has performed two pioneering measurements of the effective
lifetime $\tau_{\mu\mu}$ LHCb:2017rmj ; LHCb:2021vsc ; LHCb:2021awg .
Combining LHCb:2021vsc  and CMS:2022dbz  we obtain

$\displaystyle\tau_{\mu\mu}$ $\displaystyle=$ $\displaystyle 1.83\pm
0.21~{}\rm{ps}.$ (19)

For examples of studies on $B_{s}\rightarrow\mu^{+}\mu^{-}$ and more
generically on $b\rightarrow s\ell^{+}\ell^{-}$ processes within the context
of 2HDM see for example Crivellin:2013wna ; Crivellin:2019dun .

In the case of $B_{s}\rightarrow\tau^{+}\tau^{-}$ the current $95\%$ C.L.
bound is available LHCb:2017myy :

$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow\tau^{+}\tau^{-})<6.8\times
10^{-3}.$ (20)

Notice that according to Eq. (17), in the SM the decay ratio
$\bar{\mathcal{B}}r(B_{s}\rightarrow\tau^{+}\tau^{-})$ has the largest value
amongst all final state leptons, however the reconstruction of the $\tau$
particle is a challenging task, making the experimental extraction of the
corresponding decay ratio especially difficult.

Finally, due to the tiny mass of the electron, the transition
$B_{s}\rightarrow e^{+}e^{-}$ is rather suppressed in the SM. Indeed, as can
be seen in Eq. (15), this channel has the lowest branching fraction among all
the possible leptonic final states and lies outside the reach of current or
future particle physics experiments. However, the presence of NP scalar and
pseudoscalar mediators can drastically enhance the value of
$\bar{\mathcal{B}}r(B_{s}\rightarrow e^{+}e^{-})$ Fleischer:2017ltw . As shown
in Eq. (2), this effects boils down to the presence of the tiny factor
$m_{\ell}=m_{e}$ in the denominators of the functions $P_{\ell\ell}=P_{ee}$
and $S_{\ell\ell}=S_{ee}$ which for non-zero contributions of the differences
$\Delta C^{\ell\ell}_{P}=C^{\ell\ell}_{P}-C^{\ell\ell}_{P^{\prime}}$ and
$\Delta C^{\ell\ell}_{S}=C^{\ell\ell}_{S}-C^{\ell\ell}_{S^{\prime}}$ can
maximally lift the helicity suppression induced in the SM. In this respect,
the decay channel $B_{s}\rightarrow e^{+}e^{-}$ is special since its
measurement in any foreseeable experimental facility will be a clear and
unambiguous indication of NP.

In 2009, CDF reported the first bound on the production rate of this
particular channel at $90\%$ C.L.:

$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow e^{+}e^{-})_{\rm
Exp,CDF}<2.8\times 10^{-7}.$ (21)

This bound was updated recently by the LHCb collaboration LHCb:2020pcv
leading to the following $90\,(95)\%$ C.L. bounds:

$\displaystyle\bar{\mathcal{B}}r(B_{s}\rightarrow e^{+}e^{-})_{\rm
Exp,LHCb}<9.4\,(11.2)\times 10^{-9}.$ (22)

As described previously, the potential enhancement on
$\bar{\mathcal{B}}r(B_{s}\rightarrow e^{+}e^{-})$ as the result of NP scalar
and pseudoscalar particles was first noticed in Fleischer:2017ltw  in a model-
independent fashion. To the best of our knowledge an analysis within the
context of a renormalizable NP framework has not been performed so far. In the
following sections we take this next step and develop a NP scenario where this
effect can arise.

In order to perform the numerical calculations corresponding to the
$B$-physics processes, in this work we will make use of the flavour physics
package flavio111https://flav-io.github.io/ Straub:2018kue . This will also
allow us to combine observables in frequentist likelihood fits of experimental
data to determine constraints on the parameters of our NP model. flavio
describes NP contributions model-independently using Effective Field Theories
(EFTs) where NP enters as additions to the Wilson coefficients of the
operators of the EFT. Of interest here is the Weak Effective Theory (WET) with
five active flavours (defined at the scale $m_{b}$), where we can directly
describe the contributions from our model in the language of the relevant
Wilson coefficients as laid out below.

##  3 The General 2HDM and the process $B_{s}\rightarrow e^{+}e^{-}$

In this Section, we will focus on a mechanism that lifts the helicity
suppression in $\bar{\mathcal{B}}r(B_{s}\to e^{+}e^{-})$ leading to a large
enhancement within a minimal extension of the SM. Namely, extending the SM
with a second Higgs doublet with the same quantum numbers as the SM one; for
reviews on the Two-Higgs-Doublet Model (2HDM) we refer the reader to Refs.
Gunion:1989we ; Branco:2011iw . In the general 2HDM, both Higgs doublets are
coupled to the quarks and leptons; this scenario is commonly referred to in
the literature as the type-III 2HDM. Therefore, we can write the following
Yukawa interactions,

$\displaystyle-\mathcal{L}$
$\displaystyle\supset\bar{Q}_{L}\left(Y^{u}_{1}\widetilde{H}_{1}+Y^{u}_{2}\widetilde{H}_{2}\right)u_{R}+\bar{Q}_{L}\left(Y^{d}_{1}H_{1}+Y^{d}_{2}H_{2}\right)d_{R}$
$\displaystyle+\bar{\ell}_{L}\left(Y^{e}_{1}H_{1}+Y^{e}_{2}H_{2}\right)e_{R}+{\rm
h.c.}\,,$ (23)

with $H_{1}^{T}=(H_{1}^{+},(v_{1}+H_{1}^{0}+iA_{1}^{0})/\sqrt{2})$,
$\widetilde{H}_{1}=i\sigma_{2}H_{1}^{*}$ and correspondingly for $H_{2}$. The
vacuum expectation values (vevs) are defined by $\langle
H_{1}^{0}\rangle=v_{1}$ and $\langle H_{2}^{0}\rangle=v_{2}$.

The scalar potential for $H_{1}$ and $H_{2}$ with quantum numbers
$(\mathbf{1},\mathbf{2},1/2)$ corresponds to

$\displaystyle V(H_{1},H_{2})$
$\displaystyle=m_{11}^{2}H_{1}^{\dagger}H_{1}+m_{22}^{2}H_{2}^{\dagger}H_{2}-m_{12}^{2}\left[\left(H_{1}^{\dagger}H_{2}\right)+{\rm
h.c.}\right]$
$\displaystyle+\frac{\lambda_{1}}{2}\left(H_{1}^{\dagger}H_{1}\right)^{2}+\frac{\lambda_{2}}{2}\left(H_{2}^{\dagger}H_{2}\right)^{2}+\lambda_{3}\left(H_{1}^{\dagger}H_{1}\right)\left(H_{2}^{\dagger}H_{2}\right)+\lambda_{4}\left(H_{1}^{\dagger}H_{2}\right)\left(H_{2}^{\dagger}H_{1}\right)$
$\displaystyle+\left[\frac{\lambda_{5}}{2}\left(H_{1}^{\dagger}H_{2}\right)^{2}+\lambda_{6}\left(H_{1}^{\dagger}H_{1}\right)\left(H_{1}^{\dagger}H_{2}\right)+\lambda_{7}\left(H_{2}^{\dagger}H_{2}\right)\left(H_{1}^{\dagger}H_{2}\right)+{\rm
h.c.}\right].$ (24)

The physical Higgs fields are defined by:

$\displaystyle\begin{pmatrix}H\\\ h\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\cos\alpha&\sin\alpha\\\
-\sin\alpha&\cos\alpha\end{pmatrix}\begin{pmatrix}H_{1}^{0}\\\
H_{2}^{0}\end{pmatrix},$ (25) $\displaystyle\begin{pmatrix}G\\\
A\end{pmatrix}$ $\displaystyle=\begin{pmatrix}\cos\beta&\sin\beta\\\
-\sin\beta&\cos\beta\end{pmatrix}\begin{pmatrix}A_{1}^{0}\\\
A_{2}^{0}\end{pmatrix},$ (26) $\displaystyle\begin{pmatrix}G^{\pm}\\\
H^{\pm}\end{pmatrix}$ $\displaystyle=\begin{pmatrix}\cos\beta&\sin\beta\\\
-\sin\beta&\cos\beta\end{pmatrix}\begin{pmatrix}H_{1}^{\pm}\\\
H_{2}^{\pm}\end{pmatrix},$ (27)

where $h$ is identified as the SM-like Higgs and $H$ as an additional neutral
Higgs. In addition, $H_{i}^{0},H_{i}^{\pm},A^{0}_{i}$ are the neutral, charged
and CP-odd components of the Higgs doublets respectively. Finally
$G,\,G^{\pm}$ are the would-be Goldstone bosons. The mixing angle $\beta$ is
defined by the ratio of the vevs of the Higgs doublets,
$\tan{\beta}=v_{2}/v_{1}$ and we use $v^{2}=v_{1}^{2}+v_{2}^{2}$. The
couplings of $h$ are SM-like in the alignment limit $\sin{(\beta-\alpha)}\to
1$, which corresponds to $\alpha=\beta-\pi/2$ . Thus, the interactions between
the fermions and the neutral scalars can be written as

$\displaystyle-\mathcal{L}\supset$
$\displaystyle\,\,\bar{f}^{i}_{L}\left[\frac{M^{i}_{\rm
diag}}{v}h+\left(-\cot\beta\frac{M^{i}_{\rm
diag}}{v}+\frac{\Omega^{i}}{\sqrt{2}s_{\beta}}\right)\left(H\pm
iA\right)\right]f^{i}_{R}+{\rm h.c.}\,,$ (28)

where the super index $i$ denotes the fermion flavour for $i=u,d,e$. In the
equation above, the positive sign is assigned to the field $A$ when
considering couplings to the up-type quarks while the negative sign is
considered for couplings to the down-type quarks and charged leptons. The mass
matrices are given by

$\displaystyle
m^{i}=Y_{1}^{i}\frac{v_{1}}{\sqrt{2}}+Y_{2}^{i}\frac{v_{2}}{\sqrt{2}},$ (29)

and $M^{i}_{\rm diag}$ in Eq. (28) correspond to the diagonal mass matrices
$M^{i}_{\rm diag}=V_{L}^{i\dagger}m^{i}V_{R}^{i}$ with unitary $V_{L}^{i}$ and
$V_{R}^{i}$. Finally, the matrices $\Omega^{i}$ are given by
$\Omega^{i}=V_{L}^{i\dagger}Y_{1}^{i}V_{R}^{i}$ and are characterized by
general components.

For the charged leptons we assume the interaction with the heavy Higgs bosons
to be very close to flavour diagonal:

$\widetilde{Y}^{\ell}=-\cot\beta\frac{M^{E}_{\rm
diag}}{v}+\frac{\Omega^{\ell}}{\sqrt{2}s_{\beta}}=\begin{pmatrix}y_{ee}&\varepsilon&\varepsilon\\\
\varepsilon&y_{\mu\mu}&\varepsilon\\\
\varepsilon&\varepsilon&y_{\tau\tau}\end{pmatrix},$ (30)

where $\varepsilon\ll y_{jj}$. This allows us to evade the strong experimental
constraints from the non-observation of processes which violate lepton flavour
such as $\mu\to e\gamma$ MEG:2016leq , $\mu-e$ conversion SINDRUMII:2006dvw
and $\mu\to eee$ SINDRUM:1987nra . Since we are mostly interested in the
coupling to electrons, we assume the hierarchy $y_{\mu\mu}\ll y_{ee}$ and
$\varepsilon\ll y_{\tau\tau}$. As we will discuss in Section 5, the texture in
Eq. (30) obeying the indicated hierarchy can be motivated by embedding the
2HDM in a low-energy limit of Pati-Salam unification.

Similarly for the down-type quarks, we assume the Yukawa interaction to be
close to flavour-diagonal:

$\widetilde{Y}^{d}=-\cot\beta\frac{M^{D}_{\rm
diag}}{v}+\frac{\Omega^{d}}{\sqrt{2}s_{\beta}}=\begin{pmatrix}y_{dd}&\varepsilon&\varepsilon\\\
\varepsilon&y_{ss}&y_{bs}/2\\\ \varepsilon&y_{bs}/2&y_{bb}\end{pmatrix},$ (31)

where we write $\varepsilon$ for very small numbers obeying $\varepsilon\ll
y_{ij}$. Here, we have suppressed some off-diagonal entries in order to avoid
the strong bounds coming from measurements of neutral kaon mixing. Also, we
have kept the off-diagonal entry $y_{bs}$ since this coupling mediates the
process $B_{s}\rightarrow e^{+}e^{-}$ at tree level by coupling the NP scalar
$H$ and pseudoscalar $A$ to the quarks $b$ and $s$. Moreover, we have assumed
that $\widetilde{Y}^{d}_{sb}=\widetilde{Y}^{d}_{bs}=y_{bs}/2$, a choice that
will be motivated in Section 5. As we shall see below, the experimental
constraint from $B_{s}$ meson mixing requires this coupling to be very small.

The relevant Yukawa interactions affecting the process $B_{s}\to e^{+}e^{-}$
at tree level are

$-\mathcal{L}\supset
y_{ee}\,\bar{e}eH+y_{bs}\,\bar{b}sH-i\,y_{ee}\,\bar{e}\gamma^{5}eA-i\,y_{bs}\,\bar{b}\gamma^{5}sA,$
(32)

where the assumption of no new sources of CP violation implies all couplings
to be real. The expression in Eq. (32) will play a central role in our
subsequent discussion.

By integrating out the particles $A$ and $H$ inside Eq. (32) we can
immediately determine the Wilson coefficients $C^{ee}_{S^{(^{\prime})}}$ and
$C^{ee}_{P^{(^{\prime})}}$ given in Eq. (2) in terms of the parameters of our
model:

$\displaystyle C^{ee}_{S}$
$\displaystyle=\frac{y_{ee}y_{bs}}{M_{H}^{2}}\left(\frac{\sqrt{2}\pi}{m_{b}G_{F}V_{tb}V_{ts}^{*}\alpha}\right),\hskip
36.98866ptC^{ee}_{S^{\prime}}=C^{ee}_{S},$ (33) $\displaystyle C^{ee}_{P}$
$\displaystyle=-\frac{y_{ee}y_{bs}}{M_{A}^{2}}\left(\frac{\sqrt{2}\pi}{m_{b}G_{F}V_{tb}V_{ts}^{*}\alpha}\right),\hskip
28.45274ptC^{ee}_{P^{\prime}}=-C^{ee}_{P},$ (34)

where $M_{H}$ and $M_{A}$ are the masses of $H$ and $A$ respectively.

Given that, according to Eq. (2) and Eq. (2), the branching fraction depends
only on the difference $\Delta C^{ee}_{S}=C^{ee}_{S}-C^{ee}_{S^{\prime}}$ and
that based on Eq. (33) within our model $C^{ee}_{S}=C^{ee}_{S^{\prime}}$, we
can immediately see that the CP-even scalar $H$ does not contribute to the
observable $\bar{\mathcal{B}}r(B_{s}\to e^{+}e^{-})$. However, the CP-odd
scalar (pseudoscalar) $A$ can have large effects. A crucial point to highlight
is that in our NP scenario, the coupling between the heavy Higgs bosons and
electrons is not required to be proportional to the mass of the electron; this
is of capital importance when lifting the SM helicity suppression.

The values of the masses $M_{H}$ and $M_{A}$ depend on the parameters of the
2HDM scalar potential $\lambda_{1\to 7},m^{2}_{12},\tan\beta$ and $\alpha$
(see Eq. (24)) Gunion:2002zf . Hence $M_{H}$ and $M_{A}$ are not independent
from each other and are actually correlated. The parameters in the scalar
potential are constrained by different theoretical conditions such as
perturbativity and vacuum stability which can be combined as follows
Gunion:2002zf

$\displaystyle 0$ $\displaystyle<\lambda_{1,2}~{}~{}~{}\;<4,$ (35)
$\displaystyle-\sqrt{\lambda_{1}\lambda_{2}}$
$\displaystyle<\lambda_{3}\,~{}~{}~{}~{}~{}<4,$ (36) $\displaystyle-4$
$\displaystyle<\lambda_{4,5,6,7}<4.$ (37)

Therefore, to determine the allowed values for $M_{A}$ and $M_{H}$, we
randomly sample through the parameter space of the scalar potential which in
addition delivers the masses of the charged scalars $M_{H^{\pm}}$. During this
procedure, we use the mass of the SM Higgs boson as a constraint, i.e. we take
$M_{h}=M_{h}^{\rm SM}=125.25\pm 0.17~{}\rm{GeV}$ Workman:2022ynf  as well as
the inequality $M_{h}\leq M_{H}$, and fix $\sin(\beta-\alpha)=1$.

In Fig. 1 we present our results for the allowed values of $M_{A}$ and $M_{H}$
after imposing these constraints on the parameters in the scalar potential.
For small masses below 1 TeV, the mass splitting between $M_{A}$ and $M_{H}$
can be quite large (around 500 GeV), while for heavy masses around 10 TeV,
this difference must be rather small (around 50 GeV), and hence the heavy mass
regime satisfies $M_{A}=M_{H}$ to a very good approximation. In this limit the
NP Wilson coefficients given in Eq. (33) satisfy $C^{ee}_{P}=-C^{ee}_{S}$ and
$C^{ee}_{P^{\prime}}=C^{ee}_{S^{\prime}}$ which are two well-known
relationships obtained in SMEFT.

Figure 1:  Correlation between the masses $M_{A}$ and $M_{H}$ from the
constraints on the 2HDM scalar potential. A total of $10^{9}$ points
fulfilling the conditions on the theory were generated, in the limit of
$\sin(\beta-\alpha)=1$ and defining the perturbativity limits as
$|\lambda_{i}|<4$. The region shaded in blue shows the masses for successful
parameters of the 2HDM scalar potential.

##  4 Enhancing $B_{s}\to e^{+}e^{-}$ and Phenomenological Constraints

Our next task is to determine bounds for the couplings $y_{bs}$ and $y_{ee}$
to quarks and leptons respectively based on the phenomenological constraints
available. We first focus on the bounds on the $y_{ee}$ coupling from the
measurement of the cross-section for $e^{-}e^{+}\to e^{-}e^{+}$ performed by
the LEP collaboration. Their reported constraints on the four-electron axial-
vector interaction ALEPH:2006bhb  can be translated to the scalar and
pseudoscalar interactions; we find that the 95% confidence level lower bound
is determined by

$\frac{y_{ee}^{2}}{M_{H}^{2}}+\frac{y_{ee}^{2}}{M_{A}^{2}}<\frac{1}{(4\,\,{\rm
TeV})^{2}}\,\,.$ (38)

In the case where $M_{H}\\!=\\!M_{A}\\!$ this bound becomes

$\displaystyle\frac{y_{ee}}{M_{H}}<\frac{1}{(5.7~{}{\rm TeV})}.$ (39)

Figure 2:  Feynman diagrams for $B_{s}\to e^{-}e^{+}$ and $B_{s}-\bar{B}_{s}$
mixing induced at tree-level by the scalar $H$ and pseudoscalar $A$,
respectively.

The observable $\Delta M_{s}$ is sensitive to the presence of NP scalar and
pseudoscalar particles and thus can impose strong constraints on the coupling
$y_{bs}$; the new tree-level diagrams mediated by $H$ and $A$ are shown in
Fig. 2. In this work we will consider the following set of $\Delta B=2$
operators which contribute to $\Delta M_{s}$:

$\displaystyle\mathcal{O}^{\Delta B=2}_{V}$
$\displaystyle=\bar{s}_{i}\gamma^{\mu}(1-\gamma_{5})b_{i}~{}\bar{s}_{j}\gamma_{\mu}(1-\gamma_{5})b_{j},$
$\displaystyle\mathcal{O}^{\Delta B=2}_{LL}$
$\displaystyle=\bar{s}_{i}(1-\gamma_{5})b_{i}~{}\bar{s}_{j}(1-\gamma_{5})b_{j},$
$\displaystyle\mathcal{O}^{\Delta B=2}_{RR}$
$\displaystyle=\bar{s}_{i}(1+\gamma_{5})b_{i}~{}\bar{s}_{j}(1+\gamma_{5})b_{j},$
$\displaystyle\mathcal{O}^{\Delta B=2}_{LR}$
$\displaystyle=\bar{s}_{i}(1-\gamma_{5})b_{i}~{}\bar{s}_{j}(1+\gamma_{5})b_{j},$
(40)

where in the SM only the coefficient of $\mathcal{O}^{\Delta B=2}_{V}$ is non-
zero. In terms of the parameters of our model in Eq. (32), the coefficients of
the operators $\mathcal{O}^{\Delta B=2}_{LL},\mathcal{O}^{\Delta B=2}_{RR}$
and $\mathcal{O}^{\Delta B=2}_{LR}$ are respectively

$\displaystyle C^{\Delta B=2}_{RR}$
$\displaystyle=\frac{y^{2}_{bs}}{4}\Bigl{[}\frac{1}{m^{2}_{H}}-\frac{1}{m^{2}_{A}}\Bigl{]},$
$\displaystyle C^{\Delta B=2}_{RR}$ $\displaystyle=C^{\Delta B=2}_{LL},$
$\displaystyle C^{\Delta B=2}_{LR}$
$\displaystyle=\frac{y^{2}_{bs}}{2}\Bigl{[}\frac{1}{m^{2}_{H}}+\frac{1}{m^{2}_{A}}\Bigl{]}.$
(41)

The relevant matrix elements of the operators in Eq. (40) are given by
DiLuzio:2019jyq

$\displaystyle\langle\mathcal{O}^{\Delta B=2}_{V}\rangle$ $\displaystyle=$
$\displaystyle\frac{8}{3}M^{2}_{B_{s}}f^{2}_{B_{s}}B_{1}(\mu_{b}),$
$\displaystyle\langle\mathcal{O}^{\Delta B=2}_{LL}\rangle$ $\displaystyle=$
$\displaystyle
M^{2}_{B_{s}}f^{2}_{B_{s}}\frac{-5M^{2}_{B_{s}}}{3(\bar{m}_{b}(\mu_{b})+\bar{m}_{s}(\mu_{b}))^{2}}B_{2}(\mu_{b}),$
$\displaystyle\langle\mathcal{O}^{\Delta B=2}_{LL}\rangle$ $\displaystyle=$
$\displaystyle\langle\mathcal{O}^{\Delta B=2}_{RR}\rangle,$
$\displaystyle\langle\mathcal{O}^{\Delta B=2}_{LR}\rangle$ $\displaystyle=$
$\displaystyle
M^{2}_{B_{s}}f^{2}_{B_{s}}\Bigl{[}\frac{2M^{2}_{B_{s}}}{(\bar{m}_{b}(\mu_{b})+\bar{m}_{s}(\mu_{b}))^{2}}+\frac{1}{3}\Bigl{]}B_{4}(\mu_{b}).$

The observable $\Delta M_{s}$ is calculated according to

$\displaystyle\Delta M_{s}=2|M^{s}_{12}|,$ (42)

where

$\displaystyle M^{s}_{12}$ $\displaystyle=$
$\displaystyle\frac{G^{2}_{F}}{12\pi^{2}}\lambda^{2}_{t}M^{2}_{W}S_{0}(x_{t})\hat{\eta_{B}}f^{2}_{B_{s}}M_{B_{s}}B_{1}+\frac{1}{2M_{B_{s}}}\Bigl{[}2C^{\Delta
B=2}_{RR}\langle\mathcal{O}^{\Delta B=2}_{RR}\rangle+C^{\Delta
B=2}_{LR}\langle\mathcal{O}^{\Delta B=2}_{LR}\rangle\Bigl{]}.$

To estimate $\Delta M_{s}$ we use flavio. Our inputs are the Bag parameters
given in DiLuzio:2019jyq  and the values for
$|V_{us}|,\,|V_{cb}|,\,|V_{ub}|,\gamma$ from the CKMfitter’s Spring ‘21 update
Charles:2004jd . Thus, our determination in the SM is

$\displaystyle\Delta M^{\rm SM}_{s}$ $\displaystyle=$ $\displaystyle 17.49\pm
0.64\,{\rm ps}^{-1}.$ (44)

This result is in agreement with previous calculations, but its central value
is noticeably lower in comparison; consider for instance the result reported
in DiLuzio:2019jyq  which reads $\Delta M^{\rm
SM}_{s}=18.4^{+0.7}_{-1.2}\,$ps-1. This deviation is induced mainly by the
update on the CKM inputs, more specifically by the $\sim 1\sigma$ decrease in
$|V_{cb}|$ between the results from the CKMfitter’s Summer ‘18 report and the
one from Spring ‘21 Charles:2004jd . The experimental result for $\Delta
M_{s}$ is taken from Amhis:2022mac :

$\displaystyle\Delta M^{\rm Exp}_{s}$ $\displaystyle=$ $\displaystyle
17.765\pm 0.006\,{\rm ps}^{-1}.$ (45)

Figure 3:  Allowed parameter space of the quark coupling $y_{bs}$ and new
neutral Higgs mass $M_{A}$, from the measurement of mass mixing in the $B_{s}$
system, $\Delta M_{s}$. Left panel: In the limit of $M_{A}=M_{H}$. Right
panel: We allow $M_{A}$ and $M_{H}$ to differ within the theoretical
constraints of the model and minimize through $M_{H}$. In both plots, the
contours in dark and light blue represent the allowed space within $1\sigma$
and $2\sigma$ respectively.

In Fig. 3, we present the constraint from the measured value of neutral
$B_{s}$ meson mixing, $\Delta M_{s}$, for the allowed parameter space in the
$M_{H}$ vs $y_{bs}$ plane, both in the limit of $M_{A}=M_{H}$ and allowing the
maximum freedom between $M_{A}$ and $M_{H}$ from theory; the difference
between these two scenarios is found to be minimal. This plot shows that in
order to be in agreement with the measured value of $\Delta M_{s}$, the
coupling $y_{bs}$ has to be small; e.g. for a mass of $M_{A}\\!=\\!1\,$TeV we
find that $y_{bs}\lesssim 0.001$ at $2\sigma$.

Furthermore, we can use the processes $B\rightarrow K^{(*)}e^{+}e^{-}$ to
constrain simultaneously the couplings $y_{bs}$ and $y_{ee}$. As a matter of
fact, the NP effects in the transitions $B\rightarrow K^{(*)}e^{+}e^{-}$ can
be parameterized directly in terms of the Wilson coefficients
$C^{ee}_{S^{(^{\prime})}},\,C^{ee}_{P^{(^{\prime})}}$ in Eq. (34) which also
affect $B_{s}\to e^{+}e^{-}$. The $B\to K^{(*)}e^{+}e^{-}$ observables
considered in this work are listed in Table 1.

Observable |  $q^{2}$ bin (GeV2) | Exp. Avg. | SM Pred.
---|---|---|---
$10^{8}\times\frac{\Delta{\cal B}}{\Delta q^{2}}(B^{+}\to K^{+}e^{+}e^{-})$ | [1.0,6.0] |  $3.24\pm 0.65$ LHCb:2014vgu ; BELLE:2019xld  | $3.37\pm 0.56$
| [0.1,4.0] |  $4.70\pm 1.01$ BELLE:2019xld  | $3.40\pm 0.58$
| [4.0,8.12] |  $2.36\pm 0.79$ BELLE:2019xld  | $3.31\pm 0.54$
$10^{7}\times\frac{\Delta{\cal B}}{\Delta q^{2}}(B^{0}\to K^{*0}e^{+}e^{-})$ | [0.003,1.0] |  $3.09\pm 0.99$ LHCb:2013pra  | $2.10\pm 0.35$
$P_{4}^{\prime}(B\to K^{*}e^{+}e^{-})$ | [1.0,6.0] |  $-0.71\pm 0.40$ Belle:2016fev  | $-0.34\pm 0.04$
| [14.18,19.0] |  $-0.15\pm 0.41$ Belle:2016fev  | $-0.63\pm 0.01$
$P_{5}^{\prime}(B\to K^{*}e^{+}e^{-})$ | [1.0,6.0] |  $-0.23\pm 0.41$ Belle:2016fev  | $-0.42\pm 0.09$
| [14.18,19.0] |  $-0.86\pm 0.34$ Belle:2016fev  | $-0.63\pm 0.03$
Table 1: List of $B\to Ke^{+}e^{-}$ observables used to constrain the
couplings $y_{bs}$ and $y_{ee}$ and the masses $M_{H}$ and $M_{A}$. For
$P_{4}^{\prime}$ and $P_{5}^{\prime}$, we consider the average of the $B^{+}$
and $B^{0}$ modes.

Since the associated expressions for the observables in Table 1 are quite
lengthy, we refer the interested reader to the flavio’s documentation and code
Straub:2018kue . Here we only quote explicitly the NP components of the
helicity amplitudes for a pseudoscalar $K$ or vector $K^{*}$ final state kaon
which depend on the Wilson coefficients $C^{ee}_{P^{(^{\prime})}}$ and
$C^{ee}_{S^{(^{\prime})}}$

$\displaystyle h_{S}^{K^{*}}$
$\displaystyle=\frac{i\lambda(m_{B}^{2},m_{K^{*}}^{2},q^{2})}{2}\Big{(}C^{ee}_{S}-C^{ee}_{S^{\prime}}\Big{)}\,A_{0}(q^{2}),$
(46) $\displaystyle h_{P}^{K^{*}}$
$\displaystyle=\frac{i\lambda(m_{B}^{2},m_{K^{*}}^{2},q^{2})}{2}\Big{(}C^{ee}_{P}-C^{ee}_{P^{\prime}}+\dots\Big{)}\,A_{0}(q^{2}),$
(47) $\displaystyle h_{S}^{K}$
$\displaystyle=\frac{m_{B}^{2}-m_{K}^{2}}{2}\Big{(}C^{ee}_{S}+C^{ee}_{S^{\prime}}\Big{)}\,f_{0}(q^{2}),$
(48) $\displaystyle h_{P}^{K}$
$\displaystyle=\frac{m_{B}^{2}-m_{K}^{2}}{2}\Big{(}C^{ee}_{P}+C^{ee}_{P^{\prime}}+\dots\Big{)}\,f_{0}(q^{2}),$
(49)

where the ellipses stand for extra contributions including the purely SM ones
in $C_{10^{(^{\prime})}}^{ee}$. Moreover, $\lambda(a,b,c)$ is the Källen
function, and $A_{0}(q^{2}),\,f_{0}(q^{2})$ are each one of the $B\to K^{*}$
and $B\to K$ form factors respectively which depend on the invariant dilepton
mass squared $q^{2}$ and are constructed using Bailey:2015dka ; Horgan:2015vla
; Bharucha:2015bzk ; Gubernari:2018wyi . In Eqs. (46) and (47) we can see that
the NP contributions enter in terms of the differences of $\Delta
C^{ee}_{S}=C^{ee}_{S}-C^{ee}_{S^{\prime}}$ and $\Delta
C^{ee}_{P}=C^{ee}_{P}-C^{ee}_{P^{\prime}}$ as is the case for $B_{s}\to
e^{+}e^{-}$ and therefore the $B\to K^{*}$ modes will only be sensitive to
$M_{A}$. Conversely, from Eqs. (48) and (49), instead of the differences, the
NP effects enter in terms of the sum of the relevant Wilson coefficients and
so the $B\to K$ modes depend only on $M_{H}$.

Figure 4:  Allowed parameter space for the coupling product $y_{bs}y_{ee}$ and
the new neutral Higgs mass $M_{A}$, from the measurements of the $B\to
K^{(*)}e^{+}e^{-}$ observables in Table 1. Left panel: In the limit of
$M_{A}=M_{H}$. Right panel: We allow $M_{A}$ and $M_{H}$ to differ within the
theoretical constraints of the model and minimize through $M_{H}$. In both
plots, the contours in dark and light blue represent the possible space within
$1\sigma$ and $2\sigma$ respectively.

In Fig. 4, we show the constraints arising from the combined fit of the $B\to
K^{(*)}$ observables listed in Table 1, both in the limit of $M_{A}=M_{H}$ and
allowing the maximum freedom between $M_{A}$ and $M_{H}$ from theory. We can
see that the product of the couplings $y_{bs}y_{ee}$ is expected to be small
and is correlated with $M_{A}$ and $M_{H}$ similarly to the results drawn from
$\Delta M_{s}$ for $y_{bs}$.

Figure 5:  Allowed parameter space for the coupling product $y_{bs}y_{ee}$ and
the mass of the new neutral Higgs mass $M_{A}$. Left panel: In the limit of
$M_{A}=M_{H}$. Right panel: We allow $M_{A}$ and $M_{H}$ to differ within the
theoretical constraints of the model and minimize through $M_{H}$. In both
plots, the contours in dark and light blue represent the possible space within
$1\sigma$ and $2\sigma$ respectively. The black lines correspond to contours
for the ratio $\bar{\cal B}r(B_{s}\to e^{+}e^{-})/\bar{\cal B}r(B_{s}\to
e^{+}e^{-})_{\rm SM}$. Notice that to saturate the current experimental bound
it is required an enhancement of $10^{5}$.

In Fig. 5 we present the allowed parameter space in the $y_{bs}y_{ee}-M_{A}$
plane from all constraints considered, where once more, we take into account
two cases: first the limit $M_{A}=M_{H}$ and second the situation where the
maximum freedom between $M_{A}$ and $M_{H}$ is allowed from theory. The region
shaded in green is allowed by the LEP and the $\Delta M_{s}$ bounds. The
region shaded in blue is allowed by the bound from the $B\to
K^{(*)}e^{+}e^{-}$ processes. The black lines correspond to contours for
constant values of the ratio $\bar{\cal B}r(B_{s}\to e^{+}e^{-})/\bar{\cal
B}r(B_{s}\to e^{+}e^{-})_{\rm SM}$ which determine the enhancement in
$\bar{\cal B}r(B_{s}\to e^{+}e^{-})$ with respect to the SM prediction. Here
we can see how an enhancement by a factor as large as $10^{8}$ is allowed by
the collider and $B$ physics constraints. In fact, we can saturate the bound
imposed by the LHCb analysis of $B_{s}\to e^{+}e^{-}$ reported in LHCb:2020pcv
which is shown by the red dashed line and requires an enhancement by a factor
of $10^{5}$.

Thus Fig. 5 contains one of the main results of this work: within the context
of a type-III 2HDM, an enhancement on $\bar{\cal B}r(B_{s}\to e^{+}e^{-})$ up
to values which saturate the current experimental bounds is completely allowed
and consistent with the different phenomenological constraints known from $B$
physics and collider studies. Notice that in this Section we have focused on
determining the possible values that the coupling constants and masses of the
type-III 2HDM affecting directly $B_{s}\to e^{+}e^{-}$ can assume, although so
far we have not discussed how a model with such properties could arise from a
UV-complete theory. This is precisely the task we undertake in the next
section.

##  5 2HDM and Quark-Lepton Unification

In this Section, we demonstrate how the coupling structure we have considered
for the 2HDM can be obtained from a UV theory of quark-lepton unification that
can live at a low energy scale. We focus on the NP framework proposed in Ref.
Perez:2013osa  which is based on the gauge group

$\,{\rm SU}(4)_{C}\otimes\,{\rm SU}(2)_{L}\otimes\,{\rm U}(1)_{R}.$

Moreover, it implements the inverse seesaw mechanism in order to generate
neutrino masses and can be seen as a low energy limit of the Pati-Salam theory
Pati:1974yy . The phenomenology of the leptoquarks in this NP framework has
been studied in FileviezPerez:2021lkq ; FileviezPerez:2022rbk , while the
phenomenology of its scalar sector, which corresponds to a special case of the
type-III 2HDM, has been recently analyzed in FileviezPerez:2022fni . For
further details we refer the reader to those references.

Within our framework, the SM matter fields are unified in the following
representations,

$\displaystyle F_{QL}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}u_{r}&u_{g}&u_{b}&\nu\\\
d_{r}&d_{g}&d_{b}&e\end{array}\right)\sim(\mathbf{4},\mathbf{2},0),$ (52)
$\displaystyle F_{u}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}u^{c}_{r}&u^{c}_{g}&u^{c}_{b}&\nu^{c}\end{array}\right)\sim(\mathbf{\bar{4}},\mathbf{1},-1/2),$
(54) $\displaystyle F_{d}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}d^{c}_{r}&d^{c}_{g}&d^{c}_{b}&e^{c}\end{array}\right)\sim(\mathbf{\bar{4}},\mathbf{1},1/2).$
(56)

and hence, the leptons can be interpreted as the fourth colour of the
fermions. The Yukawa interactions for the charged fermions can be written as

$\displaystyle-\mathcal{L}_{Y}\supset
Y_{1}F_{QL}F_{u}H_{1}+Y_{2}F_{QL}F_{u}\Phi+Y_{3}H_{1}^{\dagger}F_{QL}F_{d}+Y_{4}\Phi^{\dagger}F_{QL}F_{d}+{\rm
h.c.},$ (57)

where $H_{1}\sim(\mathbf{1},\mathbf{2},1/2)$ and
$\Phi\sim(\mathbf{15},\mathbf{2},1/2)$ are required to generate fermion masses
in a consistent manner. The $\Phi$ field contains a second Higgs doublet
$H_{2}$ that is coupled to all the SM fermions

$\Phi=\left(\begin{array}[]{cc}\Phi_{8}&\Phi_{3}\\\ \Phi_{4}&0\\\
\end{array}\right)+\sqrt{2}\,T_{4}\ H_{2}\sim(\mathbf{15},\mathbf{2},1/2),$
(58)

where $T_{4}$ is one of the generators of $\,{\rm SU}(4)_{C}$ and it is
normalized as $T_{4}=\frac{1}{2\sqrt{6}}\rm{diag}(1,1,1,-3).$

Since the NP framework under consideration arises from quark-lepton
unification there are only four independent Yukawa couplings (instead of
eight) defining the interactions between the Higgs doublets and the SM
fermions:

$\displaystyle-\mathcal{L}$
$\displaystyle=\bar{u}_{R}\left(Y_{1}^{T}\widetilde{H}_{1}+\frac{1}{2\sqrt{3}}Y_{2}^{T}\widetilde{H}_{2}\right)Q_{L}+\bar{N}_{R}\left(Y_{1}^{T}\widetilde{H}_{1}-\frac{\sqrt{3}}{2}Y_{2}^{T}\widetilde{H}_{2}\right)\ell_{L}$
$\displaystyle+\bar{d}_{R}\left(Y_{3}^{T}H_{1}^{\dagger}+\frac{1}{2\sqrt{3}}Y_{4}^{T}H_{2}^{\dagger}\right)Q_{L}+\bar{e}_{R}\left(Y_{3}^{T}H_{1}^{\dagger}-\frac{\sqrt{3}}{2}Y_{4}^{T}H_{2}^{\dagger}\right)\ell_{L}+{\rm
h.c.}\,,$ (59)

and the vevs are defined by $\langle H_{1}^{0}\rangle=v_{1}/\sqrt{2}$ and
$\langle H_{2}^{0}\rangle=v_{2}/\sqrt{2}$.

As it was shown in Ref. FileviezPerez:2022fni , the interactions between the
physical Higgs bosons and the SM down-type quarks and charged leptons are
given respectively by

$\displaystyle\tilde{Y}^{\ell}$
$\displaystyle=\left(\tan{\beta}-3\cot{\beta}\right)\frac{M^{E}_{\rm
diag}}{4v}+3\left(\tan{\beta}+\cot{\beta}\right)\frac{V_{c}^{T}M^{D}_{\rm
diag}V}{4v},$ (60) $\displaystyle\tilde{Y}^{d}$
$\displaystyle=\left(3\tan{\beta}-\cot{\beta}\right)\frac{M^{D}_{\rm
diag}}{4v}+\left(\tan{\beta}+\cot{\beta}\right)\frac{V_{c}^{*}M^{E}_{\rm
diag}V^{\dagger}}{4v},$ (61)

where $V$ and $V_{c}$ are unitary matrices which contain information about the
unknown mixing between quarks and leptons. In addition, $M^{D}_{\rm diag}$ and
$M^{E}_{\rm diag}$ are the diagonal mass matrices for down-type quarks and
charged leptons. From Eqs. (60) and (61) above we can see that the theory
predicts a correlation between the couplings to quarks and leptons. As it was
demonstrated in Ref. FileviezPerez:2021arx , in the regimes with $\tan\beta\gg
1$ or $\tan\beta\ll 1$ the theory predicts unique relations among the decay
widths of heavy Higgs bosons that can be probed at the LHC. Consequently, we
focus on these two limits.

If we assume the complex phases to vanish, the $3\times 3$ unitary matrix $V$
can be parameterized in terms of three mixing angles, which here we denote as
$\theta_{12}$, $\theta_{13}$ and $\theta_{23}$, as follows

$V=\begin{pmatrix}c_{12}c_{13}&s_{12}c_{13}&s_{13}\\\\[4.30554pt]
-s_{12}c_{23}-c_{12}s_{23}s_{13}\,\,\,&c_{12}c_{23}-s_{12}s_{23}s_{13}&s_{23}c_{13}\\\\[4.30554pt]
s_{12}s_{23}-c_{12}c_{23}s_{13}&-c_{12}s_{23}-s_{12}c_{23}s_{13}\,\,\,&c_{23}c_{13}\end{pmatrix},$
(62)

where we have used $s_{ij}$ and $c_{ij}$ as short notation for
$\sin\theta_{ij}$ and $\cos\theta_{ij}$ respectively. An analogous expression
can then also be written for $V_{c}$ but with primed mixing angles
$s^{\prime}_{ij}$ and $c^{\prime}_{ij}$. For large $\tan\beta$ and in the
limit where $s_{ij}\to 1$ and $s^{\prime}_{ij}\to 1$ the interactions with the
charged leptons are simplified to

$\tilde{Y}^{\ell}=\frac{\tan\beta}{4v}\begin{pmatrix}m_{e}+3m_{b}&\varepsilon&\varepsilon\\\
\varepsilon&m_{\mu}+3m_{s}&\varepsilon\\\
\varepsilon&\varepsilon&m_{\tau}+3m_{d}\end{pmatrix},$ (63)

which gives us the flavour-diagonal couplings with the hierarchy $y_{ee}\gg
y_{\mu\mu},y_{\tau\tau}$. This motivates our choice for the couplings in
Section 3. The same conclusions hold for intermediate and small values of
$\tan\beta$.

For the down-type quarks and large $\tan\beta$, we get the following
interaction matrix

$\tilde{Y}^{d}=\frac{\tan\beta}{4v}\begin{pmatrix}3m_{d}+m_{\tau}&\varepsilon&\varepsilon\\\
\varepsilon&3m_{s}+m_{\mu}&\varepsilon\\\
\varepsilon&\varepsilon&3m_{b}+m_{e}\end{pmatrix},$ (64)

which gives us the hierarchy $y_{dd}\simeq y_{bb}\gg y_{ss}$. The same
conclusion holds for intermediate and small values of $\tan\beta$.
Unfortunately, due to the freedom in the coupling to the right-handed
neutrinos, the theory does not predict the coupling of the Higgs bosons to the
up-type quarks.

Figure 6:  Left panel: Predicted correlation between
$\bar{\mathcal{B}}r(B_{s}\to e^{+}e^{-})$ and
$\mathcal{B}r(\tau^{-}\to\mu^{-}e^{+}e^{-})$. The regions shaded in red
correspond to the experimental bounds for each decay respectively. The
different bands correspond to different values for the mass splitting $\Delta
M\equiv M_{A}-M_{H}$. We have fixed $s_{23}=s^{\prime}_{23}=0.98$ and
$\tan\beta=10$. Right panel: Same as the left panel for the predicted
correlation between $\bar{\mathcal{B}}r(B_{s}\to e^{+}e^{-})$ and
$\mathcal{B}r(\tau\to\mu\gamma)$.

Since we require a non-zero $y_{bs}$ coupling, we set all
$s_{ij}\\!=\\!s^{\prime}_{ij}\\!=\\!1$ except for $s_{23}$ and
$s^{\prime}_{23}$; this choice is motivated by requiring the first-generation
off-diagonal couplings to be vanishing. From quark-lepton unification, non-
zero entries for $\tilde{Y}^{d}_{bs}$ and $\tilde{Y}^{d}_{sb}$ imply non-zero
values for $\tilde{Y}^{\ell}_{\mu\tau}$ and $\tilde{Y}^{\ell}_{\tau\mu}$.
Namely,

$\displaystyle\tilde{Y}^{d}_{sb}$
$\displaystyle=\frac{1}{4v}\left(\tan\beta+\cot\beta\right)\left(m_{\mu}s^{\prime}_{23}c_{23}-m_{e}s_{23}c^{\prime}_{23}\right),$
(65) $\displaystyle\tilde{Y}^{d}_{bs}$
$\displaystyle=\frac{1}{4v}\left(\tan\beta+\cot\beta\right)\left(m_{\mu}s_{23}c^{\prime}_{23}-m_{e}s^{\prime}_{23}c_{23}\right),$
(66) $\displaystyle\tilde{Y}^{\ell}_{\mu\tau}$
$\displaystyle=\frac{3}{4v}\left(\tan\beta+\cot\beta\right)\left(m_{s}s^{\prime}_{23}c_{23}-m_{d}s_{23}c^{\prime}_{23}\right),$
(67) $\displaystyle\tilde{Y}^{\ell}_{\tau\mu}$
$\displaystyle=\frac{3}{4v}\left(\tan\beta+\cot\beta\right)\left(m_{s}s_{23}c^{\prime}_{23}-m_{d}s^{\prime}_{23}c_{23}\right),$
(68)

whenever $s_{23}=s^{\prime}_{23}$ then we have that
$\tilde{Y}^{d}_{sb}=\tilde{Y}^{d}_{bs}=y_{bs}/2$ which motivates the choice
made in Eq. (31); we also have that
$\tilde{Y}^{\ell}_{\mu\tau}=\tilde{Y}^{\ell}_{\tau\mu}=y_{\tau\mu}/2$. The
$\tau\mu$ couplings will generate the following dimension-six operators

$\displaystyle\mathcal{H}_{\rm eff}\supset$
$\displaystyle-\frac{y_{ee}y_{\mu\tau}}{M_{H}^{2}}(\bar{\tau}\mu)(\bar{e}e)-\frac{y_{\tau\tau}y_{\mu\tau}}{M_{H}^{2}}(\bar{\tau}\mu)(\bar{\tau}\tau)$
$\displaystyle+\frac{y_{ee}y_{\mu\tau}}{M_{A}^{2}}(\bar{\tau}\gamma^{5}\mu)(\bar{e}\gamma^{5}e)+\frac{y_{\tau\tau}y_{\mu\tau}}{M_{A}^{2}}(\bar{\tau}\gamma^{5}\mu)(\bar{\tau}\gamma^{5}\tau),$
(69)

which will induce the lepton-flavour-violating decay
$\tau^{\pm}\to\mu^{\pm}e^{-}e^{+}$ at tree-level, and
$\tau^{\pm}\to\mu^{\pm}\gamma$ at one-loop. The current experimental bounds on
these decay channels are ${\cal B}r(\tau\to\mu\gamma)<4.4\times 10^{-8}$
BaBar:2009hkt  and ${\cal B}r(\tau^{-}\to\mu^{-}\,e^{+}e^{-})<1.8\times
10^{-8}$ Hayasaka:2010np . These bounds are expected to be improved by future
$B$ factories SuperB:2010cqs .

The effective operators
$\bigl{\\{}(\bar{\tau}\mu)(\bar{e}e),(\bar{\tau}\gamma^{5}\mu)(\bar{e}\gamma^{5}e)\bigl{\\}}$
in Eq. (5) can be mapped to operators which are analogous to
$\mathcal{O}^{\Delta B=2}_{LL},\mathcal{O}^{\Delta B=2}_{RR}$ and
$\mathcal{O}^{\Delta B=2}_{LR}$ in Eq. (40), and hence, the corresponding
Wilson coefficients have also analogous expressions to the coefficients shown
in Eqs. (41). More precisely, the new set of operators to be considered is

$\displaystyle\mathcal{O}^{ee}_{LL}$
$\displaystyle=\bar{\tau}(1-\gamma_{5})\mu~{}\bar{e}(1-\gamma_{5})e,$
$\displaystyle\mathcal{O}^{ee}_{RR}$
$\displaystyle=\bar{\tau}(1+\gamma_{5})\mu~{}\bar{e}(1+\gamma_{5})e,$
$\displaystyle\mathcal{O}^{ee}_{LR}$
$\displaystyle=\bar{\tau}(1-\gamma_{5})\mu~{}\bar{e}(1+\gamma_{5})e,$
$\displaystyle\mathcal{O}^{ee}_{RL}$
$\displaystyle=\bar{\tau}(1+\gamma_{5})\mu~{}\bar{e}(1-\gamma_{5})e,$ (70)

with coefficients

$\displaystyle C^{ee}_{LL}$
$\displaystyle=C^{ee}_{RR}=\frac{y_{ee}y_{\mu\tau}}{4}\left[\frac{1}{M_{H}^{2}}-\frac{1}{M_{A}^{2}}\right],$
$\displaystyle\ C^{ee}_{LR}$
$\displaystyle=C^{ee}_{RL}=\frac{y_{ee}y_{\mu\tau}}{4}\left[\frac{1}{M_{H}^{2}}+\frac{1}{M_{A}^{2}}\right].$
(71)

Similar expressions follow for operators constructed from the set
$\bigl{\\{}(\bar{\tau}\mu)(\tau\bar{\tau}),(\bar{\tau}\gamma^{5}\mu)(\tau\gamma^{5}\bar{\tau})\bigl{\\}}$
with corresponding Wilson coefficients
$C^{\tau\tau}_{LL},C^{\tau\tau}_{RR},C^{\tau\tau}_{LR}$ and
$C^{\tau\tau}_{RL}$. After implementing these effective operators, we proceed
to compute the $\tau$-decays in flavio. In order to obtain a large
contribution to $\tau\to\mu\gamma$, we require that $C^{\tau\tau}_{RR}$ and
$C^{\tau\tau}_{LL}$ be non-zero, and hence a mass splitting between $H$ and
$A$ is needed.

Our theoretical setup has effects on the decay channel
$B_{s}\to\mu^{+}\mu^{-}$, which depends on the Wilson coefficients
$C^{\mu\mu}_{S}$ and $C^{\mu\mu}_{P}$ which are analogous to $C^{ee}_{S}$ and
$C^{ee}_{P}$ for $B_{s}\to e^{+}e^{-}$ with the leptonic coupling $y_{ee}$
replaced by $y_{\mu\mu}$, see Eq. (33). Since we have $y_{\mu\mu}\ll y_{ee}$
this implies that the NP contribution will be much smaller for muons than for
electrons. Nonetheless, we have checked that for each point in the allowed
parameter space the prediction for $\bar{\cal B}r(B_{s}\to\mu^{+}\mu^{-})$ is
in agreement with the experimental measurement given in Eq. (18) within
$2\sigma$.

Now, let us analyze the effects on $\bar{\mathcal{B}}r(B_{s}\rightarrow
e^{+}e^{-})$ and the lepton-flavour-violating decays. Firstly, we provide a
concrete example on how large the enhancement in ${\cal B}r(B_{s}\to
e^{+}e^{-})$ can be within the theory under consideration for concrete values
of the input parameters. Thus, we fix the mass of the scalar and pseudoscalar
to $M_{A}\\!=\\!800$ GeV, $M_{H}\\!=\\!400$ GeV and the rest of the parameters
to $\tan\beta=10$ and $s_{23}=s^{\prime}_{23}=0.98$. This implies an electron
coupling of $y_{ee}\simeq 0.13$ which is in agreement with the bound from LEP.
For the off-diagonal quark coupling we obtain $y_{bs}\simeq 4.2\times 10^{-4}$
which gives $\bar{\cal B}r(B_{s}\to e^{+}e^{-})\simeq 8.4\times 10^{-9}$. For
the choice of mixing angles discussed above, the off-diagonal lepton coupling
is $y_{\tau\mu}\simeq 1.1\times 10^{-3}$ and predicts ${\cal
B}r(\tau^{-}\to\mu^{-}\,e^{+}e^{-})\simeq 1.4\times 10^{-10}$ and ${\cal
B}r(\tau\to\mu\gamma)\simeq 6.6\times 10^{-13}$.

Finally, we can generalize the previous exercise while at the same time
assessing the impact on lepton-flavour-violating decays. Then, in the left
panel of Fig. 6 we show the predicted correlation between the observables
$\bar{\cal B}r(B_{s}\to e^{+}e^{-})$ and $\bar{\cal
B}r(\tau^{-}\to\mu^{-}\,e^{+}e^{-})$. The different bands correspond to
different values for the mass splitting $\Delta M\equiv M_{A}-M_{H}$. We fix
$s_{23}=s^{\prime}_{23}=0.98$ and $\tan\beta=10$. The region shaded in red
corresponds to the current experimental limits on these observables. In the
right panel in Fig. 6 we show the correlation between $\bar{\cal B}r(B_{s}\to
e^{+}e^{-})$ and ${\cal B}r(\tau\to\mu\gamma)$. From Fig. 6 we can see that it
is possible to saturate the current experimental bounds in $\bar{\cal
B}r(B_{s}\to e^{+}e^{-})$ while at the same time obeying the constraints on
the lepton-flavour-violating decays. This is the second result that we want to
highlight in this work. These plots provide a set of correlations between
different channels which make this framework phenomenologically testable.

##  6 Summary

The leptonic decay $B_{s}\to e^{+}e^{-}$ is a decay channel with interesting
properties and it can be used as smoking gun in the search for New Physics.
For instance, it is exceptionally clean. Moreover if this process takes place
as predicted by the Standard Model, due to the helicity suppression effect,
its tiny decay probability places it outside the reach of current or
forthcoming particle physics experiments. Therefore, any observation of this
channel in the near future would represent conclusive evidence for physics
beyond the Standard Model.

In this article and to the best of our knowledge, we have presented for the
first time, a concrete New Physics scenario which can provide a large
enhancement on the decay width for the channel $B_{s}\to e^{+}e^{-}$. More
specifically, by studying the general 2HDM in which both doublets are coupled
to the quarks and leptons of the Standard Model, we have demonstrated that
when the CP-odd scalar $A$ is mostly coupled to electrons, it can give a
contribution to the transition $B_{s}\to e^{+}e^{-}$ which enhances its decay
probability by up to five orders of magnitude above the Standard Model
prediction, saturating the most recent experimental upper bound established by
the LHCb collaboration. We have identified regions in the corresponding
parameter space where this potential enhancement respects all known
constraints from flavour and collider physics, including for instance neutral
$B_{s}$ mixing as well as LEP measurements of the $e^{-}e^{+}\to e^{-}e^{+}$
cross-section.

Furthermore, we have shown how the required coupling structure for the 2HDM
can arise from a UV theory of quark-lepton unification that can be realized at
a low energy scale. This framework predicts a correlation between the decay
channel $B_{s}\to e^{+}e^{-}$ and the lepton-flavour-violating decays
$\tau^{-}\to\mu^{-}e^{+}e^{-}$ and $\tau\to\mu\gamma$. We have worked out
quantitatively the interplay between these channels for different values of
the relevant free parameters. If the decay process $B_{s}\to e^{+}e^{-}$ is
observed in the near future, the presence of heavy (pseudo)scalars can be
further confirmed by searches for a heavy resonance decaying into an electron-
positron pair at the LHC. Our results show that the channel $B_{s}\to
e^{+}e^{-}$ is indeed a very interesting candidate to probe for New Physics
effects and provide additional justification to pursue further experimental
searches for it in the current and foreseeable experiments.

###### Acknowledgements.

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 945422. A.D.P. is supported by the INFN “Iniziativa Specifica”
Theoretical Astroparticle Physics (TAsP-LNF) and by the Frascati National
Laboratories (LNF) through a Cabibbo Fellowship, call 2020. M.B. and G.TX. are
supported by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
through grant 396021762 - TRR 257 “Particle Physics Phenomenology after the
Higgs Discovery”. Parts of the computations carried out for this work made use
of the OMNI cluster of the University of Siegen. We acknowledge useful
communication with Alexander Lenz and Matthew Kirk on the current updates of
$B_{s}\rightarrow\mu^{+}\mu-$.

## References

  * (1) Y. Aoki et al., _FLAG Review 2021_ , 2111.09849.
  * (2) A. Bazavov et al., _$B$ \- and $D$-meson leptonic decay constants from four-flavor lattice QCD_, _Phys. Rev. D_ 98 (2018) 074512, [1712.09262].
  * (3) ETM collaboration, A. Bussone et al., _Mass of the b quark and B -meson decay constants from N f=2+1+1 twisted-mass lattice QCD_, _Phys. Rev. D_ 93 (2016) 114505, [1603.04306].
  * (4) HPQCD collaboration, R. J. Dowdall, C. T. H. Davies, R. R. Horgan, C. J. Monahan and J. Shigemitsu, _B-Meson Decay Constants from Improved Lattice Nonrelativistic QCD with Physical u, d, s, and c Quarks_ , _Phys. Rev. Lett._ 110 (2013) 222003, [1302.2644].
  * (5) C. Hughes, C. T. H. Davies and C. J. Monahan, _New methods for B meson decay constants and form factors from lattice NRQCD_ , _Phys. Rev. D_ 97 (2018) 054509, [1711.09981].
  * (6) LHCb collaboration, R. Aaij et al., _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].
  * (7) LHCb collaboration, R. Aaij et al., _Analysis of Neutral B-Meson Decays into Two Muons_ , _Phys. Rev. Lett._ 128 (2022) 041801, [2108.09284].
  * (8) ATLAS collaboration, M. Aaboud et al., _Study of the rare decays of $B^{0}_{s}$ and $B^{0}$ mesons into muon pairs using data collected during 2015 and 2016 with the ATLAS detector_, _JHEP_ 04 (2019) 098, [1812.03017].
  * (9) 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}$, CMS-PAS-BPH-21-006_, .
  * (10) CDF collaboration, T. Aaltonen et al., _Search for the Decays $B^{0}_{s}\to e^{+}\mu^{-}$ and $B^{0}_{s}\to e^{+}e^{-}$ in CDF Run II_, _Phys. Rev. Lett._ 102 (2009) 201801, [0901.3803].
  * (11) LHCb collaboration, R. Aaij et al., _Search for the Rare Decays $B^{0}_{s}\to e^{+}e^{-}$ and $B^{0}\to e^{+}e^{-}$_, _Phys. Rev. Lett._ 124 (2020) 211802, [2003.03999].
  * (12) R. Fleischer, R. Jaarsma and G. Tetlalmatzi-Xolocotzi, _In Pursuit of New Physics with $B^{0}_{s,d}\to\ell^{+}\ell^{-}$_, _JHEP_ 05 (2017) 156, [1703.10160].
  * (13) J. F. Gunion, H. E. Haber, G. L. Kane and S. Dawson, _The Higgs Hunter’s Guide_ , vol. 80.  2000\.
  * (14) 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–102, [1106.0034].
  * (15) J. C. Pati and A. Salam, _Lepton Number as the Fourth Color_ , _Phys. Rev. D_ 10 (1974) 275–289.
  * (16) P. Minkowski, _$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?_, _Phys. Lett. B_ 67 (1977) 421–428.
  * (17) T. Yanagida, _Horizontal gauge symmetry and masses of neutrinos_ , _Conf. Proc. C_ 7902131 (1979) 95–99.
  * (18) M. Gell-Mann, P. Ramond and R. Slansky, _Complex Spinors and Unified Theories_ , _Conf. Proc. C_ 790927 (1979) 315–321, [1306.4669].
  * (19) R. N. Mohapatra and G. Senjanovic, _Neutrino Mass and Spontaneous Parity Nonconservation_ , _Phys. Rev. Lett._ 44 (1980) 912.
  * (20) P. Fileviez Perez and M. B. Wise, _Low Scale Quark-Lepton Unification_ , _Phys. Rev. D_ 88 (2013) 057703, [1307.6213].
  * (21) R. N. Mohapatra, _Mechanism for Understanding Small Neutrino Mass in Superstring Theories_ , _Phys. Rev. Lett._ 56 (1986) 561–563.
  * (22) R. N. Mohapatra and J. W. F. Valle, _Neutrino Mass and Baryon Number Nonconservation in Superstring Models_ , _Phys. Rev. D_ 34 (1986) 1642.
  * (23) 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].
  * (24) W. Altmannshofer and P. Stangl, _New physics in rare B decays after Moriond 2021_ , _Eur. Phys. J. C_ 81 (2021) 952, [2103.13370].
  * (25) LHCb collaboration, R. Aaij et al., _Measurement of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ branching fraction and effective lifetime and search for $B^{0}\to\mu^{+}\mu^{-}$ decays_, _Phys. Rev. Lett._ 118 (2017) 191801, [1703.05747].
  * (26) 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].
  * (27) A. Crivellin, D. Müller and C. Wiegand, _$b\to s\ell^{+}\ell^{-}$ transitions in two-Higgs-doublet models_, _JHEP_ 06 (2019) 119, [1903.10440].
  * (28) LHCb collaboration, R. Aaij et al., _Search for the decays $B_{s}^{0}\to\tau^{+}\tau^{-}$ and $B^{0}\to\tau^{+}\tau^{-}$_, _Phys. Rev. Lett._ 118 (2017) 251802, [1703.02508].
  * (29) D. M. Straub, _flavio: a Python package for flavour and precision phenomenology in the Standard Model and beyond_ , 1810.08132.
  * (30) MEG collaboration, A. M. Baldini et al., _Search for the lepton flavour violating decay $\mu^{+}\rightarrow\mathrm{e}^{+}\gamma$ with the full dataset of the MEG experiment_, _Eur. Phys. J. C_ 76 (2016) 434, [1605.05081].
  * (31) SINDRUM II collaboration, W. H. Bertl et al., _A Search for muon to electron conversion in muonic gold_ , _Eur. Phys. J. C_ 47 (2006) 337–346.
  * (32) SINDRUM collaboration, U. Bellgardt et al., _Search for the Decay mu+ — $>$ e+ e+ e-_, _Nucl. Phys. B_ 299 (1988) 1–6.
  * (33) J. F. Gunion and H. E. Haber, _The CP conserving two Higgs doublet model: The Approach to the decoupling limit_ , _Phys. Rev. D_ 67 (2003) 075019, [hep-ph/0207010].
  * (34) Particle Data Group collaboration, R. L. Workman and Others, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * (35) ALEPH, DELPHI, L3, OPAL, LEP Electroweak Working Group collaboration, J. Alcaraz et al., _A Combination of preliminary electroweak measurements and constraints on the standard model_ , hep-ex/0612034.
  * (36) L. Di Luzio, M. Kirk, A. Lenz and T. Rauh, _$\Delta M_{s}$ theory precision confronts flavour anomalies_, _JHEP_ 12 (2019) 009, [1909.11087].
  * (37) CKMfitter Group collaboration, J. Charles, A. Hocker, H. Lacker, S. Laplace, F. R. Le Diberder, J. Malcles et al., _updated results and plots available at:http://ckmfitter.in2p3.fr_, _Eur. Phys. J. C_ 41 (2005) 1–131, [hep-ph/0406184].
  * (38) HFLAV collaboration, Y. Amhis et al., _Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of 2021_, 2206.07501.
  * (39) LHCb collaboration, R. Aaij et al., _Test of lepton universality using $B^{+}\rightarrow K^{+}\ell^{+}\ell^{-}$ decays_, _Phys. Rev. Lett._ 113 (2014) 151601, [1406.6482].
  * (40) BELLE collaboration, S. Choudhury et al., _Test of lepton flavor universality and search for lepton flavor violation in $B\rightarrow K\ell\ell$ decays_, _JHEP_ 03 (2021) 105, [1908.01848].
  * (41) LHCb collaboration, R. Aaij et al., _Measurement of the $B^{0}\rightarrow K^{*0}e^{+}e^{-}$ branching fraction at low dilepton mass_, _JHEP_ 05 (2013) 159, [1304.3035].
  * (42) Belle collaboration, S. Wehle et al., _Lepton-Flavor-Dependent Angular Analysis of $B\to K^{\ast}\ell^{+}\ell^{-}$_, _Phys. Rev. Lett._ 118 (2017) 111801, [1612.05014].
  * (43) J. A. Bailey et al., _$B\to Kl^{+}l^{-}$ Decay Form Factors from Three-Flavor Lattice QCD_, _Phys. Rev. D_ 93 (2016) 025026, [1509.06235].
  * (44) R. R. Horgan, Z. Liu, S. Meinel and M. Wingate, _Rare $B$ decays using lattice QCD form factors_, _PoS_ LATTICE2014 (2015) 372, [1501.00367].
  * (45) A. Bharucha, D. M. Straub and R. Zwicky, _$B\to V\ell^{+}\ell^{-}$ in the Standard Model from light-cone sum rules_, _JHEP_ 08 (2016) 098, [1503.05534].
  * (46) N. Gubernari, A. Kokulu and D. van Dyk, _$B\to P$ and $B\to V$ Form Factors from $B$-Meson Light-Cone Sum Rules beyond Leading Twist_, _JHEP_ 01 (2019) 150, [1811.00983].
  * (47) P. Fileviez Perez, C. Murgui and A. D. Plascencia, _Leptoquarks and matter unification: Flavor anomalies and the muon g-2_ , _Phys. Rev. D_ 104 (2021) 035041, [2104.11229].
  * (48) P. Fileviez Perez and C. Murgui, _Flavor anomalies and quark-lepton unification_ , _Phys. Rev. D_ 106 (2022) 035033, [2203.07381].
  * (49) P. Fileviez Perez, E. Golias and A. D. Plascencia, _Two-Higgs-doublet model and quark-lepton unification_ , _JHEP_ 08 (2022) 293, [2205.02235].
  * (50) P. Fileviez Perez, E. Golias and A. D. Plascencia, _Probing quark-lepton unification with leptoquark and Higgs boson decays_ , _Phys. Rev. D_ 105 (2022) 075011, [2107.06895].
  * (51) BaBar collaboration, B. Aubert et al., _Searches for Lepton Flavor Violation in the Decays tau+- — $>$ e+- gamma and tau+- —$>$ mu+- gamma_, _Phys. Rev. Lett._ 104 (2010) 021802, [0908.2381].
  * (52) K. Hayasaka et al., _Search for Lepton Flavor Violating Tau Decays into Three Leptons with 719 Million Produced Tau+Tau- Pairs_ , _Phys. Lett. B_ 687 (2010) 139–143, [1001.3221].
  * (53) SuperB collaboration, B. O’Leary et al., _SuperB Progress Reports – Physics_ , 1008.1541.