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 δ = a0/D, where a0 is the Bohr radius and D is the screening length. The eigenvalues, nl(δ), are given in the form of Taylor series in δ which can be systematically calculated to the desired order δ. Coulomb l-degeneracy is broken by the screening effects and, for a given n, nl(δ) is larger for higher values of l which causes the crossing of levels for n ≥ 4. The convergence radius of the Taylor series can be enlarged up to the critical values using the Padé approximants technique which allows us to calculate the eigenvalues with high precision in the whole rage of values of δ where bound states exist, and to reach a precise determination of the critical screening lengths, δnl. Eigenstates have a form similar to the solutions of the Coulomb potential, with the associated Laguerre polynomials replaced by new polynomials of order δ 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 δ, which enter the phenomenology of dark matter bound states in dark gauge theories with a light dark mediator.


I. INTRODUCTION
The Yukawa potential is the effective non-relativistic description of the interaction of two particles due to the exchange of a massive particle of mass M . It was proposed in Ref. [1] by H. Yukawa as a low energy description of the strong interactions between nucleons, due to the exchange of massive particles, now known as pions. The potential is given by where α g = g 2 /4π denotes fine-structure constant of the interaction with coupling g. The range of the interaction is given by D = 1/M , also named screening length. For M = 0 we obtain the Coulomb potential and for large M we have effective short-range interactions. The first estimate of the pion mass was done in [1] based on this potential and the first experimental results for the range of the nucleon-nucleon interactions. The Yukawa potential appears in different reserch areas of physics like plasma physics at low density and high temperatures [2][3] [4][5] [6] nuclear physics [7], astrophysics [8] and solid state physics [9] [10] [11]; and is known as Debye-Huckel potential in plasma physics or Thomas-Fermi potential in solid state physics. In these applications, the interpretation of the effective parameters α g and M in Eq. (1) is different. In a cloud of charged ions and electrons at a temperature T , the Coulomb field produced by an ion is obtained from the potential in Eq. (1) with α g = Ze 2 /4π and D = 1/M describes the Debye screening distance in the system given by D = [k B T /n e e 2 (1 + Z 2 )] 1/2 , where n e is the electron density. For dopped semiconductors with injected carriers, α g = e 2 /4πκ where κ stands for the dielectric constant and λ = 1/M is the Thomas-Fermi screening length due to the injected carriers [9] [10].
The problem of the bound states of the Yukawa potential has been also considered using numerical methods [31] [32] [33] [34]. Numerical results and available approximate analytical solutions for the lowest lying states shows that Coulomb l-degeneracy is broken, and for a given n, states with higher l have a higher energy than lower l states. Also, numerical solutions exhibit the phenomena of cross-over, states of a given n, l having a higher energy than states with n + 1, l , which occurs close to the critical screening values (those for which a given state goes to the continuous) and are out of the reach of perturbative calculations. The value of the critical screening lengths are important for some applications, specially the ground state critical screening, and have been calculated using non-perturbative methods [29] and numerically solving the Yukawa potential for n = 0 to n = 9 [31].
In the past few years, the Yukawa interaction has been also proposed as a solution to the core-cusp problem on the density profiles of dark matter in dwarf galaxies. The N -body simulations of collisionless cold dark matter predicts halo distributions singular at the center [35] which is not observed in the data [36]. Self-interacting dark matter has been proposed as a possibility to solve this problem [37] but there is some tension with data [38] [39] which is alleviated if we consider that self-interaction is induced by the exchange of a massive mediator which, for non-relativistic dark matter, yields a Yukawa potential [40] [41]. The calculation of the dark matter profiles requires to solve the problem of classical scattering by the Yukawa potential in the strong coupling regime, a work done in [42]. Yukawa interaction between two dark matter particles can also give rise to dark matter bound states (darkonium) which drives us to study of the Yukawa potential at the quantum level.
We became interested in this problem in the search for gauge theory for spin-one dark matter (tensor dark matter in a spinor-like formalism) [43][44] [45], whose simplest version is a hidden U (1) DM theory containing a new massive gauge boson Z . The formation of darkonium is an interesting possibility when the mediator can also couple to standard model particles. In the hidden scenario this can occur through the renormalizable kinetic mixing with the U (1) Y hypercharge [46][47] [48][49] [50]. On the experimental side, there are searches at the LHC for signatures of dark matter [51][52] looking for monojets events with a large missing energy which can be attributed to the production of a Z which decay later into a pair of dark matter particle-antiparticle system. The possibility of the creation of darkonium has been studied recently and the corresponding phenomenology depends crucially on the bound state wave function or its derivatives at the origin [ [58]. The value of the squared wave function at the origin for the ground state, |ψ(0)| 2 , has been estimated from variational methods [19] and the obtained value used in phenomenological analysis [57].
In a previous letter, we proposed a procedure to solve this long standing problem [59]. Here, we give further details of the calculations and present a complete study of the phenomenology of the bound states of the Yukawa potential. The solution is based on the hidden supersymmetry of the Yukawa potential and on a perturbative expansion of the corresponding superpotentials in powers of δ. We find that the quantum Yukawa problem is factorizable according to [60], up to order δ 2 . At the supersymmetry level this means that the supersymmetric partner of the Yukawa Hamiltonian H l belongs to the same family with different l and we have "shape invariance" [61] of the effective potential v l (r). The complete spectrum can be obtained to this order following the well known supersymmetric quantum mechanics techniques. Beyond O(δ 2 ) we loose shape invariance. However, hidden supersymmetry is still present and we can use it to completely solve the problem.
Our paper is organized as follows. In the next section we reduce the Yukawa potential to a unidimensional problem and review the basics of supersymmetric quantum mechanics in order to set our conventions. Section III is devoted to the solution of the Yukawa potential to leading order using supersymmetry and shape invariance and to set the main recurrence relations for the solution to higher orders. In section IV we use the same techniques to obtain the complete solution of the Yukawa potential up to O(δ 2 ). Section V is dedicated to obtain the complete analytical solution to O(δ 3 ) constructing a family of supersymmetric Hamitonians. In section VI we generalize these results to arbitrary order in the expansion obtaining a systematic solution to the desired order. Section VII is devoted to the analysis of the results for the main observables for the Yukawa potential. Our conclusions are given in section VIII.

II. FACTORIZATION AND HIDDEN SUPERSYMMETRY OF THE YUKAWA POTENTIAL
In the following we will work in conventional units recovering the and c factors. The radial Schrodinger equation for a particle of reduced mass µ in a central potential is where l = 0, 1, 2, .. are the angular momentum eigenvalues. Using the dimensionless variable x = r/a where a is a typical scale of the system it can be rewritten as with the following shorthand notation In terms of R(r) = u(x)/x this equation is reduced to the simple form with the effective potential v l (x) ≡ l(l + 1) For the Yukawa potential in Eq.(1) we get v l (x) = l(l + 1) with the typical distance given by the Bohr radius of the system, a = µcαg ≡ a 0 and we use the dimensionless ratio δ = a 0 /D, which in the following will be also named screening length. The energy levels can also be written in the simple form Now we factorize the Yukawa potential and construct the hidden supersymmetry of the system. First we write our problem as where This Hamiltonian can be factorized in terms of the following operators A short calculation yields The factorization of H l in terms of these operators requires in which case H l = a l a † l + C(l, δ).
The partner Hamiltonian is defined asH whereṽ Now we construct the 2 × 2 Hamiltonian and charges where p = −i d dx . The charge operators satisfy which correspond to the N = 2 supersymmetry (SUSY) algebra [62] [63]. The supersymmetric Hamiltonian is related to H l andH l as follows with the associated potentials

III. EXPANSION AND SOLUTION TO LEADING ORDER
A complete analytical solution to Eq. (14) for the Yukawa potential can be build, using an expansion of the superpotential in powers of δ [59]. First we decompose the superpotential into δ-independent and δ-dependent parts Next, we expanding the effective potential in powers of δ v l (x) = l(l + 1) To leading order Eq.(24) involves only the δ-independent part of the superpotential and the Yukawa problem reduces to the well known Coulomb problem. Writing the superpotential as and inserting this expression in Eq. (27) we find the conditions g(l)(g(l) + 1) = l(l + 1), There are two solutions for g(l) g(l) = l, and g(l) = −(l + 1).
To O(δ 0 ) the partner Hamiltonian has the potential For g(l) = l we getH l = H l−1 [60] [64]. Similarly, for g(l) = −(l + 1) we obtainH l = H l+1 . We find it convenient to use g(l) = −(l + 1) and to work with the following super-potential for the Coulomb problem With this choice, the associated potentials satisfy From these relations we can see that these potentials satisfy the "shape invariance" condition [61] where Following [61] we define H 0 (l) ≡ H l and construct the family of Hamiltonians {H r (l)}, defined by with where f (k) denotes the k-times composition of f . Using shape invariance we get ) are SUSY partners and have a common spectrum. In consequence, the whole family {H r (l)} has a common spectrum whose levels are given by In terms of the principal quantum number, n = l + r + 1, we obtain where we change the label r → n. The angular momentum quantum number takes the values l = n − 1, n − 2, ..., 1, 0 and all the corresponding states have the same energy.
The eigenstate u n,l (x) with the highest value of l satisfies H n−1 u c n,n−1 (x) = a n−1 a † n−1 + c(n − 1) u c n,n−1 (x) = n,n−1 u c n,n−1 (x).
But from Eq. (32) we get thus, this state must satisfy a † n−1 u c n,n−1 = 0, a condition that can be used to obtain its explicit form. Indeed, has the solution where N n,n−1 is a normalization factor. For a given n, the eigenstates for lower values of l can be obtained recursively from u c n,n−1 (x) using the lowering operator a l This completes the solution of the Coulomb problem which is the leading term in the expansion of the Yukawa problem in powers of δ. The explicit form of the wave functions is given in Section VI where we address the general solution to the Yukawa potential and with a convenient choice in the phases of the normalization factors we recover the conventional solutions to the Coulomb problem in terms of the Laguerre associated polynomials. Once solved the leading order, we use the decomposition (25) in Eq. (24), to get the following equation for the complementary function w l (x, δ) Notice that the expansion in the right hand side (r.h.s.) of this equation starts at order δ and it is also an expansion in powers of x. The O(δ k+1 ) term on the r.h.s. is O(x k ). There is always a polynomial solution with δ-dependent coefficients for w l (x, δ) when we work to O(δ k ), with the advantage that powers of δ and x are correlated. Indeed, the general solution can be written as The coefficients required in the calculation to a given order in δ can be fixed matching powers of x on both sides of this equation.
The solution to order O(δ) can be written as Inserting these expressions in Eq (47) and comparing powers of x we find the solution a 1 = 0, y 1 (l) = 2, thus the solution to this order is given by To O(δ 2 ) the solution must be of the form Inserting these relations in Eq. (47), matching powers of x and keeping only up to O(δ 2 ) terms, we obtain the unique solution a 2 = −(l + 1)/2 and y 2 (l) = −(l + 1)(l + 3/2), thus the solution to this order is The solution to O(δ) is straightforward since in this case we just add a constant 2δ to the Coulomb potential, thus we get the Coulomb eigenstates with the corresponding energy levels shifted by 2δ.
The first non-trivial case appears at O(δ 2 ). In this case the solution is The Hamiltonian is factorized as The partner Hamiltonian reads such thatH The partner Hamiltonian is related to the original one although a new (δ-dependent) constant term is generated whereC In this case, the associated potentials satisfy and still satisfy the "shape invariance" condition Now we follow a similar procedure to the Coulomb case defining the analogous family {H r } of SUSY partners in Eq. (37) to obtain the spectrum of H 0 = H l as which when written in terms of the principal quantum number, n = l + r + 1, reads The eigenstate u n,l (x) for l = n − 1 satisfies H n−1 u n,n−1 (x) = a n−1 a † n−1 + C(n − 1, δ) u n,n−1 (x) = n,n−1 u n,n−1 (x).
which can be solved to obtain Using Eqs. (58,59) we obtain recursively states with lower l acting with the operator a l u n,n−k (x) = N n,n−k (δ) a n−k u n,n−k+1 .
This yields the complete solution of the quantum Yukawa problem to order δ 2 . We remark that actually the systematic calculation of the eigenstates to order δ 2 requires the expansion of both the exponential and the normalization factors in Eqs. (7071). We will address the details of this expansion in the general case considered in Section VI, but it is clear that this procedure will replace these factors by polynomials of order δ 2 with x-dependent coefficients.
To O(δ 3 ) the superpotential is given by The Hamiltonian is factorized as where The partner Hamiltonian now is given bỹ At O(δ 3 ), the Hamiltonians H l and H l−1 are not longer connected by the factorization process and we loose the shape invariance of the potentials of the family {H l } which allowed us solve easily the problem at lower orders. Notice however that shape invariance is not necessary for supersymmetry to exist. Indeed, shape invariance holds when we have supersymmetry and the Hamiltonian family is factorizable as defined in [60]. This is not the case for the Yukawa potential beyond O(δ 2 ) but the hidden supersymetry of the system still holds and we will use it to find a complete analytical solution to the problem. First, following the pattern obtained at lower orders, we expect that the state with the highest value of l still satisfy the condition a † n−1 u n,n−1 = 0 which yields A direct calculation shows that this function indeed satisfies Eq. (5) with eigenvalue The state u n,n−2 (x, δ) is no longer obtained as a n−2 u n,n−1 (x, δ) as a direct calculation shows. This is due to the fact that H l and H l−1 are not longer connected by the factorization process. However, H l andH l are SUSY partners, thus they are isospectral, and we can try to solveH l in whose case we would be finding also the eigenvalues of H l . With this aim, let us take a closer look to the SUSY partner of H l . This will be the first step towards the complete solution to this order, thus we denote the partner asH (1) l in the following. It is given bỹ Since H l andH (1) l are SUSY partners they must have a common spectrum. We can solveH (1) l at least for the highest allowed value of l following the same procedure used to solve H l for l = n − 1. First we re-factorizeH The new superpotentialW andw (1) l (x, δ) must satisfy an equation similar to Eq.(47), but with the corresponding expansion ofṽ (1) l on the right hand side, namely This equation can be used to obtainw The solution of this potential for l = n − 2 is and the corresponding energy is Now we can find the eigenstate of H l for l = n − 2 using the SUSY connection and the double factorizatioñ n−2 ) † +C (1) (n − 2, δ) = a † n−2 a n−2 + C(n − 2, δ).
The stateũ (1) n,n−2 (x) satisfies [a † n−2 a n−2 + C(n − 2, δ)]ũ (1) n,n−2 =˜ (1) n,n−2ũ Acting on the last equation with a n−2 we get H n−2 (a n−2ũn,n−2 ) = [a n−2 a † n−2 + C(n − 2, δ)](a n−2ũn,n−2 ) =˜ (1) n,n−2 (a n−2ũ and we obtain the same energy level for H n−2 andH while the eigenstate is given by u n,n−2 = N n,n−2 (δ) a n−2ũ Eigenstates and eigenvalues for l = n−3 can be calculated applying now this procedure toH (1) l . The SUSY partner, denotedH We re-factorize this Hamiltonian asH withã Following the same procedure we obtaiñ v (2) l (x) = w c (l + 2) + − 1 2 (l + 3)δ 2 + 1 6 (l + 2)(l + 3)(l + 9)δ 3 x + 1 6 (l + 3)δ 3 x 2 , The common eigenvalue of the three Hamiltonians {H n−3 ,H The corresponding eigenstates are given bỹ n,n−3 (x) = N n,n−3 (δ)a n−3ã Continuing this process we will eventually reach the lowest l = 0 level, completely solving the Yukawa problem to order O(δ 3 ). The complete set of eigenvalues to O(δ 3 ) is given by The eigenstate for a given l = n − k is obtained constructing the family of supersymmetric Hamiltonians {H (r) l }, r = 0, 1, ..., k − 1 following the procedure outlined above where H (0) l = H l . It is given by u n,n−k (x) = N n,n−k (δ)a n−kã This procedure yields the complete set of eigenvalues and eigenstates to the Yukawa problem to order δ 3 . Concerning eigenstates, a systematic calculation to order δ 3 requires to expand the solution to this order, which replaces the product of the normalization factor and the exponential inũ (k−1) n,n−k (x) with a polynomial of order δ k with x-dependent coefficients. In the next section we address this procedure.

A. Energy levels
The algorithm outlined in the previous section can be applied to any order of the expansion of the Yukawa potential in powers of δ. This yields a complete analytical solution to the quantum Yukawa problem. Our claculations show that the shape invariance property of the potentials of some supersymmetric systems is a convenient property which when it holds make results based on supersymmetry easier to obtain, but it is not a necessary condition for a supersymmetric systems. Indeed, it is a condition for supersymmetric and factorizable potentials as defined in Ref. [60]. Yukawa potential is not factorizable in this sense, but expanding the potential we can build and use supersymmetry recursively, which allows us to completely solve the problem to any order in δ. We find that the energy levels depend in general of n 2 and l(l + 1). The Taylor series for the eigenvalues can be written as n,l (δ) = ∞ k=0 ε k (n 2 , l(l + 1))δ k .
The coefficients ε k (a, b) can be recursively calculated. In the analysis below we will use this expansion to order δ 52 for the calculation of some observables. However, expressions for the coefficients are rather long and for future reference we list them here only up to k = 10, given by

B. Eigenstates
The wave functions are obtained according to Eq.(104) with the superpotentialsW (m) l (x, δ) calculated to the desired order. For a given n, the tower of states starts with u n,n−1 simply given by the condition a † n−1 u n,n−1 = 0, with the solution u n,n−1 (x) = N n,n−1 (δ)x n e − x n e − wn−1(x,δ)dx with the complementary superpotential w l (x, δ) calculated to the desired order k. Strictly speaking, this solution is valid to order k, thus we must expand the exponential to obtain u n,n−1 (x) = N n,n−1 (δ)x n e − x n P k n,n−1 (x, δ), where P k n,n−1 (x, δ) is a polynomial of order k in δ with x-dependent coefficients, which in turn, are polynomials in x. This structure is preserved when we go to the next levels. In general the solutions to order k have the following structure u n,l (x) = N n,l (δ)x l+1 e − x n P k n,l (x, δ).
The states must be normalized as Performing this integration we get ∞ 0 |u n,l (x, δ)| 2 = |N n,l (δ)| 2 n 2 where K nl (δ) is a polynomial of order δ 2k . The normalized states are where η nl are phases. A systematic calculation of the eigenstates to order δ k requires to expand also the normalization factor in the above expressions to this order which yields where M k n,l (x, δ) are new polynomials of order k in δ which incorporate additional δ terms from the normalization factors. Choosing the phases as η nl = (−1) n−l−1 , in the δ → 0 limit these polynomials reduce to M k n,l (x, 0) = L 2l+1 where L 2l+1 n−l−1 (ρ) stands for the Laguerre associated polynomials. With the above choice of phases we recover the well known solutions to the Coulomb problem for δ = 0.
Sumarizing, in terms of ρ = 2x/n = 2r/na 0 , the eigenstates of the Yukawa potential to order δ k are given by where The explicit form of these polynomials are easily calculated to the desired order k. For future reference, in the Appendix we list these polynomials to order δ 5 for the lowest lying levels, n = 1, 2, 3, 4.

VII. DISCUSSION
A. Comparison with existing analytical results.
There are many calculations with partial results of the energy levels of the Yukawa potentials in the literature but most of them focus on the ground state or the first few levels. Also, only a few of them attempt a systematic expansion and yield results to higher order with which we can compare our results. We skip comparison with calculations based on variational principles for this reason. Analytical results to order δ 2 for the Yukawa potential were obtained in [16] using a variation of perturbation theory named analytic perturbation theory. Our results to order δ 2 reproduce results in that work. Using first order conventional perturbation theory with the Coulomb potential as the unperturbed system and approximate calculation of the involved matrix elements, the energy levels of the Yukawa potential were calculated in [5] up to order δ 5 . Our results agree with results in that work up to order δ 3 . The calculation of the ground state energy has been addressed by many authors. The most elaborated perturbative calculation for n = 1 has been done in [21] [22] using Logarithmic Perturbation Theory. Our results for the ground state agree with these calculations. A more comprehensive calculation of levels was done in [28] using related techniques, where calculations for n = 1, 2, 3 levels are done up to order δ 4 . Our calculation also reproduce these results.

B. Analysis of the energy levels
The coefficients in Eqs. (106) grow with k. This raises the concern on the convergence radii of the series for the energy levels in Eq.(105). The convergence of Taylor series arising in perturbative calculations has been studied in [65][66] [67][68] [69] and several methods have been devised to construct analytical extensions beyond the corresponding convergence radii [70]. We will work here with the Padé approximants technique [71]. For the energy level expansion in Eq.(105) to order k = M + N it is always possible to construct a rational function with P N (δ) and Q M (δ) polynomials of order N and M respectively. The Taylor expansion of the [N/M ] Padé approximant coincides with the expansion of nl (δ) to order k = M + N . The coefficients of the P N and Q M polynomials, are fixed by the coefficients ε k (n 2 , l(l + 1)) and these polynomials are unique. The value of the series is bounded from above and below by the values of the [(N + 1)/N ] and [N/N ] approximants [70]. As we go to high N , these approximants converge and we can estimate the the uncertainty in the calculation of the levels as the difference in their values for a given δ. In Fig. 1 we plot the energy levels for the lower states as functions of the screening length D (in units of the Bohr radius) calculated to order k = 0, 3, 6, 9, and the [(N + 1)/N ] Padé approximant for N = 6. We use the Mathematica package for the calculation of the Padé approximants in this paper. We can see in these plots that the Taylor  In general, whenever we are far from the critical region, results to order δ 3 yield a good description of the system. The approximation at this order fails as we approach the critical region and results for the Taylor series at higher orders do not improve this description. Here, the use of the Padé approximants is necessary and a more precise calculation is reached as we increase the N of the approximant, which requires a calculation to higher order of the Taylor series. For N = 6 we reach the precision for the energy levels in the critical region obtained in the numerical solution of the Yukawa potential reported in [31], and our calculations reproduce the complete set of numerical results for n = 1, .., 9 given there. The precision in the calculation of the energy levels can be further improved taking higher N approximants, e.g., for the ground state and N = 10, which requires to calculate the Taylor series to order k = 21, we find [11/10](1) − [10/10](1) = 9 × 10 −8 and this uncertainty in the calculation is considerably lower in the low δ (large screening) region.
In general we find that, for a given n, the higher l states have higher energy. This is shown in Fig. (2), where we use the [7/6](δ) Padé approximant in the calculation of the energy levels and plot them as functions of the screening length measured in units of the Bohr radius. Also, as we go to high values of n, the gap between the energies of the l = 0 and l = n − 1 increases and eventually there is a cross-over of the energy levels, i.e., n,l > n+1,l . We find that the lowest n for which this phenomena occurs is n = 4. For low n the crossing of levels takes place near the critical regions and to clearly see it is necessary to go to higher orders in the perturbative expansion. In Fig. 3  nl and δ N nl which yields the uncertainty in the calculation. There are two important considerations in this procedure. First, each level has its own Taylor series and consequently its own set of Padé approximants in the calculation to a given order k = 2N + 1. In principle we can use a common value of N in the calculation of all nl , but this yields different uncertainties for the different levels. Second, the Padé approximants are rational functions of two polynomials, thus they have poles, some of them on the real axis. It sometimes happens that, for a given N , some of these real poles are in the physical region and close to the critical values. The reliable determination of critical screening lengths, requires to use approximants with no poles in the physical regions. With these considerations we obtain the critical screening lengths listed in Tabel I. The crossing of energy levels can also be seen in this table, where we notice that δ 43 < δ 50 . Whenever δ n,l < δ n+1,l we have this phenomena and as we go to higher n more and more n + 1 levels are crossed by the n levels.

D. Wavefunctions
The wave functions for the Yukawa potential calculated to a given order k in the expansion of the potential in powers of δ are given in Eq.(115). Notice that we incorporated the δ dependent normalization factor in the N 2l+1 n−l−1 (x, δ, k) polynomials and calculate observables to O(δ k ), thus these states are normalized as This requires that the N 2l+1 n−l−1 (x, δ, k) polynomials, even if they are δ-dependent, satisfy the following normalization relations ∞ 0 dρ ρ g+1 e −ρ (N g r (ρ, δ, k)) 2 = (2r + g + 1)(r + g)! r! + O(δ k+1 ).
We expect that eigenstates in Eq.(115) yield a good description of the system whenever we are not close to the critical screening region, where terms of O(δ k+1 ) become important and we loose the appropriate normalization. In Fig. 4 we show the radial probabilities for the lowest eigenstates, where we can see that proper normalization is lost for δ ≈ 3δ nl /4. Nevertheless, similarly to the case of the energy levels, we can construct the Padé approximants for the solutions to enlarge their convergence radii. In this case, the coefficients of the polynomials in δ entering the Padé approximant are in turn polynomials in x which makes the calculation of the Padé approximants more resource consuming and with conventional computing capabilities we cannot go very high in the order of the approximants. However, we can go high enough (N = 8) to reach the beginning of the critical screening region (δ ≈ 3δ nl /4) as shown in Fig. 4. As it is clear in these plots, a delocalization phenomenon starts taking place for these values of δ. It is interesting to notice that even for not so small δ the radial probabilities are quite similar to those of the Coulomb potential, screening becoming important only near the critical region. This is in contrast with screening effects in the energy levels which depart substantially for the Coulomb values even for small values of δ.
In cold dark matter phenomenology, dark matter is non-relativistic and the exchange of massive dark gauge fields in the non-relativistic domain produces a Yukawa potential with α D = g 2 D /4π and δ = a 0 /D = 2M GB /mα g , where g D stands for the dark matter-dark matter-gauge boson coupling, M GB denotes the dark gauge boson mass and m is the dark matter mass (the reduced mass is µ = m/2). Whenever δ = 2M GB /mα g < δ 10 , where δ 10 = 1.1906124207 (2) is the ground state critical screening length, there will be non-relativistic dark matter bound states (darkonium). This requires a dark gauge boson with a mass M GB < mα g δ 10 /2, which considering perturbative interactions yields a light gauge boson compared with the dark matter mass.
Darkonium phenomenology for s-wave bound states involves the value of the wave function at the origin ψ n00 (0) = R n0 (0)/ √ 4π, while for p-wave, observables depend on the value of the radial derivative of the wave function evaluated  [31](see also [33]. n l N δ nl Numeric n l N δ nl Numeric n l N δ nl Numeric 1 0 21 1.1906124207 (2)  at the the origin ψ n10 (0) = 3/4πR n1 (0) [72]. A calculation to order δ 10 for the lower states yields The large coefficients in these Taylor series raises again the concern on the corresponding convergence radii. The specific value of δ depends on the values of M GB , the fine structure constant α D and the dark matter mass m, and it is well possible to have large values of δ. In order to have a reliable estimate in the large δ region we must resort again to analytic methods to enlarge the convergence radii, like the Padé approximants method. In Figure 5 we show result to different orders in the perturbative expansion as well as results using the [5/5](δ) Padé approximant. It has been shown in [59] that variational methods underestimate the value of the squared wave function at the origin for the ground state when compared with the same order as calculated in our formalism. In Fig. 5 we can see the effect of higher orders terms summed up by the Padé approximants. The actual value of this observable decreases with increasing δ and vanishes at the critical value δ 10 . Similar results hold for s-wave higher states and the derivatives for p-waves. This could be important for darkonium phenomenology since, as mentioned above, for not so light mediator it is well possible that δ be in the critical region in whose case we will have weak transition matrix elements for darkonium.

VIII. CONCLUSIONS
In this work we give a detailed phenomenological analysis of the complete solution to the bound states of the Yukawa potential introduced in [59]. Eigenstates, ψ nlm , and eigenvalues, nl , are obtained using the hidden supersymmetry of the system and a Taylor expansion of the Yukawa potential in terms of the parameter δ = a 0 /D. The solutions for the eigenvalues nl are given as Taylor series in δ whose coefficients depend on n 2 and l(l + 1). The eigenstates have a similar form to the solutions of the Coulomb potential but with the Laguerre Associated polynomials, L 2l+1 n−l−1 (ρ), replaced with new polynomials of order k in δ, N 2l+1 n−l−1 (ρ, δ, k), whose coefficients are polynomial in ρ = 2r/na 0 . We provide analytical expressions up to order δ 10 for the eigenvalues and to order δ 5 for the eigenstates but both can be easily calculated to the desired order k.
Our results are in principle valid for small values of δ (large values of the screening) but since our procedure allows to easily calculate eigenvalues to high orders, the convergence radii of the Taylor series can be enlarged up to the critical screening values using the Padé approximants technique. This procedure permit us to calculate eigenvalues with high precision, which we estimate from the properties of the Padé approximants, in the whole range of values of δ where bound states exist. We find that, for a given n, states with higher l are more energetic such that for n ≥ 4 energy levels exhibit a cross-over phenomenon. The critical screening lengths, δ nl , can also be precisely calculated for every state and we give numerical results them up to n = 9 using the expansion to very high order( from δ 43 to δ 53 , depending of the state).
Similar considerations on the convergence radii apply to the eigenstates of the Yukawa potential. In this case, the information on the screening effects is contained in the new polynomials N 2l+1 n−l−1 (ρ, δ, k). We provide analytical expressions up to k = 5 for these polynomials. Higher order can be easily calculated but yields long expressions and are not shown here. We calculate the radial probabilities, finding that screening has sizable effects only for values of δ close to the critical ones. The convergence radii of the perturbative expansion of the new polynomials, can also be improved using the corresponding Padé approximants. In this case, the coefficients of the rational functions entering the Padé approximants are ρ-dependent which makes numerical computations more resource consuming. We provide results for N = 8 (meaning a calculation to order δ 16 ) for the radial probabilities which allows us to reach screenings of the order δ ≈ 3δ nl /4. For these values, a delocalization effect is already visible in the radial probabilities.
The squared absolute value at r = 0 of the wave function for s-waves and its derivative for p-waves, enter the phenomenology of the decay and production of Yukawa bound states. In particular, for cold dark matter with a gauge structure and light dark mediators, Yukawa type darkonium could exist and these quantities are relevant to asses the possibilities of detecting darkonium in hadron colliders, direct or indirect detection experiments. We calculate |ψ(0)| 2 and |ψ (0)| 2 for l = 0 and l = 1 respectively, for the lowest energy levels and in the whole range of values of the screening lengths, using the Padé approximants technique.
We wrote a Mathematica code for the calculation of eigenvalues and eigenstates of the Yukawa potential to the desired order, available upon request.

IX. APPENDIX
The polynomials entering the solution of the Yukawa potential to order k = 5 for n = 1, 2, 3, 4 are given by