Charmonium production allows us to probe QCD at its interplay between the perturbative and non-perturbative regimes. It is thus interesting to assess the convergence of the perturbative expansion in the strong coupling alpha_s for charmonium production and decay. Before going to NNLO accuracy, we will first outline quarkonium phenomenology at NLO accuracy and discuss the appearance of the negative NLO cross-section results in charmonium hadro-production. We present a new scale-fixing criterion which cures the aforementioned negative NLO cross-sections and discuss the interplay with PDFs which are rather unconstrained at scales close to the mass of the charm quark. In order to improve the renormalisation scale uncertainties for charmonium production or decay, it is essential to go to NNLO accuracy. These require the evaluation of the double-virtual contributions which involve two-loop massive Feynman integrals. We outline cutting-edge techniques and methods to compute these two-loop master integrals in a complete analytical fashion and will show some new equivalence relations among master integrals over cross topologies due to special threshold kinematics. We then discuss the techniques to produce high-precision numerics for all master integrals involved. We proceed in presenting the analytical structure for the form-factors and will assess the convergence of the perturbative series in the case of the exclusive decay of the pseudo-scalar state to di-photon.