# Combined explanations of $B$-physics anomalies, $(g-2)_{e,\mu}$ and neutrino
masses by scalar leptoquarks

Shao-Long Chen<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton
Physics (MoE) and Institute of Particle Physics, Central China Normal
University, Wuhan 430079, China  Center for High Energy Physics, Peking
University, Beijing 100871, China  Wen-wen Jiang<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton Physics (MoE) and Institute of Particle
Physics, Central China Normal University, Wuhan 430079, China  Ze-Kun Liu
<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton Physics (MoE)
and Institute of Particle Physics, Central China Normal University, Wuhan
430079, China

###### Abstract

We extend the contents of the standard model (SM) by introducing TeV-scale
scalar leptoquarks to generate neutrino masses and explain some current
observed deviations from the SM predictions, including the anomalous magnetic
moments of charged leptons (electron and muon) and $B$-physics anomalies
($R_{K^{(*)}}$ and $R_{D^{(*)}}$). The model consists of $\text{SU}(2)_{L}$
singlet leptoquark $S_{1}\sim(\bar{3},1,1/3)$, doublet leptoquark
$\tilde{R}_{2}\sim(3,2,1/6)$ and triplet leptoquark
$S_{3}\sim(\bar{3},3,1/3)$. We combine the constraints arising from the low-
energy lepton flavor violation, meson decay and mixing observables. We perform
a detailed phenomenological analysis and identify the minimized texture of
leptoquark Yukawa matrices to accommodate a unified explanation of the
anomalies and neutrino oscillation data.

##  I Introduction

The neutrino oscillation experiments have firmly established that neutrinos
are massive and have non-trivial mixing between different generations Super-
Kamiokande:1998kpq ; SNO:2002tuh ; KamLAND:2002uet ; T2K:2011ypd . The
experiments also indicate that the neutrino masses are much smaller than that
of charged fermions, which suggests that neutrinos may have specific sources
of mass generation. In the recent decades, a plethora of models have been
proposed to explain the neutrino mass and the natural way is the so called
seesaw mechanism Ma:1998dn . Type-I seasaw model Minkowski:1977sc ;
Mohapatra:1979ia ; Gell-Mann:1979vob ; Glashow:1979nm ; Schechter:1980gr
provides neutrino masses at the tree-level by extending the particle content
of the SM with three $\text{SU}(2)_{L}$-singlet right-handed neutrino fields,
while type-II Schechter:1980gr ; Mohapatra:1980yp ; Lazarides:1980nt  and
type-III Foot:1988aq  models introduce $\text{SU}(2)_{L}$-triplet scalar and
$\text{SU}(2)_{L}$-triplet fermions, respectively. Beyond tree level, the tiny
neutrino masses could radiatively originate from loop levels Zee:1980ai ;
Cheng:1980qt ; Zee:1985id ; Babu:1988ki ; Cai:2017jrq .

Extending the SM to include the source of the origin of neutrino mass and
mixing brings new physics, especially to the flavor sector. The intensity
frontier precision measurements may pin down the possible connections between
neutrino physics and flavor physics. Such as the anomalous magnetic moments of
electron and muon, there are long-standing discrepancies between the
theoretical predictions and measured values Aoyama:2017uqe ;
Keshavarzi:2018mgv ; Davier:2019can ; Aoyama:2020ynm ; Aoyama:2012wk ;
Aoyama:2019ryr ; Czarnecki:2002nt ; Gnendiger:2013pva ; Davier:2017zfy ;
Keshavarzi:2018mgv ; Colangelo:2018mtw ; Hoferichter:2019mqg ; Davier:2019can
; Keshavarzi:2019abf ; Kurz:2014wya ; Melnikov:2003xd ; Masjuan:2017tvw ;
Colangelo:2017fiz ; Hoferichter:2018kwz ; Gerardin:2019vio ; Bijnens:2019ghy ;
Colangelo:2019uex ; Blum:2019ugy ; Colangelo:2014qya ; Abi:2021gix . The
anomalies also include the ratios $R_{K^{(*)}}$ and $R_{D^{(*)}}$ in
$B$-decays, pointing towards the lepton flavor universality violation,
measured by BaBar BaBar:2012obs ; BaBar:2013mob , Belle Belle:2015qfa ;
Belle:2016dyj ; Belle:2016kgw  and LHCb LHCb:2017vlu ; LHCb:2017smo ;
LHCb:2019hip ; LHCb:2017avl ; LHCb:2021trn  collaborations. In this work, we
propose a model with scalar leptoquarks to provide a common explanation of
neutrino mass and these flavor anomalies.

Leptoquarks (LQs) have been introduced in many new physics models beyond the
SM and are very popular to explain $B$-physics anomalies with one or more
leptoquark states Dorsner:2016wpm ; Crivellin:2020mjs ; Carvunis:2021dss . The
unified solution to both $R_{K^{(*)}}$ and $R_{D^{(*)}}$ anomalies seems rule
out single scalar leptoquark models Angelescu:2018tyl . Among the scalar
leptoquarks, triplet $S_{3}\sim(\bar{3},3,1/3)$ can accommodate the
$R_{K^{(*)}}$ anomalies, while the $R_{D^{(*)}}$ anomalies can be resolved by
introducing either a singlet $S_{1}\sim(\bar{3},1,1/3)$ or a doublet
$R_{2}\sim(3,2,7/6)$ leptoquark. The double leptoquarks models were proposed
to explain both $R_{K^{(*)}}$ and $R_{D^{(*)}}$ anomalies, involving $S_{1}$
and $S_{3}$ combination Bigaran:2019bqv ; Saad:2020ihm ; Gherardi:2020qhc ;
Greljo:2021xmg ; Lee:2021jdr ; Bhaskar:2022vgk  or $R_{2}$ and $S_{3}$
combination Chen:2016dip ; Saad:2020ucl ; Babu:2020hun . Extending with
leptoquarks will give contribution to the anomalous magnetic moment of charged
lepton at one-loop level and the no-chiral scalar leptoquarks $S_{1}$ or
$R_{2}$, which have both left-chiral and right-chiral couplings, can provide
good explanations to the $a_{\mu}$ and $a_{e}$ deviations Dorsner:2020aaz ;
Bigaran:2020jil  simultaneously. The mixing between different type leptoquarks
can also generate non-trivial Majorana neutrino mass terms at one-loop level.
The minimal model to generate neutrino mass by the scalar leptoquark mixing
requires a pair of leptoquarks and the possible combinations are
$S_{1}-\tilde{R}_{2}(3,2,1/6)$, $S_{3}-\tilde{R}_{2}$ and $S_{3}-R_{2}$
AristizabalSierra:2007nf ; Dorsner:2017wwn ; Julio:2022ton ; Julio:2022bue ;
Chowdhury:2022dps . Motivated by the leptoquark abundant phenomenologies, we
attempt to extend the SM contents by scalar leptoquarks to generate neutrino
mass and explain the flavor anomalies mentioned above.

This paper is organized as follow: In Section II, we briefly introduce the
model set-up and the neutrino mass generation mechanism. In Section III, we
show how to explain the flavor anomalies in the model, including
$R_{K^{(*)}}$, $R_{D^{(*)}}$, $a_{\mu}$ and $a_{e}$. We discuss the
observables constraints on the leptoquark couplings in Section IV and then we
perform a detailed analysis of model parameter space and identify two
benchmark points in Section V and we conclude in the final section.

##  II The model and neutrino mass generation

###  II.1 The model

In addition to the SM fields, we introduce three scalar leptoquarks, including
an $\text{SU}(2)_{L}$ singlet $S_{1}\sim(\bar{3},1,1/3)$, a doublet
$\tilde{R}_{2}\sim(3,2,1/6)$ and a triplet $S_{3}\sim(\bar{3},3,1/3)$. The
scalar leptoquarks are denoted as

$\displaystyle
S_{1}(\bar{3},1,1/3)=S_{1}^{1/3}\,,\quad\tilde{R}_{2}(3,2,1/6)=(\tilde{R}_{2}^{2/3},\tilde{R}_{2}^{-1/3})^{\text{T}}\,,$
$\displaystyle
S_{3}(\bar{3},3,1/3)=\tau^{i}S_{3}^{i}=\begin{pmatrix}S_{3}^{1/3}&\sqrt{2}\,S_{3}^{4/3}\\\
\sqrt{2}\,S_{3}^{-2/3}&-S_{3}^{1/3}\end{pmatrix}\,,$ (1)

where $\tau^{i}$ ($i=1,2,3$) are the Pauli matrices and we define
$S_{3}^{4/3}=(S_{3}^{1}-\text{i}S_{3}^{2})/\sqrt{2}$,
$S_{3}^{-2/3}=(S_{3}^{1}+\text{i}S_{3}^{2})/\sqrt{2}$ and
$S_{3}^{1/3}=S_{3}^{3}$. The corresponding Yukawa terms that describe the
interactions between leptoquarks and fermions are given by

$\displaystyle\mathcal{L}_{Y}=$ $\displaystyle-
y_{1R}^{ij}\overline{u^{iC}_{R}}e^{j}_{R}S_{1}-y_{1L}^{ij}\overline{Q_{L}^{iC}}\,i\tau^{2}\,L_{L}^{j}S_{1}-y_{2L}^{ij}\
\overline{d^{i}_{R}}\tilde{R}_{2}^{T}\,i\tau^{2}\,L_{L}^{j}$ $\displaystyle-
y_{3L}^{ij}\ \overline{Q_{L}^{iC}}\,i\tau^{2}\,S_{3}L_{L}^{j}+\text{h.c.}\,,$
(2)

where $Q$ and $L$ denote the $\text{SU}(2)_{L}$ doublet left-handed quarks and
leptons, $u_{R}$, $d_{R}$ and $e_{R}$ denote the $\text{SU}(2)_{L}$ singlet
right-handed up-type quarks, down-type quarks and charged leptons,
respectively. All fields in Eq. (2) are represented in the flavor basis. For
phenomenological analysis, it is more convenient that we re-parametrize the
couplings in the fermion mass basis. The Yukawa coupling terms are then
rewritten in the mass basis of fermions as the following form,

$\displaystyle\mathcal{L}_{Y}=$ $\displaystyle-
y_{1R}^{ij}\overline{u^{iC}_{R}}e^{j}_{R}S_{1}^{1/3}+(V^{T}y_{1L})^{ij}\overline{d^{iC}_{L}}\nu_{L}^{j}S_{1}^{1/3}-y_{1L}^{ij}\overline{u^{iC}_{L}}e_{L}^{j}S_{1}^{1/3}+y_{2L}^{ij}\
\overline{d^{i}_{R}}\nu^{j}_{L}\tilde{R}_{2}^{-1/3}$ $\displaystyle-
y_{2L}^{ij}\
\overline{d^{i}_{R}}e_{L}^{j}\tilde{R}_{2}^{2/3}+(V^{T}y_{3L})^{ij}\
\overline{d^{iC}_{L}}\nu^{j}_{L}S_{3}^{1/3}+\sqrt{2}(V^{T}y_{3L})^{ij}\
\overline{d^{iC}_{L}}e^{j}_{L}S_{3}^{4/3}$ $\displaystyle-\sqrt{2}y_{3L}^{ij}\
\overline{u^{iC}_{L}}\nu^{j}_{L}S_{3}^{-2/3}+y_{3L}^{ij}\
\overline{u^{iC}_{L}}e^{j}_{L}S_{3}^{1/3}+\text{h.c.}\,.$ (3)

where $V$ is the CKM matrix. Since in our analysis of $(g-2)_{e,\mu}$ and
$B$-physics anomalies, the choice of neutrino mass or flavor basis has
negligible effect, the neutrino states in the above equation are kept in
flavor basis.

The renormalizable and gauge invariant scalar potential involving $H$,
$S_{1}$, $\tilde{R}_{2}$ and $S_{3}$ is described by

$\displaystyle V\supset\,$ $\displaystyle
m_{H}^{2}\,H^{\dagger}H+m_{1}^{2}\,S_{1}^{\dagger}S_{1}+m_{2}^{2}\,\tilde{R}_{2}^{\dagger}\tilde{R}_{2}+\frac{1}{2}m_{3}^{2}\,\text{Tr}\big{(}S_{3}^{\dagger}S_{3}\big{)}+\lambda_{H}\left(H^{\dagger}H\right)^{2}+\lambda_{1}\left(S_{1}^{\dagger}S_{1}\right)^{2}$
$\displaystyle+\lambda_{2}\left(\tilde{R}_{2}^{\dagger}\tilde{R}_{2}\right)^{2}+\lambda_{3}\left[\text{Tr}(S_{3}^{\dagger}S_{3})\right]^{2}+\lambda_{3}^{\prime}\,\text{Tr}(S_{3}^{\dagger}S_{3}^{\dagger})\text{Tr}(S_{3}S_{3})+\lambda_{H1}\,H^{\dagger}HS_{1}^{\dagger}S_{1}$
$\displaystyle+\lambda_{H2}\,H^{\dagger}H\tilde{R}_{2}^{\dagger}\tilde{R}_{2}+\frac{1}{2}\lambda_{H3}\,H^{\dagger}H\,\text{Tr}\big{(}S_{3}^{\dagger}S_{3}\big{)}+\big{(}\lambda_{13}\,H^{\dagger}S_{3}^{\dagger}HS_{1}+\mu_{1}\,\tilde{R}_{2}^{\dagger}HS_{1}^{*}$
$\displaystyle+\mu_{2}\,\tilde{R}_{2}^{\dagger}S_{3}^{\dagger}H+\text{h.c.}\big{)}\,,$
(4)

where $H$ is the SM Higgs doublet. More general interactions of leptoquarks
and SM Higgs can be found in Ref. Dorsner:2022twk . After the spontaneous
electroweak symmetry breaking, the Higgs field $H$ acquires a vacuum expecting
value (VEV) with $\langle H\rangle=v/\sqrt{2},v=246\,\text{GeV}$. The physical
scalar particles include one electric neutral Higgs boson $h$, three
$1/3$-charged leptoquarks, two $2/3$-charged leptoquarks, and one
$4/3$-charged leptoquark. In the basis of
$\rho^{1/3}\equiv(S_{1}^{1/3},\tilde{R}_{2}^{1/3},S_{3}^{1/3})^{T}$ and
$\rho^{2/3}\equiv(\tilde{R}_{2}^{2/3},S_{3}^{2/3})^{T}$, the mass matrices for
the two groups of charged scalar particles are given by

$\displaystyle
M_{1/3}^{2}=\begin{pmatrix}m_{1}^{2}+\frac{1}{2}\lambda_{H1}v^{2}&\frac{1}{\sqrt{2}}\mu_{1}v&-\frac{1}{2}\lambda_{13}v^{2}\\\\[5.69046pt]
\frac{1}{\sqrt{2}}\mu_{1}v&m_{2}^{2}+\frac{1}{2}\lambda_{H2}v^{2}&-\frac{1}{\sqrt{2}}\mu_{2}v\\\\[5.69046pt]
-\frac{1}{2}\lambda_{13}v^{2}&-\frac{1}{\sqrt{2}}\mu_{2}v&m_{3}^{2}+\frac{1}{2}\lambda_{H3}v^{2}\end{pmatrix}\,,$
(5) $\displaystyle
M_{2/3}^{2}=\begin{pmatrix}m_{2}^{2}+\frac{1}{2}\lambda_{H2}v^{2}&\mu_{2}v\\\\[5.69046pt]
\mu_{2}v&m_{3}^{2}+\frac{1}{2}\lambda_{H3}v^{2}\end{pmatrix}\,.$ (6)

After diagonalization of the above mass matrices, we obtain the physical
scalar fields: charge-1/3 leptoquarks $(\phi_{1},\phi_{2},\phi_{3})$ and
charge-2/3 leptoquarks $(\omega_{1},\omega_{2})$, which satisfy

$\displaystyle\phi_{i}=R^{1/3}_{ij}\,\rho^{1/3}_{j}\,,$ (7)
$\displaystyle\omega_{i}=R^{2/3}_{ij}\,\rho^{2/3}_{j}\,,$ (8)

where $R^{1/3}$ and $R^{2/3}$ are the corresponding rotation matrices. The
rotation matrix $R^{2/3}$ can be parametrized as

$R^{2/3}=\begin{pmatrix}\cos\alpha&\sin\alpha\\\
-\sin\alpha&\cos\alpha\end{pmatrix}\,,$ (9)

where the mixing angle is given by

$\tan
2\alpha=\frac{2\mu_{2}v}{m_{2}^{2}-m_{3}^{2}+(\lambda_{H2}-\lambda_{H3})v^{2}/2}\,.$
(10)

The rotation matrix $R^{1/3}$ need three rotation angles to be parametrized,

$R^{1/3}=R(\theta_{12})R(\theta_{13})R(\theta_{23})\,.$ (11)

In the limit where off-diagonal elements are much smaller than the diagonal
elements, the mixing angle in the rotation matrix $R^{1/3}$ can be
approximatively calculated by

