Bound states of the Yukawa potential from hidden supersymmetry

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.