A Neumann series of Bessel functions representation for solutions of the radial Dirac system
Vladislav V. Kravchenko, Elina L. Shishkina, Sergii M. Torba
aa r X i v : . [ m a t h - ph ] A ug A Neumann series of Bessel functions representation for solutions ofthe radial Dirac system
Vladislav V. Kravchenko , Elina L. Shishkina and Sergii M. Torba Departamento de Matem´aticas, CINVESTAV del IPN, Unidad Quer´etaro,Libramiento Norponiente Voronezh State University.e-mail: [email protected], [email protected], ∗ August 13, 2020
Abstract
A new representation for a regular solution of the radial Dirac system of a special formis obtained. The solution is represented as a Neumann series of Bessel functions uniformlyconvergent with respect to the spectral parameter. For the coefficients of the series convenientfor numerical computation recurrent integration formulas are given. Numerical examples arepresented.
We consider the one-dimensional radial Dirac system of the form ω − ddr + κr − p ( r ) ddr + κr − p ( r ) ω g ( r ) f ( r ) = 0 , (1)where p is absolutely continuous complex-valued function on some interval [0 , b ], ω , ω ∈ C , κ isthe spin-orbit quantum number, g and f are lower and upper radial wave functions, respectively.The system (1) with ω = mc + V s ( r )+ E − V v ( r ) ~ c , ω = mc + V s ( r ) − E + V v ( r ) ~ c and p ( r ) = V ps ( r ) ~ c arisesin quantum mechanics when studying the radial Dirac equation. Here V s is a scalar potential, V v is the time component of a vector potential and V ps is a pseudoscalar or tensor potential, see, forexample, formula 2.1 in [1], formulas (21)–(22) in [6], system (13)–(14) in [20] and (1) in [19]. Thesystem (1) is a special case of the radial Dirac equation in the presence of a tensor or a pseudoscalarpotential, and when both scalar and vector potentials are constant. System (1) appears in the recentJackiw-Pi model of the bilayer graphene [8], [10]. There is a considerable number of publicationsin which Dirac-type equations (1) are examined, but mostly either exactly solvable potentials aresought (see, e.g., [1], [6], [20]), or an approximate solution is constructed for a concrete potential(see, for example, [7]).In the present work for an arbitrary potential p ( r ) we obtain an analytical representation fora regular solution of (1) in the form of a functional series with a simple recurrent integration ∗ Research was supported by CONACYT, Mexico via the projects 222478 and 284470. Research of VladislavKravchenko was supported by the Regional mathematical center of the Southern Federal University, Russia. ω = ω , the estimates are independent of their values, while in the case ω = ω the estimates involve the factor (cid:12)(cid:12)(cid:12)p ω /ω (cid:12)(cid:12)(cid:12) , and thus depends on how much the couplingconstants differ from each other.The NSBF representations for solutions of Sturm-Liouville type equations proved to be usefulfor solving both direct and inverse spectral problems [4], [9], [11], [12], [13], [14], [15], [16], [17],[18]. In [13] an NSBF representation was obtained for solutions of the one-dimensional stationarySchr¨odinger equation. In [16] that result was generalized onto the case of an arbitrary regular Sturm-Liouville equation. Recently in [17] an NSBF representation was obtained for regular solutions ofperturbed Bessel equations. In [4], [9], [11], [12], [15] NSBF representations for solutions were usedfor solving inverse spectral problems.In the present paper an NSBF representation for regular solutions of (1) is obtained by trans-forming the system into a couple of perturbed Bessel equations and using results from [17]. Weprove the above mentioned error estimates for partial sums of the series representations and discussthe numerical implementation of the NSBF representation. We show that the spectral parameterindependent error estimates are evident, indeed, in numerical experiments and show the applica-bility of the obtained NSBF representation for solving spectral problems for (1).The paper is organized as follows. In Section 2 we obtain the NSBF representation for theregular solution of (1) and prove a convergence result for the approximate solution. In Section 3 wesummarize the steps required for numerical solution of equation (1) and related spectral problemsusing the proposed representation and show numerical results for the Dirac oscillator. Consider the following two component radial Dirac system (cid:18) ddr − κr + p ( r ) (cid:19) f = ω g, (2) (cid:18) ddr + κr − p ( r ) (cid:19) g = − ω f, (3)where ω , ω ∈ C , κ ≥ , and p ∈ AC[0 , b ] is in general a complex valued function.
Definition 1
A pair of functions ( f κ , g κ ) is called a regular solution of the system (2) – (3) if itsatisfies the system as well as the following asymptotic conditions f κ ( r ) ∼ C f r κ , g κ ( r ) ∼ C g r κ +1 , when r → where C f and C g are some constants. Together with the potential p the following functions will be considered q ( r ) = p ′ ( r ) − κr p ( r ) + p ( r ) and q ( r ) = − p ′ ( r ) − κr p ( r ) + p ( r ) . (4)2ote that if ( f κ , g κ ) is a regular solution of (2)–(3), the functions f κ and g κ are necessarily regularsolutions of the equations − f ′′ + (cid:18) κ ( κ − r + q ( r ) (cid:19) f = ω f, r ∈ (0 , b ] (5)and − g ′′ + (cid:18) κ ( κ + 1) r + q ( r ) (cid:19) g = ω g, r ∈ (0 , b ] , (6)respectively with ω = ω ω .Note that for p ∈ AC[0 , b ] both potentials q and q are such that r ε q , ( r ) ∈ L (0 , b ) for anysmall ε >
0, hence the conditions on the potential from [17] are satisfied. In order to apply theresults of [17] to equations (5) and (6) we need two solutions of the equations − f ′′ + (cid:18) κ ( κ − r + q ( r ) (cid:19) f = 0 and − g ′′ + (cid:18) κ ( κ + 1) r + q ( r ) (cid:19) g = 0 , (7)non-vanishing on (0 , b ] and satisfying the following asymptotics at zero f ( r ) ∼ r κ and g ( r ) ∼ r κ +1 , when r → . (8)The solution f can be directly obtained by taking ω = 0 in (2) and is given by f ( r ) = r κ exp (cid:18) − Z r p ( s ) ds (cid:19) . (9)To obtain the solution g note that the function 1 /f is a solution of (3) with ω = 0, and henceit is the solution of the second equation in (7) satisfying the asymptotic relation 1 /f ( r ) ∼ r − κ , r →
0. A second linearly independent solution of (3) with ω = 0 can be chosen in the form Cf ( r ) R r f ( s ) ds . Chosing C = 2 κ + 1, i.e., taking g ( r ) = (2 κ + 1) r − κ exp (cid:18)Z r p ( s ) ds (cid:19) Z r t κ exp (cid:18) − Z t p ( s ) ds (cid:19) dt (10)we obtain the solution of the second equation in (7) satisfying (8).It can be seen from (9) and (10) that the derivatives of the solutions f and g are given by f ′ ( r ) = (cid:16) κr − p ( r ) (cid:17) f ( r ) and g ′ ( r ) = (cid:16) p − κr (cid:17) g ( r ) + (2 κ + 1) f ( r ) . (11)The solution f given by (9) is always non-vanishing on (0 , b ]. The solution g given by (10)is definitely non-vanishing for real valued potentials p but may possess zeros for complex valuedfunctions p . For this reason the following assumption (A) concerning the potential p will be madethroughout this paper. We assume that the second equation in (7) admits a regular solution g which does not vanish on (0 , b ]. This assumption does not imply any additional restriction on p forthe following reason. In [17, Proposition B.1] we show that one can always choose such a constant c that the second equation in (7) with the potential e q ( x ) := q ( x ) + c possesses a non-vanishingsolution. Equation (6) can then be written as − g ′′ + (cid:16) κ ( κ +1) r + e q ( r ) (cid:17) g = e ω g with e ω = ω + c which leads to the same results and conclusions as below. Also, one may construct the regularsolution of the system (2)–(3) using only the solution f and its derivative (see Remark 5), howeverloosing an attractive possibility to verify the accuracy of approximate solutions (see Remark 4 andSubsection 3.1).Thus, without loss of generality we assume that the regular solution g of the second equationin (7) satisfying (8) does not have zeros in (0 , b ].3 heorem 2 Let p ∈ AC[0 , b ] , and the assumption (A) be fulfilled. Then a regular solution of thesystem (2) – (3) satisfying the asymptotic relations (here ω = ω ω ) f κ ( r ) ∼ − ω κ +1 ω d ( κ − r κ and g κ ( r ) ∼ ω κ +1 d ( κ ) r κ +1 , r → , has the form f κ ( r ) = − ω ω rj κ − ( ωr ) − ωω ∞ X n =0 β ,n ( r ) j κ +2 n ( ωr ) , (12) g κ ( r ) = ωrj κ ( ωr ) + ∞ X n =0 β ,n ( r ) j κ +2 n +1 ( ωr ) , (13) where j ν ( r ) = p π r J ν + ( r ) is the spherical Bessel function of the first kind, d ( κ ) := √ π κ +1 Γ( κ + 3 / . Denote u := g and u := f , where f and g are solutions of (7) satisfying (8) . Then thefunctions β j,n , j ∈ { , } , n ≥ , can be found from the recurrent formulas β j, ( r ) = (2 κ − j + 5) (cid:18) u j ( r ) r κ +2 − j − (cid:19) , j ∈ { , } , (14) β j,n ( r ) = − n + 2 κ − j + 54 n + 2 κ − j + 1 (cid:20) β j,n − ( r ) + 2(4 n + 2 κ − j + 3) u j ( r ) θ j,n ( r ) r n + κ − j +2 (cid:21) , (15) θ j,n ( r ) = Z r η j,n ( t ) − t n + κ − j +1 β j,n − ( t ) u j ( t ) u j ( t ) dt, (16) η j,n ( r ) = Z r (cid:2) tu ′ j ( t ) + (2 n + κ − j + 1) u j ( t ) (cid:3) t n + κ − j β j,n − ( t ) dt. (17) Proof.
From (3) we have that f κ = − ω (cid:0) g ′ κ + κg κ /r − pg κ (cid:1) . Hence if g κ ( r ) ∼ d ( κ )( ωr ) κ +1 , g ′ κ ( r ) ∼ ( κ +1) ωd ( κ )( ωr ) κ when r →
0, then f κ ( r ) ∼ − ωω (2 κ + 1) d ( κ )( ωr ) κ .Now we apply Theorem 5.2 from [17] in order to find out that a solution g κ of (6) satisfying therelation g κ ( r ) ∼ d ( κ )( ωr ) κ +1 , has the form (13), meanwhile a solution e f κ of (5) satisfying therelation e f κ ( r ) ∼ d ( κ − ωr ) κ , when r →
0, can be written as e f κ ( r ) = ωrj κ − ( ωr ) + ∞ X n =0 β ,n ( r ) j κ +2 n ( ωr ) . And f κ = − ωω (2 κ + 1) d ( κ ) d ( κ − e f κ = − ωω e f κ , which leads to (12).For practical use of the representation (12), (13) the estimates of the difference between theexact solution and its approximation defined as f κ,N ( r ) = − − ω ω rj κ − ( ωr ) − ωω N X n =0 β ,n ( r ) j κ +2 n ( ωr ) , (18) g κ,N ( r ) = ωrj κ ( ωr ) + N X n =0 β ,n ( r ) j κ +2 n +1 ( ωr ) (19)are needed. From Theorem 2 using [17, Theorem 5.2] the following result follows immediately.4 roposition 3 Under the conditions of Theorem 2 the following inequalities are valid | g κ ( r ) − g κ,N ( r ) | ≤ √ rε N ( r ) and | f κ ( r ) − f κ,N ( r ) | ≤ (cid:12)(cid:12)(cid:12)(cid:12) ωω (cid:12)(cid:12)(cid:12)(cid:12) √ rε N ( r ) for all ω ∈ R , ω ∈ R \ { } , where ε N is a nonnegative function independent on ω and ω , suchthat max r ∈ [0 ,b ] ε N ( r ) → as N → ∞ . Similar result holds for ω belonging to a strip | Im ω | ≤ C ,with addition of a multiplicative constant dependent only on the value of C .Suppose additionally that p ∈ W k [0 , b ] , p (0) = 0 and p ( r ) r ∈ W k − [0 , b ] for some k ∈ N . Here W k [0 , b ] denotes the class of functions having k derivatives, the last one belonging to L [0 , b ] space,and p ( r ) /r is assumed to have a finite limit as r → . Then there exists a constant c , such that ε N ( r ) ≤ cN k , N > κ + k + 1 . The independence of ε N of ω implies that the approximate solution ( f κ,N , g κ,N ) remains goodeven for very large values of | Re ω | . Remark 4
Even though the derivatives of the regular solutions ( f κ , g κ ) are readily available from (2) and (3) , an independent representation (useful, e.g., for verification of accuracy of approximatesolutions) for them can be obtained from [17, Theorem 6.3]. Under the conditions and notations ofTheorem 2 let Q j ( r ) := R r q j ( t ) dt , j ∈ { , } and let the functions γ j,n , j ∈ { , } , n ≥ be definedas γ j, ( r ) = (2 κ − j + 5) (cid:18) u ′ j ( r ) r κ − j +2 − κ − j + 2 r − Q j ( r )2 (cid:19) ,γ j,n ( r ) = − n + 2 κ − j + 54 n + 2 κ − j + 1 (cid:20) γ j,n − ( r )+ (4 n + 2 κ − j + 3) (cid:18) u ′ j ( r ) θ j,n ( r ) r n + κ − j +2 + 2 η j,n ( r ) u j ( r ) r n + κ − j +2 − β j,n − ( r ) r (cid:19)(cid:21) . Then f ′ κ ( r ) = − ω ω rj κ − ( ωr ) − (cid:18) rQ ( r )2 − κ + 1 (cid:19) ω ω j κ − ( ωr ) − ωω ∞ X n =0 γ ,n ( r ) j n + κ ( ωr ) , (20) g ′ κ ( r ) = ω rj κ − ( ωr ) + (cid:18) rQ ( r )2 − κ (cid:19) ωj κ ( ωr ) + ∞ X n =0 γ ,n ( r ) j n + κ +1 ( ωr ) . (21) Remark 5
The regular solution of the system (2) – (3) can be obtained using only the particularsolution f and related functions β ,n and γ ,n , n ≥ , without the need of the functions g , β ,n and γ ,n at all. Indeed, g κ = 1 ω (cid:16) f ′ κ − κr f κ + p ( r ) f κ (cid:17) and g ′ κ = − ω f κ − κr g κ + p ( r ) g κ , and the representations for f κ and f ′ κ are given by (12) and (20) . A numerical method based on the representation (12)–(13) of the regular solution of the system(2)–(3) consists in the following steps. 5. Compute a pair ( f , g ) of regular solutions of (7) satisfying (8) using (9) and (10). Computealso their derivatives ( f ′ , g ′ ). In the case that the coefficient p is complex valued, check if theassumption (A) holds, and if not, proceed as described in Remark 5 or look for a spectral shift(see Appendix B in [17, (8.1)]) such that a pair of solutions ( f , g ) becomes non-vanishing.2. Compute the coefficients β j,n , j ∈ { , } , n ∈ { , , . . . , N } using the formulas (14)–(17).Note that the coefficients β j,n satisfy [17] ∞ X n =0 ( − n β j,n ( r ) = rQ j ( r )2 , r ∈ [0 , b ] , j ∈ { , } (22)and decay to zero (however, not necessary monotonously) as n → ∞ . The equality (22) canbe used to estimate an optimal number of the coefficients N , as a value where the truncatedsums cease to decrease when N increases.3. Compute approximate solutions f κ,N and g κ,N using (12) and (13).4. The accuracy of the obtained approximations can be estimated by calculating the discrepan-cies f ′ κ,N − κr f κ,N + p ( r ) f κ,N − ω g κ,N and g ′ κ,N + κr f κ,N − p ( r ) f κ,N + ω g κ,N , (23)where f ′ κ,N and g ′ κ,N are computed from the truncated series (20) and (21).We refer the reader to [13] and [18] for implementation details of the proposed algorithm. As a test example for the proposed algorithm we consider the Dirac oscillator [21, 3, 5].The large radial component F ( r ) and the small radial component G ( r ) of the Dirac wavefunction are solutions of the following system (cid:18) − ddr + (cid:18) ε ( j + 1 / r + mωr (cid:19)(cid:19) G ( r ) = ( E − m ) F ( r ) , (24) (cid:18) ddr + (cid:18) ε ( j + 1 / r + mωr (cid:19)(cid:19) F ( r ) = ( E + m ) G ( r ) , (25)where j is the total angular momentum quantum number, ε = ± m is the mass of the particleand ω is the frequency. Note that the number l := j + ε j is always equal to 1 / E − m = mω (cid:0) N + 1) + ε (2 j + 1) (cid:1) for the positive-energy states, and from E − m = mω (cid:0) N + 2) + ε (2 j + 1) (cid:1) G -15 -10 -5 Abs. error of F
Abs. error of G G -15 -10 -5 Abs. error of F
Abs. error of G G -15 -10 -5 Abs. error of F
Abs. error of G
Figure 1: On the left: components F n, and G n, of the eigenfunction of the Dirac oscillator for n ∈ { , , } with the parameters j = 5 / ε = − m = ω = 1. On the right: absoluteerrors of these components.for the negative-energy states. Here N = 2 n + l , n = 0 , , , . . . , is the principal quantum number.The corresponding eigenfunctions are given by F n,l ( r ) = A (cid:0) r √ mω (cid:1) l +1 exp( − mωr / L l +1 / n ( mωr ) , (26) G n,l − ε ( r ) = A (cid:0) r √ mω (cid:1) l +1 − ε exp( − mωr / L l − ε +1 / n + ε/ − / ( mωr ) , (27)where L sk ( x ) is an associated Laguerre polynomial.The system (24)–(25) is of the type considered in this paper. Since the potential of the problemis increasing, we approximated the semiaxis spectral problem (of finding the values of E for whichthe regular solution belongs to L (0 , ∞ )) by truncating the potential and considering the Dirichlet7oundary condition. For any non-trivial solution both f κ and g κ can not be equal to zero at onepoint, so we choose the function f κ and considered f κ ( B ) = 0as the boundary condition for the problem truncated onto [0 , B ] segment. We refer the reader to[22, Section 7.4] for additional details on the convergence of the eigenvalues of truncated problemsto the exact ones. ε = 1 ε = − -5 e (r)r Q (r)/100 -5 e (r)r Q (r)/100 Figure 2: Application of formula (22) to determine optimal truncation interval for the Dirac oscil-lator with the parameters j = 5 / m = ω = 1. On both plots black dashed line shows the valueof | rQ ( r ) / | , and the blue solid line shows e ( r ) := min N ( r ) ≤ (cid:12)(cid:12)P N ( r ) n =0 ( − n β ,n ( r ) − rQ ( r ) / (cid:12)(cid:12) for computed coefficients β ,n .All the computations were performed in machine precision using Matlab 2017. We refer thereader to [18] for the details of the numerical realization. We considered two sets of parameters,having ε = ± j = 5 / m = ω = 1. For ε = 1 the corresponding potential is p = − mωr in the notations of (2), (3), and for ε = − p = mωr . In thefirst case the corresponding particular solution f given by (9) is rapidly increasing, for the secondcase f is rapidly decreasing. We decided to not implement interval subdivision techniques, andutilize the proposed representation directly to illustrate that even straightforward implementationcan deliver highly accurate results.First, we compare approximate solutions with the exact ones for the case ε = − E − m ∈ { , , } , corresponding to n ∈ { , , } . In terms of the system (2),(3) we have taken ω = 2, ω ∈ { , , } . On Figure 1 we present the solutions and correspondingabsolute errors. As one can observe, the error does not increase for large values of ω (correspondingto higher index eigenfunctions) and only increases for large values of r due to machine precisionlimitations. The approach presented in Remark 5 delivered a more accurate solution component G .This is due to the error near r = 0 in the particular solution g computed by (10). For that reasonon the plots we present the absolute errors obtained with the aid of the formula from Remark 5.Approximate solution of the spectral problem requires truncating the interval. A larger intervalallows one to compute more eigenvalues and more accurately. However this leads to larger errorsin all the coefficients β j,n computed, due to machine precision limitations. The equality (22) canbe utilized to estimate automatically a truncation parameter B . We took the segment [0 , ,
20] points and computed100 coefficients β ,n . After that we checked for each r the convergence of partial sums in (22)to rQ ( r ) /
2. Due to machine precision limitations, the difference between P Nn =0 ( − n β ,n ( r ) and rQ ( r ) / N ( r ), meaning that the difference essentiallydoes not decrease anymore when N increases. Let e ( r ) := (cid:12)(cid:12)P N ( r ) n =0 ( − n β ,n ( r ) − rQ ( r ) / (cid:12)(cid:12) . Wechose as the truncation parameter B the value 0 . · r , where r is such that for all r < r thevalue e ( r ) is small in comparison with rQ ( r ) (to be more precise, e ( r ) < | rQ ( r ) | / r > r the error e ( r ) can be larger than | rQ ( r ) | / B = 7 . ε = 1, and B = 9 . ε = −
1. See Figure 2 illustrating this procedure.In Table 1 we present approximate eigenvalues E − m for the parameters ε = ± j = 5 / m = ω = 1 computed on the truncated intervals [0 , B ]. ε = 1, on [0 , . ε = −
1, on [0 , . E − m Approximate14 13.99999999998718 17.999999999818322 21.999999998282826 25.999999964287130 29.999999427668234 33.999994269405738 37.999961665356442 42.000100568104446 46.00004836671550 50.012533027532354 54.037836743132658 58.2112436119225Number of β ,n used 24 Exact E − m Approximate0 7 . · − β ,n used 29Table 1: The eigenvalues for the Dirac oscillator problem (24), (25) truncated onto the segment[0 , B ]. Parameters used: ε = ± j = 5 / m = ω = 1. The last line shows the number of termsused in approximate solution (18). References [1] A. D. Alhaidari,
Solution of the Dirac equation for potential interaction , Int. J. Mod. Phys. A18 (2003), No. 27, 4955–4973.[2] A. Baricz, D. Jankov, T. K. Pog´any,
Series of Bessel and Kummer-type functions.
LectureNotes in Mathematics, 2207. Springer, Cham, 2017.[3] J. Ben´ıtez, R. P. Mart´ınez y Romero, H. N. N´u˜nez-Y´epez and A. L. Salas-Brito,
Solution andhidden supersymmetry of a Dirac oscillator , Phys. Rev. Lett. 64 (1990), no. 14, 1643–1645.[4] B. B. Delgado, K. V. Khmelnytskaya and V. V. Kravchenko,
The transmutation operatormethod for efficient solution of the inverse Sturm-Liouville problem on a half-line , Math. Meth.Appl. Sci. 42 (2019), 7359–7366.[5] F. Dom´ınguez-Adame and M. A. Gonz´alez,
Solvable linear potentials in the Dirac equation ,Europhys. Lett. 13 (1990), no. 3, 193–198. 96] M. Eshghi and H. Mehraban,
Eigen spectra in the Dirac-hyperbolic problem with tensor cou-pling , Chin. J. Phys. 50 (2012), No. 4, 533–543.[7] A. N. Ikot, H. Hassanabadi, E. Maghsoodi and V. Zarrinkamar,
Approximate solutions ofDirac equation for Tietz and general Manning-Rosen potentials using SUSYQM , Letters toElementary Particles and Atomic Nuclei, 11 (2014), No. 4(188), 673–687.[8] R. Jackiw and S.-Y. Pi,
Persistence of zero modes in a gauged Dirac model for bilayer graphene ,Phys. Rev. B, 78 (2008), 132104, 3pp.[9] A. N. Karapetyants, K. V. Khmelnytskaya and V. V. Kravchenko,
A practical method forsolving the inverse quantum scattering problem on a half line , J. Phys.: Conf. Ser. 1540 (2020),012007, 7pp.[10] K. V. Khmelnytskaya and H. C. Rosu,
An amplitude-phase (Ermakov–Lewis) approach for theJackiw–Pi model of bilayer graphene , J. Phys. A: Math. Theor. 42 (2009) 042004 (11pp).[11] V. V. Kravchenko,
On a method for solving the inverse Sturm–Liouville problem , J. InverseIll-pose P. 27 (2019), 401–407.[12] V. V. Kravchenko,
Direct and inverse Sturm-Liouville problems: A method of solution ,Birkh¨auser, Cham, 2020.[13] V. V. Kravchenko, L. J. Navarro and S. M. Torba,
Representation of solutions to the one-dimensional Schr¨odinger equation in terms of Neumann series of Bessel functions , Appl. Math.Comput. 314 (2017), 173–192.[14] V. V. Kravchenko, E. L. Shishkina and S. M. Torba,
On a series representation for integralkernels of transmutation operators for perturbed Bessel equations , Math. Notes 104 (2018),552–570.[15] V. V. Kravchenko, E. L. Shishkina and S. M. Torba,
A transmutation operator method forsolving the inverse quantum scattering problem . Submitted. Available at arXiv:2007.13039.[16] V. V. Kravchenko and S. M. Torba,
A Neumann series of Bessel functions representation forsolutions of Sturm-Liouville equations , Calcolo 55 (2018), article 11, 23pp.[17] V. V. Kravchenko and S. M. Torba,
An improved Neumann series of Bessel functions represen-tation for solutions of perturbed Bessel equations . Submitted. Available at arXiv:2005.10403v3.[18] V. V. Kravchenko, S. M. Torba and R. Castillo-P´erez,
A Neumann series of Bessel functionsrepresentation for solutions of perturbed Bessel equations , Appl. Anal. 97 (2018), 677–704.[19] S. Linnaeus,
Phase-integral solution of the radial Dirac equation , J. Math. Phys. 51 (2010),032304, 13pp.[20] R. Lisboa, M. Malheiro, A. S. de Castro, P. Alberto and M. Fiolhais,
Pseudospin symmetryand the relativistic harmonic oscillator , Phys. Rev. C, 69 (2004) 024319, 19pp .[21] M. Moshinsky and A. Szczepaniak,
The Dirac oscillator , J. Phys. A: Math. Gen. 22 (1989)L817–L819.[22] J. D. Pryce,
Numerical solution of Sturm-Liouville problems , Oxford: Clarendon Press, 1993.[23] G. N. Watson,
A Treatise on the theory of Bessel functions, 2nd ed., reprinted , CambridgeUniversity Press, Cambridge, UK, 1996, vi+804 pp.[24] J. E. Wilkins,