\section{Introduction}\label{intro}
Flavour-changing neutral-current (FCNC) processes involving $b\to s \ell^+\ell^-$, quark level transitions can play
a pivotal role in the indirect searches of physics beyond the Standard Model (SM). These transitions are CKM and
loop suppressed within the SM and therefore have high sensitivity to potential new physics (NP) effects.
Interestingly, recent experimental data on neutral current decays induced by $b\to s \ell^+\ell^-$
transitions have pointed towards several observables in tension with the SM predictions. Due to this fact, these
transitions currently stand among the most promising indications of NP.


The reported observables can be grouped into two sets: $b\to s \mu^+\mu^-$ observables that include only muons,
called as lepton-flavour dependent (LFD) observables, and the other known as lepton-flavour-universality violating
(LFUV) observables that involve both muons and electrons. The set of LFD observables contains several angular
observables including the so-called $P_5^{\prime}$ anomaly \cite{Descotes-Genon:2013vna,Descotes-Genon:2013wba},
first observed by the LHCb collaboration in the angular analysis of the $B^0\to K^{\ast0}\mu^+\mu^-$ decay
\cite{Aaij:2013qta,Aaij:2015oid}, later measured by Belle \cite{Abdesselam:2016llu,Wehle:2016yoi}, and more recently
by ATLAS \cite{Aaboud:2018krd} as well as CMS collaboration \cite{Sirunyan:2017dhj}. LHCb \cite{Aaij:2015oid} and ATLAS
\cite{Aaboud:2018krd} measurements for the angular observable $P_5^{\prime}$ in the low dimuon invariant mass-squared
range ($4.0\le q^2\le 6.0$ GeV$^2$) indicate $\approx 3.3\sigma$ \cite{Aebischer:2018iyb} discrepancy with respect to the
SM prediction. Similarly, Belle results deviate by $2.6\sigma$ \cite{Wehle:2016yoi} from the SM in the $4.0\le q^2\le 8.0$
GeV$^2$ bin, whereas CMS measurement is found to be consistent with the SM expectation within $1\sigma$ \cite{Sirunyan:2017dhj}.
Very recently \cite{Aaij:2020nrf}, LHCb collaboration has updated results for the angular observables with 4.7 fb$^{-1}$ of data, where
for the $P_5^{\prime}$ observable, a local discrepancy of $2.5$ and $2.9\sigma$ is reported in the $[4, 6]$ and $[6, 8]$ GeV$^2$ bins, respectively.
Adding to LFD observables, LHCb has also observed the branching fractions of the $B\to K\mu^+\mu^-$ \cite{Aaij:2014pli}, $B\to K^{\ast}\mu^+\mu^-$ \cite{Aaij:2014pli,Aaij:2013iag,Aaij:2016flj}, and $B_s\to \phi\mu^+\mu^-$ \cite{Aaij:2013aln,Aaij:2015esa} decays to be on
the lower side compared to their SM estimates. These LFD observables, while being sensitive to NP
\cite{Altmannshofer:2008dz,Bobeth:2011nj,Matias:2012xw,DescotesGenon:2012zf,Matias:2014jua}, can not establish the NP
case unambiguously because of the involvement of the hadronic uncertainties originating from the different long-distance
effects, in particular from form factors, power corrections, and charm resonances \cite{Khodjamirian:2010vf,Khodjamirian:2012rm,
Lyon:2014hpa,Descotes-Genon:2014uoa,Capdevila:2017ert,Blake:2017fyh}. Therefore, without having
additional data or a complete and reliable calculations of the hadronic uncertainties there remains a possibility to explain
the currently observed LFD anomalies with more conservative assumptions on the involved hadronic contributions
\cite{Jager:2012uw,Jager:2014rwa,Ciuchini:2015qxb,Ciuchini:2016weo,Bobeth:2017vxj}.

The second set with LFUV observables includes the ratios of branching fractions involving both $b\to s \mu^+\mu^-$ and
$b\to s e^+e^-$ transitions. The LHCb collaboration, in 2014, measured the ratio $R_K\equiv \mathcal{B}(B^+\to K^+\mu^+\mu^-)/(B^+\to K^+e^+e^-)=0.745_{-0.074}^{+0.090}\pm0.036$ \cite{Aaij:2014ora}, in the low dilepton invariant mass-squared range
($1.0\le q^2\le 6.0$ GeV$^2$), which showed $2.6\sigma$ deviation from the SM expectation of $\simeq1$
\cite{Hiller:2003js,Bordone:2016gaq}. Recently, LHCb announced new value of $R_K$ \cite{Aaij:2019wad}, by first
reanalyzing the Run I data and then combining it with the Run II result,
\begin{eqnarray}\label{01}
R_K^{[1.1,\,6]}=0.846_{-0.054}^{+0.060}(\text{stat})_{-0.014}^{+0.016}(\text{syst}).
\end{eqnarray}
Although the central value is now closer to the SM prediction, the discrepancy is still at the level of $2.5\sigma$ due to the
smaller uncertainty of the new measurement. Additionally, we have more recent Belle results of $R_K$ \cite{Abdesselam:2019lab},
which are combined together for the charged and neutral decay modes, and are presented in four different bins. However, due to
large errors, these results are consistent with the SM as well as with the LHCb data in the dilepton mass-squared range $1$ to $6$ GeV$^2$.
Similarly, in 2017, LHCb reported another ratio $R_{K^{\ast}}\equiv (B^0\to K^{\ast0}\mu^+\mu^-)/(B^0\to K^{\ast0}e^+e^-)$
\cite{Aaij:2017vbb}, in two different regions of $q^2$, the dilepton invariant mass squared. These measurements in the low $q^2$
($0.045\le q^2\le 1.1$ GeV$^2$) and central $q^2$ ($1.1\le q^2\le 6.0$ GeV$^2$) bin were found to deviate from the SM
predictions $\simeq1$ \cite{Hiller:2003js,Bordone:2016gaq}, with significance of $\approx2.4\sigma$ and $\approx2.5\sigma$,
respectively. Recently, Belle collaboration has also presented their results for $R_{K^{\ast}}$ \cite{Abdesselam:2019wac}, in multiple
$q^2$ bins, by considering both the neutral and charged $B$ modes. Again having large uncertainties, their results are in
agreement with both the SM and the LHCb measurements. Moreover, additional LFUV observables, such as
$Q_{4,5}=P_{4,5}^{\mu\prime}-P_{4,5}^{e\prime}$ \cite{Capdevila:2016ivx}, have been observed by the Belle collaboration
\cite{Wehle:2016yoi}. Although currently not statistically significant, the data indicate towards potential LFUV.
Furthermore, LHCb has also performed the first test of lepton flavour universality (LFU) violation with $b$ baryons by measuring
$R_{pK}$ ratio in $\Lambda_b \to p K^-\ell^+\ell^-$ decays \cite{Aaij:2019bzx}.
Contrary to the LFD observables, SM predictions for the LFUV observables $R_K$ and $R_{K^{\ast}}$ are theoretically
clean as the hadronic uncertainties essentially cancel and therefore they hold the key to unravel NP without ambiguity.

Because all (LFD+LFUV) observables involve $b\to s \mu^+\mu^-$ transition and additionally LFUV NP contributions are mandatory
to account for the LFUV anomalies, it seems natural to address all the $b\to s \ell^+\ell^-$ data by assuming NP present only in
the $b\to s \mu^+\mu^-$ sector, which automatically breaks lepton flavour universality and is thus purely LFUV. Interestingly,
several model independent global fit analyses \cite{Descotes-Genon:2015uva,Altmannshofer:2017fio,Alok:2017sui,Altmannshofer:2017yso,
Geng:2017svp,Ciuchini:2017mik,Capdevila:2017bsm,Alguero:2019ptt,Alok:2019ufo,Ciuchini:2019usw,Datta:2019zca,Aebischer:2019mlg,
Kowalska:2019ley,Arbey:2019duh} performed with this purely LFUV NP assumption have pointed out two simple
one-dimensional (1D) NP scenarios (S1) $C_{9\mu}^{\text{NP}}$ or (S2) $C_{9\mu}^{\text{NP}}=-C_{10\mu}^{\text{NP}}$, that can provide
better fit to all the $b\to s \ell^+\ell^-$ data with preferences reaching $\approx5-6\sigma$ compared to the SM. However, performing the separate fits,
it is observed that the inclusion of the latest $R_K$ and $R_{K^{\ast}}$ data has created tensions between the separate fit to LFD and
LFUV set in both the scenarios S1 and S2 along with increasing the significance gap between the two LFUV fits of the two scenarios
\cite{Alguero:2019ptt,Datta:2019zca}. These tensions, if not statistical fluctuations, could be indications of NP also present
in $b\to s e^+e^-$. For instance, in Ref. \cite{Datta:2019zca}, it is shown that the additional LFUV NP in $b\to s e^+e^-$
along with the basic scenarios S1 and S2 leads to a number of new scenarios, which can remove tensions along with improving the
overall fit to all data. In this case, LFD set is explained by LFUV NP only in $b\to s \mu^+\mu^-$, while the LFUV processes
are constrained by the combination of LFUV NP contributions in $b\to s \mu^+\mu^-$ and $b\to s e^+e^-$. Another complementary
approach proposed in \cite{Alguero:2018nvb}, before the latest $R_K$ and $R_{K^{\ast}}$ measurement, is to consider LFUV NP which
affects only muons along with LFU NP (with equal contributions to $e$, $\mu$, and $\tau$). In this scheme, LFD observables incorporate
contributions from LFUV+LFU NP, while LFUV set is explained mainly by the LFUV NP (while also having subleading interference with LFU NP contributions).
It was shown that several scenarios with both LFU and LFUV NP contributions to the Wilson
coefficients $C_{9\ell}^{\text{NP}}$ and $C_{10\ell}^{\text{NP}}$ can improve the agreement with the overall data. An updated global fit analysis,
including the latest experimental data, with both LFU and LFUV NP contributions extended to NP scenarios with right-handed currents (RHC),
is performed in \cite{Alguero:2019ptt}. Again, various NP scenarios are found with enhanced significance and an improved preference
for the NP scenarios with RHC is suggested to emerge after the inclusion of the latest $R_K$ and $R_{K^{\ast}}$ data. Furthermore, a better
description of data can also be obtained by increasing the degrees of freedom, i.e., 2D fits, along with the assumptions
such as NP affects only muons.

As there is no unique solution and many new scenarios are piling up due to recently emerging NP patterns in global fit
analyses, it is particularly important to discriminate between different possible scenarios and to devise methods to further
confirm or constrain patterns of NP. For that, on the experimental side, one requires more precise determinations of the already
measured observables, with different experimental setups, along with checking new additional observables. While on the theory side,
one can improve the calculations of hadronic contributions or construct new observables, with limited sensitivity to hadronic
uncertainties, either in the already existing decays or in the complementary processes induced by the same quark level
$b\to s \ell^+\ell^-$ transitions. In this work, we consider theoretically clean LFUV observables in complementary
channels and kinematical regions, and analyze the implications of the different NP scenarios proposed in the recent global analyses.
The list of decay channels induced by the $b\to s \ell^+\ell^-$ transition is long, and the LFUV ratios in a number of
decay channels have been studied \cite{Wang:2017mrd,Huang:2018rys,Dutta:2019wxo} based on the previous data. In the present study, we consider the most recent experimental results and restrict to seven exclusive channels
$M_{in}\to M_{f} \ell^+\ell^-$, with $M_{in}=B, B_s, \Lambda_b$ and
$M_{f}=f_0, K_0^{\ast}, K, K^{\ast}, \phi, K_1, \Lambda$. We study these channels in the model independent effective Hamiltonian
approach by employing the helicity formalism. The theoretical analysis of LFUV observables in complementary
hadronic decays can help to distinguish and further strengthen the emerging NP patterns along with
providing interesting cross checks.

The paper is organized as follows. In section \ref{effH}, we describe the effective Hamiltonian for $b\to s \ell^+\ell^-$ transitions.
In section \ref{matelem}, we discuss form factors present in the hadronic matrix elements for different decays along with their decay
amplitudes. In section \ref{secHelicity}, we consider the helicity formalism and work out the explicit expressions of the helicity amplitudes
for the considered decays. In section \ref{physobs}, we construct LFUV observables. Section \ref{secNum}, is devoted to numerical analysis, where
we also present our choice of the NP scenarios from different global fit analyses. The results are
summarized in section \ref{summary}. Details of the SM Wilson coefficients, numerical inputs, and the technical
considerations concerning the kinematics and polarization conventions, along with the predictions of the LFUV ratios
in the SM and the NP scenarios are collected in appendices.
\section{Effective Hamiltonian}\label{effH}
The effective weak Hamiltonian for $b\to s\ell^{+}\ell^{-}$ transition is given by
\begin{align}\label{H1}
\mathcal{H}_{\text{eff}}=-\frac{4 G_{F}}{\sqrt{2}}V_{tb}V^{\ast}_{ts}&\Bigg[\sum_{i=1}^{6}C_{i}O_{i}
+\sum_{i=7}^{8}\big(C_{i}O_{i}+C_{i^{\prime}}O_{i^{\prime}}\big)\notag\\
&+\sum_{i=9,10}\Big((C_{i}+C_{i\ell}^{\text{NP}})O_{i}+C_{i^{\prime}\ell}^{\text{NP}}O_{i^{\prime}}\Big)\Bigg],
\end{align}
where we have neglected the doubly Cabibbo suppressed contribution $(\propto V_{ub}V_{us}^\ast)$, and
$G_F$ is the Fermi coupling constant. The operators $O_{i\le6}$ are the same as the $P_{1,2}^c, P_{3,...,6}^{}$,
given in Ref. \cite{Bobeth:1999mk}, and the others are
\begin{align}\label{op1}
O_{7} &=\frac{e}{16\pi ^{2}}m_{b}\left( \bar{s}\sigma _{\mu \nu }P_{R}b\right) F^{\mu \nu },
&  O_{7^{\prime}} &=\frac{e}{16\pi ^{2}}m_{b}\left( \bar{s}\sigma _{\mu \nu }P_{L}b\right) F^{\mu \nu },\notag\\
O_{8} &=\frac{g_s}{16\pi ^{2}}m_{b}\left( \bar{s}\sigma _{\mu \nu } T^a P_{R}b\right) G^{\mu \nu\,a},
&  O_{8^{\prime}} &=\frac{g_s}{16\pi ^{2}}m_{b}\left( \bar{s}\sigma _{\mu \nu } T^a P_{L}b\right) G^{\mu \nu\,a},\notag\\
O_{9} &=\frac{e^{2}}{16\pi ^{2}}(\bar{s}\gamma _{\mu }P_{L}b)(\bar{l}\gamma^{\mu }l),
&  O_{9^{\prime}} &=\frac{e^{2}}{16\pi ^{2}}(\bar{s}\gamma _{\mu }P_{R}b)(\bar{l}\gamma^{\mu }l),\notag\\
O_{10} &=\frac{e^{2}}{16\pi ^{2}}(\bar{s}\gamma _{\mu }P_{L}b)(\bar{l} \gamma ^{\mu }\gamma _{5} l),
&  O_{10^{\prime}} &=\frac{e^{2}}{16\pi ^{2}}(\bar{s}\gamma _{\mu }P_{R}b)(\bar{l} \gamma ^{\mu }\gamma _{5} l),
\end{align}
where $e$ $(g_s)$ is the electromagnetic (strong) coupling constant, and $m_b$ represents the running
$b-$quark mass in the $\overline{\text{MS}}$ scheme.

