  This study focuses on the decay of the $B_c$ meson to S-wave charmonia. Using
lattice inputs on $B_c\to J/\psi$ form factors, we have obtained the
$B_c\to\eta_c$ form factors using heavy quark spin symmetry (HQSS) relations
between the associated form factors after parametrizing and extracting the
possible symmetry breaking corrections. Using the $q^2$ shapes of these form
factors, we have extracted the branching fractions $\mathcal{B}(B_c^-\to
\eta_c\ell^-\bar{\nu})$ (with $\ell =\tau, \mu (e)$) and the decay rate
distributions and have predicted the Standard model estimate for the observable
$R(\eta_c)=\Gamma(B_c^-\to \eta_c\tau^-\bar{\nu})/\Gamma(B_c^-\to
\eta_c\mu^-\bar{\nu}) = 0.302 \pm 0.010$. In addition, we have extracted the
radial wave functions $\psi_{B_c}^R(0)$, $\psi_{J/\psi}^R(0)$ and
$\psi_{\eta_c}^R(0)$ at small quark-antiquark distances from the available
information on the form factors from lattice and experimental data on radiative
and rare decays of the $J/\psi$ and $\eta_c$ mesons. To do so, we choose the
theory framework of nonrelativistic QCD (NRQCD) effective theory. Using our
results, we have estimated the branching fractions of a few non-leptonic decays
of $B_c$ to $J/\psi$ or $\eta_c$ and other light mesons. We have also updated
the numerical estimates of the cross sections $\sigma(e^+e^- \to J/\psi \eta_c,
\eta_c\gamma)$ and predicted the branching fractions of $Z$ boson decays to
either $J/\psi$ or $\eta_c$ final states or both.
