Front. Phys. Frontiers in Physics Front. Phys. 2296-424X Frontiers Media S.A. 1393130 10.3389/fphy.2024.1393130 Physics Original Research On the calculation of irregular solutions of the Schrödinger equation for non-spherical potentials with applications to metallic alloys Zeller 10.3389/fphy.2024.1393130 Zeller Rudolf * Peter Grünberg Institut, Forschungszentrum Jülich GmbH, Jülich, Germany

Edited by: Jamal Berakdar, Martin Luther University of Halle-Wittenberg, Germany

Reviewed by: Wolfram Hergert, Martin Luther University of Halle-Wittenberg, Germany

Arthur Ernst, Johannes Kepler University of Linz, Austria

 *Correspondence: Rudolf Zeller, ru.zeller@fz-juelich.de 30 07 2024 2024 12 1393130 28 02 2024 16 05 2024 Copyright © 2024 Zeller. 2024 Zeller

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The irregular solutions of the stationary Schrödinger equation are important for the fundamental formal development of scattering theory. They are also necessary for the analytical properties of the Green function, which in practice can greatly speed up calculations. Nevertheless, they are seldom considered in numerical treatments because of their divergent behavior at origin. This divergence demands high numerical precision that is difficult to achieve, particularly for non-spherical potentials which lead to different divergence rates in the coupled angular momentum channels. Based on an unconventional treatment of boundary conditions, an integral-equation method is here developed which is capable of dealing with this problem. The available precision is illustrated by electron-density calculations for NiTi in its monoclinic B19’ structure.

 Schrödinger equation non-spherical potentials scattering theory irregular solutions boundary conditions integral-equation method Chebyshev solver 1 Introduction

The problem of the numerical solution of the stationary Schrödinger equation for a single particle has been investigated in very many scientific publications but rarely with a focus on irregular solutions. While, for spherical potentials, the separation of radial and angular variables simplifies the problem into the solution of one-dimensional radial Schrödinger equations, the situation is more complicated for non-spherical potentials. Here, the separation of variables leads to radial equations where different angular momentum components are coupled by the non-diagonal potential matrix elements. If, as usual, a cutoff is applied by restricting angular-momentum components to l ≤ l max , then l max + 1 independent second-order linear differential equations must be solved for spherical potentials, while a set of ( l m a x + 1 ) 2 coupled second-order linear differential equations must be solved for non-spherical potentials. This represents a significant complication, particularly for the irregular solutions, which diverge with different powers of the radial variable r as r − l − 1 .

It is the purpose of this paper to present an approach which is capable of treating the divergent behavior in a numerically efficient manner. The approach is based on the integral-equation method of Gonzales et al. [1], who obtained regular solutions of the radial Schrödinger by integrations using Clenshaw–Curtis quadratures [2]. The numerical solution of second-order differential equations by integral-equation methods was introduced by Greengard [3] and Greengard and Rokhlin [4], who observed that stable, high-order numerical methods exist for the solution of integral equations. For instance, evaluation of the integrals with Clenshaw–Curtis quadratures leads to spectral accuracy. Spectral accuracy means that the results converge with the inverse p th power of the number of mesh points for p -times continuously differentiable integrands and exponentially for p   → ∞ . In contrast to this, the accuracy of solutions of differential equations by finite difference methods, like the Numerov method used in [5], is limited by a small inverse power of the number of mesh points.

The paper is organized as follows. In Section 2, the mathematical approach is presented. First, the integral equations for the coupled radial equations are defined and their boundary conditions are discussed. Then, we explain how the numerical effort can be reduced by using auxiliary integral equations and how discretization at Chebyshev collocation points leads to systems of linear algebraic equations which can be solved by standard numerical techniques. In Section 3, two examples are numerically investigated: a constant potential, for which the results are compared to the analytical results derived from the expressions given in Appendix A, and a realistic potential as it appears in all-electron density-functional electronic-structure calculations. It is shown that accurate bound-state wavefunctions and energies are obtained for constant potentials and that straightforward complex-contour integrations for calculating the electron density from the Green function can be applied. For that purpose, the correct divergence of the Green function at the origin is enforced by using unconventional boundary conditions. The numerical investigations are done with the KKRnano code of the JuKKR code package. This code is based on the full-potential screened Korringa–Kohn–Rostoker Green function method [6] and was developed for density-functional calculations for systems with thousands of atoms [7].

 2 Mathematical approach 2.1 Coupled radial equations

The coupled regular and irregular solutions of the Schrödinger equation − ∇ r 2 + V r − E Ψ r ; E = 0 for the potential V r can be defined [8] by linear Fredholm integral equations of the second kind as R L ′ L r ; k = j l ′ k r δ L ′ L + ∫ 0 ∞ d r ′ r ′ 2 g l ′ r , r ′ ; k ∑ L ″ V L ′ L ″ r ′ R L ″ L r ′ ; k and as S L ′ L r ; k = − i k h l ′ 1 k r β L ′ L k + ∫ 0 ∞ d r ′ r ′ 2 g l ′ r , r ′ ; k ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k . Here V L L ′ r = ∫ d r ̂ Y L r ̂ V r Y L ′ r ̂ are matrix elements of the potential and β L ′ L k = δ L ′ L − ∫ 0 ∞ d r r 2 j l ′ k r ∑ L ″ V L ′ L ″ r S L ″ L r ; k is a matrix, which implicitly depends on the irregular solutions. The function g l r , r ′ ;   k is given by g l r , r ′ ; k = − i k j l k r h l 1 k r ′ for r ≤ r ′ h l 1 k r j l k r ′ for r ≥ r ′ . In these equations and throughout the paper, Rydberg atomic units are used. j l , h l ( 1 ) = j l + i n l , and n l are spherical Bessel, Hankel, and Neumann functions, Y L spherical harmonics, and L a combined index for the angular momentum quantum numbers l and m , respectively. Radial and angular variables are denoted by r = r and r ̂ = r / r , respectively, and k = E is the square root of the energy variable.

The important difference between the inhomogeneous integral Eq. 2 for the regular solutions and Eq. 3 for the irregular solutions is that the source term in Eq. 2 contains Bessel functions, which lead to the r l ′ behavior of the regular solutions R L ′ L r ;   k at the origin. The source term in Eq. 3 contains Hankel functions, which lead to the r − l ′ − 1 behavior of the irregular solutions S L ′ L r ;   k at the origin. The numerical solution of Eq. 3 demands high accuracy because the integrand contains functions which increase with different powers r − l ″ − 1  at the origin.