$\theta_{ij}\simeq\frac{\left(M_{1/3}^{2}\right)_{ij}}{\left(M_{1/3}^{2}\right)_{ii}-\left(M_{1/3}^{2}\right)_{jj}}\,.$
(12)

The charge-4/3 component $S_{3}^{4/3}$ has no mixing with other scalar fields
and we denote the mass by $m_{S_{3}}$. In our analysis of the low energy
processes, we assume the leptoquark multiplets to be quasi-degenerate and set
the LQ masses as $m_{S_{1}}=m_{\phi_{1}}$, $m_{R_{2}}=m_{\phi_{2}}\approx
m_{\omega_{1}}$ and $m_{S_{3}}=m_{\phi_{3}}\approx m_{\omega_{2}}$.

###  II.2 Neutrino masses

$\times$$\nu_{L}$$d$$d^{c}$$\nu_{L}$$\phi^{1/3}_{i}$$\phi^{-1/3}_{i}$$\langle
H\rangle$ Figure 1: Feynman diagram of Majorana neutrino masses generation at
one-loop level.

In our model, the neutrino masses are induced at one-loop level through the
Feynman diagrams as shown in Fig. 1, in which the loop is mediated by the
down-type quarks and 1/3-charged leptoquarks. The neutrino mass matrix is
given by AristizabalSierra:2007nf

$\displaystyle({\cal M}_{\nu})_{\alpha\beta}=$
$\displaystyle\frac{3}{16\pi^{2}}\sum_{\begin{subarray}{c}i=1,2,3\\\
k=d,s,b\end{subarray}}m_{k}B_{0}(0,m_{k}^{2},m_{\phi_{i}}^{2})\left\\{R_{i1}^{1/3}R_{i2}^{1/3}\left[(V^{T}y_{1L})^{k\alpha}y_{2L}^{k\beta}+(V^{T}y_{1L})^{k\beta}y_{2L}^{k\alpha}\right]\right.$
$\displaystyle\left.+R_{i2}^{1/3}R_{i3}^{1/3}\left[(V^{T}y_{3L})^{k\alpha}y_{2L}^{k\beta}+(V^{T}y_{3L})^{k\beta}y_{2L}^{k\alpha}\right]\right\\}\,,$
(13)

where $B_{0}(0,m_{k}^{2},m_{\phi_{i}}^{2})$ is the Passarino-Veltman function
and its finite part is given by

$B_{0}(0,m_{k}^{2},m_{\phi_{i}}^{2})=\frac{m_{k}^{2}\log(m_{k}^{2})-m_{\phi_{i}}^{2}\log(m_{\phi_{i}}^{2})}{m_{k}^{2}-m_{\phi_{i}}^{2}}\,.$
(14)

The first term in the bracket of Eq. (13) is associated with the
$S_{1}-\tilde{R}_{2}$ combination, while the second term is associated with
the $S_{3}-\tilde{R}_{2}$ combination. To simplify the analysis, we consider
one term dominates the other. For example, when $S_{1}-\tilde{R}_{2}$
contribution is dominant ($\mu_{1}\gg\mu_{2},\lambda_{13}v$), the neutrino
mass matrix can be written as

$({\cal M}_{\nu})_{\alpha\beta}=\big{(}\hat{y}_{1L}^{T}\Lambda
y_{2L}+y_{2L}^{T}\Lambda^{T}\hat{y}_{1L}\big{)}_{\alpha\beta}\,,$ (15)

where we define $\hat{y}_{1L}\equiv(V^{T}y_{1L})$ and

$\Lambda\equiv\begin{pmatrix}\Lambda_{d}&0&0\\\ 0&\Lambda_{s}&0\\\
0&0&\Lambda_{b}\\\ \end{pmatrix}\,,\qquad\text{with}\
\Lambda_{k}\simeq\frac{3}{32\pi^{2}}m_{k}\frac{\sqrt{2}\mu_{1}v}{m_{\phi_{1}}^{2}-m_{\phi_{2}}^{2}}\log\left(\frac{m_{\phi_{1}}^{2}}{m_{\phi_{2}}^{2}}\right)\,.$
(16)

Using the method of master parametrization Cordero-Carrion:2018xre ; Cordero-
Carrion:2019qtu , we parametrize the coupling matrices $\hat{y}_{1L}$ and
$y_{2L}$ as

$\displaystyle\hat{y}_{1L}=\frac{1}{\sqrt{2}}\Sigma^{-1/2}WA\hat{D}^{1/2}U^{\dagger}\,,$
(17) $\displaystyle
y_{2L}=\frac{1}{\sqrt{2}}\Sigma^{-1/2}W^{*}B\hat{D}^{1/2}U^{\dagger}\,,$ (18)

where $U$ is the $3\times 3$ unitary neutrino mixing matrix, which brings the
neutrino mass matrix to diagonal form by

$U^{T}\,\mathcal{M}_{\nu}\,U=\text{diag}(m_{1},m_{2},m_{3})\,.$ (19)

The forms of matrix $\Sigma,W,A,B$ and $\hat{D}$ depend on the ranks of
neutrino mass matrix $\mathcal{M}_{\nu}$ and matrix $\Lambda$. Neutrino
oscillation data requires that $\cal{M}_{\nu}$ should contain two or three
non-vanishing eigenvalues. In our numerical analysis of neutrino masses, for
simplicity, we neglect the $d-$quark contribution in neutrino mass loop
$(\Lambda_{d}=0)$ and consider the normal ordering neutrino mass hierarchy
with $m_{1}=0$. In this scenario, the ranks of matrices $\mathcal{M}_{\nu}$
and $\Lambda$ are both 2 and the neutrino mass matrix $\mathcal{M}_{\nu}$ only
depends on the second and third columns of couplings $\hat{y}_{1L}$ and
$y_{2L}$. In this case, $\Sigma$ takes form as
$\text{diag}(\Lambda_{s},\Lambda_{b})$ and $\hat{D}$ takes the form as
$\text{diag}(\kappa,m_{2},m_{3})$, where $\kappa$ can be arbitrary value,
since it can always be absorbed by rescaling relevant elements in matrices A
and B. Eq. (17) and Eq. (18) give the elements of second and third columns of
couplings $\hat{y}_{1L}$ and $y_{2L}$. The matrix $W$ is a $2\times 2$ unitary
complex matrix that contains 4 real degrees of freedom. The matrices $A$ and
$B$ are defined as $A=TC_{1}$ and $B=(T^{T})^{-1}(C_{1}C_{2}+KC_{2})$, where
$T$ is an upper-triangular $2\times 2$ complex matrix with positive real
values in the diagonal and contains 4 degrees of freedom, $K$ is a $2\times 2$
anti-symmetric complex matrix that contains 2 degrees of freedom. The matrices
$C_{1}$ and $C_{2}$ are given by

$C_{1}=\begin{pmatrix}z_{1}\,&0\ &0\ \\\ z_{2}\,&0\ &0\ \end{pmatrix}\,,\ \ \
\ C_{2}=\begin{pmatrix}-1\ &0\ &0\ \\\ 0\ &1\ &0\ \\\ 0\ &0\ &1\
\end{pmatrix}\,,$ (20)

where $z_{1}$ and $z_{2}$ are two complex numbers that contains 2 degrees of
freedom with the condition $z_{1}^{2}+z_{2}^{2}=0$. The possible values of the
second and third columns of matrices $\hat{y}_{1L}$ and $y_{2L}$ can be
obtained by scanning these 12 real free parameters.

###  II.3 Effective Lagrangians

The tree-level contributions of leptoquarks to the related phenomenologies can
be described by the following effective Lagrangians,

$\displaystyle\mathcal{L}_{\bar{q}q\bar{\ell}\ell}=-\frac{4G_{F}}{\sqrt{2}}\Big{[}$
$\displaystyle\left(g^{LL}_{V,q}\right)^{ij,mn}(\bar{q}_{L}^{i}\gamma^{\mu}q_{L}^{j})(\bar{\ell}_{L}^{m}\gamma_{\mu}\ell_{L}^{n})+\left(g^{RL}_{V,q}\right)^{ij,mn}(\bar{q}_{R}^{i}\gamma^{\mu}q_{R}^{j})(\bar{\ell}_{L}^{m}\gamma_{\mu}\ell_{L}^{n})$
$\displaystyle+\left(g^{LR}_{V,q}\right)^{ij,mn}(\bar{q}_{L}^{i}\gamma^{\mu}q_{L}^{j})(\bar{\ell}_{R}^{m}\gamma_{\mu}\ell_{R}^{n})+\left(g^{RR}_{V,q}\right)^{ij,mn}(\bar{q}_{R}^{i}\gamma^{\mu}q_{R}^{j})(\bar{\ell}_{R}^{m}\gamma_{\mu}\ell_{R}^{n})$
$\displaystyle+\left(g^{LL}_{S,q}\right)^{ij,mn}(\bar{q}_{R}^{i}q_{L}^{j})(\bar{\ell}_{R}^{m}\ell_{L}^{n})+\left(g^{RR}_{S,q}\right)^{ij,mn}(\bar{q}_{L}^{i}q_{R}^{j})(\bar{\ell}_{L}^{m}\ell_{R}^{n})$
$\displaystyle+\left(g^{LL}_{T,q}\right)^{ij,mn}(\bar{q}_{R}^{i}\sigma^{\mu\nu}q_{L}^{j})(\bar{\ell}_{R}^{m}\sigma_{\mu\nu}\ell_{L}^{n})+\left(g^{RR}_{T,q}\right)^{ij,mn}(\bar{q}_{L}^{i}\sigma^{\mu\nu}q_{R}^{j})(\bar{\ell}_{L}^{m}\sigma_{\mu\nu}\ell_{R}^{n})\Big{]}\,,$
(21)
$\displaystyle\mathcal{L}_{\bar{q}q\bar{\nu}\nu}=\frac{4G_{F}}{\sqrt{2}}\Big{[}$
$\displaystyle\left(h^{LL}_{V,q}\right)^{ij,mn}(\bar{q}_{L}^{i}\gamma^{\mu}q_{L}^{j})(\bar{\nu}_{L}^{m}\gamma_{\mu}\nu_{L}^{n})+\left(h^{RL}_{V,q}\right)^{ij,mn}(\bar{q}_{R}^{i}\gamma^{\mu}q_{R}^{j})(\bar{\nu}_{L}^{m}\gamma_{\mu}\nu_{L}^{n})\Big{]}\,,$
(22)
$\displaystyle\mathcal{L}_{\bar{u}d\bar{\ell}\nu}=-\frac{4G_{F}}{\sqrt{2}}\Big{[}$
$\displaystyle\left(c^{LL}_{V}\right)^{ij,mn}(\bar{u}_{L}^{i}\gamma^{\mu}d_{L}^{j})(\bar{\ell}_{L}^{m}\gamma_{\mu}\nu_{L}^{n})+\left(c^{LL}_{S}\right)^{ij,mn}(\bar{u}_{R}^{i}d_{L}^{j})(\bar{\ell}_{R}^{m}\nu_{L}^{n})$
$\displaystyle+\left(c^{LL}_{T}\right)^{ij,mn}(\bar{u}_{R}^{i}\sigma^{\mu\nu}d_{L}^{j})(\bar{\ell}_{R}^{m}\sigma_{\mu\nu}\nu^{n}_{L})\Big{]}+\text{h.c.}\,.$
(23)

The Wilson coefficients at the leptoquark mass scale are determined by the
combinations of Yukawa couplings and summarized in Table 1. To analyze the
low-energy processes, these Wilson coefficients are needed to RGE run down to
the appropriate scale. We take the low-energy scale at the bottom-quark mass
($m_{b}=4.18\ \text{GeV}$) and the Wilson coefficients at the leading
logarithm approximation can be calculated by the following form
Chetyrkin:1997dh ; Gracey:2000am ; Dorsner:2013tla ; Hiller:2016kry ,

$C_{J}(\mu=m_{b})=\left[\frac{\alpha_{s}(m_{b})}{\alpha_{s}(m_{t})}\right]^{-\gamma_{1}^{J}/\beta^{(5)}_{1}}\left[\frac{\alpha_{s}(m_{t})}{\alpha_{s}(m_{\text{LQ}})}\right]^{-\gamma_{1}^{J}/\beta^{(6)}_{1}}C_{J}(\mu=m_{\text{LQ}})\,,$
(24)

with the QCD running coefficient $\beta^{(n_{f})}_{1}=(2n_{f}-33)/6$, where
$n_{f}$ is the relevant number of quark flavors at the hadronic scale. The
coefficients $\gamma_{1}^{J}$ are the anomalous dimension and given by
$\gamma_{1}^{V}=0$, $\gamma_{1}^{S}=2$ and $\gamma_{1}^{T}=-2/3$. In our
numerical analysis, we use the package Wilson Aebischer:2018bkb  to calculate
the running of Wilson coefficients and obtain the following Wilson coefficient
correlations between the two scales,

$\displaystyle c_{V}^{LL}(m_{b})$
$\displaystyle=1.01\,c_{V}^{LL}(m_{\text{LQ}})\,,$
$\displaystyle\begin{pmatrix}c_{S}^{LL}(m_{b})\\\
c_{T}^{LL}(m_{b})\end{pmatrix}$ $\displaystyle=\begin{pmatrix}1.64&-0.275\\\
-3.87\times
10^{-3}&0.867\end{pmatrix}\begin{pmatrix}c_{S}^{LL}(m_{\text{LQ}})\\\
c_{T}^{LL}(m_{\text{LQ}})\end{pmatrix}$ (25)

where we have taken the leptoquark masses scale as $m_{\text{LQ}}=1\
\text{TeV}$. In the following discussion of the various physical processes, we
utilize the Flavio package Straub:2018kue  to get the favored region of Wilson
coefficients at the leptoquark mass scale.

| $S_{1}$ | $\tilde{R}_{2}$ | $S_{3}$
---|---|---|---
$u^{i}\to u^{j}\bar{\ell}^{m}\ell^{n}$ | $g^{LL}_{V,u}=-(y_{1L})^{im}(y_{1L}^{*})^{jn}$ |  | $g^{LL}_{V,u}=-(y_{3L})^{im}(y_{3L}^{*})^{jn}$
$g^{RR}_{V,u}=-(y_{1R})^{im}(y_{1R}^{*})^{jn}$ |  |
$g^{LL}_{S,u}=(y_{1R})^{im}(y_{1L}^{*})^{jn}$ |  |
$g^{RR}_{S,u}=(y_{1L})^{im}(y_{1R}^{*})^{jn}$ |  |
$g^{LL}_{T,u}=-\frac{1}{4}(y_{1R})^{im}(y_{1L}^{*})^{jn}$ |  |
$g^{RR}_{T,u}=-\frac{1}{4}(y_{1L})^{im}(y_{1R}^{*})^{jn}$ |  |
$d^{i}\to d^{j}\bar{\ell}^{m}\ell^{n}$ |  | $g^{RL}_{V,d}=(y_{2L})^{jm}(y_{2L}^{*})^{in}$ | $g^{LL}_{V,d}=-2(V^{T}y_{3L})^{im}(V^{\dagger}y_{3L}^{*})^{jn}$
$u^{i}\to u^{j}\bar{\nu}^{m}\nu^{n}$ |  |  | $h^{LL}_{V,u}=2(y_{3L})^{im}(y_{3L}^{*})^{jn}$
$d^{i}\to d^{j}\bar{\nu}^{m}\nu^{n}$ | $h^{LL}_{V,d}=-(V^{T}y_{1L})^{im}(V^{\dagger}y_{1L}^{*})^{jn}$ | $h^{RL}_{V,d}=-(y_{2L})^{jm}(y_{2L}^{*})^{in}$ | $h^{LL}_{V,d}=(V^{T}y_{3L})^{im}(V^{\dagger}y_{3L}^{*})^{jn}$
$d^{i}\to u^{j}\ell^{n}\bar{\nu}^{m}$ | $c^{LL}_{V}=(V^{T}y_{1L})^{im}(y_{1L}^{*})^{jn}$ |  | $c^{LL}_{V}=-(V^{T}y_{3L})^{im}(y_{3L}^{*})^{jn}$
$c^{LL}_{S}=-(V^{T}y_{1L})^{im}(y_{1R}^{*})^{jn}$ |  |
$c^{LL}_{T}=\frac{1}{4}(V^{T}y_{1L})^{im}(y_{1R}^{*})^{jn}$ |  |
Table 1: The corresponding Wilson coefficients (in units of
$v^{2}/4\,m_{\text{LQ}}^{2}$) in Eqs. (21, 22, 23) induced by the leptoquarks
at tree level. The matching scale is set at the leptoquark mass scale.

##  III The flavor anomalies

In this section, we present the observed flavor anomalies between current
experimental observations and the SM predictions, and explore how to alleviate
these tensions by introducing scalar leptoquarks in our framework.

###  III.1 $R_{K}$ and $R_{K^{*}}$

