Study of $B_s\to K^{(*)}\ell^+ \ell^-$ decays in the PQCD factorization approach with the lattice QCD input

In this paper, we studied the semileptonic decays $B_s \to K^{(*)} \ell^+ \ell^-$ 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 take the available Lattice QCD results at the end point $q^2_{max}$ as additional input to improve the extrapolation of the form factors to the high $q^2$ region. We calculated the branching ratios and other twelve kinds of physical observables. We find the following points: (a) for $ B_s \to K l^+ l^-$ decays, the PQCD and"PQCD+Lattice"predictions for branching ratios (BRs) ${\cal B}(B_s \to K l^+ l^-)$, the ratios of the BRs $R_K^{e\mu,\mu\tau }$ and the longitudinal polarization asymmetry of the leptons $P_L$ agree well within errors; (b) the"PQCD+Lattice"prediction ${\cal B}(B_s \to K^\ast \mu^+ \mu^-)=(2.58^{+0.69}_{-0.60})\times 10^{-8}$ does agree well with the LHCb measured value $ (2.9\pm 1.1 ) \times 10^{-8}$; (c) for the angular observables $P_{1,2,3}$ and $P^\prime_{4,5,6,8}$ , our theoretical predictions for each lepton $l^-$ are consistent within errors; (d) the theoretical predictions of the angular observables $P_3$ and $P'_{6,8}$ are less than $10^{-2}$ in absolute value, but the magnitude of $P_{1,2}$ and $P'_{4,5}$ are larger than $0.2$; and (e) the PQCD and"PQCD+Lattice"predictions of the binned values of all considered observables in the two $q^2$-bins $[0.1-0.98]$GeV$^2$ and $[1.1-6]$GeV$^2$ generally agree with each other and also be consistent with the LCSR results within errors. We believe that above predictions could be measured by future LHCb and Belle-II experiments.

In this paper K * 0 denotes a vector K * 0 (892) meson, which is reconstructed in the K + π − final state experimentally by selecting candidates within 100 MeV/c 2 of the mass [30]. In LHCb experiment, however, no attempt is made to separate the vector K * 0 from the S-wave or other broad contributions which may present in the selected K + π − pair [14]. Fortunately, the S-wave fraction contribution to the B 0 → K * 0 µ + µ − mode has been measured by the LHCb and found to be small [13]. For the B s case, the S-wave contamination of the B 0 s → K * 0 µ + µ − decay is also unknown now and assumed to be small to that of the B 0 → K * 0 µ + µ − decay. Specifically, the S-wave fraction of F S (B 0 → K * 0 µ + µ − ) = (3.4 ± 0.8)% in the K + π − system [20]. Theoretically, the authors of Ref. [31] found the S-wave contribution will modify differential decay widths by about 10% in the process of B 0 → K − π + ℓ + ℓ − .
Analogous to the ratios R K and R K * for B → K ( * ) l + l − decays, we can define the similar ratios of the BRs R eµ s,K and R eµ s,K * for the B s → K ( * ) ℓ + ℓ − decays: Similarly, we can also define the ratios R µτ s,K and R µτ s,K * in the following form: These new ratios for the B s → K ( * ) ℓ + ℓ − decays, together with the ratios R K ( * ) , can help us to examine the b → (s, d)ℓ + ℓ − transitions in great details and find the possible signals of the new physics (NP) beyond SM. Unlike the well studied B → K ( * ) ℓ + ℓ − decays, the semileptonic B s → K ( * ) ℓ + ℓ − decays have not caught much attention partially due to their lower branching ratios and the lack of the relevant experimental measurements. In recent years, these decays have been studied by several authors for example in Refs. [26][27][28][29], and the first measured branching ratio as listed in Eq. (2) was reported last year by LHCb Collaboration [20]. Besides the measurements for the branching ratios, a precise angular reconstructions of the polarized K * in B s → K ( * ) ℓ + ℓ − decays was discussed in Ref. [8]. Recently, the predictions of several angular observables for the B 0 s → K * ℓ + ℓ − decays were provided using the light cone sum rule (LCSR) and the Lattice QCD method in Ref. [29].
By using the perturbative QCD (PQCD) factorization approach [32][33][34], the semileptonic B s → Kℓ + ℓ + decays have been studied by us in a previous paper [26]. We considered the next-toleading order (NLO) contributions known at 2012 and presented our PQCD predictions of the branching ratios: 73 −0.58 ) × 10 −8 , l = (e, µ), In this paper, we will make a systematic study for the semileptonic decays B s → (K, K * )ℓ + ℓ − with l = (e, µ, τ ), and present the theoretical predictions of many new physical observables: (1) For B s → Kℓ + ℓ − decays, besides the branching ratios, we also calculate its forwardbackward asymmetry A F B (q 2 ), the longitudinal lepton polarization asymmetry P L (q 2 ), the ratios R eµ s,K and R µτ s,K .
(2) For B s → K * ℓ + ℓ − decays, we treat them as a four body decay B s → K * (→ Kπ)ℓ + ℓ − described by four kinematic variables: the lepton invariant mass squared q 2 and three angles (θ K * , θ ℓ , φ). We define and calculate the full angular decay distributions, the transverse amplitudes, the partially integrated decay amplitudes over the angles (θ K * , θ ℓ , φ), the forward-backward (FB) asymmetry A F B (q 2 ), the K * polarization fraction R L,T (q 2 ) and the longitudinal lepton polarization asymmetry P L (q 2 ), and the ratios R eµ s,K * and R µτ s,K * . Since we do not know how to calculate the possible S-wave or other broad contributions related with the reconstruction of Kπ pair [13,20], we add a 10% uncertainty to the PQCD predictions of the branching ratios as an additional theoretical error [31], but neglect it in the calculations for other ratios due to the strong cancellation.
(3) We use both the PQCD factorization approach and "PQCD+Lattice" approach to determine the values and their q 2 -dependence of the B s → K ( * ) transition form factors. We use the Bourrely-Caprini-Lellouch (BCL) parametrization method [35,36] to make the extrapolation for all form factors from the low q 2 region to q 2 max . We will calculate the branching ratios and all other physical observables by using the PQCD approach and "PQCD+Lattice" approach respectively, and compare the theoretical predictions obtained based on different models.
The paper is organized as follows: In Sec. II, we give a short review for the kinematics of the B s → K ( * ) ℓ + ℓ − decays including distribution amplitudes of B s and K ( * ) mesons. Sec. III is devoted to the the theoretical framework including Hamiltonian and transition form factors based on k T factorization formalism. In Sec. IV, we list all the observables for both types of decays considered in this paper. Sec. V contains the numerical results of relevant observables and some phenomenological discussions. We conclude and summarize in the last section.