If, as is often done, the coupled radial solutions are not determined from the Fredholm integral Eqs 2, 3 but from differential equations, an additional difficulty arises. The differential equations can be obtained from Eqs 2, 3 by applying the operator L r = − d 2 d r 2 − 2 r d d r + l ′ l ′ + 1 r 2 − k 2 . With L r g l ′ r ′ , r ;   k = − δ r   −   r ′ / r ′ 2 , L r j l ′ k r = 0 , and L r h l ′ ( 1 ) ( k r ) = 0 , this leads to the coupled Schrödinger equations ∑ L ″ − d 2 d r 2 − 2 r d d r + l ′ l ′ + 1 r 2 − k 2 δ L ′ L ″ + V L ′ L ″ r R L ″ L r ; k = 0 for the regular solutions R L ″ L r ;   k and to an identical equation for the irregular solutions S L ″ L r ;   k . With the cutoff l ≤ l max , the differential equation Eq. 8 has 2 ( l m a x + 1 ) 2 linearly independent solutions—one regular and one irregular solution for each value of L . The different solutions are distinguished by different boundary conditions. These conditions must be specified explicitly for the differential equation (Eq. 8) while they are naturally contained in the integral equations (Eqs 2, 3) as a consequence of the source terms. During the numerical solution of the differential equation (Eq. 8), it is essential to maintain the linear independence of the solutions. Because of the discretization error, this represents a considerable challenge already for the regular solutions, for instance, as explained in [9], and an even greater challenge for the irregular solutions because these diverge at the origin. A discretization error, of course, also occurs in numerical treatments of integral equations, but by using Clenshaw–Curtis quadrature, the error can be made substantially smaller so that accurate results can be achieved.

 2.2 Boundary conditions

To discuss the boundary conditions, it is convenient to assume finite integration limits r min and r max . This is equivalent to the assumption that the potential vanishes for r ≤ r min and r ≥ r max . For such potentials, the integral equation Eq. 3 can be written as S L ′ L r ; k = − i k h l ′ 1 k r β L ′ L k − i k h l ′ 1 k r ∫ r m i n r d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k − i k j l ′ k r ∫ r r m a x d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k , where Eq. 6 was used. With β L ′ L k = δ L ′ L − ∫ r m i n r m a x d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k , which arises from Eq. 5 for the finite integration limits, Eq. 9 for S L ′ L r ;   k can be rewritten as S L ′ L r ; k = − i k h l ′ 1 k r δ L ′ L + i k h l ′ 1 k r ∫ r r m a x d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k − i k j l ′ k r ∫ r r m a x d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k . This shows that the irregular solutions can be expressed as S L ′ L r ; k = − i k j l ′ k r C L ′ L r ; k − i k h l ′ 1 k r D L ′ L r ; k , where the matrix functions C L ′ L r ;   k and D L ′ L r ;   k are defined as C L ′ L r ; k = ∫ r r m a x d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k and as D L ′ L r ; k = δ L ′ L − ∫ r r m a x d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k . From Eq. 12, the inner and outer boundary conditions are obtained as S L ′ L r ; k = − i k j l ′ k r C L ′ L r m i n ; k + h l ′ 1 k r D L ′ L r m i n ; k for r ≤ r m i n h l ′ 1 k r δ L ′ L for r ≥ r m a x .

Like Eq. 12, the regular solutions can be expressed as R L ′ L r ; k = j l ′ k r A L ′ L r ; k − i k h l ′ 1 k r B L ′ L r ; k with matrix functions A L ′ L r ;   k and B L ′ L r ;   k defined as A L ′ L r ; k = δ L ′ L − i k ∫ r r m a x d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ R L ″ L r ′ ; k and as B L ′ L r ; k = ∫ r m i n r d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ R L ″ L r ′ ; k . This can be shown by using Eq. 6 in Eq. 2, which results in R L ′ L r ; k = j l ′ k r δ L ′ L − i k j l ′ k r ∫ r r m a x d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ R L ″ L r ′ ; k − i k h l ′ 1 k r ∫ r m i n r d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ R L ″ L r ′ ; k . From Eq. 16, the inner and outer boundary conditions are obtained as R L ′ L r ; k = j l ′ k r A L ′ L r m i n ; k for r ≤ r m i n j l ′ k r δ L ′ L − i k h l ′ 1 k r B L ′ L r m a x ; k for r ≥ r m a x

2.3 Auxiliary integral equations

The use of integral equations instead of differential equations has been hindered in the past by the much larger computational work. When the interval from r min to r max is discretized by N mesh points, the integral equations given above can be converted into systems of linear algebraic equations with the dimension N . Thus, the computing time scales as N 3 whereas the computing time to solve linear differential equations typically scales only linearly with N . This means that the effort increases with the third power of the interval length r max   −   r min for the solution of linear integral equations, but only linearly with r max   −   r min for the solution of linear differential equations.

To overcome this problem, Greengard and Rokhlin [4] observed that the cubic scaling with r max   −   r min is avoided by dividing the interval into subintervals and by solving auxiliary integral equations locally in each subinterval. If the interval r min , r max is divided into N subintervals r n − 1 , r n with r 0 = r min and r N = r max and if p discretization points are used in each subinterval, the computing time scales as N p 3 for the solution of the auxiliary integral equations. Thus, it increases only linearly with the interval length r max   −   r min . Admittedly, the prefactor p 3 can be large, but this is not a serious drawback in view of current computer capabilities. The method of subintervals is based on the property that the integral equation Eq. 11 for the coupled irregular solutions can be written as S L ′ L r ; k = − i k j l ′ k r C L ′ L r n ; k − i k h l ′ 1 k r D L ′ L r n ; k − i k j l ′ k r ∫ r r n d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k + i k h l ′ 1 k r ∫ r r n d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k where the matrix function (Eqs 13, 14) are used for the particular value r n and the integrals arise from deviations of the matrix functions at r and r n . The idea is to solve Eq. 21 separately for each subinterval with r restricted as r n − 1 ≤ r ≤ r n by introducing auxiliary integral equations Y L ′ L n r ; k = j l ′ k r δ L ′ L − h l ′ 1 k r ∫ r r n d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ Y L ″ L n r ′ ; k + j l ′ k r ∫ r r n d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ Y L ″ L n r ′ ; k and Z L ′ L n r ; k = h l ′ 1 k r δ L ′ L − h l ′ 1 k r ∫ r r n d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ Z L ″ L n r ′ ; k + j l ′ k r ∫ r r n d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ Z L ″ L n r ′ ; k . The advantage of introducing these local solutions is that they do not depend on the unknown solution S L ′ L r ;   k , unlike Eq. 21 which contains the matrix functions Eq. 13 and Eq. 14. With the local solutions, which can be obtained numerically as described in Section 2.5, the irregular solution can be expressed in the interval r n − 1 , r n as S L ′ L ″ r ; k = − i k ∑ L Z L ′ L n r ; k C L L ″ r n ; k + Y L ′ L n r ; k D L L ″ r n ; k . This can be verified by multiplying Eq. 22 with D L L ′ ′ ′ ( r n ; k ) and Eq. 23 with C L L ′ ′ ′ ( r n ; k ) , which yields Y L ′ L n r ; k D L L ‴ r n ; k = j l ′ k r D L ′ L ‴ r n ; k − h l ′ 1 k r ∫ r r n d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ Y L ″ L n r ; k D L L ‴ r n ; k + j l ′ k r ∫ r r n d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ Y L ″ L n r ; k D L L ‴ r n ; k , Z L ′ L n r ; k C L L ‴ r n ; k = h l ′ 1 k r C L ′ L ‴ r n ; k − h l ′ 1 k r ∫ r r n d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ Z L ″ L n r ; k C L L ‴ r n ; k + j l ′ k r ∫ r r n d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ Z L ″ L n r ; k C L L ‴ r n ; k . Adding Eqs 25, 26 both multiplied with − i k and summed over L , leads to an integral equation which, compared with Eq. 21, contains the same source term and kernel. Thus, the left and right sides of Eq. 24 represent the same function.

It remains to calculate the matrix functions given by Eqs 13, 14 at the subinterval boundaries r n . This can be done recursively by recognizing that these functions satisfy the expressions C L ′ L r n − 1 ; k = C L ′ L r n ; k + ∫ r n − 1 r n d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k and D L ′ L r n − 1 ; k = D L ′ L r n ; k − ∫ r n − 1 r n d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ S L ″ L r ′ ; k and by using Eq. 24, which shows that the function S L ″ L r ′ ;   k can be expressed by the local solutions Y L ′ L n ( r ; k ) and Z L ′ L n ( r ; k ) . This leads to the integrals M L ′ L h Y n ; k = − i k ∫ r n − 1 r n d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ Y L ″ L n r ′ ; k M L ′ L h Z n ; k = − i k ∫ r n − 1 r n d r ′ r ′ 2 h l ′ 1 k r ′ ∑ L ″ V L ′ L ″ r ′ Z L ″ L n r ′ ; k M L ′ L j Y n ; k = − i k ∫ r n − 1 r n d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ Y L ″ L n r ′ ; k M L ′ L j Z n ; k = − i k ∫ r n − 1 r n d r ′ r ′ 2 j l ′ k r ′ ∑ L ″ V L ′ L ″ r ′ Z L ″ L n r ′ ; k which must be evaluated numerically as described in Section 2.5. By using Eqs 29–32 the expressions (Eqs 27, 28) can be written as recursion relations C L ′ L r n − 1 ; k = C L ′ L r n ; k + ∑ L ″ M L ′ L ″ h Z n ; k C L ″ L r n ; k + ∑ L ″ M L ′ L ″ h Y n ; k D L ″ L r n ; k D L ′ L r n − 1 ; k = D L ′ L r n ; k − ∑ L ″ M L ′ L ″ j Z n ; k C L ″ L r n ; k − ∑ L ″ M L ′ L ″ j Y n ; k D L ″ L r n ; k , starting from C L ′ L r N ;   k = 0 and D L ′ L r N ;   k = δ L ′ L .

The method of subintervals is particularly advantageous for potentials with a finite number of discontinuities in the radial direction. Such discontinuities arise, for instance, in full-potential Korringa–Kohn–Rostoker calculations where the angular integration in Eq. 4 must be cut off at boundaries of the atomic cells, which leads to jumps in the radial derivative of the matrix elements of the potential. If the interval boundaries r n are adapted to the discontinuities, the discontinuous behavior is treated without numerical approximations, which means that numerical errors only depend on the smoothness of the potential within the intervals. The method of subintervals is also advantageous from a computational point of view because the auxiliary functions Eqs 22, 23 and then the integrals in Eqs 29–32 can be calculated efficiently on multi-core processors separately for each value of n .

 2.4 Modified boundary conditions

In order to understand the necessity of modified boundary conditions for accurate density calculations, it is useful to consider the complex-contour integral n r = − 2 π lim r ′ → r I m ∫ − ∞ E F + i 0 + d ϵ G r , r ′ ; ϵ . which is used to calculate the density from the Green function for the Schrödinger equation (Eq. 1). Here, E F is the Fermi energy, which determines the total charge and i 0 + a positive infinitesimal imaginary quantity, which is used to avoid the singularities which exist for real values of ϵ . Around each atomic position, the Green function can be written as G r , r ′ ; k = ∑ L L ′ Y L r ̂ Y L ′ r ̂ ′ G L L ′ r , r ′ ; k where k is given by k = ϵ . The matrix function G L L ′ r , r ′ ;   k consists of two contributions, so-called “single-scattering” and “multiple-back-scattering” parts. While the back-scattering part is determined by coupled regular solutions alone, the single-scattering part contains the divergent irregular solutions in the form [8]: G L L ′ r , r ′ ; k = ∑ L ″ S L L ″ r > ; k R L ′ L ″ r < ; k . Here, r < and r > are defined as r < = min r , r ′ and r > = max r , r ′ , respectively. The effectiveness of the contour integral Eq. 35 arises from the fact that the contour can be chosen in the upper half of the complex- ϵ plane, where the integrand is an analytical function of ϵ , such that with only 20–30 mesh points on the contour [10], highly accurate density results are obtained even for large systems with many thousands of atoms.