The first observed anomalies we consider are the lepton flavor universality
violation ratios $R_{K}$ and $R_{K^{*}}$, which are defined as

$R_{K}=\frac{\text{Br}(B^{+}\to K^{+}\mu^{+}\mu^{-})}{\text{Br}(B^{+}\to
K^{+}e^{+}e^{-})}\,,\quad R_{K^{*}}=\frac{\text{Br}(B^{0}\to
K^{*0}\mu^{+}\mu^{-})}{\text{Br}(B^{0}\to K^{*0}e^{+}e^{-})}\,.$ (26)

The SM predictions Bobeth:2007dw ; Bordone:2016gaq  for these two ratios are

$R_{K}^{\rm{SM}}=1.0003\pm 0.0001\,,\quad R_{K^{*}}^{\rm{SM}}=1.00\pm 0.01\,.$
(27)

The new measurements of $R_{K}$ and $R_{K^{*}}$ at low $q^{2}$ region
$[1.1,\,6.0]\ \rm{GeV^{2}}$ by LHCb are given by LHCb:2019hip ; LHCb:2021trn

$R_{K}^{\rm{LHCb}}=0.846^{+0.042\,+0.013}_{-0.039\,-0.012}\,,\quad
R_{K^{*}}^{\rm{LHCb}}=0.685^{+0.113}_{-0.069}\pm 0.047\,,$ (28)

which both give deviation larger than $2.5\,\sigma$ from the SM prediction
values. These two processes are determined by the neutral current, $b\to
s\ell^{+}\ell^{-}$. The effective Hamiltonian relevant to our model can be
described by Buchalla:1995vs ,

$\mathcal{H}_{\text{eff}}=-\frac{4G_{F}}{\sqrt{2}}V_{tb}V_{ts}^{*}\Big{[}\sum_{X=9,10}(C_{X}^{\ell\ell}\mathcal{O}_{X}^{\ell\ell}+C_{X^{\prime}}^{\ell\ell}\mathcal{O}_{X^{\prime}}^{\ell\ell})\Big{]}+\text{h.c.}\,,$
(29)

where the $C_{X}^{\ell\ell}$ and $C_{X^{\prime}}^{\ell\ell}$ denote the Wilson
coefficients and $\mathcal{O}_{X}^{\ell\ell}$ and
$\mathcal{O}_{X^{\prime}}^{\ell\ell}$ are the corresponding effective
operators, which take form as

$\displaystyle\mathcal{O}_{9}^{\ell\ell}=\frac{e^{2}}{(4\pi)^{2}}(\bar{s}\gamma^{\mu}P_{L}b)(\bar{\ell}\gamma_{\mu}\ell)\,,\quad\mathcal{O}_{10}^{\ell\ell}=\frac{e^{2}}{(4\pi)^{2}}(\bar{s}\gamma^{\mu}P_{L}b)(\bar{\ell}\gamma_{\mu}\gamma_{5}\ell)\,,$
$\displaystyle\mathcal{O}_{9^{\prime}}^{\ell\ell}=\frac{e^{2}}{(4\pi)^{2}}(\bar{s}\gamma^{\mu}P_{R}b)(\bar{\ell}\gamma_{\mu}\ell)\,,\quad\mathcal{O}_{10^{\prime}}^{\ell\ell}=\frac{e^{2}}{(4\pi)^{2}}(\bar{s}\gamma^{\mu}P_{R}b)(\bar{\ell}\gamma_{\mu}\gamma_{5}\ell)\,.$
(30)

Figure 2: Contour plot of the fit to $R_{K}$ and $R_{K^{*}}$ in the plane of
Wilson coefficients at the scale $\mu=1\,\text{TeV}$. Upper left panel
corresponds to $C_{9}^{\mu\mu}$ versus $C_{10}^{\mu\mu}$ as real. The dashed
orange and dotted blue contours represent the $1\sigma$ allowed regions that
explain $R_{K}$ and $R_{K^{*}}$ respectively. The dark green (yellow) region
corresponds to the $1\sigma$ allowed region that explains $R_{K}$ and
$R_{K^{*}}$ ($R_{K^{0}_{S}}$ and $R_{K^{*+}}$) simultaneously, while lighter
region corresponds to $2\sigma$ allowed region. Upper right: The plot
corresponds to the complex plane of $C_{9,10}^{\mu\mu}$, with
$C_{9}^{\mu\mu}=-C_{10}^{\mu\mu}$ assumed. The orange and blue regions
represent the $1\sigma$ allowed regions that explain $R_{K}$ and $R_{K^{*}}$,
respectively. Bottom: Left panel shows the fit to $R_{K^{(*)}}$ using
$C_{9^{\prime}}^{\mu\mu}$ and $C_{10^{\prime}}^{\mu\mu}$ as real parameters,
while right one corresponds to the fit using complex parameters with the
assumption $C_{9^{\prime}}^{\mu\mu}=-C_{10^{\prime}}^{\mu\mu}$. No overlap
region indicates that $C_{9^{\prime}}^{\mu\mu}$ and $C_{10^{\prime}}^{\mu\mu}$
can not accommodate combined explanation for $R_{K}$ and $R_{K^{*}}$.

According to the definition of $R_{K}$ and $R_{K^{*}}$, the anomalies of
$R_{K^{(*)}}$ indicate new physics contribution to $C_{9^{(\prime)}}^{ee}$,
$C_{10^{(\prime)}}^{ee}$, $C_{9^{(\prime)}}^{\mu\mu}$ and
$C_{10^{(\prime)}}^{\mu\mu}$. The solution of $R_{K^{*}}$ is favored by new
physics coupling to muon instead of electron, with the consideration from
other observables fit Capdevila:2017bsm ; Hurth:2017hxg ; Aebischer:2019mlg .
Therefore, we set the new physics contribution related to electron is
negligible (i.e., $C_{9}^{ee}=C_{10}^{ee}\sim 0$), and the new physics
contributions to $R_{K}$ and $R_{K^{*}}$ come from the $C_{9}^{\mu\mu}$ and
$C_{10}^{\mu\mu}$ in our framework. We show the fit to $R_{K^{(*)}}$ using
$C_{9^{(\prime)},10^{(\prime)}}$ at the scale $\mu=1\ \text{TeV}$ in Fig. 2.
The upper left panel gives favored regions for $C_{9}^{\mu\mu}$ versus
$C_{10}^{\mu\mu}$ as real parameters. We also consider the recent measurements
for the ratio $R_{K^{0}_{S}}$ and $R_{K^{*+}}$ LHCb:2021lvy . Note that we
combine the constraint from $B_{s}\to\mu\mu$ in the fit. The branching ratio
of $B_{s}\to\mu\mu$ is measured to be
$\text{Br}(B_{s}\to\mu\mu)^{\text{exp}}=(2.93\pm 0.35)\times 10^{-9}$
Geng:2021nhg , which is the combined result based on measurements from ATLAS,
CMS and LHCb ATLAS:2018cur ; CMS:2019bbr ; LHCb:2021awg ; LHCb:2021vsc , while
the SM prediction value is $\text{Br}(B_{s}\to\mu\mu)^{\text{SM}}=(3.63\pm
0.13)\times 10^{-9}$ Beneke:2019slt . Taking the relation
$C_{9}^{\mu\mu}=-C_{10}^{\mu\mu}$ given by our model, the best fit point of
$R_{K^{(*)}}$ is found at $C_{9}^{\mu\mu}=-C_{10}^{\mu\mu}=-0.39$, while the
best fit point for $R_{K^{0}_{S}}$ and $R_{K^{*+}}$ is found at
$C_{9}^{\mu\mu}=-C_{10}^{\mu\mu}=-0.74$. Combining these four experimental
ratios and the branching ratio of $B_{s}\to\mu\mu$, the best fit piont of
$C_{9}^{\mu\mu}=-C_{10}^{\mu\mu}$ is $-0.45$. The upper right panel shows the
favored region for the complex case with the assumption
$C_{9}^{\mu\mu}=-C_{10}^{\mu\mu}$. The bottom panels present the fit to
$R_{K}$ and $R_{K^{*}}$ using $C_{9^{\prime},10^{\prime}}^{\mu\mu}$ and we
find no common solution. Our results are comparable with the global analysis
performed in Ref. Alok:2022pjb , where some related differential branching
ratios and angular observables are included. Relevant analyses are also found
in Ref. Becirevic:2015asa ; Carvunis:2021jga ; Altmannshofer:2021qrr .

Leptoquark $S_{1}$ doesn’t contribute to $b\to s\ell^{+}\ell^{-}$ at tree-
level but provides contribution by box-diagrams. However, $R_{K}$ and
$R_{K^{*}}$ anomalies cannot be fully accommodated with leptoquark $S_{1}$
only Becirevic:2016oho ; Angelescu:2018tyl . In our model, we expect that the
contributions to solve $R_{K}$ and $R_{K^{*}}$ anomalies come dominantly from
leptoquark $S_{3}$. The corresponding Wilson coefficients are given by

$C_{9}^{\ell\ell}=-C_{10}^{\ell\ell}=\frac{\pi
v^{2}}{V_{tb}V_{ts}^{*}\alpha_{\text{em}}}\frac{(V^{T}y_{3L})^{3\ell}(V^{\dagger}y_{3L}^{*})^{2\ell}}{m_{S_{3}}^{2}}\,.$
(31)

Leptoquark $\tilde{R}_{2}$ can also generate contribution to the process $b\to
s\ell^{+}\ell^{-}$ at tree-level by $C_{9^{\prime}}$ and $C_{10^{\prime}}$
terms. The corresponding Wilson coefficients of $\tilde{R}_{2}$ contribution
are given by

$C_{9^{\prime}}^{\ell\ell}=-C_{10^{\prime}}^{\ell\ell}=-\frac{\pi
v^{2}}{2V_{tb}V_{ts}^{*}\alpha_{\text{em}}}\frac{y_{2L}^{2\ell}y_{2L}^{*3\ell}}{m_{R_{2}}^{2}}\,.$
(32)

It is noted that the parameter space to explain $R_{K}$ is incompatible with
$R_{K^{*}}$ if one only use $C_{9^{\prime}}$ and $C_{10^{\prime}}$, as shown
in the bottom panel of Fig. 2.

###  III.2 $R_{D}$ and $R_{D^{*}}$

The next lepton flavor universality violation observables we consider are
$R_{D}$ and $R_{D^{*}}$, which are induced by charged current transitions
$b\to c\ell\bar{\nu}_{\ell}$ and defined as

$R_{D}=\frac{\text{Br}(B\to D\tau\bar{\nu})}{\text{Br}(B\to
D\ell\bar{\nu})}\,,\qquad R_{D^{*}}=\frac{\text{Br}(B\to
D^{*}\tau\bar{\nu})}{\text{Br}(B\to D^{*}\ell\bar{\nu})}\,,$ (33)

where $\ell$ denotes electron $e$ or muon $\mu$. The predicted values of these
two observed quantities in the SM are Bigi:2016mdz ; Bernlochner:2017jka ;
Jaiswal:2017rve ; Bigi:2017jbd

$R_{D}^{\rm{SM}}=0.299\pm 0.003\,,\qquad R_{D^{*}}^{\rm{SM}}=0.258\pm
0.003\,.$ (34)

These two observables have been measured independently by several
collaborations, including Babar BaBar:2012obs ; BaBar:2013mob , Belle
Belle:2015qfa ; Belle:2016dyj  and LHCb LHCb:2015gmp ; LHCb:2017smo ;
LHCb:2017rln . The average values by combining these measurements are given by
HFLAV:2019otj

$R_{D}^{\rm{exp}}=0.340\pm 0.027\pm 0.013\,,\qquad
R_{D^{*}}^{\rm{exp}}=0.295\pm 0.011\pm 0.008\,,$ (35)

which exceed the SM predictions by 1.4$\sigma$ and 2.5$\sigma$ respectively.
To confront the leptoquarks contributions with the above experimental data, we
consider the following effective Hamiltonian,

$\begin{split}\mathcal{H}_{\text{eff}}=\frac{4G_{F}}{\sqrt{2}}V_{cb}&\Big{[}g_{V_{L}}(\bar{c}_{L}\gamma_{\mu}b_{L})(\bar{\tau}_{L}\gamma^{\mu}\nu_{L})+g_{S_{L}}(\bar{c}_{R}b_{L})(\bar{\tau}_{R}\nu_{L})\\\
&+g_{T}(\bar{c}_{R}\sigma_{\mu\nu}b_{L})(\bar{\tau}_{R}\sigma^{\mu\nu}\nu_{L})\Big{]}+\text{h.c.}\,.\end{split}$
(36)

Figure 3: Upper left: Contour plot of the fit for $R_{D}$ and $R_{D^{*}}$ in
the plane of $g_{S_{L}}$ versus $g_{T}$ at the scale $\mu=1\,\text{TeV}$. The
dashed orange and dotted blue contours represent the $1\sigma$ allowed regions
that explain $R_{D}$ and $R_{D^{*}}$, respectively. The dark (light) green
region corresponds to the $1\sigma$ ($2\sigma$) allowed region that explains
both simultaneously. Upper right: The $\chi^{2}$ values to fit both $R_{D}$
and $R_{D^{*}}$ when using $g_{V_{L}}$ and $g_{S_{L}}=-4g_{T}$ respectively.
Bottom: Two plots correspond to the complex planes of $g_{S_{L}}=-4g_{T}$ and
$g_{V_{L}}$, respectively. The orange and blue regions represent the $1\sigma$
allowed region of $R_{D}$ and $R_{D^{*}}$, respectively. The deeper blue
region corresponds to the overlap scenario where accommodates combined
explanation for $R_{D}$ and $R_{D^{*}}$.

In the Fig. 3, we show the fit of $g_{S_{L}},g_{T}$ and $g_{V_{L}}$ favored
region to explain $R_{D}$ and $R_{D^{*}}$ anomalies at the scale of $1$ TeV.
The upper left panel presents the fit using real parameters $g_{S_{L}}$ and
$g_{T}$. With the relation of $g_{S_{L}}=-4g_{T}$, which is in our model, the
best fit point is $g_{S_{L}}=-4g_{T}=0.12$ and the allowed $1\,\sigma$ range
is $g_{S_{L}}=-4g_{T}\in[0.08,0.16]$. If we solely consider the Wilson
coefficient $g_{V_{L}}$, the best fit point is $g_{V_{L}}=0.08(-2.07)$ and the
allowed $1\,\sigma$ range is $g_{V_{L}}\in[0.07,0.10]\cup[-2.10,-2.05]$. We
show the $\chi^{2}$ values to fit both $R_{D}$ and $R_{D^{*}}$ in the upper
right panel. We also present the fit result when the coefficients are taken as
complex numbers. Comprehensive analyses including the ratio $R_{J/\psi}$, the
longitudinal polarization of the $P_{\tau}(D^{*})$ and $F_{L}^{D^{*}}$ can be
found in Ref. Murgui:2019czp ; Shi:2019gxi ; Cheung:2020sbq . The best fit
values in this work agree with theirs in the $1\,\sigma$ allowed range.

In the model, both $S_{1}$ and $S_{3}$ give contributions to $b\to
c\tau\bar{\nu}$ at tree-level, while $\tilde{R}_{2}$ does not. After Fierz
transformation to relevant effective Lagrangian, $g_{S_{L}}$ and $g_{T}$ have
relation $g_{S_{L}}=-4g_{T}$. The corresponding Wilson coefficient of $S_{1}$
contributions are given by

$\displaystyle
g_{V_{L}}^{\ell}=\frac{v^{2}}{4V_{cb}}\frac{(V^{T}y_{1L})^{3\ell}y_{1L}^{*23}}{m_{S_{1}}^{2}}\,,$
(37) $\displaystyle
g_{S_{L}}^{\ell}=-4g_{T}^{\ell}=-\frac{v^{2}}{4V_{cb}}\frac{(V^{T}y_{1L})^{3\ell}y_{1R}^{*23}}{m_{S_{1}}^{2}}\,.$
(38)

The contribution from leptoquark $S_{3}$ gives

$g_{V_{L}}^{\ell}=-\frac{v^{2}}{4V_{cb}}\frac{(V^{T}y_{3L})^{3\ell}y_{3L}^{*23}}{m_{S_{3}}^{2}}\,.$
(39)

However the contributions of $g_{V_{L}}^{\ell}$ from both leptoquarks $S_{1}$
and $S_{3}$ can not explain the anomalies of $R_{D^{(*)}}$ since its favored
parameters space is incompatible with $B$ meson decay process $B\to
K\nu\bar{\nu}$. To explain the anomalies of $R_{D^{(*)}}$, it is required that
$y_{1L,3L}^{33}y_{1L,3L}^{23}\sim 0.1$, but the products of couplings are
strongly constrained by the process $B\to K\nu\bar{\nu}$ with
$|y_{1L,3L}^{33}y_{1L,3L}^{23}|\lesssim 0.03$. Therefore we focus on the
Wilson coefficients $g_{S_{L}}^{\ell}$ and $g_{T}^{\ell}$ contribution from
the Leptoquark $S_{1}$ to explain the anomalies of $R_{D^{(*)}}$.

