The W-exchange process B s → η c (J/ψ)D is studied with the perturbative QCD approach. Three kinds of wave functions for B s meson and two forms of wave functions for charmonium are considered. It is estimated that branching ratios for B s → η c D, η c D, J/ψD, J/ψD decays are the order of 10 −7 , 10 −8 , 10 −8 , 10 −9 , respectively, where the largest uncertainty is from wave functions. There is a possibility for measuring these decay in the near future. PACS numbers: 13.25.Hw 14.40.Nd 1 The first evidence for B s production in e + e − annihilation at the Υ(5S) resonance was found by the CLEO collaboration [1]. Belle has accumulated 121 fb −1 data at the Υ(5S) resonance, including 7.1×10 6 B s B s pairs [2]. It is estimated that some 5.9×10 8 B s mesons in the dataset of 5 ab −1 at the Υ(5S) resonance in a New Snowmass Year (about 10 7 seconds of actual annual running time [3]) will be collected at the high luminosity e + e − asymmetric SuperKEKB [4]. More and more B s decays will have subjected the Standard Model and new physics to a series of increasingly stringent tests, through observables such as branching ratios, CP-violationg asymmetries and kinematic distributions. In the naive spectator model, the general properties of the B s meson parallel those of the B u,d mesons. The close correspondences between B s and B u,d mesons allow for sensitive tests of hadronic models [2]. Hadronic B decays are complicated on account of strong interaction effects, meanwhile, they will have provided a great opportunity to study perturbative and nonperturbative QCD. For nonleptonic two-body B decays, the low-energy effective Hamil-tonian ansatz and the factorization hypothesis are commonly used. In recent years, several attractive methods have been proposed and widely used to evaluate the hadronic matrix elements (where the local operators in the effective Hamiltonian are sandwiched between initial and final hadron states considered) based on an expansion in the QCD coupling constant α s /π and in the power Λ QCD /m Q (where Λ QCD and m Q are the QCD characteristic scale and the mass of heavy quark Q, respectively), such as the QCD factorization [5], perturbative QCD method (pQCD) [6], soft and collinear effective theory (SCET) [7], etc. Using the operator product expansion and renormalization group equation, the low energy effective Hamiltonian for the B s → η c (J/ψ)D decay can be written as [8]: H ef f = G F √ 2 V * cb V us C 1 (¯ b α c α) V −A (¯ u β s β) V −A + C 2 (¯ b α c β) V −A (¯ u β s α) V −A + G F √ 2 V * ub V cs C 1 (¯ b α u α) V −A (¯ c β s β) V −A + C 2 (¯ b α u β) V −A (¯ c β s α) V −A + H.c., (1) where G F is the Fermi coupling constant. The Cabibbo-Kobayashi-Maskawa (CKM) matrix factors V * cb V us (or V * ub V cs) and the Wilson coefficients C 1,2 describe the strengths of the local four-quark operators in the effective Hamiltonian. α and β are SU(3) color indices. (¯ qq ′) V −A = ¯ qγ µ (1 − γ 5)q ′. The Wilson coefficients, which incorporate the physics contributions from high scales, have been calcula...