Unfortunately, for the contour integral in Eq. 35, the irregular solutions are absolutely necessary to maintain the analytical behavior of the Green function as discussed in [11]. In a numerical treatment with the standard boundary conditions Eqs 15–20, the divergent behavior of the matrix function G L L ′ r , r ′ ;   k is given by − i k ∑ L ″ h l 1 k r D L L ″ r m i n ; k A L ′ L ″ r m i n ; k j l ′ k r ′ for r ≤ r ′ ≤ r min . The correct behavior − i k j l k r h l 1 k r ′ δ L L ′ is obtained from Eq. 80 derived in Appendix A. Comparison of Eqs 38, 39 shows that the transpose of the matrix A r min ;   k must be equal to the inverse of the matrix D r min ;   k . Numerically, this cannot be achieved because two different integral equations (Eqs 2, 3) are used to calculate these matrices. A solution for this problem is to apply modified irregular and regular solutions defined as S ∼ L ′ L r ; k = ∑ L ″ S L ′ L ″ r ; k D L ″ L − 1 r m i n ; k R ∼ L ′ L r ; k = ∑ L ″ R L ′ L ″ r ; k A L ″ L − 1 r m i n ; k . These solutions have the inner boundary conditions S ∼ L ′ L r ; k = − i k j l ′ k r ∑ L ″ C L ′ L ″ r m i n ; k D L ″ L − 1 r m i n ; k − i k h l ′ 1 k r δ L ′ L and R ∼ L ′ L r ; k = j l ′ k r δ L ′ L for r ≤ r min such that the divergent part of the matrix function G L L ′ r , r ′ ; k = ∑ L ″ S ∼ L L ″ r > ; k R ∼ L ′ L ″ r < ; k has the correct behavior given in Eq. 39.

The modified irregular solutions can be calculated by S ∼ L ′ L ″ r ; k = − i k ∑ L Z L ′ L n r ; k C ∼ L L ″ r n ; k + Y L ′ L n r ; k D ∼ L L ″ r n ; k , which is obtained from Eq. 24 by multiplication with the matrix D − 1 r min ;   k from the right. The matrices C ∼ ( r n ; k ) and D ∼ ( r n ; k ) , which are given by C ∼ L L ′ r n ; k = ∑ L ″ C L L ″ r n ; k D L ″ L ′ − 1 r m i n ; k and by D ∼ L L ′ r n ; k = ∑ L ″ D L L ″ r n ; k D L ″ L ′ − 1 r m i n ; k , satisfy the recursion relations Eqs 33, 34 with C and D replaced by C ∼ and D ∼ . The only difference is that the starting values are changed from C L ′ L r N ;   k = 0 and D L ′ L r N ;   k = δ L ′ L to C ∼ L ′ L ( r N ; k ) = 0 and D ∼ L ′ L ( r N ; k ) = D L ′ L − 1 ( r m i n ; k ) .

The disadvantage of the recursion for C ∼ and D ∼ is that the matrix D ∼ ( r N ; k ) is known only approximately, for instance, by using the numerically obtained result for D r min ;   k . This minor problem, however, is offset by the significant advantage that the error is known, which arises from the inaccuracy of D r min ;   k , from the numerical approximations necessary to solve the auxiliary integral equations and from roundoff errors. This error is given by the difference between D ∼ L ′ L ( 1 ) ( r m i n ; k ) , which is the numerical result obtained after the recursion, and δ L ′ L , which is the exact result for D ∼ L ′ L ( r m i n ; k ) . The knowledge of D ∼ L ′ L ( 1 ) ( r m i n ; k ) can be used in a second recursion starting from a better approximation for D ∼ ( r N ; k ) given as the product of D − 1 r min ;   k and the inverse of D ∼ ( 1 ) ( r m i n ; k ) . Further improved recursions can be added. In the present study, where electron densities up to l max = 8 were considered, rapid convergence was observed, and no more than two or three passes through the recursion were necessary.

While the straightforward use of repeated recursions for C ∼ and D ∼ successfully deals with numerical approximations, the problem of roundoff errors requires modification of the recursion relations such that not D ∼ but D ̂ L ′ L r n ; k = D ∼ L ′ L r n ; k − δ L ′ L is calculated directly. The modified recursion relations given by C ∼ L ′ L r n − 1 ; k = C ∼ L ′ L r n ; k + M L ′ L h Y n ; k + ∑ L ″ M L ′ L ″ h Z n ; k C ∼ L ″ L r n ; k + ∑ L ″ M L ′ L ″ h Y n ; k D ̂ L ″ L r n ; k D ̂ L ′ L r n − 1 ; k = D ̂ L ′ L r n ; k − M L ′ L ″ j Y n ; k − ∑ L ″ M L ′ L ″ j Z n ; k C ∼ L ″ L r n ; k − ∑ L ″ M L ′ L ″ j Y n ; k D ̂ L ″ L r n ; k lead to a matrix D ̂ L ′ L ( r m i n ; k ) with norm ‖ D ̂ L ′ L ( r m i n ; k ) ‖ ≪ 1 such that the required inverse of D ∼ L ′ L ( r m i n , k ) = D ̂ ( r m i n , k ) + δ L ′ L can be calculated reliably as − D ̂ ( r m i n , k ) + δ L ′ L . If necessary, further reduction of ‖ D ̂ L ′ L ( r m i n ; k ) ‖ can be achieved by evaluating Eqs 49, 50 with extended precision.

 2.5 Numerical treatment

The integral-equation method of Greengard [3], Greengard and Rokhlin [4], and Gonzales et al. [1] is based on expansions of the potential and solutions of the local integral equations (Eqs 22, 23) in Chebyshev polynomials T m x = cos m   arccos x . The method uses the property that the integral F τ = ∫ τ 1 d τ ′ f τ ′ of a function f τ = ∑ m = 0 M f m T m τ can be evaluated at the collocation points τ m = cos 2 m + 1 π 2 M + 1 , by matrix multiplication F τ m = ∑ m ′ = 0 M T m m ′ f τ m ′ . The collocation points τ m are the zeros of the Chebyshev polynomial T M + 1 τ . The matrix T is given by the product C − 1 S C where C and S are the “discrete cosine–transform” and “right spectral integration” matrices. The discrete cosine–transform matrix, which is given by C m m ′ = T m τ m ′ , connects the coefficients f m of the Chebyshev series Eq. 52 with values of the function at the zeros Eq. 53 by f m = ∑ m ′ = 0 M C m m ′ f τ m ′ . The right spectral integration matrix given in [1] connects the coefficients F m of the Chebyshev series F τ = ∑ m = 0 M F m T m τ with the coefficients f m by F m = ∑ m ′ = 0 M S m m ′ f m ′ . In Eq. 57, the M + 1 -th coefficient is neglected, which is justified because of the fast decay of F m with increasing m for sufficiently smooth functions. The matrix S has the non-zero elements S 00 = 1 , S 01 = 1 / 4 , S 10 = − 1 , S 12 = 1 / 2 , and, for m ≥ 2 , the non-zero elements S 0 m = 1 / ( 1 − m 2 ) , S m , m + 1 = 1 / 2 m , and S m , m − 1 = − 1 / 2 m . These values can be obtained by using the integration rules for Chebyshev polynomials. By using Eq. 54, the local integral equations (Eq. 22) are approximated by Y L ′ L n τ m ; k = h l ′ 1 k τ m δ L ′ L + r n − r n − 1 2 ∑ m ′ = 0 M ∑ L ″ A L ′ L ″ m m ′ Y L ″ L n τ m ′ ; k where the factor r n   −   r n − 1 / 2 comes from the substitution x = 2 r   −   r n − 1 / r n   −   r n − 1   −   1 , which transforms the interval r n − 1 , r n into − 1 , 1 . The matrix A is given by A L ′ L ″ m m ′ = − h l ′ 1 k τ m j l ′ k τ m ′ + j l ′ k τ m h l ′ 1 k τ m ′ T m m ′ τ m ′ 2 V L ′ L ″ τ m ′ . The system Eq. 59 of linear equations can be efficiently solved by standard linear algebra software. It has dimension ( l m a x + 1 ) 2 ( M + 1 ) and requires a computing effort that scales as ( l m a x + 1 ) 6 ( M + 1 ) 3 .

While, in principle, the subintervals can be chosen arbitrarily, the choice should be adapted to the divergent behavior of the irregular solutions for r   →   0 . A suitable choice is given by the prescription r n − 1 = α r n with α = ( r 0 / r N ) 1 / N = ( r m i n / r m a x ) 1 / N . The transformation to the standard expansion interval − 1 , 1 is obtained by the substitution r = 1 2 r n ( 1 − α ) τ + 1 + α . For inverse powers of r , the substitution leads to ∫ α r n r n 1 r l d r = 2 l − 1 1 − α l − 1 r n l − 1 ∫ − 1 1 1 a + τ l d τ with a = 1 + α / 1   −   α . Here, the integrand and the integration limits for the integral over τ do not depend on n . Thus, without changing the Chebyshev-expansion order, the same relative accuracy is obtained for all intervals.

 3 Numerical investigation

The numerical performance is investigated for two examples: a constant potential, which is analytically solvable, and a realistic non-spherical potential, which is obtained by density-functional electronic-structure calculations for an ordered nickel–titanium alloy. It is shown for the constant potential that accurate bound-state energies and wavefunctions can be obtained from the irregular solutions calculated by the integral-equation approach and that the error of calculated bound-state energies decreases exponentially with the order of the Chebyshev expansion. For the NiTi alloy, it is shown that the irregular solutions obtained can be used in complex-contour integrations to calculate the density from the full-potential multiple-scattering Green function. Thus, such calculations can be performed straightforwardly for systems with many atoms in contrast to other treatments suggested in the past [11–13], which are rather elaborate and unlikely to be useful for systems with more than a few atoms.

3.1 Bound states for a constant potential

The standard method for calculating bound states is based on the property that regular solutions of the Schrödinger equation vanish at infinity if they are evaluated for the correct bound-state energy. For trial energies, the differential equation is solved from the inside, starting with the correct power r l at r = 0 , and from the outside starting with 0 at a large value of r . The bound-state energy and wavefunction are found if, for the chosen trial energy, the logarithmic derivatives of both solutions continuously match at an intermediate value of r . The alternative method suggested here is based on the property that irregular solutions of the Schrödinger equation vanish at the origin if they are evaluated for the correct bound-state energy. For a constant potential, because of its spherical symmetry, the irregular solutions are decoupled S L ′ L r ;   k = S l ′ r ;   k δ L ′ L and the inner boundary condition Eq. 15 contains diagonal matrices C L ′ L r min ;   k = C l ′ r min ;   k δ L ′ L and D L ′ L r min ;   k = D l ′ r min ;   k δ L ′ L . The condition for a bound state is then given by D l ′ r min ;   k = 0 which eliminates the diverging Hankel functions in Eq. 15 for r ≤ r min . In mathematical scattering theory, the function D l r m i n ; k = 1 − ∫ r m i n r m a x d r ′ r ′ 2 j l k r ′ V 0 S l r ′ ; k is known as the “Jost function.” Its analytical properties in the complex- k plane are comprehensively discussed in [14], where it is explained that bound states correspond to zeros of D l r min ;   k for k values on the positive imaginary axis. The determination of these zeros is a one-dimensional root-finding problem treated here by Ridders’ method [15].

Numerical results for bound-state energies and wavefunctions are shown in Table 1 and Figure 1 for an attractive potential V 0 = − 16 Ry, which is confined to a spherical shell between r min = 0.00001   a B and r max = 3 a B . This potential has bound states up to l = 8 . The energies in Table 1 were obtained using N = 10 intervals of equal length and order M = 10 for the Chebyshev expansions. They deviate by less than 2 × 10 − 9 from the exact energies determined from the zeros of D l r min ;   k using the analytical expressions given in Appendix A. For comparison with values given in the literature, for example, in [16] where the same potential is treated by sinc-interpolants for r min = 0 , it is useful to know that the same digits as in Table 1 are obtained if r min is chosen smaller than 0.00001   a B . This is a consequence of the third-order dependence on r min , which can be established from the expressions derived in Appendix A.

Bound-state energies for different l values in Rydberg units for a potential of depth −16 Ry confined to a spherical shell between r min = 0.00001   a B and r max = 3   a B .

l e 1 e 2 e 3 e 4
0 −15.067032975 −12.287216857 −7.738182446 −1.734147318
1 −14.093970355 −10.407767360 −5.037170230
2 −12.869064652 −8.2824156334 −2.179455290
3 −11.405365235 −5.9291642729
4 −9.7123791747 −3.3710840481
5 −7.7979942918
6 −5.6694973028
7 −3.3343625870
8 −0.8012457212

(A) Deviations from the exact result for the lowest bound-state energy of l = 0 . The black, red, green, and blue curves are for N = 3 , N = 5 , N = 10 , and N = 30 intervals. (B) Irregular solutions, multiplied by r and normalized to one at r = 3   a B , for selected values of l . The black, red, green, blue, and orange curves are for the highest bound-state energies of l = 0 , l = 1 , l = 3 , l = 5 and l = 8 . The constant potential used in the calculations has a depth of −16 Ry and is confined to a spherical shell between r = 0.00001   a B and r = 3   a B .

Figure 1 shows how the energy of the lowest bound state for l = 0 converges with the order of the Chebyshev expansion. The convergence is very similar for the other bound states. The deviations decrease exponentially with the order, and accurate results are already obtained with a small number of intervals by using a sufficiently large order. A precision of 10 − 10 is achieved for N = 3 intervals and order M = 14 , leading to 45 collocation points. For a larger number of intervals, where a smaller order gives the same precision, more collocation points are necessary. This is, however, not important for the computational effort, which scales as N M + 1 3 , so that the effects of an increase of N and a decrease M are practically canceled.

Figure 1 also shows the irregular solutions for the highest bound-state energies for selected values of l calculated with N = 10 and M = 8 . For the presentation, they are multiplied with r and normalized to 1 at r = 3   a B . They exponentially decay for r   >   3   a B and, as a consequence of D l r min ;   k = 0 , regular at r = 0 . Thus, they satisfy the requirements specifying bound-state wavefunctions. It should be noted that, unlike the example shown, the condition D l r min ;   k = 0 is not always obtained numerically with sufficient precision to conceal the divergent behavior at r = 0 . It is then more appropriate to determine the bound-state wavefunction not from the irregular solution but from the regular solution by using the property that irregular and regular solutions are multiples of one another at the bound-state energy, as discussed in [14].

 3.2 Electron density for NiTi

Metallic alloys of nickel and titanium have interesting and technologically important mechanical properties, like the shape memory effect. If NiTi in its low-temperature B19’ phase is deformed and heated, it returns to its original form, which persists on cooling. Because of the low symmetry of the P 21 /m space group, the spherical-harmonics expansion of the potential V r = ∑ L V L r Y L r contains non-zero terms for all values of l in contrast to high symmetry systems like copper or silicon. The potential for NiTi was determined self-consistently for l max = 3 . The exchange-correlation potential was treated in the Vosko–Wilks–Nusair parametrization [17] and a Monkhorst–Pack grid [18] with 16× 16 × 16 points applied for the Brillouin zone integrations. The experimental lattice structure as given in [19] was used with a = 289.8 pm, b = 464.6 pm, c = 410.8 pm, γ = 97.8 °, and Wyckoff (2e) positions (±0.0372, ±0.6752, 1/4) for Ni and (±0.4176, ±0.2164, 1/4) for Ti. The order for the Chebyshev expansion was chosen as M = 8 and the subintervals were chosen as follows. On the outside of the inscribed spheres of the atomic Voronoi cells, the intervals were determined by the kinks of the shape functions (for an explanation, see [20]). On the inside, 30 intervals were used between r min = 0.00001   a B and r = 1.2   a B with increasing length corresponding to α = 0.00001 / 1.2 1 / 30 = 0.677164 and eight intervals above r = 1.2   a B with equal length. The total number of intervals was 64, leading to an overall 576 radial mesh points. It should be emphasized that the non-spherical potential was used on all radial mesh points, and no cutoff of the non-spherical part near the atomic centers was applied. Previously, such cutoffs were always necessary as explained, for instance, in [21].

The self-consistent potential determined in this way was used in calculations for the electron density for l max ≤ 8 . In order to save computer resources, a reduced Monkhorst–Pack grid with 6 × 6 × 6 points was applied. Results for the density of the valence electrons near the center of an Ni atom are shown in Figure 2 for l max = 3 and l max = 8 . The insets are blowups on a 100× smaller range. The standard boundary conditions with the recursion relations Eqs 33, 34 lead to the results shown by blue curves. They deviate from the correct behavior below r = 0.001   a B for l max = 3 and below r = 0.1   a B for l max = 8 . These deviations degrade the self-consistency procedure in density-functional calculations unless they are somehow removed by extrapolation, which is cumbersome for large systems, or completely neglected near the atomic centers. Such neglect might be justified for l max = 3 , where the affected volume is a tiny part of the total volume but might be unreasonable for l max = 8 , where the affected volume is non-negligible. Straightforward use of the modified boundary conditions leads to the results shown by green curves. In the left picture, for l max = 3 , the green curve, which is hidden under the black curve, exhibits no divergence at the origin. In the right picture, for l max = 8 , the green curve begins to diverge at r = 0.004   a B , which is considerably smaller than r = 0.1   a B where the blue curve, obtained from the unmodified solutions, begins to diverge. The use of the modified boundary conditions together with one pass through the modified recursion relations (Eqs 49, 50) leads to the orange curves. In the left picture, the orange curve is again hidden under the black curve, while in the right picture, it is considerably better than the green curve by shifting the beginning of the divergence from r = 0.004   a B to r = 0.0016   a B . A second pass through the modified recursion relations (Eqs 49, 50) gives only a small improvement seen in the red curve, which begins to diverge at about r = 0.0013   a B .

Electron density near the center of an Ni atom in NiTi plotted in (1,0,0) direction for l max = 3 (A) and l max = 8 (B). The insets display the densities on 100× smaller range. The different curves are explained in the text.

The question of whether better results can be obtained by using a larger number of intervals or a higher order of the Chebyshev expansion was investigated by doubling the number of intervals below r = 1.2   a B and by increasing the order to M = 16 . The results thus obtained are practically identical to those shown in Figure 2. This indicates that the treatment of the auxiliary integral equations with double precision floating-point arithmetic is accurate enough. Whether better results can be obtained by treating the recursion relations more accurately was investigated by using quadruple precision instead of double precision. Then, instead of below r = 0.04   a B , r = 0.0016   a B , and r = 0.0013   a B , the density results diverge only below r = 0.003   a B , r = 0.0004   a B , and r = 0.00001   a B . Thus, divergence-free densities as shown by the black curves can be obtained for NiTi at least up to l max = 8 .

 4 Summary and outlook

An integral-equation approach was presented for the calculation of irregular solutions of the Schrödinger equation for non-spherical potentials. It was shown how expansions in Chebyshev polynomials can be used to convert the integral equations into systems of algebraic equations, which can be solved by standard software. For that purpose, no explicit construction of the Chebyshev series is needed but only function values at the zeros of the Chebyshev polynomials. It was explained that the numerical effort is considerably reduced by a subinterval technique suggested by Greengard and Rokhlin and that this technique, with appropriately adapted intervals, is beneficial for potentials with a finite number of radial discontinuities; this is because the numerical precision is determined by the smoothness of the potentials between the discontinuities but not by the discontinuous behavior. A numerical investigation was presented for a constant potential and for a realistic non-spherical potential, which was obtained by density-functional calculations for a nickel–titanium alloy. It was shown that accurate bound-state energies can be obtained from the calculated irregular solutions. It was explained how a precise description of the divergent behavior of the coupled irregular solutions can be obtained such that accurate density calculations by complex-contour integrations are possible.

The approach presented can be extended in several directions. It is not restricted to the non-relativistic Schrödinger equation treated here but is also useful for including scalar-relativistic and spin-orbit-coupling effects [22] and for full-relativistic calculations by the Dirac equation (Eq. 23). It can be extended to calculate bound-state energies for non-spherical potentials, although with more effort because of nearby roots caused by degeneracy splitting. It also can be extended to calculate scattering resonances which are determined by zeros of Jost functions in the complex- k plane. These zeros can be obtained by contour integrals, for which a Fortran package is available [24]. Calculations of exchange-correlation and Coulomb potentials, which require differentiations and integrations, are easily done with spectral differentiation and integration matrices without introducing additional numerical approximations. An interesting subject for further research is the question of how much numerical precision is needed to evaluate the recursion relations for higher values of l max beyond l max = 8 , which was the limit set by the present capabilities of the applied KKRnano code.

 Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

 Author contributions

RZ: methodology, software, writing–original draft and writing–review and editing.

 Funding

The author declares that no financial support was received for the research, authorship, and/or publication of this article.

I would like to thank Paul F. Baumeister for a critical reading of the manuscript and valuable discussions to clarify the mathematical ideas and their presentation.

 Conflict of interest

Author RZ was employed by Forschungszentrum Jülich GmbH.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Gonzales RA Eisert J Koltracht I Neumann M Rawitscher G . Integral equation method for the continuous spectrum radial Schrödinger equation. J Comput Phys (1997) 134:134–49. 10.1006/jcph.1997.5679 Clenshaw CW Curtis AR . A method for numerical integration on an automatic computer. Numer Math (1960) 2:197–205. 10.1007/bf01386223 Greengard L . Spectral integration and two-point boundary value problems. SIAM J Numer Anal (1991) 28:1071–80. 10.1137/0728057 Greengard L Rokhlin V . On the numerical solution of two-point boundary value problems. Commun Pure Appl Math (1991) 44:419–52. 10.1002/cpa.3160440403 Hatada K Hayakawa K Benfatto M Natoli CR . Full-potential multiple scattering theory with space-filling cells for bound and continuum states. J Phys Condens Matter (2010) 22:185501. 10.1088/0953-8984/22/18/185501 Zeller R Dederichs PH Újfalussy B Szunyogh L Weinberger P . Theory and convergence properties of the screened Korringa-Kohn-Rostoker method. Phys Rev B (1995) 52:8807–12. 10.1103/physrevb.52.8807 Thiess A Zeller R Bolten M Dederichs PH Blügel S . Massively parallel density functional calculations for thousands of atoms: KKRnano. Phys Rev B (2012) 85:235103. 10.1103/physrevb.85.235103 Zeller R . Projection potentials and angular momentum convergence of total energies in the full-potential Korringa-Kohn-Rostoker method. J Phys Condens Matter (2013) 25:105505. 10.1088/0953-8984/25/10/105505 Ershov SN Vaagen JS Zhukov MV . Modified variable phase method for the solution of coupled radial Schrödinger equations. Phys Rev C (2011) 84:064308. 10.1103/physrevc.84.064308 Zeller R Deutz J Dederichs PH . Application of complex energy integration to selfconsistent electronic structure calculations. Solid State Commun (1982) 44:993–7. 10.1016/0038-1098(82)90320-9 Kotani T Akai H . KKR-ASA method in exact exchange-potential band-structure calculations. Phys Rev B (1996) 54:16502–14. 10.1103/physrevb.54.16502 Ogura M Akai H . The full potential Korringa-Kohn-Rostoker method and its application in electric field gradient calculations. J Phys Condens Matter (2005) 17:5741–55. 10.1088/0953-8984/17/37/011 Rusanu A Stocks GM Wang Y Faulkner JS . Green’s functions in full-potential multiple-scattering theory. Phys Rev B (2011) 84:035102. 10.1103/physrevb.84.035102 Newton RG . Analytic properties of radial wave functions. J Math Phys (1960) 1:452–347. 10.1063/1.1703680 Ridders CJF . A new algorithm for computing a single root of a real continuous function. IEEE Trans Circuits Syst (1979) 26:979–80. 10.1109/tcs.1979.1084580 Annaby MH Asharabi RM . Sinc-interpolants in the energy plane for regular solution, Jost function, and its zeros of quantum scattering. J Math Phys (2018) 59:013502. 10.1063/1.5001078 Vosko SH Wilk L Nusair M . Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis. Can J Phys (1980) 58:1200–11. 10.1139/p80-159 Monkhorst HJ Pack JD . Special points for Brillouin-zone integrations. Phys Rev B (1976) 13:5188–92. 10.1103/physrevb.13.5188 Kudoh Y Tokonami M Miyazaki S Otsuka K . Crystal structure of the martensite in Ti-49.2 at.%Ni alloy analyzed by the single crystal X-ray diffraction method. Acta Mater (1985) 33:2049–56. 10.1016/0001-6160(85)90128-2 Stefanou N Akai H Zeller R . An efficient numerical method to calculate shape truncation functions for Wigner-Seitz atomic polyhedra. Comp Phys Commun (1990) 60:231–8. 10.1016/0010-4655(90)90009-p Zeller R . Large scale supercell calculations for forces around substitutional defects in NiTi. Phys Status Solidi B (2014) 251:2048–54. 10.1002/pssb.201350406 Bauer DSG Development of a relativistic full-potential first-principles multiple scattering Green function method applied to complex magnetic textures of nano structures at surfaces. Ph.D. thesis. Aachen, Germany: RWTH Aachen University (2013). Geilhufe M Achilles S Köbis MA Arnold M Mertig I Hergert W Numerical solution of the relativistic single-site scattering problem for the Coulomb and the Mathieu potential. J Phys Condens Matter (2015) 27:435202. 10.1088/0953-8984/27/43/435202 Kravanja P Van Barel M Ragos O Vrahatis MN Zafiropoulos FA . Zeal: a mathematical software package for computing zeros of analytic functions. Comp Phys Commun (2000) 124:212–32. 10.1016/s0010-4655(99)00429-4 Appendix A Analytical results for constant potentials