###  III.3 The anomalous magnetic moments of charged leptons

The last observable anomalies we consider are the anomalous magnetic moments
of charged leptons, including electron and muon, which both exist long-
standing discrepancy between the SM predictions and experimental measurements.
The recent combined result of Fermilab Abi:2021gix  and BNL Bennett:2006fi
increases the tension of muon $(g-2)$, which gives a $4.2\,\sigma$ level
deviation from the SM prediction. The precise discrepancy between the SM
predictions and experimental values reads Aoyama:2017uqe ; Keshavarzi:2018mgv
; Davier:2019can ; Aoyama:2020ynm ; Aoyama:2012wk ; Aoyama:2019ryr ;
Czarnecki:2002nt ; Gnendiger:2013pva ; Davier:2017zfy ; Colangelo:2018mtw ;
Hoferichter:2019mqg ; Keshavarzi:2019abf ; Kurz:2014wya ; Melnikov:2003xd ;
Masjuan:2017tvw ; Colangelo:2017fiz ; Hoferichter:2018kwz ; Gerardin:2019vio ;
Bijnens:2019ghy ; Colangelo:2019uex ; Blum:2019ugy ; Colangelo:2014qya ;
Abi:2021gix

$\displaystyle\Delta a_{e}=a_{e}^{\text{exp}}-a_{e}^{\text{SM}}=-(8.7\pm
3.6)\times 10^{-13}\,,$ (40) $\displaystyle\Delta
a_{\mu}=a_{\mu}^{\text{exp}}-a_{\mu}^{\text{SM}}=(2.51\pm 0.59)\times
10^{-9}\,.$ (41)

We start to discuss the contributions to $(g-2)_{\ell}$ from general scalar
Leptoquark interactions, which is described by Cheung:2001ip ; Dorsner:2017wwn

$\displaystyle\mathcal{L}^{F=0}$
$\displaystyle=\overline{q}_{i}(y_{R}^{ij}P_{R}+y_{L}^{ij}P_{L})\ell_{j}S+\text{h.c.}\,,$
(42) $\displaystyle\mathcal{L}^{|F|=2}$
$\displaystyle=\overline{q^{C}_{i}}(y_{R}^{\prime ij}P_{R}+y_{L}^{\prime
ij}P_{L})\ell_{j}S+\text{h.c.}\,.$ (43)

Here $q^{i}$ denotes quark, $\ell$ denotes charged leptons, $S$ stands for
leptoquarks and $F$ is the fermion number. The contributions to $\Delta
a_{\ell}\equiv(g-2)_{\ell}/2$ from $F=0$ terms are illustrated in Fig. 4 and
given by

$\Delta
a_{\ell}=-\frac{3m_{\ell}}{8\pi^{2}m_{S}^{2}}\sum_{q}\left[m_{\ell}(|y_{R}^{q\ell}|^{2}+|y_{L}^{q\ell}|^{2})F(x)+m_{q}\text{Re}(y_{L}^{*q\ell}y_{R}^{q\ell})G(x)\right]\,,$
(44)

where

$\begin{split}F(x)&=Q_{S}\,f_{S}(x)-f_{F}(x)\,,\\\
G(x)&=Q_{S}\,g_{S}(x)-g_{F}(x)\,,\end{split}$ (45)

and the loop functions are calculated by following formulas,

$\displaystyle f_{S}(x)=\frac{x+1}{4(1-x)^{2}}+\frac{x\ln x}{2(1-x)^{3}}\,,$
$\displaystyle f_{F}(x)=\frac{2+5x-x^{2}}{12(1-x)^{3}}+\frac{x\ln
x}{2(1-x)^{4}}\,,$ $\displaystyle g_{S}(x)=\frac{-1}{1-x}-\frac{\ln
x}{(1-x)^{2}}\,,$ $\displaystyle g_{F}(x)=\frac{x-3}{2(1-x)^{2}}-\frac{\ln
x}{(1-x)^{3}}\,,$ (46)

where $x=m_{q}^{2}/m_{S}^{2}$ and $Q_{S}$ is the charge of leptoquark $S$. The
$|F|=2$ scalar leptoquarks contribution can be obtained by changing the
couplings $y\to y^{\prime}$ in Eq. (44). It is noted that the no-chiral scalar
leptoquarks which have both left-handed and right-handed couplings to quarks
can give a chiral-enhanced contributions to $\Delta a_{\ell}$ by the quark
masses. This is revealed by Eq. (44), in which the first term is proportional
to the lepton mass, while the second term is proportional to the internal
quark mass. Besides, it is worthy to notice that only the second term in Eq.
(44) can provide different sign contribution, since the deviation $\Delta
a_{e}$ and $\Delta a_{\mu}$ have opposite sign. Thereby among all the scalar
leptoquarks, only singlet $S_{1}$ or doublet $R_{2}$ could provide solution to
explain $\Delta a_{e}$ and $\Delta a_{\mu}$ simultaneously. However, the
constraint from the branching ratio of $\mu\to e\gamma$ excludes the one
internal quark, such as the top quark, dominating solution Crivellin:2018qmi .
In our model, we choose the scenario that the contributions to $\Delta a_{e}$
and $\Delta a_{\mu}$ come from different quarks. The new contribution is
mainly coming from the leptoquark $S_{1}$ mediated loop and the contribution
to $\Delta a_{\ell}$ is given by

$\Delta
a_{\ell}\simeq-\sum_{q}\frac{3m_{\ell}m_{q}}{8\pi^{2}m_{S_{1}}^{2}}\text{Re}\Big{(}y_{1R}^{*q\ell}y_{1L}^{q\ell}\Big{)}\bigg{[}\frac{7}{6}+\frac{2}{3}\ln
x\bigg{]}\,.$ (47)

$\ell$$q\,(q^{C})$$\ell^{(\prime)}$$S\,(S^{*})$$\gamma$ Figure 4: One-loop
diagram contributing to the charged leptons anomalous magnetic moments and the
flavor changing process $\ell\to\ell^{\prime}\gamma$.

##  IV Low energy constraints

In the previous section, we have discussed the solution to the $B$-physics
anomalies, $R_{K^{(*)}}$ and $R_{D^{(*)}}$, and the anomalous charged lepton
magnetic moments, $\Delta a_{e}$ and $\Delta a_{\mu}$. The model also gives
rise to various flavor violating processes and rare meson decays, which are
severely constrained by current experiments. In this section, we summarize the
most stringent low-energy processes and give the relevant constraints to the
leptoquark couplings in the model.

###  IV.1 $\ell\to\ell^{\prime}\gamma$ processes

The lepton flavor violation $\ell\to\ell^{\prime}\gamma$ processes, such as
$\mu\to e\gamma$, $\tau\to e\gamma$ and $\tau\to\mu\gamma$, can be induced via
the one-loop diagrams shown in Fig. 4. The non-chiral leptoquark $S_{1}$
contribution to $\ell\to\ell^{\prime}\gamma$ processes is enhanced by the
quark mass. On the contrary, the chiral leptoquarks $\tilde{R}_{2}$ and
$S_{3}$ induce $\ell\to\ell^{\prime}\gamma$ processes without chiral
enhancement. The current experimental limits on the lepton flavor violation
$\ell\to\ell^{\prime}\gamma$ processes are summarized as following
TheMEG:2016wtm ; Aubert:2009ag ,

$\displaystyle\text{Br}(\mu\to e\gamma)<4.2\times 10^{-13}\,,$ (48)
$\displaystyle\text{Br}(\tau\to e\gamma)<3.3\times 10^{-8}\,,$ (49)
$\displaystyle\text{Br}(\tau\to\mu\gamma)<4.4\times 10^{-8}\,.$ (50)

The branching ratio of the process $\ell\to\ell^{\prime}\gamma$ mediated by
the leptoquarks can be calculated by the following formula,

$\text{Br}(\ell\to\ell^{\prime}\gamma)=\frac{\alpha_{\text{em}}(m_{\ell}^{2}-m_{\ell^{\prime}}^{2})^{3}}{4m_{\ell}^{3}\,\Gamma(\ell)}\Big{(}|\sigma_{R}^{\ell\ell^{\prime}}|^{2}+|\sigma_{L}^{\ell\ell^{\prime}}|^{2}\Big{)}\,,$
(51)

where $\Gamma(\ell)$ is the total decay width of the lepton $\ell$ and the
form factors $\sigma_{L}^{\ell\ell^{\prime}}$ and
$\sigma_{R}^{\ell\ell^{\prime}}$ originating from $S_{1}$, $\tilde{R}_{2}$ and
$S_{3}$ contribution are calculated as

$\displaystyle\sigma_{L,S_{1}}^{\ell\ell^{\prime}}=$
$\displaystyle\frac{3}{16\pi^{2}m_{S_{1}}^{2}}\sum_{q=u,c,t}\bigg{\\{}\Big{[}m_{\ell}\,y_{1R}^{q\ell}y_{1R}^{*q\ell^{\prime}}+m_{\ell^{\prime}}\,y_{1L}^{q\ell}y_{1L}^{*q\ell^{\prime}}\Big{]}\big{[}\frac{1}{3}f_{S}(x)-f_{F}(x)\big{]}$
$\displaystyle+m_{q}\big{(}y_{1L}^{q\ell}y_{1R}^{*q\ell^{\prime}}\big{)}\big{[}\frac{1}{3}\textsl{g}_{S}(x)-\textsl{g}_{F}(x)\big{]}\bigg{\\}}\,,$
(52) $\displaystyle\sigma_{R,S_{1}}^{\ell\ell^{\prime}}=$
$\displaystyle\frac{3}{16\pi^{2}m_{S_{1}}^{2}}\sum_{q=u,c,t}\bigg{\\{}\Big{[}m_{\ell}\,y_{1L}^{q\ell}y_{1L}^{*q\ell^{\prime}}+m_{\ell^{\prime}}\,y_{1R}^{q\ell}y_{1R}^{*q\ell^{\prime}}\Big{]}\big{[}\frac{1}{3}f_{S}(x)-f_{F}(x)\big{]}$
$\displaystyle+m_{q}\big{(}y_{1R}^{q\ell}y_{1L}^{*q\ell^{\prime}}\big{)}\big{[}\frac{1}{3}\textsl{g}_{S}(x)-\textsl{g}_{F}(x)\big{]}\bigg{\\}}\,,$
(53) $\displaystyle\sigma_{L,\tilde{R}_{2}}^{\ell\ell^{\prime}}=$
$\displaystyle\frac{3}{16\pi^{2}m_{R_{2}}^{2}}\sum_{q=d,s,b}m_{\ell^{\prime}}\,y_{2L}^{q\ell}y_{2L}^{*q\ell^{\prime}}\big{[}\frac{2}{3}f_{S}(x)-f_{F}(x)\big{]}\,,$
(54) $\displaystyle\sigma_{R,\tilde{R}_{2}}^{\ell\ell^{\prime}}=$
$\displaystyle\frac{3}{16\pi^{2}m_{R_{2}}^{2}}\sum_{q=d,s,b}m_{\ell}\,y_{2L}^{q\ell}y_{2L}^{*q\ell^{\prime}}\big{[}\frac{2}{3}f_{S}(x)-f_{F}(x)\big{]}\,,$
(55) $\displaystyle\sigma_{L,S_{3}}^{\ell\ell^{\prime}}=$
$\displaystyle\frac{3}{16\pi^{2}m_{S_{3}}^{2}}\bigg{\\{}\sum_{q=d,s,b}2m_{\ell^{\prime}}\,(V^{T}y_{3L})^{q\ell}(V^{\dagger}y_{3L}^{*})^{q\ell^{\prime}}\big{[}\frac{4}{3}f_{S}(x)-f_{F}(x)\big{]}$
$\displaystyle+\sum_{q=u,c,t}m_{\ell^{\prime}}\,y_{3L}^{q\ell}y_{3L}^{*q\ell^{\prime}}\big{[}\frac{1}{3}f_{S}(x)-f_{F}(x)\big{]}\bigg{\\}}\,,$
(56) $\displaystyle\sigma_{R,S_{3}}^{\ell\ell^{\prime}}=$
$\displaystyle\frac{3}{16\pi^{2}m_{S_{3}}^{2}}\bigg{\\{}\sum_{q=d,s,b}2m_{\ell}\,(V^{T}y_{3L})^{q\ell}(V^{\dagger}y_{3L}^{*})^{q\ell^{\prime}}\big{[}\frac{4}{3}f_{S}(x)-f_{F}(x)\big{]}$
$\displaystyle+\sum_{q=u,c,t}m_{\ell}\,y_{3L}^{q\ell}y_{3L}^{*q\ell^{\prime}}\big{[}\frac{1}{3}f_{S}(x)-f_{F}(x)\big{]}\bigg{\\}}\,,$
(57)

where the loop functions $f_{S,F}(x)$ and $g_{S,F}(x)$ are defined in Eqs.
(46). The terms above proportional to $m_{q}$ arising from the non-chiral
leptoquark $S_{1}$ give an enhancement and the corresponding couplings are
more severely limited. Whereas the chiral leptoquarks $\tilde{R}_{2}$ and
$S_{3}$ only consist of the terms proportional to $m_{\ell^{(\prime)}}$ and
get weaker limits. To get the constraints on the couplings, we assume that
only the relevant term dominates the contribution. The relevant constraints on
the leptoquark Yukawa couplings are summarized in Table 2.

Process | Constraints
---|---
$\mu\to e\gamma$ | $|y_{1L}^{12}y_{1R}^{*11}|,|y_{1L}^{*11}y_{1R}^{12}|<3.57\times 10^{-4}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{22}y_{1R}^{*21}|,|y_{1L}^{*21}y_{1R}^{22}|<1.29\times
10^{-6}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{32}y_{1R}^{*31}|,|y_{1L}^{*31}y_{1R}^{32}|<5.38\times
10^{-8}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{i2}y_{1L}^{*i1}|,|y_{1R}^{i2}y_{1R}^{*i1}|,|y_{3L}^{i2}y_{3L}^{*i1}|<1.31\times
10^{-3}\left(\frac{m_{S_{1}/S_{3}}}{\text{TeV}}\right)^{2}$
$|(V^{T}y_{3L})^{i2}(V^{\dagger}y_{3L}^{*})^{i1}|<3.98\times
10^{-4}\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$\tau\to e\gamma$ | $|y_{1L}^{13}y_{1R}^{*11}|,|y_{1L}^{*11}y_{1R}^{13}|<3.99\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{23}y_{1R}^{*21}|,|y_{1L}^{*21}y_{1R}^{23}|<1.45\times
10^{-2}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{33}y_{1R}^{*31}|,|y_{1L}^{*31}y_{1R}^{33}|<6.02\times
10^{-4}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{i3}y_{1L}^{*i1}|,|y_{1R}^{i3}y_{1R}^{*i1}|,|y_{3L}^{i3}y_{3L}^{*i1}|<0.874\left(\frac{m_{S_{1}/S_{3}}}{\text{TeV}}\right)^{2}$
$|(V^{T}y_{3L})^{i3}(V^{\dagger}y_{3L}^{*})^{i1}|<0.240\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$\tau\to\mu\gamma$ | $|y_{1L}^{13}y_{1R}^{*12}|,|y_{1L}^{*12}y_{1R}^{13}|<4.61\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{23}y_{1R}^{*22}|,|y_{1L}^{*22}y_{1R}^{23}|<1.67\times
10^{-2}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{33}y_{1R}^{*32}|,|y_{1L}^{*32}y_{1R}^{33}||<6.95\times
10^{-4}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{i3}y_{1L}^{*i2}|,|y_{1R}^{i3}y_{1R}^{*i2}|,|y_{3L}^{i3}y_{3L}^{*i2}|<1.01\left(\frac{m_{S_{1}/S_{3}}}{\text{TeV}}\right)^{2}$
$|(V^{T}y_{3L})^{i3}(V^{\dagger}y_{3L}^{*})^{i2}|<0.278\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$\mu\text{Au}\to e\text{Au}$ | $|y_{1L}^{12}y_{1L}^{*11}|,|y_{1R}^{12}y_{1R}^{*11}|<4.20\times 10^{-6}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{1L}^{12}y_{1R}^{*11}|,|y_{1R}^{12}y_{1L}^{*11}|<8.12\times
10^{-6}\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$|y_{2L}^{12}y_{2L}^{*11}|<3.40\times
10^{-6}\left(\frac{m_{R_{2}}}{\text{TeV}}\right)^{2}$
$|y_{3L}^{12}y_{3L}^{*11}|<2.14\times
10^{-6}\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$|(V^{T}y_{3L})^{12}(V^{\dagger}y_{3L}^{*})^{11}|<1.70\times
10^{-6}\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
Table 2: Upper limits on the leptoquark couplings from the processes
$\ell\to\ell^{\prime}\gamma$ and $\mu\text{Au}\to e\text{Au}$.

