A generating integral for the matrix elements of the Coulomb Green's function with the Coulomb wave functions

TL;DR

This paper analytically evaluates integrals involving the Coulomb Green's function, enabling precise calculations of energy spectra and physical properties in multi-electron atomic systems.

Contribution

It provides explicit formulas for generating integrals and moments of the Coulomb Green's function for all convergent parameter values, advancing analytical methods in quantum physics.

Findings

Derived explicit expressions for integrals involving Coulomb Green's function.

Applicable to second-order perturbation theory in atomic physics.

Facilitates analytical computation of atomic energy levels and densities.

Abstract

We analytically evaluate the generating integral $K_{nl}(\beta,\beta') = \int_{0}^{\infty}\int_{0}^{\infty} e^{-\beta r - \beta' r'}G_{nl} r^{q} r'^{q'} dr dr'$ and integral moments $J_{nl}(\beta, r) = \int_{0}^{\infty} dr' G_{nl}(r,r') r'^{q} e^{-\beta r'}$ for the reduced Coulomb Green's function $G_{nl}(r,r')$ for all values of the parameters $q$ and $q'$, when the integrals are convergent. These results can be used in second-order perturbation theory to analytically obtain the complete energy spectra and local physical characteristics such as electronic densities of multi-electron atoms or ions.

Tables

Table 1. Table 1: Comparison of computational times of the generating integral for different values of the parameters n 𝑛 n , l 𝑙 l , q 𝑞 q and q ′ superscript 𝑞 ′ q^{\prime} of analytical expressions ( 20 ) and ( 71 ) with the direct numerical integration using Mathematica 16 , when the energy is shifted by i ​ δ i 𝛿 \mathrm{i}\delta from its resonant value. The last column shows the value of δ 𝛿 \delta 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 δ 𝛿 \delta .

$n$	$l$	$q$	$q^{\prime }$	Analytic time	Numeric time	$\delta$
3	1	2	0	8$\times {10}^{-3}$ s	3.83 s	${10}^{-4}$
7	5	4	1	0.128 s	2.39 s	${10}^{-4}$
16	15	10	10	0.188 s	NaN	NaN
37	1	1	0	8.38 s	39.6 s	${10}^{-7}$

Table 2. Table 2: Comparison of computational times of the integral moments for different values of parameters n 𝑛 n , l 𝑙 l , q 𝑞 q and λ 𝜆 \lambda of analytical expressions with the direct numerical integration using Mathematica 16 at 100 different values of x 𝑥 x . The last column shows the value of δ 𝛿 \delta required to obtain the accurate results to four significant figures.

$n$	$l$	$q$	$\lambda$	Analytic time	Numeric time	$\delta$
1	0	0	0.37	0.19 s	21.7 s	${10}^{-4}$
2	4	3	0.37	0.01 s	32.2 s	${10}^{-4}$
5	4	7	0.37	3.89 s	68.1 s	${10}^{-5}$
6	8	8	0.37	0.05 s	53.5 s	${10}^{-6}$

Table 3. Table 3: Comparison of computational times of second order single-electron correction to ground state energies of some example neutral atoms.

Element	Atomic number	Analytic time	Numeric time	$\delta$
Li	3	3.26 s	29.3 s	${10}^{-3}$
F	9	14.2 s	261 s	${10}^{-4}$
Ne	10	29.3 s	321 s	${10}^{-4}$

Equations

$$K_{nl}(\beta ,{\beta }^{\prime })={\int }_0^{\infty }{\int }_0^{\infty }{e}^{-\beta r-{\beta }^{\prime }{r}^{\prime }}{G}_{nl}(r,r^{\prime })r^q{r^{\prime }}^{q^{\prime }}drd{r}^{\prime },$$

$$J_{nl}(\beta ,r)={\int }_0^{\infty }{G}_{nl}(r,r^{\prime })r^{\prime q}{e}^{-\beta r^{\prime }}d{r}^{\prime }$$

$${\int }_0^{\infty }d{r}^{\prime }{G}_{nl}(r,r^{\prime })r^{\prime k+2}{e}^{-Z{r}^{\prime }/n}$$

$$\left[-\frac{1}{2}{\nabla }^2-\frac{Z}{r}-E\right]\psi (\mathbf{r})=0.$$

$${\psi }_{nlm}(\mathbf{r})=R_{nl}(r)Y_{lm}(\Omega ),$$

$$R_{nlZ}(r)$$

$$t$$

$$R_{nlZ}(r)=Z^{3/2}{R}_{nl}(rZ)$$

$$N_{n,l}(i)$$

