  In this work we calculate the mass spectrum, weak decay constants, two photon
decay widths, and two gluon decay widths of ground (1S), and radially excited
(2S,3S,...) states of pseudoscalar charmoniuum and bottomonium such as
$\eta_{c}$ and $\eta_{b}$, as well as the mass spectrum and leptonic decay
constants of ground state (1S), excited (2S,1D,3S,2D,4S,and 3D) states of
vector charmonium and bottomonium such as $J/\psi$, and $\Upsilon$, using the
formulation of Bethe-Salpeter equation under covariant Instantaneous Ansatz
(CIA). Our results are in good agreement with data (where ever available) and
other models. In this framework, from the beginning, we employ a $4\times 4$
representation for two-body ($q\overline{q}$) BS amplitude for calculating both
the mass spectra as well as the transition amplitudes. However, the price we
have to pay is to solve a coupled set of equations for both pseudoscalar and
vector quarkonia, which we have explicitly shown get decoupled in the
heavy-quark approximation, leading to mass spectral equation with analytical
solutions for both masses, as well as eigenfunctions for all the above states,
in an approximate harmonic oscillator basis. Further, in the present framework
of BSE, the hadron-quark vertex function (and the full 4D BS wave function) for
pseudoscalar and vector quarkonia used for calculation of various processes,
accommodates all the Dirac structures from their complete set in a natural
manner.