###  IV.2 $\mu-e$ conversion in nuclei

Besides the charged lepton flavor violating radiative decay processes, $\mu-e$
conversion in nuclei is also a rare process providing stringent constraints on
the strength of leptoquark interactions. The current experimental search on
$\mu-e$ conversion using gold nucleus provides the most stringent upper limits
and the upper bound to the branching ratio is set by the SINDRUM experiment as
SINDRUMII:2006dvw

$\text{Br}(\mu-e)_{\text{Au}}=\frac{\Gamma(\mu-e)_{\text{Au}}}{\Gamma_{\text{capture}}}<7\times
10^{-13}\,,$ (58)

where the $\Gamma_{\text{capture}}=8.6\times 10^{-18}$ GeV denotes the muon
capture rate by gold nucleus Suzuki:1987jf . The $\mu-e$ conversion rate in
nuclei can be calculated by following formula Kitano:2002mt ; Arganda:2007jw

$\Gamma(\mu-e)=2G_{F}^{2}m_{\mu}^{5}\left|\tilde{g}_{LS}^{(p)}S^{(p)}+\tilde{g}_{LS}^{(n)}S^{(n)}+\tilde{g}_{LV}^{(p)}V^{(p)}+\tilde{g}_{LV}^{(n)}V^{(n)}\right|^{2}+(L\to
R)\,.$ (59)

The overlap integral values of gold nucleus are $S^{(p)}=0.0523$,
$S^{(n)}=0.0610$, $V^{(p)}=0.0859$, $V^{(n)}=0.108$ Kitano:2002mt . With the
effective Lagrangian given in Eq. (21), the coupling constants $\tilde{g}$ are
defined as

$\displaystyle\tilde{g}_{LS,RS}^{(p)}=\sum_{q}G_{S}^{(q,p)}\frac{1}{2}(g_{S,q}^{LL,RR})^{ii,12}\,,$
(60)
$\displaystyle\tilde{g}_{LS,RS}^{(n)}=\sum_{q}G_{S}^{(q,n)}\frac{1}{2}(g_{S,q}^{LL,RR})^{ii,12}\,,$
(61)
$\displaystyle\tilde{g}_{LV}^{(p)}=\left[(g_{V,u}^{LL})^{11,12}+(g_{V,u}^{RL})^{11,12}\right]+\frac{1}{2}\left[(g_{V,d}^{LL})^{11,12}+(g_{V,d}^{RL})^{11,12}\right]\,,$
(62)
$\displaystyle\tilde{g}_{RV}^{(p)}=\left[(g_{V,u}^{RR})^{11,12}+(g_{V,u}^{LR})^{11,12}\right]+\frac{1}{2}\left[(g_{V,d}^{RR})^{11,12}+(g_{V,d}^{LR})^{11,12}\right]\,,$
(63)
$\displaystyle\tilde{g}_{LV}^{(n)}=\frac{1}{2}\left[(g_{V,u}^{LL})^{11,12}+(g_{V,u}^{RL})^{11,12}\right]+\left[(g_{V,d}^{LL})^{11,12}+(g_{V,d}^{RL})^{11,12}\right]\,,$
(64)
$\displaystyle\tilde{g}_{RV}^{(n)}=\frac{1}{2}\left[(g_{V,u}^{RR})^{11,12}+(g_{V,u}^{LR})^{11,12}\right]+\left[(g_{V,d}^{RR})^{11,12}+(g_{V,d}^{LR})^{11,12}\right]\,,$
(65)

where the coefficients of scalar operators are $G_{S}^{u,p}=G_{S}^{d,n}=5.1$,
$G_{S}^{d,p}=G_{S}^{u,n}=4.3$ and $G_{S}^{s,p}=G_{S}^{s,n}=2.5$ Kosmas:2001mv
. The bounds on the leptoquark couplings from $\text{Br}(\mu-e)_{\text{Au}}$
are summarized in Table 2.

###  IV.3 Rare meson leptonic decays

Introducing leptoquarks could induce meson rare decay processes. In this
subsection, we consider the relevant $B_{s}$ meson rare leptonic decays that
include leptonic conserving decays, $B_{s}\to\mu^{+}\mu^{-}/\tau^{+}\tau^{-}$,
and leptonic flavor violation decay $B_{s}\to\mu^{\pm}\tau^{\mp}$. The
corresponding 4-fermion operators
$\mathcal{O}_{9^{(\prime)}},\mathcal{O}_{10^{(\prime)}}$ are given in the Eq.
(30). The recent experimental measurements of these processes are given by
Geng:2021nhg ; LHCb:2017myy ; LHCb:2019ujz

$\displaystyle\text{Br}(B_{s}\to\mu^{+}\mu^{-})=(2.93\pm 0.35)\times
10^{-9}\,,$ (66) $\displaystyle\text{Br}(B_{s}\to\tau^{+}\tau^{-})<6.8\times
10^{-3}\,,$ (67)
$\displaystyle\text{Br}(B_{s}\to\mu^{\pm}\tau^{\mp})<1.4\times 10^{-5}\,.$
(68)

Among these processes, only the $B_{s}\to\mu^{+}\mu^{-}$ has been observed by
the current experiments and the branching ratio agrees with SM prediction
value at a level of $4\sigma$, $\text{Br}(B_{s}\to\mu\mu)^{\text{SM}}=(3.63\pm
0.13)\times 10^{-9}$ Beneke:2019slt , while the current experiments only give
upper bounds for the other two processes. The contribution to the decay width
of a neutral meson to two charged leptons $P\to\ell^{+}\ell^{\prime-}$ can be
written as Becirevic:2016zri

$\displaystyle\Gamma_{P\to\ell^{+}\ell^{\prime-}}=$
$\displaystyle\frac{1}{64\pi^{3}}\frac{G_{F}^{2}\alpha^{2}_{\text{em}}}{m_{P}^{3}}f_{P}^{2}|V_{qj}V^{*}_{qi}|^{2}\lambda_{1}^{1/2}\lambda_{2}^{1/2}$
$\displaystyle\times\Bigg{\\{}\lambda_{1}\cdot\left|(m_{\ell}-m_{\ell^{\prime}})\left(C_{9}^{ij\ell\ell^{\prime}}-C_{9^{\prime}}^{ij\ell\ell^{\prime}}\right)+\frac{m_{P}^{2}}{m_{q}+m_{q^{\prime}}}(C_{S}-C_{S}^{\prime})\right|^{2}$
$\displaystyle+\lambda_{2}\cdot\left|(m_{\ell}+m_{\ell^{\prime}})\left(C_{10}^{ij\ell\ell^{\prime}}-C_{10^{\prime}}^{ij\ell\ell^{\prime}}\right)+\frac{m_{P}^{2}}{m_{q}+m_{q^{\prime}}}(C_{P}-C_{P}^{\prime})\right|^{2}\Bigg{\\}}\,,$
(69)

where $f_{P}$ is the meson decay constant,
$\lambda_{1,2}=m_{P}^{2}-(m_{\ell}\pm m_{\ell^{\prime}})^{2}$ and
$m_{q},m_{q^{\prime}}$ are the masses of the valence quarks in the
pseudoscalar meson $P$. It is noted that the lepton flavor conserving decay
process $P\to\ell^{+}\ell^{-}$ is independent of Wilson coefficients
$C_{9^{(\prime)}}$.

Processes | Constraints on the couplings
---|---
$B^{0}_{s}\to\mu^{+}\mu^{-}$ | $y_{2L}^{32}y_{2L}^{*22}\in[-1.1,1.1]\times 10^{-4}\times\left(\frac{m_{R_{2}}}{\text{TeV}}\right)^{2}$
$(V^{T}y_{3L})^{22}(V^{\dagger}y_{3L}^{*})^{32}\in[-5.8,5.3]\times
10^{-4}\times\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$B^{0}_{s}\to\tau^{+}\tau^{-}$ | $y_{2L}^{33}y_{2L}^{*23}\in[-1.3,1.3]\times\left(\frac{m_{R_{2}}}{\text{TeV}}\right)^{2}$
$(V^{T}y_{3L})^{23}(V^{\dagger}y_{3L}^{*})^{33}\in[-0.63,0.63]\times\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$B^{0}_{s}\to\mu^{\pm}\tau^{\mp}$ | $y_{2L}^{32}y_{2L}^{*23},y_{2L}^{33}y_{2L}^{*22}\in[-0.080,0.080]\times\left(\frac{m_{R_{2}}}{\text{TeV}}\right)^{2}$
$(V^{T}y_{3L})^{23}(V^{\dagger}y_{3L}^{*})^{32},(V^{T}y_{3L})^{22}(V^{\dagger}y_{3L}^{*})^{33}\in[-0.040,0.040]\times\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$B\to K\nu\bar{\nu}$ | $(V^{T}y_{1L})^{3i}(V^{\dagger}y_{1L}^{*})^{2j}\in[-0.070,0.029]\times\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$y_{2L}^{3i}y_{2L}^{*2j}\in[-0.032,0.061]\times\left(\frac{m_{R_{2}}}{\text{TeV}}\right)^{2}$
$(V^{T}y_{3L})^{3i}(V^{\dagger}y_{3L}^{*})^{2j}\in[-0.070,0.029]\times\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$B_{s}^{0}-\overline{B_{s}^{0}}$ | $(V^{T}y_{1L})^{2i}(V^{\dagger}y_{1L}^{*})^{3i}\in[-0.14,0.14]\times\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$y_{2L}^{2i}y_{2L}^{*3i}\in[-0.14,0.14]\times\left(\frac{m_{R_{2}}}{\text{TeV}}\right)^{2}$
$(V^{T}y_{3L})^{2i}(V^{\dagger}y_{3L}^{*})^{3i}\in[-0.061,0.061]\times\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
$K^{0}-\overline{K^{0}}$ | $(V^{T}y_{1L})^{2i}(V^{\dagger}y_{1L}^{*})^{1i}\in[-0.026,0.026]\times\left(\frac{m_{S_{1}}}{\text{TeV}}\right)^{2}$
$y_{2L}^{2i}y_{2L}^{*1i}\in[-0.026,0.026]\times\left(\frac{m_{R_{2}}}{\text{TeV}}\right)^{2}$
$(V^{T}y_{3L})^{2i}(V^{\dagger}y_{3L}^{*})^{1i}\in[-0.013,0.013]\times\left(\frac{m_{S_{3}}}{\text{TeV}}\right)^{2}$
Table 3: Bounds on the leptoquark couplings from neutral mesons mixing and
rare decay processes.

###  IV.4 Rare meson semi-leptonic decays

The meson rare semi-leptonic decays can be induced at the tree-level by the
leptoquarks and present constraints on the corresponding parameters. Here we
consider $B\to K\nu\bar{\nu}$ and $B\to K^{*}\nu\bar{\nu}$ processes, related
to the $(\bar{q}q\bar{\nu}\nu)$ interactions. The corresponding SM predictions
are $\text{Br}(B^{0}\to K^{0}\nu\nu)=(4.1\pm 0.5)\times 10^{-6}$ and
$\text{Br}(B^{0}\to K^{*0}\nu\nu)=(9.2\pm 1.0)\times 10^{-6}$ Buras:2014fpa ;
Altmannshofer:2009ma , while the current experimental upper limit bounds are
given as $2.6\times 10^{-5}$ and $1.8\times 10^{-5}$ by the Belle
collaboration Belle:2017oht  respectively. To describe the constraints on new
physics from the $B\to K\nu\bar{\nu}$ and $B\to K^{*}\nu\bar{\nu}$ processes,
the ratio $R_{K^{(*)}}^{\nu\nu}$ is introduced and defined as

$R_{K^{(*)}}^{\nu\nu}=\frac{\text{Br}^{\text{SM+NP}}(B\to
K^{(*)}\nu\bar{\nu})}{\text{Br}^{\text{SM}}(B\to K^{(*)}\nu\bar{\nu})}\,.$
(70)

The latest Belle results Belle:2017oht  imply $R_{K}^{\nu\nu}<3.9$ and
$R_{K^{*}}^{\nu\nu}<2.7\,$. As shown in Table 1, the contributions to $B\to
K^{(*)}\nu\bar{\nu}$ from leptoquarks $S_{1}$ and $S_{3}$ are represented by
the Wilson coefficients $h_{V,d}^{LL}$, while by the Wilson coefficient
$h_{V,d}^{RL}$ for the case of leptoquark $\tilde{R}_{2}$. If the new physics
contribution is dominated by the $h_{V,d}^{LL}$ term, the ratios $R_{K}$ and
$R_{K^{*}}$ can be calculated by the following formula Buras:2014fpa ,

$R_{K^{(*)}}^{\nu\nu}=\frac{2}{3}+\sum_{\nu,\nu^{\prime}}\frac{1}{3|C_{L}^{\text{SM}}|^{2}}\left|\delta^{\nu\nu^{\prime}}C_{L}^{\text{SM}}+(h_{V,d}^{LL})^{32;\nu\nu^{\prime}}\right|^{2}\,,$
(71)

where $C_{L}^{\text{SM}}$ describes the SM contribution and the value is
$C_{L}^{\text{SM}}=-6.35$. Note that since the experiments cannot detect the
neutrinos in the final state, we need sum over all the flavor. On the other
hand, if the new physics contribution is only originated from the
$h_{V,d}^{RL}$ term, one has $R_{K}^{\nu\nu}\neq R_{K^{*}}^{\nu\nu}$ and the
ratios are then presented by

$\displaystyle
R_{K}^{\nu\nu}=\frac{2}{3}+\sum_{\nu,\nu^{\prime}}\frac{1}{3|C_{L}^{\text{SM}}|^{2}}\left[\delta^{\nu\nu^{\prime}}C_{L}^{\text{SM}}+(h_{V,d}^{RL})^{32;\nu\nu^{\prime}}\right]\left[1+2\frac{\delta^{\nu\nu^{\prime}}C_{L}^{\text{SM}}\,\text{Re}(h_{V,d}^{RL})^{32;\nu\nu^{\prime}}}{|C_{L}^{\text{SM}}|^{2}+|(h_{V,d}^{RL})^{32;\nu\nu^{\prime}}|^{2}}\right]\,,$
(72) $\displaystyle
R_{K^{*}}^{\nu\nu}=\frac{2}{3}+\sum_{\nu,\nu^{\prime}}\frac{1}{3|C_{L}^{\text{SM}}|^{2}}\left[\delta^{\nu\nu^{\prime}}C_{L}^{\text{SM}}+(h_{V,d}^{RL})^{32;\nu\nu^{\prime}}\right]\left[1-1.34\frac{\delta^{\nu\nu^{\prime}}C_{L}^{\text{SM}}\,\text{Re}(h_{V,d}^{RL})^{32;\nu\nu^{\prime}}}{|C_{L}^{\text{SM}}|^{2}+|(h_{V,d}^{RL})^{32;\nu\nu^{\prime}}|^{2}}\right]\,.$
(73)

###  IV.5 Neutral meson mixing

Leptoquarks can induce neutral meson mixing via box diagrams mediated by
leptons and leptoquarks. In this subsection we study the constraints from the
$B_{s}^{0}-\overline{B_{s}^{0}}$ and $K^{0}-\overline{K^{0}}$ mixing. The
related effective Hamiltonian can be described by Bobeth:2017ecx

$\mathcal{H}_{\text{eff}}=C^{ij}_{LL}(\bar{d}_{L}^{i}\gamma^{\mu}d_{L}^{j})(\bar{d}_{L}^{i}\gamma_{\mu}d_{L}^{j})+C^{ij}_{RR}(\bar{d}_{R}^{i}\gamma^{\mu}d_{R}^{j})(\bar{d}_{R}^{i}\gamma_{\mu}d_{R}^{j})+C^{ij}_{LR}(\bar{d}_{L}^{i}\gamma^{\mu}d_{L}^{j})(\bar{d}_{R}^{i}\gamma_{\mu}d_{R}^{j})\,,$
(74)

where $i,j=3,2$ corresponding to $B_{s}^{0}-\overline{B_{s}^{0}}$ mixing and
$i,j=2,1$ related to $K^{0}-\overline{K^{0}}$ mixing. Mapping the contribution
of leptoquarks $S_{1},\tilde{R}_{2}$ and $S_{3}$, we have the Wilson
coefficients at the scale $\mu=m_{\text{LQ}}$ in the following form,

$\displaystyle S_{1}:\qquad
C^{ij}_{LL}=-\frac{1}{128\pi^{2}m_{S_{1}}^{2}}\sum_{k}\left[(V^{T}y_{1L})^{ik}(V^{\dagger}y_{1L}^{*})^{jk}\right]^{2}\,,$
(75) $\displaystyle\tilde{R}_{2}:\qquad
C^{ij}_{RR}=-\frac{1}{128\pi^{2}m_{R_{2}}^{2}}\sum_{k}2(y_{2L}^{ik}y_{2L}^{*jk})^{2}\,,$
(76) $\displaystyle S_{3}:\qquad
C^{ij}_{LL}=-\frac{1}{128\pi^{2}m_{S_{3}}^{2}}\sum_{k}5\left[(V^{T}y_{3L})^{ik}(V^{\dagger}y_{3L}^{*})^{jk}\right]^{2}\,.$
(77)

