  In this work, we present a phenomenological study of the complete analytical
solution to the bound eigenstates and eigenvalues of the Yukawa potential
obtained previously using the hidden supersymmetry of the system and a
systematic expansion of the Yukawa potential in terms of $\delta=a_{0}/D$,
where $a_{0}$ is the Bohr radius and $D$ is the screening length. The
eigenvalues, $\epsilon_{nl}(\delta)$, are given in the form of Taylor series in
$\delta$ which can be systematically calculated to the desired order
$\delta^{k}$. Coulomb $l$-degeneracy is broken by the screening effects and,
for a given $n$, $\epsilon_{nl}(\delta)$ is larger for higher values of $l$
which causes the crossing of levels for $n\ge4$. The convergence radius of the
Taylor series can be enlarged up to the critical values using the Pad\'e
approximants technique which allows us to calculate the eigenvalues with high
precision in the whole rage of values of $\delta$ where bound states exist, and
to reach a precise determination of the critical screening lengths,
$\delta_{nl}$. Eigenstates have a form similar to the solutions of the Coulomb
potential, with the associated Laguerre polynomials replaced by new polynomials
of order $\delta^{k}$ with $r$-dependent coefficients which, in turn, are
polynomials in $r$. In general we find sizable deviations from the Coulomb
radial probabilities only for screening lengths close to their critical values.
We use these solutions to find the squared absolute value at the origin of the
wave function for $l=0$, and their derivatives for $l=1$, for the lowest
states, as functions of $\delta$, which enter the phenomenology of dark matter
bound states in dark gauge theories with a light dark mediator.
