Ju n 20 19 A generating integral for the matrix elements of the Coulomb Green ’ s function with the Coulomb wave functions

The reduced Coulomb Green’s function is used in Rayleigh-Schrödinger perturbation theory for hydrogenic systems to analytically evaluate matrix elements [1–3]. In addition, it has recently been demonstrated [4] that a multi-electron atom can be effectively described via a very simple model using an effective charge. In this approach one starts from a Hamiltonian of a multi-electron atom in the secondary-quantized representation, written in the hydrogen-like basis with an effective charge Z being a free parameter. The effective charge depends only on the set of occupation numbers of a given state and a charge of an atom. Then one constructs perturbation theory (PT), where corrections to energy levels and wave functions are given in terms of the matrix elements of a reduced Coulomb Green’s function (RCGF) Gnl(r, r ) with the hydrogen like wave functions. By expanding all functions in spherical harmonics Ylm(Ω) the angular integrations can be performed in closed form and are expressed through the 3-j symbols, while the radial integration can be further reduced to the computation of the generating integral of the form


A. Outline
The reduced Coulomb Green's function is used in Rayleigh-Schrödinger perturbation theory for hydrogenic systems to analytically evaluate matrix elements [1][2][3]. In addition, it has recently been demonstrated [4] that a multi-electron atom can be effectively described via a very simple model using an effective charge. In this approach one starts from a Hamiltonian of a multi-electron atom in the secondary-quantized representation, written in the hydrogen-like basis with an effective charge Z * being a free parameter. The effective charge depends only on the set of occupation numbers of a given state and a charge of an atom. Then one constructs perturbation theory (PT), where corrections to energy levels and wave functions are given in terms of the matrix elements of a reduced Coulomb Green's function (RCGF) G nl (r, r ′ ) with the hydrogen like wave functions. By expanding all functions in spherical harmonics Y lm (Ω) the angular integrations can be performed in closed form and are expressed through the 3-j symbols, while the radial integration can be further reduced to the computation of the generating integral of the form which is convergent for q ≥ 0 and q ′ ≥ 0. In Eq. (1) the letter l denotes the index of angular expansion. In contrast to the single particle Coulomb Green's function (CGF) the index n describes the state that was subtracted from the CGF. In principle, one can use the direct numerical integration of Eq. (1). However, when the indices n and l in Eq. (1) are large, which is the typical situation for Rydberg states [5] the direct numerical integration of Eq. (1) becomes very inefficient or even impossible due to the increasing number of nodes of the integrand. In addition, if repeated evaluations of the generating integrals are required the efficient scheme of computation of Eq. (1) is desirable. Consequently, in our work we analytically evaluate the generating integral K nl (β, β ′ ) and integral moments J nl (β, r ′ ) = ∞ 0 G nl (r, r ′ )r ′q e −βr ′ dr ′ (2) for all values of the parameters q and q ′ when the integrals are convergent.
The article is organized in the following way. For the readers convenience in Sec. II B we summarize the main results without any derivations. In Sec. II C we compare the evaluation time based on our analytical expressions with the direct numerical evaluation of the integrals. In Secs. III-V we provide the details of all derivations.

B. Relationship to previous research
Our starting point is the work of Johnson and Hirschfelder [6], who derived the explicit expressions for the RCGF in terms of elementary functions and computed integral moments ∞ 0 dr ′ G nl (r, r ′ )r ′k+2 e −Zr ′ /n (3) of the RCGF. In addition, we mention the work of Hill and Huxtable [7], who evaluated the generating integral for the case of q, q ′ = l + 1 and provided the recurrence relations, which allow one to compute the generating integral with the increasing powers of r, r ′ starting from q, q ′ = l + 1. Their result is based on treating the generating integral as a Laplace transform of β and β ′ and solving the differential equation for this Laplace transform. However, the generating integral in the range 0 ≤ q, q ′ ≤ l has not been evaluated. Consequently, here we evaluate the generating integral and integral moments for this case. In addition, we provide the closed form result for values q, q ′ > l + 1, which is directly applicable without the use of recurrence relations.
In our work we follow the approach of the direct integration of the RCGF with the corresponding powers of r, r ′ and exponentials. The main complications come from the fact, that when 0 ≤ q, q ′ < l + 1, the individual terms of the RCGF possess the singularities that are explicitly cancelled only upon summing all expressions together.