The transition of the Wilson coefficients from $\mu=1\,\text{TeV}$ to
$\mu=m_{b}$ scale are evaluated by the Wilson package Aebischer:2018bkb  and
the results are given by

$C^{ij}_{LL,RR}(\mu=1\,\text{TeV})=0.78\,C^{ij}_{LL,RR}(\mu=m_{b})\,.$ (78)

The current measurements of the mass differences in
$B_{s}^{0}-\overline{B_{s}^{0}}$ and $K^{0}-\overline{K^{0}}$ mixing are
ParticleDataGroup:2020ssz ,

$\displaystyle\Delta m_{B_{s}}^{\text{exp}}=(17.741\pm 0.020)\times
10^{12}\,\text{sec}^{-1}\,,$ (79) $\displaystyle\Delta
m_{K}^{\text{exp}}=(3.484\pm 0.0009)\times 10^{10}\,\text{sec}^{-1}\,.$ (80)

For the mass difference $\Delta m_{B_{s}}^{\text{exp}}$, the SM prediction
value is $\Delta m_{B_{s}}^{\text{SM}}=(18.3\pm 2.7)\times
10^{12}\,\text{sec}^{-1}$ Artuso:2015swg ; FermilabLattice:2016ipl ;
Jubb:2016mvq . But the SM prediction for the mass difference in
$K^{0}-\overline{K^{0}}$ mixing has not been precisely estimated Brod:2011ty ;
Buras:2014maa . Thereby in our analysis, we take the new physics contribution
to $K^{0}-\overline{K^{0}}$ mixing to be compatible with the experimental
value. The bounds on the leptoquarks couplings from neutral meson mixing and
rare decay processes are summarized in Table 3.

##  V Numerical analysis

Figure 5: Allowed regions of the various leptoquark Yukawa couplings obtained
from fitting to the corresponding processes. The dark (light) green band
represent the $1\sigma$ ($2\sigma$) allowed regions that explain the
corresponding anomalies and the lines denote the constraints from the labelled
processes.

In this section, we perform a numerical analysis of the model parameter space
to supply a common explanation of $B$-physics anomalies in $R_{K^{(*)}}$,
$R_{D^{(*)}}$ and the charged leptons anomalous magnetic moment
$(g-2)_{e,\mu}$, as well as the neutrino oscillation data. Instead of
exploring the entire parameter space, we find the minimal parameters of the
model and combine the constraints from the low-energy processes given in the
Sec. IV. We fix the components of singlet leptoquark $S_{1}$ and triplet
leptoquark $S_{3}$ mass at 1 TeV ($m_{S_{1}}=m_{S_{3}}=1\,\text{TeV}$) and fix
the components of doublet leptoquark $\tilde{R}_{2}$ mass at 2 TeV
($m_{R_{2}}=2\ \text{TeV}$). We use the python package Flavio to obtain the
appropriate values of Wilson coefficients that explain the anomalies at the
scale of leptoquark masses ($\mu=1\ \text{TeV}$) and then analyze the model
parameter space.

In order to minimize the number of parameters, we adopt the following form of
the Yukawa coupling matrices in the analysis.

$\displaystyle y_{1R}=\begin{pmatrix}0\ &0\ &0\\\
{\color[rgb]{0,0.75,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.75,1}y_{1R}^{21}}&0\
&{\color[rgb]{1,0.6484375,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6484375,0}y_{1R}^{23}}\\\
0\
&{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}y_{1R}^{32}}&0\end{pmatrix}\,,$
$\displaystyle\qquad y_{1L}=\begin{pmatrix}0\ &0\ &0\ \\\
{\color[rgb]{0,0.75,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.75,1}y_{1L}^{21}}&0\
&y_{1L}^{23}\\\ 0\
&{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}y_{1L}^{32}}&{\color[rgb]{1,0.6484375,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6484375,0}y_{1L}^{33}}\end{pmatrix}\,,$
$\displaystyle y_{3L}=\begin{pmatrix}0\ &0&0\ \\\ 0\
&{\color[rgb]{1,0.26953125,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.26953125,0}y_{3L}^{22}}&y_{3L}^{23}\
\\\ y_{3L}^{31}\
&{\color[rgb]{1,0.26953125,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.26953125,0}y_{3L}^{32}}&0\
\end{pmatrix}\,,$ $\displaystyle\qquad y_{2L}=\begin{pmatrix}0\ &0\ &0\ \\\
y_{2L}^{21}\ &y_{2L}^{22}\ &y_{2L}^{23}\\\ y_{2L}^{31}\ &y_{2L}^{32}\
&y_{2L}^{33}\end{pmatrix}\,.$ (81)

The coupling combination $(y_{3L}^{22},\ y_{3L}^{32})$ can explain the
anomalies of $R_{K}$ and $R_{K^{*}}$, while the couplings $(y_{2R}^{33},\
y_{2L}^{23})$ can contribute to $R_{D}$ and $R_{D^{*}}$. The couplings
$(y_{1R}^{21},\ y_{1L}^{21})$ and $(y_{1R}^{32},\ y_{1L}^{32})$ give
contributions to $\Delta a_{e}$ and $\Delta a_{\mu}$ respectively. The other
non-zero couplings are needed to fit the neutrino masses and mixing angles.
For a simple illustration, in Fig. 5, we present the allowed parameter space
to explain these anomalies and satisfy the relevant processes constraints with
taking the coupling as real. Specifically, we provide two concrete benchmark
points of the leptoquarks Yukawa couplings. In benchmark point 1, the
couplings are chosen as complex number. While for benchmark point 2, we choose
the Dirac CP angle in the neutrino mixing matrix as $180^{\circ}$, which is
within $1\,\sigma$ allowed range ParticleDataGroup:2020ssz  and it is possible
to take all the leptoquark Yukawa coupling values as real. The corresponding
values of observables for these two benchmark points are summarized in Table
4.

Benchmark point 1:

$\displaystyle y_{1R}=\begin{pmatrix}0\,&0\,&0\,\\\ -0.37\,&0\,&-0.70\,\\\
0\,&0.054\,&0\,\end{pmatrix}\,,\qquad y_{3L}=\begin{pmatrix}0\,&0\,&0\,\\\
0\,&0.029\,&0\,\\\ 0\,&0.023\,&0\,\end{pmatrix}\,,$ $\displaystyle
y_{1L}=\begin{pmatrix}0\ &0\,&0\,\\\ 0.012+0.016i\,&0\,&-0.049-0.0042i\,\\\
0\,&0.57+0.0082i\,&0.59+0.052i\,\end{pmatrix}\,,$ $\displaystyle
y_{2L}=\begin{pmatrix}0\,&0\,&0\,\\\ 0.043-0.042i\,&0.044\,&-0.048\,\\\
-0.00013\,&0.00038\,&0.00027\,\end{pmatrix}\,.$ (82)

Benchmark point 2:

$\displaystyle y_{1R}=\begin{pmatrix}0\,&0\,&0\,\\\ -0.014\,&0\,&-0.93\,\\\
0\,&0.012\,&0\,\end{pmatrix}\,,\ \ \ \ y_{1L}=\begin{pmatrix}0\ &0\,&0\,\\\
0.37\,&0\,&0\,\\\ 0\,&0.21\,&0.38\,\end{pmatrix}\,,$ $\displaystyle
y_{3L}=\begin{pmatrix}0\,&0\,&0\,\\\ 0\,&0.070\,&0.038\,\\\
0.0019\,&0.0059\,&0\,\end{pmatrix}\,,\ \ \ \ \
y_{2L}=\begin{pmatrix}0\,&0\,&0\,\\\ -0.015\,&-0.0020\,&0.023\,\\\
0.0015\,&0.0031\,&-0.00078\,\end{pmatrix}\,.$ (83)

Observables | Allowed range | BP1 | BP2
---|---|---|---
$\Delta m_{21}^{2}(10^{-5}\text{eV}^{2})$ | [6.82, 8.04] | 7.44 | 7.40
$\Delta m_{32}^{2}(10^{-3}\text{eV}^{2})$ | [2.435, 2.598] | 2.50 | 2.51
$\sin^{2}\theta_{12}$ | [0.269, 0.343] | 0.305 | 0.301
$\sin^{2}\theta_{23}$ | [0.405, 0.620] | 0.569 | 0.570
$\sin^{2}\theta_{13}$ | [0.02064, 0.02430] | 0.0226 | 0.0225
$\delta_{\text{CP}}/^{\circ}$ | [169, 246] | 194 | 180
$R_{K}$ | [0.795, 0.901] | 0.808 | 0.812
$R_{K^{*}}$ | [0.569, 0.845] | 0.794 | 0.817
$R_{K_{S}^{0}}$ | [0.48, 0.88] | 0.808 | 0.812
$R_{K^{*+}}$ | [0.53, 0.91] | 0.825 | 0.844
$R_{D}$ | [0.310, 0.370] | 0.358 | 0.351
$R_{D^{*}}$ | [0.281, 0.309] | 0.305 | 0.292
$\Delta a_{e}\,(10^{-13})$ | [-12.3, -5.1] | -8.48 | -9.89
$\Delta a_{\mu}\,(10^{-9})$ | [1.93, 3.11] | 2.52 | 2.06
Br$(\mu\to e\gamma)$ | $<4.2\times 10^{-13}$ | $6.49\times 10^{-21}$ | $1.04\times 10^{-17}$
Br$(\tau\to e\gamma)$ | $<3.3\times 10^{-8}$ | $2.98\times 10^{-18}$ | $6.09\times 10^{-16}$
Br$(\tau\to\mu\gamma)$ | $<4.4\times 10^{-8}$ | $2.99\times 10^{-18}$ | $7.82\times 10^{-16}$
Br$(\mu-e)_{\text{Au}}$ | $<7\times 10^{-13}$ | $4.84\times 10^{-19}$ | $2.62\times 10^{-15}$
Br$(B^{0}_{s}\to\mu\mu)$ | [2.58, 3.28] $\times 10^{-9}$ | $2.93\times 10^{-9}$ | $3.42\times 10^{-9}$
Br$(B^{0}_{s}\to\tau\tau)$ | $<6.8\times 10^{-3}$ | $7.82\times 10^{-7}$ | $7.95\times 10^{-7}$
Br$(B^{0}_{s}\to\mu\tau)$ | $<1.4\times 10^{-5}$ | $2.11\times 10^{-14}$ | $1.40\times 10^{-10}$
$R_{K}^{\nu\nu}$ | $<3.9$ | 1.3 | 0.75
$R_{K^{*}}^{\nu\nu}$ | $<2.7$ | 1.4 | 0.76
$\Delta m_{B_{s}}^{\text{SM+NP}}/\Delta m_{B_{s}}^{\text{SM}}$ | [0.85,1.15] | 1.03 | 1.01
$\Delta m_{K}^{\text{NP}}(10^{10}\,s^{-1})$ | $<0.95$ | 0.0016 | 0.52
Table 4: Summary of the observable values for the benchmark points.

##  VI Conclusion

In this paper, we have proposed a simple model by extending SM with three TeV-
scale scalar leptoquarks $S_{1}$, $\tilde{R}_{2}$ and $S_{3}$, where the
source of tiny neutrino masses, the lepton flavor anomalies in $B$-meson
decays $(R_{K^{(*)}}$, $R_{D^{(*)}})$ and the tension in the charged lepton
(electron and muon) anomalous magnetic moments have a common solution. In the
model, $R_{K^{(*)}}$ anomalies are resolved by the leptoquark $S_{3}$ via the
Wilson coefficients $C_{9,10}^{\mu\mu}$. Leptoquark $S_{1}$ explains the
anomalies of $R_{D^{(*)}}$ through the Wilson coefficients $g_{S_{L}},g_{T}$,
as well as the deviations of leptonic magnetic moments $(g-2)_{e,\mu}$ by one-
loop level contribution. The small mixing of leptoquarks $S_{1}$ with
$\tilde{R}_{2}$ or $\tilde{R}_{2}$ with $S_{3}$ can generate tiny neutrino
masses. We analyze the parameter space of the leptoquark Yukawa couplings and
obtain the corresponding viable region. We study the relevant experimental
constraints and conclude there is an appropriate parameter space accommodate
to combined explanation for these anomalies and deviations.

Acknowledgements. This work is supported in part by the National Science
Foundation of China (12175082, 11775093) and the Fundamental Research Funds
for the Central Universities (CCNU22LJ004).