For a potential, which has a constant value V 0 for r min ≤ r ≤ r max and vanishes for r   <   r min and r   >   r max , the irregular solutions can be calculated analytically. Because of the spherical symmetry of the potential, the irregular solutions for different L channels are decoupled and can be written as S L ′ L r ; k = S l ′ r ; k δ L ′ L . With k 1 = k 2 − V 0 the functions are given by S l r ; k = − i k h l 1 k r for r ≥ r m a x c l j l k 1 r + d l h l 1 k 1 r for r ≤ r m a x . The constants c l and d l are determined by the conditions c l j l k 1 r m a x + d l h l 1 k 1 r m a x = − i k h l 1 k r m a x − k 1 c l j l + 1 k 1 r m a x − k 1 d l h l + 1 1 k 1 r m a x = i k 2 h l + 1 1 k r m a x which guarantee that the solutions are continuous and continuously differentiable at r max . Eq. A4 is obtained from Eq. A3 by differentiation using the formula f l ′ ( x ) = − f l + 1 ( x ) + ( l / x ) f l ( x ) for derivatives of Bessel and Hankel functions. The terms arising from l / x f l x are omitted in Eq. A4 because they are a simple multiple of Eq. A3. Solving Eq. A3 and Eq. A4 for c l and d l leads to c l = k k 1 r m a x 2 k 1 h l + 1 1 k 1 r m a x h l 1 k r m a x − k h l 1 k 1 r m a x h l + 1 1 k r m a x d l = − k k 1 r m a x 2 k 1 j l + 1 k 1 r m a x h l 1 k r m a x − k j l k 1 r m a x h l + 1 1 k r m a x Inserting Eq. A2 into Eq. 62 and using the standard results for integrals of products of Bessel and Hankel functions of the same order and different arguments leads to D l r m i n ; k = D l 1 + D l 2 with D l 1 = 1 − c l [ k r m a x 2 j l + 1 k r m a x j l k 1 r m a x − k 1 r m a x 2 j l k r m a x j l + 1 k 1 r m a x ] − d l k r m a x 2 j l + 1 k r m a x h l 1 k 1 r m a x − k 1 r m a x 2 j l k r m a x h l + 1 1 k 1 r m a x D l 2 = c l [ k r m i n 2 j l + 1 k r m i n j l k 1 r m i n − k 1 r m i n 2 j l k r m i n j l + 1 k 1 r m i n ] + d l k r m i n 2 j l + 1 k r m i n h l 1 k 1 r m i n − k 1 r m i n 2 j l k r m i n h l + 1 1 k 1 r m i n In Eq. A8 the constants c l and d l can be eliminated by using Eq. A3 in the terms proportional to k and Eq. A4 in the terms proportional to k 1 . This leads to D l 1 = 1 + i k 2 r m a x 2 j l + 1 k r m a x h l 1 k r m a x − j l k r m a x h l + 1 1 k r m a x From the Wronskian relation j l + 1 x h l 1 x − j l x h l + 1 1 x = i x 2 for spherical Bessel functions it follows that D 1 vanishes and that D l r min ;   k is given by Eq. A9.

 Green function at the origin

The behavior of the Green function for arguments smaller than r min can be determined from the Dyson equation G r , r ′ ; k = g r , r ′ ; k + ∫ d r ″ g r , r ″ ; k V r ″ G r ″ , r ′ ; k . by using Eq. 36 for G r , r ′ ;   k and the corresponding result g r , r ′ ; k = ∑ L Y L r ̂ Y L r ̂ ′ g l r , r ′ ; k for g r , r ′ ;   k . This leads to ∑ L L ′ Y L r ̂ Y L ′ r ̂ ′ G L L ′ r , r ′ ; k = ∑ L Y L r ̂ Y L r ̂ ′ g l r , r ′ ; k + ∑ L L ′ L ″ ∫ r m i n r m a x d r ″ r ″ 2 Y L r ̂ g l r , r ″ ; k V L L ″ r ″ G L ″ L ′ r ″ , r ′ ; k Y L ′ r ̂ ′ where the integration over the angles was done by using the definition Eq. 4 for the potential matrix elements. With the orthogonality of the spherical harmonics the angular coordinates in Eq. A14 can be eliminated, which yields G L L ′ r , r ′ ; k = g l r , r ′ ; k δ L L ′ + ∫ r m i n r m a x d r ″ r ″ 2 g l r , r ″ ; k ∑ L ″ V L L ″ r ″ G L ″ L ′ r ″ , r ′ ; k For r ≤ r ′ ≤ r min , the use of Eq. 6 and Eq. 37 leads to G L L ′ r , r ′ ; k = − i k j l k r h l 1 k r ′ δ L L ′ − i k j l k r ∑ L ″ L ‴ R L ′ L ‴ r ′ ; k ∫ r m i n r m a x d r ″ r ″ 2 h l 1 k r ″ × V L L ″ r ″ S L ″ L ‴ r ″ ; k By using Eq. 13 the final result is given by G L L ′ r , r ′ ; k = − i k j l k r h l 1 k r ′ δ L L ′ − i k j l k r ∑ L ‴ R L ′ L ‴ r ′ ; k C L L ‴ r m i n ; k Here the only divergent expression is the first term.