$$\left[-\frac{1}{2}{\nabla }^2-\frac{Z}{r}-E\right]G(\mathbf{r},{\mathbf{r}}^{\prime };E)=-\delta (\mathbf{r}-{\mathbf{r}}^{\prime })$$

$$G(\mathbf{r},{\mathbf{r}}^{\prime };E)$$

$$=\sum _{l=0}^{\infty }\sum _{m=-l}^l{Y}_{lm}(\Omega )Y_{lm}^{\ast }({\Omega }^{\prime })\text{sideset}{}^{\prime }{\sum }_{n^{\prime }=l+1}^{\infty }\frac{R_{n^{\prime }l}(r)R_{n^{\prime }l}^{\ast }(r^{\prime })}{E-E_{n^{\prime }}}=\sum _{l=0}^{\infty }\sum _{m=-l}^l{Y}_{lm}(\Omega )Y_{lm}^{\ast }({\Omega }^{\prime })G_l(r,r^{\prime };E).$$

$$(H_l-E)G_l(r,r^{\prime };E)=\left[-\frac{1}{2r}\frac{{\partial }^2}{\partial r^2}r+\frac{l(l+1)}{2r^2}-\frac{Z}{r}-E\right]G_l(r,r^{\prime };E)=-\frac{\delta (r-r^{\prime })}{r{r}^{\prime }}$$

$$G_l(r,r^{\prime };E)$$

$$t_{<}$$

$$G_n(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$=\sum _{l=0}^{\infty }\sum _{m=-l}^l{Y}_{lm}(\Omega )Y_{lm}^{\ast }({\Omega }^{\prime })\text{sideset}{}^{\ast }{\sum }_{\begin{array}{c}{n}^{\prime }=l+1\\ n^{\prime }\ne n\end{array}}^{\infty }\frac{R_{n^{\prime }l}(r)R_{n^{\prime }l}^{\ast }(r^{\prime })}{E_n-E_{n^{\prime }}}=\sum _{l=0}^{\infty }\sum _{m=-l}^l{Y}_{lm}(\Omega )Y_{lm}^{\ast }({\Omega }^{\prime })G_{nl}(r,r^{\prime }).$$

$$(H_l-E_n)G_{nl}(r,r^{\prime })=-\frac{\delta (r-r^{\prime })}{r{r}^{\prime }},$$

$$G_{nl}(r,r^{\prime })$$

$$\times \left[e^{t_{<}}\sum _{j=0}^{l-n}\left(\genfrac{}{}{0ex}{}{2l-j}{l+n}\right)\frac{(-t_{<}{)}^j}{j!}-\sum _{k=0}^{l+n}\left(\genfrac{}{}{0ex}{}{2l-k}{l-n}\right)\frac{t_{<}^k}{k!}\right],$$

$$(H_l-E_n)G_{nl}(r,r^{\prime })=R_{nl}(r^{\prime })R_{nl}(r)-\frac{\delta (r-r^{\prime })}{r{r}^{\prime }}.$$

$$G_{nl}(r,r^{\prime })$$

$$\times \left[\ln t+\ln t^{\prime }+\frac{t+t^{\prime }}{2n}-\psi (n-l)-\psi (n+l+1)-\frac{4l+5}{2n}-\mathrm{Ei}(t_{<})\right]$$

$$+L_{n-l-1}^{2l+1}(t)\left[\frac{t^{\prime }}{n}{L}_{n-l-2}^{2l+2}(t^{\prime })+\sum _{k=0}^{n-l-2}{A}_k{t}^{\prime k}\right]+L_{n-l-1}^{2l+1}(t^{\prime })\left[\frac{t}{n}{L}_{n-l-2}^{2l+2}(t)+\sum _{k=0}^{n-l-2}{A}_k{t}^k\right]$$

$$+L_{n-l-1}^{2l+1}(t_{>})e^{t_{<}}{\Phi }_{nl}(t_{<})-L_{n-l-1}^{2l+1}(t_{<})\sum _{k=1}^{2l+1}\left(\genfrac{}{}{0ex}{}{n+l}{n-l-1+k}\right)\frac{(k-1)!}{t_{>}^k}$$

$$\displaystyle+L_{n-l-1}^{2l+1}(t_{>})\sum_{k=1}^{2l+1}\frac{(k-1)!}{t_{<}^{k}}\left[{{n+l-k}\choose{n-l-1}}e^{t_{<}}-{{n+l}\choose{n-l-1+k}}\right]\Bigg{\}},$$

$$A_k$$

$${\Phi }_{nl}(x)$$

