  Given the unique role of gravitational form factors in unveiling the internal
mechanics of hadrons, we examine the gravitational form factors (GFFs) of
quarkonium systems within the Dyson-Schwinger and Bethe-Salpeter equations
framework. A contact interaction model is employed, along with a novel approach
to the dressed amputated meson-meson scattering amplitude that makes contact
with the energy-momentum tensor and so the GFFs. The resulting form factors
fulfill the anticipated symmetry constraints. The corresponding charge and mass
radii are also analyzed, as well as the so-called $D-$term. For pseudoscalar
mesons, we show that the $D-$term is bounded within the $(-1, -1/3)$ range; the
bounds corresponding to the massless and infinitely massive cases,
respectively. Considering the current diversity of opinions in the field, we
anticipate that lattice QCD simulations can provide a refined analysis, as
resolving the gravitational form factors and \textit{D}-term of quarkonium
holds significant physical implications.
