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