## References

  * (1) Super-Kamiokande collaboration, _Evidence for oscillation of atmospheric neutrinos_ , _Phys. Rev. Lett._ 81 (1998) 1562 [hep-ex/9807003].
  * (2) SNO collaboration, _Direct evidence for neutrino flavor transformation from neutral current interactions in the Sudbury Neutrino Observatory_ , _Phys. Rev. Lett._ 89 (2002) 011301 [nucl-ex/0204008].
  * (3) KamLAND collaboration, _First results from KamLAND: Evidence for reactor anti-neutrino disappearance_ , _Phys. Rev. Lett._ 90 (2003) 021802 [hep-ex/0212021].
  * (4) T2K collaboration, _Indication of Electron Neutrino Appearance from an Accelerator-produced Off-axis Muon Neutrino Beam_ , _Phys. Rev. Lett._ 107 (2011) 041801 [1106.2822].
  * (5) E. Ma, _Pathways to naturally small neutrino masses_ , _Phys. Rev. Lett._ 81 (1998) 1171 [hep-ph/9805219].
  * (6) P. Minkowski, _$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?_, _Phys. Lett. B_ 67 (1977) 421.
  * (7) R.N. Mohapatra and G. Senjanovic, _Neutrino Mass and Spontaneous Parity Nonconservation_ , _Phys. Rev. Lett._ 44 (1980) 912.
  * (8) M. Gell-Mann, P. Ramond and R. Slansky, _Complex Spinors and Unified Theories_ , _Conf. Proc. C_ 790927 (1979) 315 [1306.4669].
  * (9) S.L. Glashow, _The Future of Elementary Particle Physics_ , _NATO Sci. Ser. B_ 61 (1980) 687.
  * (10) J. Schechter and J.W.F. Valle, _Neutrino Masses in SU(2) x U(1) Theories_ , _Phys. Rev. D_ 22 (1980) 2227.
  * (11) R.N. Mohapatra and G. Senjanovic, _Neutrino Masses and Mixings in Gauge Models with Spontaneous Parity Violation_ , _Phys. Rev. D_ 23 (1981) 165.
  * (12) G. Lazarides, Q. Shafi and C. Wetterich, _Proton Lifetime and Fermion Masses in an SO(10) Model_ , _Nucl. Phys. B_ 181 (1981) 287.
  * (13) R. Foot, H. Lew, X.G. He and G.C. Joshi, _Seesaw Neutrino Masses Induced by a Triplet of Leptons_ , _Z. Phys. C_ 44 (1989) 441.
  * (14) A. Zee, _A Theory of Lepton Number Violation, Neutrino Majorana Mass, and Oscillation_ , _Phys. Lett. B_ 93 (1980) 389.
  * (15) T.P. Cheng and L.-F. Li, _Neutrino Masses, Mixings and Oscillations in SU(2) x U(1) Models of Electroweak Interactions_ , _Phys. Rev. D_ 22 (1980) 2860.
  * (16) A. Zee, _Quantum Numbers of Majorana Neutrino Masses_ , _Nucl. Phys. B_ 264 (1986) 99.
  * (17) K.S. Babu, _Model of ’Calculable’ Majorana Neutrino Masses_ , _Phys. Lett. B_ 203 (1988) 132.
  * (18) Y. Cai, J. Herrero-García, M.A. Schmidt, A. Vicente and R.R. Volkas, _From the trees to the forest: a review of radiative neutrino mass models_ , _Front. in Phys._ 5 (2017) 63 [1706.08524].
  * (19) T. Aoyama, T. Kinoshita and M. Nio, _Revised and Improved Value of the QED Tenth-Order Electron Anomalous Magnetic Moment_ , _Phys. Rev. D_ 97 (2018) 036001 [1712.06060].
  * (20) A. Keshavarzi, D. Nomura and T. Teubner, _Muon $g-2$ and $\alpha(M_{Z}^{2})$: a new data-based analysis_, _Phys. Rev._ D97 (2018) 114025 [1802.02995].
  * (21) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, _A new evaluation of the hadronic vacuum polarisation contributions to the muon anomalous magnetic moment and to $\mathbf{\bm{\alpha}(m_{Z}^{2})}$_, _Eur. Phys. J._ C80 (2020) 241 [1908.00921].
  * (22) T. Aoyama et al., _The anomalous magnetic moment of the muon in the Standard Model_ , _Phys. Rept._ 887 (2020) 1 [2006.04822].
  * (23) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, _Complete Tenth-Order QED Contribution to the Muon $g-2$_, _Phys. Rev. Lett._ 109 (2012) 111808 [1205.5370].
  * (24) T. Aoyama, T. Kinoshita and M. Nio, _Theory of the Anomalous Magnetic Moment of the Electron_ , _Atoms_ 7 (2019) 28.
  * (25) A. Czarnecki, W.J. Marciano and A. Vainshtein, _Refinements in electroweak contributions to the muon anomalous magnetic moment_ , _Phys. Rev._ D67 (2003) 073006 [hep-ph/0212229].
  * (26) C. Gnendiger, D. Stöckinger and H. Stöckinger-Kim, _The electroweak contributions to $(g-2)_{\mu}$ after the Higgs boson mass measurement_, _Phys. Rev._ D88 (2013) 053005 [1306.5546].
  * (27) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, _Reevaluation of the hadronic vacuum polarisation contributions to the Standard Model predictions of the muon $g-2$ and ${\alpha(m_{Z}^{2})}$ using newest hadronic cross-section data_, _Eur. Phys. J._ C77 (2017) 827 [1706.09436].
  * (28) G. Colangelo, M. Hoferichter and P. Stoffer, _Two-pion contribution to hadronic vacuum polarization_ , _JHEP_ 02 (2019) 006 [1810.00007].
  * (29) M. Hoferichter, B.-L. Hoid and B. Kubis, _Three-pion contribution to hadronic vacuum polarization_ , _JHEP_ 08 (2019) 137 [1907.01556].
  * (30) A. Keshavarzi, D. Nomura and T. Teubner, _The $g-2$ of charged leptons, $\alpha(M_{Z}^{2})$ and the hyperfine splitting of muonium_, _Phys. Rev._ D101 (2020) 014029 [1911.00367].
  * (31) A. Kurz, T. Liu, P. Marquard and M. Steinhauser, _Hadronic contribution to the muon anomalous magnetic moment to next-to-next-to-leading order_ , _Phys. Lett._ B734 (2014) 144 [1403.6400].
  * (32) K. Melnikov and A. Vainshtein, _Hadronic light-by-light scattering contribution to the muon anomalous magnetic moment revisited_ , _Phys. Rev._ D70 (2004) 113006 [hep-ph/0312226].
  * (33) P. Masjuan and P. Sánchez-Puertas, _Pseudoscalar-pole contribution to the $(g_{\mu}-2)$: a rational approach_, _Phys. Rev._ D95 (2017) 054026 [1701.05829].
  * (34) G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, _Dispersion relation for hadronic light-by-light scattering: two-pion contributions_ , _JHEP_ 04 (2017) 161 [1702.07347].
  * (35) M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold and S.P. Schneider, _Dispersion relation for hadronic light-by-light scattering: pion pole_ , _JHEP_ 10 (2018) 141 [1808.04823].
  * (36) A. Gérardin, H.B. Meyer and A. Nyffeler, _Lattice calculation of the pion transition form factor with $N_{f}=2+1$ Wilson quarks_, _Phys. Rev._ D100 (2019) 034520 [1903.09471].
  * (37) J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, _Short-distance constraints for the HLbL contribution to the muon anomalous magnetic moment_ , _Phys. Lett._ B798 (2019) 134994 [1908.03331].
  * (38) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, _Longitudinal short-distance constraints for the hadronic light-by-light contribution to $(g-2)_{\mu}$ with large-$N_{c}$ Regge models_, _JHEP_ 03 (2020) 101 [1910.13432].
  * (39) T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung et al., _The hadronic light-by-light scattering contribution to the muon anomalous magnetic moment from lattice QCD_ , _Phys. Rev. Lett._ 124 (2020) 132002 [1911.08123].
  * (40) G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera and P. Stoffer, _Remarks on higher-order hadronic corrections to the muon $g-2$_, _Phys. Lett._ B735 (2014) 90 [1403.7512].
  * (41) Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * (42) BaBar collaboration, _Evidence for an excess of $\bar{B}\to D^{(*)}\tau^{-}\bar{\nu}_{\tau}$ decays_, _Phys. Rev. Lett._ 109 (2012) 101802 [1205.5442].
  * (43) BaBar collaboration, _Measurement of an Excess of $\bar{B}\to D^{(*)}\tau^{-}\bar{\nu}_{\tau}$ Decays and Implications for Charged Higgs Bosons_, _Phys. Rev. D_ 88 (2013) 072012 [1303.0571].
  * (44) Belle collaboration, _Measurement of the branching ratio of $\bar{B}\to D^{(\ast)}\tau^{-}\bar{\nu}_{\tau}$ relative to $\bar{B}\to D^{(\ast)}\ell^{-}\bar{\nu}_{\ell}$ decays with hadronic tagging at Belle_, _Phys. Rev. D_ 92 (2015) 072014 [1507.03233].
  * (45) Belle collaboration, _Measurement of the $\tau$ lepton polarization and $R(D^{*})$ in the decay $\bar{B}\to D^{*}\tau^{-}\bar{\nu}_{\tau}$_, _Phys. Rev. Lett._ 118 (2017) 211801 [1612.00529].
  * (46) Belle collaboration, _Measurement of the branching ratio of $\bar{B}^{0}\rightarrow D^{*+}\tau^{-}\bar{\nu}_{\tau}$ relative to $\bar{B}^{0}\rightarrow D^{*+}\ell^{-}\bar{\nu}_{\ell}$ decays with a semileptonic tagging method_, in _51st Rencontres de Moriond on EW Interactions and Unified Theories_ , 3, 2016 [1603.06711].
  * (47) LHCb collaboration, _Measurement of the ratio of branching fractions $\mathcal{B}(B_{c}^{+}\,\to\,J/\psi\tau^{+}\nu_{\tau})$/$\mathcal{B}(B_{c}^{+}\,\to\,J/\psi\mu^{+}\nu_{\mu})$_, _Phys. Rev. Lett._ 120 (2018) 121801 [1711.05623].
  * (48) LHCb collaboration, _Measurement of the ratio of the $B^{0}\to D^{*-}\tau^{+}\nu_{\tau}$ and $B^{0}\to D^{*-}\mu^{+}\nu_{\mu}$ branching fractions using three-prong $\tau$-lepton decays_, _Phys. Rev. Lett._ 120 (2018) 171802 [1708.08856].
  * (49) LHCb collaboration, _Search for lepton-universality violation in $B^{+}\to K^{+}\ell^{+}\ell^{-}$ decays_, _Phys. Rev. Lett._ 122 (2019) 191801 [1903.09252].
  * (50) LHCb collaboration, _Test of lepton universality with $B^{0}\rightarrow K^{*0}\ell^{+}\ell^{-}$ decays_, _JHEP_ 08 (2017) 055 [1705.05802].
  * (51) LHCb collaboration, _Test of lepton universality in beauty-quark decays_ , _Nature Phys._ 18 (2022) 277 [2103.11769].
  * (52) I. Doršner, S. Fajfer, A. Greljo, J.F. Kamenik and N. Košnik, _Physics of leptoquarks in precision experiments and at particle colliders_ , _Phys. Rept._ 641 (2016) 1 [1603.04993].
  * (53) A. Crivellin, C. Greub, D. Müller and F. Saturnino, _Scalar Leptoquarks in Leptonic Processes_ , _JHEP_ 02 (2021) 182 [2010.06593].
  * (54) A. Carvunis, A. Crivellin, D. Guadagnoli and S. Gangal, _The Forward-Backward Asymmetry in $B\to D^{*}\ell\nu$: One more hint for Scalar Leptoquarks?_, _Phys. Rev. D_ 105 (2022) L031701 [2106.09610].
  * (55) A. Angelescu, D. Bečirević, D.A. Faroughy and O. Sumensari, _Closing the window on single leptoquark solutions to the $B$-physics anomalies_, _JHEP_ 10 (2018) 183 [1808.08179].
  * (56) I. Bigaran, J. Gargalionis and R.R. Volkas, _A near-minimal leptoquark model for reconciling flavour anomalies and generating radiative neutrino masses_ , _JHEP_ 10 (2019) 106 [1906.01870].
  * (57) S. Saad, _Combined explanations of $(g-2)_{\mu}$, $R_{D^{(*)}}$, $R_{K^{(*)}}$ anomalies in a two-loop radiative neutrino mass model_, _Phys. Rev. D_ 102 (2020) 015019 [2005.04352].
  * (58) V. Gherardi, D. Marzocca and E. Venturini, _Low-energy phenomenology of scalar leptoquarks at one-loop accuracy_ , _JHEP_ 01 (2021) 138 [2008.09548].
  * (59) A. Greljo, P. Stangl and A.E. Thomsen, _A model of muon anomalies_ , _Phys. Lett. B_ 820 (2021) 136554 [2103.13991].
  * (60) H.M. Lee, _Leptoquark option for B-meson anomalies and leptonic signatures_ , _Phys. Rev. D_ 104 (2021) 015007 [2104.02982].
  * (61) A. Bhaskar, A.A. Madathil, T. Mandal and S. Mitra, _Combined explanation of $W$-mass, muon $g-2$, $R_{K^{(*)}}$ and $R_{D^{(*)}}$ anomalies in a singlet-triplet scalar leptoquark model_, 2204.09031.
  * (62) C.-H. Chen, T. Nomura and H. Okada, _Explanation of $B\to K^{(*)}\ell^{+}\ell^{-}$ and muon $g-2$, and implications at the LHC_, _Phys. Rev. D_ 94 (2016) 115005 [1607.04857].
  * (63) S. Saad and A. Thapa, _Common origin of neutrino masses and $R_{D^{(\ast)}}$, $R_{K^{(\ast)}}$ anomalies_, _Phys. Rev. D_ 102 (2020) 015014 [2004.07880].
  * (64) K.S. Babu, P.S.B. Dev, S. Jana and A. Thapa, _Unified framework for $B$-anomalies, muon $g-2$ and neutrino masses_, _JHEP_ 03 (2021) 179 [2009.01771].
  * (65) I. Doršner, S. Fajfer and S. Saad, _$\mu\to e\gamma$ selecting scalar leptoquark solutions for the $(g-2)_{e,\mu}$ puzzles_, _Phys. Rev. D_ 102 (2020) 075007 [2006.11624].
  * (66) I. Bigaran and R.R. Volkas, _Getting chirality right: Single scalar leptoquark solutions to the $(g-2)_{e,\mu}$ puzzle_, _Phys. Rev. D_ 102 (2020) 075037 [2002.12544].
  * (67) D. Aristizabal Sierra, M. Hirsch and S.G. Kovalenko, _Leptoquarks: Neutrino masses and accelerator phenomenology_ , _Phys. Rev. D_ 77 (2008) 055011 [0710.5699].
  * (68) I. Doršner, S. Fajfer and N. Košnik, _Leptoquark mechanism of neutrino masses within the grand unification framework_ , _Eur. Phys. J. C_ 77 (2017) 417 [1701.08322].
  * (69) J. Julio, S. Saad and A. Thapa, _A flavor-inspired radiative neutrino mass model_ , _JHEP_ 08 (2022) 270 [2202.10479].
  * (70) J. Julio, S. Saad and A. Thapa, _Marriage between neutrino mass and flavor anomalies_ , _Phys. Rev. D_ 106 (2022) 055003 [2203.15499].
  * (71) T.A. Chowdhury and S. Saad, _Leptoquark-vectorlike quark model for the CDF mW, (g-2) $\mu$, RK(*) anomalies, and neutrino masses_, _Phys. Rev. D_ 106 (2022) 055017 [2205.03917].
  * (72) I. Doršner, S. Fajfer and O. Sumensari, _Triple-leptoquark interactions for tree- and loop-level proton decays_ , _JHEP_ 05 (2022) 183 [2202.08287].
  * (73) I. Cordero-Carrión, M. Hirsch and A. Vicente, _Master Majorana neutrino mass parametrization_ , _Phys. Rev. D_ 99 (2019) 075019 [1812.03896].
  * (74) I. Cordero-Carrión, M. Hirsch and A. Vicente, _General parametrization of Majorana neutrino mass models_ , _Phys. Rev. D_ 101 (2020) 075032 [1912.08858].
  * (75) K.G. Chetyrkin, _Quark mass anomalous dimension to O (alpha-s**4)_ , _Phys. Lett. B_ 404 (1997) 161 [hep-ph/9703278].
  * (76) J.A. Gracey, _Three loop MS-bar tensor current anomalous dimension in QCD_ , _Phys. Lett. B_ 488 (2000) 175 [hep-ph/0007171].
  * (77) I. Doršner, S. Fajfer, N. Košnik and I. Nišandžić, _Minimally flavored colored scalar in $\bar{B}\to D^{(*)}\tau\bar{\nu}$ and the mass matrices constraints_, _JHEP_ 11 (2013) 084 [1306.6493].
  * (78) G. Hiller, D. Loose and K. Schönwald, _Leptoquark Flavor Patterns & B Decay Anomalies_, _JHEP_ 12 (2016) 027 [1609.08895].
  * (79) J. Aebischer, J. Kumar and D.M. Straub, _Wilson: a Python package for the running and matching of Wilson coefficients above and below the electroweak scale_ , _Eur. Phys. J. C_ 78 (2018) 1026 [1804.05033].
  * (80) D.M. Straub, _flavio: a Python package for flavour and precision phenomenology in the Standard Model and beyond_ , 1810.08132.
  * (81) C. Bobeth, G. Hiller and G. Piranishvili, _Angular distributions of $\bar{B}\to\bar{K}\ell^{+}\ell^{-}$ decays_, _JHEP_ 12 (2007) 040 [0709.4174].
  * (82) M. Bordone, G. Isidori and A. Pattori, _On the Standard Model predictions for $R_{K}$ and $R_{K^{*}}$_, _Eur. Phys. J. C_ 76 (2016) 440 [1605.07633].
  * (83) G. Buchalla, A.J. Buras and M.E. Lautenbacher, _Weak decays beyond leading logarithms_ , _Rev. Mod. Phys._ 68 (1996) 1125 [hep-ph/9512380].
  * (84) B. Capdevila, A. Crivellin, S. Descotes-Genon, J. Matias and J. Virto, _Patterns of New Physics in $b\to s\ell^{+}\ell^{-}$ transitions in the light of recent data_, _JHEP_ 01 (2018) 093 [1704.05340].
  * (85) T. Hurth, F. Mahmoudi, D. Martinez Santos and S. Neshatpour, _Lepton nonuniversality in exclusive $b{\rightarrow}s{\ell}{\ell}$ decays_, _Phys. Rev. D_ 96 (2017) 095034 [1705.06274].
  * (86) J. Aebischer, W. Altmannshofer, D. Guadagnoli, M. Reboud, P. Stangl and D.M. Straub, _$B$ -decay discrepancies after Moriond 2019_, _Eur. Phys. J. C_ 80 (2020) 252 [1903.10434].
  * (87) LHCb collaboration, _Tests of lepton universality using $B^{0}\to K^{0}_{S}\ell^{+}\ell^{-}$ and $B^{+}\to K^{*+}\ell^{+}\ell^{-}$ decays_, _Phys. Rev. Lett._ 128 (2022) 191802 [2110.09501].
  * (88) L.-S. Geng, B. Grinstein, S. Jäger, S.-Y. Li, J. Martin Camalich and R.-X. Shi, _Implications of new evidence for lepton-universality violation in b→s $\ell$+$\ell$\- decays_, _Phys. Rev. D_ 104 (2021) 035029 [2103.12738].
  * (89) ATLAS collaboration, _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].
  * (90) CMS collaboration, _Measurement of properties of B ${}^{0}_{\mathrm{s}}\to\mu^{+}\mu^{-}$ decays and search for B${}^{0}\to\mu^{+}\mu^{-}$ with the CMS experiment_, _JHEP_ 04 (2020) 188 [1910.12127].
  * (91) 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].
  * (92) LHCb collaboration, _Analysis of Neutral B-Meson Decays into Two Muons_ , _Phys. Rev. Lett._ 128 (2022) 041801 [2108.09284].
  * (93) 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].
  * (94) A.K. Alok, N.R.S. Chundawat, S. Gangal and D. Kumar, _A global analysis of $b\to s\ell\ell$ data in heavy and light $Z^{\prime}$ models_, 2203.13217.
  * (95) D. Bečirević, S. Fajfer and N. Košnik, _Lepton flavor nonuniversality in b→s $\ell$+$\ell$- processes_, _Phys. Rev. D_ 92 (2015) 014016 [1503.09024].
  * (96) A. Carvunis, F. Dettori, S. Gangal, D. Guadagnoli and C. Normand, _On the effective lifetime of B s→ $\mu$$\mu$$\gamma$_, _JHEP_ 12 (2021) 078 [2102.13390].
  * (97) W. Altmannshofer and P. Stangl, _New physics in rare B decays after Moriond 2021_ , _Eur. Phys. J. C_ 81 (2021) 952 [2103.13370].
  * (98) D. Bečirević, N. Košnik, O. Sumensari and R. Zukanovich Funchal, _Palatable Leptoquark Scenarios for Lepton Flavor Violation in Exclusive $b\to s\ell_{1}\ell_{2}$ modes_, _JHEP_ 11 (2016) 035 [1608.07583].
  * (99) D. Bigi and P. Gambino, _Revisiting $B\to D\ell\nu$_, _Phys. Rev. D_ 94 (2016) 094008 [1606.08030].
  * (100) F.U. Bernlochner, Z. Ligeti, M. Papucci and D.J. Robinson, _Combined analysis of semileptonic $B$ decays to $D$ and $D^{*}$: $R(D^{(*)})$, $|V_{cb}|$, and new physics_, _Phys. Rev. D_ 95 (2017) 115008 [1703.05330].
  * (101) S. Jaiswal, S. Nandi and S.K. Patra, _Extraction of $|V_{cb}|$ from $B\to D^{(*)}\ell\nu_{\ell}$ and the Standard Model predictions of $R(D^{(*)})$_, _JHEP_ 12 (2017) 060 [1707.09977].
  * (102) D. Bigi, P. Gambino and S. Schacht, _$R(D^{*})$ , $|V_{cb}|$, and the Heavy Quark Symmetry relations between form factors_, _JHEP_ 11 (2017) 061 [1707.09509].
  * (103) LHCb collaboration, _Measurement of the ratio of branching fractions $\mathcal{B}(\bar{B}^{0}\to D^{*+}\tau^{-}\bar{\nu}_{\tau})/\mathcal{B}(\bar{B}^{0}\to D^{*+}\mu^{-}\bar{\nu}_{\mu})$_, _Phys. Rev. Lett._ 115 (2015) 111803 [1506.08614].
  * (104) LHCb collaboration, _Test of Lepton Flavor Universality by the measurement of the $B^{0}\to D^{*-}\tau^{+}\nu_{\tau}$ branching fraction using three-prong $\tau$ decays_, _Phys. Rev. D_ 97 (2018) 072013 [1711.02505].
  * (105) HFLAV collaboration, _Averages of b-hadron, c-hadron, and $\tau$-lepton properties as of 2018_, _Eur. Phys. J. C_ 81 (2021) 226 [1909.12524].
  * (106) C. Murgui, A. Peñuelas, M. Jung and A. Pich, _Global fit to $b\to c\tau\nu$ transitions_, _JHEP_ 09 (2019) 103 [1904.09311].
  * (107) R.-X. Shi, L.-S. Geng, B. Grinstein, S. Jäger and J. Martin Camalich, _Revisiting the new-physics interpretation of the $b\to c\tau\nu$ data_, _JHEP_ 12 (2019) 065 [1905.08498].
  * (108) K. Cheung, Z.-R. Huang, H.-D. Li, C.-D. Lü, Y.-N. Mao and R.-Y. Tang, _Revisit to the $b\to c\tau\nu$ transition: In and beyond the SM_, _Nucl. Phys. B_ 965 (2021) 115354 [2002.07272].
  * (109) Muon g-2 collaboration, _Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL_ , _Phys. Rev. D_ 73 (2006) 072003 [hep-ex/0602035].
  * (110) K.-m. Cheung, _Muon anomalous magnetic moment and leptoquark solutions_ , _Phys. Rev. D_ 64 (2001) 033001 [hep-ph/0102238].
  * (111) A. Crivellin, M. Hoferichter and P. Schmidt-Wellenburg, _Combined explanations of $(g-2)_{\mu,e}$ and implications for a large muon EDM_, _Phys. Rev. D_ 98 (2018) 113002 [1807.11484].
  * (112) MEG collaboration, _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].
  * (113) BaBar collaboration, _Searches for Lepton Flavor Violation in the Decays $\tau^{\pm}\to e^{\pm}\gamma$ and $\tau^{\pm}\to\mu^{\pm}\gamma$_, _Phys. Rev. Lett._ 104 (2010) 021802 [0908.2381].
  * (114) SINDRUM II collaboration, _A Search for muon to electron conversion in muonic gold_ , _Eur. Phys. J. C_ 47 (2006) 337.
  * (115) T. Suzuki, D.F. Measday and J.P. Roalsvig, _Total Nuclear Capture Rates for Negative Muons_ , _Phys. Rev. C_ 35 (1987) 2212.
  * (116) R. Kitano, M. Koike and Y. Okada, _Detailed calculation of lepton flavor violating muon electron conversion rate for various nuclei_ , _Phys. Rev. D_ 66 (2002) 096002 [hep-ph/0203110].
  * (117) E. Arganda, M.J. Herrero and A.M. Teixeira, _mu-e conversion in nuclei within the CMSSM seesaw: Universality versus non-universality_ , _JHEP_ 10 (2007) 104 [0707.2955].
  * (118) T.S. Kosmas, S. Kovalenko and I. Schmidt, _Nuclear muon- e- conversion in strange quark sea_ , _Phys. Lett. B_ 511 (2001) 203 [hep-ph/0102101].
  * (119) LHCb collaboration, _Search for the decays $B_{s}^{0}\to\tau^{+}\tau^{-}$ and $B^{0}\to\tau^{+}\tau^{-}$_, _Phys. Rev. Lett._ 118 (2017) 251802 [1703.02508].
  * (120) LHCb collaboration, _Search for the lepton-flavour-violating decays $B^{0}_{s}\to\tau^{\pm}\mu^{\mp}$ and $B^{0}\to\tau^{\pm}\mu^{\mp}$_, _Phys. Rev. Lett._ 123 (2019) 211801 [1905.06614].
  * (121) D. Bečirević, O. Sumensari and R. Zukanovich Funchal, _Lepton flavor violation in exclusive $b\rightarrow s$ decays_, _Eur. Phys. J. C_ 76 (2016) 134 [1602.00881].
  * (122) A.J. Buras, J. Girrbach-Noe, C. Niehoff and D.M. Straub, _$B\to{K}^{\left(\ast\right)}\nu\overline{\nu}$ decays in the Standard Model and beyond_, _JHEP_ 02 (2015) 184 [1409.4557].
  * (123) W. Altmannshofer, A.J. Buras, D.M. Straub and M. Wick, _New strategies for New Physics search in $B\to K^{*}\nu\bar{\nu}$, $B\to K\nu\bar{\nu}$ and $B\to X_{s}\nu\bar{\nu}$ decays_, _JHEP_ 04 (2009) 022 [0902.0160].
  * (124) Belle collaboration, _Search for $\bm{B\to h\nu\bar{\nu}}$ decays with semileptonic tagging at Belle_, _Phys. Rev. D_ 96 (2017) 091101 [1702.03224].
  * (125) C. Bobeth and A.J. Buras, _Leptoquarks meet $\varepsilon^{\prime}/\varepsilon$ and rare Kaon processes_, _JHEP_ 02 (2018) 101 [1712.01295].
  * (126) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * (127) M. Artuso, G. Borissov and A. Lenz, _CP violation in the $B_{s}^{0}$ system_, _Rev. Mod. Phys._ 88 (2016) 045002 [1511.09466].
  * (128) Fermilab Lattice, MILC collaboration, _$B^{0}_{(s)}$ -mixing matrix elements from lattice QCD for the Standard Model and beyond_, _Phys. Rev. D_ 93 (2016) 113016 [1602.03560].
  * (129) T. Jubb, M. Kirk, A. Lenz and G. Tetlalmatzi-Xolocotzi, _On the ultimate precision of meson mixing observables_ , _Nucl. Phys. B_ 915 (2017) 431 [1603.07770].
  * (130) J. Brod and M. Gorbahn, _Next-to-Next-to-Leading-Order Charm-Quark Contribution to the $CP$ Violation Parameter $\epsilon_{K}$ and $\Delta M_{K}$_, _Phys. Rev. Lett._ 108 (2012) 121801 [1108.2036].
  * (131) A.J. Buras, J.-M. Gérard and W.A. Bardeen, _Large $N$ Approach to Kaon Decays and Mixing 28 Years Later: $\Delta I=1/2$ Rule, $\hat{B}_{K}$ and $\Delta M_{K}$_, _Eur. Phys. J. C_ 74 (2014) 2871 [1401.1385].


