Exclusive charmless $B_s$ hadronic decays into $\eta'$ and $\eta$
  Using the next-to-leading order QCD-corrected effective Hamiltonian, charmless exclusive nonleptonic decays of the $B_s$ meson into $\eta$ or $\eta'$ are calculated within the generalized factorization approach. Nonfactorizable contributions, which can be parametrized in terms of the effective number of colors $N_c^{\rm eff}$ for $P P$ and $V P$ decay modes, are studied in two different schemes: (i) the one with the "homogeneous" structure in which $(N_c^{\rm eff})_1 \approx (N_c^{\rm eff})_2 \approx ... \approx (N_c^{\rm eff})_{10}$ is assumed, and (ii) the "heterogeneous" one in which the possibility of $N_c^{\rm eff}(V+A) \neq N_c^{\rm eff}(V-A)$ is considered, where $N_c^{\rm eff}(V+A)$ denotes the effective value of colors for the $(V-A)(V+A)$ penguin operators and $N_c^{\rm eff}(V-A)$ for the $(V-A)(V-A)$ ones. For processes depending on the $N_c^{\rm eff}$-stable $a_i$ such as $\bar{B_s} \to (\pi,\rho) (\eta^{'},\eta)$, the predicted branching ratios are not sensitive to the factorization approach we choose. While for the processes depending on the $N_c^{\rm eff}$-sensitive $a_i$ such as $\bar{B_s} \to \omega \eta^{(')}$, there is a wide range for the branching ratios depending on the choice of the $N_c^{\rm eff}$ involved. We have included the QCD anomaly effect in our calculations and found that it is important for $\bar{B_s} \to \eta^{(')} \eta^{(')}$. The effect of the $(c\bar c) \to \eta'$ mechanism is found to be tiny due to a possible CKM-suppression and the suppression in the decay constants except for the $\bar{B_s} \to \phi \eta$ decay within the "naive" factorization approach, where the internal $W$ diagram is CKM-suppressed and the penguin contributions are compensated. 
B. Tseng