Within the SM, major role in $b\to s\ell^{+}\ell^{-}$ transition, is played by operators $O_{7,9,10}$, whereas contributions of
primed dipole operators $O_{7^{\prime},8^{\prime}}$ are suppressed by $m_s/m_b$, and therefore we neglect them. Furthermore,
the factorizable contributions from current-current, QCD penguins and chromomagnetic dipole operators $O_{1-6,8}$ can be absorbed
into the effective Wilson coefficients $C_{7}^{\text{eff}}(q^{2})$ and $C_{9}^{\text{eff}}(q^{2})$ \cite{Du:2015tda}.
The explicit expressions of these Wilson coefficients, which we used, are presented in appendix \ref{append}. It is important to mention that, in
Eq. (\ref{H1}), we have considered NP contributions only in $O_{9^{(\prime)}}$ and $O_{10^{(\prime)}}$ operators because the emerging
viable NP solutions from the global fits of all the $b\to s\ell^{+}\ell^{-}$ data, which we consider in our study, are only
in the form of vector and axial-vector operators. The numerical values of Wilson coefficients at $\mu\sim m_{b}$
in the SM are presented in table \ref{wc table}.
\begin{table*}[t]
\setlength\tabcolsep{3.5pt}
\centering
\begin{tabular*}{.9\linewidth}{|cccccccccc|}
\hline
$C_{1}$&$C_{2}$&$C_{3}$&$C_{4}$&$C_{5}$&$C_{6}$&$C_{7}$&$C_{8}$&$C_{9}$&$C_{10}$
\\ \hline
  $-0.294$ &   $1.017$  & $-0.0059$  &   $-0.087$  &
  $0.0004$  &  $0.0011$   &   $-0.324$  &  $-0.176$  &
    $4.114$  &  $-4.193$ \\
\hline
\end{tabular*}
\caption{The SM Wilson coefficients $C_{i}^{\mu}$ up to NNLL accuracy given at the scale $\mu\sim m_{b}$.}
\label{wc table}
\end{table*}

Using the effective Hamiltonian given in Eq. (\ref{H1}), the decay amplitude for the process $M_{in}\to M_{f}\ell^{+}\ell^{-}$,
including the SM and the NP contributions, can be written as
\begin{eqnarray}\label{Amp1}
\mathcal{M}(M_{in}\to M_{f}\ell^{+}\ell^{-})&=&\frac{G_{F}\alpha}{2\sqrt{2}\pi}V_{tb}V^{\ast}_{ts}
\Bigg[\Big\{(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}})\langle M_{f}(k)|\bar s\gamma_{\mu}(1-\gamma_{5})b|M_{in}(p)\rangle\notag
\\
&+&C_{9^\prime\ell}^{\text{NP}}\langle M_{f}(k)|\bar s\gamma_{\mu}(1+\gamma_{5})b|M_{in}(p)\rangle\notag\\
&-&\frac{2m_{b}}{q^{2}}C_{7}^{\text{eff}}
\langle M_{f}(k)|\bar s i\sigma_{\mu\nu}q^{\nu}(1+\gamma_{5})b|M_{in}(p)\rangle\Big\}\big(\bar\ell\gamma^{\mu}\ell\big)\notag\\
&+&\Big\{(C_{10}+C_{10\ell}^{\text{NP}})\langle M_{f}(k)|\bar s\gamma_{\mu}(1-\gamma_{5})b|M_{in}(p)\rangle\notag\\
&+&C_{10^\prime\ell}^{\text{NP}}\langle M_{f}(k)|\bar s\gamma_{\mu}(1+\gamma_{5})b|M_{in}(p)\rangle\Big\}
\big(\bar\ell\gamma^{\mu}\gamma_{5}\ell\big)\Bigg].
\end{eqnarray}
\section{Hadronic matrix elements and decay amplitudes}\label{matelem}
The hadronic matrix elements in Eq. (\ref{Amp1}) can be parameterized in terms of different form factors. As we consider various decay
channels with final state hadron $(M_f)$ having scalar, pseudoscalar, vector, axial vector and spin half $(\Lambda)$ hadron,
so we consider transition form factors and decay amplitudes for each case separately.
\subsection{$M_{in}\to S\ell^{+}\ell^{-}$ case}
The matrix elements for the process $M_{in}\to S\ell^{+}\ell^{-}$, where the parent particle $M_{in}=B_s$ or $B$, and
the daughter particle $S$ is a scalar meson $0^{+}$, such as $S=f_{0}(980)$ or $K^{\ast}_0(1430)$, are given by
\begin{eqnarray}
\langle f_0\big(K^{\ast}_0\big)(k)|\bar s\gamma_{\mu}\gamma_{5}b|B_s\big(B\big)(p)\rangle&=&
-i\Big[f_{+}^{{f_0}({K^{\ast}_0})}(q^{2})P_{\mu}+f_{-}^{{f_0}({K^{\ast}_0})}(q^{2})q_{\mu}\Big],\label{2.2.1}\\
\langle f_0\big(K^{\ast}_0\big)(k)|\bar s i \sigma_{\mu\nu}q^{\nu}\gamma_{5}b|B_s\big(B\big)(p)\rangle&=& -i\frac{f_{T}^{{f_0}({K^{\ast}_0})}(q^{2})}{\big(m_{{B_s}(B)}+m_{{f_0}({K^{\ast}_0})}\big)}\notag
\\
&&\times\Big[q^{2}P_{\mu}-\Big(m^{2}_{{B_s}({B})}-m^{2}_{{f_0}({K^{\ast}_0})}\Big)q_{\mu}\Big],\label{2.2.2}
\end{eqnarray}
where $P_{\mu}= p_{\mu}+k_{\mu}$, and $q_{\mu}=p_{\mu}-k_{\mu}$. Using these matrix elements in Eq. (\ref{Amp1}), the decay amplitude for $M_{in}\to S\ell^{+}\ell^{-}$, can be written as
\begin{eqnarray}
\mathcal{M}\left(M_{in}\to S\ell^{+}\ell^{-}\right)=\frac{G_{F}\alpha}{2\sqrt{2}\pi}V_{tb}V^{\ast}_{ts}\{T^{1,S}_{\mu}(\bar{\ell}\gamma^{\mu}\ell)
+T^{2,S}_{\mu}(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell)\},\label{H1a}
\end{eqnarray}
where $T^{i,S}_{\mu}$ $(i=1, 2)$, incorporate both Wilson coefficients and transition form factors
\begin{eqnarray}
T^{i,{f_0}({K^{\ast}_0})}_{\mu}=iP_{\mu}\mathcal{F}^{i,{f_0}({K^{\ast}_0})}_{+}(q^2)+iq_{\mu}\mathcal{F}^{i,{f_0}
({K^{\ast}_0})}_{-}(q^2),\label{He6}
\end{eqnarray}
\begin{align}
\mathcal{F}^{1,{f_0}({K^{\ast}_0})}_{+}(q^2)&=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
f_{+}^{{f_0}({K^{\ast}_0})}(q^{2})\notag\\
&\qquad\qquad+\frac{2m_{b}}{\big(m_{{B_s}(B)}+m_{{f_0}({K^{\ast}_0})}\big)}
C_{7}^{\text{eff}}f_{T}^{{f_0}({K^{\ast}_0})}(q^{2}),\label{He7}\\
\mathcal{F}^{1,{f_0}({K^{\ast}_0})}_{-}(q^2)&=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
f_{-}^{{f_0}({K^{\ast}_0})}(q^{2})\notag\\
&\qquad\qquad-\frac{2m_{b}}{q^{2}}\big(m_{{B_s}(B)}-m_{{f_0}({K^{\ast}_0})}\big)
C_{7}^{\text{eff}}f_{T}^{{f_0}({K^{\ast}_0})}(q^{2}),\label{He8}\\
\mathcal{F}^{2,{f_0}({K^{\ast}_0})}_{\pm}(q^2)&=(C_{10}+C_{10\ell}^{\text{NP}}
-C_{10^{\prime}\ell}^{\text{NP}})f_{\pm}^{{f_0}({K^{\ast}_0})}(q^{2}).\label{He2}
\end{align}
For the $B_s\to f_{0}(980)\ell^{+}\ell^{-}$ decay, $f_+^{f_0}, f_0^{f_0}, f_T^{f_0}$ form factors are used in the numerical analysis.
For that $f_{-}^{f_0}(q^{2})$ form factor can be expressed as
\begin{eqnarray}
f_{-}^{f_0}(q^{2})=\frac{m^{2}_{B_s}-m^{2}_{f_0}}{q^{2}}(f_{0}^{f_0}(q^{2})-f_{+}^{f_0}(q^{2})).\label{2.2.3}
\end{eqnarray}

