results [@Li:2015zda] by as much as $\sim$40% ($\sim$20%). See Table \[tab:improvements\] for further comparisons. Our spectroscopy is competitive with those obtained from other relativistic models [@Crater:2010fc; @McNeile:2012qf; @Hilger:2014nma; @Leitao:2016bqq]. Not only are the mass spectra improved, the spread of the mass eigenvalues $\overline{\delta_j M}$ due to the violation of rotational symmetry, is also significantly reduced as mentioned. A related issue is the quenching of the hyperfine splitting found within the fixed coupling results. With the running coupling, this issue is resolved and the hyperfine splittings are consistent with the experimental values, as shown in Fig. \[fig:spectroscopy\] (cf. Fig. \[fig:conv\]). Therefore, the violation of the rotational symmetry is significantly reduced.

Figure \[fig:conv\] shows the trends of the charmonium mass eigenvalues as functions of $N_{\max}^{-1}$ (with $N_{\max}=L_{\max}$). The left panel presents the convergence trends of the ground-state masses ($\eta_c$ and $J/\psi$). The right panel presents the convergence trends for the hyperfine splittings between 1S ($M_{J/\psi}-M_{\eta_c}$) and 2S ($M_{\psi'}-M_{\eta'_c}$) states. Two sets of parameters are used: the fix-parameter results use model parameters from the $N_{\max}=L_{\max}=32$ fit; the refit-parameter calculation refits the model parameters for each $N_{\max}=L_{\max}$. Smooth extrapolations are made using three types of functions: $a+b/N_{\max}+c/N_{\max}^2$ (solid), $a+b \exp(-c N_{\max})$ (dashed), $a+b \exp(-c