$$K_{nl}(\lambda ,{\lambda }^{\prime })=\int e^{-\lambda x-{\lambda }^{\prime }{x}^{\prime }}{G}_{nl}(x,x^{\prime })x^q{x}^{\prime q^{\prime }}dxd{x}^{\prime }$$

$$\displaystyle=2Z\sum_{i_{1}=l}^{n-1}\sum_{i_{2}=1-l}^{1+l}{{n-l-1}\choose{n-i_{1}-1}}(-1)^{i_{1}+l+1}\frac{(i_{2}-1+l)!}{(i_{1}+l+1)!}\Big{(}\frac{2}{n}\Big{)}^{i_{1}-i_{2}+1}$$

$$\displaystyle\mspace{50.0mu}\times\Bigg{\{}{{n-i_{2}}\choose{n-l-1}}\left[f_{q-i_{2}+1}^{q^{\prime}+i_{1}+1}\left(\alpha^{\prime},\alpha-\frac{2}{n}\right)+f_{q^{\prime}-i_{2}+1}^{q+i_{1}+1}\left(\alpha,\alpha^{\prime}-\frac{2}{n}\right)\right]$$

$$\displaystyle\mspace{190.0mu}-{{n+l}\choose{n-1+i_{2}}}\left[{g}_{q-i_{2}+1}^{q^{\prime}+i_{1}+1}(\alpha,\alpha^{\prime})+{g}_{q^{\prime}-i_{2}+1}^{q+i_{1}+1}(\alpha^{\prime},\alpha)\right]\Bigg{\}}$$

$$\text{mspace}30.0mu+2Z\sum _{i_1=l+1}^n\sum _{i_2=l+1}^n\left(\genfrac{}{}{0ex}{}{n-l-1}{n-i_1}\right)\left(\genfrac{}{}{0ex}{}{n+l}{n-i_2}\right)\frac{(-2/n{)}^{i_1+i_2-1}}{(i_1+l)!(i_2-l-1)!}$$

$$\displaystyle\mspace{70.0mu}\times\Bigg{\{}\frac{(q+i_{1}-1)!}{\alpha^{q+i_{1}}}\frac{(q^{\prime}+i_{2}-1)!}{\alpha^{q^{\prime}+i_{2}}}\Bigg{(}\Psi(q+i_{1})+\Psi(q^{\prime}+i_{2})-\ln\left(\frac{n^{2}}{4}\alpha\alpha^{\prime}\right)+A_{i_{1},i_{2}}$$

$$\displaystyle\mspace{250.0mu}+\frac{(2n+l-i_{1}+1)(q+i_{1})}{(i_{1}+l+1)\alpha n^{2}}+\frac{(2n+l-i_{2}+1)(q+i_{2})}{(i_{2}+l+1)\alpha^{\prime}n^{2}}\Bigg{)}$$

$$\displaystyle\mspace{70.0mu}-\left(f_{q+i_{1}}^{q^{\prime}+i_{2}}\left(\alpha^{\prime},\alpha-\frac{2}{n}\right)+f_{q^{\prime}+i_{1}}^{q+i_{2}}\left(\alpha,\alpha^{\prime}-\frac{2}{n}\right)\right)B_{i_{2}}+\left(\frac{n}{2}\right)^{q+q^{\prime}+i_{1}+i_{2}}I_{q+i_{1}-1}^{q^{\prime}+i_{2}-1}\left(\frac{n}{2}\alpha,\frac{n}{2}\alpha^{\prime}\right)\Bigg{\}},$$

$$A_{i_1,i_2}$$

$$\text{mspace}30.0mu-\Psi (1+i_1+l)-\Psi (i_2-l)-\Psi (1+i_2+l),$$

$$B_i$$

$$=\sum _{k=1}^{n-i}\frac{(i-l-1)!(i+l)!(n-i)!(-1{)}^k(k-1)!}{(k+i-l-1)!(k+i+l)!(n-i-k)!},$$

$$f_a^b(x,y)=\frac{\Gamma (a+b)}{a{x}^b}{}_2{F}_1\left(a,a+b,a+1,-\frac{y}{x}\right),$$

Full text

atomic perturbation theory, Coulomb Green function, matrix elements, effective charge

I Introduction

It has recently been demonstrated 1 that a multi-electron atom can be effectively described via a 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 1; 2; *hameka_use_1967; 4; 5; 6; 7 (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}(\bm{r},\bm{r}^{\prime})$ with the hydrogen like wave functions.

Due to the nature of the Coulomb field the radial and angular variables are separated, which allows one to integrate out the angular parts in matrix elements. The remaining radial integrals are further reduced to the computation of the generating integral from RCGF of a form

[TABLE]

which is convergent for $q\geq 0$ and $q^{\prime}\geq 0$. In Eq. (1) the letter $l$ denotes the index of angular expansion.

In principle, one can use the direct numerical integration of Eq. (1). However, if the matrix elements can be evaluated analytically then the model will provide analytical expressions for the observable characteristics of multi-electron atoms with Hartree-Fock accuracy 1. Moreover, when the indices $n$ and $l$ in Eq. (1) are large, which is the typical situation for Rydberg states 8 the direct numerical integration of Eq. (1) becomes very inefficient or even impossible due to the increasing number of nodes of the integrand. Furthermore, 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}(\beta,\beta^{\prime})$ and integral moments

[TABLE]

for all values of the parameters $q$ and $q^{\prime}$ when the integrals are convergent.

A starting point for the evaluation of (1) and (2) is the work of Johnson and Hirschfelder 9, who derived the explicit expressions for the RCGF in terms of elementary functions and computed integral moments

[TABLE]

of the RCGF, where $k\geq-l-2$. This expression is a special case of Eq. (2) when $\beta=Z/n$.

In addition, we mention the work of Hill and Huxtable 10, who evaluated the generating integral (1) for the case of $q,q^{\prime}=l+1$ and provided recurrence relations, which allow one to compute the generating integral with the increasing powers of $r,r^{\prime}$ starting from $q,q^{\prime}=l+1$. Their result is based on treating the generating integral as a Laplace transform of $\beta$ and $\beta^{\prime}$ and solving the differential equation for this Laplace transform. However, the generating integral in the range $0\leq q,q^{\prime}\leq l$ has not been evaluated. Consequently, here we evaluate the generating integral and integral moments for this case. In addition, we provide a closed form result for values $q,q^{\prime}>l+1$, which is directly applicable without the use of recurrence relations. Finally, we extend $J_{nl}(\beta,r)$ to values of $\beta$ when the integral is convergent.

In our work we follow the approach of the direct integration of the RCGF with the corresponding powers of $r,r^{\prime}$ and exponentials. The main complications come from the fact, that when $0\leq q,q^{\prime}<l+1$, the individual terms of the RCGF possess singularities that are explicitly cancelled only upon summing all expressions together.

The article is organized in the following way. For the readers convenience in Sec. I.2 we summarize the main results without any derivations. In Sec. I.3 we compare the evaluation time based on our analytical expressions with the direct numerical evaluation of the integrals. In Secs. II-IV we provide the details of all derivations. Finally, a Mathematica notebook has been prepared as a supplementary information, where we have programmed the main results of the article.

I.1 Definitions

The Hydrogen-like wave functions $\psi(\bm{r})$ satisfy the Schrödinger equation

[TABLE]

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 11. Consequently, the eigenfunctions of a discrete spectrum are represented as a product

[TABLE]

where

[TABLE]

with $L_{n-l-1}^{2l+1}(t)$ being the associated Laguerre polynomials 12 and $Y_{lm}(\Omega)$ the spherical harmonics 11. An analogous expression can be obtained for the continuous spectrum as well 11.

The associated Laguerre polynomials can be expanded according to their definition, thus leading to

[TABLE]

where the binomial is defined as ${{n}\choose{k}}=n!/(k!(n-k)!)$.

In Eqs. (4)–(8) $n$ is an integer, $n>0$, $l$ is an integer, $0\leq l\leq n-1$ and $m$ is an integer, $-l\leq m\leq 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 9

[TABLE]

and can be expanded over the eigenfunctions of Eq. (4)

[TABLE]

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^{\prime};E)$ satisfies the equation 9

[TABLE]

and is expressed in terms of the Whittaker functions 13

[TABLE]

with $\Gamma(t)$ being the gamma function 12.

The RCGF is defined as a CGF from which a state with the principal quantum number $n$ is subtracted 4

[TABLE]

The summation in RCGF $G_{nl}(r,r^{\prime})$ starts from $l+1$. Therefore, one needs to distinguish the two distinct cases. When $l\leq n-1$ the term with $n^{\prime}=n$ must be explicitly excluded. For the opposite situation of $l\geq n$ no such term appears anyway, thus making the restriction $n\neq n^{\prime}$ redundant.

As shown in 9, for the case of $l\geq n$ the RCGF, satisfying the equation

[TABLE]

is given by ($t=2Zr/n$, $t^{\prime}=2Zr^{\prime}/n$):

[TABLE]

while for $l\leq n-1$, the RCGF satisfying an equation

[TABLE]

reads:

[TABLE]

where

[TABLE]

I.2 The main result

I.2.1 Generating integrals

When $l+1\leq n$ the closed form main result is given via111We work in the rescaled variables, i.e., $x=Zr$, $x^{\prime}=Zr^{\prime}$ and $\lambda=\beta/Z$, $\lambda^{\prime}=\beta^{\prime}/Z$. Therefore, the factor from the Jacobian $Z^{-q-q^{\prime}-2}$ should be included when relating $K_{nl}(\beta,\beta^{\prime})$ and $K_{nl}(\lambda,\lambda^{\prime})$.

[TABLE]

where $\alpha=\lambda+\frac{1}{n}$ and $\alpha^{\prime}=\lambda^{\prime}+\frac{1}{n}$ and the constants $A_{i_{1},i_{2}}$ and $B_{i}$ are given by:

[TABLE]

where $\Psi(n)$ is the logarithmic derivative of the Gamma function 12 (also known as the zeroth polygamma function) and ${}_{p}F_{q}(a_{1},\ldots,a_{p};b_{1},\ldots,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_{a}^{b}(x,y)$ is defined as follows:

[TABLE]

• •

$a>0$ and $a+b>0$

[TABLE]

• •

$a=-n$ is a non-positive integer and $a+b>0$:

[TABLE]

• •

$a+b=-m$ is a non-positive integer and $a>0$

[TABLE]

• •

$a+b=-m$ is a non-positive integer and $b>0$

[TABLE]

• •

$a=-m$ and $b=-n$ are both non-positive integers

[TABLE]

When $y=0$ the function $f_{a}^{b}(x,0)$ becomes:

[TABLE]

• •

$a+b>0$ and $a\neq 0$:

[TABLE]

• •

$b>0$ and $a=0$

[TABLE]

• •

$a+b=-m$ is a non-positive integer and $a\neq 0$:

[TABLE]

• •

$b=-m$ is a negative integer and $a=0$

[TABLE]

The function ${g}_{a}^{b}(x,y)$ is defined as follows:

[TABLE]

• •

$a>0$ and $b>0$

[TABLE]

• •

$a=-n$ is a non-positive integer

[TABLE]

The function $I_{m}^{n}$ is defined as:

[TABLE]

In addition, in Eqs. (20)–(25b) ${{}_{3}\tilde{F}_{2}(a,b,c;d,e;z)}$ is the generalized regularized hypergeometric function 12 and in Eq. (26b) $H_{n}$ is a harmonic number: $H_{n}=1+\frac{1}{2}+\frac{1}{3}+...+\frac{1}{n}$.

The closed form result for $l\geq n$ reads

[TABLE]

I.2.2 Integral moments

The main result for the integral moments for the case $n\geq l+1$ reads ($x=Zr$, $\lambda=\beta/Z$)

[TABLE]

where the function $\mathcal{F}$ is defined as:

[TABLE]

In addition, the function $\mathcal{I}$ is defined as:

[TABLE]

while the function $\mathcal{G}$ is:

[TABLE]

The closed form result for the case $n\leq l$ is given via ($x=Zr$, $\lambda=\beta/Z$)

[TABLE]

I.3 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 (38) 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) (See also supplementary information). Since this expression is divergent, we introduce a parameter $\delta$ and substitute the energy as $E_{n}-\mathrm{i}\delta$, effectively regularizing the numerical integration. That is, we evaluate the integrals from the expression $G_{l}(r,r^{\prime};E_{n}-\mathrm{i}\delta)-R_{nl}(r)R_{nl}^{*}(r^{\prime})/(-\mathrm{i}\delta)$. Since the resulting answer for the generating integral is $\delta$ dependent we always ensure that the numerical value is independent of the $\delta$ by varying the parameter $\delta$. For all comparisons we employed the standard desktop Intel 2600k 3.4GHz processor 15.

In Tab. 1 we provide the summary of the evaluation times of the generating integral for different values of parameters $n$, $l$, $q$ and $q^{\prime}$. As can be observed from the table the numerical evaluation is several orders of magnitude slower then our analytical expressions. As follows from Tab. 1 the evaluation time increases significantly when $n\gg 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. 2 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 $\lambda$. 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 1. Therefore, we compare a few cases of evaluation times between our analytical and numerical approaches, which is given in Tab. 3. Typically a fully sequential version of an effz program from Ref. 1 implemented in Mathematica requires one order of magnitude larger times for the direct numerical evaluation as compared to the analytical expressions.

