  In this paper, we calculate the branching ratios and CP violating asymmetries
of the five $\bar{B}^0_s \to (\pi^0\eta^{(*)},\eta^{(*)}\eta^{(*)})$ decays, by
employing the perturbative QCD (pQCD) factorization approach and with the
inclusion of all currently known next-to-leading order (NLO) contributions. We
find that (a) the NLO contributions can provide about 100% enhancements to the
LO pQCD predictions for the decay rates of $\bar{B}_s^0 \to \eta\eta^\prime$
and $\eta^\prime \eta^\prime$ decays, but result in small changes to
$Br(\bar{B}_s \to \pi^0 \eta^{(*)})$ and $Br(\bar{B}_s \to \eta\eta)$; (b) the
newly known NLO twist-2 and twist-3 contributions to the relevant form factors
can provide about 10% enhancements to the decay rates of the considered decays;
(c) for $\bar{B}_s \to \pi^0 \eta^{(*)}$ decays, their direct CP-violating
asymmetries $\cala_f^{dir}$ could be enhanced significantly by the inclusion of
the NLO contributions; and (d) the pQCD predictions for $Br(\bar{B}_s \to \eta
\eta^{(*)})$ and $Br(\bar{B}_s \to \eta^\prime \eta^\prime)$ can be as large as
$4\times 10^{-5}$, which may be measurable at LHCb or the forthcoming super-B
experiments.
