We analytically evaluate the generating integral Knl(β,β)=00eβrβrGnl(r,r)rqrqdrdr and integral moments Jnl(β,r)=0drGnl(r,r)rqeβr for the reduced Coulomb Green’s function Gnl(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.

It has recently been demonstrated1 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 theory1–6 (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.

Due to the nature of the Coulomb field, the radial and angular variables are separated, which allows one to integrate 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

Knl(β,β)=00eβrβrGnl(r,r)rqrqdrdr,
(1)

which is convergent for q ≥ 0 and q′ ≥ 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,7 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 Knl(β, β′) and integral moments

Jnl(β,r)=0Gnl(r,r)rqeβrdr,
(2)

for all values of the parameters q and q′ when the integrals are convergent.

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

0drGnl(r,r)rk+2eZr/n,
(3)

of the RCGF, where k ≥ −l − 2. This expression is a special case of Eq. (2) when β = Z/n.

In addition, we mention the work of Hill and Huxtable,9 who evaluated the generating integral (1) for the case of q, q′ = l + 1 and provided 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 a closed form result for values q, q′ > l + 1, which is directly applicable without the use of recurrence relations. Finally, we extend Jnl(β, r) to values of β 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′ and exponentials. The main complications come from the fact, that when 0 ≤ q, q′ < l + 1, the individual terms of the RCGF possess singularities that are explicitly canceled only upon summing all expressions together.

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

The hydrogen-like wave functions ψ(r) satisfy the Schrödinger equation,

122ZrEψ(r)=0.
(4)

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

ψnlm(r)=Rnl(r)Ylm(Ω),
(5)

where

RnlZ(r)=Z3(nl1)!(n+l)!2n2tlet/2Lnl12l+1(t),t=2Zrn,
(6)

with Lnl12l+1(t) being the associated Laguerre polynomials11 and Ylm(Ω) being the spherical harmonics.10 An analogous expression can be obtained for the continuous spectrum as well.10 

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

RnlZ(r)=Z3/2Rnl(rZ)=Z3/2eZrni=l+1nNn,l(i)(rZ)i1,
(7)
Nn,l(i)=(nl1)!(n+l)!1n2ni(1)l+1(i1l)!n+lni,
(8)

where the binomial is defined as nk=n!/(k!(nk)!).

In Eqs. (4)–(8), n is an integer, n > 0, l is an integer, 0 ≤ ln – 1, and m is an integer, − lml.

The bound states of Eq. (4) are described by energy eigenvalues En = −Z2/(2n2).

Green’s function of the hydrogen-like atom satisfies the equation8 

122ZrEG(r,r;E)=δ(rr)
(9)

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

G(r,r;E)=*n=1l=0n1m=llψnlm(r)ψnlm*(r)EEn=l=0m=llYlm(Ω)Ylm*(Ω)n=l+1Rnl(r)Rnl*(r)EEn=l=0m=llYlm(Ω)Ylm*(Ω)Gl(r,r;E).
(10)

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 Gl(r, r′; E) satisfies the equation8 

(HlE)Gl(r,r;E)=12r2r2r+l(l+1)2r2ZrEGl(r,r;E)=δ(rr)rr
(11)

and is expressed in terms of the Whittaker functions,12 

Gl(r,r;E)=4ZνΓ(l+1ν)t<t>Mν,l+1/2(t<)Wν,l+1/2(t>),t<=min(t,t),t>=max(t,t),ν=Z22E,
(12)

with Γ(t) being the gamma function.11 

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

Gn(r,r)=*n=1nnl=0n1m=llψnlm(r)ψnlm*(r)EnEn=l=0m=llYlm(Ω)Ylm*(Ω)*n=l+1nnRnl(r)Rnl*(r)EnEn=l=0m=llYlm(Ω)Ylm*(Ω)Gnl(r,r).
(13)

The summation in RCGF Gnl(r, r′) starts from l + 1. Therefore, one needs to distinguish the two distinct cases. When ln – 1, the term with n′ = n must be explicitly excluded. For the opposite situation of ln, no such term appears anyway, thus making the restriction nn′ redundant.

As shown in Ref. 8, for the case of ln, the RCGF, satisfying the equation

(HlEn)Gnl(r,r)=δ(rr)rr,
(14)

is given by (t = 2Zr/n, t′ = 2Zr′/n),

Gnl(r,r)=(1)l+1n4Zn(ln)!(l+n)!(tt)l1e(t+t)/2i=0l+n2lilnt>ii!×et<j=0ln2ljl+n(t<)jj!k=0l+n2lklnt<kk!,
(15)

while for ln − 1, the RCGF satisfying an equation

(HlEn)Gnl(r,r)=Rnl(r)Rnl(r)δ(rr)rr
(16)

reads

Gnl(r,r)=4Zn(nl1)!(n+l)!e(t+t)/2Lnl12l+1(t)Lnl12l+1(t)×lnt+lnt+t+t2nψ(nl)ψ(n+l+1)4l+52nEi(t<)+Lnl12l+1(t)tnLnl22l+2(t)+k=0nl2Aktk+Lnl12l+1(t)tnLnl22l+2(t)+k=0nl2Aktk+Lnl12l+1(t>)et<Φnl(t<)Lnl12l+1(t<)k=12l+1n+lnl1+k(k1)!t>k+Lnl12l+1(t>)k=12l+1(k1)!t<kn+lknl1et<n+lnl1+k,
(17)

where

Ak=(1)kk!n+lnl1kj=k+1nl12j+2l+1j(j+2l+1),
(18)
Φnl(x)=j=1nl1(x)jj!n+lnl1jk=1j(k1)!xk.
(19)

1. Generating integrals

When l + 1 ≤ n, the closed form main result is given via13 

Knl(λ,λ)=eλxλxGnl(x,x)xqxqdxdx=2Zi1=ln1i2=1l1+lnl1ni11(1)i1+l+1(i21+l)!(i1+l+1)!2ni1i2+1  ×ni2nl1fqi2+1q+i1+1α,α2n+fqi2+1q+i1+1α,α2n        n+ln1+i2gqi2+1q+i1+1(α,α)+gqi2+1q+i1+1(α,α)+2Zi1=l+1ni2=l+1nnl1ni1n+lni2(2/n)i1+i21(i1+l)!(i2l1)!  ×(q+i11)!αq+i1(q+i21)!αq+i2Ψ(q+i1)+Ψ(q+i2)lnn24αα+Ai1,i2       +(2n+li1+1)(q+i1)(i1+l+1)αn2+(2n+li2+1)(q+i2)(i2+l+1)αn2  fq+i1q+i2α,α2n+fq+i1q+i2α,α2nBi2+n2q+q+i1+i2Iq+i11q+i21n2α,n2α,
(20)

where α=λ+1n and α=λ+1n and the constants Ai1,i2 and Bi are given by

Ai1,i2=4l+52n+Ψ(n+l+1)+Ψ(nl)Ψ(i1l)Ψ(1+i1+l)Ψ(i2l)Ψ(1+i2+l),
(21)
Bi=(ni)(i+l+1)(il)F32(1,1,1n+i,1+il,2+i+l,1)=k=1ni(il1)!(i+l)!(ni)!(1)k(k1)!(k+il1)!(k+i+l)!(nik)!,
(22)

where Ψ(n) is the logarithmic derivative of the Gamma function11 (also known as the zeroth polygamma function) and Fpq(a1,,ap;b1,,bq;z) is the generalized hypergeometric function, which in this case can be evaluated as a finite series [see Eq. (22)].

The function fab(x,y) is defined as follows:

  • a > 0 and a + b > 0:

fab(x,y)=Γ(a+b)axbF21a,a+b,a+1,yx.
(23a)
  • a = −n is a non-positive integer and a + b > 0:

fnb(x,y)=2Γ(b)xb(y)nn!Ψ(b)lnx+i=0n1Γ(bn+i)(in)xbn+i(y)ii!+b!(y)1+nxb+1F̃321,1,b+1,n+2,2,yx.
(23b)
  • a + b = −m is a non-positive integer and a > 0:

faam(x,y)=i=0m(x)mi(mi)!(i+a)(y)ii!Ψ(m+1i)lnx12(i+a)+(a+m)!(y)m+1xF̃321,1,a+m+1,m+2,2+a+m,yx.
(23c)
  • a + b = −m is a non-positive integer and b > 0:

fmbb(x,y)=i=0m(x)mi(mi)!(imb)(y)ii!Ψ(m+1i)lnx12(imb)+i=1+mm+b1Γ(im)(imb)xim(y)ii!+2Γ(b)xb(y)m+b(m+b)!(Ψ(b)lnx)+b!(y)m+b+1xb+1F̃321,1,1+b,2,2+m+b,yx.
(23d)
  • a = −m and b = −n are both non-positive integers:

fmn(x,y)=i=0,imm+n(x)m+ni(m+ni)!(im)(y)ii!Ψ(m+n+1i)ln(x)12(im)+(x)nn!(y)mm!(Ψ(n+1)ln(x))2Ψ(1)(n+1)+π23+n!(y)m+n+1xF̃321,1,1+n,m+n+2,n+2,yx.
(23e)

When y = 0, the function fab(x,0) becomes

  • a + b > 0 and a ≠ 0:

fab(x,0)=Γ(a+b)axa+b.
(24a)
  • b > 0 and a = 0:

f0b(x,0)=2Γ(b)xb(Ψ(b)lnx).
(24b)
  • a + b = −m is a non-positive integer and a ≠ 0:

faam(x,0)=(x)mam!Ψ(1+m)ln(x)12a.
(24c)
  • b = −m is a negative integer and a = 0:

f0m(x,0)=(x)mm!π23+[Ψ(1+m)ln(x)]2Ψ(1)(1+m).
(24d)

The function gab(x,y) is defined as follows:

  • a > 0 and b > 0:

gab(x,y)=Γ(a)xaΓ(b)yb.
(25a)
  • a = −n is a non-positive integer:

gnb(x,y)=Γ(b)yb(x)nn!(Ψ(1+n)+Ψ(b)ln(xy)).
(25b)

The function Imn is defined as

Imn(x,y)=n!yn+1i=0nyii!fm+1iy1,x+m!xm+1fn+10(1,y),
(26a)

and when y = 1,

Imn(x,1)=n!t=1n1t!ftm+1(x,0)+m!xm+1(Hmln(x)).
(26b)

In addition, in Eqs. (20)–(25b), F̃32(a,b,c;d,e;z) is the generalized regularized hypergeometric function,11 and in Eq. (26b), Hn is a harmonic number Hn=1+12+13++1n. The closed form result for ln reads

Knl(λ,λ)=eλxλxGnl(x,x)xqxqdxdx=(1)n+lZni=lnj=ln(li)!(l+j)!l+ni+l2ni+j×ljl+n(1)j+lfj+qi+qα,α2n+fj+qi+qα,α2nljlnfj+qi+qα,α+fj+qi+qα,α.
(27)

2. Integral moments

The main result for the integral moments for the case nl + 1 reads (x = Zr, λ = β/Z)

Jnl(λ,x)=Zq1eλyGn,l(x,y)yqdy=2Zqexni1=ln1i2=1l1+lnl1ni11(1)i1+l+1(i21+l)!(i1+l+1)!2ni1i2+1×ni2nl1xi1Fq+1i2(2/nα,x)+xi2e2x/nα1+i1+qΓ(1+i1+q,xα)n+ln1+i2xi1Gq+1i2(α)+(q+i1)!xi2α1+q+i1+i1=l+1ni2=l+1nnl1ni1(2/n)i1+i21(i1+l)!(i2l1)!×xi11Fq+i22nα,x+xi21αi1+qexnΓ(i1+q,xα)Bi2xi11(q+i21)!αq+i2ln4n2xα+Ψ(q+i2)+2n+li1+1i1+l+1xn2+2n+li2+1i2+l+1i2+qn2α+Ai1,i2xi11n2q+i2Iq+i21(αn/2,2x/n),
(28)

where the function F is defined as

Fq(x,y)=yqqF11(q,q+1,xy),
(29a)

and in the case of q being a non-positive integer,

Fn(x,y)=i=0n1xiinyini!+xnn!lny+yx1+n(1+n)!F22(1,1,2,2+n,xy),
(29b)

while in the case of x = 0,

Fq(0,y)=yqqF0(0,y)=lny.
(29c)

In addition, the function I is defined as

Iq(x,y)=q!xq+1Ei(yxy)ln(1x)+i=1q(xy)iii!F11(i,i+1,yyx),
(30a)

and for x = 1,

Iq(1,y)=q!lny+γ+i=1qyiii!,
(30b)

while the function G is

Gq(x)=Γ(q)xq,
(31a)

and for negative integers,

Gn(x)=(x)nn!(Ψ(1+n)ln(x)).
(31b)

The closed form result for the case nl is given via (x = Zr, λ = β/Z)

Jnl(λ,x)=Zq1eλyGnl(x,y)yqdy=(1)n+lZqnexni=lnj=ln(li)!(l+j)!n+li+l2ni+j×ljl+n(1)j+lxi1Fj+q(2/nα,x)+xj1αi+qe2x/nΓ(i+q,xα)ljlnxi1Fj+q(α,x)+xj1αi+qΓ(i+q,xα).
(32)

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 En = −Z2/(2n2) has been explicitly subtracted from Green’s function Eq. (12) (see also the supplementary material). Since this expression is divergent, we introduce a parameter δ and substitute the energy as En − iδ, effectively regularizing the numerical integration. That is, we evaluate the integrals from the expression Gl(r,r;Eniδ)Rnl(r)Rnl*(r)/(iδ). Since the resulting answer for the generating integral is δ dependent, we always ensure that the numerical value is independent of δ by varying the parameter δ. For all comparisons, we employed the standard desktop Intel 2600K 3.4 GHz processor.14 

In Table 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 than our analytical expressions. As follows from Table I, the evaluation time increases significantly when nl. 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.

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 (71) with the direct numerical integration using Mathematica,15 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 δ.

nlqqAnalytic time (s)Numeric time (s)δ
8 × 10−3 3.83 10−4 
0.128 2.39 10−4 
16 15 10 10 0.188 NaN NaN 
37 8.38 39.6 10−7 
nlqqAnalytic time (s)Numeric time (s)δ
8 × 10−3 3.83 10−4 
0.128 2.39 10−4 
16 15 10 10 0.188 NaN NaN 
37 8.38 39.6 10−7 

In Table 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.

TABLE II.

Comparison of computational times of the integral moments for different values of parameters n, l, q, and λ of analytical expressions with the direct numerical integration using Mathematica15 at 100 different values of x. The last column shows the value of δ required to obtain the accurate results to four significant figures.

nlqλAnalytic time (s)Numeric time (s)δ
0.37 0.19 21.7 10−4 
0.37 0.01 32.2 10−4 
0.37 3.89 68.1 10−5 
0.37 0.05 53.5 10−6 
nlqλAnalytic time (s)Numeric time (s)δ
0.37 0.19 21.7 10−4 
0.37 0.01 32.2 10−4 
0.37 3.89 68.1 10−5 
0.37 0.05 53.5 10−6 

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 Table III. 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.

TABLE III.

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

AtomicAnalyticNumeric
Element number time (s) time (s)δ
Li 3.26 29.3 10−3 
14.2 261 10−4 
Ne 10 29.3 321 10−4 
AtomicAnalyticNumeric
Element number time (s) time (s)δ
Li 3.26 29.3 10−3 
14.2 261 10−4 
Ne 10 29.3 321 10−4 

As a final test (see the supplementary material), we evaluated generating integrals for the large number of input parameters, which corresponds to the repeated evaluation of integrals. We considered 81 examples of Knl(λ, λ′) and Jnl(λ, x) each (n ∈ ⟨2, 4⟩, l ∈ ⟨0, 2⟩, q ∈ ⟨0, 2⟩, λ ∈ ⟨1, 2⟩) for a total computation 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 oscillatory behavior of Rydberg states.7 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 δ to keep the constant accuracy since in many cases, large positive values are almost completely canceled 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 calculation, 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 s on Intel 2600K 3.4 GHz processor.

We want to evaluate the following integrals:

Knl(λ,λ)=eλxλxGnl(x,x)xqxqdxdx,
(33)
Jnl(λ,x)=eλyGnl(x,y)yqdy.
(34)

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′ 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′ and then to find the Laurent series of (33) and (34) around q = m + δ, q′ = m′ + δ, where m and m′ are non-negative integers and finally show that for q, q′ ≥ 0, the divergent parts (terms proportional to δ−1 and δ−2) always vanish.

First, we evaluate the following auxiliary integral:

0reλrλrra1rb1drdr=0eλrra1λbΓ(b,λr)dr=Γ(a+b)aλa+bF21a,a+b,a+1,λλ,
(35)

which is convergent whenever Re a > 0, Re(a + b) > 0, and Re(λ + λ′) > 0. Here, F21(a,b,c,z) is the Gauss Hypergeometric function.11 Let us introduce the following notation:

uab(x,y)=Γ(a+b)axa+bF21a,a+b,a+1,yx=i=0Γ(a+b+i)(a+i)xa+b+i(y)ii!.
(36)

In the case when a and b are positive integers, the expression (36) simplifies if one uses the continuous relations11 for F21(a,b,c,z). This allows us to rewrite it as a finite sum,

uab(x,y)=Γ(a)yaΓ(b)xbx+ybi=0a1Γ(b+i)i!yy+xi.
(37)

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

ua0(x,y)=Γ(a)yaln1+yxi=1a11iyy+xi.
(38)

Equation (36) is useful when working with modern computer algebra software such as Mathematica15 that can automatically simplify F21(a,b,c,z) for given integer parameters, while (37) and (38) are useful for evaluating uab with finite precision arithmetic, as they require evaluating an explicitly finite amount of terms.

Finally, the introduced functions uab allow us to evaluate the following integral:

0eμr<λrλrrqrqr<pr>pdrdr=0re(μλ)rλrrq+prq+pdrdr+0re(μλ)rλrrq+prq+pdrdr=uq+p+1q+p+1(λ,λμ)+uq+p+1q+p+1(λ,λμ).
(39)

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

Gnl(x,y)=(1)n+lni=lnj=ln(li)!(l+j)!n+li+l2ni+jexynmax[x,y]i1min[x,y]j1×ljl+n(1)j+le2/nmin[x,y]ljln.
(40)

Then, we simply use Eq. (39) to write the answer in terms of uba(x,y),

eλxλyGnl(x,y)xqyqdxdy=(1)n+lZni=lnj=ln(li)!(l+j)!l+ni+l2ni+j×ljl+n(1)j+luj+qi+qα,α2n+uj+qi+qα,α2nljlnuj+qi+qα,α+uj+qi+qα,α.
(41)

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

Gnl(x,y)=4Zn(nl1)!(n+l)!Gnl(sg)(x,y)+Gnl(nsg)(x,y),
(42)

where Gnl(sg) contains all terms that are singular around the origin and Gnl(nsg) those that are not. In terms of explicit powers of x and y, we then get

Gnl(sg)(x,y)=exyni1=ln1i2=1l1+ln+l1ni11(1)i1+l(i21+l)!(i1l)!2ni1i2×ni2nl1e2/nmin[x,y]max[x,y]i1min[x,y]i2n+ln1+i2xi1yi2+yi1xi2,
(43)
Gnl(nsg)(x,y)=exyni1=l+1ni2=l+1nn+lni1n+lni2(1)i1+i2+1(i1l1)!(i2l1)!2ni1+i22×max[x,y]i11min[x,y]i21e2/nmin[x,y]Bi2xi11yi21ln4n2xyEi2nmin[x,y]+Ai1,i2+2n+li1+1i1+l+1xn2+2n+li2+1i2+l+1yn2,
(44)

where Ei(x) is the exponential integral function11 and the constants Ai1,i2 and Bi are defined in Eqs. (21) and (22).

We can use Eq. (39) to immediately evaluate the elements of Gnl(sg), as

eλxλyGnl(sg)(x,y)xqyqdxdy=i1=ln1i2=1l1+ln+lni11(1)i1+l(i21+l)!(i1l)!2ni1i2×ni2nl1uqi2+1q+i1+1α,α2n+uqi2+1q+i1+1α,α2nn+ln1+i2(q+i1)!α1+q+i1Γ(qi2+1)αqi2+1+Γ(qi2+1)αqi2+1(q+i1)!α1+q+i1.
(45)

Before we can evaluate Gnl(nsg), we need expressions for the terms containing logarithm and exponential integral functions. The evaluation of the term containing the logarithmic function is straightforward,

eλxλyln(μxy)xqyqdxdy=q!λq+1q!λq+1Ψ(q+1)+Ψ(q+1)lnλλμ,
(46)

but the exponential integral function is non-trivial,

Iq,q(λ,λ)=eλxλyxqyqEimin[x,y]dxdy.
(47)

Therefore, we first evaluate a simple case,

I0,0(λ,λ)=eλxλyEimin[x,y]dxdy=ln(λ+λ1)λλ,
(48)

and take derivatives, with respect to λ and λ′,

Iq,q(λ,λ)=(1)q+qλ(q)λ(q)ln(λ+λ1)λλ=(1)q+qsqtqqsqtλ(s)λ(t)ln(λ+λ1)λ(qs)λ(qt)1λλ=(1)q+qs,t,s+t0q,qqsqt(1)s+t+1Γ(s+t)(λ+λ1)s+t(1)s+q(qs)!λqs(1)q+t(qt)!λqt+(1)q(q)!λq(1)q(q)!λqlog(λ+λ1)=(q)!λq(q)!λqs=0qt=1qΓ(s+t)s!t!λsλt(λ+λ1)s+t
(49)
+s=1qΓ(s)s!λs(λ+λ1)sln(λ+λ1)
(50)
=q!λq+1i=0qλii!uq+1iλ1,λ+q!λq+1uq+10(1,λ),
(51)

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

eλxλyEiμmin[x,y]xqyqdxdy=μ2qqIq,qλμ,λμ=q!λq+1i=0qλii!uq+1iλμ,λq!λq+1uq+10(μ,λ).
(52)

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

λ2qqIq,qλλ,1=(q)!λq+1(q)!λq+1s=0qt=1qΓ(s+t)s!t!λt(λ)t+s=1qΓ(s)s!ln(λλ)=q!λq+1λq+1t=1q1t!utq+1λλ,0+q!Hqlnλλ.
(53)

Here, Hq is the harmonic number.

Finally, we evaluate elements of Gnl(nsg) that come out as

eλxλyGnl(nsg)(x,y)xqyqdxdy=i1=l+1ni2=l+1nn+lni1n+lni2(1)i1+i2+1(i1l1)!(i2l1)!2ni1+i22×uq+i1q+i2(α,α2n)+uq+i1q+i2(α,α2n)Bi2(q+i11)!λq+i1(q+i21)!λq+i2Ψ(q+i1)+Ψ(q+i2)lnn24αα+Ai1,i2+(2n+li1+1)(q+i1)(i1+l+1)λn2+(2n+li2+1)(q+i2)(i2+l+1)λn2+n2q+q+i1+i2Iq+i11,q+i21n2λ,n2λ,
(54)

finalizing the derivation.

We start from the evaluation of the following auxiliary integral:

0eμmin[r,r]λrrqmin[r,r]amax[r,r]bdr=r1+a+b+q1+a+qF11(1+a+q,2+a+q,r(μλ))+raλ1+b+qeμrΓ(1+b+q,rλ)=rbF1+a+q(μλ,r)+raλ1+b+qeμrΓ(1+b+q,rλ),
(55)

where for the purpose of dealing with singularities, we define

Fa(x,y)=yaaF11[a,a+1,xy]=i=0xia+iya+ii!.
(56)

Integrating (43), gives

0eλyyqGnl(sg)(x,y)dy=exni1=ln1i2=1l1+ln+lni11(1)i1+l(i21+l)!(i1l)!2ni1i2×ni2nl1xi1Fq+1i2(α2/n,x)+xi2α1+i1+qe2x/nΓ(1+i1+q,xα)n+ln1+i2xi1α1+qi2Γ(qi2+1)+1xi2α1+q+i1Γ(q+i1+1).
(57)

To integrate (44), we first need

0eλyyqln(xy)dy=q!λq+1lnxλ+Ψ(1+q),
(58)

as well as

Iq(λ,x)=0eλyyqEi[μmin(x,y)]dy,
(59)

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

I0(λ,x)=0eλyEi[min(x,y)]dy,
(60)

and then differentiate with respect to λ.

For this, we note that

λ(Ei(x(1λ))ln(1λ))=ex(1λ)λ1=(x)F11(1,2,x(1λ)),
(61)

and using the formula for derivatives of F11, we get that for i > 0,

λ(i)(Ei(x(1λ))log(1λ))=(x)iiF11(i,i+1,x(1λ)).
(62)

Finally, we can express

Iq(λ,x)=(1)qλ(q)I0(λ,x)=(1)qi=0qqiλ(i)(Ei(x(1λ))log(1λ))λ(qi)λ1=q!λqI0(λ,x)+i=1qq!i!xiiF11[i,i+1,x(1λ)]λ1+qi.
(63)

The generalization for the case μ ≠ 1 is performed via a change of variables

0yqeλyEi[μmin(x,y)]dy=μ1qIqλμ,xμ.
(64)

As in Sec. II C, when μ = λ, there is an indeterminacy in the expression. In order to handle this case, we need to take the limit of λ → 1 in the above expression (64) to get

Iq(1,x)=q!I0(1,x)+i=1qq!i!xii=q!lnx+γ+i=1qxiii!.
(65)

Finally, we integrate (44) to get

0eλyyqGnl(nsg)(x,y)dy=exni1=l+1ni2=l+1nn+lni1(1)i1+i2+1(i1l1)!(i2l1)!2ni1+i22×xi11Fq+i2(2/nα,x)+xi21αi1+qexnΓ(i1+q,xα)Bi2xi11(q+i21)!αq+i2ln4n2xμ+Ψ(q+i2)+2n+li1+1i1+l+1xn2+2n+li2+1i2+l+1i2+qn2α+Ai1,i2xi11(n/2)q+i2Iq+i21(αn/2,2x/n).
(66)

For the case of ln, we can immediately integrate the expression for Gnl(x, y) to get

Jnl(λ,x)=0eλyyqGnl(x,y)dy=(1)n+lni=lnj=ln(li)!(l+j)!n+li+l2ni+jexn×ljl+n(1)j+lxi1Fj+q(2/nα,x)+xj1αi+qe2/nxΓ(i+q,xα)ljlnxi1Fj+q(α,x)+xj1αi+qΓ(i+q,xα).
(67)

In order to integrate the Green function with q, q′ < l, we will need an expansion of uab functions in the Laurent series in the neighborhood of the singularities when a = −m + δ and a + b = −n + 2δ. 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. The latter part is finite and is reduced to

ũab(x,y,n)=i=nΓ(a+b+i)(a+i)xa+b+i(y)ii!=Γ(a+b+n)xa+b+nΓ(a+n)(y)nF̃321,a+n,a+b+n,1+n,1+a+n,yx,
(68)

where F̃32(a1,a2,a3,b1,b2,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

Γ(δ+n)=Γ(n)(1+δΨ(n))+O(δ2),
(69)

while for negative,

Γ(δn)=(1)nn!1δ+Ψ(n+1)+δ2Ψ(n+1)2Ψ(1)(n+1)+π23+O(δ2).
(70)

Other relevant functions are

1xδ+n=xn1δlnx+δ2ln2x2+O(δ3)
(71)

and

1x+δ=1x1δx+O(δ2).
(72)

In addition, we use the following expansions:

An(x)=Γ(n+2δ)δxn+2δ=Γ(n)xn1δ+2Ψ(n)2lnx+O(δ),
(73)
Bmn(x)=Γ(n+2δ)(m+δ)xn+2δ=(x)nn!m12δ+Ψ(n+1)lnx12m+O(δ),
(74)
Cn(x)=Γ(n+2δ)δxn+2δ=(x)nn!12δ2+Ψ(n+1)lnxδ+Ψ(n+1)2Ψ(1)(n+1)+π232ln(x)Ψ(n+1)+ln2x+O(δ).
(75)

This allows us to find the Laurent series of uab 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),

un+δb+δ(x,y)=i=0n1Γ(bn+i)(in)xb(y)ii!+Ab(x)(y)n(n)!+ũnb(x,y,n+1).
(76)
  • When a + b is a non-positive integer and a is positive, the first −ab terms diverge

uδ+aδan(x,y)=i=0nBnia+i(x)(y)ii!+ũaan(x,y,n+1).
(77)
  • When a + b is a non-positive integer and b is positive, the first −ab terms diverge as Γ(−m) and independently the term with i = −a diverges as 1/(aa),

uδmbδ+b(x,y)=i=0mBimbmi(x)(y)ii!+i=1+mm+b1Γ(im)(imb)xim(y)ii!+Ab(x)ym+b(m+b)!+ũmbb(x,y,m+b+1).
(78)
  • When a and b are both non-positive integers, the first −ab terms diverge as Γ(−mn) with the i = −a term diverging faster, as Γ(−mn)/(aa),

uδmδn(x,y)=i=0,imm+nBimm+ni(x)(y)ii!+Cn(x)(y)mm!+ũmn(x,y,m+n+1).
(79)

Finally, examining (45), it is clear that we also need to find the Laurent series representation of the function uab defined as

uab(x,y)=Γ(a)xaΓ(b)yb.
(80)

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

uδnδ+b(x,y)=Γ(b)yb(x)nn!1δ+Ψ(1+n)+Ψ(b)log(xy).
(81)

We first consider the case of ln. 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

eλxλyGnl(x,y)xqyqdxdy,=(1)n+lni=lmin[n,lq](li)!l+ni+l2niαqCiq(α)hn,lq,nα2+i=lmin[n,lq](li)!l+ni+l2niαqCiq(α)hn,lq,nα2+i=qn(li)!l+ni+l2niAi+q(α)αqhn,lq,nα2+i=qn(li)!l+ni+l2niAi+q(α)αqhn,lq,nα2+i=lmin[n,lqq](li)!l+ni+l2nibi(q,q,α,α)+bi(q,q,α,α)+O(δ0),
(82)

where

hn,l(q,x)=i=qlxi(li)!(iq)!l+il+n(1)i+l1+1xiql+iln
(83)

and

bi(q,q,α,α)=j=liqqk=0,kjqjiqqBj+q+kijqqk(α)(l+j)!2nj1k!×ljl+n(1)j+l2nαkljln(α)k=(α)lqqiδdn,llqqi,lq,nα2,nα2+O(δ0),
(84)

where

dn,l(p,r,x,y)=i=0pj=0,jripixijyj(pji)!(j+ir)i!j!2lil+n(1)i1y+1j2liln.
(85)

In Sec. IV, we show that both hn,l(q, x) and dn,l(p, r, x, y) vanish for every integer value of 0 ≤ ql and 0 ≤ pl, respectively, proving that (41) converges in these cases.

For the case of n > l, we only need to consider Gnl(sg)(x,y) as Gnl(nsg)(x,y) never contains any singularities. In the case of l = 0, we get

eλxλyGn0(sg)(x,y)xq+δyq+δdxdy=i1=0n1nni1(1)i1(i1)!2ni11uq+δq+i1+1+δα,α2n+uq+δq+i1+1+δα,α2nuq+δ1+q+i1+δ(α,α)uq+δ1+q+i1+δ(α,α),
(86)

and by the use of Eq. (76),

eλxλyGn0(sg)(x,y)xq+δyq+δdxdy=i1=0n1nni1(1)i1(i1)!2ni11Aq+i1+1(α)(2/nα)q(q)!+Aq+i1+1(α)(2/nα)q(q)!(q+i1)!α1+q+i1(2/nα)qδ(q)!(2/nα)qδ(q)!(q+i1)!α1+q+i1+O[δ],
(87)

in which the divergent part vanishes by the definition of A.

For the case of l > 0, we get

eλxλyGnl(sg)(x,y)xq+δyq+δdxdy=i1=ln1i2=1l1+ln+lni11(1)i1+l(i21+l)!(i1l)!2ni1i2×ni2nl1uqi2+1+δq+i1+1+δα,α2n+uqi2+1+δq+i1+1+δα,α2nn+ln1+i2uqi2+1+δ1+q+i1+δ(α,α)+uqi2+1+δ1+q+i1+δ(α,α),
(88)

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

eλxλyGnl(sg)(x,y)xq+δyq+δdxdy=i1=ln1i2=1+q1+ln+lni11(1)i1+l(i21+l)!(i1l)!2ni1i2×1δni2nl1(q+i1)!αq+i1+1(2/nα)i2q1(i2q1)!n+ln1+i2(q+i1)!α1+q+i1(α)i2q1(i2q1)!+i1=ln1i2=1+q1+ln+lni11(1)i1+l(i21+l)!(i1l)!2ni1i2×1δni2nl1(q+i1)!αq+i1+1(2/nα)i2q1(i2q1)!n+ln1+i2(q+i1)!α1+q+i1(α)i2q1(i2q1)!+O(δ0),
(89)

we can perform the sum over i1 and isolate the sum over i2 to get

eλxλyGnl(sg)(x,y)xq+δyq+δdxdy=2nln+lnl11δ(q+l)!α1+l+qF211+ln,1+l+q,2l+2,2nαcn,l(q,2nα)αq+1
(90)
+(q+l)!α1+l+qF211+ln,1+l+q,2l+2,2nαcn,l(q,2nα)αq+1+O(δ0),

where

cn,l(q,x)=i=1+q1+lΓ(i+l)(iq1)!xinil+1i1+1xiq1n+ll+1i.
(91)

In Sec. IV, we show that cn,l(q, x) vanishes for every integer value of 0 ≤ ql, provided that n > l, proving that (45), indeed, converges.

We now define the functions fab to be the constant term, i.e., ∼δ0 of the Laurent series of uab, given in Eqs. (76)–(79), thus arriving at our main result for integer values of q and q′.

For the case of non-positive integer, we can write the Laurent series of Fn as

Fδn(x,y)=i=0n1xiin+δyin+δi!+xnδyδn!+i=n+1xiin+δyin+δi!=i=0n1xiinyini!+xnn!1δ+lny+yx1+n(1+n)!F22(1,1,2,2a,xy)+O(δ),
(92)

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

0eλyyq+δGnl(sg)(x,y)dy=exni1=ln1i2=1+q1+ln+l1ni11(1)i1+l(i21+l)!(i1l)!2ni1i2×ni2nl1xi1Fq+1i2+δ(α2/n,x)n+ln1+i2xi11α1+qi2+δΓ[qi2+1+δ]+O(δ0)=exni1=ln1i2=1+q1+ln+l1ni11(1)i1+l(i21+l)!(i1l)!2ni1i2×ni2nl1xi1(α2/n)i21q(i21q)!1δn+ln1+i2xi1α1+qi2(1)qi2+1(i2q1)!1δ+O(δ0)=n+lnl1n2xMn,l+1/22xnα1qcn,lq,nα21δ+O(δ0),
(93)

where M is the first Whittaker function.

On the other hand, the case of ln comes out as

0eλyyq+δGnl(x,y)dy=(1)n+lni=lnj=lq(li)!(l+j)!n+li+l2ni+jexn×ljl+n(1)j+lxi1Fj+q+δ(2/nα,x)ljln(xi1Fj+q+δ(α,x))+O(δ0)=(1)n+lni=lnj=lq(li)!(l+j)!n+li+l2ni+jexnxi1×ljl+n(1)j+l(2/nα)jq(jq)!1δljln(α)jq(jq)!1δ+O(δ0)=(1)n+l2(2l)!n2xl+1F11nl,2l,2nx(α)qhn,lq,nα21δ+O(δ0).
(94)

Since cn,l(q, x) and hn,l(q, x) both vanish for every integer value of 0 ≤ ql, we can redefine Fq(x,y) to the constant term of its Laurent series completing the derivation of J.

Here, we show that cn,l vanishes for all integers ql in the case of n > l, while hn,l vanishes for all ql in the case ln.

The function cn,l is given by a difference of two finite series. Let us call them cn,l(1) and cn,l(2). Using binomial theorem, we get

cn,l(1)(q,x)=i=q+11+lΓ(i+l)(iq1)!xinil+1i1+1xiq1=i=q+11+lk=0iq1Γ(i+l)(iq1)!iq1kxiknil+1i.
(95)

Relabeling indices and using properties of binomials, this transforms into

cn,l(1)(q,x)=i=q+1l+1xiΓ(i+l)(iq1)!k=01+lik+l+i1l+i1nkinl1=i=q+1l+1xiΓ(i+l)(iq1)!n+lli+1=cn,l(2)(q,x),
(96)

where in the last step, we have used16 

i=0bca+kabkc=a+b+1bc.
(97)

Proceeding along those same lines for hn,l, we can similarly use the binomial theorem to write

hn,l(1)(q,x)=i=ql(1)l+ixi(li)!(iq)!l+il+n1+1xiq=i=qlk=0iq(1)l+ixik(li)!(iq)!l+il+niqk
(98)

and redefine summation indices to arrive to

hn,l(1)(q,x)=i=qlk=0li(1)l+i+kxi(lki)!(kq+i)!l+k+il+nkq+ik.
(99)

Finally, using properties of binomials, we get

hn,l(1)(q,x)=i=qlk=0li(1)i+l+kxi(li)!(iq)!l+k+il+nlik=i=qlxi(li)!(iq)!l+iln=hn,l(2)(q,x),
(100)

where in the last step, we have used the following identity valid whenever b < a,16 

k=0ba+kb+cbk(1)b+k=ac.
(101)

Here, we show that dn,l(p, r, x, y) vanishes for all integer values of 0 ≤ p ≤ 2l and 0 < r < l, assuming that ln.

First, dn,l is given as a difference of two series, that we denote, as dn,l(1) and dn,l(2). Thus, we have

dn,l(1)(p,r,x,y)=i=0pj=0,jripixijyj(pij)!(j+ir)i!j!2lil+n(1)i1y+1j.
(102)

We can use the binomial theorem to transform it to

dn,l(1)(p,r,x,y)=i=0pj=0,jripik=0jxijyjr(pij)!(j+ir)i!j!2lil+n(1)ijk.
(103)

Redefinition of the summation variables, by kj, ik, and ji + jk and changing the order of sums yield

dn,l(1)(p,r,x,y)=i=0pj=0,jripik=0ixijyj(1)k(pij)!(j+ir)(i+jk)!k!2lkl+ni+jkj.
(104)

Using properties of binomials, we can write it as

dn,l(1)(p,r,x,y)=i=0pj=0,jripik=0ixijyj(pij)!(j+ir)i!j!2lkl+n(1)kik.
(105)

Finally, we use

k=0caibci(1)i=acab
(106)

to arrive at

dn,l(1)(p,r,x,y)=i=0pj=0,jripixijyj(pij)!(j+ir)i!j!2liln=dn,l(2)(p,r,x,y),
(107)

completing the proof.

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′ and exponentials. The main complications came from the fact that when 0 ≤ q, q′ < l + 1, the integrals from the individual terms of the RCGF possess singularities that are explicitly canceled 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′ and then find the Laurent series of (33) and (34) around q = m + δ, q′ = m′ + δ, where m and m′ are non-negative integers. After this, we demonstrated that for q, q′ ≥ 0, the divergent parts, that is, the terms proportional to δ−1 and δ−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.

See the Mathematica notebook (supplementary material) where we have programmed the main results of the article and provided the comparisons of analytical vs numerical evaluations of the generating integral and integral moments.

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 Ph.D. thesis work of Kamil Dzikowski to be submitted to the Heidelberg University, Germany.

1.
O. D.
Skoromnik
,
I. D.
Feranchuk
,
A. U.
Leonau
, and
C. H.
Keitel
,
J. Phys. B: At., Mol. Opt. Phys.
50
,
245007
(
2017
).
2.
H. F.
Hameka
,
J. Chem. Phys.
48
,
4810
(
1968
);
H. F.
Hameka
,
J. Chem. Phys.
47
,
2728
(
1967
).
3.
L. C.
Hostler
,
Phys. Rev.
178
,
126
(
1969
).
4.
L. C.
Hostler
and
W. W.
Repko
,
Phys. Rev. A
23
,
2425
(
1981
).
5.
G. S.
Adkins
,
J. Phys.: Conf. Ser.
1138
,
012005
(
2018
).
6.
7.
N.
Sibalic
,
Rydberg Physics
(
IOP Publishing
,
2018
).
8.
B. R.
Johnson
and
J. O.
Hirschfelder
,
J. Math. Phys.
20
,
2484
(
1979
).
9.
R. N.
Hill
and
B. D.
Huxtable
,
J. Math. Phys.
23
,
2365
(
1982
).
10.
L. D.
Landau
and
E. M.
Lifshitz
,
Quantum Mechanics: Non-Relativistic Theory
(
Butterworth Heinemann
,
1977
).
11.
I. S.
Gradštejn
,
J. M.
Ryžik
,
A.
Jeffrey
, and
D.
Zwillinger
,
Table of Integrals, Series and Products
, 7th ed. (
Elsevier Academic Press
,
Amsterdam
,
2009
), oCLC: 846182468.
12.
Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables
, Dover Books on Mathematics, 9th ed., edited by
M.
Abramowitz
and
I. A.
Stegun
(
Dover Publication
,
New York, NY
,
2013
), oCLC: 935935300.
13.

We work in the rescaled variables, i.e., x = Zr, x′ = Zr′ and λ = β/Z, λ′ = β′/Z. Therefore, the factor from the Jacobian Zqq′−2 should be included when relating Knl(β, β′) and Knl(λ, λ′).

14.
Intel Corporation
, Intel Core i7-2600K Processor, https://ark.intel.com/content/www/us/en/ark/products/52214/intel-core-i7-2600k-processor-8m-cache-up-to-3-80-ghz.html (2020); accessed 13 January 2020.
15.
Wolfram Research, Inc.
, Mathematica, Version 11.3,
Champaign, IL
,
2018
.
16.
A. P.
Prudnikov
,
J. A.
Bryčkov
, and
O. I.
Maričev
,
Integrals and Series
(
Gordon and Breach Science Publication
,
New York, NY
,
1986
).

Supplementary Material