As a final test (see supplementary information), we evaluated generating integrals for the large number of input parameters, which corresponds to the repeated evaluation of integrals. We considered 81 examples of $K_{nl}(\lambda,\lambda^{\prime})$ and $J_{nl}(\lambda,x)$ each ($n\in\langle 2,4\rangle$, $l\in\langle 0,2\rangle$, $q\in\langle 0,2\rangle$, $\lambda\in\langle 1,2\rangle$) for a total compuaion time of 1.43 s and 0.192 s respectively, as compared to 348 s and 4.62 s numerically.

Evaluation time of numerical integrals of the RCGF depends primarily on the principal quantum number $n$ due to the oscilatory behaviour of Rydberg states 8. This 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 $\delta$ 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.

On the other hand, evaluation time of our analytical results is harder to investigate as it depends strongly on the methods used for computation of hypergeometric functions, appearing in the main results. This is further complicated by the fact that hypergeometric functions with integer coefficients can be in most cases expressed by elementary functions. This suggests that optimization of associated algorithms could reduce the evaluation time of the analytical calculaion, which is however beyond the scope of this paper. Nevertheless, 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.

II Derivation of the main results

We want to evaluate the following integrals:

[TABLE]

The main difficulty arises from the fact that the integrals of the individual terms of the RCGF contain divergences in the case when $q$ or $q^{\prime}$ are integers smaller than $l$ and these divergences cancel out upon summing all terms together. For this reason, the strategy employed in our work is to first derive expressions valid for non-integer values of $q$ and $q^{\prime}$, then find the Laurent series of (33) and (34) around $q=m+\delta$, $q^{\prime}=m^{\prime}+\delta$, where $m$ and $m^{\prime}$ are non-negative integers and finally show that for $q,q^{\prime}\geq 0$ the divergent parts (terms proportional to $\delta^{-1}$ and $\delta^{-2}$) always vanish.

II.1 Evaluation of the auxiliary integrals

First, we evaluate the following auxiliary integral

[TABLE]

convergent whenever $\operatorname{Re}a>0$, $\operatorname{Re}(a+b)>0$ and $\operatorname{Re}(\lambda+\lambda^{\prime})>0$. Here ${{}_{2}F_{1}}(a,b,c,z)$ is the Gauss Hypergeometric function 12. Let us introduce the following notation

[TABLE]

In the case when $a$ and $b$ are positive integers the expression (36) simplifies if one uses the contingous relations 12 for ${{}_{2}F_{1}}(a,b,c,z)$. This allows us to rewrite it as a finite sum

[TABLE]

Furthermore, for $b=0$ and $a$ a positive integer we get:

[TABLE]

Eq. (36) is useful when working with modern computer algebra software such as Mathematica 16 that can automatically simplify ${{}_{2}F_{1}}(a,b,c,z)$ for given integer parameters, while (37) and (38) are useful for evaluating $u_{a}^{b}$ with finite precision arithmetic, as they require evaluating an explicitly finite amount of terms.

Finally, the introduced functions $u_{a}^{b}$ allow us to evaluate the following integral:

[TABLE]

II.2 Derivation for $l\geq n$

In this case the derivation is straightforward. We employ the definition of the RCGF (15) and shift the index of summation to get:

[TABLE]

Then we simply use Eq. (39) to write the answer in terms of $u^{a}_{b}(x,y)$:

[TABLE]

II.3 Derivation for $n>l$

In the case of $n>l$, we split the Green’s function into two parts:

[TABLE]

where $G^{(\mathrm{sg})}_{nl}$ contains all terms that are singular around the origin and $G^{(\mathrm{nsg})}_{nl}$ those that are not. In terms of explicit powers of $x$ and $y$, we then get:

[TABLE]

where $\mathrm{Ei}(x)$ is the exponential integral function 12 and the constants $A_{i_{1},i_{2}}$ and $B_{i}$ are defined in Eqs. (21)-(22).

We can use Eq. (39) to immediately evaluate the elements of $G^{(\mathrm{sg})}_{nl}$, as:

[TABLE]

Before we can evaluate $G^{(\mathrm{nsg})}_{nl}$ we need expressions for the terms containing logarithm and exponential integral functions. The evaluation of the term containing the logarithmic function is straightforward:

[TABLE]

but the exponential integral function is non-trivial:

[TABLE]

Therefore, we first evaluate a simple case:

[TABLE]

and take derivatives, with respect to $\lambda$ and $\lambda^{\prime}$:

[TABLE]

where in the last step we have made use of Eqs. (37) and (38). This gives us the integrals of the exponential integral function as:

[TABLE]

