# Study of $B_{s}\to\phi\ell^{+}\ell^{-}$ decays in the PQCD factorization
approach with lattice QCD input

Su-Ping Jin1<EMAIL_ADDRESS>Zhen-Jun Xiao1,2<EMAIL_ADDRESS>Corresponding author 1\. Department of Physics and Institute of Theoretical
Physics, Nanjing Normal University, Nanjing, Jiangsu 210023, People’s Republic
of China,  2\. Jiangsu Key Laboratory for Numerical Simulation of Large Scale
Complex Systems, Nanjing Normal University, Nanjing 210023, People’s Republic
of China

###### Abstract

In this paper, we studied systematically the semileptonic decays $B_{s}\to\phi
l^{+}l^{-}$ with $l^{-}=(e^{-},\mu^{-},\tau^{-})$ by using the perturbative
QCD (PQCD) and the “PQCD+Lattice” factorization approach, respectively. We
first evaluated all relevant form factors $F_{i}(q^{2})$ in the low $q^{2}$
region using the PQCD approach, and we also took the available lattice QCD
results at the high-$q^{2}$ region as additional input to improve the
extrapolation of $F_{i}(q^{2})$ from the low-$q^{2}$ region to the endpoint
$q^{2}_{max}$. We then calculated the branching ratios and many other physical
observables: $A_{FB}^{l}$, $F_{L}^{\phi}$, $S_{3,4,7}$, $A_{5,6,8,9}$ and the
clean angular observables $P_{1,2,3}$ and $P^{\prime}_{4,5,6,8}$. From our
studies, we find the following points: (a) the PQCD and “PQCD+Lattice”
predictions of ${\cal B}(B_{s}\to\phi\mu^{+}\mu^{-})$ are about $7\times
10^{-7}$, which agree well with the LHCb measured values and the QCD sum rule
prediction within still large errors; (b) we defined and calculated the ratios
of the branching ratios $R_{\phi}^{e\mu}$ and $R_{\phi}^{\mu\tau}$; (c) the
PQCD and ”PQCD+Lattice” predictions of the longitudinal polarization $F_{L}$,
the CP averaged angular coefficients $S_{3,4,7}$ and the CP asymmetry angular
coefficients $A_{5,6,8,9}$, agree with the LHCb measurements in all considered
bins within the still large experimental errors; and (d) for those currently
still unknown observables
$R_{\phi}^{e\mu},R_{\phi}^{\mu\tau},A_{FB}^{l},P_{1,2,3}$ and
$P^{\prime}_{4,5,6,8}$, we suggest LHCb and Belle-II Collaboration to measure
them in their experiments.

###### pacs:

13.20.He, 12.38.Bx, 14.40.Nd

##  I Introduction

In the Standard Model (SM) of particle physics, one treats these three
generations of the charged leptons $\ell^{-}=(e^{-},\mu^{-},\tau^{-})$ as
exact copies of each other. These charged leptons behave in the same way but
differ only in the masses determined by their Yukawa coupling to the Higgs
boson. The lepton flavor universality (LFU), i.e. the equality of the coupling
to the all electroweak gauge bosons among three families of leptons, has been
regarded as an exact symmetry for quite a long timeDey:2018arx . In recent
years, however, some physics observables associated with the flavor-changing
neutral current (FCNC) transitions $b\to s\ell\ell$ have exhibited deviations
from the SM expectations. These include the LFU-violating(LFUV) ratios $R_{K}$
and $R_{K^{\ast}}$ Aaij:2014ora ; Aaij:2017vbb , whose measurements deviates
from $\mu-e$ universality Bordone:2016gaq ; Hiller:2003js ; Hiller2015  by
around $2.5\sigma$. More notably, the measurements of the angular observable
$P^{\prime}_{5}$ of $B\to K^{\ast}\mu^{+}\mu^{-}$ decay in the large recoil
region Descotes-Genon:2013wba ; Descotes-Genon:2013vna ; Matias:2012xw ;
DescotesGenon:2012zf ; Matias:2014jua  as reported by the LHCb Aaij:2013qta ;
Aaij:2015oid  and Belle Collaboration Abdesselam:2016llu  point to a deviation
of about $3\sigma$ with respect to the SM prediction Descotes-Genon:2014uoa .

As is well known, the FCNC $b\to s$ transition is forbidden at tree-level, but
proceeds by way of loop diagrams with a very low rate. Due to the strong
suppression within SM, such kinds of FCNC decays may be sensitive to the
possible new physics (NP) effects. Therefore, the semileptonic $b\to
s\ell\ell$ decay has received striking attentions by means of measurements of
the inclusive $B\to X_{s}\ell^{+}\ell^{-}$ and/or the exclusive $B\to
K^{(\ast)}\ell^{+}\ell^{-}$ decays and their comparison with the SM
predictions. Besides the decay rates, many angular observables of the
semileptonic $B\to K^{\ast}\mu^{+}\mu^{-}$ decays have also been measured
previously Aaij:2013qta ; Aaij:2015oid ; Abdesselam:2016llu . The precision of
the experimental measurements will also be expected to upgrade remarkably in
the forthcoming year.

The semileptonic decay $B_{s}\to\phi\mu^{+}\mu^{-}$, which is closely relevant
to the decay $B\to K^{\ast}\mu^{+}\mu^{-}$, offers an alternative scene to
check out the same fundamental quark process, in a different hadronic
background. On the theoretical side, various studies on the quark level $b\to
s$ transition and the exclusive $B_{(s)}\to V\ell^{+}\ell^{-}$ decays by using
rather different theories or models have been performed within the SM, such as
the constituent quark model or covariant quark model Deandrea:2001qs ;
Dubnicka:2016nyy , the light front quark model Geng:2003su , the QCD
factorization (QCDF) Bobeth:2008ij  and the light-cone sum rule (LCSR)
Ball:2004rg ; Altmannshofer:2014rta ; Wu:2006rd ; Straub:2015ica ; Gao2020 ;
Descotes-Genon:2015uva , and beyond the SM, such as the universal extra
dimension Mohanta:2006ae ; Li:2011yn  and the supersymmetric theory Yuan-
Guo:2013vpa . On the experimental side, the $B_{s}\to\phi\mu^{+}\mu^{-}$ decay
mode was first observed and studied by the CDF collaboration Aaltonen:2011cn
and subsequently by the LHCb collaboration Aaij:2013aln ; Aaij:2015esa ;
Aaij:2103 ; Aaij:2021pkz ; LHCb-2107 . Beyond the measurement of the branching
ratio, a rich phenomenology of various kinematical distributions can be
presented. While the angular distributions was found to be consistent with the
SM expectations obtained in Refs. Altmannshofer:2014rta ; Straub:2015ica ,
however, LHCb also observed a deficit with respect to the SM prediction for
the branching ratio $B^{0}_{s}\to\phi\mu^{+}\mu^{-}$ in the low-$q^{2}$
region: the tension between the theory and experiment is about $3\sigma$ in
the region $1.0\leq q^{2}\leq 6.0GeV^{2}$, where the form factors are
evaluated by using the combined fit of lattice and the LCSR results
Altmannshofer:2014rta ; Straub:2015ica .

In a previous paper Jin:2020jtu , the semileptonic $B_{s}\to
K^{(\ast)}\ell^{+}\ell^{+}$ decays have been studied by us using the
perturbative QCD (PQCD) factorization approach Keum:2000wi ; Lu:2000em ;
li2003 ; Li:2009pr ; Fan:2012kn ; Fan:2013qz ; Wang:2012ab ; Wang:2013ix ;
Xiao:2011tx ; Ali:2007ff ; xiao18a ; xiao18b . In this paper, we will make
systematic studies for the the semileptonic $B_{s}\to\phi\ell^{+}\ell^{-}$ and
present the theoretical predictions for many physical observables:

  1. (1)

For $B_{s}\to\phi\ell^{+}\ell^{-}$ decays, we treat them as a four body decay
$B_{s}\to\phi(\to K^{-}K^{+})\ell^{+}\ell^{-}$ described by four kinematic
variables: the lepton invariant mass squared $q^{2}$ and three angles
$({\theta_{K}},{\theta_{\ell}},\Phi)$. We defined and calculated the full
angular decay distribution, the transverse amplitudes, the partially
integrated decay amplitudes over the angles
$({\theta_{K}},{\theta_{\ell}},\Phi)$, the CP averaged differential branching,
the ratios $R^{e\mu}(\phi)$ and $R^{\mu\tau}(\phi)$ of the branching ratios,
the forward-backward asymmetry $A_{FB}(q^{2})$, the $\phi$ polarization
fraction $F_{L}(q^{2})$, the CP averaged (asymmetry) angular coefficients
$S_{i}$ ($A_{i}$) and the optimized observables $P_{i}$ and $P_{i}^{\prime}$.
Following Ref. Doring2013  , where the authors approved that the possible
S-wave correction to the branching fractions of $B_{s}\to\phi(\to
K^{-}K^{+})\ell^{+}\ell^{-}$ decays is small and may modify the differential
decay widths by about $5\%$ only, we therefore will take the S-wave correction
to the branching fractions as an additional uncertainty of $5\%$ in magnitude.

  2. (2)

We used both the PQCD factorization approach and the “PQCD+Lattice” approach
to determine the values and their $q^{2}$-dependence of the $B_{s}\to\phi$
transition form factors. We used the $z$-series parametrization to make the
extrapolation for all form factors from the low $q^{2}$ region to the endpoint
$q^{2}_{max}$. We will calculate the branching ratios and all other physical
observables by using the PQCD approach itself and the “PQCD+Lattice” approach
respectively, and compare their predictions with those currently available
experimental measurements.

The paper is organized as follows: In Sec. II, we give a short review for the
kinematics of the $B_{s}\to\phi\ell^{+}\ell^{-}$ decays including distribution
amplitudes of $B_{s}$ and $\phi$ mesons, and the effective Hamiltonian for the
quark level $b\to s\ell^{+}\ell^{-}$. In Sec. III, we define explicitly all
physical observables for $B_{s}\to\phi\ell^{+}\ell^{-}$ decays. In Sec. IV we
present our theoretical predictions of all relevant physical observables of
the considered decay modes, compare these predictions with those currently
available experimental measurements and make some phenomenological analysis. A
short summary is given in the last section.

##  II Kinematics and theoretical framework

###  II.1  Kinematics and wave functions

We treat the $B_{s}$ meson at rest as a heavy-light system. The kinematics of
the semileptonic $B_{s}\to\phi\ell^{+}\ell^{-}$ decays in the large-recoil
(low $q^{2}$) region will be discussed below, where the PQCD factorization
approach is applicable to the considered decays. In the rest frame of
$\bar{B}^{0}_{s}$ meson, we define the $\bar{B}^{0}_{s}$ meson momentum
$p_{1}$, the $\phi$ momentum $p_{2}$ in the light-cone coordinates as Ref.
Fan:2013qz . We also use $x_{i}$ to denote the momentum fraction of light
anti-quark in each meson and set the momentum $p_{i}$ and $k_{i}$ ( the
momenta carried by the spectator quark in $B_{s}$ and $\phi$ meson ) in the
following forms:

$\displaystyle p_{1}$ $\displaystyle=$
$\displaystyle\frac{m_{B_{s}}}{\sqrt{2}}(1,1,0_{\bot}),\quad
p_{2}=\frac{rm_{B_{s}}}{\sqrt{2}}(\eta^{+},\eta^{-},0_{\bot}),$ $\displaystyle
k_{1}$ $\displaystyle=$
$\displaystyle(0,x_{1}\frac{m_{B_{s}}}{\sqrt{2}},k_{1\bot}),\quad
k_{2}=\frac{m_{B_{s}}}{\sqrt{2}}(x_{2}r\eta^{+},x_{2}r\eta^{-},k_{2\bot}).$
(1)

where the mass ratio $r=m_{\phi}/m_{B_{s}}$ , and the factor $\eta^{\pm}$ is
defined in the following form:

$\displaystyle\eta^{\pm}=\eta\pm\sqrt{\eta^{2}-1},\quad{\rm
with}\quad\eta=\frac{1}{2r}\left[1+r^{2}-\frac{q^{2}}{m_{B_{s}}^{2}}\right],$
(2)

where $q=p_{1}-p_{2}$ is the lepton-pair four-momentum. For the final state
$\phi$ meson, its longitudinal and transverse polarization vector
$\epsilon_{L,T}$ can be written in the form of
$\epsilon_{L}=(\eta^{+},-\eta^{-},0_{\bot})/\sqrt{2}$ and
$\epsilon_{T}=(0,0,1)$.

For the $B_{s}$ meson wave function, we use the same kind of parameterizations
as in Refs. Ali:2007ff ; Xiao:2011tx ; Wang:2012ab

$\displaystyle\Phi_{B_{s}}=\frac{i}{\sqrt{2N_{c}}}(p/_{B_{s}}+m_{B_{s}})\gamma_{5}\phi_{B_{s}}({k_{1}}).$
(3)

Here only the contribution of the Lorentz structure $\phi_{B_{s}}(k_{1})$ is
taken into account, since the contribution of the second Lorentz structure
$\bar{\phi}_{B_{s}}$ is numerically small Kurimoto:2001zj ; Lu2003a  and has
been neglected. We adopted the distribution amplitude of the $B_{s}$ meson in
the similar form as that of $B$-meson in the $SU(3)_{f}$ limit being widely
used in the PQCD approach Ali:2007ff ; Xiao:2011tx ; Wang:2012ab

$\displaystyle\phi_{B_{s}}(x,b)$ $\displaystyle=$ $\displaystyle
N_{B_{s}}x^{2}(1-x)^{2}\mathrm{\exp}\left[-\frac{m_{B_{s}}^{2}\;x^{2}}{2\omega_{B_{s}}^{2}}-\frac{1}{2}(\omega_{B_{s}}b)^{2}\right].$
(4)

In order to estimate the theoretical uncertainties induced by the variations
of $\phi_{B_{s}}(x,b)$, one usually take $\omega_{B_{s}}=0.50\pm 0.05$ GeV for
$B_{s}^{0}$ meson li2003 ; Ali:2007ff . The normalization factor $N_{B_{s}}$
depends on the values of the shape parameter $\omega_{B_{s}}$ and the decay
constant $f_{B_{s}}$ and defined through the normalization relation :
$\int_{0}^{1}dx\;\phi_{B_{s}}(x,0)=f_{B_{s}}/(2\sqrt{6})$ Ali:2007ff ;
Wang:2012ab .

For the vector meson $\phi$, the longitudinal and transverse polarization
components can both provide the contribution. Here we adopt the wave functions
of the vector $\phi$ as in Ref. Ali:2007ff :

$\displaystyle\Phi_{\phi}^{||}(p,\epsilon_{L})$ $\displaystyle=$
$\displaystyle\frac{i}{\sqrt{6}}\left[\not\\!\epsilon_{L}m_{\phi}\phi_{\phi}(x)+\not\\!\epsilon_{L}p/\phi^{t}_{\phi}(x)+m_{\phi}\phi^{s}_{\phi}(x)\right],$
(5) $\displaystyle\Phi_{\phi}^{\perp}(p,\epsilon_{T})$ $\displaystyle=$
$\displaystyle\frac{i}{\sqrt{6}}\left[\not\\!\epsilon_{T}m_{\phi}\phi^{v}_{\phi}(x)+\not\\!\epsilon_{T}p/\phi^{T}_{\phi}(x)+m_{\phi}i\epsilon_{\omega\upsilon\rho\sigma}\gamma_{5}\gamma^{\omega}\epsilon^{v}_{T}n^{\rho}v^{\sigma}\phi^{a}_{\phi}(x)\right],$
(6)

where $p$ and $m_{\phi}$ are the momentum and the mass of the $\phi$ meson,
$\epsilon_{L}$ and $\epsilon_{T}$ correspond to the longitudinal and
transverse polarization vectors of the vector meson $\phi$, respectively. The
twist-2 DAs $\phi_{\phi}$ and $\phi_{\phi}^{T}$ in Eqs. (5,6) can be
reconstructed as a Gegenbauer expansion Ali:2007ff :

$\displaystyle\phi_{\phi}(x)$ $\displaystyle=$
$\displaystyle\frac{3f_{\phi}}{\sqrt{6}}x(1-x)\left[1+\sum^{2}_{n=1}a_{n\phi}^{||}C_{n}^{3/2}(t)\right],$
$\displaystyle\phi^{T}_{\phi}(x)$ $\displaystyle=$
$\displaystyle\frac{3f^{T}_{\phi}}{\sqrt{6}}x(1-x)\left[1+\sum^{2}_{n=1}a_{n\phi}^{\perp}C_{n}^{3/2}(t)\right],$
(7)

where $t=2x-1$, $a^{||,\perp}_{1,2}$ are the Gegenbauer moments, while
$C^{3/2}_{1,2}$ are the Gegenbauer polynomials as given in Ref. Ali:2007ff .
$f_{\phi}$ and $f_{\phi}^{T}$ are the longitudinal and transverse components
of the decay constants of the vector meson $\phi$ with $f_{\phi}=0.231\pm
0.004$ GeV and $f^{T}_{\phi}=0.20\pm 0.01$ GeV as given in Ref. Ali:2007ff .
For the relevant Gegenbauer moments we use the same ones as those in Refs.
Ball:2006wn ; Ali:2007ff ; xiao18a ; xiao18b .

$\displaystyle a^{||,\perp}_{1}=0,\quad a_{2\phi}^{||}=0.18\pm 0.08,\quad
a_{2\phi}^{\perp}=0.14\pm 0.07.$ (8)

The twist-3 DAs $\phi^{s,t}_{\phi}$ and $\phi^{v,a}_{\phi}$ in Eqs. (5,6) are
the same ones as those defined in Ref. Ali:2007ff :

$\displaystyle\phi^{t}_{\phi}=\frac{3f^{T}_{\phi}}{2\sqrt{6}}t^{2},\quad\phi^{s}_{\phi}=\frac{3f_{\phi}^{T}}{2\sqrt{6}}(-t)~{},\quad\phi^{v}_{\phi}=\frac{3f_{\phi}}{8\sqrt{6}}(1+t^{2}),\quad\phi^{a}_{\phi}=\frac{3f_{\phi}}{4\sqrt{6}}(-t)~{},$
(9)

where $t=2x-1$.

###  II.2 Effective Hamiltonian for $b\to s\ell^{+}\ell^{-}$ decays

The effective Hamiltonian for the considered semileptonic decay
$B_{s}\to\phi\ell^{+}\ell^{-}$ is defined by the same one as in Refs.
Buchalla:1995vs ; Li:2008tk ; Kindra:2018ayz ; Singh:2019hvj ; Nayek:2018rcq :

$\displaystyle{\cal H}_{\text{eff}}$ $\displaystyle=$
$\displaystyle-\frac{4G_{F}}{\sqrt{2}}\Big{\\{}V_{tb}V^{*}_{ts}\left[C_{1}(\mu)\mathcal{O}_{1}^{c}(\mu)+C_{2}(\mu)\mathcal{O}_{2}^{c}(\mu)+\sum_{i=3}^{10}{C}_{i}(\mu){\mathcal{O}}_{i}(\mu)\right]$
(10)
$\displaystyle+V_{ub}V^{*}_{us}\Big{[}C_{1}(\mu)\left[\mathcal{O}_{1}^{c}(\mu)-\mathcal{O}^{u}_{1}(\mu)\right]+C_{2}(\mu)\left[\mathcal{O}_{2}^{c}(\mu)-\mathcal{O}^{u}_{2}(\mu)\right]\Big{]}\Big{\\}}+{\rm
h.c}.,$

where $G_{F}=1.16638\times 10^{-5}{\rm GeV}^{-2}$ is the Fermi constant,
$V_{ij}$ are the CKM matrix elements. For the operators ${\mathcal{O}}_{i}$ we
adopt those as defined in the so-called $\gamma_{5}$-free basis
Chetyrkin:1996vx ; Chetyrkin:1997gb . Following Ref. Gambino:2003zm , the
operators $\mathcal{O}_{i}$ can be written in the following form:

$\displaystyle{\mathcal{O}}^{c}_{1}$
$\displaystyle=(\bar{s}\gamma_{\mu}T^{a}P_{L}c)(\bar{c}\gamma_{\mu}T^{a}P_{L}b),$
$\displaystyle{\mathcal{O}}^{c}_{2}$
$\displaystyle=(\bar{s}\gamma_{\mu}P_{L}c)(\bar{c}\gamma_{\mu}P_{L}b),$
$\displaystyle{\mathcal{O}}^{u}_{1}$
$\displaystyle=(\bar{s}\gamma_{\mu}T^{a}P_{L}u)(\bar{u}\gamma_{\mu}T^{a}P_{L}b),$
$\displaystyle{\mathcal{O}}^{u}_{2}$
$\displaystyle=(\bar{s}\gamma_{\mu}P_{L}u)(\bar{u}\gamma_{\mu}P_{L}b),$
$\displaystyle{\mathcal{O}}_{3}$
$\displaystyle=(\bar{s}\gamma_{\mu}P_{L}b){\textstyle\sum_{q}}(\bar{q}\gamma_{\mu}q),$
$\displaystyle{\mathcal{O}}_{4}$
$\displaystyle=(\bar{s}\gamma_{\mu}T^{a}P_{L}b){\textstyle\sum_{q}}(\bar{q}\gamma_{\mu}T^{a}q),$
$\displaystyle{\mathcal{O}}_{5}$
$\displaystyle=(\bar{s}\gamma_{\mu}\gamma_{\nu}\gamma_{\rho}P_{L}b){\textstyle\sum_{q}}(\bar{q}\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}q),$
$\displaystyle{\mathcal{O}}_{6}$
$\displaystyle=(\bar{s}\gamma_{\mu}\gamma_{\nu}\gamma_{\rho}T^{a}P_{L}b){\textstyle\sum_{q}}(\bar{q}\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}T^{a}q),$
$\displaystyle{\mathcal{O}}_{7}$
$\displaystyle=\frac{e}{g^{2}}m_{b}\bar{(}s\sigma^{\mu\nu}P_{R}b)F_{\mu\nu},$
$\displaystyle{\mathcal{O}}_{8}$
$\displaystyle=\frac{1}{g}m_{b}\bar{(}s\sigma^{\mu\nu}T^{a}P_{R}b)G^{a}_{\mu\nu},$
$\displaystyle{\mathcal{O}}_{9}$
$\displaystyle=\frac{e^{2}}{g^{2}}(\bar{s}\gamma_{\mu}P_{L}b){\textstyle\sum_{\ell}}(\bar{\ell}\gamma^{\mu}\ell),$
$\displaystyle{\mathcal{O}}_{10}$
$\displaystyle=\frac{e^{2}}{g^{2}}(\bar{s}\gamma_{\mu}P_{L}b){\textstyle\sum_{\ell}}(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell).$
(11)

where ${\mathcal{O}}^{c,u}_{1,2}$ are the current-current operators,
${\mathcal{O}}_{3-6}$ are the QCD penguin operators, ${\mathcal{O}}_{7,8}$ are
the electromagnetic and chromomagnetic penguin operators respectively, and
finally ${\mathcal{O}}_{9,10}$ are the semileptonic operators. The inclusion
of the factors $4\pi/g^{2}=1/{\alpha_{s}}$ in the definition of the operators
${\mathcal{O}}_{7,8,9,10}$ serves to allow a more transparent organisation of
the expansion of the relevant Wilson coefficients as defined in Refs.
Bobeth:1999mk ; Gambino:2003zm  up to next-to-next-to leading order (NNLO).
They are then evolved from the scale $\mu=m_{W}$ down to the scale $\mu=m_{b}$
using the renormalization group equations.

Since the contributions from the subleading chromomagnetic penguin, quark-loop
and annihilation diagrams are highly suppressed for the considered $b\to
s\ell^{+}\ell^{-}$ decays Kindra:2018ayz , we will neglect them in our
calculations. Using the effective Hamiltonian in Eq. (10), the decay amplitude
for $b\to s\ell^{+}\ell^{-}$ loop transition can be decomposed as a product of
a short-distance contributions through Wilson coefficients and long-distance
contribution which is further expressed in terms of form factors,

$\displaystyle{\cal A}(b\to
s\ell^{+}\ell^{-})=\frac{G_{F}}{\sqrt{2}}\frac{\alpha_{\rm{em}}}{\pi}V_{tb}V^{*}_{ts}\bigg{\\{}$
$\displaystyle
C_{9}^{\rm{eff}}(q^{2})[\bar{s}\gamma_{\mu}P_{L}b][\bar{\ell}\gamma^{\mu}\ell]+C_{10}[\bar{s}\gamma_{\mu}P_{L}b][\bar{\ell}\gamma^{\mu}\gamma_{5}\ell]$
(12)
$\displaystyle-2m_{b}C_{7}^{\rm{eff}}\big{[}\bar{s}i\sigma_{\mu\nu}\frac{q^{\nu}}{q^{2}}P_{R}b\big{]}[\bar{\ell}\gamma^{\mu}\ell]\bigg{\\}},$

where $C_{7}^{eff}(\mu)$ and $C_{9}^{eff}(\mu)$ are the effective Wilson
coefficients, defined as in Refs. Chen:2001zc ; Wang:2012ab

$\displaystyle C_{7}^{\rm{eff}}(\mu)$ $\displaystyle=$ $\displaystyle
C_{7}(\mu)+C^{\prime}_{b\to s\gamma}(\mu),$ (13) $\displaystyle
C_{9}^{\rm{eff}}(\mu,q^{2})$ $\displaystyle=$ $\displaystyle
C_{9}(\mu)+Y_{\rm{pert}}(q^{2})+Y_{\rm{res}}(q^{2}).$ (14)

The term $C^{\prime}_{b\to s\gamma}$ in Eq. (13) is the absorptive part of
$b\to s\gamma$ transition and was given in Ref. Chen:2001zc

$\displaystyle C^{\prime}_{b\to
s\gamma}(\mu)=i\alpha_{s}\left\\{\frac{2}{9}\eta^{14/23}\left[\frac{x_{t}\left(x_{t}^{2}-5x_{t}-2\right)}{8\left(x_{t}-1\right)^{3}}+\frac{3x_{t}^{2}\ln
x_{t}}{4(x_{t}-1)^{4}}-0.1687\right]-0.03C_{2}(\mu)\right\\},$ (15)

where $x_{t}=m_{t}^{2}/m_{W}^{2}$ and
$\eta=\alpha_{s}(m_{W})/\alpha_{s}(\mu)$. The explicit expressions of the term
$Y_{\rm{pert}}(q^{2})$ and $Y_{\rm{res}}(q^{2})$ in Eq. (14) are of the
following form Ali:1991is ; Lim:1988yu ; Deshpande:1988bd ; ODonnell:1991cdx ;
Nayek:2018rcq

$\displaystyle Y_{\rm{pert}}(q^{2})$ $\displaystyle=$ $\displaystyle
0.124\,\omega(\hat{s})+g(\hat{m}_{c},\hat{s})C_{0}+\lambda_{u}\left[g(\hat{m}_{c},\hat{s})-g(\hat{m}_{u},\hat{s})\right](3C_{1}+C_{2})$
(16)
$\displaystyle-\frac{1}{2}g(\hat{m}_{b},\hat{s})(C_{3}+3C_{4})-\frac{1}{2}g(\hat{m}_{b},\hat{s})(4C_{3}+4C_{4}+3C_{5}+C_{6})$
$\displaystyle+\frac{2}{9}(3C_{3}+C_{4}+3C_{5}+C_{6}),$ $\displaystyle
Y_{\rm{res}}(q^{2})$ $\displaystyle=$
$\displaystyle-\frac{3\pi}{\alpha_{\rm{em}}^{2}}\Big{[}C_{0}\cdot\\!\sum_{V=J/\Psi,\Psi^{\prime}...}\frac{m_{V}{\cal
B}(V\to
l^{+}l^{-})\Gamma_{\rm{tot}}^{V}}{q^{2}-m_{V}^{2}+im_{V}\Gamma_{\rm{tot}}^{V}}$
(17)
$\displaystyle-\lambda_{u}\;g(\hat{m}_{u},\hat{s})(3C_{1}+C_{2})\cdot\sum_{V=\rho,\omega,\phi}\frac{m_{V}{\cal
B}(V\to
l^{+}l^{-})\Gamma_{\rm{tot}}^{V}}{q^{2}-m_{V}^{2}+im_{V}\Gamma_{\rm{tot}}^{V}}\Big{]},$

where $C_{0}=3C_{1}+C_{2}+3C_{3}+C_{4}+3C_{5}+C_{6}$,
$\hat{s}\\!=\\!q^{2}/m^{2}_{b}$, $\hat{m}_{q}\\!=\\!{m_{q}}/{m_{b}}$ and the
CKM ratio $\lambda_{u}=V_{ub}V^{*}_{us}/(V_{tb}V^{*}_{ts})$. In Eq. (16), the
function $\omega(\hat{s})$ is the soft-gluon correction to the matrix element
of operator $\mathcal{O}_{9}$. The function $g(\hat{m}_{q},\hat{s})$ in Eqs.
(16,17) is related to the basic fermion loop. The contributions from four-
quark operators ${\mathcal{O}}_{1}-{\mathcal{O}}_{6}$ are usually combined
with coefficient $C_{9}$ into an ”effective” one. One can find the explicit
expressions of the function $\omega(\hat{s})$ and $g(\hat{m}_{q},\hat{s})$
easily for example in Ref. Jin:2020jtu  and references therein.

The term $Y_{\rm{pert}}(q^{2})$ in Eqs. (14,16) defines the short distance
perturbative part that involves the indirect contributions from the matrix
element of the four quark operators
$\sum^{10}_{i=1}\langle\ell^{+}\ell^{-}s|\mathcal{O}_{i}|b\rangle$ Ali:1991is
; Lim:1988yu ; Deshpande:1988bd ; ODonnell:1991cdx  and lies at the place far
away from $c\bar{c}$ resonance regions.

The term $Y_{\rm{res}}(q^{2})$ in Eqs. (14,17) describes the long distance
resonant contributions related with the $B_{s}\to\phi V\to\phi(V\to
l^{+}l^{-})$ transitions in the resonance regions, where
$V=(\rho,\omega,\phi,J/\Psi,\psi^{\prime},\cdots)$ are the light vector mesons
and $c\bar{c}$ charmonium states. Up to now the term $Y_{\rm{res}}(q^{2})$ can
not be calculated from the first principle of QCD and may also introduce the
double-counting problem with the term $Y_{\rm{pert}}(q^{2})$. For more details
about such kinds of double-counting problem, one can see the discussions as
given in Refs. Khod2010 ; Khod2013 . In this paper, we checked the possible
effects on the theoretical predictions for the branching ratios and other
considered physical observables by including the term $Y_{\rm{res}}(q^{2})$ or
not in our numerical calculations, and we found that the resulted variations
of the theoretical predictions are less than $5\%$. It is much smaller than
the total theoretical errors: say around $30-40\%$. According to the argument
in Ref. Ahmed:2013lba , the term $Y_{\rm{res}}(q^{2})$ is also generally
small. Because of its smallness and the possible double-counting problem we
here simply drop the term $Y_{\rm{res}}(q^{2})$ out in our numerical
evaluations for all physical observables considered in this paper.

###  II.3  $B_{s}\to\phi$ transition form factors

For the vector meson $\phi$ with polarization vector $\epsilon^{*}$, as usual,
the relevant form factors for $B_{s}\to\phi$ transitions are $V(q^{2})$ and
$A_{0,1,2}(q^{2})$ of the vector and axial-vector currents, and $T_{1,2,3}$ of
the tensor currents. Between the form factors $A_{0,1,2}(q^{2})$ at the point
$q^{2}=0$, there is an exact relation
$2m_{\phi}A_{0}(0)=(m_{B_{s}}+m_{\phi})A_{1}(0)-(m_{B_{s}}-m_{\phi})A_{2}(0)$
in order to avoid the kinematical singularity. Between the form factor
$T_{1,2}$, there also exist a relation $T_{1}(0)=T_{2}(0)$ in an algebraic
manner which is implied by the identity
$\sigma^{\mu\nu}\gamma_{5}=-\frac{i}{2}\epsilon^{\mu\nu\alpha\beta}\sigma_{\alpha\beta}$
with the $\epsilon^{0123}=+1$ convention for the the Levi-Civita tensor.

Using the well-studied wave functions as given in previous subsection, the
PQCD factorization formulas for the relevant form factors of
$B_{s}\to\phi\ell^{+}\ell^{-}$ decays can be calculated and written in the
following form:

$\displaystyle V(q^{2})$ $\displaystyle=$ $\displaystyle 8\pi
m_{B_{s}}^{2}C_{F}(1\\!+\\!r)\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\phi_{B_{s}}(x_{1})$ (18)
$\displaystyle\times\Bigl{\\{}\Big{[}\\!-\\!x_{2}r\phi^{v}_{\phi}(x_{2})\\!+\\!\phi^{T}_{\phi}(x_{2})\\!+\\!\frac{1\\!+\\!x_{2}r\eta}{\sqrt{\eta^{2}\\!-\\!1}}\phi^{a}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{1}(t_{1})$
$\displaystyle\\!+\\!\Big{[}\left(r\\!+\\!\frac{x_{1}}{2\sqrt{\eta^{2}\\!-\\!1}}\right)\phi^{v}_{\phi}(x_{2})\\!-\\!\frac{x_{1}-2r\eta}{2\sqrt{\eta^{2}\\!-\\!1}}\phi^{a}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{2}(t_{2})\Bigr{\\}},$
$\displaystyle A_{0}(q^{2})$ $\displaystyle=$ $\displaystyle 8\pi
m_{B_{s}}^{2}C_{F}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\phi_{B_{s}}(x_{1})\times\Bigl{\\{}\Big{[}\left(1\\!+\\!x_{2}r(2\eta\\!-\\!r)\right)\phi_{\phi}(x_{2})$
(19)
$\displaystyle+\\!(1\\!-\\!2x_{2})r\phi^{t}_{\phi}(x_{2})\\!+\\!\frac{(1\\!-\\!r\eta)-2x_{2}r(\eta\\!-\\!r)}{\sqrt{\eta^{2}\\!-\\!1}}\phi^{s}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{1}(t_{1})$
$\displaystyle\\!+\\!\Big{[}\Big{[}\frac{x_{1}}{\sqrt{\eta^{2}\\!-\\!1}}\left(\frac{\eta\\!+\\!r}{2}\\!-\\!r\eta^{2}\right)\\!+\\!\left(\frac{x_{1}}{2}\\!-\\!x_{1}r\eta\\!+\\!r^{2}\right)\Big{]}\phi_{\phi}(x_{2})$
$\displaystyle-\\!\Big{[}\frac{x_{1}(1\\!-\\!r\eta)+2r(r\\!-\\!\eta)}{\sqrt{\eta^{2}\\!-\\!1}}\\!-\\!x_{1}r\Big{]}\phi^{s}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{2}(t_{2})\Bigr{\\}},$
$\displaystyle A_{1}(q^{2})$ $\displaystyle=$ $\displaystyle 16\pi
m_{B_{s}}^{2}C_{F}\frac{r}{1\\!+\\!r}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\phi_{B_{s}}(x_{1})$ (20)
$\displaystyle\times\Bigl{\\{}\Big{[}(1\\!+\\!x_{2}r\eta)\phi^{v}_{\phi}(x_{2})\\!+\\!(\eta\\!-\\!2x_{2}r)\phi^{T}_{\phi}(x_{2})\\!+\\!x_{2}r\sqrt{\eta^{2}\\!-\\!1}\phi^{a}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{1}(t_{1})$
$\displaystyle\\!+\\!\Big{[}\left(r\eta\\!-\\!\frac{x_{1}}{2}\right)\phi^{v}_{\phi}(x_{2})\\!+\\!\left(r\sqrt{\eta^{2}\\!-\\!1}\\!+\\!\frac{x_{1}}{2}\right)\phi^{a}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{2}(t_{2})\Bigr{\\}},$
$\displaystyle A_{2}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{(1\\!+\\!r)^{2}(\eta\\!-\\!r)}{2r(\eta^{2}\\!-\\!1)}A_{1}(q^{2})-8\pi
m_{B_{s}}^{2}C_{F}\frac{1\\!+\\!r}{\eta^{2}\\!-\\!r}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\phi_{B_{s}}(x_{1})$ (21)
$\displaystyle\times\Bigl{\\{}\Big{[}\big{[}\eta\left(1\\!-\\!x_{2}r^{2}\right)\\!+\\!r\left(x_{2}(2\eta^{2}\\!-\\!1)\\!-\\!1\right)\big{]}\phi_{\phi}(x_{2})\\!+\\!\Big{[}1\\!+\\!2x_{2}r^{2}\\!-\\!(1\\!+\\!2x_{2})r\eta\Big{]}\phi^{t}_{\phi}(x_{2})$
$\displaystyle+r(1\\!-\\!2x_{2})\sqrt{\eta^{2}\\!-\\!1}\phi^{s}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{1}(t_{1})$
$\displaystyle+\\!\Big{[}\Big{[}\left(r\eta\\!-\\!\frac{1}{2}\right)x_{1}\sqrt{\eta^{2}\\!-\\!1}\\!-\\!\Big{[}r\left(r\eta\\!-\\!1\\!-\\!x_{1}\eta^{2}\right)\\!+\\!\frac{x_{1}(r\\!+\\!\eta)}{2}\Big{]}\Big{]}\phi_{\phi}(x_{2})$
$\displaystyle+\\!\Big{[}x_{1}(r\eta\\!-\\!1)\\!+\\!(x_{2}\\!-\\!2)r\sqrt{\eta^{2}\\!-\\!1}\Big{]}\phi^{s}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{2}(t_{2})\Bigr{\\}},$
$\displaystyle T_{1}(q^{2})$ $\displaystyle=$ $\displaystyle 8\pi
m_{B_{s}}^{2}C_{F}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\phi_{B_{s}}(x_{1})\times\Bigl{\\{}\Big{[}(1\\!-\\!2x_{2})r\phi^{v}_{\phi}(x_{2})\\!$
(22)
$\displaystyle+\left(1\\!+\\!2x_{2}r\eta\\!-\\!x_{2}r^{2}\right)\phi^{T}_{\phi}(x_{2})\\!+\\!\frac{1\\!+\\!2x_{2}r^{2}\\!-\\!(1\\!+\\!2x_{2})r\eta}{\sqrt{\eta^{2}\\!-\\!1}}\phi^{a}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{1}(t_{1})$
$\displaystyle+\Big{[}\Big{[}\left(1\\!-\\!\frac{x_{1}}{2}\right)r-\frac{x_{1}(r\eta\\!-\\!1)}{2\sqrt{\eta^{2}\\!-\\!1}}\Big{]}\phi^{v}_{\phi}(x_{2})$
$\displaystyle+\Big{[}\frac{r(\eta\\!-\\!r)}{\sqrt{\eta^{2}\\!-\\!1}}\\!+\\!\frac{x_{1}}{2}\left(r+\frac{r\eta\\!-\\!1}{\sqrt{\eta^{2}\\!-\\!1}}\right)\Big{]}\phi^{a}_{\phi}(x_{2})\Big{]}\cdot
H_{2}(t_{2})\Bigr{\\}},\quad$ $\displaystyle T_{2}(q^{2})$ $\displaystyle=$
$\displaystyle 16\pi m_{B_{s}}^{2}C_{F}\frac{r}{1\\!-\\!r^{2}}\int
dx_{1}dx_{2}\int b_{1}db_{1}b_{2}db_{2}\phi_{B_{s}}(x_{1})$ (23)
$\displaystyle\times\Bigl{\\{}\Big{[}(1\\!-\\!(1+2x_{2})r\eta\\!+\\!2x_{2}r^{2})\phi^{v}_{\phi}(x_{2})\\!$
$\displaystyle+\\!\Big{(}x_{2}r\eta(2\eta\\!-\\!r)\\!-\\!x_{2}r\\!+\\!\eta\\!-\\!r\Big{)}\phi^{T}_{\phi}(x_{2})+\\!(1\\!-\\!2x_{2})r\sqrt{\eta^{2}\\!-\\!1}\phi^{a}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{1}(t_{1})\\!$
$\displaystyle+\\!\Big{[}\Big{[}\frac{x_{2}}{2}\Big{(}1\\!+\\!\frac{\eta}{\sqrt{\eta^{2}\\!-\\!1}}\Big{)}(r\eta\\!-\\!1)\\!+\\!\Big{(}r\\!+\\!\frac{x_{1}}{2\sqrt{\eta^{2}\\!-\\!1}}\Big{)}(\eta\\!-\\!r)\Big{]}\phi^{v}_{\phi}(x_{2})$
$\displaystyle\\!+\\!\Big{[}\Big{(}1\\!-\\!\frac{x_{1}}{2}\Big{)}r\sqrt{\eta^{2}-1}\\!+\\!\frac{x_{1}}{2}(1\\!-\\!r\eta)\Big{]}\phi^{a}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{2}(t_{2})\Bigr{\\}},$
$\displaystyle T_{3}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{(1\\!-\\!r)^{2}(\eta\\!+\\!r)}{2r(\eta^{2}\\!-\\!1)}T_{2}(q^{2})-8\pi
m_{B_{s}}^{2}C_{F}\frac{1\\!-\\!r^{2}}{\eta^{2}\\!-\\!1}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\phi_{B_{s}}(x_{1})$ (24)
$\displaystyle\times\Bigl{\\{}\Big{[}\frac{\eta^{2}\\!-\\!(1\\!+\\!2x_{2})r\eta\\!+\\!2x_{2}r^{2}}{\eta\\!-\\!r}\phi_{\phi}(x_{2})\\!+\\!(1\\!+\\!x_{2}r\eta)\phi^{t}_{\phi}(x_{2})\\!+\\!x_{2}r\sqrt{\eta^{2}\\!-\\!1}\phi^{s}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{1}(t_{1})$
$\displaystyle+\\!\Big{[}\big{[}r\\!-\\!\frac{x_{1}}{2}(\eta\\!+\\!\sqrt{\eta^{2}\\!-\\!1})\big{]}\phi_{\phi}(x_{2})\\!+\\!(x_{1}\\!+\\!2r\sqrt{\eta^{2}\\!-\\!1})\phi^{s}_{\phi}(x_{2})\Big{]}\\!\cdot\\!H_{2}(t_{2})\Bigr{\\}},$

where $r=m_{\phi}/m_{B_{s}}$, the twist-2 DAs $(\phi_{\phi},\phi^{T}_{\phi})$
and the twist-3 DAs $(\phi_{\phi}^{s,t},\phi^{v,a}_{\phi})$ have been defined
in Eqs. (7,9). The function $H_{i}(t_{i})$ in above equations are of the
following form

$\displaystyle
H_{i}(t_{i})=h_{i}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{i})\cdot
S_{t}(x_{2})\exp\left[-S_{ab}(t_{i})\right],\quad for\quad i=(1,2).$ (25)

The hard functions $h_{1,2}(x_{1},x_{2},b_{1},b_{2})$ come form the Fourier
transform of virtual quark and gluon propagators and they can be defined by

$\displaystyle h_{1}$ $\displaystyle=$ $\displaystyle
K_{0}(\beta_{1}b_{1})\left[\theta(b_{1}-b_{2})I_{0}(\alpha_{1}b_{2})K_{0}(\alpha_{1}b_{1})+\theta(b_{2}-b_{1})I_{0}(\alpha_{1}b_{1})K_{0}(\alpha_{1}b_{2})\right],$
$\displaystyle h_{2}$ $\displaystyle=$ $\displaystyle
K_{0}(\beta_{2}b_{1})\left[\theta(b_{1}-b_{2})I_{0}(\alpha_{2}b_{2})K_{0}(\alpha_{2}b_{1})+\theta(b_{2}-b_{1})I_{0}(\alpha_{2}b_{1})K_{0}(\alpha_{2}b_{2})\right],$
(26)

where $K_{0}$ and $I_{0}$ are modified Bessel functions, and

$\displaystyle\alpha_{1}=m_{B_{s}}\sqrt{x_{2}r\eta^{+}},\quad\alpha_{2}=m_{B_{s}}\sqrt{x_{1}r\eta^{+}-r^{2}+r_{s}^{2}},\quad\beta_{1}=\beta_{2}=m_{B_{s}}\sqrt{x_{1}x_{2}r\eta^{+}},\quad$
(27)

where $r=m_{\phi}/m_{B_{s}},r_{s}=m_{s}/m_{B_{s}}$. The hard scales $t_{i}$ in
Eq. (25) are chosen as the largest scale of the virtuality of the internal
particles in the hard $b$-quark decay diagram, including $1/b_{i}(i=1,2)$:

$\displaystyle t_{1}=\max\\{\alpha_{1},1/b_{1},1/b_{2}\\},\quad
t_{2}=\max\\{\alpha_{2},1/b_{1},1/b_{2}\\}.$ (28)

The threshold resummation factor $S_{t}(x)$ in Eq. (25) is adopted from
Kurimoto:2001zj ,

$\displaystyle
S_{t}=\frac{2^{1+2c}\Gamma(3/2+c)}{\sqrt{\pi}\Gamma(1+c)}[x(1-x)]^{c},$ (29)

with a fitted parameter $c(Q^{2})=0.04Q^{2}-0.51Q+1.87$ Li:2009pr  and
$Q^{2}=m^{2}_{B_{s}}(1-r^{2})$ Wang:2015uea . The function $S_{t}(x)$ is
normalized to unity. The function $\exp[-S_{ab}(t)]$ in Eq. (25) contains the
Sudakov logarithmic corrections and the renormalization group evolution
effects of both the wave functions and the hard scattering amplitude, for more
details of function $\exp[-S_{ab}(t)]$ one can see Refs. Lu:2000em ;
Kurimoto:2001zj .

##  III Observables for $B_{s}\to\phi\ell^{+}\ell^{-}$ decays

In experimental analysis, the $\bar{B}_{s}\to\phi\ell^{+}\ell^{-}$ decay is
treated as the four body differential decay distribution
$\bar{B}_{s}\to\phi(\to K^{+}K^{-})\ell^{+}\ell^{-}$, and has been described
in terms of the four kinematic variables Aaij:2015oid ; DescotesGenon:2012zf ;
Descotes-Genon:2013vna ; Bobeth:2008ij : the lepton invariant mass squared
$q^{2}$ and the three decay angles
$\vec{\Omega}=(\cos{\theta_{K}},\cos{\theta_{\ell}},\Phi)$. The angle
${\theta_{K}}$ is the angle between the direction of flight of $K^{+}$ and
$B_{s}$ meson in the rest frame of $\phi$, ${\theta_{\ell}}$ is the angle made
by $\ell^{-}$ with respect to the $B_{s}$ meson in the dilepton rest frame and
$\Phi$ is the azimuthal angle between the two planes formed by dilepton and
$K^{+}K^{-}$.

With the hadronic and leptonic amplitudes defined in Eq. (12), we write down
the four fold differential distribution of four-body $\bar{B}_{s}\to\phi(\to
K^{+}K^{-})\ell^{+}\ell^{-}$ decay Kindra:2018ayz ; Aaij:2015oid ;
Altmannshofer:2008dz ; Becirevic:2011bp ,

$\frac{d^{4}\Gamma}{dq^{2}\,d\vec{\Omega}}=\frac{9}{32\pi}I(q^{2},\vec{\Omega}),\quad
d\vec{\Omega}=d\cos{\theta_{K}}\,d\cos{\theta_{\ell}}\,d\Phi,$ (30)

where the functions $I(q^{2},\vec{\Omega})$ can be written in terms of a set
of angular coefficients and trigonometric functions Altmannshofer:2008dz :

$\displaystyle I(q^{2},\vec{\Omega})$ $\displaystyle=$
$\displaystyle{\textstyle\sum_{i}}I_{i}(q^{2}){f_{i}(\vec{\Omega})}$ (31)
$\displaystyle=$ $\displaystyle
I_{1s}\sin^{2}{\theta_{K}}+I_{1c}\cos^{2}{\theta_{K}}+(I_{2s}\sin^{2}{\theta_{K}}+I_{2c}\cos^{2}{\theta_{K}})\cos
2{\theta_{\ell}}$ $\displaystyle+$ $\displaystyle
I_{3}\sin^{2}{\theta_{K}}\sin^{2}{\theta_{\ell}}\cos 2\Phi+I_{4}\sin
2{\theta_{K}}\sin 2{\theta_{\ell}}\cos\Phi$ $\displaystyle+$ $\displaystyle
I_{5}\sin
2{\theta_{K}}\sin{\theta_{\ell}}\cos\Phi+I_{6s}\sin^{2}{\theta_{K}}\cos{\theta_{\ell}}+I_{7}\sin
2{\theta_{K}}\sin{\theta_{\ell}}\sin\Phi$ $\displaystyle+$ $\displaystyle
I_{8}\sin 2{\theta_{K}}\sin
2{\theta_{\ell}}\sin\Phi+I_{9}\sin^{2}{\theta_{K}}\sin^{2}{\theta_{\ell}}\sin
2\Phi\,.$

For the CP-conjugated mode $B_{s}\to\phi(\to K^{-}K^{+})\ell^{+}\ell^{-}$, the
corresponding expression of the angular decay distribution is

$\displaystyle\frac{d^{4}\bar{\Gamma}}{dq^{2}\,d\vec{\Omega}}=\frac{9}{32\pi}\bar{I}(q^{2},\vec{\Omega})\,,$
(32)

where the function $\bar{I}(q^{2},\vec{\Omega})$ is obtained from
$I(q^{2},\vec{\Omega})$ in Eq. (31) by making the complex conjugation for all
weak phases in $I_{i}$ Altmannshofer:2008dz , and numerically by the following
substitution:

$\displaystyle I_{1(c,s),2(c,s),3,4,7}\to\bar{I}_{1(c,s),2(c,s),3,4,7},\quad
I_{5,6s,8,9}\to-\bar{I}_{5,6s,8,9}.$ (33)

The minus sign in Eq. (33) is a result of the convention that, under the
previous definitions of three angles $(\theta_{K},\theta_{l},\Phi)$, a CP
transformation interchanges the lepton and anti-lepton, i.e., leading to the
transformation ${\theta_{\ell}}\to{\theta_{\ell}}-\pi$ and $\Phi\to-\Phi$.

The angular coefficients $I_{i}$, which are functions of $q^{2}$ only, are
usually expressed in terms of the transverse amplitudes Matias:2012xw ;
Aaij:2015oid . In the limit of massless leptons, there are six such complex
amplitudes: ${\cal A}_{0}^{L,R}$, ${\cal A}_{\|}^{L,R}$ and ${\cal
A}_{\perp}^{L,R}$, where $L$ and $R$ refer to the chirality of the leptonic
current. For the massive case, an additional complex amplitude ${\cal A}_{t}$
is required, where the timelike component of the virtual gauge boson (which
can later decay into dilepton) couple to an axial-vector current.

In Table 1, we summarize the treatment of the angular distribution by
decomposition of the angular coefficients $I_{i}(q^{2})$ into seven transverse
amplitude ${\cal A}^{L,R}_{\perp,\|,0}$ and ${\cal A}_{t}$ as well as the
corresponding trigonometric factor ${f_{i}(\vec{\Omega})}$. Here we will not
consider scalar contribution to facilitate the comparison with Ref.
Jin:2020jtu . Notice that the distribution including lepton masses (but
neglecting scalar $I_{6c}=0$) contains eleven $I_{i}$ where only 10 of them
are independent Egede:2010zc ; Matias:2012xw . In the limit of massless
leptons, it is easy to obtain the relations $I_{1s}=3I_{2s}$ and
$I_{1c}=-I_{2c}$ Altmannshofer:2008dz .

Table 1: The explicit expressions of the angular coefficients $I_{i}(q^{2})$ and $f_{i}(\vec{\Omega})$ appeared in Eq. (31).  $i$ | $I_{i}(q^{2})$ | $f_{i}(\vec{\Omega})$
---|---|---
$1s$ | $(\frac{3}{4}\\!-\\!\hat{m}^{2}_{\ell})\left[|{\cal A}_{\parallel}^{\rm L}|^{2}\\!+\\!|{\cal A}_{\perp}^{\rm L}|^{2}\\!+\\!|{\cal A}_{\parallel}^{\rm R}|^{2}\\!+\\!|{\cal A}_{\perp}^{\rm R}|^{2}\right]\\!+\\!4\hat{m}^{2}_{\ell}{\rm Re}\left[{\cal A}_{\perp}^{\rm L}{\cal A}_{\perp}^{\rm R\ast}\\!+\\!{\cal A}_{\parallel}^{\rm L}{\cal A}_{\parallel}^{\rm R\ast}\right]$ | $\sin^{2}{\theta_{K}}$
$1c$ | $|{\cal A}_{0}^{\rm L}|^{2}+|{\cal A}_{0}^{\rm R}|^{2}+4\hat{m}^{2}_{\ell}\left[|{\cal A}_{t}|^{2}+2{\rm Re}[{\cal A}_{0}^{\rm L}{\cal A}_{0}^{\rm R\ast}]\right]$ | $\cos^{2}{\theta_{K}}$
$2s$ | $\frac{1}{4}\beta^{2}_{\ell}\left[|{\cal A}_{\parallel}^{\rm L}|^{2}+|{\cal A}_{\perp}^{\rm L}|^{2}+|{\cal A}_{\parallel}^{\rm R}|^{2}+|{\cal A}_{\perp}^{\rm R}|^{2}\right]$ | $\sin^{2}{\theta_{K}}\cos 2{\theta_{\ell}}$
$2c$ | $-\beta^{2}_{\ell}\left[|{\cal A}_{0}^{\rm L}|^{2}+|{\cal A}_{0}^{\rm R}|^{2}\right]$ | $\cos^{2}{\theta_{K}}\cos 2{\theta_{\ell}}$
3 | $\frac{1}{2}\beta^{2}_{\ell}\left[|{\cal A}_{\perp}^{\rm L}|^{2}-|{\cal A}_{\parallel}^{\rm L}|^{2}+|{\cal A}_{\perp}^{\rm R}|^{2}-|{\cal A}_{\parallel}^{\rm R}|^{2}\right]$ | $\sin^{2}{\theta_{K}}\sin^{2}\theta_{\ell}\cos 2\Phi$
4 | $\sqrt{\frac{1}{2}}\beta^{2}_{\ell}{\rm Re}({\cal A}_{0}^{\rm L}{\cal A}_{\parallel}^{{\rm L}\ast}+{\cal A}_{0}^{\rm R}{\cal A}_{\parallel}^{{\rm R}\ast})$ | $\sin 2{\theta_{K}}\sin 2{\theta_{\ell}}\cos\Phi$
5 | $\sqrt{2}\beta_{\ell}{\rm Re}({\cal A}_{0}^{\rm L}{\cal A}_{\perp}^{{\rm L}\ast}-{\cal A}_{0}^{\rm R}{\cal A}_{\perp}^{{\rm R}\ast})$ | $\sin 2{\theta_{K}}\sin{\theta_{\ell}}\cos\Phi$
$6s$ | $2\beta_{\ell}{\rm Re}({\cal A}_{\parallel}^{\rm L}{\cal A}_{\perp}^{{\rm L}\ast}-{\cal A}_{\parallel}^{\rm R}{\cal A}_{\perp}^{{\rm R}\ast})$ | $\sin^{2}{\theta_{K}}\cos{\theta_{\ell}}$
7 | $\sqrt{2}\beta_{\ell}{\rm Im}({\cal A}_{0}^{\rm L}{\cal A}_{\parallel}^{{\rm L}\ast}-{\cal A}_{0}^{\rm R}{\cal A}_{\parallel}^{{\rm R}\ast})$ | $\sin 2{\theta_{K}}\sin{\theta_{\ell}}\sin\Phi$
8 | $\sqrt{\frac{1}{2}}\beta^{2}_{\ell}{\rm Im}({\cal A}_{0}^{\rm L}{\cal A}_{\perp}^{{\rm L}\ast}+{\cal A}_{0}^{\rm R}{\cal A}_{\perp}^{{\rm R}\ast})$ | $\sin 2{\theta_{K}}\sin 2{\theta_{\ell}}\sin\Phi$
9 | $\beta^{2}_{\ell}{\rm Im}({\cal A}_{\parallel}^{{\rm L}\ast}{\cal A}_{\perp}^{\rm L}+{\cal A}_{\parallel}^{{\rm R}\ast}{\cal A}_{\perp}^{\rm R})$ | $\sin^{2}{\theta_{K}}\sin^{2}{\theta_{\ell}}\sin 2\Phi$

The seven transverse amplitudes ${\cal A}_{0}^{L,R}$, ${\cal A}_{\|}^{L,R}$,
${\cal A}_{\perp}^{L,R}$ and ${\cal A}_{t}$ of $B_{s}\to\phi\ell^{+}\ell^{-}$
decay, in turn, can be parameterized by means of the relevant form factors
Egede:2008uy ; Altmannshofer:2008dz :

$\displaystyle{\cal A}_{\perp}^{\rm L,R}$ $\displaystyle=$ $\displaystyle-
N_{\ell}\sqrt{2N_{\phi}}\sqrt{\lambda}\left[(C_{9}^{\rm{eff}}\mp
C_{10})\frac{V(q^{2})}{m_{B_{s}}+m_{\phi}}+2\hat{m}_{b}C_{7}^{\rm{eff}}T_{1}(q^{2})\right],$
(34) $\displaystyle{\cal A}_{\parallel}^{\rm L,R}$ $\displaystyle=$
$\displaystyle N_{\ell}\sqrt{2N_{\phi}}\Big{[}(C_{9}^{\rm{eff}}\mp
C_{10})(m_{B_{s}}+m_{\phi})A_{1}(q^{2})+2\hat{m}_{b}C_{7}^{\rm{eff}}(m^{2}_{B_{s}}-m^{2}_{\phi})T_{2}(q^{2})\Big{]},$
(35) $\displaystyle{\cal A}_{0}^{\rm L,R}$ $\displaystyle=$
$\displaystyle\frac{N_{\ell}\sqrt{N_{\phi}}}{2m_{\phi}\sqrt{q^{2}}}\Big{\\{}(C_{9}^{\rm{eff}}\\!\mp\\!C_{10})\left[(m^{2}_{B_{s}}\\!-\\!m^{2}_{\phi}\\!-\\!q^{2})(m_{B_{s}}\\!+\\!m_{\phi})A_{1}(q^{2})\\!-\\!\frac{\lambda}{m_{B_{s}}\\!+\\!m_{\phi}}A_{2}(q^{2})\right]$
(36)
$\displaystyle+2m_{b}C_{7}^{\rm{eff}}\Big{[}(m^{2}_{B_{s}}\\!+\\!3m^{2}_{\phi}\\!-\\!q^{2})T_{2}(q^{2})\\!-\\!\frac{\lambda}{m^{2}_{B_{s}}\\!-\\!m^{2}_{\phi}}T_{3}(q^{2})\Big{]}\Big{\\}},$
$\displaystyle{\cal A}_{t}$ $\displaystyle=$ $\displaystyle
2N_{\ell}\sqrt{N_{\phi}}\frac{\sqrt{\lambda}}{\sqrt{q^{2}}}C_{10}A_{0}(q^{2}),$
(37)

where
$\lambda=(m^{2}_{B_{s}}\\!-\\!m^{2}_{\phi}\\!-\\!q^{2})^{2}\\!-\\!4m^{2}_{\phi}q^{2}$,
$\hat{m}_{b}=m_{b}/q^{2}$ and the normalization constants are given as:

$\displaystyle
N_{\ell}=\frac{i\alpha_{em}G_{F}}{4\sqrt{2}\pi}V_{tb}V^{\ast}_{ts},\quad
N_{\phi}=\frac{8\sqrt{\lambda}q^{2}}{3\times
256\pi^{3}m^{3}_{B_{s}}}\sqrt{1\\!-\frac{4m_{\ell}^{2}}{q^{2}}}{\cal
B}(\phi\to K^{+}K^{-}).$ (38)

In numerical calculations, we take ${\cal B}(\phi\to K^{+}K^{-})=0.492$ from
PDG 2018 pdg2018 . It is easy to see that the narrow width approximation works
well in the case of $\phi$ meson since $\Gamma_{\phi}/{m_{\phi}}=4.17\times
10^{-3}\sim 0$.

Analogous to Ref. Altmannshofer:2008dz , to separate CP-conserving and CP-
violating effects, one can define the CP averaged angular coefficients $S_{i}$
and CP asymmetry angular coefficients $A_{i}$ normalized by the differential
(CP-averaged) decay rate to reduce the theoretical uncertainties,

$\displaystyle
S_{i}=\frac{I_{i}+\bar{I}_{i}}{d(\Gamma+\bar{\Gamma})/dq^{2}},\quad
A_{i}=\frac{I_{i}-\bar{I}_{i}}{d(\Gamma+\bar{\Gamma})/dq^{2}},$ (39)

where $I_{i}$ and $\bar{I}_{i}$ have been defined in Eqs. (31,32,33) and Table
1, and the differential decay rate reads (analogously for $\bar{\Gamma}$),

$\displaystyle{\frac{d\Gamma}{dq^{2}}}={\frac{1}{4}}\left(3I_{1c}+6I_{1s}-I_{2c}-2I_{2s}\right)\,.$
(40)

Based on the definition of $S_{i}$ , one can find the relation
$3S_{1c}+6S_{1s}-S_{2c}-2S_{2s}=4$. Consequently, all established observables
can be expressed in terms of $S_{i}$ and $A_{i}$:

  * (1)

the CP asymmetry

$\displaystyle A_{\rm
CP}(q^{2})=\frac{d\Gamma/dq^{2}-d\bar{\Gamma}/dq^{2}}{d\Gamma/dq^{2}+d\bar{\Gamma}/dq^{2}}={\frac{1}{4}}\left(3A_{1c}+6A_{1s}-A_{2c}-2A_{2s}\right).$
(41)

  * (2)

The lepton forward-backward (CP) asymmetry:

$\displaystyle{A}_{\mathrm{FB}}(q^{2})=\frac{\left[\int_{0}^{1}-\int_{-1}^{0}\right]d\cos\theta_{\ell}\frac{d^{2}(\Gamma-\bar{\Gamma})}{dq^{2}d\cos\theta_{\ell}}}{d(\Gamma+\bar{\Gamma})/dq^{2}}=\frac{3}{4}S_{6s},$
(42) $\displaystyle{A}^{\rm
CP}_{\mathrm{FB}}(q^{2})=\frac{\left[\int_{0}^{1}-\int_{-1}^{0}\right]d\cos\theta_{\ell}\frac{d^{2}(\Gamma+\bar{\Gamma})}{dq^{2}d\cos\theta_{\ell}}}{d(\Gamma+\bar{\Gamma})/dq^{2}}=\frac{3}{4}A_{6s}.$
(43)

  * (3)

The $\phi$ polarization fractions:

$\displaystyle F_{L}(q^{2})=\frac{1}{4}(3S_{1c}-S_{2c}),\quad
F_{T}(q^{2})=\frac{1}{2}(3S_{1s}-S_{2s}).$ (44)

In the massless limit, since the CP-averaged observable $S_{1(c,s),2(c,s)}$
obey the relations $S_{1s}=3S_{2s}$ and $S_{1c}=-S_{2c}$, the definitions of
the polarization fractions can be simplified directly as:

$\displaystyle F_{L}(q^{2})=S_{1c}=-S_{2c},\quad
F_{T}(q^{2})=\frac{4}{3}S_{1s}=4S_{2s}.$ (45)

  * (4)

The clean (no S-wave pollution) observables $P_{1,2,3}$ and
$P^{\prime}_{4,5,6}$ in the natural basis can be defined in terms of the
coefficients $S_{i}$ through the following relations Kruger:2005ep ;
Becirevic:2011bp ; Wang2015 :

$\displaystyle P_{1}$ $\displaystyle=$
$\displaystyle\frac{S_{3}}{2S_{2s}},\quad
P_{2}=\beta_{\ell}\frac{S_{6s}}{8S_{2s}},\quad P_{3}=-\frac{S_{9}}{4S_{2s}}\
,$ (46) $\displaystyle P^{\prime}_{4}$ $\displaystyle=$
$\displaystyle\frac{S_{4}}{\sqrt{S_{1c}S_{2s}}},\quad
P^{\prime}_{5}=\frac{\beta_{\ell}S_{5}}{2\sqrt{S_{1c}S_{2s}}},\quad
P^{\prime}_{6}=-\frac{\beta_{\ell}S_{7}}{2\sqrt{S_{1c}S_{2s}}},\quad
P^{\prime}_{8}=-\frac{S_{8}}{\sqrt{S_{1c}S_{2s}}},$ (47)

where $\beta_{\ell}=\sqrt{1-4m_{\ell}^{2}/q^{2}}$.

  * (5)

In the massless limit of leptons, the optimised observables $P^{(\prime)}_{i}$
Matias:2012xw  can be transformed as the following form:

$\displaystyle P_{1}$ $\displaystyle=$ $\displaystyle\frac{2S_{3}}{F_{T}}\quad
P_{2}=\frac{S_{6s}}{2F_{T}},\quad P_{3}=\frac{-S_{9}}{F_{T}},$ (48)
$\displaystyle P^{\prime}_{4}$ $\displaystyle=$
$\displaystyle\frac{2S_{4}}{\sqrt{F_{L}(1-F_{L})}},\quad
P^{\prime}_{5}=\frac{S_{5}}{\sqrt{F_{L}(1-F_{L})}},$ $\displaystyle
P^{\prime}_{6}$ $\displaystyle=$
$\displaystyle-\frac{S_{7}}{\sqrt{F_{L}(1-F_{L})}},\quad
P^{\prime}_{8}=-\frac{2S_{8}}{\sqrt{F_{L}(1-F_{L})}}.$ (49)

One should know that our definitions of the CP averaged angular coefficients
$S_{i}$ , the CP asymmetry angular coefficients $A_{i}$ and the clean
observable $P_{1,2,3}$ and$P^{\prime}_{4,5,6}$ differ from those adopted by
the LHCb collaboration. To be specific, the reasons are the following:

  * (1)

Our conventions for the angles to define the $B_{s}\to\phi\ell^{+}\ell^{-}$
kinematics are identical to the Ref. Altmannshofer:2008dz  but different from
the LHCb choices Aaij:2015oid ; Aaij:2015esa . The corresponding relations are
the following :

$\displaystyle{\theta_{K}}^{\rm LHCb}={\theta_{K}},\quad{\theta_{\ell}}^{\rm
LHCb}=\pi-{\theta_{\ell}},\quad\Phi^{\rm LHCb}=-\Phi.$ (50)

Some angular coefficients $I_{i}$ , $S_{i}$ and $A_{i}$, consequently, will
have different signs:

$\displaystyle I^{\rm LHCb}_{4,6,7,9}=-I_{4,6s,7,9},\quad S^{\rm
LHCb}_{4,6,7,9}=-S_{4,6s,7,9},\quad A^{\rm LHCb}_{4,6,7,9}=-A_{4,6s,7,9}.$
(51)

Other remaining coefficients $I_{i}$ ($S_{i}$ and $A_{i}$), however, have the
same sign in both conventions.

  * (2)

Our definitions of the clean observables $P_{1,2,3}$ and
$P^{\prime}_{4,5,6,8}$ in Eq. (47) in terms of $S_{i}$ may be different from
those defined and used by the LHCb Collaboration for example in Ref.
Aaij:2015oid . The resultant differences of the sign and normalization are of
the following:

$\displaystyle P^{\rm LHCb}_{1}=P_{1},\quad P^{\rm LHCb}_{2,3}=-P_{2,3},\quad
P^{\prime\rm LHCb}_{4,8}=-\frac{1}{2}P^{\prime}_{4,8},\quad P^{\prime\rm
LHCb}_{5,6}=P^{\prime}_{5,6}.$ (52)

For more details about the angular conventions of the angular observables of
the semileptonic decays $B_{(s)}\to Vl^{+}l^{-}$, one can see Ref.
prd93:2016ghz .

##  IV Numerical results and discussions

In the numerical calculations we use the following input parameters (here
masses and decay constants are in units of GeV) pdg2018 ; Ali:2007ff :

$\displaystyle\Lambda^{f=4}_{\overline{\rm MS}}$ $\displaystyle=$
$\displaystyle 0.250,\;\quad\tau_{B_{s}^{0}}=1.509{\rm ps},\;\quad
m_{b}=4.8,\;\quad m_{W}=80.38,\;\quad m_{\phi}=1.019$ $\displaystyle
m_{B_{s}}$ $\displaystyle=$ $\displaystyle 5.367,\;\quad
m_{e}=0.000511,\;\quad m_{\mu}=0.105,\;\quad m_{\tau}=1.777,$ $\displaystyle
f_{B_{s}}$ $\displaystyle=$ $\displaystyle 0.23,\quad f_{\phi}=0.231(4),\quad
f^{\perp}_{\phi}=0.20(1),\quad a^{||}_{2\phi}=0.18(8),\quad
a^{\perp}_{2\phi}=0.14(7).$ (53)

For the CKM matrix elements and angles, we adopt the following values as given
in Ref. pdg2018 :

$\displaystyle V_{tb}$ $\displaystyle=$ $\displaystyle 1.019(25),\quad
V_{us}=0.2243(5),\quad|V_{ts}|=(39.4\pm 2.3)\times 10^{-3},$
$\displaystyle|V_{ub}|$ $\displaystyle=$ $\displaystyle(3.94\pm 0.36)\times
10^{-3},\quad
2\beta_{s}=0.021(31),\quad\gamma=\left({73.5}^{+4.2}_{-5.1}\right)^{\circ}.$
(54)

###  IV.1 The form factors

For the considered semileptonic decays, the differential decay rates and other
physical observables strongly rely on the value and the shape of the relevant
form factors $V(q^{2}),$ $A_{0,1,2}(q^{2})$ and $T_{1,2,3}(q^{2})$ for
$B_{s}\to\phi\ell^{+}\ell^{-}$ decays. These form factors have been calculated
in rather different theories or models Straub:2015ica ; Keum:2000wi ;
Lu:2000em ; li2003 ; Ali:2007ff . Since the PQCD predictions for the
considered form factors are valid only at the large hadronic recoil (
low-$q^{\rm 2}$) region, we usually calculate explicitly the values of the
relevant form factors at the low-$q^{2}$ region, say $0\leq q^{\rm 2}\leq
m_{\rm\tau}^{\rm 2}$, and then make an extrapolation for all relevant form
factors from the low-$q^{\rm 2}$ towards the high-$q^{\rm 2}$ region by using
the pole model parametrizationcheng2004 ; wang2009  or other different
methods.

In Refs. Fan:2015 ; Hu:2019 ; Hu:2020 , we developed a new method: the so-
called “PQCD+Lattice” approach. Here we still use the PQCD approach to
evaluate the form factors at the low $q^{2}$ region, but take those currently
available lattice QCD results for the relevant form factors at the
high-$q^{2}$ region as the lattice QCD input to improve the extrapolation of
the form factors up to $q^{2}_{max}$. In Refs. Hu:2019 ; Hu:2020 , we used the
Bourrely-Caprini-Lellouch (BCL) parametrization method bcl09 ; jhep1905-094
instead of the traditional pole model parametrization.

Table 2: The PQCD predictions of the seven form factors $F_{B_{s}\to\phi}^{i}(0)$ with the theoretical uncertainties from the variations of the parameters $\omega_{B_{s}}$, $a^{||,\perp}_{2\phi}$, $f_{\phi}$ and$f^{T}_{\phi}$.  $F_{B_{s}\to\phi}^{i}(0)$ | PQCD predictions
---|---
$V(0)$ |  ${0.311}$ ${}^{+0.063}_{-0.051}(\omega_{B_{s}})$ ${}^{+0.008}_{-0.007}(a^{\perp}_{2\phi})$ $\pm 0.004(f_{\phi})$ $\pm 0.004(f^{T}_{\phi})$
$A_{0}(0)$ |  ${0.262}$ ${}^{+0.046}_{-0.038}(\omega_{B_{s}})$ ${}^{+0.009}_{-0.008}(a^{||}_{2\phi})$ $\pm 0.002(f_{\phi})$ $\pm 0.008(f^{T}_{\phi})$
$A_{1}(0)$ |  ${0.247}$ ${}^{+0.051}_{-0.041}(\omega_{B_{s}})$ ${}^{+0.005}_{-0.006}(a^{\perp}_{2\phi})$ $\pm 0.003(f_{\phi})$ $\pm 0.003(f^{T}_{\phi})$
$A_{2}(0)$ |  ${0.239}$ ${}^{+0.054}_{-0.042}(\omega_{B_{s}})$ $\pm 0.009(a^{||,\perp}_{2\phi})$ $\pm 0.004(f_{\phi})$ $\pm 0.001(f^{T}_{\phi})$
$T_{1,2}(0)$ |  ${0.264}$ ${}^{+0.052}_{-0.042}(\omega_{B_{s}})$ $\pm 0.006(a^{\perp}_{2\phi})$ $\pm 0.003(f_{\phi})$ $\pm 0.005(f^{T}_{\phi})$
$T_{3}(0)$ |  ${0.196}$ ${}^{+0.041}_{-0.033}(\omega_{B_{s}})$ $\pm 0.008(a^{||,\perp}_{2\phi})$ $\pm 0.003(f_{\phi})$ $\pm 0.001(f^{T}_{\phi})$
Table 3: The central values of the theoretical predictions for the form factors of $B_{s}\to\phi$ transitions at $q^{2}=0$ from different theories or models Ali:2007ff ; Li:2009tx ; Ball:2004rg ; Yilmaz:2008pa ; Straub:2015ica ; Peng:2020ivn ; Melikhov:2000yu ; Faustov:2013pca ; Lu:2007sg ; Wu:2006rd ; Su:2011eq ; Dubnicka:2016nyy  .  | $V(0)$ | $A_{0}(0)$ | $A_{1}(0)$ | $A_{2}(0)$ | $T_{1,2}(0)$ | $T_{3}(0)$
---|---|---|---|---|---|---
This work | $0.311$ | $0.262$ | $0.247$ | $0.239$ | $0.264$ | $0.196$
PQCDAli:2007ff  | $0.25$ | $0.30$ | $0.19$ | $-$ | $-$ |
PQCDLi:2009tx  | 0.26 | $0.31$ | $0.18$ | $0.12$ | $0.23$ | $0.19$
LCSRStraub:2015ica  | $0.387$ | $0.389$ | $0.296$ | $-$ | $0.309$ | $-$
LCSRBall:2004rg  | $0.434$ | $0.474$ | $0.311$ | $0.234$ | $0.349$ | $0.175$
LCSRYilmaz:2008pa  | $0.433$ | $0.382$ | $0.296$ | $0.255$ | $0.348$ | $0.254$
QCDSRPeng:2020ivn  | $0.45$ | $0.30$ | $0.32$ | $0.30$ | $0.33$ | $0.22$
RDAMelikhov:2000yu  | $0.44$ | $0.42$ | $0.34$ | $0.31$ | $0.38$ | $0.26$
RQMFaustov:2013pca  | $0.406$ | $0.322$ | $0.320$ | $0.318$ | $0.275$ | $0.133$
SCETLu:2007sg  | $0.329$ | $0.279$ | $0.232$ | $0.210$ | $0.276$ | $0.170$
HQEFTWu:2006rd  | $0.339$ | $0.269$ | $0.271$ | $0.212$ | $0.299$ | $0.191$
SQEHSu:2011eq  | $0.259$ | $0.311$ | $0.194$ | $-$ | $-$ | $-$
CQMDubnicka:2016nyy  | $0.31$ | $0.28$ | $0.27$ | $0.27$ | $0.27$ | $0.18$

In Table 2, we list the PQCD predictions for all seven relevant form factors
$V(0)$, $A_{0,1,2}(0)$ and $T_{1,2,3}(0)$ for $B_{s}\to\phi$ transition at the
point $q^{2}=0$. The dominant theoretical errors come from the uncertainties
of the parameter $\omega_{B_{s}}=0.50\pm 0.05$ Ali:2007ff  , the Gegenbauer
moments $a_{2\phi}^{||}=0.18\pm 0.08$ and $a_{2\phi}^{\perp}=0.14\pm 0.07$
Ball:2006wn ; Ali:2007ff ; xiao18a ; xiao18b , as well as the decay constants
$f_{\phi}=0.231\pm 0.004$ GeV and $f^{T}_{\phi}=0.20\pm 0.01$ GeV Ali:2007ff .

In Table 3, as a comparison, we also list the central values of the
theoretical predictions for the form factors $F_{B_{s}\to\phi}^{i}(0)$ at
$q^{2}=0$ evaluated in the PQCD approach Ali:2007ff ; Li:2009tx , and in other
different theories or models Ball:2004rg ; Yilmaz:2008pa ; Straub:2015ica ;
Peng:2020ivn ; Melikhov:2000yu ; Faustov:2013pca ; Lu:2007sg ; Wu:2006rd ;
Su:2011eq ; Dubnicka:2016nyy . One can see that there exist always some
differences between different authors, even among the authors using the same
approach. Taking the calculations based on the LCSR method as an example, the
authors of Ref. Straub:2015ica  introduced the hadronic input parameters, Ball
and Zwicky considered the one-loop radiative corrections Ball:2004rg , Yilmaz
included the radiative and higher twist corrections and SU(3) breaking effects
Yilmaz:2008pa .

In Table 4, for the $B_{s}\to\phi$ transition form factors
$(V,A_{0,1,2},T_{1,2,3})$, we quote directly the values of the lattice QCD
results at two or three reference points of the high $q^{2}$ region, say
$q^{2}=12,16$ GeV2 and $q^{2}_{max}=(m_{B_{s}}-m_{\phi})^{2}\approx 18.9$
GeV2, as listed in Table XXXI of Ref. Horgan:2013hoa . In Ref. Horgan:2013hoa
, the authors defined the helicity form factors $A_{12}(q^{2})$ and
$T_{23}(q^{2})$ from the ordinary form factors $A_{1,2}(q^{2})$ and
$T_{2,3}(q^{2})$:

$\displaystyle A_{12}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{(m_{B_{s}}+m_{\phi})^{2}(m^{2}_{B_{s}}-m^{2}_{\phi}-q^{2})A_{1}(q^{2})-\lambda(q^{2})A_{2}(q^{2})}{16m_{B_{s}}m^{2}_{\phi}(m_{B_{s}}+m_{\phi})},$
$\displaystyle T_{23}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{m_{B_{s}}+m_{\phi}}{8m_{B_{s}}m^{2}_{\phi}}\left[(m^{2}_{B_{s}}+3m^{2}_{\phi}-q^{2})T_{2}(q^{2})-\frac{\lambda(q^{2})T_{3}(q^{2})}{m^{2}_{B_{s}}-m^{2}_{\phi}}\right],$
(55)

where the kinematic variable $\lambda(q^{2})=(t_{+}-q^{2})(t_{-}-q^{2})$ with
$t_{\pm}=(m_{B_{s}}\pm m_{\phi})^{2}$. From above two equations and the
numerical values of $(A_{1}(q^{2}),T_{2}(q^{2}),A_{12}(q^{2}),T_{23}(q^{2}))$
as given in Ref. Horgan:2013hoa , we can find the corresponding lattice QCD
results of $A_{2}(q^{2})$ and $T_{3}(q^{2})$ at the two points $q^{2}=(12,16)$
GeV2 by direct numerical calculations. When $q^{2}\to
q^{2}_{max}=(m_{B_{s}}-m_{\phi})^{2}$, however, the parameter
$\lambda(q^{2}_{max})$ in Eq. (55) is also approaching zero simultaneously,
one therefore can not determine $A_{2}(q^{2}_{max})$ and $T_{3}(q^{2}_{max})$
reliably from the values of $A_{12}(q^{2}_{max})$ and $T_{23}(q^{2}_{max})$ as
given in Ref. Horgan:2013hoa . Consequently, $A_{2}(q^{2}_{max})$ and
$T_{3}(q^{2}_{max})$ are absent in Table 4.

Table 4: The values for the lattice QCD results of the relevant $B_{s}\to\phi$ transition form factors at two or three reference points Horgan:2013hoa .  $q^{2}$ | $V(q^{2})$ | $A_{0}(q^{2})$ | $A_{1}(q^{2})$ | $A_{2}(q^{2})$ | $T_{1}(q^{2})$ | $T_{2}(q^{2})$ | $T_{3}(q^{2})$
---|---|---|---|---|---|---|---
12 | $0.77(6)$ | $0.90(6)$ | $0.44(3)$ | $0.48(4)$ | $0.69(4)$ | $0.45(3)$ | $0.46(4)$
16 | $1.19(7)$ | $1.32(7)$ | $0.52(3)$ | $0.54(4)$ | $0.99(5)$ | $0.53(3)$ | $0.70(5)$
18.9 | $1.74(10)$ | $1.85(10)$ | $0.62(3)$ | $-$ | $1.36(8)$ | $0.62(3)$ | $-$

In this work, we will use both the PQCD factorization approach and the
“PQCD+Lattice” approach to evaluate all relevant form factors over the whole
range of $q^{2}$.

  1. (1)

In the PQCD approach, we use the definitions and formulae to calculate the
values of all relevant form factors $V(q^{2})$ , $A_{0,1,2}(q^{2})$ and
$T_{1,2,3}(q^{2})$ at some points in the region of $0\leq q^{2}\leq
m_{\tau}^{2}$. We then make the extrapolation for these form factors to the
large $q^{2}$ region up to $q_{max}^{2}$ by using the selected parametrization
method directly.

  2. (2)

In the “PQCD+Lattice ” approach, we take the lattice QCD results of the form
factors at the points of $q^{2}=(12,16,18.9)$ GeV2 as the new input in the
high $q^{2}$ region and then make a combined fit of the PQCD results in low
$q^{2}$ region and the lattice QCD results in high $q^{2}$ region to determine
the relevant parameters $b_{k}^{i}$ in the $z$-series expansion and then
complete the extrapolation.

  3. (3)

For both approaches, we always use the model-independent $z$-series
parametrization, which is based on a rapidly converging series in the
parameter $z$, as in Refs. Kindra:2018ayz ; Straub:2015ica  to make the
extrapolation. The entire cut $q^{2}$-plane will be mapped onto the unit disc
$|z(q^{2})|\leq 1$ under the conformal transformation as Lu:2018cfc

$\displaystyle
z(q^{2})=\frac{\sqrt{t_{+}-q^{2}}-\sqrt{t_{+}-t_{0}}}{\sqrt{t_{+}-q^{2}}+\sqrt{t_{+}-t_{0}}}$
(56)

where $t_{\pm}=(m_{B_{s}}\pm m_{\phi})^{2}$ and $0\leq t_{0}<t_{-}$ is a
auxiliary parameter which can be optimised to reduce the maximum value of
$|z(q^{2})|$ in the physical range of the form factors and will be taken in
the same way as in Ref. Bharucha:2010im :
$t_{0}=t_{+}(1-\sqrt{1-t_{-}/t_{+}})$. The form factors are finally
parameterized in the BCL version of the $z$-series expansion bcl09

$\displaystyle F_{B_{s}\to\phi}^{i}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{F_{B_{s}\to\phi}^{i}(0)}{1-q^{2}/m_{i,\,\rm
pole}^{2}}\,\left\\{1+\sum_{k=1}^{N}\,b_{k}^{i}\,\left[z(q^{2},t_{0})^{k}-z(0,t_{0})^{k}\right]\right\\}$
(57) $\displaystyle=$
$\displaystyle\frac{F_{B_{s}\to\phi}^{i}(0)}{1-q^{2}/m_{i,\,\rm
pole}^{2}}\,\left\\{1+b_{1}^{i}\,\left[z(q^{2},t_{0})-z(0,t_{0})\right]\right\\}+\cdots.$

Since the term $|z(q^{2},t_{0})|^{2}\leq 0.04$ in the whole considered $q^{2}$
region, the high order $N\geq 2$ terms in Eq. (57) should be very small in
magnitude and therefore can be neglected. After the truncation at $N=1$, the
coefficient $b^{i}_{1}$ for the corresponding form factor
$F_{B_{s}\to\phi}^{i}(q^{2})$ can be determined by fitting to the PQCD
predictions at low $q^{2}$ region and the lattice QCD results in the high
$q^{2}$ region. Taking the form factor $A_{2}(q^{2})$ as an example, we
calculate $A_{2}(q^{2})$ first by employing the PQCD approach in the sixteen
points in the low $0\leq q^{2}\leq m_{\tau}^{2}$ region, take the lattice QCD
results $A_{2}(12)=0.48\pm 0.04$ and $A_{2}(16)=0.54\pm 0.04$ as additional
input, and finally make the fitting for the parameter $b_{1}^{A_{2}}$ and find
that:

$\displaystyle b_{1}^{A_{2}}=-1.820\pm 0.148{\rm(standard\ \ error)},$ (58)

with the goodness-of-fit $R^{2}=0.9978$. For other form factors $F^{i}(q^{2})$
we find the results by following the same kind of procedure. The input values
of the various $\bar{s}b$-resonance mass $m_{i,pole}$ in Eq. (57) can be found
from Ref. pdg2018  and are collected in Table 5. For further discussions on
the systematic uncertainties due to the dependence of truncation schemes and
on the implementation of the strong unitary constraints one can see Refs.
Bharucha:2010im ; Bigi:2016mdz .

Table 5: The masses $m_{i,{\rm pole}}$ in Eq. (57) for the form factors $F^{i}_{B_{s}\to\phi}(q^{2})$ Straub:2015ica .  $F_{B_{s}\to\phi}^{i}(q^{2})$ | $\bar{s}b(J^{P})$ |  $m_{i,{\rm pole}}$(GeV)
---|---|---
${A}_{0}(q^{2})$ | $B_{s}(0^{-})$ | 5.366
${V}(q^{2})$, ${T}_{1}(q^{2})$ | $B_{s}^{\ast}(1^{-})$ | 5.415
${A}_{1,2}(q^{2})$, ${T}_{2,3}(q^{2})$ | $B_{s1}(1^{+})$ | 5.829
  4. (4)

In Fig. 1, we show the theoretical predictions of the form factors
$V(q^{2})$,$A_{0,1,2}(q^{2})$ and $T_{1,2,3}(q^{2})$ for $B_{s}\to\phi$
transition based on the PQCD approach ( red curves) and the “PQCD+Lattice”
approach (blue curves) with the extrapolation from $q^{2}=0$ to
$q^{2}_{max}=(m_{B_{s}}-m_{\phi})^{2}$ by applying the $z$-series
parameterizations. The shaded bands represent the total theoretical error
obtained by adding in quadrature of the separate errors from the uncertainty
of the parameter $\omega_{B_{s}}$, $a^{||,\perp}_{2\phi}$, $f_{\phi}$ and
$f_{\phi}^{\perp}$. The black error bars in the low-$q^{2}$ region correspond
to the PQCD predictions of the corresponding form factors, while the error
bars in the high-$q^{2}$ region are the currently known lattice QCD results as
collected in Table 4

Figure 1: Theoretical predictions of the relevant form factors for
$B_{s}\to\phi$ transition in the PQCD approach ( red curves) and the
“PQCD+Lattice” approach (blue curves). The red (blue) shaded band represents
the theory uncertainties. The black error bars in the low-$q^{2}$ region
correspond to the PQCD predictions. The error bars in high-$q^{2}$ region are
the lattice QCD results.

###  IV.2 Numerical results

We now proceed to explore the phenomenological aspects of the cascade decays
$B_{s}\to\phi(\to K^{-}K^{+})\ell^{+}\ell^{-}$, which allow us to define and
compute a number of physical observables and compare them with those measured
by experiments. We first compare our results for the branching ratios and
angular observables with the experimental data reported by the LHCb
CollaborationAaij:2015esa ; Aaij:2021pkz . As studied systematically in last
section, the physical observables accessible in the semileptonic decays
$B_{s}\to\phi\ell^{+}\ell^{-}$ are the CP averaged differential branching
fraction $d{\cal B}/dq^{2}$ Aaij:2015esa ; Aaij:2021pkz , the CP-averaged
$\phi$ meson longitudinal polarization fraction $F_{L}$, the forward-backward
asymmetry ${\cal A}_{FB}$, the angular coefficients $S_{i}$ and $A_{i}$, and
the optimized observables $P_{i}$ and $P_{j}^{\prime}$ Descotes-Genon:2015uva
. The CP asymmetry angular coefficients $A_{5,6,8,9}$ in the SM are induced by
the weak phase from the CKM matrix. For the $b\to s$ transition, the CP
asymmetries proportional to ${\rm
Im}(V_{ub}V^{\ast}_{us}/V_{tb}V^{\ast}_{ts})$, which is of order $10^{-2}$
Bobeth:2008ij  as measured by the LHCb Collaboration (see Table 3 in Ref.
Aaij:2015esa ), but the statistical uncertainties are still large. For these
reasons, we will focus on the CP averaged quantities when taking the binned
observables into consideration.

Figure 2: Theoretical predictions for the $q^{2}$-dependence of differential
branching fraction $d{\cal B}/dq^{2}$ for the semileptonic decays
$B_{s}\to\phi\mu^{+}\mu^{-}$ in the PQCD (red band) and “PQCD+Lattice” (blue
band) approach, respectively. The crosses show the LHCb measurements in
different bins as given in Ref. Aaij:2015esa  ( orange ones ) and in Ref.
Aaij:2021pkz  (black ones). The vertical grey blocks are the experimental veto
regions.

We begin with the branching ratios of the decays
$B_{s}\to\phi\ell^{+}\ell^{-}$. From the differential decay rates as defined
in Eq. (40), it is straightforward to make the integration over the range of
$4m_{\ell}^{2}\\!\leq\\!q^{2}\\!\leq\\!(m_{B_{s}}-m_{\phi})^{2}$. In order to
be consistent with the choices made by LHCb Collaboration in their data
analysis, we here also cut off the regions of dilepton-mass squared around the
charmonium resonances $J/\psi(1S)$ and $\psi(2S)$: i.e.,
$8.0\\!<\\!q^{2}\\!<\\!11.0\,GeV^{2}$ and
$12.5\\!<\\!q^{2}\\!<\\!15.0\,GeV^{2}$ for $\ell=(e,\mu,\tau)$ cases. We
display the PQCD and “PQCD+Lattice” predictions for the differential branching
ratios $d{\cal B}/dq^{2}$ in Fig. 2 for the cases of $l=(\mu,\tau)$, including
currently available LHCb results in six or eight bins of $q^{2}$ Aaij:2015esa
; Aaij:2021pkz  indicated by the crosses for $B_{s}\to\phi\mu^{+}\mu^{-}$
decay. From Fig. 2 one can see that both PQCD and “PQCD+Lattice” predictions
for the differential branching ratios do agree well with the LHCb results
within the still large errors. Since the theoretical prediction for the
differential branching ratio of the electron mode is almost identical with the
one of the muon, we do not draw the figure of $d{\cal B}(B_{s}\to\phi
e^{+}e^{-})/dq^{2}$ in Fig. 2.

Table 6: Theoretical predictions for the total branching fractions ${\cal B}(B_{s}\to\phi\ell^{+}\ell^{-})$ ( in units of $10^{-7}$ ) in the PQCD (the first row) and “PQCD+Lattice” (the second row) approaches. As a comparison, we also list the LHCb measured value for muon channel corresponding to an integrated luminosity of 3$\rm fb^{-1}$Aaij:2015esa  and 9$\rm fb^{-1}$Aaij:2021pkz  and the QCDSR predictions for all three channels Peng:2020ivn .  BFs | PQCD / “PQCD+Lattice” | QCDSR Peng:2020ivn  | LHCb
---|---|---|---
${\cal B}(B_{s}\to\phi e^{+}e^{-})$ | ${8.55}^{+4.02}_{-2.69}(\rm FFs)\pm{0.15}(\mu)\pm{0.42}(V_{tb})\pm{0.65}(V_{ts})$ | $7.12\pm 1.40$ |
| ${8.24}^{+2.03}_{-1.53}(\rm FFs)\pm{0.14}(\mu)\pm{0.41}(V_{tb})\pm{0.63}(V_{ts})$ |  |
${\cal B}(B_{s}\to\phi\mu^{+}\mu^{-})$ | ${7.07}^{+3.37}_{-2.25}(\rm FFs)\pm{0.12}(\mu)\pm{0.38}(V_{tb})\pm{0.53}(V_{ts})$ | $7.06\pm 1.59$ |  $7.97^{+0.81}_{-0.80}$Aaij:2015esa
| ${6.76}^{+1.39}_{-1.09}(\rm FFs)\pm{0.11}(\mu)\pm{0.33}(V_{tb})\pm{0.52}(V_{ts})$ |  |  $8.14\pm 0.47$Aaij:2021pkz
${\cal B}(B_{s}\to\phi\tau^{+}\tau^{-})$ | ${0.81}^{+0.42}_{-0.27}(\rm FFs)\pm{0.02}(\mu)\pm{0.04}(V_{tb})\pm{0.06}(V_{ts})$ | $0.35\pm 0.17$ |
| ${0.68}^{+0.06}_{-0.06}(\rm FFs)\pm{0.02}(\mu)\pm{0.03}(V_{tb})\pm{0.05}(V_{ts})$ |  |

In Table 6 we present the theoretical predictions of the total branching
fractions for $B_{s}\to\phi\ell^{+}\ell^{-}$ with $\ell=(e,\mu,\tau)$ obtained
by the integration over the six $q^{2}$ bins using the PQCD ( the first row)
and “PQCD+Lattice” approach ( the second row), respectively. The major
theoretical errors from different sources, such as the form factors (FFs) as
listed in Table 2, the scale $\mu$, the CKM matrix element $V_{tb}$ and
$V_{ts}$, are also listed. As in Ref. Aaij:2015esa , a correction factor
$f_{veto}=1.52$ is applied to account for the contribution in the veto $q^{2}$
bins for $\ell=(e,\mu)$ cases. As a comparison, we also show the LHCb measured
value ${\cal B}(B_{s}\to\phi\mu^{+}\mu^{-})=(7.97^{+0.81}_{-0.80})\times
10^{-7}$ Aaij:2015esa and $(8.14^{+0.47}_{-0.47})\times 10^{-7}$Aaij:2021pkz
and the QCDSR predictions ${\cal B}(B_{s}\to\phi\ell^{+}\ell^{-})$ for all
three decay modes Peng:2020ivn . For $B_{s}\to\phi\mu^{+}\mu^{-}$ decay, for
instance, the theoretical predictions and the LHCb measurement Aaij:2015esa ;
Aaij:2021pkz  (in unit of $10^{-7}$) are the following:

$\displaystyle{\cal
B}(B_{s}\to\phi\mu^{+}\mu^{-})=\left\\{\begin{array}[]{ll}7.07^{+3.43}_{-2.34},&{\rm
in\ \ PQCD},\\\ 6.76^{+1.52}_{-1.25},&{\rm in\ \ PQCD+Lattice},\\\
7.06^{+1.59}_{-1.59},&{\rm in\ \ QCDSR}\ [87],\\\ 7.97^{+0.81}_{-0.80},&{\rm
LHCb}\ \ [31].\\\ 8.14^{+0.47}_{-0.47},&{\rm LHCb}\ \ [33].\\\
\end{array}\right.$ (64)

From the numerical results in above equation and Table 6, one can see that

  1. (1)

The PQCD and “PQCD+Lattice” predictions for the branching ratio ${\cal
B}(B_{s}\to\phi\ell^{+}\ell^{-})$ with $\ell=(e,\mu,\tau)$ do agree well with
each other within the errors, while the “PQCD+ Lattice” predictions of ${\cal
B}(B_{s}\to\phi\ell^{+}\ell^{-})$ have smaller errors than those of the PQCD
predictions.

  2. (2)

Both PQCD and “PQCD+Lattice” predictions of ${\cal
B}(B_{s}\to\phi\mu^{+}\mu^{-})$ do agree well with currently available LHCb
measured values Aaij:2015esa ; Aaij:2021pkz  within errors. For the electron
and tau mode, however, we have to wait for the future experimental
measurements.

  3. (3)

For all three decay modes, our theoretical predictions of the branching ratios
do agree well with the theoretical predictions obtained from the QCD sum rule
Peng:2020ivn .

Since the large theoretical uncertainties of the branching ratios could be
largely canceled in the ratio of the branching ratios of
$B_{s}\to\phi\ell^{+}\ell^{-}$ decays, one can define and check the physical
observables $R^{e\mu}_{\phi}$ and $R^{\mu\tau}_{\phi}$ Hiller:2003js . In the
region $q^{2}<4m^{2}_{\mu}$, where only the $e^{+}e^{-}$ modes are allowed,
there is a large enhancement due to the $1/q^{2}$ scaling of the photon
penguin contribution Aubert:2008ps . In order to remove the phase space
effects in the ratio $R^{e\mu}_{\phi}$ and keep consistent with other analysis
Hiller:2003js , we here also use the lower cut of $4m^{2}_{\mu}$ for both the
electron and muon modes in the definition of the ratio $R^{e\mu}_{\phi}$ as in
Ref. Hiller:2003js :

$\displaystyle
R^{e\mu}_{\phi}=\frac{\int^{q^{2}_{max}}_{4m^{2}_{\mu}}dq^{2}\frac{d{\cal
B}(B_{s}\to\phi\mu^{+}\mu^{-})}{dq^{2}}}{\int^{q^{2}_{max}}_{4m^{2}_{\mu}}dq^{2}\frac{d{\cal
B}(B_{s}\to\phi e^{+}e^{-})}{dq^{2}}}=\left\\{\begin{array}[]{ll}0.992\pm
0.002,&{\rm in\ \ PQCD},\\\ 0.991\pm 0.002,&{\rm in\ \ PQCD+Lattice}.\\\
\end{array}\right.$ (67)

For the case of the ratio $R^{\mu\tau}_{\phi}$ we have

$\displaystyle
R^{\mu\tau}_{\phi}=\frac{\int^{q^{2}_{max}}_{4m^{2}_{\tau}}dq^{2}\frac{d{\cal
B}(B_{s}\to\phi\tau^{+}\tau^{-})}{dq^{2}}}{\int^{q^{2}_{max}}_{4m^{2}_{\mu}}dq^{2}\frac{d{\cal
B}(B_{s}\to\phi\mu^{+}\mu^{-})}{dq^{2}}}=\left\\{\begin{array}[]{ll}0.115\pm
0.004,&{\rm in\ \ PQCD},\\\ 0.100\pm 0.009,&{\rm in\ \ PQCD+Lattice},\\\
\end{array}\right.$ (70)

where the total error is the combination of the individual errors in
quadrature. We suggest the LHCb and Belle-II to measure these two ratios.

Figure 3:  Theoretical predictions for the $q^{2}$-dependence of the
observables $F_{L}(q^{2})$, $S_{3,4,7}(q^{2})$ and $A_{5,6,8,9}(q^{2})$ for
the decay $B_{s}\to\phi\mu^{+}\mu^{-}$ in the PQCD ( red lines ) and
“PQCD+Lattice” ( blue lines ) approach. The crosses in each figure represent
the LHCb measurements in six $q^{2}$ bins Aaij:2015esa . The vertical grey
blocks are the two experimental veto regions.

Figure 4: Theoretical predictions for the $q^{2}$-dependence of the
observables $A_{FB}(q^{2})$, $P_{1,2,3}(q^{2})$ and
$P^{\prime}_{4,5,6,8}(q^{2})$ for $B_{s}\to\phi\mu^{+}\mu^{-}$ decay in the
PQCD (red lines) and “PQCD+Lattice” (blue lines )approach. The vertical grey
blocks are the two experimental veto regions.

For $B_{s}\to\phi\mu^{+}\mu^{-}$ decay, we show our theoretical predictions
for the $q^{2}$-dependence of the longitudinal polarization $F_{L}(q^{2})$,
the CP averaged angular coefficients $S_{3,4,7}(q^{2})$ and the CP asymmetry
angular coefficients $A_{5,6,8,9}(q^{2})$ in Fig. 3. As a comparison, the
currently available LHCb measurements for these observables of
$B_{s}\to\phi\mu^{+}\mu^{-}$ decay in the six $q^{2}$ bins Aaij:2015esa  are
also shown by those crosses explicitly. One can see from the Fig. 3 that:

  1. (1)

For the longitudinal polarization $F_{L}(q^{2})$, although both PQCD and “PQCD
+ Lattice” predictions all agree well with the LHCb measurements in the six
bins, our theoretical predictions in the region of the fourth and fifth bin
are little larger than the measured ones.

  2. (2)

For the CP averaged angular coefficients $S_{3,4,7}(q^{2})$, the PQCD and
“PQCD + Lattice” predictions agree very well with each other, and are
consistent with the LHCb results within the still large experimental errors.
For the last two high $q^{2}$ bins, the LHCb results of $S_{3}$ ($S_{7}$) is a
little larger (smaller) than our theoretical predictions.

  3. (3)

For the CP asymmetry angular coefficients $A_{5,6,8,9}(q^{2})$ , the PQCD and
“PQCD + Lattice” predictions are very small: in the range of $10^{-4}$ to
$10^{-2}$. For the LHCb measurements in the six bins, they are clearly
consistent with our theoretical predictions due to still large experimental
errors.

In Fig. 4, we show our theoretical predictions for the $q^{2}$-dependence of
the forward-backward asymmetry $A_{FB}(q^{2})$, the optimized observables
$P_{1,2,3}(q^{2})$ and $P^{\prime}_{4,5,6,8}(q^{2})$ for
$B_{s}\to\phi\mu^{+}\mu^{-}$ decay. Unfortunately, there exist no any
experimental measurements for these observables. We have to wait for future
LHCb and Belle-II measurements. Analogous to Fig. 3, the vertical grey blocks
in Fig. 4 also denote the two experimental veto regions of $q^{2}$:
$8.0\\!<\\!q^{2}\\!<\\!11.0\,GeV^{2}$ and
$12.5\\!<\\!q^{2}\\!<\\!15.0\,GeV^{2}$.

Table 7: Theoretical predictions for the observables $F^{\phi}_{L}$, $A_{FB}$, $S_{3,4,7}$, $A_{5,6,8,9}$, $P_{1,2,3}$ and $P^{\prime}_{4,5,6,8}$ integrated over the whole kinematic region for $B_{s}\to\phi\ell^{+}\ell^{-}$ decays in the PQCD (the first row) and “PQCD+Lattice” (the second row) approaches, respectively.  ${\rm Obs.}$ | $\ell=e$ | $\ell=\mu$ | $\ell=\tau$ | ${\rm Obs.}$ | $\ell=e$ | $\ell=\mu$ | $\ell=\tau$
---|---|---|---|---|---|---|---
$F^{\phi}_{L}$ | ${0.383}^{+0.003}_{-0.004}$ | ${0.454}^{+0.006}_{-0.007}$ | ${0.396}^{+0.002}_{-0.003}$ | $-A_{FB}$ | ${0.192}^{+0.003}_{-0.004}$ | ${0.233}^{+0.004}_{-0.004}$ | ${0.173}^{+0.002}_{-0.002}$
| ${0.446}^{+0.012}_{-0.013}$ | ${0.533}^{+0.001}_{-0.001}$ | ${0.442}^{+0.002}_{-0.002}$ |  | ${0.152}^{+0.010}_{-0.008}$ | ${0.186}^{+0.005}_{-0.004}$ | ${0.151}^{+0.002}_{-0.003}$
$-S_{3}$ | ${0.120}^{+0.005}_{-0.004}$ | ${0.144}^{+0.005}_{-0.005}$ | ${0.080}^{+0.001}_{-0.001}$ | $-P_{1}$ | ${0.399}^{+0.014}_{-0.011}$ | ${0.555}^{+0.010}_{-0.012}$ | ${0.795}^{+0.007}_{-0.004}$
| ${0.102}^{+0.009}_{-0.008}$ | ${0.124}^{+0.007}_{-0.006}$ | ${0.075}^{+0.001}_{-0.001}$ |  | ${0.381}^{+0.041}_{-0.039}$ | ${0.564}^{+0.026}_{-0.024}$ | ${0.817}^{+0.007}_{-0.005}$
$S_{4}$ | ${0.210}^{+0.004}_{-0.004}$ | ${0.258}^{+0.003}_{-0.004}$ | ${0.100}^{+0.001}_{-0.001}$ | $-P_{2}$ | ${0.213}^{+0.002}_{-0.003}$ | ${0.299}^{+0.001}_{-0.001}$ | ${0.281}^{+0.003}_{-0.004}$
| ${0.201}^{+0.012}_{-0.013}$ | ${0.249}^{+0.007}_{-0.006}$ | ${0.098}^{+0.001}_{-0.001}$ |  | ${0.189}^{+0.015}_{-0.013}$ | ${0.282}^{+0.005}_{-0.004}$ | ${0.268}^{+0.004}_{-0.004}$
$10^{3}S_{7}$ | ${0.350}^{+0.010}_{-0.011}$ | ${0.371}^{+0.012}_{-0.013}$ | ${0.022}^{+0.001}_{-0.000}$ | $10^{2}P_{3}$ | ${0.103}^{+0.003}_{-0.006}$ | ${0.142}^{+0.002}_{-0.007}$ | ${0.122}^{+0.003}_{-0.006}$
| ${0.392}^{+0.026}_{-0.024}$ | ${0.423}^{+0.040}_{-0.035}$ | ${0.030}^{+0.001}_{-0.000}$ |  | ${0.118}^{+0.012}_{-0.012}$ | ${0.174}^{+0.007}_{-0.008}$ | ${0.158}^{+0.002}_{-0.003}$
$-10^{3}A_{5}$ | ${0.489}^{+0.010}_{-0.009}$ | ${0.585}^{+0.012}_{-0.013}$ | ${0.040}^{+0.000}_{-0.000}$ | $P^{\prime}_{4}$ | ${1.049}^{+0.013}_{-0.015}$ | ${1.111}^{+0.011}_{-0.012}$ | ${1.338}^{+0.002}_{-0.002}$
| ${0.518}^{+0.020}_{-0.020}$ | ${0.623}^{+0.041}_{-0.046}$ | ${0.038}^{+0.001}_{-0.001}$ |  | ${1.033}^{+0.033}_{-0.037}$ | ${1.098}^{+0.024}_{-0.027}$ | ${1.345}^{+0.002}_{-0.003}$
$-10^{3}A_{6}$ | ${0.482}^{+0.001}_{-0.001}$ | ${0.577}^{+0.002}_{-0.004}$ | ${0.068}^{+0.000}_{-0.001}$ | $-P^{\prime}_{5}$ | ${0.482}^{+0.002}_{-0.004}$ | ${0.524}^{+0.003}_{-0.005}$ | ${0.427}^{+0.005}_{-0.006}$
| ${0.429}^{+0.020}_{-0.021}$ | ${0.518}^{+0.038}_{-0.043}$ | ${0.058}^{+0.001}_{-0.001}$ |  | ${0.470}^{+0.008}_{-0.010}$ | ${0.515}^{+0.002}_{-0.003}$ | ${0.408}^{+0.006}_{-0.006}$
$10^{4}A_{8}$ | ${0.541}^{+0.030}_{-0.027}$ | ${0.241}^{+0.017}_{-0.020}$ | ${0.011}^{+0.000}_{-0.000}$ | $-10^{3}P^{\prime}_{6}$ | ${0.875}^{+0.031}_{-0.030}$ | ${0.780}^{+0.028}_{-0.029}$ | ${0.069}^{+0.001}_{-0.001}$
| ${0.550}^{+0.088}_{-0.095}$ | ${0.237}^{+0.060}_{-0.070}$ | ${0.014}^{+0.000}_{-0.000}$ |  | ${1.007}^{+0.089}_{-0.101}$ | ${0.910}^{+0.078}_{-0.089}$ | ${0.097}^{+0.001}_{-0.001}$
$10^{4}A_{9}$ | ${0.040}^{+0.002}_{-0.004}$ | ${0.054}^{+0.002}_{-0.004}$ | ${0.013}^{+0.001}_{-0.001}$ | $-10^{2}P^{\prime}_{8}$ | ${0.643}^{+0.009}_{-0.003}$ | ${0.640}^{+0.004}_{-0.008}$ | ${0.287}^{+0.002}_{-0.002}$
| ${0.041}^{+0.006}_{-0.008}$ | ${0.056}^{+0.004}_{-0.007}$ | ${0.015}^{+0.001}_{-0.001}$ |  | ${0.749}^{+0.004}_{-0.004}$ | ${0.748}^{+0.006}_{-0.006}$ | ${0.374}^{+0.001}_{-0.001}$
Table 8: Theoretical predictions for the $q^{2}$-binned observables ${\cal B}(B_{s}\to\phi\mu^{+}\mu^{-})$ (in unit of $10^{-7}$) in the PQCD (the first low) and “PQCD+Lattice” (the second row) approach, respectively. For a comparison, we also list the new LHCb measurements Aaij:2021pkz . $q^{2}$ bin $(\rm GeV^{2})$ | ${\cal B}(\ell=\mu)$ | LHCb Aaij:2021pkz  |  $q^{2}$ bin $(\rm GeV^{2})$ | ${\cal B}(\ell=\mu)$ | LHCb Aaij:2021pkz
---|---|---|---|---|---
$[0.10,0.98]$ | ${0.25}^{+0.10}_{-0.07}$ | ${0.68}\pm{0.06}$ | $[1.10,2.50]$ | ${0.19}^{+0.08}_{-0.05}$ | ${0.44}\pm{0.05}$
| ${0.25}^{+0.10}_{-0.07}$ |  |  | ${0.21}^{+0.07}_{-0.05}$ |
$[2.50,4.00]$ | ${0.23}^{+0.09}_{-0.07}$ | ${0.35}\pm{0.04}$ | $[4.00,6.00]$ | ${0.41}^{+0.18}_{-0.12}$ | ${0.62}\pm{0.06}$
| ${0.27}^{+0.09}_{-0.06}$ |  |  | ${0.50}^{+0.14}_{-0.11}$ |
$[6.00,8.00]$ | ${0.61}^{+0.28}_{-0.19}$ | ${0.63}\pm{0.06}$ | $[11.0,12.5]$ | ${0.78}^{+0.38}_{-0.25}$ | ${0.72}\pm{0.06}$
| ${0.69}^{+0.18}_{-0.13}$ |  |  | ${0.76}^{+0.14}_{-0.11}$ |
$[15.00,17.00]$ | ${1.26}^{+0.64}_{-0.42}$ | ${1.05}\pm{0.08}$ | $[17.0,19.0]$ | ${0.86}^{+0.45}_{-0.30}$ | ${0.84}\pm{0.07}$
| ${1.04}^{+0.11}_{-0.10}$ |  |  | ${0.67}^{+0.05}_{-0.06}$ |
$[1.10,6.00]$ | ${0.83}^{+0.35}_{-0.24}$ | ${1.41}\pm{0.10}$ | $[15.0,19.0]$ | ${2.12}^{+1.09}_{-0.71}$ | ${1.85}\pm{0.13}$
| ${0.98}^{+0.30}_{-0.22}$ |  |  | ${1.70}^{+0.17}_{-0.16}$ |
Table 9: Theoretical predictions for the $q^{2}$-binned observables ${\cal B}(B_{s}\to\phi\ell^{+}\ell^{-})$ (in unit of $10^{-7}$) and $F^{\phi}_{L}$ with $\ell=(\mu,\tau)$ in the PQCD (the first low) and “PQCD+Lattice” (the second row) approach, respectively. For a comparison, we also list the LHCb measurements as given in Refs. Aaij:2015esa ; Aaij:2021pkz .  $q^{2}$ bin $(\rm GeV^{2})$ | ${\cal B}(\ell=\mu)$ | LHCb | ${\cal B}(\ell=\tau)$ | $F_{L}^{\phi}(\ell=\mu)$ | LHCb Aaij:2015esa  | $F_{L}^{\phi}(\ell=\tau)$
---|---|---|---|---|---|---
$[0.1,2.0]$ | ${0.39}^{+0.16}_{-0.11}$ |  ${1.11}\pm{0.16}$ Aaij:2015esa  | $-$ | ${0.441}^{+0.007}_{-0.008}$ | ${0.20}\pm{0.09}$ | $-$
| ${0.40}^{+0.16}_{-0.11}$ |  | $-$ | ${0.472}^{+0.011}_{-0.012}$ |  | $-$
$[2.0,5.0]$ | ${0.48}^{+0.20}_{-0.14}$ |  ${0.77}\pm{0.14}$ Aaij:2015esa  | $-$ | ${0.738}^{+0.008}_{-0.009}$ | ${0.68}\pm{0.15}$ | $-$
| ${0.58}^{+0.18}_{-0.13}$ |  | $-$ | ${0.796}^{+0.007}_{-0.007}$ |  | $-$
$[5.0,8.0]$ | ${0.84}^{+0.37}_{-0.25}$ |  ${0.96}\pm{0.15}$ Aaij:2015esa  | $-$ | ${0.584}^{+0.008}_{-0.009}$ | ${0.54}\pm{0.10}$ | $-$
| ${0.96}^{+0.25}_{-0.19}$ |  | $-$ | ${0.682}^{+0.010}_{-0.008}$ |  | $-$
$[11.0,12.5]$ | ${0.78}^{+0.38}_{-0.25}$ |  ${0.71}\pm{0.12}$ Aaij:2015esa  | $-$ | ${0.433}^{+0.003}_{-0.005}$ | ${0.29}\pm{0.11}$ | $-$
| ${0.76}^{+0.14}_{-0.11}$ |  $0.72\pm 0.06$Aaij:2021pkz  | $-$ | ${0.524}^{+0.008}_{-0.008}$ |  | $-$
$[15.0,17.0]$ | ${1.12}^{+0.61}_{-0.40}$ |  ${0.90}\pm{0.13}$ Aaij:2015esa  | ${0.45}^{+0.22}_{-0.14}$ | ${0.368}^{+0.001}_{-0.001}$ | ${0.23}\pm{0.09}$ | ${0.421}^{+0.003}_{-0.004}$
| ${0.99}^{+0.11}_{-0.10}$ |  $1.05\pm 0.08$Aaij:2021pkz  | ${0.39}^{+0.04}_{-0.04}$ | ${0.412}^{+0.004}_{-0.004}$ |  | ${0.479}^{+0.003}_{-0.003}$
$[17.0,19.0]$ | ${0.80}^{+0.41}_{-0.27}$ |  ${0.75}\pm{0.13}$ Aaij:2015esa  | ${0.36}^{+0.18}_{-0.12}$ | ${0.346}^{+0.001}_{-0.001}$ | ${0.40}\pm{0.14}$ | ${0.365}^{+0.001}_{-0.002}$
| ${0.67}^{+0.05}_{-0.05}$ |  $0.82\pm 0.07$Aaij:2021pkz  | ${0.29}^{+0.02}_{-0.02}$ | ${0.363}^{+0.001}_{-0.001}$ |  | ${0.391}^{+0.001}_{-0.001}$
$[1.0,6.0]$ | ${0.93}^{+0.29}_{-0.27}$ |  ${1.29}\pm{0.19}$ Aaij:2015esa  | $-$ | ${0.708}^{+0.007}_{-0.009}$ | ${0.63}\pm{0.09}$ | $-$
| ${1.10}^{+0.34}_{-0.25}$ |  | $-$ | ${0.777}^{+0.008}_{-0.006}$ |  | $-$
$[15.0,19.0]$ | ${1.99}^{+1.02}_{-0.71}$ |  ${1.62}\pm{0.20}$ Aaij:2015esa  | ${0.82}^{+0.40}_{-0.26}$ | ${0.359}^{+0.002}_{-0.001}$ | ${0.29}\pm{0.07}$ | ${0.396}^{+0.003}_{-0.003}$
| ${1.60}^{+0.16}_{-0.16}$ |  $1.85\pm 0.13$Aaij:2021pkz  | ${0.68}^{+0.06}_{-0.06}$ | ${0.394}^{+0.003}_{-0.003}$ |  | ${0.442}^{+0.002}_{-0.001}$

In Table 7, we list the theoretical predictions for the values of the
observables $F^{\phi}_{L}$, $A_{FB}$, $S_{3,4,7}$, $A_{5,6,8,9}$, $P_{1,2,3}$
and $P^{\prime}_{4,5,6,8}$, obtained after the integrations over the whole
kinematic region of $q^{2}$ for the semileptonic decays
$B_{s}\to\phi\ell^{+}\ell^{-}$ with $\ell=(e,\mu,\tau)$ in the PQCD (the first
row) and “PQCD+Lattice” (the second row) approaches, respectively. Of course,
the regions corresponding to resonance $J/\psi(1S)$ and $\psi(2S)$, say
$8.0\\!<\\!q^{2}\\!<\\!11.0\,GeV^{2}$ and
$12.5\\!<\\!q^{2}\\!<\\!15.0\,GeV^{2}$ numerically, are also cut off here. The
total errors are the combinations of the individual errors from the form
factors, the renormalization scales and the relevant CKM matrix elements. The
above theoretical predictions should be tested in the near future LHCb and
Belle-II experiments. For the considered $B_{s}$ meson decays, one should
consider the effects from the $B_{s}$-$\bar{B}_{s}$ mixing Hiller2015 . The
theoretical framework for examining the time-dependent decays with the
inclusion of such mixing effects can been found in Ref. Descotes-Genon:2015hea
. The authors of Ref. Descotes-Genon:2015hea  proved that the mixing effects
on the values of decay rates and CP averaged observables are generally within
a few percent and could be neglected.

###  IV.3 The $q^{2}$-binned predictions

Table 10: Theoretical predictions for the $q^{2}$-binned observables $S_{3,4,7}$ of the decays $B_{s}\to\phi\mu^{+}\mu^{-}$ in the PQCD (the first low) and “PQCD+Lattice” (the second row) approaches. For a comparison, we also list the LHCb measured values Aaij:2015esa . | $S_{3}$ | $S_{4}$ | $S_{7}$
---|---|---|---
$q^{2}$ bin $(\rm GeV^{2})$ | Theor. | LHCb | Theor. | LHCb | Theor. ($10^{-3}$) | LHCb
$[0.1,2.0]$ | ${0.003}^{+0.000}_{-0.000}$ | ${-0.05}\pm{0.13}$ | ${-0.054}^{+0.001}_{-0.001}$ | ${-0.27}\pm{0.23}$ | $1.571^{+0.001}_{-0.002}$ | ${-0.04}\pm{0.12}$
| ${0.002}^{+0.000}_{-0.001}$ |  | ${-0.053}^{+0.000}_{-0.000}$ |  | $1.551^{+0.003}_{-0.006}$ |
$[2.0,5.0]$ | ${-0.024}^{+0.002}_{-0.002}$ | ${-0.06}\pm{0.21}$ | ${0.191}^{+0.004}_{-0.004}$ | ${0.47}\pm{0.37}$ | $1.065^{+0.003}_{-0.003}$ | ${0.03}\pm{0.21}$
| ${-0.021}^{+0.001}_{-0.001}$ |  | ${0.177}^{+0.003}_{-0.003}$ |  | ${0.979}^{+0.017}_{-0.018}$ |
$[5.0,8.0]$ | ${-0.057}^{+0.003}_{-0.003}$ | ${-0.10}\pm{0.25}$ | ${0.270}^{+0.002}_{-0.002}$ | ${0.10}\pm{0.17}$ | ${0.453}^{+0.003}_{-0.003}$ | ${-0.04}\pm{0.18}$
| ${-0.050}^{+0.004}_{-0.003}$ |  | ${0.259}^{+0.003}_{-0.003}$ |  | ${0.453}^{+0.004}_{-0.005}$ |
$[11.0,12.5]$ | ${-0.124}^{+0.004}_{-0.003}$ | ${-0.19}\pm{0.21}$ | ${0.296}^{+0.001}_{-0.001}$ | ${0.47}\pm{0.25}$ | ${0.153}^{+0.002}_{-0.002}$ | ${0.00}\pm{0.16}$
| ${-0.115}^{+0.002}_{-0.001}$ |  | ${0.303}^{+0.001}_{-0.001}$ |  | ${0.185}^{+0.001}_{-0.001}$ |
$[15.0,17.0]$ | ${-0.219}^{+0.003}_{-0.003}$ | ${-0.06}\pm{0.18}$ | ${0.314}^{+0.001}_{-0.001}$ | ${0.03}\pm{0.15}$ | ${0.052}^{+0.001}_{-0.001}$ | ${-0.12}\pm{0.15}$
| ${-0.213}^{+0.001}_{-0.001}$ |  | ${0.323}^{+0.001}_{-0.001}$ |  | ${0.071}^{+0.001}_{-0.001}$ |
$[17.0,19.0]$ | ${-0.283}^{+0.001}_{-0.001}$ | ${-0.07}\pm{0.25}$ | ${0.325}^{+0.001}_{-0.001}$ | ${0.39}\pm{0.30}$ | ${0.019}^{+0.001}_{-0.001}$ | ${-0.20}\pm{0.26}$
| ${-0.281}^{+0.001}_{-0.001}$ |  | ${0.329}^{+0.001}_{-0.001}$ |  | ${0.027}^{+0.001}_{-0.001}$ |
$[1.0,6.0]$ | ${-0.026}^{+0.002}_{-0.002}$ | ${-0.02}\pm{0.13}$ | ${0.180}^{+0.004}_{-0.004}$ | ${0.19}\pm{0.14}$ | ${1.063}^{+0.001}_{-0.001}$ | ${0.03}\pm{0.14}$
| ${-0.023}^{+0.001}_{-0.001}$ |  | ${0.169}^{+0.001}_{-0.001}$ |  | ${0.985}^{+0.020}_{-0.022}$ |
$[15.0,19.0]$ | ${-0.245}^{+0.002}_{-0.002}$ | ${-0.09}\pm{0.12}$ | ${0.318}^{+0.001}_{-0.001}$ | ${0.14}\pm{0.11}$ | ${0.038}^{+0.002}_{-0.002}$ | ${-0.13}\pm{0.11}$
| ${-0.239}^{+0.001}_{-0.001}$ |  | ${0.325}^{+0.001}_{-0.001}$ |  | ${0.054}^{+0.001}_{-0.001}$ |

For $B_{s}\to\phi\mu^{+}\mu^{-}$ decay mode, the LHCb Collaboration has
reported their experimental measurements for many physical observables in
several $q^{2}$ bins Aaij:2015esa ; Aaij:2021pkz . In order to compare our
theoretical predictions with the LHCb results bin by bin, we make the same
choices of the $q^{2}$ bins as LHCb did, calculate and show our theoretical
predictions for the branching ratio ${\cal B}(B_{s}\to\phi\ell^{+}\ell^{-})$
and the asymmetry $F_{L}^{\phi}$ with $\ell=(\mu,\tau)$ in Table 8 and 9, and
the observables $S_{3,4,7}$ with $\ell=\mu$ in Table 10. For observables
$S_{7}$ and $A_{5,6,8,9}$, in fact, our theoretical predictions for their
values are very small, say in the range of $10^{-3}-10^{-4}$ in magnitude, but
still agree with the LHCb measurements in different bins Aaij:2015esa  due to
still large experimental errors. For the observables $P_{3}$ and
$P^{\prime}_{6,8}$, they are also very small in size: in the range of
$10^{-3}-10^{-4}$ and there exist no corresponding data at present. For
observables $P_{1,2}$ and $P^{\prime}_{4,5}$, on the other hand, although
there exist no experimental measurements for them at present, they are
relatively large in size and may be measured in the near future LHCb and
Belle-II experiments, so we calculate and list the theoretical predictions of
these observables bin by bin for the cases of $\ell=(\mu,\tau)$ in Table 11
and 12. Very recently, LHCb reported some new measurements for the angular
observables of $B_{s}\to\phi\mu^{+}\mu^{-}$ decay LHCb-2107  in the $q^{2}$
bins different from those in their previous work Aaij:2015esa , which will be
studied in our next work.

The definitions of the $q^{2}$-binned observables are the following:

$\displaystyle{\cal
B}(q^{2}_{1},q^{2}_{2})=\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}\frac{d{\cal
B}(B_{s}\to\phi\ell^{+}\ell^{-})}{dq^{2}},$ (71) $\displaystyle
F_{L}^{\phi}(q^{2}_{1},q^{2}_{2})=\frac{\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}[3(I^{c}_{1}+\bar{I}^{c}_{1})-(I^{c}_{2}+\bar{I}^{c}_{2})]}{4\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}[d(\Gamma+\bar{\Gamma})/dq^{2}]},\quad$
(72) $\displaystyle{A}_{\rm
FB}(q^{2}_{1},q^{2}_{2})=\frac{3\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(I^{s}_{6}+\bar{I}^{s}_{6})}{4\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}[d(\Gamma+\bar{\Gamma})/dq^{2}]},$
(73) $\displaystyle
S_{3,4,7}(q^{2}_{1},q^{2}_{2})=\frac{\int_{q_{1}^{2}}^{q_{2}^{2}}dq^{2}(I_{3,4,7}+\bar{I}_{3,4,7})}{\int_{q_{1}^{2}}^{q_{2}^{2}}dq^{2}[d(\Gamma+\bar{\Gamma})/dq^{2}]},\quad$
(74) $\displaystyle P_{1}(q^{2}_{1},q^{2}_{2})$ $\displaystyle=$
$\displaystyle\frac{\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(S_{3})}{2\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(S^{s}_{2})},\quad
P_{2}(q^{2}_{1},q^{2}_{2})=\frac{\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(\beta_{\ell}S^{s}_{6})}{8\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(S^{s}_{2})},$
(75) $\displaystyle P^{\prime}_{4}(q^{2}_{1},q^{2}_{2})$ $\displaystyle=$
$\displaystyle\frac{\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(S_{4})}{\sqrt{-\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(S^{c}_{2}S^{s}_{2})}},\quad
P^{\prime}_{5}(q^{2}_{1},q^{2}_{2})=\frac{\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(\beta_{\ell}S_{5})}{2\sqrt{-\int^{q^{2}_{2}}_{q^{2}_{1}}dq^{2}(S^{c}_{2}S^{s}_{2})}}.$
(76)

Table 11: Theoretical predictions for the $q^{2}$-binned observables $A^{\mu}_{FB}$, $P_{1,2}$ of the decays $B_{s}\to\phi\ell^{+}\ell^{-}$ with $\ell=(\mu,\tau)$ in the PQCD (the first low) and “PQCD+Lattice” (the second row) approaches. $q^{2}$ bin $(\rm GeV^{2})$ | $A^{\mu}_{FB}$ | $A^{\tau}_{FB}$ | $P_{1}(\ell=\mu)$ | $P_{1}(\ell=\tau)$ | $P_{2}(\ell=\mu)$ | $P_{2}(\ell=\tau)$
---|---|---|---|---|---|---
$[0.1,2.0]$ | ${0.131}^{+0.003}_{-0.003}$ | $-$ | ${0.015}^{+0.001}_{-0.001}$ | $-$ | ${0.206}^{+0.001}_{-0.001}$ | $-$
| ${0.122}^{+0.003}_{-0.002}$ | $-$ | ${0.013}^{+0.001}_{-0.001}$ | $-$ | ${0.204}^{+0.001}_{-0.001}$ | $-$
$[2.0,5.0]$ | ${-0.047}^{+0.004}_{-0.003}$ | $-$ | ${-0.197}^{+0.007}_{-0.007}$ | $-$ | ${-0.128}^{+0.003}_{-0.003}$ | $-$
| ${-0.038}^{+0.001}_{-0.002}$ | $-$ | ${-0.225}^{+0.005}_{-0.006}$ | $-$ | ${-0.132}^{+0.002}_{-0.002}$ | $-$
$[5.0,8.0]$ | ${-0.265}^{+0.006}_{-0.005}$ | $-$ | ${-0.280}^{+0.009}_{-0.008}$ | $-$ | ${-0.431}^{+0.001}_{-0.001}$ | $-$
| ${-0.200}^{+0.006}_{-0.006}$ | $-$ | ${-0.323}^{+0.003}_{-0.004}$ | $-$ | ${-0.427}^{+0.001}_{-0.001}$ | $-$
$[11.0,12.5]$ | ${-0.365}^{+0.001}_{-0.001}$ | $-$ | ${-0.440}^{+0.010}_{-0.008}$ | $-$ | ${-0.431}^{+0.002}_{-0.003}$ | $-$
| ${-0.298}^{+0.005}_{-0.005}$ | $-$ | ${-0.489}^{+0.002}_{-0.003}$ | $-$ | ${-0.420}^{+0.001}_{-0.001}$ | $-$
$[15.0,17.0]$ | ${-0.326}^{+0.001}_{-0.001}$ | ${-0.188}^{+0.001}_{-0.001}$ | ${-0.697}^{+0.009}_{-0.006}$ | ${-0.707}^{+0.009}_{-0.005}$ | ${-0.345}^{+0.003}_{-0.004}$ | ${-0.341}^{+0.003}_{-0.004}$
| ${-0.290}^{+0.005}_{-0.005}$ | ${-0.160}^{+0.003}_{-0.003}$ | ${-0.730}^{+0.007}_{-0.005}$ | ${-0.738}^{+0.006}_{-0.006}$ | ${-0.330}^{+0.003}_{-0.003}$ | ${-0.326}^{+0.003}_{-0.004}$
$[17.0,19.0]$ | ${-0.226}^{+0.002}_{-0.003}$ | ${-0.153}^{+0.001}_{-0.002}$ | ${-0.869}^{+0.005}_{-0.003}$ | ${-0.875}^{+0.005}_{-0.003}$ | ${-0.231}^{+0.002}_{-0.003}$ | ${-0.225}^{+0.003}_{-0.004}$
| ${-0.208}^{+0.004}_{-0.003}$ | ${-0.137}^{+0.003}_{-0.002}$ | ${-0.884}^{+0.004}_{-0.004}$ | ${-0.890}^{+0.005}_{-0.003}$ | ${-0.219}^{+0.003}_{-0.004}$ | ${-0.213}^{+0.004}_{-0.004}$
Table 12: Theoretical predictions for the $q^{2}$-binned optimized observables $P^{\prime}_{4}$ and $P^{\prime}_{5}$ of the decays $B_{s}\to\phi\ell^{+}\ell^{-}$ with $\ell=(e,\mu,\tau)$ in the PQCD (the first low) and “PQCD+Lattice” (the second row) approach.  $q^{2}$ bin $(\rm GeV^{2})$ | $P^{\prime}_{4}(\ell=e)$ | $P^{\prime}_{4}(\ell=\mu)$ | $P^{\prime}_{4}(\ell=\tau)$ | $P^{\prime}_{5}(\ell=e)$ | $P^{\prime}_{5}(\ell=\mu)$ | $P^{\prime}_{5}(\ell=\tau)$
---|---|---|---|---|---|---
$[0.1,2.0]$ | ${-0.322}^{+0.005}_{-0.005}$ | ${-0.289}^{+0.005}_{-0.005}$ | $-$ | ${0.484}^{+0.001}_{-0.002}$ | ${0.459}^{+0.001}_{-0.002}$ | $-$
| ${-0.318}^{+0.001}_{-0.001}$ | ${-0.285}^{+0.001}_{-0.001}$ | $-$ | ${0.475}^{+0.002}_{-0.001}$ | ${0.449}^{+0.001}_{-0.001}$ | $-$
$[2.0,5.0]$ | ${0.903}^{+0.007}_{-0.007}$ | ${0.903}^{+0.007}_{-0.007}$ | $-$ | ${-0.606}^{+0.001}_{-0.002}$ | ${-0.607}^{+0.001}_{-0.002}$ | $-$
| ${0.916}^{+0.001}_{-0.001}$ | ${0.916}^{+0.001}_{-0.001}$ | $-$ | ${-0.608}^{+0.004}_{-0.003}$ | ${-0.608}^{+0.004}_{-0.004}$ | $-$
$[5.0,8.0]$ | ${1.110}^{+0.003}_{-0.005}$ | ${1.110}^{+0.004}_{-0.005}$ | $-$ | ${-0.816}^{+0.005}_{-0.006}$ | ${-0.816}^{+0.005}_{-0.006}$ | $-$
| ${1.127}^{+0.002}_{-0.001}$ | ${1.127}^{+0.002}_{-0.001}$ | $-$ | ${-0.796}^{+0.003}_{-0.001}$ | ${-0.796}^{+0.002}_{-0.002}$ | $-$
$[11.0,12.5]$ | ${1.198}^{+0.002}_{-0.005}$ | ${1.198}^{+0.002}_{-0.005}$ | $-$ | ${-0.727}^{+0.005}_{-0.007}$ | ${-0.727}^{+0.005}_{-0.007}$ | $-$
| ${1.217}^{+0.001}_{-0.001}$ | ${1.217}^{+0.001}_{-0.001}$ | $-$ | ${-0.698}^{+0.001}_{-0.001}$ | ${-0.698}^{+0.002}_{-0.001}$ | $-$
$[15.0,17.0]$ | ${1.302}^{+0.002}_{-0.004}$ | ${1.302}^{+0.002}_{-0.004}$ | ${1.306}^{+0.002}_{-0.004}$ | ${-0.533}^{+0.005}_{-0.007}$ | ${-0.533}^{+0.005}_{-0.007}$ | ${-0.525}^{+0.006}_{-0.008}$
| ${1.314}^{+0.002}_{-0.003}$ | ${1.314}^{+0.002}_{-0.003}$ | ${1.317}^{+0.003}_{-0.002}$ | ${-0.506}^{+0.005}_{-0.006}$ | ${-0.506}^{+0.005}_{-0.006}$ | ${-0.498}^{+0.006}_{-0.006}$
$[17.0,19.0]$ | ${1.366}^{+0.001}_{-0.002}$ | ${1.366}^{+0.001}_{-0.002}$ | ${1.369}^{+0.001}_{-0.002}$ | ${-0.341}^{+0.004}_{-0.006}$ | ${-0.341}^{+0.004}_{-0.006}$ | ${-0.331}^{+0.004}_{-0.006}$
| ${1.371}^{+0.001}_{-0.001}$ | ${1.371}^{+0.001}_{-0.001}$ | ${1.374}^{+0.001}_{-0.002}$ | ${-0.323}^{+0.005}_{-0.006}$ | ${-0.323}^{+0.005}_{-0.006}$ | ${-0.314}^{+0.006}_{-0.006}$

From the numerical values as shown in Fig. 3 and in Table 8, 9 and 10, we find
the following points about the relevant physical observables of the considered
$B_{s}\to\phi\mu^{+}\mu^{-}$ decays in bins:

  1. (1)

For $B_{s}\to\phi\mu^{+}\mu^{-}$ decay, besides the good consistency between
the theory and the LHCb data for the integrated total branching ratio ${\cal
B}(B_{s}\to\phi\mu^{+}\mu^{-})$ as listed in Eq. (64), the PQCD and
“PQCD+Lattice” predictions for ${\cal B}(B_{s}\to\phi\mu^{+}\mu^{-})$ in most
bins do agree well with the measured ones within $2\sigma$ errors. For the
first low-$q^{2}$ bin $0.1<q^{2}<2$ $({\rm GeV^{2}})$, however, the central
value of the LHCb result $1.11\pm 0.16$ is larger than the theoretical ones by
roughly a factor of three. The LHCb results of ${\cal
B}(B_{s}\to\phi\mu^{+}\mu^{-})$ in different bins of $q^{2}$ as listed in the
third column of Table 9 are obtained from the results as given in Refs.
Aaij:2015esa ; Aaij:2021pkz  by multiplying the LHCb measured values of
differential decay rate $d{\cal B}(B_{s}\to\phi\mu^{+}\mu^{-})/dq^{2}$ with
the width of the corresponding bin $(q_{2}^{2}-q_{1}^{2})$. The theoretical
errors of our theoretical predictions of the branching ratios in bins are
still relatively large, while the differences between the PQCD and
“PQCD+Lattice” predictions for ${\cal B}(B_{s}\to\phi\ell^{+}\ell^{-})$ with
$\ell=(\mu,\tau)$ are small.

  2. (2)

In the first low-$q^{2}$ bin $0.1<q^{2}<2$ $({\rm GeV^{2}})$, both the PQCD
and ”PQCD+Lattice” predictions for $F_{L}^{\phi}(\ell=\mu)$ are larger than
the LHCb measured results $F_{L}^{\phi}(\ell=\mu)|_{{\rm
LHCb}}=0.20^{+0.08}_{-0.09}\pm 0.02$ Aaij:2015esa . For other bins, both PQCD
and “PQCD+Lattice” predictions of $F^{\phi}_{L}$ for muon mode do agree very
well with currently available LHCb measured values Aaij:2015esa  within
$2\sigma$ errors. It is worth of remaining that our theoretical predictions of
$F_{L}^{\phi}$ have a little error of $\sim 2\%$ due to the strong
cancellation of the theoretical errors in the ratios. The theoretical
predictions for ${\cal B}(B_{s}\to\phi\tau^{+}\tau^{-})$ and
$F_{L}^{\phi}(\ell=\tau)$ in different bins of $q^{2}$ as listed in Table 9
will be tested by future experimental measurements,

  3. (3)

For the observables $S_{3,4,7}$, as listed in Table 10, the PQCD and
“PQCD+Lattice” predictions for their values in all bins are in the range of
$10^{-3}-10^{-1}$, and show a good agreement with the LHCb measured values
Aaij:2015esa  . The errors of the theoretical predictions are also very small,
$\sim 2\%$ in magnitude, because of their nature of the ratios. In all bins,
the LHCb measured values of $S_{3,4,7}$ are still consistent with zero due to
their still large errors, which is a clear feature as can be seen easily from
the numerical values in Table 10 and the crosses in Fig. 3.

In Table 11 and 12, we show the PQCD and “PQCD+Lattice” predictions for the
physical observables $A_{FB}^{\mu,\tau},P_{1,2}^{\mu,\tau}$ and
$P^{\prime}_{4,5}(\ell=e,\mu,\tau)$ in six bins. These physical observables
could be tested in the near future LHCb and Belle-II experiments.

##  V Summary

In this paper, we made a systematic study of the semileptonic decays
$B_{s}\to\phi\ell^{+}\ell^{-}$ with $\ell^{-}=(e^{-},\mu^{-},\tau^{-})$ using
the PQCD and the “PQCD+Lattice” factorization approach respectively. We first
evaluated all relevant form factors in the low $q^{2}$ region using the PQCD
approach, and we also took currently available lattice QCD results at the
high-$q^{2}$ points $q^{2}=(12,16,18.9)$ GeV2 as additional input to improve
the extrapolation of the form factors from the low to the high-$q^{2}$ region.
We calculated the branching ratios ${\cal B}(B_{s}\to\phi\ell^{+}\ell^{-})$,
the CP averaged $\phi$ longitudinal polarization fraction $F_{L}(q^{2})$, the
forward-backward asymmetry ${\cal A}_{FB}(q^{2})$, the CP averaged angular
coefficients $S_{3,4,7}(q^{2})$, the CP asymmetry angular coefficients
$A_{5,6,8,9}(q^{2})$, the optimized observables $P_{1,2,3}(q^{2})$ and
$P^{\prime}_{4,5,6,8}(q^{2})$. For $B_{s}\to\phi\mu^{+}\mu^{-}$ decay mode, we
calculated the binned values of the branching ratio ${\cal
B}(B_{s}\to\phi\mu^{+}\mu^{-})$, the observables $F_{L}^{\phi}$ and
$S_{3,4,7}$ in the same bins as defined by LHCb Collaboration Aaij:2015esa  in
order to compare our theoretical predictions with those currently available
LHCb measurements bin by bin directly.

Based on the analytical evaluations, the numerical results and the
phenomenological analysis, we found the following main points:

  1. (1)

For the branching ratio ${\cal B}(B_{s}\to\phi\mu^{+}\mu^{-})$, the PQCD and
“PQCD+Lattice” prediction are $(7.07^{+3.43}_{-2.34})\times 10^{-7}$ and
$(6.76^{+1.52}_{-1.25})\times 10^{-7}$ respectively, which agree well with the
LHCb measured values Aaij:2015esa ; Aaij:2021pkz  and the QCDSR prediction
within still large errors. For the electron and tau mode, our theoretical
predictions for their decay rates are also well consistent with the
corresponding QCDSR predictions and to be tested by future experimental
measurements.

  2. (2)

For the ratios of the branching ratios $R_{\phi}^{e\mu}$ and
$R_{\phi}^{\mu\tau}$, the PQCD and “PQCD+Lattice” predictions agree with each
other and with small theoretical errors because of the strong cancellation of
the theoretical errors in such ratios. We suggest the LHCb and Belle-II
collaboration to measure these ratios.

  3. (3)

For the longitudinal polarization $F_{L}$, both PQCD and ”PQCD+Lattice”
predictions agree with the LHCb measurements in the considered bins within the
errors. For the CP averaged angular coefficients $S_{3,4,7}$, the PQCD and
”PQCD+Lattice” predictions in all bins are small in magnitude, in the range of
$10^{-3}-10^{-1}$, and agree well with the LHCb results within the still large
experimental errors. For the CP asymmetry angular coefficients $A_{5,6,8,9}$ ,
the PQCD and ”PQCD+Lattice” predictions are very small, in the range of
$10^{-4}-10^{-2}$, and clearly consistent with the LHCb measurements in the
six bins.

  4. (4)

For the physical observables $A_{FB}^{l}$, $P_{1,2,3}$ and
$P^{\prime}_{4,5,6,8}$, the experimental measurements are still absent now, we
think that the PQCD and “PQCD+Lattice” predictions for these physical
observables will be tested in the near future LHCb and Belle-II experiments.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
under Grant No. 11775117 and 11235005.

## References

  * (1) B. Dey [LHCb Collaboration], _Lepton Flavor Universality tests in $b\to sl^{+}l^{-}$ decays at LHCb_, PoS ICHEP 2018, 069 (2019).
  * (2) R. Aaij et al. [LHCb Collaboration], _Test of lepton universality using $B^{+}\to K^{+}l^{+}l^{-}$ decays_, Phys. Rev. Lett. 113, 151601 (2014).
  * (3) R. Aaij et al. [LHCb Collaboration], _Test of lepton universality with $B^{0}\to K^{*0}l^{+}l^{-}$ decays_, JHEP 1708, 055 (2017).
  * (4) M. Bordone, G. Isidori and A. Pattori, _On the Standard Model predictions for $R_{K}$ and $R_{K^{*}}$ _, Eur. Phys. J. C 76, 440 (2016).
  * (5) G. Hiller and F. Kruger, _More model-independent analysis of $b\to s$ processes _, Phys. Rev. D 69, 074020 (2004).
  * (6) G. Hiller and M. Schmaltz, _Diagnosing lepton-nonuniversalityin $b\to sll$ test _, JHEP 02, 055 (2015).
  * (7) S. Descotes-Genon, J. Matias and J. Virto, _Understanding the $B\to K^{*}\mu^{+}\mu^{-}$ anomaly_, Phys. Rev. D 88, 074002 (2013).
  * (8) J. Matias, F. Mescia, M. Ramon and J. Virto, _Complete Anatomy of $\bar{B}_{d}\to\bar{K}^{*0}(\to K\pi)l^{+}l^{-}$ and its angular distribution_, JHEP 1204, 104 (2012).
  * (9) J. Matias and N. Serra, _Symmetry relations between angular observables in $B^{0}\to K^{*}\mu^{+}\mu^{-}$ and the LHCb $P_{5}^{\prime}$ anomaly_, Phys. Rev. D 90, 034002 (2014).
  * (10) S. Descotes-Genon, J. Matias, M. Ramon and J. Virto, _Implications from clean observables for the binned analysis of $B\to K^{*}\mu^{+}\mu^{-}$ at large recoil_, JHEP 1301, 048 (2013).
  * (11) S. Descotes-Genon, T. Hurth, J. Matias and J. Virto, _Optimizing the basis of $B\to K^{*}ll$ observables in the full kinematic range _, JHEP 1305, 137 (2013).
  * (12) R. Aaij et al. [LHCb Collaboration], _Measurement of Form-Factor-Independent Observables in the Decay $B^{0}\to K^{*0}\mu^{+}\mu^{-}$_, Phys. Rev. Lett. 111, 191801 (2013).
  * (13) R. Aaij et al. [LHCb Collaboration], _Angular analysis of the $B^{0}\to K^{*0}\mu^{+}\mu^{-}$ decay using 3 fb-1 of integrated luminosity_, JHEP 1602, 104 (2016).
  * (14) A. Abdesselam et al. [Belle], _Angular analysis of $B^{0}\to K^{*}(892)^{0}\ell^{+}\ell^{-}$_, [arXiv:1604.04042 [hep-ex]].
  * (15) S. Descotes-Genon, L. Hofer, J. Matias and J. Virto, _On the impact of power corrections in the prediction of $B\to K^{*}\mu^{+}\mu^{-}$ observables_, JHEP 1412, 125 (2014).
  * (16) A. Deandrea and A. D. Polosa, _The Exclusive $B_{s}\to\phi\mu^{+}\mu^{-}$ process in a constituent quark model_, Phys. Rev. D 64, 074012 (2001).
  * (17) S. Dubnicka, A. Z. Dubnickova, A. Issadykov, M. A. Ivanov, A. Liptaj and S. K. Sakhiyev, _Decay $B_{s}\to\phi\ell^{+}\ell^{-}$ in covariant quark model_, Phys. Rev. D 93, 094022 (2016).
  * (18) C. Q. Geng and C. C. Liu, _Study of $B_{s}\to(\eta,\eta^{\prime},\phi)\ell\bar{\ell}$ decays_, J. Phys. G 29, 1103 (2003).
  * (19) C. Bobeth, G. Hiller and G. Piranishvili, _CP Asymmetries in $\bar{B}\to\bar{K}^{*}(\to\bar{K}\pi)\bar{\ell}\ell$ and Untagged $\bar{B}_{s}$, $B_{s}\to\phi(\to K^{+}K^{-})\bar{\ell}\ell$ Decays at NLO_, JHEP 0807, 106 (2008).
  * (20) P. Ball and R. Zwicky, _$B_{d,s}\to\rho,\omega,K^{*},\phi$ decay form-factors from light-cone sum rules revisited_, Phys. Rev. D 71, 014029 (2005).
  * (21) Y. L. Wu, M. Zhong and Y. B. Zuo, _$B_{s},D_{s}\to\pi,K,\eta,\rho,K^{*},\omega,\phi$ Transition Form Factors and Decay Rates with Extraction of the CKM parame ters$|V_{ub}|,|V_{cs}|,|V_{cd}|$_, Int. J. Mod. Phys. A 21, 6125 (2006).
  * (22) W. Altmannshofer and D. M. Straub, _New physics in $b\to s$ transitions after LHC run I _, Eur. Phys. J. C 75, 382 (2015).
  * (23) A. Bharucha, D. M. Straub and R. Zwicky, _ $B\to V\ell^{+}\ell^{-}$ in the Standard Model from light-cone sum rules_, JHEP 1608, 098 (2016).
  * (24) J. Gao, C.D. Lü, Y.L. Shen, Y.M. Wang and Y.B. Wei, _Precision calculations of $B\to V$ form factors from soft-collinear effective theory sum rules on the light-cone_, Phys. Rev. D 101, 074035 (2020).
  * (25) S. Descotes-Genon, L. Hofer, J. Matias and J. Virto, _Global analysis of $b\to s\ell\ell$ anomalies_, JHEP 1606, 092 (2016).
  * (26) R. Mohanta and A. K. Giri, _Study of FCNC mediated rare $B_{s}$ decays in a single universal extra dimension scenario_, Phys. Rev. D 75, 035008 (2007).
  * (27) Y. Li and J. Hua, _Study of $B_{s}\to\phi\ell^{+}\ell^{-}$ Decay in a Single Universal Extra Dimension_, Eur. Phys. J. C 71, 1764 (2011).
  * (28) Y. G. Xu, L. H. Zhou, B. Z. Li and R. M. Wang, _Analysis of $B_{s}\to\phi\mu^{+}\mu^{-}$ decay within supersymmetry_, Chin. Phys. C 37, 063104 (2013).
  * (29) T. Aaltonen et al. [CDF Collaboration], _Measurement of the Forward-Backward Asymmetry in the $B\to K^{(*)}\mu^{+}\mu^{-}$ Decay and First Observation of the $B^{0}_{s}\to\phi\mu^{+}\mu^{-}$ Decay_, Phys. Rev. Lett. 106, 161801 (2011).
  * (30) R. Aaij et al. [LHCb Collaboration], _Differential branching fraction and angular analysis of the decay $B_{s}^{0}\to\phi\mu^{+}\mu^{-}$ _, JHEP 1307, 084 (2013).
  * (31) R. Aaij et al. [LHCb Collaboration], _Angular analysis and differential branching fraction of the decay $B^{0}_{s}\to\phi\mu^{+}\mu^{-}$_, JHEP 1509, 179 (2015).
  * (32) R. Aaij et al. [LHCb Collaboration], _Precise measurement of the $f_{s}/f_{d}$ ratio of fragmentation fractions and of $B^{0}_{s}$ decay branching fractions_, Phys. Rev. D 104, 032005 (2021).
  * (33) R. Aaij et al. [LHCb Collaboration], _Branching fraction measurements of the rare $B^{0}_{s}\rightarrow\phi\mu^{+}\mu^{-}$ and $B^{0}_{s}\rightarrow f_{2}^{\prime}(1525)\mu^{+}\mu^{-}$ decays_, arXiv:2105.14007 [hep-ex].
  * (34) R. Aaij et al. [LHCb Collaboration], _Angular analysis of the rare decay $B^{0}_{s}\rightarrow\phi\mu^{+}\mu^{-}$_, arXiv:2107.13428 [hep-ex].
  * (35) S. P. Jin, X. Q. Hu and Z. J. Xiao, _Study of $B_{s}\to K^{(*)}\ell^{+}\ell^{-}$ decays in the PQCD factorization approach with lattice QCD input_, Phys. Rev. D 102, 013001 (2020).
  * (36) Y. Y. Kim, H.N. Li and A.I. Sanda, _Penguin enhancement and $B\to K\pi$ decays in perturbative QCD_, Phys. Rev. D 63, 054008 (2001).
  * (37) C. D. Lu, K. Ukai and M. Z. Yang, _Branching ratio and CP violation of $B\to\pi\pi$ decays in perturbative QCD approach_, Phys. Rev. D 63, 074009 (2001).
  * (38) H. N. Li, _QCD Aspects of Exclusive B Meson Decays_ , Prog.Part. $\&$ Nucl. Phys. 51, 85 (2003) and references therein.
  * (39) H. N. Li and S. Mishima, _Pion transition form factor in $k_{T}$ factorization_, Phys. Rev. D 80, 074024 (2009).
  * (40) Y. Y. Fan, W. F. Wang, S. Cheng and Z. J. Xiao, _Anatomy of $B\to K\eta^{(^{\prime})}$ decays in different mixing schemes and effects of next-to-leading order contributions in the perturbative QCD approach_, Phys. Rev. D 87, 094003 (2013).
  * (41) Y. Y. Fan, W. F. Wang, S. Cheng and Z. J. Xiao, _Semileptonic decays $B\to D^{(*)}l\nu$ in the perturbative QCD factorization approach_, Chin. Sci. Bull. 59, 125 (2014).
  * (42) A. Ali, G. Kramer, Y. Li, C. D. Lu, Y. L. Shen, W. Wang and Y. M. Wang, _Charmless non-leptonic $B_{s}$ decays to $PP$, $PV$ and $VV$ final states in the pQCD approach_, Phys. Rev. D 76, 074018 (2007).
  * (43) Z. J. Xiao, W. F. Wang and Y. Y. Fan, _Revisiting the pure annihilation decays $B_{s}\to\pi^{+}\pi^{-}$ and $B^{0}\to K^{+}K^{-}$: the data and the pQCD predictions_, Phys. Rev. D 85, 094003 (2012).
  * (44) W.F. Wang and Z.J. Xiao, _The semileptonic decays $B/B_{s}\to(\pi,K)(l^{+}l^{-},l\nu,\nu\bar{\nu})$ in the perturbative QCD approach beyond the leading-order_, Phys. Rev. D 86, 114025 (2012).
  * (45) W.F. Wang, Y.Y. Fan, M. Liu and Z.J. Xiao, _Semileptonic decays $B/B_{s}\to(\eta,\eta^{\prime},G)$ $(l^{+}l^{-},l\bar{\nu},\nu\bar{\nu})$ in the perturbative QCD approach beyond the leading order _, Phys. Rev. D 87, 097501 (2013).
  * (46) D.C. Yan, P. Yang, X. Liu, and Z.J. Xiao, _Anatomy of $B_{s}\to PV$ decays and effects of next-to-leading order contributions in the perturbative QCD factorization approach_, Nucl. Phys. B 931, 79 (2018).
  * (47) D.C. Yan, X. Liu, and Z.J. Xiao, _Anatomy of $B_{s}\to VV$ decays and effects of next-to-leading order contributions in the perturbative QCD factorization approach_, Nucl. Phys. B 935, 17 (2018).
  * (48) M. Döring, Ulf-G. Meißner and Wei Wang, _Chiral dynamics and S-wave contributions in semileptonic B decays_ , JHEP 10, 011 (2013).
  * (49) T. Kurimoto, H.N. Li and A.I. Sanda, _Leading power contributions to $B\to\pi,\rho$ transition form-factors_, Phys. Rev. D 65, 014007 (2002).
  * (50) C.D. Lu and M.Z. Yang, _B to light meson transition form factors calculated in perturbative QCD approach_ , Eur. Phys. J. C 28, 515 (2003).
  * (51) P. Ball, V. M. Braun and A. Lenz, _Higher-twist distribution amplitudes of the K meson in QCD_ , JHEP 0605, 004 (2006).
  * (52) G. Buchalla, A. J. Buras and M. E. Lautenbacher, _Weak decays beyond leading logarithms_ , Rev. Mod. Phys. 68, 1125 (1996).
  * (53) R. H. Li, C. D. Lu, W. Wang and X. X. Wang, _$B\to S$ Transition Form Factors in the PQCD approach_, Phys. Rev. D 79, 014013 (2009).
  * (54) S. R. Singh and B. Mawlong, _331- $Z^{\prime}$ mediated FCNC effects on $b\to d\mu^{+}\mu^{-}$ transitions_, Int. J. Mod. Phys. A 33, 1850225 (2019).
  * (55) B. Kindra and N. Mahajan, _Predictions of angular observables for $\bar{B}_{s}\to K^{*}\ell\ell$ and $\bar{B}\to\rho\ell\ell$ in the standard model_, Phys. Rev. D 98, 094012 (2018).
  * (56) P. Nayek, P. Maji and S. Sahoo, _Study of semileptonic decays $B\to\pi\ell^{+}\ell^{-}$ and $B\to\rho\ell^{+}\ell^{-}$ in nonuniversal Z’ model _, Phys. Rev. D 99, 013005 (2019).
  * (57) K. G. Chetyrkin, M. Misiak and M. Munz, _Weak radiative B meson decay beyond leading logarithms_ , Phys. Lett. B 400, 206 (1997).
  * (58) K. G. Chetyrkin, M. Misiak and M. Munz, _$|\Delta F|=1$ nonleptonic effective Hamiltonian in a simpler scheme_, Nucl. Phys. B 520, 279 (1998).
  * (59) P. Gambino, M. Gorbahn and U. Haisch, _Anomalous dimension matrix for radiative and rare semileptonic B decays up to three loops_ , Nucl. Phys. B 673, 238 (2003).
  * (60) C. Bobeth, M. Misiak and J. Urban, _Photonic penguins at two loops and $m_{t}$ dependence of $BR(B\to X_{s}l^{+}l^{-})$ _, Nucl. Phys. B 574, 291 (2000).
  * (61) C. H. Chen and C. Q. Geng, _Baryonic rare decays of $\Lambda(b)\to\Lambda l^{+}l^{-}$_, Phys. Rev. D 64, 074001 (2001).
  * (62) C. S. Lim, T. Morozumi and A. I. Sanda, _A Prediction for $d\Gamma(b\to s\ell\bar{\ell})/dq^{2}$ Including the Long Distance Effects_, Phys. Lett. B 218, 343 (1989).
  * (63) N. G. Deshpande, J. Trampetic and K. Panose, _Resonance Background to the Decays $b\to s\ell^{+}\ell^{-}$, $B\to K^{*}\ell^{+}\ell^{-}$ and $B\to K\ell^{+}\ell^{-}$ _, Phys. Rev. D 39, 1461 (1989).
  * (64) A. Ali, T. Mannel and T. Morozumi, _Forward backward asymmetry of dilepton angular distribution in the decay $b\to s\ell^{+}\ell^{-}$_, Phys. Lett. B 273, 505 (1991).
  * (65) P. J. O’Donnell and H. K. K. Tung, _Resonance contributions to the decay $b\to s\ell^{+}\ell^{-}$_, Phys. Rev. D 43, 2067 (1991).
  * (66) A. Khodjamirian, Th. Mannel, A.A. Pivovarov and Y.M. Wang, _Charm-loop effect in $B\to K^{(*)}l^{+}l^{-}$ and $B\to K^{*}\gamma$_, JHEP 09, 089 (2010).
  * (67) A. Khodjamirian, Th. Mannel, A.A. Pivovarov and Y.M. Wang, _ $B\to Kl^{+}l^{-}$ decay at large hadronic recoil_, JHEP 02, 010 (2013).
  * (68) I. Ahmed, M. J. Aslam and M. Ali Paracha, _Impact of $Z^{\prime}$ and universal extra dimension parameters on different asymmetries in $B_{s}\to\phi\ell^{+}\ell^{-}$ decays_, Phys. Rev. D 88, 014019 (2013).
  * (69) W. F. Wang, H. N. Li, W. Wang and C. D. Lu, _$S$ -wave resonance contributions to the $B^{0}_{(s)}\to J/\psi\pi^{+}\pi^{-}$ and $B_{s}\to\pi^{+}\pi^{-}\mu^{+}\mu^{-}$ decays_, Phys. Rev. D 91, 094024 (2015).
  * (70) W. Altmannshofer, P. Ball, A. Bharucha, A. J. Buras, D. M. Straub and M. Wick, _Symmetries and Asymmetries of $B\to K^{*}\mu^{+}\mu^{-}$ Decays in the Standard Model and Beyond_, JHEP 0901, 019 (2009).
  * (71) D. Becirevic and E. Schneider, _On transverse asymmetries in $B\to K^{*}l^{+}l^{-}$_, Nucl. Phys. B 854, 321 (2012).
  * (72) U. Egede, T. Hurth, J. Matias, M. Ramon and W. Reece, _New physics reach of the decay mode $\bar{B}\to\bar{K}^{*0}\ell^{+}\ell^{-}$_, JHEP 1010, 056 (2010).
  * (73) U. Egede, T. Hurth, J. Matias, M. Ramon and W. Reece, _New observables in the decay mode $\bar{B}_{d}\to\bar{K}^{*0}\ell^{+}\ell^{-}$_, JHEP 0811, 032 (2008).
  * (74) M. Tanabashi et al. [Particle Data Group], _Review of Particle Physics_ , Phys. Rev. D 98, 030001 (2018).
  * (75) Y.M. Wang and Y.L. Shen, _QCD corrections to $B\to\pi$ form factors from light-cone sum rules_, Nucl. Phys. B 898, 563 (2015).
  * (76) F. Kruger and J. Matias, _Probing new physics via the transverse amplitudes of $B^{0}\to K^{*0}(\to K^{-}\pi^{+})l^{+}l^{-}$ at large recoil_, Phys. Rev. D 71, 094009 (2005)
  * (77) J. Gratrexa, M. Hopfera and R. Zwicky, _Generalised helicity formalism, higher moments and the $B\to K_{J_{K}}(\to K\pi)\bar{l}_{1}l_{2}$ angular distributions_, Phys. Rev. D 93, 054008 (2016).
  * (78) H. Y. Cheng, C. K. Chua and C. W. Hwang, _Covariant light-front approach for s-wave and p-wave mesons: its application to decay constants and form factors_ , Phys. Rev. D 69, 074025 (2004).
  * (79) W. Wang , Y. L. Shen and C. D. Lu, _Covariant light-front approach for $B_{c}$ transition form factors_, Phys. Rev. D 79, 054012 (2009).
  * (80) Y. Y. Fan, Z. J. Xiao, R. M. Wang and B. Z. Li, _The $B\to D^{(*)}l\nu_{l}$ decays in the pQCD approach with the Lattice QCD input_, Sci. Bull. 60 (2015) 2009-2015.
  * (81) X. Q. Hu, S. P. Jin and Z. J. Xiao, _Semileptonic decays $B_{c}\to(\eta_{c},J/\psi)l\bar{\nu}_{l}$ in the “PQCD + Lattice” approach_, Chin. Phys. C 44  (2020) 023104.
  * (82) X. Q. Hu, S. P. Jin and Z. J. Xiao, _Semileptonic decays $B/B_{s}\to(D^{(*)},D_{s}^{(*)})l\nu_{l}$ in the PQCD factorization approach with the lattice QCD input_, Chin. Phys. C 44  (2020) 053102.
  * (83) C. Bourrely, I. Caprini and L. Lellouch, Model-independent description of $B\to\pi l\nu$ decays and a determination of $|V_{ub}|$, Phys. Rev. D 79, 013008 (2009). [Erratum ibid. 82 , 099902 (2010)].
  * (84) D. Leljak, B. Melic and M. Patra, On lepton flavour universality in semileptonic $B_{c}\to\eta_{c},J/\Psi$ decays, JHEP 05, 094 (2019).
  * (85) R. H. Li, C. D. Lu and W. Wang, _Transition form factors of B decays into p-wave axial-vector mesons in the perturbative QCD approach_ , Phys. Rev. D 79, 034014 (2009).
  * (86) U. O. Yilmaz, _Analysis of $B_{s}\to\phi\ell^{+}\ell^{-}$ decay with new physics effects_, Eur. Phys. J. C 58, 555 (2008).
  * (87) Y. Q. Peng and M. Z. Yang, _Study of semileptonic decay of $\bar{B}_{s}^{0}\to\phi l^{+}l^{-}$ in QCD sum rule_, arXiv:2001.08459 [hep-ph].
  * (88) D. Melikhov and B. Stech, _Weak form-factors for heavy meson decays: An Update_ , Phys. Rev. D 62, 014006 (2000).
  * (89) R. N. Faustov and V. O. Galkin, _Rare $B_{s}$ decays in the relativistic quark model_, Eur. Phys. J. C 73, 2593 (2013)
  * (90) C. D. Lu, W. Wang and Z. T. Wei, _Heavy-to-light form factors on the light cone_ , Phys. Rev. D 76, 014013 (2007).
  * (91) F. Su, Y. L. Wu, C. Zhuang and Y. B. Yang, _Charmless $B_{s}\to PP,PV,VV$ Decays Based on the Six-Quark Effective Hamiltonian with Strong Phase Effects II_, Eur. Phys. J. C 72, 1914 (2012).
  * (92) R. R. Horgan, Z. Liu, S. Meinel and M. Wingate, _Lattice QCD calculation of form factors describing the rare decays $B\to K^{*}\ell^{+}\ell^{-}$ and $B_{s}\to\phi\ell^{+}\ell^{-}$_, Phys. Rev. D 89, 094501 (2014).
  * (93) C. D. Lu, Y. L. Shen, Y. M. Wang and Y. B. Wei, _QCD calculations of $B\to\pi,K$ form factors with higher-twist corrections_, JHEP 1901, 024 (2019).
  * (94) A. Bharucha, T. Feldmann and M. Wick, _Theoretical and Phenomenological Constraints on Form Factors for Radiative and Semi-Leptonic B-Meson Decays_ , JHEP 1009, 090 (2010).
  * (95) D. Bigi and P. Gambino, _Revisiting $B\to D\ell\nu$_, Phys. Rev. D 94, 094008 (2016)
  * (96) B. Aubert et al. [BaBar Collaboration], _Direct CP, Lepton Flavor and Isospin Asymmetries in the Decays $B\to K^{(*)}\ell^{+}\ell^{-}$_, Phys. Rev. Lett. 102, 091803 (2009)
  * (97) S. Descotes-Genon and J. Virto, _Time dependence in $B\to Vll$ decays_, JHEP 1504, 045 (2015), Erratum: [ JHEP 1507, 049 (2015) ].