A. Definitions
The Hydrogen-like wave functions ψ(r) satisfy the Schrödinger equation Here E is the energy of the system, Z is the charge of the nucleus and the atomic units are employed. The variables in Eq. (4) are separated in spherical coordinates [8]. Consequently, the eigenfunctions of a discrete spectrum are represented as a product where with L 2l+1 n−l−1 (t) being the associated Laguerre polynomials [9] and Y lm (Ω) the spherical harmonics [8]. An analogous expression can be obtained for the continuous spectrum as well [8].
The associated Laguerre polynomials can be expanded according to their definition, thus leading to where the binomial is defined as n k = n!/(k!(n − k)!). In Eqs. (4)-(8) n is an integer, n > 0, l is an integer, 0 ≤ l ≤ n − 1 and m is an integer, −l ≤ m ≤ l. The bound states of Eq. (4) are described by energy eigenvalues E n = −Z 2 /(2n 2 ). The Green's function of the Hydrogen like atom satisfies the equation [ and can be expanded over the eigenfunctions of Eq. (4) Here the sum with asterisk denotes the summation over both the discrete and continuous spectra. When E equals the energy of the bound state the CGF has a pole. The radial CGF G l (r, r ′ ; E) satisfies the equation [6] (H l − E)G l (r, r ′ ; E) = − 1 2r and is expressed in terms of the Whittaker functions [10] G l (r, with Γ(t) being the gamma function [9]. The RCGF is defined as a CGF from which a state with the principal quantum number n is subtracted [2] G n (r, The summation in RCGF G nl (r, r ′ ) starts from l + 1. Therefore, one needs to distinguish the two distinct cases. When l ≤ n − 1 the term with n ′ = n must be explicitly excluded. For the opposite situation of l ≥ n no such term appears anyway, thus making the restriction n = n ′ redundant.
As shown in [6], for the case of l ≥ n the RCGF, satisfying the equation is given by (t = 2Zr/n, t ′ = 2Zr ′ /n): while for l ≤ n − 1, the RCGF satisfying an equation reads: where B. The main result

Generating integrals
When l + 1 ≤ n the closed form main result is given via 1 where α = λ + 1 n and α ′ = λ ′ + 1 n and the constants A i1,i2 and B i are given by: where Ψ(n) is the logarithmic derivative of the Gamma function [9] (also known as the zeroth polygamma function) and p F q (a 1 , . . . , a p ; b 1 , . . . , b q ; z) is the generalized hypergeometric function, which in this case can be evaluated as a finite series (see Eq. (22)). The function f b a (x, y) is defined as follows: • a = −n is a non-positive integer and a + b > 0: • a + b = −m is a non-positive integer and a > 0 • a = −m and b = −n are both non-positive integers When y = 0 the function f b a (x, 0) becomes: • a + b > 0 and a = 0: • b > 0 and a = 0 • a + b = −m is a non-positive integer and a = 0: • b = −m is a negative integer and a = 0 The function g b a (x, y) is defined as follows: • a > 0 and b > 0 • a = −n is a non-positive integer The function I q,q ′ is defined as: and when y = 1 In addition, in Eqs. (20)-(25b) 3F2 (a, b, c; d, e; z) is the generalized regularized hypergeometric function [9] and in Eq. (26b) H q is a harmonic number:

Integral moments
The main result for the integral moments for the case n ≥ l + 1 reads (x = Zr, λ = β/Z) where the function F is defined as: and in the case of q being a non-positive integer: In addition, the function I is defined as: and for x = 1 The closed form result for the case n ≤ l is given via (

C. Comparison of evaluation time of an analytical calculation with the numerical one
In this section we compare the evaluation time of the generating integral via analytical expressions (20) and (35) with the direct numerical evaluation of the integral when the term with E n = −Z 2 /(2n 2 ) has been explicitly subtracted from the Green's function Eq. (12). Since this expression is divergent, we introduce a parameter δ and substitute the energy as E n − iδ, effectively regularizing the numerical integration. That is, we evaluate the integrals from the expression G l (r, Since the resulting answer for the generating integral is δ dependent we always ensure that the numerical value is independent of the δ by varying the parameter δ. We have found out that in all relevant cases the total evaluation time through analytical expressions is of the order of 0.001-0.1 seconds on Intel 2600k 3.4GHz processor. On the other hand the direct numerical integration of the Green's function for some finite values of δ becomes very inefficient for large values of n and l, which is a relevant situation in the description of Rydberg states [5]. This  Table I. Comparison of computational times of the generating integral for different values of the parameters n, l, q and q ′ of analytical expressions (20) and (39) with the direct numerical integration using Mathematica [11], when the energy is shifted by iδ from its resonant value. The last column shows the value of δ required to obtain the accurate results to four significant figures. NaN means that the integral did not converge to the correct value for any values of δ.  happens due to the increasing number of nodes of the integrand, resulting in oscillatory behavior. Consequently, the accurate evaluation of the integral demands the smaller and smaller values of δ to keep the constant accuracy, since in many cases large positive values are almost completely cancelled by large negative values. Therefore, the precise result would require a forbiddingly accurate evaluation of the integrand at every point. In Tab. I we provide the summary of the evaluation times of the generating integral for different values of parameters n, l, q and q ′ . As can be observed from the table the numerical evaluation is several orders of magnitude slower then our analytical expressions. As follows from Tab. I the evaluation time increases significantly when n ≫ l. Consequently, the evaluation time of some high-n Rydberg states is still large and requires further optimization. For example, our analytic results allow one to derive asymptotic expressions for large values of n.
In Tab. II we compare the evaluation time of analytically computed integral moments with the numerically computed ones for some values of the parameters q, n, l and λ. In this simulation, we numerically evaluated the integral moments at one hundred different values of r in order to compare with an analytically calculated curve. As in the situation of the generating integral the evaluation time of the direct numerical integration is a few orders of magnitude slower.
Moreover, as was explained in the introduction the generating integral can be used to calculate the second order single-electron corrections to energies of neutral atoms and ions in the effective charge model of a multi-electron atom [4]. Therefore, we compare a few cases of evaluation times between our analytical and numerical approaches, which is given in Tab. III. Typically a fully sequential version of an effz program from Ref. [4] implemented in Mathematica requires one order of magnitude larger times for the direct numerical evaluation as compared to the analytical expressions.

A. Evaluation of the auxiliary integrals
First, we evaluate the following auxiliary integral convergent whenever Re a > 0, Re(a + b) > 0 and Re(λ + λ ′ ) > 0. Here 2 F 1 (a, b, c, z) is the Gauss Hypergeometric function [9]. Let us introduce the following notation In the case when a and b are positive integers the expression (33) simplifies if one uses the contingous relations [9] for 2 F 1 (a, b, c, z). This allows us to rewrite it as a finite sum Furthermore, for b = 0 and a a positive integer we get: Eq. (33) is useful when working with modern computer algebra software such as Mathematica [11] that can automatically simplify 2 F 1 (a, b, c, z) for given integer parameters, while (34) and (35) are useful for evaluating u b a with finite precision arithmetic, as they require evaluating an explicitly finite amount of terms.

B. Laurent series of u b a
During the calculation of the Green function we will need an expansion of u b a functions in the Laurent series in the neighbourhood of the singularities when a = −m + δ and a + b = −n + δ. In order to obtain the Laurent series it is most convenient to start from the infinite series representation (33), and separate singular terms. For this we split the series into two parts, namely i = 0..n − 1 and n..∞. The latter part is finite and is reduced to: where 3F2 (a 1 , a 2 , a 3 , b 1 , b 2 , z) is the reduced generalized hypergeometric function. Therefore, the divergent terms appear only in the first part. We now introduce the following useful relations. The expansion of the gamma function at positive integers gives: while for negative: Other relevant functions are and In addition, we use the following expansions This allows us to find the Laurent series of u b a in four different cases depending on the signs of a and b: • When a is a non-positive integer, we get exactly one divergent term (with index i = −a): • When a + b is a non-positive integer and a is positive, the first −a − b terms diverge: • When a + b is a non-positive integer and b is positive, the first −a − b temrs diverge as Γ(−m) and independently the term with i = −a diverges as 1/(a − a): • When a and b are both non-positive integers, the first −a − b terms diverge as Γ(−m − n) with the i = −a term diverging faster, as Γ(−m − n)/(a − a): Lastly, it is required to find the Laurent series representation of the function u ′b a defined as In the case of a being a non-positive integer and b > 0, we get: Finally, the introduced functions u b a allow us to evaluate the following integral: C. Derivation for l ≥ n In this case the derivation is straightforward. We employ the definition of the RCGF (15) and shift the index of summation Then we use Eq. (50), which for the case q, q ′ ≥ l + 1 does not contain any singularities. The case of q, q ′ ≤ n requires the usage of the Laurent series expansions derived in the previous subsection III B, as it contains different divergences depending on the values of q and q ′ .
Whenever q + q ′ + i + j ≤ 0 we get B divergences, whenever j + q ≤ 0 and i + q ′ ≥ 0 we get A divergences and for j + q ≤ 0 and i + q ′ ≤ 0, we get the mixed divergence C. Separating the sums over i and j in those three cases we get: where In section V we show that both h n,l (q, x) and d n,l (p, r, x, y) vanish for every integer value of 0 ≤ q ≤ l and 0 ≤ p ≤ l respectively. We now redefine the functions u b a into the functions f b a to be the constant term, i.e., ∼ δ 0 of its Laurent series, given in equations (44), (45), (46) and (47), completing the proof of our main result.

D. Derivation for n > l
In this case of n > l, we split the Green function into two parts: where G (sg) nl contains all terms that are singular around the origin and G (nsg) nl those that are not. In terms of explicit powers of x and y, we then get: where Ei(x) is the exponential integral function [9] and the constants A i1,i2 and B i are defined in Eqs. (21)-(22). In the case of q > l and q ′ > l, we can use Eq. (50) to immediately evaluate the elements of G (sg) nl , as: Before we can evaluate G (nsg) nl we need expressions for the terms containing logarithm and exponential integral functions. The term containing the logarithmic function is straightforward: but the exponential integral function is non-trivial: Therefore, we first evaluate a simple case: and take derivatives, with respect to λ and λ ′ : where in the last step we have made use of Eqs. (34) and (35). This gives us the integrals of the exponential integral function as: The usage of contiguous relations for the hypergeometric function introduced the indeterminacy into the expression (66) when µ = λ ′ . However, Eqs. (63)-(64) are well defined. Therefore, we start from Eqs. (63)-(64) in which we plug-in λ/λ ′ for the first argument, one for the second and apply the contiguous relations again. This yields Here H q is the harmonic number. Finally, we evaluate elements of G (nsg) nl that come out as: finalizing the derivation for the case of q, q ′ > l. Whenever q ≤ l or q ′ ≤ l the expression (59) for G (sg) nl (x, y) contains singularities and cannot be evaluated directly. To calculate this case we again employ the Laurent series of u b+δ a+δ and then show that the divergent part, that is terms with negative powers of δ, vanishes.
For this we firstly consider the case of l = 0, when we get: and by the use of Eq. (44): in which the divergent part vanishes by the definition of A.
For the case of l > 0, we get: and expanding the divergent part with the help of Eq. (44): we can perform the sum over i 1 and isolate the sum over i 2 to get: nl (x, y)x q+δ y q ′ +δ dxdy In section V we show that c n,l (q, x) vanishes for every integer value of 0 ≤ q ≤ l, provided that n > l, proving that (59) indeed converges.
To integrate (58) we first need: ∞ 0 e −λy y q ln(xy)dy = q! λ q+1 ln as well as: which is non-trivial. Therefore, we first evaluate Eq. (79) when q = 0, µ = 1 and then differentiate with respect to λ. For this we note that: and using the formula for derivatives of 1 F 1 , we get that for i > 0: Finally, we can express: The generalization for the case µ = 1 is performed via a change of variables ∞ 0 y q e −λy Ei[µ min(x, y)]dy = µ −1−q I q λ µ , xµ .
For the case of l ≥ n, we can immediately integrate the expression for G nl (x, y), to get: For the case of q non-positive integer, we can write the Laurent series of F −n , as: