  We study the $B\to \phi K_S$ and $B\to \ep K_S$ decays in MSSM by calculating
hadronic matrix elements of operators with QCD factorization approach and
including neutral Higgs boson (NHB) contributions. We calculate the Wilson
coefficients of operators including the new operators which are induced by NHB
penguins at LO using the MIA with double insertions. We calculate the
$\alpha_s$ order hadronic matrix elements of the new operators for $B\to \phi
K_s$ and $B\to \eta^\prime K_s$. It is shown that the recent experimental
results on the time-dependent CP asymmetries in $B\to \phi K_S$ and $B\to \ep
K_S$, $S_{\phi K}$ is negative and $S_{\eta^{\prime} K}$ is positive, which can
not be explained in SM, can be explained in MSSM if there are flavor
non-diagonal squark mass matrix elements of 2nd and 3rd generations whose size
satisfies all relevant constraints from known experiments ($B\to X_S\gamma,
B_s\to \mu^+\mu^-, B\to X_s \mu^+\mu^-, B\to X_s g, \Delta M_s$, etc.). In
particular, we find that one can explain the experimental results with a flavor
non-diagonal mass insertion of any chirality when $\alpha_s$ corrections of
hadronic matrix elements of operators are included, in contrast with the claim
in the literature. At the same time, the branching ratios for the two decays
can also be in agreement with experimental measurements.
