Bound states of the Yukawa potential from hidden supersymmetry
BBound states of the Yukawa potential from hidden supersymmetry
M. Napsuciale ∗ Departamento de F´ısica, Universidad de Guanajuato, Lomas del Campestre 103,Fraccionamiento Lomas del Campestre, Le´on, Guanajuato, M´exico, 37150.
S. Rodr´ıguez † Facultad de Ciencias F´ısico-Matem´aticas, Universidad Aut´onoma de Coahuila,Edificio A, Unidad Camporredondo, 25000, Saltillo, Coahuila, M´exico.
In this work, we present a phenomenological study of the complete analytical solution to thebound eigenstates and eigenvalues of the Yukawa potential obtained previously using the hiddensupersymmetry of the system and a systematic expansion of the Yukawa potential in terms of δ = a /D , where a is the Bohr radius and D is the screening length. The eigenvalues, (cid:15) nl ( δ ), aregiven in the form of Taylor series in δ which can be systematically calculated to the desired order δ k .Coulomb l -degeneracy is broken by the screening effects and, for a given n , (cid:15) nl ( δ ) is larger for highervalues of l which causes the crossing of levels for n ≥
4. The convergence radius of the Taylor seriescan be enlarged up to the critical values using the Pad´e approximants technique which allows us tocalculate 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 formsimilar to the solutions of the Coulomb potential, with the associated Laguerre polynomials replacedby new polynomials of order δ 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 lengthsclose to their critical values. We use these solutions to find the squared absolute value at the originof 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 lightdark mediator. I. INTRODUCTION
The Yukawa potential is the effective non-relativistic description of the interaction of two particles due to theexchange of a massive particle of mass M . It was proposed in Ref. [1] by H. Yukawa as a low energy description of thestrong interactions between nucleons, due to the exchange of massive particles, now known as pions. The potential isgiven by V ( r ) = − α g e − Mr r , (1)where α g = g / π denotes fine-structure constant of the interaction with coupling g . The range of the interaction isgiven by D = 1 /M , also named screening length. For M = 0 we obtain the Coulomb potential and for large M wehave effective short-range interactions. The first estimate of the pion mass was done in [1] based on this potentialand 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 hightemperatures [2][3][4][5][6] nuclear physics [7], astrophysics [8] and solid state physics [9] [10][11]; and is known asDebye-Huckel potential in plasma physics or Thomas-Fermi potential in solid state physics. In these applications, theinterpretation of the effective parameters α g and M in Eq. (1) is different. In a cloud of charged ions and electrons ata temperature T , the Coulomb field produced by an ion is obtained from the potential in Eq. (1) with α g = Ze / π and D = 1 /M describes the Debye screening distance in the system given by D = [ k B T /n e e (1 + Z )] / , where n e isthe electron density. For dopped semiconductors with injected carriers, α g = e / πκ where κ stands for the dielectricconstant and λ = 1 /M is the Thomas-Fermi screening length due to the injected carriers [9] [10].The importance of the quantum Yukawa problem for many research areas of physics, motivated the search forsolutions under diverse methods. It is well known that for a finite screening length there is a finite number ofbound states [5][7] [12]. Approximate calculations for some of the energy levels exist, using variational methods [4] ∗ Electronic address: mauro@fisica.ugto.mx † Electronic address: [email protected] a r X i v : . [ h e p - ph ] F e b [8][9][13][14][15], Rayleigh-Schrodinger perturbation theory [5] [16][17] [18] [19] [20]; new perturbation schemes [21][22], and other methods [23][24][25][26][27] [28][29] [30]. The ground state energy has been calculated at very highorders in the expansion on the parameter δ = a /D = M/α g µ , where a = 1 /α g µ denotes the Bohr radius and µ stands for the reduced mass of the bounded system.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 showsthat 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 stateswith n + 1 , l (cid:48) , 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 someapplications, 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 thedensity profiles of dark matter in dwarf galaxies. The N -body simulations of collisionless cold dark matter predictshalo distributions singular at the center [35] which is not observed in the data [36]. Self-interacting dark matter hasbeen proposed as a possibility to solve this problem [37] but there is some tension with data [38][39] which is alleviatedif we consider that self-interaction is induced by the exchange of a massive mediator which, for non-relativistic darkmatter, yields a Yukawa potential [40][41]. The calculation of the dark matter profiles requires to solve the problemof classical scattering by the Yukawa potential in the strong coupling regime, a work done in [42]. Yukawa interactionbetween two dark matter particles can also give rise to dark matter bound states (darkonium) which drives us tostudy 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 matterin a spinor-like formalism) [43][44][45], whose simplest version is a hidden U (1) DM theory containing a new massivegauge boson Z (cid:48) . The formation of darkonium is an interesting possibility when the mediator can also couple tostandard 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 darkmatter [51][52] looking for monojets events with a large missing energy which can be attributed to the productionof a Z (cid:48) which decay later into a pair of dark matter particle-antiparticle system. The possibility of the creation ofdarkonium has been studied recently and the corresponding phenomenology depends crucially on the bound statewave function or its derivatives at the origin [53][54][55][56][57] [58]. The value of the squared wave function at theorigin for the ground state, | ψ (0) | , has been estimated from variational methods [19] and the obtained value used inphenomenological analysis [57].In a previous letter, we proposed a procedure to solve this long standing problem [59]. Here, we give further detailsof 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 thecorresponding superpotentials in powers of δ . We find that the quantum Yukawa problem is factorizable accordingto [60], up to order δ . At the supersymmetry level this means that the supersymmetric partner of the YukawaHamiltonian H l belongs to the same family with different l and we have ”shape invariance” [61] of the effectivepotential v l ( r ). The complete spectrum can be obtained to this order following the well known supersymmetricquantum mechanics techniques. Beyond O ( δ ) we loose shape invariance. However, hidden supersymmetry is stillpresent 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 problemand review the basics of supersymmetric quantum mechanics in order to set our conventions. Section III is devoted tothe solution of the Yukawa potential to leading order using supersymmetry and shape invariance and to set the mainrecurrence relations for the solution to higher orders. In section IV we use the same techniques to obtain the completesolution of the Yukawa potential up to O ( δ ). Section V is dedicated to obtain the complete analytical solution to O ( δ ) constructing a family of supersymmetric Hamitonians. In section VI we generalize these results to arbitraryorder in the expansion obtaining a systematic solution to the desired order. Section VII is devoted to the analysis ofthe 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 (cid:126) and c factors. The radial Schrodinger equationfor a particle of reduced mass µ in a central potential is (cid:20) − (cid:126) µ (cid:18) r ddr ( r ddr ) − l ( l + 1) r (cid:19) + V ( r ) (cid:21) R ( r ) = E R ( r ) , (2)where l = 0 , , , .. are the angular momentum eigenvalues. Using the dimensionless variable x = r/a where a is atypical scale of the system it can be rewritten as (cid:20) − x ddx ( x ddx ) + l ( l + 1) x + U ( x ) (cid:21) R ( ax ) = (cid:15) R ( ax ) , (3)with the following shorthand notation U ( x ) = 2 µa (cid:126) V ( ax ) , (cid:15) = 2 µa (cid:126) E. (4)In terms of R ( r ) = u ( x ) /x this equation is reduced to the simple form (cid:20) − d dx + v l ( x ) (cid:21) u l ( x ) = (cid:15) l u l ( x ) , (5)with the effective potential v l ( x ) ≡ l ( l + 1) x + U ( x ) . (6)For the Yukawa potential in Eq.(1) we get v l ( x ) = l ( l + 1) x − x e − δx , (7)with the typical distance given by the Bohr radius of the system, a = (cid:126) µcα g ≡ a and we use the dimensionless ratio δ = a /D , which in the following will be also named screening length. The energy levels can also be written in thesimple form E l = 12 µc α g (cid:15) l . (8)Now we factorize the Yukawa potential and construct the hidden supersymmetry of the system. First we write ourproblem as H l u l = (cid:15) l u l (9)where H l = − d dx + v l ( x ) . (10)This Hamiltonian can be factorized in terms of the following operators a l = − ddx + W l ( x ) , a † l = ddx + W l ( x ) . (11)A short calculation yields a l a † l = − d dx + W l ( x ) − W (cid:48) l ( x ) , (12) a † l a l = − d dx + W l ( x ) + W (cid:48) l ( x ) . (13)The factorization of H l in terms of these operators requires W l ( x ) − W (cid:48) l ( x ) = v l ( x ) − C ( l, δ ) , (14)in which case H l = a l a † l + C ( l, δ ) . (15)The partner Hamiltonian is defined as ˜ H l = a † l a l + C ( l, δ ) = − d dx + ˜ v l ( x ) . (16)where ˜ v l ( x ) = W l ( x ) + W (cid:48) l ( x ) + C ( l, δ ) = v l ( x ) + 2 W (cid:48) l ( x ) . (17)Now we construct the 2 × H = (cid:18) a † l a l a l a † l (cid:19) = (cid:18) − d dx + W l ( x ) (cid:19) + W (cid:48) l ( x ) σ , (18)and charges Q = 12 ( σ p + σ W l ) = (cid:18) − ia l ia † l (cid:19) (19) Q = 12 ( σ p + σ W l ) = (cid:18) a l a † l (cid:19) , (20)where p = − i ddx . The charge operators satisfy { Q i , Q j } = 2 δ ij H, [ Q i , H ] = 0 , (21)which correspond to the N = 2 supersymmetry (SUSY) algebra [62][63]. The supersymmetric Hamiltonian is relatedto H l and ˜ H l as follows H = (cid:18) ˜ H l + C ( l, δ ) 00 H l + C ( l, δ ) (cid:19) = (cid:32) − d dx + U + ( x, l ) 00 − d dx + U − ( x, l ) (cid:33) , (22)with the associated potentials U ± ( x, l ) = W l ( x ) ± W (cid:48) l ( x ) . (23) III. EXPANSION AND SOLUTION TO LEADING ORDER
A complete analytical solution to Eq.(14) for the Yukawa potential W l ( x, δ ) − W (cid:48) l ( x, δ ) + C ( l, δ ) = l ( l + 1) x − x e − δx . (24)can be build, using an expansion of the superpotential in powers of δ [59]. First we decompose the superpotentialinto δ -independent and δ -dependent parts W l ( x, δ ) = w c ( x, l ) + w l ( x, δ ) , C ( l, δ ) = c ( l ) + y ( l, δ ) . (25)Next, we expanding the effective potential in powers of δv l ( x ) = l ( l + 1) x − x + 2 δ − δ x + 13 δ x + .... (26)To leading order Eq.(24) involves only the δ -independent part of the superpotential w c ( x, l ) − w (cid:48) c ( x, l ) + c ( l ) = l ( l + 1) x − x , (27)and the Yukawa problem reduces to the well known Coulomb problem. Writing the superpotential as w c ( x, l ) = g ( l ) x − g ( l ) , (28)and inserting this expression in Eq. (27) we find the conditions g ( l )( g ( l ) + 1) = l ( l + 1) , c ( l ) = − g ( l ) . (29)There are two solutions for g ( l ) g ( l ) = l, and g ( l ) = − ( l + 1) . (30)To O ( δ ) the partner Hamiltonian has the potential˜ v l ( x ) = v l ( x ) + 2 w (cid:48) c ( x, l ) = l ( l + 1) − g ( l ) x − x . (31)For g ( l ) = l we get ˜ H l = H l − [60][64]. Similarly, for g ( l ) = − ( l + 1) we obtain ˜ H l = H l +1 . We find it convenient touse g ( l ) = − ( l + 1) and to work with the following super-potential for the Coulomb problem w c ( x, l ) = 1 l + 1 − l + 1 x , c ( l ) = − l + 1) . (32)With this choice, the associated potentials satisfy U + ( x, l ) ≡ w c ( x, l ) + w (cid:48) c ( x, l ) = v l +1 ( x ) − c ( l ) , (33) U − ( x, l + 1) ≡ w c ( x, l + 1) − w (cid:48) c ( x, l + 1) = v l +1 ( x ) − c ( l + 1) . (34)From these relations we can see that these potentials satisfy the ”shape invariance” condition [61] U + ( x, l ) = U − ( x, f ( l )) + R ( f ( l )) (35)where f ( l ) = l + 1 , R ( f ( l )) = c ( f ( l )) − c ( f ( l − . (36)Following [61] we define H ( l ) ≡ H l and construct the family of Hamiltonians { H r ( l ) } , defined by H r ( l ) = − d dx + U − ( x, f ( r ) ( l )) + r (cid:88) k =1 R ( f ( k ) ( l )) (37)with f ( k ) ( l ) = f ( f ( k − ( l )) = f ( f ( f ( k − ( l ))) = ... = f ( k ) ( l ) = l + k, (38)where f ( k ) denotes the k -times composition of f . Using shape invariance we get H r +1 ( l ) = − d dx + U − ( x, f ( f ( r ) ( l )) + R ( f ( f ( r ) ( l )) + r (cid:88) k =1 R ( f ( k ) ( l )) = − d dx + U + ( x, f ( r ) ( l )) + r (cid:88) k =1 R ( f ( k ) ( l )) , (39)thus H r ( l ) − (cid:80) rk =1 R ( f ( k ) ( l )) and H r +1 ( l ) − (cid:80) rk =1 R ( f ( k ) ( l )) are SUSY partners and have a common spectrum. Inconsequence, the whole family { H r ( l ) } has a common spectrum whose levels are given by (cid:15) r,l = c ( l ) + r (cid:88) k =1 R ( f ( k ) ( l )) = − l + 1) + r (cid:88) k =1 (cid:18) l + k ) − l + k + 1) (cid:19) = − l + r + 1) . (40)In terms of the principal quantum number, n = l + r + 1, we obtain (cid:15) n,l = − n , (41)where we change the label r → n . The angular momentum quantum number takes the values l = n − , n − , ..., , u n,l ( x ) with the highest value of l satisfies H n − u cn,n − ( x ) = (cid:16) a n − a † n − + c ( n − (cid:17) u cn,n − ( x ) = (cid:15) n,n − u cn,n − ( x ) . (42)But from Eq. (32) we get c ( n −
1) = − n = (cid:15) n,n − , (43)thus, this state must satisfy a † n − u cn,n − = 0, a condition that can be used to obtain its explicit form. Indeed, (cid:20) ddx + w c ( n − , x ) (cid:21) u cn,n − ( x ) = 0 , (44)has the solution u cn,n − ( x ) = N n,n − e − (cid:82) w cn − ( x ) dx = N n,n − x n e − xn , (45)where N n,n − is a normalization factor. For a given n , the eigenstates for lower values of l can be obtained recursivelyfrom u cn,n − ( x ) using the lowering operator a l u cn,n − k ( x ) = N n,n − k (cid:18) − ddx + w c ( x, n − k ) (cid:19) u cn,n − k +1 ( x ) . (46)This completes the solution of the Coulomb problem which is the leading term in the expansion of the Yukawaproblem in powers of δ . The explicit form of the wave functions is given in Section VI where we address the generalsolution to the Yukawa potential and with a convenient choice in the phases of the normalization factors we recoverthe 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 thecomplementary function w l ( x, δ ) w l ( x, δ ) − w (cid:48) l ( x, δ ) + 2 w c ( x, l ) w l ( x, δ ) = 2 δ − δ x + 13 δ x + ... + y ( l, δ ) . (47)Notice that the expansion in the right hand side (r.h.s.) of this equation starts at order δ and it is also an expansionin powers of x . The O ( δ k +1 ) term on the r.h.s. is O ( x k ). There is always a polynomial solution with δ -dependentcoefficients 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 w l ( x, δ ) = a δ + ( a δ + a δ + a δ ... ) x + ( b δ + b δ + b δ + ... ) x + ( c δ + c δ + c δ + ... ) x + ..., (48) y ( l, δ ) = y ( l ) δ + y ( l ) δ + y ( l ) δ + ... (49)The coefficients required in the calculation to a given order in δ can be fixed matching powers of x on both sides ofthis equation. IV. SOLUTION OF THE YUKAWA PROBLEM TO ORDER δ . The solution to order O ( δ ) can be written as w l ( x, δ ) = a δ, y ( l, δ ) = y ( l ) δ. (50)Inserting these expressions in Eq(47) and comparing powers of x we find the solution a = 0, y ( l ) = 2, thus thesolution to this order is given by O ( δ ) : w l ( x, δ ) = 0 , y ( l, δ ) = 2 δ. (51)To O ( δ ) the solution must be of the form w l ( x, δ ) = a δ x, y ( l, δ ) = 2 δ + y ( l ) δ . (52)Inserting these relations in Eq.(47), matching powers of x and keeping only up to O ( δ ) terms, we obtain the uniquesolution a = − ( l + 1) / y ( l ) = − ( l + 1)( l + 3 / O ( δ ) : w l ( x, δ ) = −
12 ( l + 1) δ x, y ( l, δ ) = 2 δ − ( l + 1)( l + 32 ) δ . (53) A. Energy levels at O ( δ ) The solution to O ( δ ) is straightforward since in this case we just add a constant 2 δ to the Coulomb potential, thuswe get the Coulomb eigenstates with the corresponding energy levels shifted by 2 δ .The first non-trivial case appears at O ( δ ). In this case the solution is W l ( x, δ ) = − l + 1 − l + 1 x −
12 ( l + 1) δ x, C ( l, δ ) = − l + 1) + 2 δ − ( l + 1)( l + 32 ) δ . (54)The Hamiltonian is factorized as H l = − d dx + l ( l + 1) x − x + 2 δ − δ x = a l a † l + C ( l, δ ) . (55)The partner Hamiltonian reads˜ H l = a † l a l + C ( l, δ ) = − d dx + ( l + 1)( l + 2) x − x + 2 δ − δ x − ( l + 1) δ , (56)such that ˜ H l − = a † l − a l − + C ( l − , δ ) = − d dx + l ( l + 1) x − x + 2 δ − δ x − lδ = H l − lδ . (57)The partner Hamiltonian is related to the original one although a new ( δ -dependent) constant term is generated a l a † l = H l − C ( l, δ ) , (58) a † l − a l − = H l − C ( l − , δ ) − lδ = H l − ˜ C ( l − , δ ) , (59)where ˜ C ( l, δ ) = C ( l, δ ) + ( l + 1) δ . (60)In this case, the associated potentials satisfy U + ( x, l ) ≡ W l ( x, δ ) + W (cid:48) l ( x, δ ) = v l +1 ( x, δ ) − ˜ C ( l, δ ) , (61) U − ( x, l + 1) ≡ W l +1 ( x, δ ) − W (cid:48) l +1 ( x, δ ) = v l +1 ( x, δ ) − C ( l + 1 , δ ) , (62)and still satisfy the ”shape invariance” condition U + ( x, l ) = U − ( x, f ( l )) + R ( f ( l )) (63)with f ( l ) = l + 1 , R ( f ( l )) = C ( f ( l ) , δ ) − ˜ C ( f ( l − , δ ) . (64)Now we follow a similar procedure to the Coulomb case defining the analogous family { H r } of SUSY partners inEq.(37) to obtain the spectrum of H = H l as (cid:15) r,l = C ( l, δ ) + r (cid:88) k =1 R ( f ( k ) ( l )) = − l + r + 1) + 2 δ − (cid:20) ( l + 1)( l + 32 ) + 3 r ( r + 2( l + 1)) (cid:21) δ , (65)which when written in terms of the principal quantum number, n = l + r + 1, reads (cid:15) n,l = − n + 2 δ −
12 [3 n − l ( l + 1)] δ . (66) B. Eigenstates to O ( δ ) The eigenstate u n,l ( x ) for l = n − H n − u n,n − ( x ) = (cid:16) a n − a † n − + C ( n − , δ ) (cid:17) u n,n − ( x ) = (cid:15) n,n − u n,n − ( x ) . (67)Comparing Eqs.(54,66) we obtain C ( n − , δ ) = − n + 2 δ − n ( n + 12 ) = (cid:15) n,n − . (68)This state must satisfy a † n − u n,n − = 0, i.e. (cid:20) ddx + W n − ( x, δ ) (cid:21) u n,n − ( x ) = 0 , (69)which can be solved to obtain u n,n − ( x, δ ) = N n,n − ( δ ) e − (cid:82) W n − ( x,δ ) dx = N n,n − ( δ ) x n e − xn + nδ x . (70)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 . (71)This yields the complete solution of the quantum Yukawa problem to order δ . We remark that actually the sys-tematic calculation of the eigenstates to order δ requires the expansion of both the exponential and the normalizationfactors in Eqs. (7071). We will address the details of this expansion in the general case considered in Section VI, butit is clear that this procedure will replace these factors by polynomials of order δ with x -dependent coefficients. V. SOLUTION TO THE YUKAWA POTENTIAL TO O ( δ ) . To O ( δ ) the superpotential is given by W l ( x ) = 1 l + 1 − l + 1 x − (cid:20)
12 ( l + 1) δ + 16 ( l + 1) ( l + 2) δ (cid:21) x + l + 16 δ x . (72)The Hamiltonian is factorized as H l = − d dx + l ( l + 1) x − x + 2 δ − δ x + 13 δ x = a l a † l + C ( l, δ ) . (73)where C ( l, δ ) = − l + 1) + 2 δ − ( l + 1)( l + 32 ) δ + 13 ( l + 1) ( l + 2)( l + 32 ) δ . (74)The partner Hamiltonian now is given by˜ H l − = a † l − a l − + C ( l − , δ ) = H l − lδ + 13 l [ l ( l + 1) + 2 x ] δ . (75)At O ( δ ), the Hamiltonians H l and H l − are not longer connected by the factorization process and we loose the shapeinvariance of the potentials of the family { H l } which allowed us solve easily the problem at lower orders. Noticehowever that shape invariance is not necessary for supersymmetry to exist. Indeed, shape invariance holds when wehave supersymmetry and the Hamiltonian family is factorizable as defined in [60]. This is not the case for the Yukawapotential beyond O ( δ ) but the hidden supersymetry of the system still holds and we will use it to find a completeanalytical solution to the problem. First, following the pattern obtained at lower orders, we expect that the statewith the highest value of l still satisfy the condition a † n − u n,n − = 0 which yields u n,n − ( x, δ ) = N n,n − ( δ ) e − (cid:82) W n − ( x,δ ) dx = N n,n − ( δ ) x n e − xn e [ n δ − n ( n +1) δ ] x − n δ x . (76)A direct calculation shows that this function indeed satisfies Eq. (5) with eigenvalue (cid:15) n,n − = − n + 2 δ − n ( n + 12 ) δ + 13 n ( n + 1)( n + 12 ) δ = C ( n − , δ ) . (77)The state u n,n − ( x, δ ) is no longer obtained as a n − u n,n − ( x, δ ) as a direct calculation shows. This is due to the factthat H l and H l − are not longer connected by the factorization process. However, H l and ˜ H l are SUSY partners,thus they are isospectral, and we can try to solve ˜ H 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 completesolution to this order, thus we denote the partner as ˜ H (1) l in the following. It is given by˜ H (1) l ≡ a † l a l + C ( l, δ ) = − d dx + U + ( x, l ) + C ( l, δ ) = d dx + ˜ v (1) l ( x ) , (78)where to O ( δ ) ˜ v (1) l ( x ) = v l +1 ( x ) − ( l + 1) δ + 13 ( l + 1)[( l + 1)( l + 2) + 2 x ] δ . (79)Since H l and ˜ H (1) l are SUSY partners they must have a common spectrum. We can solve ˜ H (1) l at least for the highestallowed value of l following the same procedure used to solve H l for l = n −
1. First we re-factorize ˜ H (1) l as follows˜ H (1) l = ˜ a (1) l (˜ a (1) l ) † + ˜ C (1) ( l, δ ) (80)where ˜ a (1) l = − ddx + ˜ W (1) l ( x ) , (˜ a (1) l ) † = ddx + ˜ W (1) l ( x ) . (81)The new superpotential ˜ W (1) l ( x ) is decomposed as˜ W (1) l ( x, δ ) = w c ( x, l + 1) + ˜ w (1) l ( x, δ ) , (82)and ˜ w (1) l ( x, δ ) must satisfy an equation similar to Eq.(47), but with the corresponding expansion of ˜ v (1) l on the righthand side, namely( ˜ w (1) l ( x, δ )) − ( ˜ w (1) l ( x, δ )) (cid:48) + 2 w c ( x, l + 1) ˜ w (1) l ( x, δ ) = 2 δ − ( l + 1) δ + 13 ( l + 1) ( l + 2) δ − ( δ −
23 ( l + 1) δ ) x + 13 δ x . (83)This equation can be used to obtain ˜ w (1) l ( x, δ ) and ˜ C (1) ( l, δ ) to O ( δ ) as we did for w l ( x, δ ). We get˜ w (1) l ( x, δ ) = (cid:20) −
12 ( l + 2) δ + 16 ( l + 2)( l + 7 l + 8) δ (cid:21) x + 16 ( l + 2) δ x , (84)˜ C (1) ( l, δ ) = − l + 2) + 2 δ − ( l + 32 )( l + 4) δ + 13 ( l + 32 )( l + 2) ( l + 7) δ . (85)The solution of this potential for l = n − u (1) n,n − ( x ) = e (cid:82) ˜ W (1) n − ( x ) dx = x n e − xn e nδ − n ( n +3 n − δ x − nδ x (86)and the corresponding energy is˜ (cid:15) (1) n,n − = ˜ C (1) ( n − , δ ) = − n + 2 δ − ( n −
12 )( n + 2) δ + 13 ( n −
12 ) n ( n + 5) δ . (87)Now we can find the eigenstate of H l for l = n − H (1) n − = ˜ a (1) n − (˜ a (1) n − ) † + ˜ C (1) ( n − , δ ) = a † n − a n − + C ( n − , δ ) . (88)0The state ˜ u (1) n,n − ( x ) satisfies [ a † n − a n − + C ( n − , δ )]˜ u (1) n,n − = ˜ (cid:15) (1) n,n − ˜ u (1) n,n − . (89)Acting on the last equation with a n − we get H n − ( a n − ˜ u n,n − ) = [ a n − a † n − + C ( n − , δ )]( a n − ˜ u n,n − ) = ˜ (cid:15) (1) n,n − ( a n − ˜ u (1) n,n − ) , (90)and we obtain the same energy level for H n − and ˜ H (1) n − as expected (cid:15) n,n − = ˜ (cid:15) (1) n,n − = − n + 2 δ − ( n −
12 )( n + 2) δ + 13 ( n −
12 ) n ( n + 5) δ , (91)while the eigenstate is given by u n,n − = N n,n − ( δ ) a n − ˜ u (1) n,n − . (92)Eigenstates and eigenvalues for l = n − H (1) l . The SUSY partner,denoted ˜ H (2) l is ˜ H (2) l ≡ (˜ a (1) l ) † ˜ a (1) l + ˜ C (1) ( l, δ ) ≡ d dx + ˜ v (2) l ( x ) . (93)We re-factorize this Hamiltonian as ˜ H (2) l ≡ ˜ a (2) l (˜ a (2) l ) † + ˜ C (2) ( l, δ ) , (94)with ˜ a (2) l = − ddx + ˜ W (2) l ( x ) , (˜ a (2) l ) † = ddx + ˜ W (2) l ( x ) . (95)Following the same procedure we obtain˜ v (2) l ( x ) = v c ( l + 2) + 2 δ − (2 l + 3) δ + 23 ( l + 32 )( l + 2)( l + 3) − ( δ −
43 ( l + 32 ) δ ) x + 13 δ x , (96)˜ W (2) l ( x ) = w c ( l + 2) + (cid:20) −
12 ( l + 3) δ + 16 ( l + 2)( l + 3)( l + 9) δ (cid:21) x + 16 ( l + 3) δ x , (97)˜ C (2) ( l, δ ) = − l + 3) + 2 δ −
12 (2 l + 17 l + 27) δ + 13 ( l + 3) ( l + 2)( l + 232 ) δ . (98)The common eigenvalue of the three Hamiltonians { H n − , ˜ H (1) n − , ˜ H (2) n − } is (cid:15) n,n − = − n + 2 δ − ( n + 52 n − δ + 13 n ( n − n + 172 ) δ . (99)The corresponding eigenstates are given by˜ u (2) n,n − ( x ) = N (2) n,n − ( δ ) e (cid:82) ˜ W (3) n − ( x ) dx , (100)˜ u (1) n,n − ( x ) = N (1) n,n − ( δ )˜ a (1) n − ˜ u (2) n,n − ( x ) , (101) u n,n − ( x ) = N (0) n,n − ( δ ) a n − ˜ u (1) n,n − ( x ) = N n,n − ( δ ) a n − ˜ a (1) n − ˜ u (2) n,n − ( x ) . (102)Continuing this process we will eventually reach the lowest l = 0 level, completely solving the Yukawa problem toorder O ( δ ). The complete set of eigenvalues to O ( δ ) is given by (cid:15) n,l ( δ ) = − n + 2 δ −
12 [3 n − l ( l + 1)] δ + n n + 1 − l ( l + 1)) δ . (103)The eigenstate for a given l = n − k is obtained constructing the family of supersymmetric Hamiltonians { ˜ H ( r ) l } , r = 0 , , ..., k − H (0) l = H l . It is given by u n,n − k ( x ) = N n,n − k ( δ ) a n − k ˜ a (1) n − k ... ˜ a ( k − n − k ˜ u ( k − n,n − k ( x ) . (104)This procedure yields the complete set of eigenvalues and eigenstates to the Yukawa problem to order δ . Concerningeigenstates, a systematic calculation to order δ requires to expand the solution to this order, which replaces theproduct of the normalization factor and the exponential in ˜ u ( k − n,n − k ( x ) with a polynomial of order δ k with x -dependentcoefficients. In the next section we address this procedure.1 VI. COMPLETE SOLUTION TO THE YUKAWA POTENTIAL TO ANY ORDER IN δ .A. Energy levels The algorithm outlined in the previous section can be applied to any order of the expansion of the Yukawa potentialin powers of δ . This yields a complete analytical solution to the quantum Yukawa problem. Our claculations show thatthe shape invariance property of the potentials of some supersymmetric systems is a convenient property which whenit holds make results based on supersymmetry easier to obtain, but it is not a necessary condition for a supersymmetricsystems. Indeed, it is a condition for supersymmetric and factorizable potentials as defined in Ref. [60]. Yukawapotential 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 generalof n and l ( l + 1). The Taylor series for the eigenvalues can be written as (cid:15) n,l ( δ ) = ∞ (cid:88) k =0 ε k ( n , l ( l + 1)) δ k . (105)The coefficients ε k ( a, b ) can be recursively calculated. In the analysis below we will use this expansion to order δ forthe calculation of some observables. However, expressions for the coefficients are rather long and for future referencewe list them here only up to k = 10, given by ε = − a ,ε = 2 ,ε = −
12 (3 a − b ) ,ε = a a − b + 1) ,ε = − a
96 (77 a + 55 a − ab − b − b ) ,ε = a
160 (171 a + 245 a − ab − b − b + 4) ,ε = − a (cid:2) a − a (69 b − − a (cid:0) b + 420 b − (cid:1) − b (cid:0) b + 41 b + 6 (cid:1)(cid:3) ,ε = a (cid:2) a − a (141 b − − a (cid:0) b + 1245 b − (cid:1) − (cid:0) b + 1281 b + 686 b − (cid:1)(cid:3) ,ε = − a a − a (38034 b − − a (cid:0) b + 66312 b − (cid:1) − a (cid:0) b + 237265 b + 303534 b − (cid:1) − b (cid:0) b + 2228 b + 580 b + 48 (cid:1) ] ,ε = a a − a (619482 b − − a (cid:0) b + 359448 b − (cid:1) − a (cid:0) b + 148791 b + 381258 b − (cid:1) − (cid:0) b + 767060 b + 622580 b + 211632 b − (cid:1) ] ,ε = − a a − a (2749521 b − − a (cid:0) b + 961480800 b − (cid:1) − a (cid:0) b + 23985045 b + 113412222 b − (cid:1) − a (cid:0) b + 17424260 b + 31739120 b + 26224200 b − (cid:1) − b (cid:0) b + 223865 b + 82252 b + 12924 b + 720 (cid:1) ] . (106) B. Eigenstates
The wave functions are obtained according to Eq.(104) with the superpotentials ˜ W ( m ) l ( x, δ ) calculated to the desiredorder. For a given n , the tower of states starts with u n,n − simply given by the condition a † n − u n,n − = 0, with thesolution u n,n − ( x ) = N n,n − ( δ ) x n e − xn e − (cid:82) w n − ( x,δ ) dx (107)2with the complementary superpotential w l ( x, δ ) calculated to the desired order k . Strictly speaking, this solution isvalid to order k , thus we must expand the exponential to obtain u n,n − ( x ) = N n,n − ( δ ) x n e − xn P kn,n − ( x, δ ) , (108)where P kn,n − ( 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 followingstructure u n,l ( x ) = N n,l ( δ ) x l +1 e − xn P kn,l ( x, δ ) . (109)The states must be normalized as (cid:90) ∞ dr r | R n,l ( r, δ ) | = a (cid:90) ∞ dx | u n,l ( x, δ ) | = 1 . (110)Performing this integration we get (cid:90) ∞ | u n,l ( x, δ ) | = | N n,l ( δ ) | (cid:16) n (cid:17) ( n + l )!2 n ( n − l − (cid:16) n (cid:17) l K nl ( δ ) (111)where K nl ( δ ) is a polynomial of order δ k . The normalized states are u n,l ( x, δ ) = η nl (cid:115)(cid:18) na (cid:19) ( n − l − n + l )!2 n (cid:18) n (cid:19) l (cid:112) K nl ( δ ) x l +1 e − xn P kn,l ( x, δ ) , (112)where η nl are phases. A systematic calculation of the eigenstates to order δ k requires to expand also the normalizationfactor in the above expressions to this order which yields u n,l ( x, δ, k ) = (cid:115)(cid:18) na (cid:19) ( n − l − n + l )!2 n (cid:18) n (cid:19) l x l +1 e − xn M kn,l ( x, δ ) , (113)where M kn,l ( x, δ ) are new polynomials of order k in δ which incorporate additional δ terms from the normalizationfactors. Choosing the phases as η nl = ( − n − l − , in the δ → M kn,l ( x,
0) = L l +1 n − l − (cid:18) xn (cid:19) , (114)where L l +1 n − l − ( ρ ) stands for the Laguerre associated polynomials. With the above choice of phases we recover the wellknown solutions to the Coulomb problem for δ = 0.Sumarizing, in terms of ρ = 2 x/n = 2 r/na , the eigenstates of the Yukawa potential to order δ k are given by ψ nlm ( r , δ, k ) = (cid:115)(cid:18) na (cid:19) ( n − l − n + l )!2 n ρ l e − ρ N l +1 n − l − ( ρ, δ, k ) Y ml ( θ, φ ) , (115)where N l +1 n − l − ( ρ, δ, k ) = M kn,l ( nρ , δ ) . (116)The explicit form of these polynomials are easily calculated to the desired order k . For future reference, in theAppendix we list these polynomials to order δ for the lowest lying levels, n = 1 , , , VII. DISCUSSIONA. Comparison with existing analytical results.
There are many calculations with partial results of the energy levels of the Yukawa potentials in the literature butmost of them focus on the ground state or the first few levels. Also, only a few of them attempt a systematic expansion3and yield results to higher order with which we can compare our results. We skip comparison with calculations basedon variational principles for this reason. Analytical results to order δ for the Yukawa potential were obtained in[16] using a variation of perturbation theory named analytic perturbation theory. Our results to order δ reproduceresults in that work. Using first order conventional perturbation theory with the Coulomb potential as the unperturbedsystem and approximate calculation of the involved matrix elements, the energy levels of the Yukawa potential werecalculated in [5] up to order δ . Our results agree with results in that work up to order δ . The calculation of theground state energy has been addressed by many authors. The most elaborated perturbative calculation for n = 1has been done in [21] [22] using Logarithmic Perturbation Theory. Our results for the ground state agree with thesecalculations. A more comprehensive calculation of levels was done in [28] using related techniques, where calculationsfor n = 1 , , δ . Our calculation also reproduce these results. k = = = = [ / ]( D / a ) - - - - - / a ϵ k = = = = [ / ]( D / a ) - - - - - / a ϵ k = = = = [ / ]( D / a ) - - - - - / a ϵ k = = = = [ / ]( D / a )
10 15 20 25 30 35 40 - - - - - - / a ϵ FIG. 1: Lowest energy levels of the Yukawa potential (in Rydberg units) as functions of the screening length D , calculated atdifferent order k in the perturbative expansion, and the [( N + 1) /N ] Pad´e approximant for N = 6. 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 theenergy 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 correspondingconvergence radii [70]. We will work here with the Pad´e approximants technique [71]. For the energy level expansion4in Eq.(105) to order k = M + N it is always possible to construct a rational function[ N/M ]( δ ) = P N ( δ ) Q M ( δ ) (117)with P N ( δ ) and Q M ( δ ) polynomials of order N and M respectively. The Taylor expansion of the [ N/M ] Pad´eapproximant coincides with the expansion of (cid:15) nl ( δ ) to order k = M + N . The coefficients of the P N and Q M polynomials, are fixed by the coefficients ε k ( n , l ( l + 1)) and these polynomials are unique. The value of the series isbounded 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 differencein 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 , , ,
9, and the [( N + 1) /N ] Pad´e approximant for N = 6.We use the Mathematica package for the calculation of the Pad´e approximants in this paper. We can see in these plotsthat the Taylor series yields a good description everywhere except close to the critical screening, while the [( N + 1) /N ]and [ N/N ] approximants, enlarge the convergence radii up to the critical region. The precision in the calculation ofthe ground state in the region close to the critical screening, for N = 6 and D/a = 1 is [7 / − [6 / . × − .In general, whenever we are far from the critical region, results to order δ 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 higherorders do not improve this description. Here, the use of the Pad´e approximants is necessary and a more precisecalculation is reached as we increase the N of the approximant, which requires a calculation to higher order of theTaylor series. For N = 6 we reach the precision for the energy levels in the critical region obtained in the numericalsolution of the Yukawa potential reported in [31], and our calculations reproduce the complete set of numerical resultsfor n = 1 , .., 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 / × − 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 weuse the [7 / δ ) Pad´e approximant in the calculation of the energy levels and plot them as functions of the screeninglength 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 − (cid:15) n,l > (cid:15) n +1 ,l (cid:48) . We findthat the lowest n for which this phenomena occurs is n = 4. For low n the crossing of levels takes place near thecritical regions and to clearly see it is necessary to go to higher orders in the perturbative expansion. In Fig. 3 weshow results of the calculation of the energy levels for n = 4 , , , , N + 1) /N ] Pad´e approximant with N = 10 which requires a calculation of the series to order k = 21. C. Critical screening lengths
The critical screening lengths δ nl , defined as the values of δ at which the levels (cid:15) n,l go to the continuum, are veryimportant for practical applications. They are in the large δ region and we must use the [( N + 1) /N ] and [ N/N ]Pad´e approximants in their calculation. For a given level (cid:15) nl , we find the values of δ ≡ δ N +1 nl and δ ≡ δ Nnl at which[( N + 1) /N ]( δ N +1 nl ) = 0 and [ N/N ]( δ Nnl ) = 0. For large enough values of N , these values of δ coincide up to a givenfigure and the actual value of δ nl lies between δ N +1 nl and δ Nnl which yields the uncertainty in the calculation.There are two important considerations in this procedure. First, each level has its own Taylor series and consequentlyits own set of Pad´e approximants in the calculation to a given order k = 2 N + 1. In principle we can use a commonvalue of N in the calculation of all (cid:15) nl , but this yields different uncertainties for the different levels. Second, thePad´e approximants are rational functions of two polynomials, thus they have poles, some of them on the real axis. Itsometimes happens that, for a given N , some of these real poles are in the physical region and close to the criticalvalues. The reliable determination of critical screening lengths, requires to use approximants with no poles in thephysical regions. With these considerations we obtain the critical screening lengths listed in Tabel I. The crossing ofenergy levels can also be seen in this table, where we notice that δ < δ . Whenever δ n,l < δ n +1 ,l (cid:48) we have thisphenomena 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 powersof δ are given in Eq.(115). Notice that we incorporated the δ dependent normalization factor in the N l +1 n − l − ( x, δ, k )5 l = =
15 10 15 20 25 30 - - - - / a ϵ l l = = =
210 15 20 25 30 35 40 - - - - - - - / a ϵ l l = = = =
320 30 40 50 60 70 80 - - - - / a ϵ l l = = = = =
430 40 50 60 70 80 90 - - - - - / a ϵ l FIG. 2: Lowest energy levels of the Yukawa potential (in Rydberg units), (cid:15) nl , as functions of the screening length D (in Bohrradius units) calculated using the [( N + 1) /N ] Pad´e approximant with N = 6. polynomials and calculate observables to O ( δ k ), thus these states are normalized as (cid:90) d r | ψ nlm ( r , δ, k ) | = 1 + O ( δ k +1 ) . (118)This requires that the N l +1 n − l − ( x, δ, k ) polynomials, even if they are δ -dependent, satisfy the following normalizationrelations (cid:90) ∞ dρ ρ g +1 e − ρ ( N gr ( ρ, δ, k )) = (2 r + g + 1)( r + g )! r ! + O ( δ k +1 ) . (119)We expect that eigenstates in Eq.(115) yield a good description of the system whenever we are not close to the criticalscreening region, where terms of O ( δ k +1 ) become important and we loose the appropriate normalization. In Fig. 4we show the radial probabilities for the lowest eigenstates, where we can see that proper normalization is lost for δ ≈ δ nl /
4. Nevertheless, similarly to the case of the energy levels, we can construct the Pad´e approximants forthe solutions to enlarge their convergence radii. In this case, the coefficients of the polynomials in δ entering thePad´e approximant are in turn polynomials in x which makes the calculation of the Pad´e approximants more resourceconsuming 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 ( δ ≈ δ nl /
4) as shown6 ϵ ϵ ϵ ϵ
20 22 24 26 28 30 32 34 - - - - - / a ϵ n l ϵ ϵ ϵ ϵ ϵ
30 35 40 45 50 - - - - / a ϵ n l ϵ ϵ ϵ ϵ ϵ ϵ
40 50 60 70 80 - - - - / a ϵ n l ϵ ϵ ϵ ϵ ϵ ϵ ϵ
60 70 80 90 100 110 - - - - - - / a ϵ n l FIG. 3: Crossing of energy levels of the Yukawa potential. Curves correspond to energy levels calculated using the [( N + 1) /N ]Pad´e approximant with N = 10. in Fig. 4. As it is clear in these plots, a delocalization phenomenon starts taking place for these values of δ . It isinteresting to notice that even for not so small δ the radial probabilities are quite similar to those of the Coulombpotential, screening becoming important only near the critical region. This is in contrast with screening effects in theenergy 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 inthe non-relativistic domain produces a Yukawa potential with α D = g D / π and δ = a /D = 2 M 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 thedark matter mass (the reduced mass is µ = m/ δ = 2 M GB /mα g < δ , where δ = 1 . M GB < mα g δ /
2, which considering perturbative interactions yields a lightgauge boson compared with the dark matter mass.Darkonium phenomenology for s -wave bound states involves the value of the wave function at the origin ψ n (0) = R n (0) / √ π , while for p -wave, observables depend on the value of the radial derivative of the wave function evaluated7 TABLE I: Critical screening lengths, δ nl , for n = 1 , ...,
9, calculated with the [ N + 1 /N ] and [ N/N ] Pad´e approximants. Forcomparison we also show results from the numerical solutions given in Ref.[31](see also [33]. n l N δ nl Numeric n l N δ nl Numeric n l N δ nl Numeric1 0 21 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . at the the origin ψ (cid:48) n (0) = (cid:112) / πR (cid:48) n (0) [72]. A calculation to order δ for the lower states yields | ψ (0) | = 1 πa (cid:18) − δ + 116 δ − δ + 1427160 δ − δ + 3194474032 δ − δ + 6390853769120 δ − δ (cid:19) , (120) | ψ (0) | = 18 πa (cid:18) − δ + 3283 δ − δ + 277845 δ − δ + 3166386463 δ − δ + 7950329914135 δ − δ (cid:19) , (121) | ψ (cid:48) (0) | = 132 πa (cid:18) − δ + 3203 δ − δ + 208245 δ − δ + 1335711235 δ − δ + 276918997063 δ − δ (cid:19) . (122) | ψ (cid:48) (0) | = 8729 πa (cid:18) − δ + 1215 δ − δ + 1979307980 δ − δ + 53785680507560 δ − δ + 3345713402181037168 δ − δ (cid:19) . (123)The large coefficients in these Taylor series raises again the concern on the corresponding convergence radii. Thespecific value of δ depends on the values of M GB , the fine structure constant α D and the dark matter mass m , andit is well possible to have large values of δ . In order to have a reliable estimate in the large δ region we must resortagain to analytic methods to enlarge the convergence radii, like the Pad´e approximants method. In Figure 5 we showresult to different orders in the perturbative expansion as well as results using the [5 / δ ) Pad´e approximant. It hasbeen shown in [59] that variational methods underestimate the value of the squared wave function at the origin forthe ground state when compared with the same order as calculated in our formalism. In Fig. 5 we can see the effectof higher orders terms summed up by the Pad´e approximants. The actual value of this observable decreases withincreasing δ and vanishes at the critical value δ . 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 mediatorit is well possible that δ be in the critical region in whose case we will have weak transition matrix elements for8 FIG. 4: Radial probabilities for the lowest states as functions of δ calculated to order k = 5 in the perturbative expansion.Curves labelled as [ N/N ] correspond to the Pad´e approximant with N = 8 and δ close to the critical region. darkonium. VIII. CONCLUSIONS
In this work we give a detailed phenomenological analysis of the complete solution to the bound states of the Yukawapotential introduced in [59]. Eigenstates, ψ nlm , and eigenvalues, (cid:15) nl , are obtained using the hidden supersymmetryof the system and a Taylor expansion of the Yukawa potential in terms of the parameter δ = a /D . The solutions forthe eigenvalues (cid:15) nl are given as Taylor series in δ whose coefficients depend on n and l ( l + 1). The eigenstates havea similar form to the solutions of the Coulomb potential but with the Laguerre Associated polynomials, L l +1 n − l − ( ρ ),replaced with new polynomials of order k in δ , N l +1 n − l − ( ρ, δ, k ), whose coefficients are polynomial in ρ = 2 r/na . Weprovide analytical expressions up to order δ for the eigenvalues and to order δ for the eigenstates but both can beeasily 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 allowsto easily calculate eigenvalues to high orders, the convergence radii of the Taylor series can be enlarged up to thecritical screening values using the Pad´e approximants technique. This procedure permit us to calculate eigenvalueswith high precision, which we estimate from the properties of the Pad´e 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 ≥ [ / ]( δ ) k = = = δ π a | ψ ( ) [ / ]( δ ) k = = = δ π a | ψ ( ) [ / ]( δ ) k = = = δ π a | ψ ' ( ) [ / ]( δ ) k = = = δ π a | ψ ' ( ) / FIG. 5: Squared wave functions (for l = 0) and radial derivatives (for l = 1) evaluated at the origin for the lowest lying statesas functions of δ calculated to order k = 10 in the perturbative expansion. Curves labelled as [5 / δ ) correspond to the Pad´eapproximant calculation with N = 5. energy levels exhibit a cross-over phenomenon. The critical screening lengths, δ nl , can also be precisely calculated forevery state and we give numerical results them up to n = 9 using the expansion to very high order( from δ to δ ,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 l +1 n − l − ( ρ, δ, k ). We provide analyticalexpressions up to k = 5 for these polynomials. Higher order can be easily calculated but yields long expressions andare 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 beimproved using the corresponding Pad´e approximants. In this case, the coefficients of the rational functions enteringthe Pad´e approximants are ρ -dependent which makes numerical computations more resource consuming. We provideresults for N = 8 (meaning a calculation to order δ ) for the radial probabilities which allows us to reach screeningsof the order δ ≈ δ 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 thephenomenology of the decay and production of Yukawa bound states. In particular, for cold dark matter with a gaugestructure and light dark mediators, Yukawa type darkonium could exist and these quantities are relevant to asses thepossibilities of detecting darkonium in hadron colliders, direct or indirect detection experiments. We calculate | ψ (0) | and | ψ (cid:48) (0) | for l = 0 and l = 1 respectively, for the lowest energy levels and in the whole range of values of thescreening lengths, using the Pad´e approximants technique.We wrote a Mathematica code for the calculation of eigenvalues and eigenstates of the Yukawa potential to thedesired order, available upon request.0
IX. APPENDIX
The polynomials entering the solution of the Yukawa potential to order k = 5 for n = 1 , , , N ( ρ, δ,
5) = M , ( ρ , δ ) = 1 − (cid:18) − ρ (cid:19) δ + (cid:18) − ρ − ρ (cid:19) δ − (cid:18) − ρ − ρ − ρ (cid:19) δ + (cid:18) − ρ − ρ − ρ − ρ (cid:19) δ , (124) N ( ρ, δ,
5) = M , ( ρ, δ ) = 2 − ρ + (cid:18) −
24 + 12 ρ + 2 ρ − ρ (cid:19) δ + (cid:18) − ρ − ρ ρ ρ (cid:19) δ + (cid:18) − ρ − ρ ρ
36 + 19 ρ − ρ (cid:19) δ + (cid:18) − ρ ρ − ρ − ρ
15 + 106 ρ
225 + 53 ρ (cid:19) δ , (125) N ( ρ, δ,
5) = M , ( ρ, δ ) = 1 + (cid:18) −
15 + ρ (cid:19) δ + (cid:18) − ρ − ρ (cid:19) δ + (cid:18) − − ρ ρ
36 + 7 ρ (cid:19) δ + (cid:18) ρ − ρ − ρ − ρ (cid:19) δ , (126) N ( ρ, δ,
5) = M , ( 3 ρ , δ ) = 3 − ρ + ρ (cid:18) − ρ − ρ − ρ
16 + 27 ρ (cid:19) δ + (cid:18) − ρ ρ ρ − ρ − ρ (cid:19) δ + (cid:18) − ρ − ρ
128 + 20169 ρ − ρ − ρ
32 + 405 ρ (cid:19) δ + (cid:18) − ρ
320 + 60293403 ρ − ρ
128 + 400221 ρ ρ − ρ − ρ (cid:19) δ , (127) N ( ρ, δ,
5) = M , ( 3 ρ , δ ) = 4 − ρ + (cid:18) −
270 + 135 ρ ρ − ρ (cid:19) δ + (cid:18) − ρ − ρ + 9 ρ ρ (cid:19) δ + (cid:18) − ρ − ρ + 567 ρ ρ − ρ (cid:19) δ + (cid:18) − ρ
160 + 38637 ρ − ρ − ρ
320 + 105381 ρ ρ (cid:19) δ , (128) N ( ρ, δ,
5) = M , ( 3 ρ , δ ) = 1 + (cid:18) − ρ (cid:19) δ + (cid:18) − ρ − ρ (cid:19) δ + (cid:18) − − ρ ρ
128 + 405 ρ (cid:19) δ + (cid:18) ρ − ρ − ρ − ρ (cid:19) δ , (129) N ( ρ, δ,
5) = M , (2 ρ, δ ) = 4 − ρ + 2 ρ − ρ (cid:18) −
768 + 1152 ρ − ρ − ρ + 12 ρ − ρ (cid:19) δ + (cid:18) − ρ + 18880 ρ − ρ − ρ ρ ρ (cid:19) δ + (cid:18) − ρ − ρ ρ − ρ − ρ ρ − ρ (cid:19) δ + (cid:18) − ρ ρ − ρ
15 + 107712 ρ ρ − ρ ρ
75 + 824 ρ (cid:19) δ , (130)1 N ( ρ, δ,
5) = M , (2 ρ, δ ) = 10 − ρ + ρ (cid:0) − ρ − ρ − ρ + 2 ρ (cid:1) δ + (cid:18) − ρ ρ ρ − ρ − ρ (cid:19) δ + (cid:18) − ρ − ρ ρ − ρ − ρ ρ (cid:19) δ + (cid:18) − ρ + 24671104 ρ − ρ ρ ρ − ρ − ρ (cid:19) δ , (131) N ( ρ, δ,
5) = M , (2 ρ, δ ) = 6 − ρ + (cid:0) − ρ + 36 ρ − ρ (cid:1) δ + (cid:18) − ρ − ρ + 16 ρ + 16 ρ (cid:19) δ + (cid:18) − ρ − ρ + 736 ρ + 620 ρ − ρ (cid:19) δ + (cid:18) − ρ ρ − ρ − ρ ρ
25 + 1648 ρ (cid:19) δ , (132) N ( ρ, δ,
5) = M , (2 ρ, δ ) = 1 + (cid:0) −
360 + 4 ρ (cid:1) δ + (cid:18) − ρ − ρ (cid:19) δ + (cid:18) − − ρ ρ ρ (cid:19) δ + (cid:18) ρ − ρ − ρ − ρ (cid:19) δ . (133) [1] H. Yukawa, On the Interaction of Elementary Particles I , Proc. Phys. Math. Soc. Jap. (1935) 48.[2] V. Debye and E. Huckel, Zur Theorie der Elektrolyte , Physikalische Zeitschrift (1923) 185.[3] H. Margenau and M. Lewis, Structure of spectral lines from plasmas , Rev. Mod. Phys. (1959) 569.[4] G. M. Harris, Attractive two-body interactions in partially ionized plasmas , Phys. Rev. (1962) 1131.[5] C. R. Smith,
Bound states in a debye-h¨uckel potential , Phys. Rev. (1964) A1235.[6] C. S. Lam and Y. P. Varshni,
Ionization energy of the helium atom in a plasma , Phys. Rev. A (1983) 418.[7] H. M. Schey and J. L. Schwartz, Counting the bound states in short-range central potentials , Phys. Rev. (1965)B1428.[8] K. M. Roussel and R. F. O’Connell,
Variational solution of schr¨odinger’s equation for the static screened coulombpotential , Phys. Rev. A (1974) 52.[9] J. B. Krieger, Electron shielding in heavily doped semiconductors , Phys. Rev. (1969) 1337.[10] B. Zee,
Models and method of calculation of doping and injection-dependent impurity density of states in gaas , Phys. Rev.B (1979) 3167.[11] A. Ferraz, N. March and F. Flores, Metal-insulator transition in hydrogen and in expanded alkali metals , Journal ofPhysics and Chemistry of Solids (1984) 627 .[12] V. Bargmann, On the number of bound states in a central field of force , Proceedings of the National Academy of Sciences (1952) 961 [ ].[13] G. J. Iafrate, Critical screening parameters for the generalized screened coulomb potential , Phys. Rev. A (1973) 1138.[14] C. Stubbins, Bound states of the hulth´en and yukawa potentials , Phys. Rev. A (1993) 220.[15] O. A. Gomes, H. Chacham and J. R. Mohallem, Variational calculations for the bound-unbound transition of the yukawapotential , Phys. Rev. A (1994) 228.[16] J. McEnnan, L. Kissel and R. H. Pratt, Analytic perturbation theory for screened coulomb potentials: Nonrelativistic case , Phys. Rev. A (1976) 532.[17] J. P. Edwards, U. Gerber, C. Schubert, M. A. Trejo and A. Weber, The Yukawa potential: ground state energy andcritical screening , PTEP (2017) 083A01 [ ].[18] E. R. Vrscay,
Hydrogen atom with a yukawa potential: Perturbation theory and continued-fractions–pad´e approximants atlarge order , Phys. Rev. A (1986) 1433.[19] C. S. Lam and Y. P. Varshni, Energies of s eigenstates in a static screened coulomb potential , Phys. Rev. A (1971) 1875.[20] R. Dutt, K. Chowdhury and Y. P. Varshni, An improved calculation for screened coulomb potentials inrayleigh-schrodinger perturbation theory , Journal of Physics A: Mathematical and General (1985) 1379.[21] V. Eletsky, V. Popov and V. Weinberg, Logarithmic Perturbation Theory for the Screened Coulomb Potential , Phys. Lett.A (1981) 235.[22] V. Vainberg, V. Eletskii and V. Popov, Logarithmic Perturbation Theory for a Screened Coulomb Potential and a Charmonium Potential , Sov. Phys. JETP (1981) 833.[23] A. E. S. Green, Energy eigenvalues for yukawa potentials , Phys. Rev. A (1982) 1759.[24] C. S. Lai, Pad´e approximants and perturbation theory for screened coulomb potentials , Phys. Rev. A (1981) 455.[25] S. L. Garavelli and F. A. Oliveira, Analytical solution for a yukawa-type potential , Phys. Rev. Lett. (1991) 1310.[26] C. S. Lai and B. Suen, Alternative approach to perturbation theory for screened coulomb potentials , Phys. Rev. A (1980) 1100.[27] G. Moreno and A. Zepeda, N Expansion for a Yukawa Potential , J. Phys. B (1984) 21.[28] B. Gonul, K. Koksal and E. Bakir, An alternative treatment for yukawa-type potentials , Physica Scripta (2006) 279.[29] S. H. Patil, Energy levels of screened coulomb and hulthen potentials , Journal of Physics A: Mathematical and General (1984) 575.[30] H. E. Montgomery, K. D. Sen and J. Katriel, Critical screening in the one- and two-electron yukawa atoms , Phys. Rev. A (2018) 022503.[31] F. J. Rogers, H. C. Graboske and D. J. Harwood, Bound eigenstates of the static screened coulomb potential , Phys. Rev.A (1970) 1577.[32] C. A. Rouse, Screened coulomb solutions of the schr¨odinger equation , Phys. Rev. (1967) 41.[33] C. G. Diaz, F. M. Fernandez and E. A. Castro,
Critical screening parameters for screened coulomb potentials , Journal ofPhysics A: Mathematical and General (1991) 2061.[34] Y. Li, X. Luo and H. Kroger, Bound states and critical behavior of the yukawa potential , Science in China Series G (2006) 60.[35] J. F. Navarro, A. Ludlow, V. Springel, J. Wang, M. Vogelsberger, S. D. White et al., The Diversity and Similarity ofCold Dark Matter Halos , Mon. Not. Roy. Astron. Soc. (2010) 21 [ ].[36] W. de Blok,
The Core-Cusp Problem , Adv. Astron. (2010) 789293 [ ].[37] D. N. Spergel and P. J. Steinhardt,
Observational evidence for selfinteracting cold dark matter , Phys. Rev. Lett. (2000)3760 [ astro-ph/9909386 ].[38] N. Yoshida, V. Springel, S. D. White and G. Tormen, Collisional dark matter and the structure of dark halos , Astrophys.J. Lett. (2000) L103 [ astro-ph/0002362 ].[39] J. Miralda-Escude,
A test of the collisional dark matter hypothesis from cluster lensing , Astrophys. J. (2002) 60[ astro-ph/0002050 ].[40] A. Loeb and N. Weiner,
Cores in Dwarf Galaxies from Dark Matter with a Yukawa Potential , Phys. Rev. Lett. (2011) 171302 [ ].[41] M. H. Chan,
Observational evidences of the Yukawa Potential Interacting Dark Matter , Astrophys. J. Lett. (2013) L2[ ].[42] S. Khrapak, A. Ivlev, G. Morfill and S. Zhdanov,
Scattering in the Attractive Yukawa Potential in the Limit of StrongInteraction , Phys. Rev. Lett. (2003) 225002.[43] H. Hern´andez-Arellano, M. Napsuciale and S. Rodr´ıguez, Spin portal to dark matter , Phys. Rev.
D98 (2018) 015001[ ].[44] H. Hern´andez-Arellano, M. Napsuciale and S. Rodr´ıguez,
Spin-one dark matter and gamma ray signals from the galacticcenter , JHEP (2020) 106 [ ].[45] M. Napsuciale, S. Rodr´ıguez and H. Hern´andez-Arellano, Kinetic mixing, custodial symmetry and a lower bound on thedark Z (cid:48) mass , .[46] K. Babu, C. F. Kolda and J. March-Russell, Implications of generalized Z - Z-prime mixing , Phys. Rev. D (1998) 6788[ hep-ph/9710441 ].[47] B. Holdom, Two U(1)’s and Epsilon Charge Shifts , Phys. Lett. B (1986) 196.[48] J. L. Hewett and T. G. Rizzo,
Low-Energy Phenomenology of Superstring Inspired E(6) Models , Phys. Rept. (1989)193.[49] K. R. Dienes, C. F. Kolda and J. March-Russell,
Kinetic mixing and the supersymmetric gauge hierarchy , Nucl. Phys. B (1997) 104 [ hep-ph/9610479 ].[50] P. Langacker,
The Physics of Heavy Z (cid:48) Gauge Bosons , Rev. Mod. Phys. (2009) 1199 [ ].[51] ATLAS collaboration,
Search for dark matter and other new phenomena in events with an energetic jet and large missingtransverse momentum using the ATLAS detector , JHEP (2018) 126 [ ].[52] CMS collaboration,
Search for dark matter produced with an energetic jet or a hadronically decaying W or Z boson at √ s = 13 TeV , JHEP (2017) 014 [ ].[53] J. D. March-Russell and S. M. West, WIMPonium and Boost Factors for Indirect Dark Matter Detection , Phys. Lett. B (2009) 133 [ ].[54] H. An, M. B. Wise and Y. Zhang,
Effects of Bound States on Dark Matter Annihilation , Phys. Rev. D (2016) 115020[ ].[55] M. Cirelli, P. Panci, K. Petraki, F. Sala and M. Taoso, Dark Matter’s secret liaisons: phenomenology of a dark U(1)sector with bound states , JCAP (2017) 036 [ ].[56] K. Petraki, M. Postma and J. de Vries, Radiative bound-state-formation cross-sections for dark matter interacting via aYukawa potential , JHEP (2017) 077 [ ].[57] A. Krovi, I. Low and Y. Zhang, Broadening Dark Matter Searches at the LHC: Mono-X versus Darkonium Channels , JHEP (2018) 026 [ ].[58] J. Harz and K. Petraki, Higgs Enhancement for the Dark Matter Relic Density , Phys. Rev. D (2018) 075041[ ]. [59] M. Napsuciale and S. Rodriguez, Complete analytical solution to the quantum Yukawa potential , .[60] L. Infeld and T. Hull, The factorization method , Rev. Mod. Phys. (1951) 21.[61] L. Gendenshtein, Derivation of Exact Spectra of the Schrodinger Equation by Means of Supersymmetry , JETP Lett. (1983) 356.[62] E. Witten, Dynamical Breaking of Supersymmetry , Nucl. Phys. B (1981) 513.[63] F. Cooper and B. Freedman,
Aspects of supersymmetric quantum mechanics , Annals of Physics (1983) 262 .[64] D. J. Fernandez,
New hydrogen - like potentials , Lett. Math. Phys. (1984) 337 [ quant-ph/0006119 ].[65] C. M. Bender and T. T. Wu, Anharmonic oscillator , Phys. Rev. (1969) 1231.[66] C. M. Bender and T. Wu,
Anharmonic oscillator. 2: A Study of perturbation theory in large order , Phys. Rev. D (1973)1620.[67] J. Zinn-Justin, Perturbation Series at Large Orders in Quantum Mechanics and Field Theories: Application to theProblem of Resummation , Phys. Rept. (1981) 109.[68] P. M. Stevenson, Optimized perturbation theory , Phys. Rev. D (1981) 2916.[69] A. Okopinska, Nonstandard expansion techniques for the effective potential in lambda phi**4 quantum field theory , Phys.Rev. D (1987) 1835.[70] G. Arteca, F. Fernandez and E. Castro, Large order perturbation theory and summation methods in quantum mechanics .Springer-Verlag, Berlin, 1990.[71] G. A. Baker, Jr,
The theory and application of the pade approximant method . , pp 1-58 of Advances in TheoreticalPhysics. Vol. I. Brueckner, Keith A. (ed.). New York, Academic Press, 1965. .[72] B. Guberina, J. H. Kuhn, R. D. Peccei and R. Ruckl, Rare Decays of the Z0 , Nucl. Phys. B174