The usage of contiguous relations for the hypergeometric function introduced the indeterminacy into the expression (52) when $\mu=\lambda^{\prime}$. However, Eqs. (49)-(50) are well defined. Therefore, we start from Eqs. (49)-(50) in which we plug-in $\lambda/\lambda^{\prime}$ for the first argument, one for the second and apply the contiguous relations again. This yields

[TABLE]

Here $H_{q}$ is the harmonic number.

Finally, we evaluate elements of $G^{(\mathrm{nsg})}_{nl}$ that come out as:

[TABLE]

finalizing the derivation.

III Derivation of the integral moments

We start from the evaluation of the following auxiliary integral

[TABLE]

where for the purpose of dealing with singularities we define:

[TABLE]

Integrating (43), gives:

[TABLE]

To integrate (44) we first need:

[TABLE]

as well as:

[TABLE]

which is non-trivial. Therefore, we first evaluate Eq. (59) when $q=0$, $\mu=1$

[TABLE]

and then differentiate with respect to $\lambda$.

For this we note that:

[TABLE]

and using the formula for derivatives of ${{}_{1}F_{1}}$, we get that for $i>0$:

[TABLE]

Finally, we can express:

[TABLE]

The generalization for the case $\mu\neq 1$ is performed via a change of variables

[TABLE]

As in the previous section, when $\mu=\lambda$ there is an indeterminacy in the expression. In order to handle this case, we need to take the limit of $\lambda\to 1$ in the above expression (64), to get:

[TABLE]

Finally we integrate (44) to get:

[TABLE]

For the case of $l\geq n$, we can immediately integrate the expression for $G_{nl}(x,y)$, to get:

[TABLE]

III.1 Laurent series of $u_{a}^{b}$

In order to integrate the Green function with $q,q^{\prime}<l$ we will need an expansion of $u_{a}^{b}$ functions in the Laurent series in the neighbourhood of the singularities when $a=-m+\delta$ and $a+b=-n+2\delta$. In order to obtain the Laurent series it is most convenient to start from the infinite series representation (36), and separate singular terms. For this we split the series into two parts, namely $i=0..n-1$ and $n..\infty$. The latter part is finite and is reduced to:

[TABLE]

where ${{}_{3}\tilde{F}_{2}}(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:

[TABLE]

while for negative:

[TABLE]

Other relevant functions are

[TABLE]

and

[TABLE]

In addition, we use the following expansions

[TABLE]

This allows us to find the Laurent series of $u_{a}^{b}$ 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$):

[TABLE]

• •

When $a+b$ is a non-positive integer and $a$ is positive, the first $-a-b$ terms diverge:

[TABLE]

• •

When $a+b$ is a non-positive integer and $b$ is positive, the first $-a-b$ temrs diverge as $\Gamma(-m)$ and independently the term with $i=-a$ diverges as $1/(a-a)$:

[TABLE]

• •

When $a$ and $b$ are both non-positive integers, the first $-a-b$ terms diverge as $\Gamma(-m-n)$ with the $i=-a$ term diverging faster, as $\Gamma(-m-n)/(a-a)$:

[TABLE]

Lastly, examining (45) it is clear that we also need to find the Laurent series representation of the function $u_{a}^{\prime b}$ defined as

[TABLE]

In the case of $a$ being a non-positive integer and $b>0$, we get:

[TABLE]

III.2 Laurent series of $K$

We first consider the case of $l\geq n$. Whenever $q+q^{\prime}+i+j\leq 0$ we get $\mathcal{B}$ divergences, whenever $j+q\leq 0$ and $i+q^{\prime}\geq 0$ we get $\mathcal{A}$ divergences and for $j+q\leq 0$ and $i+q^{\prime}\leq 0$, we get the mixed divergence $\mathcal{C}$. Separating the sums over $i$ and $j$ in those three cases we get:

[TABLE]

where

[TABLE]

and

[TABLE]

where

[TABLE]

In section IV we show that both $h_{n,l}(q,x)$ and $d_{n,l}(p,r,x,y)$ vanish for every integer value of $0\leq q\leq l$ and $0\leq p\leq l$ respectively, proving that (41) converges in these cases.

For the case of $n>l$, we only need to consider $G_{nl}^{(\mathrm{sg})}(x,y)$, as $G_{nl}^{(\mathrm{nsg})}(x,y)$ never contains any singularities. In the case of $l=0$, we get:

[TABLE]

and by the use of Eq. (76):

[TABLE]

in which the divergent part vanishes by the definition of $\mathcal{A}$.

For the case of $l>0$, we get:

[TABLE]

and expanding the divergent part with the help of Eq. (76):

[TABLE]