The form factors, for the decay $B_{s}\to f_{0}(980)\ell^{+}\ell^{-}$, and $B\to K^{\ast}_0(1430)\ell^{+}\ell^{-}$
can be calculated using the light cone QCD sum rule approach \cite{Colangelo:2010bg}, and three-point QCD sum rules \cite{Aliev:2007rq}.
For whole $q^{2}$ the form factors for the decay can be parameterized in the form of double pole as follows
\begin{eqnarray}
f_{i}^{{f_0}({K^{\ast}_0})}(q^{2})=\frac{f_{i}^{{f_0}
({K^{\ast}_0})}(0)}{1-a_{i}\frac{q^{2}}{m^{2}_{{B_s}(B)}}+b_{i}\frac{q^{4}}{m^{4}_{{B_s}(B)}}}.\label{2.8}
\end{eqnarray}
The numerical values of the form factors at $q^{2}=0$, and parameters $a_{i}$, $b_{i}$, for the decay $B_{s}\to f_{0}(980)\ell^{+}\ell^{-}$, and $B\to K^{\ast}_0(1430)\ell^{+}\ell^{-}$ are presented in table-II of Ref. \cite{Colangelo:2010bg}, and table-I of Ref. \cite{Aliev:2007rq}.
\subsection{$M_{in}\to P\ell^{+}\ell^{-}$ case}
The matrix elements for the process $M_{in}\to P\ell^{+}\ell^{-}$, where both initial $M_{in}=B$, and
final state meson $P=K$, are pseudoscalar in nature, can be expressed as
\begin{eqnarray}
\langle K(k)|\bar s\gamma_{\mu}b|B(p)\rangle&=&f_{+}^K(q^{2})P_{\mu}+f_{-}^K(q^{2})q_{\mu},\label{2.3}\\
\langle K(k)|\bar s\sigma_{\mu\nu}q^{\nu}b|B(p)\rangle&=&-\frac{f_{T}^K(q^{2})}{m_B+m_K}
\Big[q^{2}P_{\mu}-\Big(m^{2}_{B}-m^{2}_{K}\Big)q_{\mu}\Big].\label{2.4}
\end{eqnarray}
Using the above matrix elements in Eq. (\ref{Amp1}), the decay amplitude for $M_{in}\to P\ell^{+}\ell^{-}$, is given as
\begin{eqnarray}
\mathcal{M}\left(M_{in}\to P\ell^{+}\ell^{-}\right)=\frac{G_{F}\alpha}{2\sqrt{2}\pi}V_{tb}V^{\ast}_{ts}\{T^{1,P}_{\mu}(\bar{\ell}\gamma^{\mu}\ell)
+T^{2,P}_{\mu}(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell)\},\label{H1aa}
\end{eqnarray}
where $T^{i,P}_{\mu}$ $(i=1, 2)$, are again given in terms of Wilson coefficients and transition form factors
\begin{eqnarray}
T^{i,K}_{\mu}=P_{\mu}\mathcal{F}^{i,K}_{+}(q^2)+q_{\mu}\mathcal{F}^{i,K}_{-}(q^2),\label{He66}
\end{eqnarray}
\begin{align}
\mathcal{F}^{1,K}_{+}(q^2)&=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
f_{+}^{K}(q^{2})+\frac{2m_{b}}{\big(m_{B}+m_{K}\big)}
C_{7}^{\text{eff}}f_{T}^{K}(q^{2}),\label{He77}\\
\mathcal{F}^{1,K}_{-}(q^2)&=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
f_{-}^{K}(q^{2})-\frac{2m_{b}}{q^{2}}\big(m_{B}-m_{K}\big)
C_{7}^{\text{eff}}f_{T}^{K}(q^{2}),\label{He88}\\
\mathcal{F}^{2,K}_{\pm}(q^2)&=(C_{10}+C_{10\ell}^{\text{NP}}
+C_{10^{\prime}\ell}^{\text{NP}})f_{\pm}^{K}(q^{2}).\label{He22}
\end{align}
For the $B\to K\ell^{+}\ell^{-}$ decay, $f_+^{K}, f_0^{K}, f_T^{K}$ form factors are used in the numerical analysis. Therefore,
$f_{-}^{K}(q^{2})$ form factor is decomposed using a similar expression to Eq. (\ref{2.2.3}). For $B\to K$ transition
form factors, light cone sum rules (LCSR)
predictions can be extrapolated at $q^{2}\leq 8\text{GeV}^{2}$ to the whole
kinematical region by applying  $z$-series expansion \cite{Lu:2018cfc}.
\begin{eqnarray}
z(q^{2},t_{0})=\frac{\sqrt{t_{+}-q^{2}}-\sqrt{t_{+}-t_{0}}}{\sqrt{t_{+}-q^{2}}+\sqrt{t_{+}-t_{0}}},\label{2.5}
\end{eqnarray}
where $t_{+}=(m_{B}+m_{K})^{2}$ and $t_{0}=(m_{B}+m_{K})(\sqrt{m_{B}}-\sqrt{m_{K}})^{2}$. The simplified series
expansion for $B\to P$ form factors has been adopted which was originally proposed in \cite{Bourrely:2008za}.
\begin{eqnarray}
f_{+,T}^K\equiv f^{+,T}_{B\to K}(q^{2})=\frac{f^{+,T}_{B\to K}(0)}{1-q^{2}/m^{2}_{B}}\Big\{1+\sum_{k=1}^{N-1}b^{+,T}_{k,K}\Big(z(q^{2},t_{0})^{k}-z(0,t_{0})^{k}\notag\\
-(-1)^{N-k}\frac{k}{N}[z(q^{2},t_{0})^{N}-z(0,t_{0})^{N}]\Big)\Big\},\label{2.6}\\
f_0^K\equiv f^{0}_{B\to K}(q^{2})=f^{0}_{B\to K}(0)\Big\{1+\sum_{k=1}^{N}b^{0}_{k,K}\big(z(q^{2},t_{0})^{k}-z(0,t_{k})^{k}\big)\Big\}.\label{2.7}
\end{eqnarray}
The numerical values of the form factors at $q^{2}=0$, and parameters given in Eqs. (\ref{2.6}) and (\ref{2.7}) are presented in table-2 of Ref. \cite{Lu:2018cfc}.
\subsection{$M_{in}\to V\ell^{+}\ell^{-}$ case}
For the process $M_{in}\to V\ell^{+}\ell^{-}$, where the parent particle $M_{in}=B$ or $B_s$, and
the daughter particle $V$ is a vector meson $1^{-}$, such as $V=K^{\ast}$ or $\phi$, the matrix elements
for such decays can be parameterized in terms of the form factors as
\begin{align}
\left\langle K^\ast\big(\phi\big)(k,\overline\epsilon)\left\vert \bar{s}\gamma
_{\mu }b\right\vert B\big(B_s\big)(p)\right\rangle &=\frac{2\epsilon_{\mu\nu\alpha\beta}}
{m_{B(B_s)}+m_{K^\ast(\phi)}}\overline\epsilon^{\,\ast\nu}p^{\alpha}k^{\beta}V^{K^\ast(\phi)}(q^{2}),\label{2.13a}
\\
\left\langle K^\ast\big(\phi\big)(k,\overline\epsilon)\left\vert \bar{s}\gamma_{\mu}\gamma_{5}b\right\vert
B(B_s\big)(p)\right\rangle &=i\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)g_{\mu\nu}\overline\epsilon^{\,\ast\nu}A_{1}^{K^\ast(\phi)}(q^{2})
\notag\\
&-iP_{\mu}(\overline\epsilon^{\,\ast}\cdot q)\frac{A_{2}^{K^\ast(\phi)}(q^{2})}{\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)}\notag\\
&-i\frac{2m_{K^\ast(\phi)}}{q^{2}}q_{\mu}(\overline\epsilon^{\,\ast}\cdot q)
\left[A_{3}^{K^\ast(\phi)}(q^{2})-A_{0}^{K^\ast(\phi)}(q^{2})\right],\label{2.13b}
\end{align}
where
\begin{eqnarray}
A_{3}^{K^\ast(\phi)}(q^{2})&=&\frac{m_{B(B_s)}+m_{K^\ast(\phi)}}{2m_{K^\ast(\phi)}}A_{1}^{K^\ast(\phi)}(q^{2})
-\frac{m_{B(B_s)}-m_{K^\ast(\phi)}}{2m_{K^\ast(\phi)}}A_{2}^{K^\ast(\phi)}(q^{2}),\label{A3}
\end{eqnarray}
with $A_3(0)=A_0(0)$. Here and throughout the whole study, we use $\epsilon_{0123}=+1$ convention for the Levi-Civita tensor.
The additional form factors are the tensor form factors which can be expressed as
\begin{align}
\left\langle K^\ast\big(\phi\big)(k,\overline\epsilon)\left\vert \bar{s}i\sigma
_{\mu \nu }q^{\nu }b\right\vert B\big(B_s\big)(p)\right\rangle
&=-2\epsilon _{\mu\nu\alpha\beta}\overline\epsilon^{\,\ast\nu}p^{\alpha}k^{\beta}T_{1}^{K^\ast(\phi)}(q^{2}),\label{FF11}\\
\left\langle K^\ast\big(\phi\big)(k,\overline\epsilon )\left\vert \bar{s}i\sigma
_{\mu \nu }q^{\nu}\gamma_{5}b\right\vert B\big(B_s\big)(p)\right\rangle
&=i\Big[\left(m^2_{B(B_s)}-m^2_{K^\ast(\phi)}\right)g_{\mu\nu}\overline\epsilon^{\,\ast\nu}\notag\\
&-(\overline\epsilon^{\,\ast }\cdot q)P_{\mu}\Big]T_{2}^{K^\ast(\phi)}(q^{2})+i(\overline\epsilon^{\,\ast}\cdot q)\notag\\
&\times\left[q_{\mu}-\frac{q^{2}}{m^2_{B(B_s)}-m^2_{K^\ast(\phi)}}P_{\mu}
\right]T_{3}^{K^\ast(\phi)}(q^{2}),\label{F3}
\end{align}
with $T_{1}(0)=T_{2}(0)$. With the help of above matrix elements, for the
case of $M_{in}\to V\ell^{+}\ell^{-}$, generalized decay amplitude in Eq. (\ref{Amp1}), takes the form
\begin{eqnarray}
\mathcal{M}\left(M_{in}\to V\ell^{+}\ell^{-}\right)=\frac{G_{F}\alpha}{2\sqrt{2}\pi}V_{tb}V^{\ast}_{ts}\{T^{1,V}_{\mu}(\overline\epsilon)(\bar{\ell}\gamma^{\mu}\ell)
+T^{2,V}_{\mu}(\overline\epsilon)(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell)\},\label{H1bb1}
\end{eqnarray}
where
\begin{eqnarray}
T^{1,K^\ast(\phi)}_{\mu}(\overline\epsilon)=\overline\epsilon^{\,\ast\nu}(r)T^{1,K^\ast(\phi)}_{\mu\nu},\label{V11}
\end{eqnarray}
\begin{eqnarray}
T^{1,K^\ast(\phi)}_{\mu\nu}=\epsilon_{\mu\nu\alpha\beta}p^{\alpha}k^{\beta}\mathcal{F}_{1}^{K^\ast(\phi)}
-ig_{\mu\nu}\mathcal{F}_{2}^{K^\ast(\phi)}+iP_{\mu}q_{\nu}\mathcal{F}_{3}^{K^\ast(\phi)}
+iq_{\mu}q_{\nu}\mathcal{F}_{4}^{K^\ast(\phi)},\label{V12}
\end{eqnarray}
\begin{eqnarray}
\mathcal{F}_{1}^{K^\ast(\phi)}&=&2(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})\frac{V^{K^\ast(\phi)}}
{m_{B(B_s)}+m_{K^\ast(\phi)}}
+\frac{4m_{b}}{q^{2}}C_{7}^\text{eff}T_{1}^{K^\ast(\phi)},\label{v13}\\
\mathcal{F}_{2}^{K^\ast(\phi)}&=&(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)A_{1}^{K^\ast(\phi)}\notag\\
&&+\frac{2m_{b}}{q^{2}}C_{7}^\text{eff}\left(m^2_{B(B_s)}-m^2_{K^\ast(\phi)}\right)T_{2}^{K^\ast(\phi)},\label{v14}
\end{eqnarray}
\begin{eqnarray}
\mathcal{F}_{3}^{K^\ast(\phi)}&=&(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
\frac{A_{2}^{K^\ast(\phi)}}{m_{B(B_s)}+m_{K^\ast(\phi)}}\notag\\
&&+\frac{2m_{b}}{q^{2}}C_{7}^\text{eff}\left(T_{2}^{K^\ast(\phi)}+\frac{q^{2}}{\left(m^2_{B(B_s)}-m^2_{K^\ast(\phi)}\right)}
T_{3}^{K^\ast(\phi)}\right),\label{v15}\\
\mathcal{F}_{4}^{K^\ast(\phi)}&=&(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
\frac{2m_{K^\ast(\phi)}}{q^{2}}[A_{3}^{K^\ast(\phi)}-A_{0}^{K^\ast(\phi)}]
-\frac{2m_{b}}{q^{2}}C_{7}^\text{eff}T_{3}^{K^\ast(\phi)}.\label{v16}
\end{eqnarray}
The hadronic amplitude associated with $\bar\ell\gamma^{\mu}\gamma^{5}\ell$ is
\begin{eqnarray}
T^{2,K^\ast(\phi)}_{\mu}(\overline\epsilon)=\overline\epsilon^{\,\ast\nu}(r)T^{2,K^\ast(\phi)}_{\mu\nu},\label{v17}
\end{eqnarray}
where
\begin{eqnarray}
T^{2,K^\ast(\phi)}_{\mu\nu}=\epsilon_{\mu\nu\alpha\beta}p^{\alpha}k^{\beta}\mathcal{F}_{5}^{K^\ast(\phi)}
-ig_{\mu\nu}\mathcal{F}_{6}^{K^\ast(\phi)}+iP_{\mu}q_{\nu}\mathcal{F}_{7}^{K^\ast(\phi)}
+iq_{\mu}q_{\nu}\mathcal{F}_{8}^{K^\ast(\phi)},\label{V18}
\end{eqnarray}
\begin{eqnarray}
\mathcal{F}_{5}^{K^\ast(\phi)}&=&2(C_{10}+C_{10\ell}^{\text{NP}}
+C_{10^{\prime}\ell}^{\text{NP}})\frac{V^{K^\ast(\phi)}}{m_{B(B_s)}+m_{K^\ast(\phi)}},\label{v19}\\
\mathcal{F}_{6}^{K^\ast(\phi)}&=&(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})
\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)A_{1}^{K^\ast(\phi)},\label{v20}\\
\mathcal{F}_{7}^{K^\ast(\phi)}&=&(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})
\frac{A_{2}^{K^\ast(\phi)}}{m_{B(B_s)}+m_{K^\ast(\phi)}},\label{V21}\\
\mathcal{F}_{8}^{K^\ast(\phi)}&=&(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})
\frac{2m_{K^\ast(\phi)}}{q^{2}}[A_{3}^{K^\ast(\phi)}-A_{0}^{K^\ast(\phi)}].\label{v22}
\end{eqnarray}
For the decays $B\to K^{\ast}\ell^{+}\ell^{-}$ and $B_{s}\to\phi\ell^{+}\ell^{-}$,
we use the series expansion fits to LCSR and lattice form factors \cite{Straub:2015ica}. The fits are based on a rapidly converging series parameter
\begin{eqnarray}
z^{\prime}(q^2)=\frac{\sqrt{t^{\prime}_{+}-q^2}-\sqrt{t^{\prime}_{+}-t^{\prime}_{0}}}{\sqrt{t^{\prime}_{+}-q^2}
+\sqrt{t^{\prime}_{+}-t^{\prime}_{0}}},\label{fit}
\end{eqnarray}
where $t^{\prime}_{\pm}=\left(m_{B(B_s)}\pm m_{K^\ast(\phi)}\right)^2$ and $t^{\prime}_{0}=t^{\prime}_{+}\left(1-\sqrt{1-\frac{t^{\prime}_{-}}{t^{\prime}_{+}}}\right)$. The transition
form factors in terms of rapidly converging series parameter can be expressed as \cite{Straub:2015ica}
\begin{eqnarray}
F_{i}(q^{2})=P_{i}(q^{2})\sum_{k}\alpha^{i}_{k}[z^{\prime}(q^{2})-z^{\prime}(0)]^{k},\label{FF}
\end{eqnarray}
where $P_{i}(q^{2})=\frac{1}{(1-q^{2}/m^{2}_{R,i})}$ is simple pole representing the first resonance in the spectrum. In order
to evaluate the transition form factors for $B\to K^{\ast}\ell^{+}\ell^{-}$ and $B_{s}\to\phi\ell^{+}\ell^{-}$ decays, the
form factors are evaluated at $q^{2}=0$ and the numerical values of these form factors are presented in table-4 of Ref. \cite{Straub:2015ica}.
The values of the masses of resonances $m_{R,i}$ and parameters $\alpha^{i}_{k}$ are given in table-3 and table-15, respectively \cite{Straub:2015ica}.
\subsection{$M_{in}\to A\ell^{+}\ell^{-}$ case}
For $M_{in}\to A\ell^{+}\ell^{-}$ decay, where $M_{in}=B$, and $M_{f}=A$, is a final state axial vector
meson $1^{+}$, such as $K_1(1270, 1400)$ meson. For this decay the matrix element can be parameterized in terms
of transition form factors as follows
\begin{eqnarray}
\langle K_1(k,\overline\epsilon)|\bar{s}\gamma_{\mu}b|B(p)\rangle&=&-\left(m_{B}+m_{K_1}\right)g_{\mu\nu}\overline\epsilon^{\,\ast\nu}V_{1}^{K_1}(q^{2})
+P_{\mu}(\overline\epsilon^{\,\ast}\cdot q)\frac{V_{2}^{K_1}(q^{2})}{\left(m_{B}+m_{K_1}\right)}\notag\\
&+&\frac{2m_{K_1}}{q^{2}}q_{\mu}(\overline\epsilon^{\,\ast}\cdot q)[V_{3}^{K_1}(q^{2})-V_{0}^{K_1}(q^{2})],\label{F1}\\
\langle K_1(k,\overline\epsilon)|\bar{s}\gamma_{\mu}\gamma_5 b|B(p)\rangle&=& \frac{2i\epsilon_{\mu\nu\alpha\beta}}{m_{B}+m_{K_1}}
\overline\epsilon^{\,\ast\nu}p^{\alpha}k^{\beta}A^{K_1}(q^{2}),\label{F2}
\end{eqnarray}
where
\begin{eqnarray}
V_{3}^{K_1}(q^{2})&=&\frac{m_{B}+m_{K_1}}{2m_{K_1}}V_{1}^{K_1}(q^{2})-\frac{m_{B}-m_{K_1}}{2m_{K_1}}V_{2}^{K_1}(q^{2}),\label{V3}\\
V_{3}^{K_1}(0)&=&V_{0}^{K_1}(0).\notag
\end{eqnarray}
The other contributions from the tensor form factors are
\begin{eqnarray}
\langle K_1(k,\overline\epsilon)|\bar si\sigma_{\mu\nu}q^{\nu}b|B(p)\rangle&=&\Big[\Big(m_{B}^{2}-m_{K_1}^{2}\Big)g_{\mu\nu}\overline\epsilon^{\,\ast\nu}
-(\overline\epsilon^{\,\ast}\cdot q)P_{\mu}\Big]T_{2}^{K_1}(q^{2})\notag\\
&+&(\overline\epsilon^{\,\ast}\cdot q)\left[q_{\mu}-\frac{q^{2}}{m^2_{B}-m^2_{K_1}}P_{\mu}\right]T_{3}^{K_1}(q^{2}),\label{T1}\\
\langle K_1(k,\overline\epsilon)|\bar si\sigma_{\mu\nu}q^{\nu}\gamma_{5}b|B(p)\rangle&=&2i\epsilon_{\mu\nu\alpha\beta}\overline\epsilon^{\,\ast\nu}p^{\alpha}k^{\beta}T_{1}^{K_1}(q^{2})\,,\label{T2}
\end{eqnarray}
with $T_{1}^{K_1}(0)=T_{2}^{K_1}(0)$. Including the above matrix elements in Eq. (\ref{Amp1}), leads to the decay amplitude
for the case of $M_{in}\to A\ell^{+}\ell^{-}$, to be
\begin{eqnarray}
\mathcal{M}\left(M_{in}\to A\ell^{+}\ell^{-}\right)=\frac{G_{F}\alpha}{2\sqrt{2}\pi}V_{tb}V^{\ast}_{ts}\{T^{1,A}_{\mu}(\overline\epsilon)(\bar{\ell}\gamma^{\mu}\ell)
+T^{2,A}_{\mu}(\overline\epsilon)(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell)\},\label{HA1bb1}
\end{eqnarray}
where
\begin{eqnarray}
T^{1,K_1}_{\mu}(\overline\epsilon)=\overline\epsilon^{\,\ast\nu}(r)T^{1,K_1}_{\mu\nu},\label{VA11}
\end{eqnarray}
\begin{eqnarray}
T^{1,K_1}_{\mu\nu}=-i\epsilon_{\mu\nu\alpha\beta}p^{\alpha}k^{\beta}\mathcal{F}_{1}^{K_1}
-g_{\mu\nu}\mathcal{F}_{2}^{K_1}+P_{\mu}q_{\nu}\mathcal{F}_{3}^{K_1}
+q_{\mu}q_{\nu}\mathcal{F}_{4}^{K_1},\label{VA12}
\end{eqnarray}
\begin{align}
\mathcal{F}_{1}^{K_1}&=2(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})\frac{A^{K_1}}
{m_{B}+m_{K_1}}
+\frac{4m_{b}}{q^{2}}C_{7}^\text{eff}T_{1}^{K_1},\label{vA13}\\
\mathcal{F}_{2}^{K_1}&=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
\left(m_{B}+m_{K_1}\right)V_{1}^{K_1}
+\frac{2m_{b}}{q^{2}}C_{7}^\text{eff}\left(m^2_{B}-m^2_{K_1}\right)T_{2}^{K_1},\label{vA14}\\
\mathcal{F}_{3}^{K_1}&=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
\frac{V_{2}^{K_1}}{m_{B}+m_{K_1}}
+\frac{2m_{b}}{q^{2}}C_{7}^\text{eff}\left(T_{2}^{K_1}+\frac{q^{2}}{(m^2_{B}-m^2_{K_1})}
T_{3}^{K_1}\right),\label{vA15}\\
\mathcal{F}_{4}^{K_1}&=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
\frac{2m_{K_1}}{q^{2}}[V_{3}^{K_1}-V_{0}^{K_1}]
-\frac{2m_{b}}{q^{2}}C_{7}^{\text{eff}}T_{3}^{K_1}.\label{vA16}
\end{align}
Similarly, hadronic amplitude associated with leptonic axial-vector current is
\begin{eqnarray}
T^{2,K_1}_{\mu}(\overline\epsilon)=\overline\epsilon^{\,\ast\nu}(r)T^{2,K_1}_{\mu\nu},\label{vA17}
\end{eqnarray}
where
\begin{eqnarray}
T^{2,K_1}_{\mu\nu}=-i\epsilon_{\mu\nu\alpha\beta}p^{\alpha}k^{\beta}\mathcal{F}_{5}^{K_1}
-g_{\mu\nu}\mathcal{F}_{6}^{K_1}+P_{\mu}q_{\nu}\mathcal{F}_{7}^{K_1}
+q_{\mu}q_{\nu}\mathcal{F}_{8}^{K_1},\label{VA18}
\end{eqnarray}
\begin{eqnarray}
\mathcal{F}_{5}^{K_1}&=&2(C_{10}+C_{10\ell}^{\text{NP}}
-C_{10^{\prime}\ell}^{\text{NP}})\frac{A^{K_1}}{m_{B}+m_{K_1}},\label{vA19}\\
\mathcal{F}_{6}^{K_1}&=&(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})
\left(m_{B}+m_{K_1}\right)V_{1}^{K_1},\label{vA20}\\
\mathcal{F}_{7}^{K_1}&=&(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})
\frac{V_{2}^{K_1}}{m_{B}+m_{K_1}},\label{VA21}\\
\mathcal{F}_{8}^{K_1}&=&(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})
\frac{2m_{K_1}}{q^{2}}[V_{3}^{K_1}-V_{0}^{K_1}].\label{vA22}
\end{eqnarray}
The physical states $K_{1}(1270)$ and $K_{1}(1400)$ are mixed
states of $K_{1A}$ and $K_{1B}$ with mixing angle $\theta_{K_1}$ defined as
\begin{eqnarray}
|K_{1}(1270)\rangle&=&|K_{1A}\rangle\sin\theta_{K_{1}}+|K_{1B}\rangle\cos\theta_{K_{1}}\,,\label{theta}\\
|K_{1}(1400)\rangle&=&|K_{1A}\rangle\cos\theta_{K_{1}}-|K_{1B}\rangle\sin\theta_{K_{1}}\,.\label{theta1}
\end{eqnarray}
The corresponding mixing relations among different matrix elements and for the form factors are explicitly
given in \cite{Hatanaka:2008gu,Ishaq:2013toa,Huang:2018rys,Munir:2015gsp}. For numerical analysis we use the light-cone QCD sum rule form factors \cite{Hatanaka:2008gu}
\begin{equation}
\mathcal{T}^{X}_{i}(q^{2})=\frac{\mathcal{T}^{X}_{i}(0)}{1-a_{i}^{X}\left(q^{2}/m^{2}_{B}\right)+b_{i}^{X}\left(q^{2}/m^{2}_{B}\right)^{2}}\,,\label{m11}
\end{equation}
where $\mathcal{T}$ is $A$, $V$ or $T$ form factors and the subscript
$i$ can take a value 0, 1, 2 or 3, and the superscript $X$ denotes the
$K_{1A}$ or $K_{1B}$ states.
\subsection{$\Lambda_b\to \Lambda\ell^{+}\ell^{-}$ case}
The matrix elements for the process $\Lambda_b\to \Lambda\ell^{+}\ell^{-}$, where both initial $M_{in}=\Lambda_b$, and
final state baryon $\Lambda$, are spin half particles, can be conveniently written in the helicity basis
\begin{align}
\left\langle \Lambda(k, s_{\Lambda})\left\vert \bar{s}\gamma_{\mu}b\right\vert \Lambda_b(p, s_{\Lambda_b})\right\rangle
&=\bar u_{\Lambda}(k, s_{\Lambda})\Big[f_t^V(q^{2})(m_{\Lambda_b}-m_{\Lambda})\frac{q_{\mu}}{q^2}\notag\\
&+f_0^V(q^2)\frac{m_{\Lambda_b}+m_{\Lambda}}{s_+}\big\{p_{\mu}+k_{\mu}-\frac{q_{\mu}}
{q^2}(m^2_{\Lambda_b}-m^2_{\Lambda})\big\}\notag
\\
&+f_{\perp}^V(q^2)\big\{\gamma_{\mu}-\frac{2m_{\Lambda}}{s_+}p_{\mu}-
\frac{2m_{\Lambda_b}}{s_+}k_{\mu}\big\}\Big]u_{\Lambda_b}(p, s_{\Lambda_b}),\label{HM1}
\\
\left\langle \Lambda(k, s_{\Lambda})\left\vert \bar{s}\gamma_{\mu}\gamma_5 b\right\vert \Lambda_b(p, s_{\Lambda_b})\right\rangle
&=-\bar u_{\Lambda}(k, s_{\Lambda})\gamma_5\Big[f_t^A(q^{2})(m_{\Lambda_b}+m_{\Lambda})\frac{q_{\mu}}{q^2}\notag\\
&+f_0^A(q^2)\frac{m_{\Lambda_b}-m_{\Lambda}}{s_-}\big\{p_{\mu}+k_{\mu}-\frac{q_{\mu}}
{q^2}(m^2_{\Lambda_b}-m^2_{\Lambda})\big\}\notag
\\
&+f_{\perp}^A(q^2)\big\{\gamma_{\mu}+\frac{2m_{\Lambda}}{s_-}p_{\mu}-
\frac{2m_{\Lambda_b}}{s_-}k_{\mu}\big\}\Big]u_{\Lambda_b}(p, s_{\Lambda_b}),\label{HM2}
\end{align}
where we have $s\pm=(m_{\Lambda_b}\pm m_\Lambda)^2-q^2$. Additionally,
\begin{align}
\langle \Lambda(k, s_{\Lambda})| \bar{s}i\sigma_{\mu\nu}\gamma^{\nu}&b | \Lambda_b(p, s_{\Lambda_b})\rangle
=-\bar u_{\Lambda}(k, s_{\Lambda})\Big[f_0^T(q^2)\frac{q^2}{s_+}\big\{p_{\mu}+k_{\mu}-\frac{q_{\mu}}
{q^2}(m^2_{\Lambda_b}-m^2_{\Lambda})\big\}\notag
\\
&+f_{\perp}^T(q^2)(m_{\Lambda_b}+m_{\Lambda})\big\{\gamma_{\mu}-\frac{2m_{\Lambda}}{s_+}p_{\mu}-
\frac{2m_{\Lambda_b}}{s_+}k_{\mu}\big\}\Big]u_{\Lambda_b}(p, s_{\Lambda_b}),\label{HM3}
\\
\langle \Lambda(k, s_{\Lambda})| \bar{s}i\sigma_{\mu\nu}\gamma^{\nu}\gamma_5 &b| \Lambda_b(p, s_{\Lambda_b})\rangle
=-\bar u_{\Lambda}(k, s_{\Lambda})\gamma_5\Big[f_0^{T_5}(q^2)\frac{q^2}{s_-}\big\{p_{\mu}+k_{\mu}-\frac{q_{\mu}}
{q^2}(m^2_{\Lambda_b}-m^2_{\Lambda})\big\}\notag
\\
&+f_{\perp}^{T_5}(q^2)(m_{\Lambda_b}-m_{\Lambda})\big\{\gamma_{\mu}+\frac{2m_{\Lambda}}{s_-}p_{\mu}-
\frac{2m_{\Lambda_b}}{s_-}k_{\mu}\big\}\Big]u_{\Lambda_b}(p, s_{\Lambda_b}),\label{HM4}
\end{align}
After including these matrix elements in Eq. (\ref{Amp1}), the decay amplitude for $\Lambda_b\to \Lambda\ell^{+}\ell^{-}$ decay can be given as
\begin{align}
\mathcal{M}^{s_{\Lambda_b},s_{\Lambda}}\left(\Lambda_b\to \Lambda\ell^{+}\ell^{-}\right)=\frac{G_{F}\alpha}{2\sqrt{2}\pi}
V_{tb}V^{\ast}_{ts}&\big\{T^{1}_{\mu}(s_{\Lambda_b}, s_{\Lambda})(\bar{\ell}\gamma^{\mu}\ell)\notag\\
&+T^{2}_{\mu}(s_{\Lambda_b}, s_{\Lambda})(\bar{\ell}\gamma^{\mu}\gamma_{5}\ell)\big\},\label{Lamamp}
\end{align}
where
\begin{align}
T^{1}_{\mu}(s_{\Lambda_b}, s_{\Lambda})&=\bar u_{\Lambda}(k, s_{\Lambda})\Big[q_{\mu}(\mathcal{R}_1+\mathcal{R}_2 \gamma_5)
+p_{\mu}(\mathcal{R}_3+\mathcal{R}_4 \gamma_5)\notag\\
&+k_{\mu}(\mathcal{R}_5+\mathcal{R}_6 \gamma_5)+\gamma_{\mu}
(\mathcal{R}_7-\mathcal{R}_8 \gamma_5)\Big]u_{\Lambda_b}(p, s_{\Lambda_b}),\label{lamamp1}
\end{align}
\begin{align}
\mathcal{R}_{1(2)}=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}&\pm C_{9^{\prime}\ell}^{\text{NP}})\frac{\left(m_{\Lambda_b}\mp m_{\Lambda}\right)}{q^2}
\left(f^{V(A)}_t-\frac{\left(m_{\Lambda_b}\pm m_{\Lambda}\right)^2}{s_{\pm}}f^{V(A)}_0\right)\notag\\
&-\frac{2m_b}{q^2}C_{7}^{\text{eff}}\frac{\left(m^2_{\Lambda_b}-m^2_{\Lambda}\right)}{s_{\pm}}f^{T(T_5)}_0,\label{auxL1}
\end{align}
\begin{align}
\mathcal{R}_{3(4)}=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}&\pm C_{9^{\prime}\ell}^{\text{NP}})
\left(\frac{m_{\Lambda_b}\pm m_{\Lambda}}{s_{\pm}}f^{V(A)}_0\mp\frac{2m_{\Lambda}}{s_{\pm}}f^{V(A)}_{\perp}\right)\notag\\
&+\frac{2m_b}{q^2}C_{7}^{\text{eff}}\left(\frac{q^2}{s_{\pm}}f^{T(T_5)}_0\mp\frac{2m_{\Lambda}}{s_{\pm}}
\left(m_{\Lambda_b}\pm m_{\Lambda}\right)f^{T(T_5)}_{\perp}\right),\label{auxL3}
\\
\mathcal{R}_{5(6)}=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}&\pm C_{9^{\prime}\ell}^{\text{NP}})
\left(\frac{m_{\Lambda_b}\pm m_{\Lambda}}{s_{\pm}}f^{V(A)}_0-\frac{2m_{\Lambda_b}}{s_{\pm}}f^{V(A)}_{\perp}\right)\notag\\
&+\frac{2m_b}{q^2}C_{7}^{\text{eff}}\left(\frac{q^2}{s_{\pm}}f^{T(T_5)}_0-\frac{2m_{\Lambda_b}}{s_{\pm}}
\left(m_{\Lambda_b}\pm m_{\Lambda}\right)f^{T(T_5)}_{\perp}\right),\label{auxL4}
\\
\mathcal{R}_{7(8)}=(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}&\pm C_{9^{\prime}\ell}^{\text{NP}})f^{V(A)}_{\perp}+\frac{2m_b}{q^2}C_{7}^{\text{eff}}
\left(m_{\Lambda_b}\pm m_{\Lambda}\right)f^{T(T_5)}_{\perp},\label{auxL5}
\end{align}
Similarly,
\begin{align}
T^{2}_{\mu}(s_{\Lambda_b}, s_{\Lambda})&=\bar u_{\Lambda}(k, s_{\Lambda})\Big[q_{\mu}(\mathcal{R}_9+\mathcal{R}_{10} \gamma_5)
+p_{\mu}(\mathcal{R}_{11}+\mathcal{R}_{12} \gamma_5)\notag\\
&+k_{\mu}(\mathcal{R}_{13}+\mathcal{R}_{14} \gamma_5)+\gamma_{\mu}
(\mathcal{R}_{15}-\mathcal{R}_{16} \gamma_5)\Big]u_{\Lambda_b}(p, s_{\Lambda_b}),\label{lamamp2}
\end{align}
\begin{align}
&\mathcal{R}_{9(10)}=(C_{10}+C_{10\ell}^{\text{NP}}\pm C_{10^{\prime}\ell}^{\text{NP}})\frac{\left(m_{\Lambda_b}\mp m_{\Lambda}\right)}{q^2}
\left(f^{V(A)}_t-\frac{\left(m_{\Lambda_b}\pm m_{\Lambda}\right)^2}{s_{\pm}}f^{V(A)}_0\right),\label{auxL9}
\\
&\mathcal{R}_{11(12)}=(C_{10}+C_{10\ell}^{\text{NP}}\pm C_{10^{\prime}\ell}^{\text{NP}})
\left(\frac{m_{\Lambda_b}\pm m_{\Lambda}}{s_{\pm}}f^{V(A)}_0\mp \frac{2m_{\Lambda}}{s_{\pm}}f^{V(A)}_{\perp}\right),\label{auxL11}
\\
&\mathcal{R}_{13(14)}=(C_{10}+C_{10\ell}^{\text{NP}}\pm C_{10^{\prime}\ell}^{\text{NP}})
\left(\frac{m_{\Lambda_b}\pm m_{\Lambda}}{s_{\pm}}f^{V(A)}_0-\frac{2m_{\Lambda_b}}{s_{\pm}}f^{V(A)}_{\perp}\right),\label{auxL12}
\\
&\mathcal{R}_{15(16)}=(C_{10}+C_{10\ell}^{\text{NP}}\pm C_{10^{\prime}\ell}^{\text{NP}})f^{V(A)}_{\perp}.\label{auxL13}
\end{align}
For the $\Lambda_b\to \Lambda\ell^{+}\ell^{-}$ decay, we use the lattice QCD results of the form factors for whole $q^2$ range \cite{Detmold:2016pkz}.
The form factors used in \cite{Detmold:2016pkz} are related to our notation of the form factors as $f^V_{t, 0, \perp}=f_{0, +, \perp}$,
$f^A_{t, 0, \perp}=g_{0, +, \perp}$, $f^T_{0, \perp}=h_{+, \perp}$, and $f^{T_5}_{0, \perp}=\tilde h_{+, \perp}$. In the numerical analysis,
to estimate the central value of the LFUV ratio and the total uncertainty due to the form factors, we employ the steps given in Eqs. (50)-(55), in
Ref. \cite{Detmold:2016pkz}. For that we use the ``nominal'' fit and ``higher-order'' fit form factors which correspond to the following $z^{\prime\prime}$-parametrization
\begin{eqnarray}
f(q^2)=\frac{1}{1-q^2/(m^f_{\text{pole}})^2}[a^f_0+a^f_1 z^{\prime\prime}(q^2)],\label{formzzlam}
\end{eqnarray}
\begin{eqnarray}
f(q^2)=\frac{1}{1-q^2/(m^f_{\text{pole}})^2}[a^f_0+a^f_1 z^{\prime\prime}(q^2)+a^f_2 [z^{\prime\prime}(q^2)]^2],\label{formzzlammm}
\end{eqnarray}
where
\begin{eqnarray}
z^{\prime\prime}(q^2)=\frac{\sqrt{t^{\prime\prime}_{+}-q^2}-\sqrt{t^{\prime\prime}_{+}-t^{\prime\prime}_{0}}}
{\sqrt{t^{\prime\prime}_{+}-q^2}+\sqrt{t^{\prime\prime}_{+}-t^{\prime\prime}_{0}}},\label{fitlam}
\end{eqnarray}
with $t^{\prime\prime}_0=(m_{\Lambda_b}-m_{\Lambda})^2$ and $t^{\prime\prime}_+=(m_B+m_K)^2$. The values of the fit parameters
are taken from table-IV and table-V, while the correlation matrices are given in appendix of Ref. \cite{Detmold:2016pkz}.
\section{Helicity framework and helicity amplitudes}\label{secHelicity}
The decay amplitudes can be expressed in terms of helicity basis as described in \cite{Faessler:2002ut}
and references therein. The orthonormality and completeness properties of helicity basis $\varepsilon^{\alpha}(n=t, +, -, 0)$, with three
spin $1$ components orthogonal to momentum transfer i.e., $q\cdot\varepsilon(\pm)=q\cdot\varepsilon=0$, can be expressed as follows
\begin{eqnarray}
\varepsilon^{\ast\alpha}(n)\varepsilon_{\alpha}(l)=g_{nl}, \qquad\quad \sum_{n, l=t, +, -, 0}\varepsilon^{\ast\alpha}(n)\varepsilon^{\beta}(l)g_{nl}=g^{\alpha\beta},\label{C22}
\end{eqnarray}
with $g_{nl}=\text{diag}(+, -, -, -)$. Using the completeness property given in Eq. (\ref{C22}), the contraction of leptonic tensors $L^{(k)\alpha\beta}$ and hadronic tensors $H^{ij}_{\alpha\beta}=T^{i,M_f}_{\alpha}\overline{T}^{\,j,M_f}_{\beta}$ $(i, j=1, 2)$, can be written as
\begin{eqnarray}
L^{(k)\alpha\beta}H^{ij}_{\alpha\beta}=\sum_{n, n^{\prime}, l, l^{\prime}}L^{(k)}_{nl}g_{nn^{\prime}}g_{ll^{\prime}}H^{ij}_{n^{\prime}l^{\prime}},\label{LH}
\end{eqnarray}
where the leptonic and hadronic tensors are expressed in the helicity basis as follows
\begin{eqnarray}
L^{(k)}_{nl}=\varepsilon^{\alpha}(n)\varepsilon^{\ast\beta}(l)L^{(k)}_{\alpha\beta}, &&\qquad H^{ij}_{nl}=\varepsilon^{\ast\alpha}(n)\varepsilon^{\beta}(l)H^{ij}_{\alpha\beta}.\label{LHT}
\end{eqnarray}
Both leptonic and hadronic tensors given in Eq. (\ref{LHT}), will be evaluated in two different frame of references. The lepton
tensor $L^{(k)}_{nl}$ will be evaluated in $l\bar{l}$ CM frame. However the hadron tensor $H^{ij}_{nl}$ will be evaluated
in the rest frame of decaying hadron.
\subsection{Helicity amplitudes for $M_{in}\to S\ell^{+}\ell^{-}$ decays}
\begin{eqnarray}
H^{ij}_{nl}=\big(\varepsilon^{\ast\alpha}(n)T^{i,S}_{\alpha}\big)\cdot\big(\overline{\varepsilon^{\ast\beta}(l)T^{j,S}_{\beta}}\big)
\equiv H^{i,S}_n \,\overline{H}^{\,j,S}_l,\label{HA1}
\end{eqnarray}
where, for $M_{in}=B_s$, and $S=f_0(980)$, explicit helicity amplitudes are obtained as
\begin{eqnarray}
H^{1,f_0}_t&=&i\frac{m^2_{B_s}-m^2_{f_0}}{\sqrt{q^2}}(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})f^{f_0}_0(q^2),\notag
\\
H^{2,f_0}_t&=&i\frac{m^2_{B_s}-m^2_{f_0}}{\sqrt{q^2}}(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})f^{f_0}_0(q^2),\notag
\\
H^{i,f_0}_{\pm}&=&0,\notag
\\
H^{1,f_0}_0&=&i\sqrt{\frac{\lambda}{q^2}}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})f^{f_0}_+(q^2)
+\frac{2m_b}{m_{B_s}+m_{f_0}}C_{7}^{\text{eff}}f^{f_0}_T(q^2)\Big],\notag
\\
H^{2,f_0}_0&=&i\sqrt{\frac{\lambda}{q^2}}\Big[(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})f^{f_0}_+(q^2)\Big],\label{HA2}
\end{eqnarray}
similarly, for $M_{in}=B$, and $S=K^{\ast}_0(1430)$,
\begin{eqnarray}
H^{1,K^{\ast}_0}_t&=&i(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
\Big[\frac{m^2_{B}-m^2_{K^{\ast}_0}}{\sqrt{q^2}}f^{K^{\ast}_0}_+(q^2)+\sqrt{q^2}f^{K^{\ast}_0}_-(q^2)\Big],\notag
\\
H^{2,K^{\ast}_0}_t&=&i(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})
\Big[\frac{m^2_{B}-m^2_{K^{\ast}_0}}{\sqrt{q^2}}f^{K^{\ast}_0}_+(q^2)+\sqrt{q^2}f^{K^{\ast}_0}_-(q^2)\Big],\notag
\\
H^{i,K^{\ast}_0}_{\pm}&=&0,\notag
\\
H^{1,K^{\ast}_0}_0&=&i\sqrt{\frac{\lambda}{q^2}}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})f^{K^{\ast}_0}_+(q^2)
+\frac{2m_b}{m_{B}+m_{K^{\ast}_0}}C_{7}^{\text{eff}}f^{K^{\ast}_0}_T(q^2)\Big],\notag
\\
H^{2,K^{\ast}_0}_0&=&i\sqrt{\frac{\lambda}{q^2}}\Big[(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})f^{K^{\ast}_0}_+(q^2)\Big].\label{HA3}
\end{eqnarray}
Here $\lambda\equiv \lambda(m^2_{B_s(B)}, m^2_{f_0(K^{\ast}_0)}, q^2)$.
\subsection{Helicity amplitudes for $M_{in}\to P\ell^{+}\ell^{-}$ decays}
\begin{eqnarray}
H^{ij}_{nl}=\big(\varepsilon^{\ast\alpha}(n)T^{i,P}_{\alpha}\big)\cdot\big(\overline{\varepsilon^{\ast\beta}(l)T^{j,P}_{\beta}}\big)
\equiv H^{i,P}_n \,\overline{H}^{\,j,P}_l,\label{HAA11}
\end{eqnarray}
where, for $M_{in}=B$, and $P=K$, explicit helicity amplitudes are calculated as
\begin{eqnarray}
H^{1,K}_t&=&\frac{m^2_{B}-m^2_{K}}{\sqrt{q^2}}(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})f^{K}_0(q^2),\notag
\\
H^{2,K}_t&=&\frac{m^2_{B}-m^2_{K}}{\sqrt{q^2}}(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})f^{K}_0(q^2),\notag
\\
H^{i,K}_{\pm}&=&0,\notag
\\
H^{1,K}_0&=&\sqrt{\frac{\lambda}{q^2}}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})f^{K}_+(q^2)
+\frac{2m_b}{m_{B}+m_{K}}C_{7}^{\text{eff}}f^{K}_T(q^2)\Big],\notag
\\
H^{2,K}_0&=&\sqrt{\frac{\lambda}{q^2}}\Big[(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})f^{K}_+(q^2)\Big].\label{HA4}
\end{eqnarray}
Here $\lambda\equiv \lambda(m^2_{B}, m^2_{K}, q^2)$.
\subsection{Helicity amplitudes for $M_{in}\to V\ell^{+}\ell^{-}$ decays}
\begin{eqnarray}
H^{ij}_{nl}&=&\big(\varepsilon^{\ast\alpha}(n)T^{i,V}_{\alpha}\big)\cdot\big(\overline{\varepsilon^{\ast\beta}(l)T^{j,V}_{\beta}}\big)\notag
\\
&=&\big(\varepsilon^{\ast\alpha}(n)\overline\epsilon^{\ast\mu}(r)T^{i,V}_{\alpha,\mu}\big)\cdot\big(\overline{\varepsilon^{\ast\beta}(l)
\overline\epsilon^{\ast\nu}(s)T^{j,V}_{\beta,\nu}}\big)\delta_{rs}\equiv H^{i,V}_n \,\overline{H}^{\, j,V}_l,\label{HA5}
\end{eqnarray}
where, from angular momentum conservation, $r=n$ and $s=l$ for $n, l=\pm,0$ and $r, s=0$ for $n, l=t$.
The explicit helicity amplitudes for $M_{in}=B(B_s)$, and $V=K^{\ast}(\phi)$, are derived in terms of the Wilson coefficients and the form factors as
\begin{align}
H^{1,K^{\ast}(\phi)}_t&=-i\sqrt{\frac{\lambda}{q^2}}(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})A^{K^{\ast}(\phi)}_0,\notag
\\
H^{2,K^{\ast}(\phi)}_t&=-i\sqrt{\frac{\lambda}{q^2}}(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})A^{K^{\ast}(\phi)}_0,\notag
\\
H^{1,K^\ast(\phi)}_{\pm}&=-i\left(m^2_{B(B_s)}-m^2_{K^\ast(\phi)}\right)\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
\frac{A_{1}^{K^\ast(\phi)}}{\left(m_{B(B_s)}-m_{K^\ast(\phi)}\right)}\notag
\\
&+\frac{2m_{b}}{q^{2}}C_{7}^{\text{eff}}T_{2}^{K^\ast(\phi)}\Big]
\pm i\sqrt{\lambda}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
\frac{V^{K^\ast(\phi)}}{\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)}\notag
\\
&\qquad\quad+\frac{2m_{b}}{q^{2}}C_{7}^{\text{eff}}T_{1}^{K^\ast(\phi)}\Big],\notag
\\
H^{2,K^\ast(\phi)}_{\pm}&=-i(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)
A_{1}^{K^\ast(\phi)}\notag
\\
&\pm i\sqrt{\lambda}(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})
\frac{V^{K^\ast(\phi)}}{\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)},\notag
\\
H^{1,K^\ast(\phi)}_0&=-\frac{8im_{B(B_s)}m_{K^\ast(\phi)}}{\sqrt{q^2}}\Bigg[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
A_{12}^{K^\ast(\phi)}+m_b C_{7}^{\text{eff}}\frac{T_{23}^{K^\ast(\phi)}}{m_{B(B_s)}+m_{K^\ast(\phi)}}\Bigg],\notag
\\
H^{2,K^\ast(\phi)}_0&=-\frac{8im_{B(B_s)}m_{K^\ast(\phi)}}{\sqrt{q^2}}\Bigg[(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})
A_{12}^{K^\ast(\phi)}\Bigg],\label{HA6}
\end{align}
where $\lambda\equiv \lambda(m^2_{B(B_s)}, m^2_{K^\ast(\phi)}, q^2)$, and we have used relations between the form factors \cite{Straub:2015ica}
\begin{align}
A_{12}^{K^\ast(\phi)}&=\frac{\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)^2(m^2_{B(B_s)}-m^2_{K^\ast(\phi)}-q^2)
A_{1}^{K^\ast(\phi)}-\lambda A_{2}^{K^\ast(\phi)}}{16m_{B(B_s)}m^2_{K^\ast(\phi)}\left(m_{B(B_s)}+m_{K^\ast(\phi)}\right)},\notag
\\
T_{23}^{K^\ast(\phi)}&=\frac{\left(m^2_{B(B_s)}-m^2_{K^\ast(\phi)}\right)(m^2_{B(B_s)}+3m^2_{K^\ast(\phi)}-q^2)
T_{2}^{K^\ast(\phi)}-\lambda T_{3}^{K^\ast(\phi)}}{8m_{B(B_s)}m^2_{K^\ast(\phi)}\left(m_{B(B_s)}-m_{K^\ast(\phi)}\right)}.\label{relation}
\end{align}
\subsection{Helicity amplitudes for $M_{in}\to A\ell^{+}\ell^{-}$ decays}
\begin{eqnarray}
H^{ij}_{nl}&=&\big(\varepsilon^{\ast\alpha}(n)T^{i,A}_{\alpha}\big)\cdot\big(\overline{\varepsilon^{\ast\beta}(l)T^{j,A}_{\beta}}\big)\notag
\\
&=&\big(\varepsilon^{\ast\alpha}(n)\overline\epsilon^{\ast\mu}(r)T^{i,A}_{\alpha,\mu}\big)\cdot\big(\overline{\varepsilon^{\ast\beta}(l)
\overline\epsilon^{\ast\nu}(s)T^{j,A}_{\beta,\nu}}\big)\delta_{rs}\equiv H^{i,A}_n \,\overline{H}^{\, j,A}_l,\label{HA7}
\end{eqnarray}
where, from angular momentum conservation, $r=n$ and $s=l$ for $n, l=\pm,0$ and $r, s=0$ for $n, l=t$.
The explicit helicity amplitudes for $M_{in}=B$, and $A=K_1$, are given as
\begin{align}
H^{1,K_1}_t&=-\sqrt{\frac{\lambda}{q^2}}(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})V^{K_1}_0,\notag
\\
H^{2,K_1}_t&=-\sqrt{\frac{\lambda}{q^2}}(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})V^{K_1}_0,\notag
\end{align}
\begin{align}
H^{1,K_1}_{\pm}&=-\left(m^2_{B}-m^2_{K_1}\right)\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
\frac{V_{1}^{K_1}}{m_{B}-m_{K_1}}+\frac{2m_{b}}{q^{2}}C_{7}^{\text{eff}}T_{2}^{K_1}\Big]\notag
\\
&\pm \sqrt{\lambda}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
\frac{A^{K_1}}{m_{B}+m_{K_1}}+\frac{2m_{b}}{q^{2}}C_{7}^{\text{eff}}T_{1}^{K_1}\Big],\notag
\\
H^{2,K_1}_{\pm}&=-(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})\left(m_{B}+m_{K_1}\right)
V_{1}^{K_1}\pm \sqrt{\lambda}(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})
\frac{A^{K_1}}{m_{B}+m_{K_1}},\notag
\\
H^{1,K_1}_0&=-\frac{1}{2m_{K_1}\sqrt{q^2}}\Bigg[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
\Big\{(m^2_{B}-m^2_{K_1}-q^2)\left(m_{B}+m_{K_1}\right)V_{1}^{K_1}\notag
\\
&-\frac{\lambda}{m_{B}+m_{K_1}}V_{2}^{K_1}\Big\}+2m_b C_{7}^{\text{eff}}\Big\{(m^2_{B}+3m^2_{K_1}-q^2)T_{2}^{K_1}
-\frac{\lambda}{m^2_{B}-m^2_{K_1}}T_{3}^{K_1}\Big\}
\Bigg],\notag
\\
H^{2,K_1}_0&=-\frac{1}{2m_{K_1}\sqrt{q^2}}(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})
\Bigg[(m^2_{B}-m^2_{K_1}-q^2)\left(m_{B}+m_{K_1}\right)V_{1}^{K_1}\notag
\\
&-\frac{\lambda}{m_{B}+m_{K_1}}V_{2}^{K_1}\Bigg].\label{HA8}
\end{align}
Here $\lambda\equiv \lambda(m^2_{B}, m^2_{K_1}, q^2)$.
\subsection{Helicity amplitudes for $\Lambda_b \to \Lambda\ell^{+}\ell^{-}$ decay}
\begin{eqnarray}
H^{ij}_{nl}&=&\sum_{s_{\Lambda_b}, s_{\Lambda}}\big(\varepsilon^{\ast\alpha}(n)T^{i}_{\alpha}(s_{\Lambda_b}, s_{\Lambda})\big)\cdot\big(\overline{\varepsilon^{\ast\beta}(l)
T^{j}_{\beta}(s_{\Lambda_b}, s_{\Lambda})}\big)\notag
\\
&\equiv&\sum_{s_{\Lambda_b}, s_{\Lambda}}H^{i}_n(s_{\Lambda_b}, s_{\Lambda})\, \overline{H}^{\, j}_l(s_{\Lambda_b}, s_{\Lambda}).\label{HA9}
\end{eqnarray}
The helicity $s_{\Lambda_b}$ of the parent baryon is fixed by angular momentum conservation relation, $s_{\Lambda_b}=-s_{\Lambda}+\lambda_{j_{\text{eff}}}$.
The possible helicity configurations are shown in table \ref{table2}. Using the explicit results of the spinor matrix elements for different combinations of spin orientations, represented in appendix \ref{spinorele}, we work out the expressions of the non-vanishing helicity amplitudes
\begin{align}
H^{1}_t(\pm1/2, \mp1/2)&=\mp\frac{m_{\Lambda_b}-m_{\Lambda}}{\sqrt{q^2}}\sqrt{s_{+}}   (C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})f^V_t\notag
\\
&-\frac{m_{\Lambda_b}+m_{\Lambda}}{\sqrt{q^2}}\sqrt{s_{-}}
(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})f^A_t,\notag
\\
H^{2}_t(\pm1/2, \mp1/2)&=\mp\frac{m_{\Lambda_b}-m_{\Lambda}}{\sqrt{q^2}}\sqrt{s_{+}}
(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})f^V_t\notag
\\
&-\frac{m_{\Lambda_b}+m_{\Lambda}}{\sqrt{q^2}}\sqrt{s_{-}}
(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})f^A_t,\notag
\\
H^{1}_{\pm}(\pm1/2, \pm1/2)&=\pm\sqrt{2s_-}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})f^V_{\perp}
+\frac{2m_b}{q^2} C_{7}^{\text{eff}}(m_{\Lambda_b}+m_{\Lambda})f^T_{\perp}\Big]\notag
\\
&-\sqrt{2s_+}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})f^A_{\perp}
+\frac{2m_b}{q^2} C_{7}^{\text{eff}}(m_{\Lambda_b}-m_{\Lambda})f^{T_5}_{\perp}\Big],\notag
\\
H^{2}_{\pm}(\pm1/2, \pm1/2)&=\pm\sqrt{2s_-}(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})f^V_{\perp}
-\sqrt{2s_+}(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})f^A_{\perp},\notag
\end{align}
\begin{align}
H^{1}_{0}(\pm1/2, \mp1/2)&=\mp\sqrt{\frac{s_-}{q^2}}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}+C_{9^{\prime}\ell}^{\text{NP}})
(m_{\Lambda_b}+m_{\Lambda})f^V_{0}+2m_b C_{7}^{\text{eff}}f^T_{0}\Big]\notag
\\
&-\sqrt{\frac{s_+}{q^2}}\Big[(C_{9}^{\text{eff}}+C_{9\ell}^{\text{NP}}-C_{9^{\prime}\ell}^{\text{NP}})
(m_{\Lambda_b}-m_{\Lambda})f^A_{0}+2m_b C_{7}^{\text{eff}}f^{T_5}_{0}\Big],\notag
\\
H^{2}_{0}(\pm1/2, \mp1/2)&=\mp\sqrt{\frac{s_-}{q^2}}(C_{10}+C_{10\ell}^{\text{NP}}+C_{10^{\prime}\ell}^{\text{NP}})(m_{\Lambda_b}+m_{\Lambda})f^V_{0}\notag
\\
&-\sqrt{\frac{s_+}{q^2}}(C_{10}+C_{10\ell}^{\text{NP}}-C_{10^{\prime}\ell}^{\text{NP}})(m_{\Lambda_b}-m_{\Lambda})f^A_{0}.\label{HA10}
\end{align}
It is important here to mention that the expressions of the helicity amplitudes correspond to intermediate results and depend upon the
kinematics and polarization vectors convention. For the $\Lambda_b \to \Lambda\ell^{+}\ell^{-}$ decay, our conventions are consistent with
that of used in Ref. \cite{Gutsche:2013pp}. However, the final decay observables remain same and are independent of the conventions used.
\begin{table}
\centering
\begin{tabular}{|c|c|c|}
\hline $\hspace{0.5cm}s_{\Lambda_b}\hspace{0.5cm}$ & $\hspace{0.5cm}s_{\Lambda}\hspace{0.5cm}$ & $\hspace{0.5cm}\lambda_{j_{\text{eff}}}\hspace{0.5cm}$ \\
\hline $+\frac{1}{2}$ & $-\frac{1}{2}$ & $0(t)$ \\
$-\frac{1}{2}$ & $+\frac{1}{2}$ & $0(t)$ \\
$+\frac{1}{2}$ & $+\frac{1}{2}$ & $1$ \\
$-\frac{1}{2}$ & $-\frac{1}{2}$ & $-1$ \\
\hline
\end{tabular}
\caption{The possible helicity configurations for $\Lambda_b \to \Lambda\ell^{+}\ell^{-}$ decay.}\label{table2}
\end{table}
\section{Formulation of physical observables}\label{physobs}
The differential decay rate in terms of helicity amplitudes for $M_{in}\to M_{f} \ell^+\ell^-$ transitions, with
$M_{in}=B, B_s$ and $M_{f}=f_0, K_0^{\ast}, K, K^{\ast}, \phi, K_1$, can be expressed as \cite{Faessler:2002ut}
\begin{eqnarray}
\frac{d\Gamma\left(M_{in}\to M_{f} \ell^+\ell^-\right)}{dq^{2}}&=&\frac{G_{F}^2\alpha^2|V_{tb}V^{\ast}_{ts}|^{2}q^2\sqrt{\lambda}\beta_l}{3.2^9 m_{in}^{3}\pi^5}\Bigg[\frac{2m_{\ell}^{2}}{q^{2}}
3\mathcal{R}e\left(H^{2,M_f}_{t}\overline{H}^{\,2,M_f}_{t}\right)\notag
\\
&+&\left(1+\frac{2m_{\ell}^{2}}{q^2}\right)\Big[H^{1,M_f^T}\overline{H}^{\,1,M_f^T}+\mathcal{R}e\left(H^{1,M_f}_{0}\overline{H}^{\,1,M_f}_{0}\right)\Big]\notag
\\
&+&\left(1-\frac{4m_{\ell}^{2}}{q^2}\right)\Big[H^{2,M_f^T}\overline{H}^{\,2,M_f^T}
+\mathcal{R}e\left(H^{2,M_f}_{0}\overline{H}^{\,2,M_f}_{0}\right)\Big]\Bigg],\label{Drateori}
\end{eqnarray}
where
\begin{eqnarray}
H^{i,M_f^T}\overline{H}^{\,i,M_f^T}\equiv \mathcal{R}e\left(H^{i,M_f}_{+}\overline{H}^{\,i,M_f}_{+}\right)
+\mathcal{R}e\left(H^{i,M_f}_{-}\overline{H}^{\,i,M_f}_{-}\right).\label{Drate1}
\end{eqnarray}
When the final state $(M_f)$, is a vector or axial-vector, the longitudinal and transverse
polarizations can be separated and labeled as $L$ and $T$, respectively. The corresponding decay rates are written as
\begin{eqnarray}
\frac{d\Gamma(M_{in}\to M_{f}^L\ell^+\ell^-)}{dq^{2}}&=&\frac{G_{F}^2\alpha^2|V_{tb}V^{\ast}_{ts}|^{2}q^2\sqrt{\lambda}\beta_l}{3.2^9 m_{in}^{3}\pi^5}\Bigg[\frac{2m_{\ell}^{2}}{q^{2}}
3\mathcal{R}e\left(H^{2,M_f}_{t}\overline{H}^{\,2,M_f}_{t}\right)\notag
\\
&+&\left(1+\frac{2m_{\ell}^{2}}{q^2}\right)\mathcal{R}e\left(H^{1,M_f}_{0}\overline{H}^{\, 1,M_f}_{0}\right)\notag
\\
&+&\left(1-\frac{4m_{\ell}^{2}}{q^2}\right)\mathcal{R}e\left(H^{2,M_f}_{0}\overline{H}^{\,2,M_f}_{0}\right)\Bigg],\label{DrateL}
\end{eqnarray}
\begin{eqnarray}
\frac{d\Gamma(M_{in}\to M_{f}^T \ell^+\ell^-)}{dq^{2}}&=&\frac{G_{F}^2\alpha^2|V_{tb}V^{\ast}_{ts}|^{2}q^2\sqrt{\lambda}\beta_l}{3.2^9 m_{in}^{3}\pi^5}\Bigg[\left(1+\frac{2m_{\ell}^{2}}{q^2}\right)H^{1,M_f^T}\overline{H}^{\,1,M_f^T}\notag
\\
&+&\left(1-\frac{4m_{\ell}^{2}}{q^2}\right)H^{2,M_f^T}\overline{H}^{\,2,M_f^T}\Bigg],\label{DrateT}
\end{eqnarray}
Similarly, differential decay rate for $\Lambda_b \to \Lambda\ell^{+}\ell^{-}$ decay is calculated as
\begin{align}
\frac{d\Gamma\left(\Lambda_b \to \Lambda\ell^{+}\ell^{-}\right)}{dq^{2}}&=\frac{G_{F}^2\alpha^2|V_{tb}V^{\ast}_{ts}|^{2}q^2\sqrt{\lambda}\beta_l}{3.2^{10} m_{\Lambda_b}^{3}\pi^5}\Bigg[\frac{2m_{\ell}^{2}}{q^{2}}3\Big\{\left|H^{2}_{t}(+1/2, -1/2)\right|^2\notag
\\
&+\left|H^{2}_{t}(-1/2, +1/2)\right|^2\Big\}+\left(1+\frac{2m_{\ell}^{2}}{q^2}\right)
\Big\{\left|H^{1}_{+}(+1/2, +1/2)\right|^2\notag
\\
&+\left|H^{1}_{-}(-1/2, -1/2)\right|^2+\left|H^{1}_{0}(+1/2, -1/2)\right|^2
+\left|H^{1}_{0}(-1/2, +1/2)\right|^2\Big\}\notag
\\
&+\left(1-\frac{4m_{\ell}^{2}}{q^2}\right)\Big\{\left|H^{2}_{+}(+1/2, +1/2)\right|^2+\left|H^{2}_{-}(-1/2, -1/2)\right|^2\notag
\\
&+\left|H^{2}_{0}(+1/2, -1/2)\right|^2+\left|H^{2}_{0}(-1/2, +1/2)\right|^2\Big\}\Bigg].\label{Drate3}
\end{align}
The decay rate in Eq. (\ref{Drate3}) can be separated into two parts. The first part corresponding to $\Lambda_b$ and $\Lambda$ having opposite spins
is denoted as $d\Gamma\left(\Lambda_b \to \Lambda^0\ell^{+}\ell^{-}\right)/dq^2$, while the other part with $\Lambda_b$ and $\Lambda$ having same spins
is labeled as $d\Gamma\left(\Lambda_b \to \Lambda^1\ell^{+}\ell^{-}\right)/dq^2$.
The LFUV observables are constructed by taking the ratio of decay rates for $M_{in}\to M_{f}\mu^{+}\mu^{-}$
and $M_{in}\to M_{f}e^{+}e^{-}$,
\begin{eqnarray}
R_{M_{f}}\left[q^2_{\text{min}}, q^2_{\text{max}}\right]=\frac{\displaystyle\int^{q^2_{\text{max}}}_{q^2_{\text{min}}}dq^2d\Gamma(M_{in}\to M_{f}\mu^{+}\mu^{-})/dq^{2}}{\displaystyle\int^{q^2_{\text{max}}}_{q^2_{\text{min}}}dq^2d\Gamma(M_{in}\to M_{f}e^{+}e^{-})/dq^{2}}.\label{RK1}
\end{eqnarray}
For the vector and axial-vector final states, polarized LFUV ratios are defined as
\begin{eqnarray}
R_{M_{f}^{L,\,T}}\left[q^2_{\text{min}}, q^2_{\text{max}}\right]=\frac{\displaystyle\int^{q^2_{\text{max}}}_{q^2_{\text{min}}}dq^2d\Gamma(M_{in}\to
M^{L,\,T}_{f}\mu^{+}\mu^{-})/dq^{2}}{\displaystyle\int^{q^2_{\text{max}}}_{q^2_{\text{min}}}dq^2d\Gamma(M_{in}\to M^{L,\,T}_{f}e^{+}e^{-})/dq^{2}}.\label{RK2}
\end{eqnarray}
Similarly, for $\Lambda_b \to \Lambda\ell^{+}\ell^{-}$ decay
\begin{eqnarray}
R_{\Lambda}\left[q^2_{\text{min}}, q^2_{\text{max}}\right]=\frac{\displaystyle\int^{q^2_{\text{max}}}_{q^2_{\text{min}}}dq^2d\Gamma(\Lambda_b \to \Lambda\mu^{+}\mu^{-})/dq^{2}}{\displaystyle\int^{q^2_{\text{max}}}_{q^2_{\text{min}}}dq^2d\Gamma(\Lambda_b \to \Lambda e^{+}e^{-})/dq^{2}},\label{RK3}
\end{eqnarray}
\begin{eqnarray}
R_{\Lambda^{0,\,1}}\left[q^2_{\text{min}}, q^2_{\text{max}}\right]=\frac{\displaystyle\int^{q^2_{\text{max}}}_{q^2_{\text{min}}}dq^2d\Gamma(\Lambda_b \to \Lambda^{0,\,1}\mu^{+}\mu^{-})/dq^{2}}{\displaystyle\int^{q^2_{\text{max}}}_{q^2_{\text{min}}}dq^2d\Gamma(\Lambda_b \to \Lambda^{0,\,1} e^{+}e^{-})/dq^{2}}.\label{RK4}
\end{eqnarray}
\section{Predictions for LFUV ratios in the SM and the NP scenarios}\label{secNum}
\subsection{NP scenarios}
To give predictions and perform numerical analysis of the LFUV ratios, we first specify our choice of the NP scenarios, from the recent global fit analyses,
which could be easily realized in the specific simple NP models.
\begin{enumerate}[label=\arabic*)]
	\item Assuming LFUV NP in $b\to s\mu^+\mu^-$ only, two basic (1D) NP scenarios (S1) $C_{9\mu}^{\text{NP}}$ and (S2) $C_{9\mu}^{\text{NP}}=-C_{10\mu}^{\text{NP}}$, continue to provide better fit to all data, including the latest experimental inputs \cite{Alguero:2019ptt,Datta:2019zca}. Therefore, for the S1 and S2 scenarios, we consider the best-fit values of the Wilson coefficients from table-1
of Ref. \cite{Datta:2019zca}, and collect them in table \ref{tab:bestfitWC}, for the sake of completeness. S1 and S2 can be realized in the simplest
NP models involving the tree-level exchange of a leptoquark (LQ) or a $Z^{\prime}$ boson. While S1 is only possible with a $Z^{\prime}$, S2 can appear
in both LQ and $Z^{\prime}$ models \cite{Alok:2017sui}.
\begin{table*}[!htbp]
	\begin{center}
			\begin{tabular}{|clcc|}
				\hline
				Scenario & & Best-fit value & $1\sigma$ \\
				\hline
                S1 & $C_{9\mu}^{\text{NP}}$                          &$-1.10$       & $[-1.26, -0.94]$\\
                S2 & $C_{9\mu}^{\text{NP}}=-C_{10\mu}^{\text{NP}}$   &$-0.53$       & $[-0.61, -0.45]$ \\
                S3 & $C_{9\mu}^{\text{NP}}=-C_{10\mu}^{\text{NP}}$   &$-0.67$       & $[-0.82, -0.52]$ \\
                   & $C_{9e}^{\text{NP}}=-C_{10e}^{\text{NP}}$       &$-0.28$       & $[-0.48, -0.08]$ \\
                S4 & $C_{9\mu}^{\text{NP}}=-C_{10\mu}^{\text{NP}}$   &$-0.64$       & $[-0.78, -0.50]$ \\
                   & $C_{9e}^{\text{NP}}$                            &$-0.65$       & $[-1.09, -0.21]$ \\
                S5 & $C_{9\mu}^{\text{V}}=-C_{10\mu}^{\text{V}}$     &$-0.30$       & $[-0.42, -0.20]$  \\
                   & $C_{9}^{\text{U}}$                              &$-0.74$       & $[-0.96, -0.51]$  \\
                S6 & $C_{9\mu}^{\text{V}}$                           &$-1.03$       & $[-1.22, -0.84]$ \\
                   & $C_{10^{\prime}}^{\text{U}}$                    &$-0.29$       & $[-0.47, -0.12]$  \\
                \hline
	\end{tabular}		
			\caption{Best-fit values of the Wilson coefficients, and the $1\sigma$ ranges of different NP scenarios with assumptions, such as, purely LFUV NP
        in $b\to s\mu^+\mu^-$, additional arbitrary LFUV NP in $b\to s e^+e^-$ along with LFUV NP in $b\to s\mu^+\mu^-$, and both LFU and
        LFUV NP.}\label{tab:bestfitWC}
	\end{center}
\end{table*}
	\item  Motivated by removing the tensions between the separate fits to LFD and LFUV observables, we consider the NP scenarios extending S1 and S2, with additional arbitrary LFUV NP in $b\to se^+e^-$, which affects only LFUV observables, leading to improved pulls with respect to the SM. While several
scenarios extending S1 and S2, with the addition of one nonzero NP WC in $b\to se^+e^-$ are reported in \cite{Datta:2019zca}, we pick only those scenarios,
which can be realized in the context of the LQ and $Z^{\prime}$ models, and have improved pulls with respect to the SM, compared to the ones obtained in S1 and S2. Therefore, we consider S3 and S4 from table-4 of Ref. \cite{Datta:2019zca}, that can be generated in $Z^{\prime}$ model, whereas only S3 can be realized in the LQ models due to the fact that leptoquarks can only contribute to $C_{9\ell}^{\text{NP}}=-C_{10\ell}^{\text{NP}}$, $\ell=e, \mu$. The best-fit values of the WCs for S3 and S4 are listed
in table \ref{tab:bestfitWC}.
    \item Next, we consider the NP hypothesis which allows LFU NP (equal contributions for all the lepton flavours), in addition to LFUV contributions to
muons only. NP Wilson coefficients in this case can be represented as
\begin{eqnarray}
C_{i^{(\prime)}e}^{\text{NP}}=C_{i^{(\prime)}}^{\text{U}},\qquad\quad C_{i^{(\prime)}\mu}^{\text{NP}}=C_{i^{(\prime)}}^{\text{U}}
+C_{i^{(\prime)}\mu}^{\text{V}},\label{WC910}
\end{eqnarray}
with $i=9, 10$, for the $b\to se^+e^-$ and $b\to s\mu^+\mu^-$ transitions, respectively. The superscript ``U'' and ``V'' represents the
LFU and LFUV contribution, respectively. Several NP scenarios
with both LFU and LFUV NP contributions are presented in table-5 of Ref. \cite{Alguero:2019ptt}. For the sake of simplicity,
we restrict to the NP scenarios, which only extend S1 and S2, yielding equal or improved pulls compared to the corresponding ones for the S1 and S2 scenarios, given in table-1 of \cite{Alguero:2019ptt}, and can be fairly easily realized in specific NP models. It is important to mention that one
should be very careful while comparing pulls found in different analyses, as they strongly depend on the choice of observables, treatment of the
theoretical errors, and the fact that how the analysis is performed. Therefore, we only consider comparison of pulls between scenarios
obtained within a single analysis. Based on the above criteria, we consider S8 and S11 given in table-5 of \cite{Alguero:2019ptt}, and label them
as S5 and S6, as shown in table \ref{tab:bestfitWC}. Scenario S5 can be generated via off-shell
photon penguins \cite{Crivellin:2018yvo} in a LQ model, while S6 can be generated in $Z^{\prime}$ model with vector couplings
to muons and additional Vectorlike quarks with the quantum numbers of left-handed quarks doublets \cite{Bobeth:2016llm}.
\end{enumerate}
It is worth mentioning that the above two methods of considering additional $b\to se^+e^-$ NP are complementary and each NP scenario
in one method can be translated into the other, and vice versa \cite{Kumar:2019qbv}, however, they offer distinct fitting mechanism to
LFD and LFUV observables and therefore may correspond to unique NP predictions. For example, it is suggested in Ref. \cite{Alguero:2018nvb},
that assuming both LFU and LFUV NP provides a different mechanism to obey the constraint from the LFD observable $\mathcal{B}(B_s\to\mu^+\mu^-)$, with large
value of $C_{10\mu}^{\text{V}}$ WC with opposite sign $C_{10}^{\text{U}}$ WC value, and hence allows the possibility of new class of NP models
with large LFU and LFUV contributions to $C_{10\mu}$ at the same time, to account for the combined LFD+LFUV observables. This result is not obtained with only LFUV NP contributions to both $b\to s\mu^+\mu^-$ and $b\to se^+e^-$, as the additional LFUV NP in $b\to se^+e^-$,
affects only LFUV observables and the LFD observables, in this case, explained only by the LFUV NP contributions to $b\to s\mu^+\mu^-$, lead to the other
favoured NP scenarios with large pulls.
\subsection{Predictions for the LFUV ratios}
In this section, we give predictions for the various LFUV ratios, including (pseudo-)scalar final states, $R_{f_0}$, $R_{K^*_0}$, $R_{K}$,
unpolarized and polarized (axial-)vector final states, $R_{\phi^{(L,\,T)}}$, $R_{{K^*}^{(L,\,T)}}$, $R_{K_1^{(L,\,T)}(1270,\,1400)}$, and
for $\Lambda$ baryon with different spin orientations $R_{\Lambda^{(0,\,1)}}$. Experimentally, $R_{K^{(*)}}$ has already been measured by LHCb in the kinematical region $q^2\leq 6$ GeV$^2$, and by Belle in the low and high $q^2$ regions with large errors. Future precision measurements of high $q^2$ bins at Belle II and LHCb will be complementary and important for testing LFU, therefore, in our analysis, we give predictions for only high $q^2$ bins of $R_{K^{(*)}}$.

In figures~\ref{Fig:SPSPlot}-\ref{Fig:RLamda}, we show the predicted values in the SM and the NP scenarios (using the best-fit Wilson coefficients) for the LFUV ratios in the low $q^2$ bin, [0.045, 1] GeV$^2$, the central $q^2$ bin, [1, 6] GeV$^2$, and the high $q^2$ bin, [14, $q^2_{\text{max}}$] GeV$^2$, which include the errors due to the hadronic form factors. We perform our numerical analysis mostly based on these results. Additionally, tables~\ref{tab:Rf0bin}-\ref{tab:RLambdabin}, are given in appendix~\ref{app:LFUV}, where we explicitly list the central values and the errors due to form factors for the LFUV ratios in the SM and the NP scenarios. In order to be more conservative, we also give the deviations of LFUV ratios due to the $1\sigma$ ranges of the best-fit Wilson coefficients in these tables.
\subsubsection{SM and NP predictions for \texorpdfstring{$R_{S}$}{} and \texorpdfstring{$R_{P}$}{}}
\begin{figure*}[ht]
\centering
\includegraphics[scale=0.55]{SPSPlot}
\caption{Predictions for the LFUV ratios involving decays with scalar or pseudoscalar final state particles, $R_{f_0}$, $R_{K^*_0}$, and $R_{K}$. Three kinematical regions, low $[0.045, 1]$ GeV$^2$, central $[1, 6]$ GeV$^2$, and high $[14, q^2_{\text{max}}]$ GeV$^2$, are chosen, where
$q^2_{\text{max}}= 19.2, 14.9$, and $22.9$ GeV$^2$, for $R_{f_0}$, $R_{K^*_0}$, and $R_{K}$, respectively. In each case, predictions from left to right,
correspond to the SM and scenarios S1 to S6, depicted with different colors.}
\label{Fig:SPSPlot}
\end{figure*}
In this section, we consider the LFUV ratios involving decays with scalar or pseudoscalar final states, $R_{f_0}$, $R_{K^*_0}$, and $R_{K}$. SM and NP predictions for these ratios are shown in figure~\ref{Fig:SPSPlot}. It is clear from figure~\ref{Fig:SPSPlot}, that in all $q^2$ bins, NP predictions for
these ratios are considerably lower than the corresponding SM predictions. Considering, $R_{f_0}$ first, the SM predictions of $R_{f_0}$ in the central and high $q^2$ bins are relatively clean and these bins are also very useful to distinguish among the different NP scenarios except between S1 and S6 scenario.
Next, for $R_{K^*_0}$, SM values in the low and central $q^2$ bins are very clean therefore future measurements of $R_{K^*_0}$ in these bins have
the potential to reveal NP unambiguously, however in order to distinguish among the different NP scenarios very precise measurements of $R_{K^*_0}$
will be required in these two bins as the form factor uncertainties in the NP scenarios also largely cancel out. On the other hand, the low $q^2$ bin of $R_{f_0}$ is not very useful for distinguishing the NP scenarios as the NP predictions overlap frequently due to large uncertainties, and the high $q^2$ bin of $R_{K^*_0}$ can only help to differentiate S3 scenario because it does not overlap with the ranges of any other NP scenario. In any case, NP predictions for $R_{f_0}$ and $R_{K^*_0}$ should differ from the SM predictions, therefore it would be very useful for testing LFU by measuring them. In addition, very interestingly, the measurement of $R_K$ at high $q^2$, which can be accessible at Belle II \cite{Kou:2018nap}, can almost distinguish all NP scenarios, making such measurement very anticipated.
\subsubsection{SM and NP predictions for \texorpdfstring{$R_{V^{(L,\,T)}}$}{} and \texorpdfstring{$R_{A^{(L,\,T)}}$}{}}
\begin{figure*}[ht]
\centering
\includegraphics[scale=0.50]{VAPlot}
\caption{Predictions for the LFUV ratios involving decays with vector or axial-vector final state particles, $R_{\phi}$, $R_{K^*}$, $R_{K_1(1270)}$, and $R_{K_1(1400)}$, where only $R_{K_1(1270)}$ values with $\theta_{K_1}=-34^{\circ}$, and $R_{K_1(1400)}$ values with $\theta_{K_1}=34^{\circ}$ are presented.
Three kinematical regions, low $[0.045, 1]$ GeV$^2$, central $[1, 6]$ GeV$^2$, and high $[14, q^2_{\text{max}}]$ GeV$^2$, are chosen, where
$q^2_{\text{max}}= 18.9, 19.2, 16$, and $15$ GeV$^2$, for $R_{\phi}$, $R_{K^*}$, $R_{K_1(1270)}$, and $R_{K_1(1400)}$ respectively. In each case, predictions from left to right, correspond to the SM and scenarios S1 to S6, depicted with different colors.}
\label{Fig:VAPlot}
\end{figure*}
In this section, we consider the LFUV ratios involving decays with unpolarized and polarized vector or axial-vector final states,
$R_{\phi^{(L,\,T)}}$, $R_{{K^*}^{(L,\,T)}}$, $R_{K_1^{(L,\,T)}(1270)}$, and $R_{K_1^{(L,\,T)}(1400)}$. Before presenting our predictions, we need to specify what value of the $K_1$ mixing angle $\theta_{K_1}$, we adopt. In fact, there are two widely used values, i.e., $\theta_{K_1}\sim-34^{\circ}$ \cite{Hatanaka:2008xj}, from $B\to K_1\gamma$ and $\tau\to K_1(1270)\nu$, and $\theta_{K_1}\sim34^{\circ}$ \cite{Cheng:2011pb,Tanabashi:2018oca}, from the study of the $f_1(1285)-f_1(1420)$ and $h_1(1170)-h_1(1380)$ mixing. These different possibilities of $\theta_{K_1}$ lead to different predictions for the observables. In the case of $\theta_{K_1}=-34^{\circ}$, the branching ratio of $\mathcal B(B\to K_1(1400)\ell^+\ell^-)$ is suppressed by one to two orders of magnitude with respect to $\mathcal B(B\to K_1(1270)\ell^+\ell^-)$, while in the case of $\theta_{K_1}=34^{\circ}$, the situation is reversed, $\mathcal B(B\to K_1(1270)\ell^+\ell^-)$ is more suppressed. Given the highly suppressed decay modes are difficult to measure experimentally, we only present the enhanced mode for each possibility of $\theta_{K_1}$, i.e., $B\to K_1(1270)\ell^+\ell^-$, for $\theta=-34^\circ$, and $B\to K_1(1400)\ell^+\ell^-$, for $\theta=34^\circ$. In fact, these two cases have very analogous predictions for $R_{K_1}$ as can be seen in the subsequent analysis.

In figure~\ref{Fig:VAPlot}, we present the SM and the NP predictions for the unpolarized LFUV ratios, $R_{\phi}$, $R_{{K^*}}$, and $R_{K_1(1270,\,1400)}$, where we observe that in the low $q^2$ region [0.045, 1] GeV$^2$, $R_\phi$ is able to discriminate between the SM and the NP values although it cannot distinguish any specific NP scenario, and on the contrary $R_{K_1(1270,\,1400)}$, in the same $q^2$ bin, do not have good sensitivity to NP as the NP predictions overlap with the SM ranges, which also have relatively large uncertainties. With the increase of the momentum transfer, $R_\phi$ and $R_{K_1}$ become more sensitive to NP: in the central $q^2$ region [1, 6] GeV$^2$, $R_\phi$ can distinguish S2 and S5, and $R_{K_1(1270,\,1400)}$ are also able to discriminate S2. The measurement of $R_\phi$ in this region is very useful given the statistical uncertainty can be less than 0.05 after 50 fb$^{-1}$ data is accumulated at LHCb~\cite{Bediaga:2018lhg}. Furthermore, in the high $q^2$ region, sensitivity to NP becomes even more clear as both $R_\phi$ and $R_{K_1(1270,\,1400)}$ have very small errors for the SM and NP predictions, and thus should be able to distinguish among most NP scenarios except between S1 and S6 in the case of $R_\phi$, and between S3 and S6 in case of $R_{K_1(1270,\,1400)}$. Besides, the high $q^2$ bin of $R_{K^*}$ is also very useful for differentiating among the NP scenarios except between S1 and S6, therefore future measurements of both $R_K$ and $R_{K^*}$ in high $q^2$ region would be crucial for probing NP given the Belle II sensitivities are less than $4\%$ with 50 ab$^{-1}$ data~\cite{Kou:2018nap}.
\begin{figure*}[ht]
\centering
\includegraphics[scale=0.50]{VALPlot}
\caption{Same as figure \ref{Fig:VAPlot}, where final state particles are longitudinally polarized, giving polarized LFUV ratios,
$R_{\phi^L}$, $R_{{K^*}^L}$, $R_{K_1^L(1270)}$, and $R_{K_1^L(1400)}$.}
\label{Fig:VALPlot}
\end{figure*}
\begin{figure*}[ht]
\centering
\includegraphics[scale=0.50]{VATPlot}
\caption{Same as figure \ref{Fig:VAPlot}, where final state particles are transversely polarized, giving polarized LFUV ratios,
$R_{\phi^T}$, $R_{{K^*}^T}$, $R_{K_1^T(1270)}$, and $R_{K_1^T(1400)}$.}
\label{Fig:VATPlot}
\end{figure*}

Additionally, in figure~\ref{Fig:VALPlot}, we show the results for the LFUV ratios, with the final vector and axial-vector states longitudinally polarized, $R_{\phi^L}$, $R_{K_1^L(1270,\,1400)}$ and $R_{{K^*}^L}$. These ratios can provide complementary information for testing the
lepton flavour universality. Although the ratios corresponding to longitudinally polarized final states, i.e., $R_{\phi^L}$, $R_{K_1^L(1270,\,1400)}$ and $R_{{K^*}^L}$ have similar behaviours with respect to $R_{\phi}$, $R_{K_1(1270,\,1400)}$ and $R_{K^*}$, $R_{\phi^L}$ and $R_{K_1^L}$ are more sensitive to NP in the central $q^2$ region, especially $R_{\phi^L}$ which can distinguish almost all NP scenarios, if measured precisely. In contrast, the LFUV ratios for transversely polarized final state mesons, as shown in figure~\ref{Fig:VATPlot}, have even more interesting behaviours in the low $q^2$ region: they are sensitive to effects from the NP scenarios except S4 because in these scenarios they are greater than the SM predictions and with small errors. In the central $q^2$ region, $R_{\phi^T}$ in different NP scenarios except S1 and S6 are distinguishable, but in the same $q^2$ bin $R_{K_1^{T}(1270,\,1400)}$ in different NP scenarios have relatively large errors, making it hard to discriminate among the NP scenarios. In the high $q^2$ region, analogous to $R_{\phi^{L}}$, $R_{K_1^{L}(1270,\,1400)}$ and $R_{{K^*}^L}$, the ratios for transverse polarization $R_{\phi^{T}}$, $R_{K_1^{T}(1270,\,1400)}$ and $R_{{K^*}^T}$ in both the SM and the NP scenarios have relatively small errors and most of them can be differentiated from each other.
\subsubsection{SM and NP predictions for \texorpdfstring{$R_{\Lambda^{(0,\,1)}}$}{}}
\begin{figure*}[!t]
\centering
\includegraphics[scale=0.55]{RLamda}
\caption{Predictions for the LFUV ratios, $R_{\Lambda}$, $R_{\Lambda^0}$, and $R_{\Lambda^1}$, involving baryonic final state. Three kinematical regions, low $[0.045, 1]$ GeV$^2$, central $[1, 6]$ GeV$^2$, and high $[14, q^2_{\text{max}}]$ GeV$^2$, are chosen, where
$q^2_{\text{max}}= 20.3$ GeV$^2$. In each case, predictions from left to right, correspond to the SM and scenarios S1 to S6, depicted with different colors.}
\label{Fig:RLamda}
\end{figure*}
In this section, we consider the LFUV ratios $R_{\Lambda^{(0,\,1)}}$. SM and NP predictions for the scenarios S1-S6, for the ratios $R_{\Lambda}$, $R_{\Lambda^0}$, and $R_{\Lambda^1}$ are presented in figure~\ref{Fig:RLamda}. It is observed that the behaviour of $R_\Lambda$ and $R_{\Lambda^{0}}$ is analogous in the sense that the large-recoil bins having low sensitivity to NP cannot distinguish among the different NP scenarios, the central $q^2$ bins with increased sensitivity to NP scenarios can partially distinguish the NP scenarios (S2 for $R_\Lambda$ and S2, S4, and S5 for $R_{\Lambda^{0}}$), and the low-recoil bins with very small errors in the SM and NP predictions can well distinguish all the NP scenarios. In contrast, for $R_{\Lambda^1}$, most NP scenarios are non-distinguishable by using the central $q^2$ bins, partially distinguishable by using the high-recoil bins and almost fully distinguishable by using the low-recoil bins. Therefore, the most remarkable conclusion on $R_{\Lambda^{(0,\,1)}}$ is that it would be most helpful to measure the high $q^2$ bins of $R_{\Lambda^{(0,\,1)}}$ because these bins have very small uncertainties. Lastly, similar to $R_{\phi^T}$, $R_{\Lambda^1}$ corresponding to S1 and S6, in central $q^2$ region may exceed 1, which can be an interesting characteristic for these scenarios, although they are not distinguishable from each other.
\newpage
\section{Summary and conclusions}\label{summary}
In recent years, a number of experimental measurements for the $b\to s\ell^+\ell^-$ transitions have shown deviations from the SM expectations. Such measurements include the branching ratios $\mathcal B(B\to K^{(*)}\mu^+\mu^-)$ and $\mathcal B(B_s\to \phi\mu^+\mu^-)$, the angular observables in $B\to K^*\mu^+\mu^-$ decay including the famous $P_5^{\prime}$ anomaly, and very importantly, the LFUV ratios $R_{K^{(*)}}$ which are ``clean'' probe for LFUV/NP. On the other hand, experimental measurements of the LFUV ratios for either more $b\to s\ell^+\ell^-$ channels or more kinematical regions at Belle II and LHCb have been put on the agenda~\cite{Kou:2018nap,Bediaga:2018lhg}. In light of the current stage, we have studied the LFUV ratios for various $b\to s\ell^+\ell^-$ channels with (pseudo-)scalar and (axial-)vector final state mesons including $R_{f_0}$, $R_{K_0^{\ast}}$, $R_K$, $R_{K^{\ast}}$, $R_\phi$, $R_{K_1}$ as well as  $R_\Lambda$ for $\Lambda_b\to\Lambda\ell^+\ell^-$. In particular, for the cases when spin-1 meson or the $\Lambda$ baryon is the final state, we have also considered the LFUV ratios with the final state hadron longitudinally and transversally polarized.

In our calculation, we have adopted the recent results of hadronic form factors calculated in lattice QCD or/and QCD light-cone sum rules (LCSR). Within the framework of the effective field theory, we have studied various decay channels by employing the helicity formalism, and give the expressions of
the physical observables in terms of the helicity amplitudes by keeping lepton mass effects. Further, we have explicitly worked out
the expressions of the helicity amplitudes in terms of the (NP) Wilson coefficients and the general hadronic form factors, in a consistent
manner, by using the same kinematical configuration and polarization conventions for all the decay channels, which allows others
to easily check and use our expressions.

In the numerical analysis, we have made predictions and performed analysis for the SM and the selected NP scenarios. Given
the updated measurements of $R_{K^{(*)}}$ suggest NP also present in $b\to se^+e^-$, therefore, besides considering
the two basic $b\to s\mu^+\mu^-$ NP scenarios S1 ($C_{9\mu}^{\text{NP}}$ only) or S2 ($C_{9\mu}^{\text{NP}}=-C_{10\mu}^{\text{NP}}$),
we have also considered four NP scenarios which extend S1 and S2, assuming additional arbitrary LFUV NP in $b\to se^+e^-$ or both LFU and LFUV NP,
and have explicit model interpretations. These scenarios are two $b\to s\mu^+\mu^-$ plus $b\to se^+e^-$  NP scenarios S3 and S4, and two LFU plus LFUV NP scenarios S5 and S6. The conclusions on the numerical results can be summarized as follows:
\begin{itemize}
  \item $R_{K^{(*)}}$ in the high $q^2$ region have quite good sensitivity to NP, therefore future precision measurements on the high $q^2$ bins will be important complement to the measurements of low $q^2$ bins in probing LFUV/NP in the $b\to s\ell^+\ell^-$ transition.
  \item Measurements of the LFUV ratios with scalar mesons in final states are also very helpful for distinguishing NP scenarios. In particular, $R_{f_0}$ in the bins [1, 6] GeV$^2$ and $[14,\,q^2_{\text{max}}]$ GeV$^2$, and $R_{K^{*}_0}$ in the bins [0.045, 1] GeV$^2$ and [1, 6] GeV$^2$ are useful because the theoretical uncertainties in these bins are relatively small compared with other bins.
  \item $R_\phi$ is useful for testing LFUV/NP in all kinematical regions and especially in the high $q^2$ region where theoretical predictions have small errors and the different NP scenarios are distinguishable from each other. In contrast, $R_{K_1}$ in the SM and NP scenarios have larger errors
       in the low and central $q^2$ bins and it has good sensitivity to NP in the central and high $q^2$ region, with NP scenarios more distinguishable
       in the high $q^2$ region.
  \item $R_{\phi^L}$ and $R_{K_1^L}$ corresponding to longitudinally polarized final state meson have similar behaviours with respect to the unpolarized ratios $R_\phi$ and $R_{K_1}$ in all $q^2$ bins, but for $R_{\phi^T}$ and $R_{K_1^T}$ corresponding to transversely polarized $\phi$ and $K_1$, the low and high $q^2$ bins are more sensitive to NP, while the central $q^2$ bins in the NP scenarios have large uncertainties.
  \item $R_{\Lambda}$, $R_{\Lambda^0}$ and $R_{\Lambda^1}$ are all very sensitive to NP with tiny theoretical errors in high $q^2$ region, which can be used to distinguish among different NP scenarios. $R_{\Lambda^{(0,\,1)}}$, in the region [0.045, 1] GeV$^2$ do not have good sensitivity to NP, while in [1, 6] GeV$^2$, the measurements of $R_{\Lambda^{(0)}}$ can partially distinguish the NP scenarios, e.g., the central $q^2$ bins of $R_{\Lambda^{(0)}}$ can distinguish S2.
\end{itemize}
In conclusion, similar to $R_{K^{(*)}}$, SM predictions for the LFUV ratios in the various complementary decay channels are
theoretically clean in different kinematical regions and have high sensitivity to NP effects. Therefore, the future measurements on the LFUV ratios for these additional channels, along with the more precise $R_{K^{(*)}}$ measurements, can provide critical information on testing NP/LFUV in the $b\to s\ell^+\ell^-$ FCNC transitions. In addition, LFUV ratios with polarized final state particles are also found to be sensitive to different new physics
scenarios, and therefore can provide additional complementary probe of NP. With the running of Belle II and future upgrade of LHCb,
the measurements of many LFUV ratios studied in this work will be accessible, especially $R_{K^{(*)}}$ in high $q^2$ region
and $R_\phi$ which have already been planned~\cite{Kou:2018nap,Bediaga:2018lhg}. We hope upcoming experimental and theoretical
studies on the LFUV ratios in the $b\to s\ell^+\ell^-$ transitions, along with giving crucial evidence of possible NP behind
the $b\to s\ell^+\ell^-$ anomalies, will also help to identify the true structure of the underlying NP, by differentiating among
the emerging NP scenarios.
\section*{Acknowledgments}
F.M.B would like to acknowledge the financial support from the provincial postdoctoral talent introduction fund under grant No. 2019YJ-01. The work is partly supported by National Science Foundation of China under the Grants 11775012, 11521505 and 11621131001. M.A.P would like to thank the hospitality provided by IHEP when this project was initiated.