II. KINEMATICS AND THE WAVE FUNCTIONS
We discuss kinematics of these decays in the large-recoil (low q 2 ) region, where the PQCD factorization approach is applicable to the considered semileptonic decays involving K ( * ) as the The typical Feynman diagrams for the semileptonic decaysB 0 s → K ( * ) ℓ + ℓ − in PQCD approach with the flavor-changing neutral current (FCNC) contributions due to the operators O i denoted as black squares.
final state meson. In the rest frame ofB 0 s meson, we define theB 0 s meson momentum p 1 , the K ( * ) momentum p 2 in the light-cone coordinates as Ref. [37] where the mass ratio r = m K /m Bs or m * K /m Bs , and the factor η ± is defined in the following form: where q = p 1 −p 2 is the lepton-pair four-momentum. For the final state K * meson, its longitudinal and transverse polarization vector ǫ L,T can be written as The momenta of the spectator quarks in B s and K ( * ) mesons are parameterized as we make the approximation in the small k ⊥ . For theB 0 s meson wave function, we use the same parameterizations as in Refs. [26,38] Here only the contribution of the Lorentz structure φ Bs (k 1 ) is taken into account, since the contribution of the second Lorentz structureφ Bs is numerically small and has been neglected. We adopted the B s -meson distribution amplitude the same as B-meson in the SU(3) f limit widely used in the PQCD approach In order to analyze the uncertainties of theoretical predictions induced by the inputs, one usually take ω Bs = 0.50 ± 0.05 GeV for B 0 s meson. The normalization factor N Bs depends on the values of the shape parameter ω Bs and the decay constant f Bs and defined through the normalization relation : [26]. For the pseudoscalar K meson , the wave function can be chosen as the same one in Ref. [39]: where m K 0 and p is the chiral mass and the momentum of the meson K. The parameter ζ = 1 or −1 when the momentum fraction of the quark (anti-quark) of the meson is set to be x. The distribution amplitudes (DA's) of the kaon meson can be found easily in Refs. [40][41][42][43]: where t = 2x − 1, f K is the decay constant of kaon meson and ρ K = m K /m 0 K is the mass ratio. The Gegenbauer moments and other parameters are [40][41][42][43]: The Gegenbauer polynomials appeared in Eqs. (14,15) are of the following form [40][41][42][43]: For the light vector meson K * , the longitudinal and transverse polarization components can provide the contribution. Here we adopt the wave functions of the vector K * as in Ref. [43]: where p and m K * are the momentum and the mass of the K * meson, ǫ L and ǫ T correspond to the longitudinal and transverse polarization vectors of the vector meson, respectively. The φ K * and φ T K * in Eqs. (19,20) are the twist-2 DAs [43]: where f K * and f T K * are the longitudinal and transverse components of the decay constants. The Gegenbauer moments in Eqs. (19,20) are the same ones as those in Ref. [43]: The twist-3 DAs φ s,t K * and φ v,a K * in Eqs. (19,20) are defined with the asymptotic form as in Ref. [43]: For the considered b → dℓ + ℓ − transitions, the effective Hamiltonian in the framework of the SM can be written in the following form [29,[44][45][46]: where is a ratio of the CKM elements, C i (µ) and O i (µ) are the Wilson coefficients and the 4-fermion operators at the renormalization scale µ. In SM, a suitable basis of the operators O i (µ) for b → dℓ + ℓ − transition is given by the current-current operators O u,c 1,2 , the QCD penguin operators O 3−6 , the electromagnetic penguin operator O 7 and the chromomagnetic penguin operator O 8 , as well as the semileptonic operators O 9,10 : where T a denotes the generators of the SU(3) C group and m b is the running b quark mass in the MS scheme; F µν and G a µν are the electromagnetic and chromomagnetic tensors, respectively. The labels V ± A refers to the Lorentz structure γ µ (1 ± γ 5 ). The dominant contribution to b → dℓ + ℓ − transitions are given by O 7 and O 9,10 , as well as O u,c 1,2 . The operator O 7 corresponds to the γpenguin diagram, as shown in Fig. 2(a). The operators O 9,10 describe the sum of the contributions from the Z penguin in Fig. 2(a) and the W box diagrams in Fig 2(b). The current-current operators O u,c 1,2 involve a long-distance (LD) contribution, which origins in the real uū,dd and cc intermediate states, namely the (ρ, ω, φ) and J/ψ family in Fig 2(c), coupled to the lepton pair via the virtual photon. This contribution is proportional to C 9 and can be absorbed into an effective Wilson coefficient C eff 9 [47]. Here we neglect the contribution from subleading chromomagnetic penguin, quark-loop and annihilation diagrams because these effects are highly suppressed [29]. Hence the decay amplitude for b → dl + l − loop transition can be decomposed as   where C ef f 7 (µ) and C ef f 9 (µ) are the effective Wilson coefficients, defined as in Refs. [26,50] The analytic expressions for all Wilson coefficients in the NLO approximation can be found easily in Ref. [48]. The numerical values of the NLO Wilson coeffients C i (µ) at three different renormalization scales µ = (m b /2, m b , 3m b /2) are listed in Table I. Note that the Wilson coefficient C 10 is independent of the µ scale and C 9 (µ) is relatively sensitive to the choice of µ.
The term C ′ b→dγ in Eq. (29) is the absorptive part of b → dγ and was given in Ref. [50] C ′ b→dγ (µ) = iα s where η = α s (m W )/α s (µ), x t = m 2 t /m 2 W and Besides the ordinary Wilson coefficient C 9 (µ), the effective Wilson coefficient C eff 9 (q 2 ) in Eq. (30) also contains two additional effective terms Y pert (ŝ) and Y res (q 2 ). The term Y pert (ŝ) describes the short distance contribution from the soft-gluon emission and the one-loop contribution of the four-quark operators O 1 − O 6 . The term Y res (q 2 ) includes the contributions of the virtual resonances described by the Breit-Wigner form prescribed in Refs. [46,49,[51][52][53].
In the above expressions, ω(ŝ) is the soft-gluon correction to the matrix element of operator O 9 and was given in Refs. [46,54] .
The functions g(m q ,ŝ) in Eqs. (33,34) describe the one-loop (qq) contributions to the four-quark operators O 1 − O 6 , and can be written as where x ≡ 4m 2 q /ŝ. The term Y res in Eq. (34) denotes the long-distance resonance contributions from those B s → K ( * ) V → K ( * ) (V → ℓ + ℓ − ) transitions , where V stands for the possible intermediate resonance states decaying to lepton pairs: (1) The charmless light vector mesons V = (ρ, ω, φ). The kinematic region where the light resonances (ρ, ω, φ) contribute is typically not excluded from the experimental analyses because their effects on branching fractions and other physical observables might be substantial [55].
(2) The cc charmonia V cc = ψ(1S, 2S, 3770, 4040, 4160, 4415). The two lowest charmonium states ψ(1S) and ψ(2S) (i.e. J/ψ and ψ ′ ), whose masses are below the open charm threshold (DD), have tiny width and can induce large breaking of quark-hadron duality. Hence, the narrow charmonia resonance regions are routinely rejected in the theoretical and experimental analysis. For the four higher charmonium resonances, however, they are broad and overlapping throughout the high-q 2 regions. One usually make the integration over the full high-q 2 range.
As reported in Ref. [56], a resonance above ψ(2S) compatible with the ψ(4160) has been observed by LHCb in B → Kµ + µ − decay. Consequently, nearly all available contribution about the J P C = 1 −− charmonium resonances above the open charm threshold should be taken into account [57]. In Table II, we list the properties of all considered intermediate resonance states: their mass, width, and branching fractions of the leptonic decay channel V → l + l − [30]. For the case ℓ = τ , only the fraction of J/ψ(2S) → τ + τ − does not vanish , which equals 3.1 × 10 −3 from Ref. [30]. TABLE II. The masses, decay widths and branching fractions of the dilepton decays of the vector charmonium states [30]. The B s → K transition can be induced by the vector current V µ and the tensor currents T µν : where V µ =dγ µ b and T µν =dσ µν b, and q = p 1 − p 2 is the momentum carried off by the lepton pairs and σ µν = i[γ µ , γ ν ]/2. The B s → K transition form factors F + (q 2 ) and F 0 (q 2 ) can be written as a combination of the auxiliary form factors f 1 (q 2 ) and f 2 (q 2 ) in Eq. (37): We also have the relation F + (0) = F 0 (0) in order to smear the pole at q 2 = 0.
Using the well-studied wave functions as given in Sec. II, we calculated the three B s → K form factors f 1 (q 2 ), f 2 (q 2 ) and F T (q 2 ) in the PQCD factorization approach: where C F = 4/3 is a color factor, r 0 = m 0 K /m Bs , r = m k /m Bs , η is defined in Eq. (8), and the function H i (t i ) in the following form The explicit expressions of the hard functions h 1,2 (x 1 , x 2 , b 1 , b 2 ) , the hard scales t 1,2 and the Sudakov factors S ab (t i ) will be given in Appendix A. For the vector meson K * with polarization vector ǫ * , the relevant form factors for B s → K * transitions are V (q 2 ) and A 0,1,2 (q 2 ) of the vector and axial-vector currents, and T 1,2,3 (q 2 ) of the tensor currents. In the PQCD factorization approach, these seven form factors of B s → K * ℓ + ℓ − decays can be calculated and written in the following form: where r = m K * /m Bs , the twist-2 DAs (φ K * , φ T K * ) and other four twist-3 DAs are defined in Eqs. (19,20,24), the functions H 1,2 (t 1,2 ) are the same ones as those defined in Eq. (44) for B s → K transition, but with a replacement of r = m K /m Bs by r = m K * /m Bs .
Within the SM operator basis, the decay amplitude of B s → Kℓ + ℓ − decay can be written in the following form [58]: where p µ 1 denotes the four-momentum of the B s -meson and m ℓ is the lepton mass. Based on the matrix element of the operators in terms of the form factors, we obtain the double differential decay rate forB s → Kℓ + ℓ − with respect to q 2 and θ ℓ with lepton flavor ℓ [59], The angle θ ℓ is defined as the angle between theB s -direction and the ℓ − -direction in the ℓ + ℓ − rest frame. The corresponding angular coefficients a ℓ , b ℓ and c ℓ can be written as [58,59] with the factor N , where β l = 1 − 4m 2 l withm l = m l / q 2 , α em = 1/137 is the fine structure constant, m ℓ means the lepton mass and λ = λ(m 2 Bs , m 2 K , q 2 ) is the Källen function: λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + bc + ca).
Integration over the polar angle θ ℓ leads to the expression for the differential decay rate, We see that the linear dependence on cos θ ℓ is lost after integration over θ ℓ , consequently , the lepton forward-backward asymmetry A FB will also become zero, Another observable of interest that we calculate is the longitudinal polarization asymmetry P L (q 2 ) of the leptons defined as [45]: where h ℓ = +1(−1) implies a right(left)-handed charged lepton ℓ − in the final state. For thē B s → Kℓ + ℓ − decay, the lepton polarization is given by [45]: For a four body decay, B s → K * (→ πK)ℓ + ℓ − , the decay distribution can be completely described in terms of four kinematic variables [7,84,85]: the lepton invariant mass squared (q 2 ) and three angles θ K * , θ ℓ , and φ.The angle θ K * is the angle between direction of flight of K with respect to B s meson in the rest frame of θ K * , θ ℓ is the angle made by ℓ − with respect to the B s meson in the dilepton rest frame and φ is the azimuthal angle between the two planes formed by dilepton and πK. The full angular decay distribution of B s → K * (→ πK)ℓ + ℓ − is given by [8,29,60], where the functions I(q 2 , θ K * , θ ℓ , φ) are of the following form [8]: = I s 1 sin 2 θ K * + I c 1 cos 2 θ K * + (I s 2 sin 2 θ K * + I c 2 cos 2 θ K * ) cos 2θ ℓ + I 3 sin 2 θ K * sin 2 θ ℓ cos 2φ + I 4 sin 2θ K * sin 2θ ℓ cos φ + I 5 sin 2θ K * sin θ ℓ cos φ + I s 6 sin 2 θ K * cos θ ℓ + I 7 sin 2θ K * sin θ ℓ sin φ + I 8 sin 2θ K * sin 2θ ℓ sin φ + I 9 sin 2 θ K * sin 2 θ ℓ sin 2φ .
The angular coefficients I i of the distributions in above equation can be written in terms of the transverse amplitudes [16,61]. For the massless case there are six such complex amplitudes: For the massive case an additional complex amplitude A t is required. In Table III, we show the expressions for those angular coefficients I i (q 2 ) and the corresponding angular factor f i (θ K * , θ ℓ , φ) as those defined in Refs. [16,61].
The seven transversity amplitudes A R,L 0 , A R,L , A R,L ⊥ and A t , in turn, can be expressed in terms of the relevant B s → K * ℓ + ℓ − form factors [6,8]: where the factors N ℓ and N K * are of the following form: Analogous to Ref. [62], one can write down three partially integrated decay distributions, integrating all but one angle at a time: (1) The θ K * distribution : (2) The θ ℓ distribution : (3) The φ distribution : From the full angular distributions as defined in Eq. (64), we set various coefficients apart and combine them into diverse quantities normalized to the differential decay rate and other observables [62].
(1) The differential decay rate: (2) The lepton forward-backward asymmetry: (3) The K * polarization fraction: where Γ L and Γ T represent the longitudinal and transverse K * polarization decay rates, Alternatively, one can define the quantity F K * L which is a measure of the longitudinally polarized K * 's in the whole ensemble of B s → K * ℓℓ decays, which is linked to R L,T (q 2 ) as: where F K * L is a number obtained by integrating F K * L (q 2 ) over the proper phase space.
The above observables are constructed from Eqs. (63,64) by integrating over the angles in various ranges. These observables which have a form factor dependence in the leading order are called form factor dependent (FFD) observables and generally plagued by the large uncertainties of the form factors. To avoid this problem, a lot of works have been done to construct observables which are theoretically clean in low-q 2 region. Such observables are free from this dependence at the leading order and are called form factor independent (FFI) observables. In this paper, we study both kinds of observables.
As a necessary and sufficient condition, such FFI observable must be invariant under the symmetry transformations of the transverse amplitudes A's; we then say that the observable respects the symmetries of the angular distribution. Fortunately, there exists a systematic procedure to construct all such possible observables as discussed in Ref. [61].
We start defining the following complex vectors [63], With these vectors we can construct the products |n i | 2 = n † i n i and n † i n j , We examine the following (clean) FFI observables [61]: The primed observables are also defined in the following form [29]: These primed observables P ′ 4,5,6,8 are clean and a good approximation to P 4,5,6,8 due to the fact that P 1 ≃ 0 in the SM. From the experimental perspective, fitting the primed observables will be simpler and more efficient despite the whole analysis can be performed directly in terms of the observables P 4,5,6,8 .
Since the most observables are written in terms of the ratios, O ℓ i (q 2 ) = N ℓ i (q 2 )/D ℓ i (q 2 ) with N and D being generically a numerator and a denominator, the integrated quantities are then defined as in Ref. [62]: We also check the physical observables R eµ K,K * and R µτ K,K * , as defined in Eqs. (3,4), since the theoretical uncertainties are largely canceled in the ratio of the branching ratios of B s → K ( * ) ℓ + ℓ − decays.
In the region q 2 < 4m 2 µ , where only the e + e − modes are allowed, there is a large enhancement of B s → K * e + e − due to the 1/q 2 scaling of the photon penguin contribution [64]. If we take 4m 2 e (2m 2 µ ) as the lower limit of the integration for B s → K * e + e − ( B s → K * µ + µ − ) decay mode, we will find the values of the ratios: at the scale µ = m b . In order to remove the phase space effects in the ratio R eµ K * and keep consistent with other analysis [65], we here also use the lower cut of 4m 2 µ for both the electron and muon modes in the definition of the ratio R eµ K * [65]:

V. NUMERICAL RESULTS AND DISCUSSIONS
In the numerical calculations we use the following input parameters (here masses and decay constants are in units of GeV) [30]: 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 F 0,+ (q 2 ) and F T (q 2 ) for B s → Kℓ + ℓ − decays, and the form factors V (q 2 ), A 0,1,2 (q 2 ) and T 1,2,3 (q 2 ) for B s → K * ℓ + ℓ − decays. These form factors have been calculated in rather different theories or models, such as the relativistic quark model (RQM) [66], the light cone sum rule (LCSR) [28,67] and the covariant confined quark model (CCQM) [68]. For the heavy B/B s to light meson (such as K, π, η ′ , ρ, K * , etc) transitions, on the other hand, the relevant form factors at the low q 2 region have been evaluated successfully by employing the PQCD factorization approach for example in Refs. [26, 27, 32-34, 40, 41, 69].
Since the PQCD predictions for the considered form factors are reliable only at the low q 2 region, we usually calculate explicitly the values of the relevant form factors at the low q 2 region, say 0 ≤ q 2 ≤ m 2 τ , and then make an extrapolation for all relevant form factors from the low q 2 region to the large q 2 region by using the pole model parametrization [70,71] or other different methods.
In Refs. [72][73][74], 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. [73,74], we used the Bourrely-Caprini-Lellouch (BCL) parametrization method [35,36] instead of the traditional pole model parametrization since the BCL method has better convergence.
In Table IV and V, we list the values of the lattice QCD results for the relevant B s → K * transition form factors at three reference points of q 2 [81,82] used in this paper. The systematic uncertainties are included.
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 .
12 0.56(9) 0.84 (9) (7) in the low q 2 region: 0 ≤ q 2 ≤ m 2 τ . We then make the extrapolation for these form factors to the large q 2 region up to q 2 max by using proper parametrization method.
(2) In the "PQCD+Lattice " approach, we take the lattice QCD results for the form factors at some large q 2 points as input and then make a combined fit to the PQCD and the lattice QCD results at the low and high q 2 region.
In Table VI, as a comparison, we show the centre values of all relevant form factors in this work and other theoretical predictions as given in Refs. [26,43,66,[75][76][77][78][79][80] at the scale q 2 = 0. The PQCD factorization approach is applied in Refs. [43,77] and in Ref. [26] with the inclusion of the NLO corrections. Covariant confined quark model (CCQM) is used in Ref. [68]. Calculations based on Light-cone sum rules (LCSR) in Ref. [67] with hadronic input parameters, TABLE VI. The theoretical predictions for the form factors of the B s → K ( * ) transitions at q 2 = 0 obtained by using rather different theories or models. and in Ref. [75] with the inclusion of the one-loop radiative corrections. In Ref. [66], the authors used the relativistic quark model based on the quasi-potential approach. In Ref. [76], the authors used the quark model and relativistic dispersion approach. In Ref. [78], the light-cone quark model (LCQM) is utilized based on the basis of the soft collinear effective theory. The authors of Ref. [79] employed the LCSR in the framework of the heavy quark effective theory. In Ref. [80], the authors evaluated the transition form factors in the six-quark effective Hamiltonian approach. One can see that there is no significant difference between the theoretical predictions for the B s → K ( * ) transition form factors evaluated at q 2 = 0 in various models or approaches. transitions. Form factors F +,0,T (q 2 ),V (q 2 ),A 0,1,2 (q 2 ) and T 1,2,3 (q 2 ) are fitted by using Eq. (95). In Table VII, we list the PQCD predictions for the form factors F +,0,T (q 2 ),V (q 2 ),A 0,1,2 (q 2 ) and T 1,2,3 (q 2 ) with the corresponding pole and resonance masses, the fitting parametrization constants "α 0 " , "α 1 " and "α 2 " in Eq. form factors as shown in Table VII are the two major errors from the uncertainties of the parameter ω Bs = 0.50 ± 0.05 GeV and f Bs = 0.23 ± 0.02 GeV.  In Table VIII, we list the "PQCD+Lattice" predictions for the form factors F +,0 (q 2 ), V (q 2 ), A 0,1,2 (q 2 ) and T 1,2,3 (q 2 ) by taking into account the lattice QCD results for the form factors at some points of q 2 as listed in Table (IV and V) from Refs. [81,82], in a similar way as what we did in Refs. [72][73][74]. The additional form factors A 12 (q 2 ) and T 23 (q 2 ) can be defined as the linear combinations of A 1 (q 2 ) and A 2 (q 2 ), T 2 (q 2 ) and T 3 (q 2 ), together with kinematic variable λ as given in Eqs. (10) and (11) from Ref. [82]. In Figs. 3 and 4, we show the q 2 -dependence of the form factors F +,0,T (q 2 ), V (q 2 ), A 0,1,2 (q 2 ) and T 1,2,3 (q 2 ) in the PQCD ( the red curves) and PQCD+Lattice ( the blue curves ) approaches for B s → K ( * ) transitions. In Figs. 3 and 4, the error bars of the initial PQCD and relevant lattice QCD results as listed in Table IV and V are blackened, in order to show them clearly.
B. Observables for B s → Kℓ + ℓ − From the differential decay rates as given in Eq. (59), it is conventional to make the integration over the range of 4m 2 ℓ ≤ q 2 ≤ (m Bs − m K ) 2 . In order to be consistent with the choices made by experiment collaborations in their data analysis, however, we have to cut off the regions of dilepton-mass squared around the charmonium resonances J/ψ(1S) and ψ(2S): i.e., 8.0 < q 2 < 11.0 GeV 2 and 12.5 < q 2 < 15.0 GeV 2 for ℓ = (e, µ, τ ) cases. The PQCD predictions for the branching ratios and the longitudinal polarization asymmetry P L of the semileptonic decaysB s → Kℓ + ℓ − at three different renormalization scales are listed in Table IX and X, where the theoretical errors are the combinations of the two major theoretical errors from the uncertainties of ω Bs = 0.50 ± 0.05 GeV and f Bs = 0.24 ± 0.02 GeV. To reduce the large theoretical uncertainties, we also check the physical observables R eµ K and R µτ K , as defined in Eqs. (3,4), i.e., the ratio of the branching ratios of B s → Kℓ + ℓ − decays.
In Fig. 5, we show the q 2 -dependence of the theoretical predictions of the differential branching fraction dB/dq 2 and the longitudinal lepton polarization P L (q 2 ) for B s → Kℓ + ℓ − decays calculated by using the PQCD (the red solid curves ) and " PQCD+Lattice" (the blue dashed curves) approach, with the choice of the scale µ = m b and q 2 max = 23.71 GeV 2 . The shaded bands indicate the theoretical error of our predictions due to the uncertainties of the input parameters. The two vertical grey blocks are the experimental veto regions [1] to remove contributions from B s → J/ψ(1S)(→ ℓ + ℓ − )K ( the left-hand band) and B s → ψ ′ (2S)(→ ℓ + ℓ − )K (the right-hand band) for the q 2 -dependence of dB/dq 2 and P L .
From the numerical results as listed in Table IX and X, one can see the following points: (1) The theoretical predictions from both PQCD and " PQCD+Lattice" approaches have a relatively weak dependence on the choice of the renormalization scale µ. The variations of the central values due to the µ-dependence are about 10% in magnitude and much smaller than the combined errors from the uncertainties of other input parameters.
(2) For B s → Kµ + µ − decay mode, for example, the theoretical prediction B(B s → Kµ + µ − ) in "PQCD+Lattice" approach is smaller than the corresponding decay rate in the PQCD approach: 1.13 +0.28 −0.36 × 10 −8 against 1.36 +0.73 −0.64 × 10 −8 . The difference of their central values is still smaller than the errors. For other two decay modes with final electron or muon leptons, we find similar results.
(3) For the ratio R eµ K , we find R eµ K = 0.997(1) in the PQCD and "PQCD+Lattice" approaches and for µ = (0.5m b , m b , 1.5m b ). For the ratio R µτ K , it has a weak µ-dependence: less than 8% in magnitude for 0.5m b ≤ µ ≤ 1.5m b . The PQCD prediction for R µτ K is a little similar than the one in "PQCD+Lattice " approach: 0.305±0.046 against 0.353±0.045 numerically for µ = m b .
(4) For the case of l − = (e − , µ − ), the PQCD and "PQCD+Lattice " predictions for the values of P L are very similar and close to −1 in value. For the τ − lepton, however, the PQCD and "PQCD+Lattice " prediction is P L ∼ −0.39 and P L ∼ −0.24, respectively. The difference is relatively large in size.   XI. The PQCD and "PQCD+Lattice" predictions for the branching ratios (in unit of 10 −8 ) of the semileptonic decays B s → K * ℓ + ℓ − , the lepton forward-backward asymmetry A F B , the K * polarization fraction F K * L and the ratios of the branching ratios R eµ K * and R µτ K * .
In Table XI, we listed the PQCD and " PQCD+Lattice" predictions for the branching ratios B(B s → K * ℓ + ℓ − ) with l = (e, µ, τ ) , the lepton forward-backward asymmetry, the longitudinal polarization asymmetry of the leptons, and the ratios of the branching ratios R eµ K * and R µτ K * with the choice of the scale µ = m b . In Table XII, we listed the the PQCD and " PQCD+Lattice" predictions for the values of those angular observables P i (i = 1, 2, 3) and P ′ j (j = 4, 5, 6, 8) in e, µ, τ mode. The total error for each physical observable is the combination of the major theoretical errors from the uncertainties of the input parameters ω Bs = 0.50 ± 0.05 GeV and f Bs = 0.24 ± 0.02 GeV. For the branching ratios, the error from the S-wave pollution up to 10% has been added additionally [31]. As a comparison, we also list the LCSR prediction of B(B s → K * µ + µ − ) = (2.897 ± 0.732) × 10 −8 as given the Ref. [29] and the measured value B(B s → K * µ + µ − ) = (2.9 ± 1.1) × 10 −8 as reported by LHCb collaboration [20]. In numerical calculations, we here use the mean value of decay rate Γ B 0 s = 1.509 × 10 12 s −1 [30]. In Fig. 6, we show the PQCD and the "PQCD+Lattice" predictions of q 2 -dependence of the differential decay rate dB/dq 2 , the forward-backward asymmetry A F B (q 2 ), the longitudinal polarization F K * L (q 2 ) for B s → K * ℓ + ℓ − decays with ℓ = (e, µ, τ ) and q 2 max = 20.02 GeV 2 . The red (blue) lines (dashed lines) correspond to the predictions obtained using the PQCD ("PQCD+Lattice") approach, while the shaded narrow bands ( red and blue) indicate the uncertainty of our predictions due to the variations of the input parameters. For the cases of the decays B s → K * ℓ + ℓ − with l = (e, µ), the two vertical grey blocks show the experimental veto regions [1] in order to remove the contributions from the resonance J/ψ(1S) (left-hand band) and ψ ′ (2S) (right-hand band) to the form factor dependent(FFD) observables. For the case of the decay B s → K * τ + τ − , on the other hand, there is one vertical grey block which shows the experimental veto region [1] for the resonance ψ ′ (2S) only.
In Fig. 7 and 8, we show the q 2 -dependence of the angular observables P i (i = 1, 2, 3) and P ′ j (j = 4, 5, 6, 8) for the considered semileptonic decays B s → K * l + l − with ℓ = (e, µ, τ ), respectively. The symbols in these two figures have the same meaning with the ones in Fig. 6. It is easy to see that the PQCD and "PQCD+Lattice" predictions agree very well in the whole range of q 2 .
From the numerical predictions as given in Tables XI and XII and in Figs. 6-8, we find the following points about the physical observables of B s → K * l + l − (ℓ = e, µ, τ ) decays: (1) For the considered decay modes, the PQCD and "PQCD+Lattice" predictions for B(B s → K * l + l − ) with ℓ = (e, µ, τ ) do agree well with each other within the errors. The "PQCD+ Lattice" predictions of B(B s → K * l + l − ) have much smaller errors than those of the PQCD predictions. Both PQCD and "PQCD+Lattice" predictions of B(B s → K * µ + µ − ) do agree well with the only available LHCb measured value [20] within errors. It is also easy to see that the central value of the "PQCD+Lattice" prediction B(B s → K * µ + µ − ) = (2.58 +0.69 −0.60 )× 10 −8 shows a better agreement with the measured one: B(B s → K * µ + µ − ) LHCb = (2.9 ± 1.1) × 10 −8 . For the l − = (e − , τ − ) mode, however, we have to wait for the future experimental measurements.
(2) For the ratio R eµ K * , the theoretical predictions from both PQCD and "PQCD+Latatice" approach are almost the same one, with a tiny ∼ 1% error because of the great cancellation of the errors in the ratio of the branching ratios. For the ratio R µτ K * , however, the remaining error of the theoretical predictions from both PQCD and "PQCD+Latatice" approach are still around 10%. These two ratios should be measured in the future experiments.
(3) For physical observables A F B and F K * L , the differences between the central values of the PQCD and "PQCD+Lattice" are about (20 ∼ 30)% in magnitude, while the error of the theoretical predictions are less than 10%. Bs→K * τ + τ -FIG. 6. The theoretical predictions for the q 2 -dependence of dB/dq 2 , A F B (q 2 ) and F K * L (q 2 ) for the semileptonic decays B s → K * ℓ + ℓ − with ℓ = (e, µ, τ ) in the PQCD and "PQCD+Lattice" approaches. For more details see the text.
(5) One can see from the curves in Figs. 7 and 8 that all angular observables P i and P ′ j have weak q 2 -dependence in the major region of q 2 due to the large cancellation of q 2 -dependence in the ratios.
For the semileptonic decays B → K * ℓ + ℓ − (ℓ = e, µ), some regions of q 2 corresponding to some resonance states, such as the charmonium J/Ψ, ψ(2S), etc, , are usually removed for the sake of date analysis. Following Ref. [29], we here also present the binned value of the observables as a function of lepton-pair momentum q 2 covering two q 2 regions:[0.1-0.98]GeV 2 and [1.1-6]GeV 2 and consider the mass effect in the final state. We employ the PQCD and "PQCD+Lattice" approach to evaluate the form factors and compare the resultant results. The q 2 -binned observables are defined in the following form: , P 2 bin = bin dq 2 (β ℓ I s 6 ) 8 bin dq 2 (I s 2 ) , P 3 bin = − bin dq 2 (I 9 ) 4 bin dq 2 (I s 2 ) , , , In Table XIII, we listed the PQCD and "PQCD+Lattice" predictions for the binned values of all eleven physical observables considered in this paper including the branching ratios. The uncertainties are due to the errors in determination of the form factors . For the decay mode B s → K * µ + µ − , as a comparison, we also insert an extra row of the results from the LCSR approach [29] into the table, for each physical observable. It is necessary to tell the reader that there exist three differences between our predictions and the LCSR results as given in Ref. [29]: (1) The sign definition of the forward-backward asymmetry A F B in Ref. [29] is opposite to ours as given in Eq. (74).  [29]. Because we try to remove the possible contribution from the light resonance φ(1020).
(3) The authors in Ref. [29] considered the nonfactorizable corrections like weak annihilation and spectator scattering in the bin [1 − 6] GeV 2 while these effects in our analysis are very small and have been neglected.
On the theoretical side, from the numerical results as listed in Table XIII, one can find the following points: (1) For observables P 1,2 and A F B (ℓ) , the differences between the PQCD and "PQCD+Lattice predictions are around (20 − 40)% of the central values. For other eight physical observales, however, the PQCD and "PQCD+Lattice predictions agree very well within errors. The source of the difference come from the variations of the form factors of these two factorization approaches.
(2) The differences between our results and LCSR predictions are generally not large in magnitude and could be understood if one takes the three differences between our approaches and the LCSR as specified in last paragraph. Current difference will be tested in the future when the experimental measurements become available.
(3) For observables P 3 and P ′ 6,8 , their SM values are tiny, say about 10 −3 to 10 −2 in magnitude because they are basically driven by the NLO contributions. It is noted that the observable P ′ 6 stems from the absorptive part of b → dγ, a small imaginary one and thus turns to be tiny. Since these observables are not protected from hadronic uncertainties in general, their values are more sensitive to the choice of the method of calculating the form factors or to the variations of the input parameters being used in calculations.
(4) In this paper, the possible long-distance charm loop effects has been taken into account. The modification induced to C 9 was encoded in a shift where the factorizable charm loop and nonfactorizable soft gluon are taken into account. We also use a phenomenological model to account for light resonances like ρ(770) and ω(782) in the low-q 2 region. It is interesting to note that such particular effect is difficult to estimate and can be large in size, casting some doubts on the possibility to exploit the bins between J/ψ and ψ(2S) for comparison with experiments.

VI. SUMMARY AND CONCLUSIONS
In the framework of the SM, we here studied the rare semileptonic decays B s → K ( * ) ℓ + ℓ − with l − = (e − , µ − , τ − ) by using the PQCD and "PQCD+Lattice" factorization approaches and provided the theoretical predictions for the thirteen kinds of physical observables: the branching ratios B(B s → K ( * ) ℓ + ℓ − ), the ratios of the branching ratios R eµ K,K * and R µτ K,K * , the lepton FB asymmetry A F B (l), the longitudinal polarization asymmetry of the leptons P L and the quantity F K * L , the angular observables P i with (i = 1, 2, 3) and P ′ j with (j = 4, 5, 6, 8). In the PQCD factorization approach, specifically, we first evaluated the relevant form factors F 0,+,T (q 2 ), V (q 2 ), A 0,1,2 (q 2 ) and T 1,2,3 (q 2 ) in the low q 2 region and then extrapolate them to the whole q 2 region using the BCL parametrization method. In the "PQCD+Lattice" approach, we also take those currently available Lattice QCD results for the relevant form factors at the end point q 2 max as additional input to improve the extrapolation of the form factors to the whole range of q 2 .
(2) The PQCD and "PQCD+Lattice" predictions for B(B s → K * l + l − ) with ℓ = (e, µ, τ ) do agree well with each other within the errors. Our theoretical predictions of B(B s → K * µ + µ − ) , for example, also agree well with the LHCb measured value [20] as well as the LCSR prediction [29]: (3) For the ratios R eµ K * and R µτ K * , the PQCD and "PQCD+Lattice" predictions agree very well and have a small error less than 10% due to the cancellation of the theoretical uncertainties in the ratios of the branching ratios. For physical observables A F B and F K * L , the differences between the central values of the PQCD and "PQCD+Lattice" are about (20 ∼ 30)% in magnitude, while the error of the theoretical predictions are less than 10%.
(4) For the angular observables P 1,2,3 and P ′ 4,5,6,8 , the PQCD and "PQCD+Lattice " predictions for each lepton l − are consistent within errors. The theoretical predictions of P 3 and P ′ 6,8 are tiny, say less than 10 −2 in absolute value, and thus hard to be measured. For the remaining P 1,2 and P ′ 4,5 , on the other hand, their magnitudes are larger than 0.2 and therefore could be measured by future LHCb and Belle-II experiments. In general, we believe that most above physical observables could be measured in the future LHCb or Belle-II experiments. Any clear deviations from above SM predictions might be a signal of new physics beyond the SM. TABLE XIII. The binned values of observables for the process B s → K * e + e − and B s → K * µ + µ − at µ = m b scale using the PQCD and "PQCD+Lattice" factorization approaches. The uncertainties shown are due to errors in determination of form factors. The LCSR predictions forB s → K * µ + µ − decay as given in Ref. [29] were added as a comparison. For more details, see the text.   (38)