we can perform the sum over $i_{1}$ and isolate the sum over $i_{2}$ to get:

[TABLE]

where

[TABLE]

In section IV we show that $c_{n,l}(q,x)$ vanishes for every integer value of $0\leq q\leq l$, provided that $n>l$, proving that (45) indeed converges.

We now define the functions $f_{a}^{b}$ to be the constant term, i.e., $\sim\delta^{0}$ of the Laurent series of $u_{a}^{b}$, given in equations (76), (77), (78) and (79), thus arriving at our main result for integer values of $q$ and $q^{\prime}$.

III.3 Laurent series of J

For the case of non-positive integer, we can write the Laurent series of $\mathcal{F}_{-n}$, as:

[TABLE]

which gives the Laurent series in the case of $n>l$, as:

[TABLE]

where $M$ is the first Whittaker function.

On the other hand, the case of $l\geq n$ comes out as:

[TABLE]

Since $c_{n,l}(q,x)$ and $h_{n,l}(q,x)$ both vanish for every integer value of $0\leq q\leq l$, we can redefine $\mathcal{F}_{q}(x,y)$ to the constant term of it’s Laurent series completing the derivation of $J$.

IV Proofs

IV.1 $c_{n,l}=0$

Here we show that $c_{n,l}$ vanishes for all integers $q\leq l$ in the case of $n>l$, while $h_{n,l}$ vanishes for all $q\leq l$ in the case $l\geq n$.

The function $c_{n,l}$ is given by a difference of two finite series. Let us call them $c^{(1)}_{n,l}$ and $c^{(2)}_{n,l}$. Using binomial theorem we get:

[TABLE]

Relabeling indices and using properties of binomials, this transforms into:

[TABLE]

where in the last step we have used 17:

[TABLE]

IV.2 $h_{n,l}=0$

Proceeding along those same lines for $h_{n,l}$, we can similarly use the binomial theorem to write:

[TABLE]

and redefine summation indices to arrive to:

[TABLE]

Finally, using properties of binomials we get:

[TABLE]

where in the last step we have used the following identity valid whenever $b<a$ 17:

[TABLE]

IV.3 $d_{n,l}=0$

Here we show that $d_{n,l}(p,r,x,y)$ vanishes for all integer values of $0\leq p\leq 2l$ and $0<r<l$, assuming that $l\geq n$.

Firstly, $d_{n,l}$ is given as a difference of two series, that we denote, as $d_{n,l}^{(1)}$ and $d_{n,l}^{(2)}$. Thus we have:

[TABLE]

We can use the binomial theorem to transform it to:

[TABLE]

Redefinition of the summation variables, by $k\rightarrow j$, $i\rightarrow k$ and $j\rightarrow i+j-k$ and changing the order of sums yield:

[TABLE]

Using properties of binomials we can write it as:

[TABLE]

Finally, we use:

[TABLE]

to arrive at:

[TABLE]

completing the proof.

V Conclusion

In our work we calculated the generating integral and integral moments of the RCGF that appear during the calculation of the matrix elements in second-order perturbation theory of multi-electron atoms and ions. Our closed form analytical results allow one to effectively compute the generating integrals and are expressed through elementary and special functions.

In contrast to previous works, we followed the approach of the direct integration of the RCGF with the corresponding powers of $r,r^{\prime}$ and exponentials. The main complications came from the fact, that when $0\leq q,q^{\prime}<l+1$, the integrals from the individual terms of the RCGF possess singularities that are explicitly cancelled only upon summing all parts of the integrals of RCGF together. For this reason, we employed the strategy of first to derive expressions valid for non-integer values of $q$ and $q^{\prime}$ and then find the Laurent series of (33) and (34) around $q=m+\delta$, $q^{\prime}=m^{\prime}+\delta$, where $m$ and $m^{\prime}$ are non-negative integers. After this we demonstrated that for $q,q^{\prime}\geq 0$ the divergent parts, that is the terms proportional to $\delta^{-1}$ and $\delta^{-2}$ always vanish.

We validated the evaluation times of the generating integrals for different input parameters and found out that the results via our analytical expressions are typically several orders of magnitude faster than the corresponding ones via direct numerical integrations. Moreover, our analytical expressions provide correct results even in cases where the direct numerical integration typically fails to converge.

Acknowledgements.

The authors are grateful to I. D. Feranchuk, C. H. Keitel, A. U. Leonau, S. Bragin, N. Oreshkina and Z. Harman for useful discussions. This article comprises parts of the PhD thesis work of Kamil Dzikowski to be submitted to the Heidelberg University, Germany.