  A continuum approach to Quantum Chromodynamics (QCD), based upon
Schwinger-Dyson (SD) and Bethe-Salpeter (BS) equations, is employed to provide
a tightly constrained prediction for the $\gamma^{*} \gamma^{*} \rightarrow \{
\pi^0, \eta, \eta', \eta_c, \eta_b \}$ transition form factors (TFFs) and their
corresponding pole contribution to the hadronic light-by-light (HLbL) piece of
the anomalous magnetic moment of the muon ($a_\mu$). This work relies on a
practical and well-tested quark-photon vertex Ansatz approach to evaluate the
TFFs for arbitrary space-like photon virtualities, in the impulse
approximation. The numerical results are parametrized meticulously, ensuring a
reliable evaluation of the HLbL contributions to $a_\mu$. We obtain:
$a_{\mu}^{\pi^0-\textrm{pole}} = (6.14 \pm 0.21) \times 10^{-10}$,
$a_{\mu}^{\eta-\textrm{pole}} = (1.47 \pm 0.19) \times 10^{-10}$,
$a_{\mu}^{\eta'-\textrm{pole}} = (1.36 \pm 0.08) \times 10^{-10}$, yielding a
total value of $a_{\mu}^{\pi^0+\eta+\eta'-\textrm{pole}} = (8.97 \pm 0.48)
\times 10^{-10}$, compatible with contemporary determinations. Notably, we find
that $a_{\mu}^{\eta_c+\eta_b-\textrm{pole}} \approx
a_{\mu}^{\eta_c-\textrm{pole}} = (0.09 \pm 0.01) \times 10^{-10}$, which might
not be negligible once the percent precision in the computation of the light
pseudoscalars is reached.