Submitted to the 28${}^{\text{th}}$ ACM SIGPLAN International Conference on
Functional Programming, Haskell Symposium.
Accepted for oral presentation.
Presented Sep 8, 2023.

parallel computing, syntax tree parsing, work-stealing scheduler

††copyright: none††journalyear: 2023††doi: XXXXXXX.XXXXXXX††conference: ACM
SIGPLAN Haskell Symposium 2023; Sep 4–9, 2023; Seattle, Washington,
USA††booktitle: Proceedings of the 16th ACM SIGPLAN International Haskell
Symposium (Haskell ’23), Sep 4–9, 2023, Seattle, Washington, USA††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Computing methodologies Parallel
computing methodologies

##  1\. Introduction

In this paper, we designed, implemented, and benchmarked a Haskell auto-
parallelizer with a simple yet powerful interface by taking advantage of the
default purity of Haskell functions. One of the main challenges in distributed
computing is building interfaces and APIs that allow programmers with limited
background in distributed systems to write scalable, performant, and fault-
tolerant applications on large clusters. While accounting for the complex
side-effects induced by the separate pieces of code running on different
computers, designers of these frameworks often have to find the compromise
between more limited API or looser consistencies and semantics. For example,
MapReduce (Dean and Ghemawat, 2004) provides the programming models of Map and
Reduce that completely abstracts the distributed nature of the program away
from the programmer, with the caveat that side-effects have to be atomic and
idempotent. Ray, on the other hand, asks the programmer to decide whether the
execution should be stateful or not. In large projects, it could be very
difficult for a programmer to reason about these properties because some
library deep down in the nested function calls could introduce unexpected side
effects, such as writing to the same temporary file across all calls, which
would result in incorrect execution of the program.

By design, the purity of Haskell mitigates these problems drastically and
allows us to provide a much simpler interface than traditionally possible with
other object-oriented languages. Since the purity of a function call can be
directly inferred from its type signature at compile time, we can parallelize
the pure functions across the workers without worrying about the side-effects
while also ensuring that the impure functions are executed in order. More
specifically, given a Haskell program and a section of the code to
parallelize, a scheduler can parse the program’s data dependencies between
function calls and greedily schedules tasks to worker nodes as their inputs
are ready, all without requiring the user to understand the way in which
underlying libraries execute.

We have built a prototype implementing the ideas above. While our prototype
does not yet have industrial-scale scalability or fault-tolerance, it has
shown promises as a simple way to achieve parallel speed-up on existing
workloads. Hence, in this paper, we benchmarked our implementation to
illustrate the potential for future work in this direction.

##  2\. Design

The user specifies which section of the code to parallelize. This should be a
section with high levels of abstraction (in other words, each function call
takes some amount of time to execute). For example, a user working on natural
language processing could write something like the following:

    data Summary = ...  -- A custom data type

    clean_files :: IO Summary
    clean_files = ...

    complex_evaluation :: Summary -> Int
    complex_evaluation x = ...

    semantic_analysis :: IO Int
    semantic_analysis = ...

    main :: IO()
    main = do
        x <- clean_files
        let y = complex_evaluation x
        z <- semantic_analysis
        print (y, z)

In this case, the user should specify that the `main` function is the function
they want to parallelize, because it calls on other high-level functions. From
there, a parser could infer the following data dependency graph:

Figure 1. Data dependency graph generated from example Haskell program

Notice that `RealWorld` is considered an input and output by each `IO`
function. In this case, once `clean_files` is done, both `complex_evaluation`
and `semantic_analysis` can be scheduled for execution.

Another example is a deep learning project, in which the user specify the
forward and backward passes of the neural network. In practice, a user is
usually writing high-level code, so the section worth parallelizing usually
coincides with the section that the user writes, which is very convenient. In
our prototype, only the main function is parallelized, but we could see a
future implementation where the user can specify any arbitrary function.

While we only implemented a shallow parser in our implementation, we think
that incorporating a more powerful parser such as Graph Trace that can parse
arbitrary depth could further allow the user to specify the granularity of
distribution.

##  3\. Related Works

This section provides both works that tread similar paths as us and works that
a future implementation of this idea could build on.

There are a number of automatic parallelization compilers for other languages.
These compilers parallelize at instruction level on a shared memory machine.
Many such research compilers exist for Fortran and C, such as the Vienna
Fortran Compiler (Benkner, 1999) and the Intel C++ Compiler (Intel, 2021). A
recent work by Google (Zheng et al., 2022) enables automatic parallelization
of models written under the JAX framework (Bradbury et al., 2018).

There are a number of existing packages for shared-memory parallel in Haskell.
Built-in to Haskell are the Control.concurrent package provides basic APIs to
threads and forks, and the parallel library provides primitives such as par
and pseq. Projects such as monad-par (Marlow et al., 2011) and lvish (Kuper
and Newton, 2013) provide additional interfaces through monads to enable work-
stealing schedulers.

A few libraries exist as backbones for distributed computing. Cloud Haskell
(Epstein et al., 2011) provides APIs for serializing functional closures and
channels for network communications. The network package provides a low-level
networking interface. A recent advancement is Sparkle (I/O, 2023), which
provides a Haskell interface to Spark (Zaharia et al., 2012), a unified
framework for resilient distributed applications written in Scala. We believe
future implementations could benefit from Sparkle for its robustness.

Finally, while we wrote a custom script for generating dependency graph based
on a shallow parsing of the program in our implementation for simplicity,
there are robust libraries that can generate functional graphs for complex
Haskell programs. Two that we found are Graph Trace and SourceGraph, which
future implementations could benefit from.

##  4\. Experiments

We performed experiments on matrix operations (generation and multiplication
of large random matrices) with different numbers of workers simulated using
Cloud Haskell. Although this may seem simple, it is the foundation for modern
deep learning computations and good performance and scalability would indicate
great potential for future generalization. For reference, we used single-
thread and Haskell’s built-in SMP parallelism as our baselines. The task size
is the number of times that the matrix operations are performed, and all
results are rounded to the nearest second.

Figure 2. Benchmark results on large matrix multiplication tasks

###### Acknowledgements.

We would like to thank Chengrui Hou, Danyang Zhuo and Simon Peyton-Jones for
their generous feedbacks.

## References

  * (1)
  * Benkner (1999) Siegfried Benkner. 1999\.  VFC: The Vienna Fortran Compiler.  _Scientific Programming_ 7, 1 (1999), 67–81.
  * Bradbury et al. (2018) James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. 2018\.  _JAX: composable transformations of Python+NumPy programs_.  http://github.com/google/jax
  * Dean and Ghemawat (2004) Jeffrey Dean and Sanjay Ghemawat. 2004.  MapReduce: Simplified Data Processing on Large Clusters. In _USENIX Symposium on Operating Systems Design and Implementation (OSDI)_ _(OSDI’04)_. 137–150.  https://www.usenix.org/legacy/publications/library/proceedings/osdi04/tech/full_papers/dean/dean.pdf
  * Epstein et al. (2011) Jeff Epstein, Andrew P Black, and Simon Peyton-Jones. 2011\.  Towards Haskell in the cloud. In _Proceedings of the 4th ACM symposium on Haskell_. 118–129.
  * Intel (2021) Intel. 2021.  _Intel® C++ Compiler Classic Developer Guide and Reference_.  Intel.
  * I/O (2023) Tweag I/O. 2023\.  sparkle: Apache Spark applications in Haskell.  https://github.com/charlespwd/project-title.
  * Kuper and Newton (2013) Lindsey Kuper and Ryan R Newton. 2013.  LVars: lattice-based data structures for deterministic parallelism. In _Proceedings of the 2nd ACM SIGPLAN workshop on Functional high-performance computing_. 71–84.
  * Marlow et al. (2011) Simon Marlow, Ryan Newton, and Simon Peyton Jones. 2011\.  A monad for deterministic parallelism.  _ACM SIGPLAN Notices_ 46, 12 (2011), 71–82.
  * Zaharia et al. (2012) Matei Zaharia, Mosharaf Chowdhury, Tathagata Das, Ankur Dave, Justin Ma, Murphy McCauly, Michael J Franklin, Scott Shenker, and Ion Stoica. 2012.  Resilient distributed datasets: A fault-tolerant abstraction for in-memory cluster computing. In _Presented as part of the 9th $\\{$USENIX$\\}$ Symposium on Networked Systems Design and Implementation ($\\{$NSDI$\\}$ 12)_. 15–28.
  * Zheng et al. (2022) Lianmin Zheng, Zhuohan Li, Hao Zhang, Yonghao Zhuang, Zhifeng Chen, Yanping Huang, Yida Wang, Yuanzhong Xu, Danyang Zhuo, Eric P Xing, et al. 2022\.  Alpa: Automating Inter-and $\\{$Intra-Operator$\\}$ Parallelism for Distributed Deep Learning. In _16th USENIX Symposium on Operating Systems Design and Implementation (OSDI 22)_. 559–578.