A new expression for the Coulomb Potential corresponding to the product of two Exponential Functions based on the properties of the Integral Representations of the Bessel Functions

ABSTRACT

The calculation of the Coulomb Potential corresponding to the product of two Exponential Type Functions, inherently has numerical challenges that must be resolved. In order to address these problems, in this paper it is presented a new partition of the Coulomb Potential. The proposed partition involves two terms. One of the terms is a one-dimensional integral, which allows geometrical and statistical interpretations. The other term is proportional to a Modified Bessel Function and it is obtained from a two-step procedure. As a first step, a Non-Rational Function is used for approximating one of the two integrals involved. Then, the remaining improper integral can be identified with an integral representation of an appropriate Modified Bessel Function. The existence of such a Non-Rational Approximant is proved and its numerical performance is shown through some examples

Keywords:
Bessel functions; non-rational functions; integral representation; improper integrals; oscillating integrand; exponential type functions; Coulomb potential

1 INTRODUCTION

In this paper, it is presented an alternative treatment of the corresponding presented in References ⁵ and ⁴ for the Coulomb Potential (CP). Explicitly, the CP generated by the product of two Exponential Type Function (ETF) ⁵, denoted by V, is given by

$$\begin{gathered} V\left(\overrightarrow{R}\right)\begin{array}{ll}=& \int d\overrightarrow{r}\frac{\exp (-\alpha |\overrightarrow{r}-\overrightarrow{A}\left|\right)\exp (-\beta |\overrightarrow{r}-\overrightarrow{B}\left|\right)}{|\overrightarrow{r}-\overrightarrow{R}|}\end{array} \\ =C\frac{2^{5/2}{\pi }^{1/2}}{D^5}{\int }_0^{1}duu(1-u){\int }_0^{\infty }dt\cos \left(Dt\right)\frac{\left[1-\exp (-p\gamma )\left\{1+\frac{5p\gamma }{8}+\frac{(p\gamma {)}^2}{8}\right\}\right]}{p{\gamma }^6\left[u\right(1-u){]}^3}, \end{gathered}$$ (1.1)

where α and β are positive constants; $\overrightarrow{A}$, $\overrightarrow{B}$ and $\overrightarrow{R}$ are fixed three dimensional vectors and

$\gamma =\frac{[t^2+u{\alpha }^2+(1-u){\beta }^2{]}^{1/2}}{\left[u\right(1-u){]}^{1/2}}$ (1.2)

with $D=|\overrightarrow{A}-\overrightarrow{B}|,p=p\left(u\right)=|u\overrightarrow{A}+(1-u)\overrightarrow{B}-\overrightarrow{R}|$ and $C=(\frac{32}{\pi }{)}^{1/2}\alpha \beta D^5$.

This work proposes a new partition of equation (1.1), given by

$$\begin{gathered} V\left(\overrightarrow{R}\right)=C\frac{2^{5/2}{\pi }^{1/2}}{D^5}{\int }_0^{1}du\frac{u(1-u)}{p}{\int }_0^{\infty }dt\frac{\cos \left(Dt\right)}{[t^2+u{\alpha }^2+(1-u){\beta }^2{]}^3}- \\ -C\frac{2^{5/2}{\pi }^{1/2}}{D^5}{\int }_0^{\infty }dt\cos \left(Dt\right)g\left(t\right), \end{gathered}$$ (1.3)

where

$g\left(t\right)={\int }_0^{1}du\frac{\left[u\right(1-u){]}^{-2}}{p{\gamma }^6}\left(\exp (-p\gamma )\left(\frac{(p\gamma {)}^2}{8}+\frac{5p\gamma }{8}+1\right)\right).$ (1.4)

Analysing the equation (1.3), it could be seen that the first term can be transformed into a one-dimensional integral by using the following integral representation ⁶ for the K _ν (λz) Modified Bessel Function (KMBF),

$K_{\nu }\left(\Lambda z\right)=\frac{\Gamma (\nu +1/2)}{{\pi }^{1/2}{\Lambda }^{\nu }}{\int }_0^{\infty }dt\frac{cos\left(\Lambda t\right)}{(t^2+z^2{)}^{\nu +1/2}}.$ (1.5)

In other words, denoting the first term in equation (1.3) by T ₁, it is obtained

$$\begin{gathered} T_1=C\frac{2^{5/2}{\pi }^{1/2}}{D^5}{\int }_0^{1}du\frac{u(1-u)}{p}{\int }_0^{\infty }dt\frac{\cos \left(Dt\right)}{[t^2+u{\alpha }^2+(1-u){\beta }^2{]}^3} \\ =C\frac{\pi }{D^5\Gamma \left(3\right)}{\int }_0^{1}du\frac{u(1-u)}{p\left(u\right)}\frac{{\pi }^{1/2}{D}^{5/2}{K}_{5/2}\left(D(u{\alpha }^2+(1-u){\beta }^2{)}^{1/2}\right)}{(u{\alpha }^2+(1-u){\beta }^2{)}^{5/2}} \end{gathered}$$ (1.6)

Now, considering the second term in the equation (1.3) and denoting it by T ₂, it can also be transformed into a simpler expression, approximating the function g(t) by an appropriate Non-Rational Approximant (NRA), which allows to compute the remaining improper integration using again the equation (1.5). Then, the final result that is achieved is

$T_2\simeq C\frac{2^{5/2}{\pi }^{1/2}}{D^5}g\left(0\right){\left(\frac{n}{a}\right)}^n\frac{\Gamma \left(\frac{1}{2}\right)}{\Gamma \left(n\right)}{\left(\frac{D{a}^{1/2}}{2n^{1/2}}\right)}^{n-\frac{1}{2}}{K}_{n-\frac{1}{2}}\left(D{\left(\frac{n}{a}\right)}^{1/2}\right)$ (1.7)

Therefore, the current approach is expressed in the following equation

$V\simeq T_1-T_2$ (1.8)

It is worth mentioning that in spite of using the CP ² in equation (1.1) as an electrostatic object, it is important because it can provide approximated models in many other situations. In our multi-disciplinary research group, we are interested in introducing Exponential Type Functions to perform Molecular Calculations. In this context, the goal is to find alternative evaluations of the equation (1.1), i.e. expressions that can avoid improper integrals with oscillating integrands, since thousands of them must be evaluated routinely.

In the next sections we present the details of the method. Section 2 gives the explicit construction of the NRA, i.e. the approximation of the function g(t) defined in equation (1.4) and the parameters involved in the calculation are described. Section 3 shows the numerical performance of the new expression of V, by the experimentation with some examples. Finally, Section 4 expresses the concluding remarks.

2 METHOD

The method proposed is motivated by the behaviour of the non-oscillating part in the integrand of equation (1.3), i.e. g(t) is a decreasing non-negative function with slope zero at the origin, going to zero at the infinite. Consequently, these facts together with the asymptotic behaviour of the approximant, which is given in Appendix B (see ¹, ³ for analogies), make it possible to approximate the function g(t) by the following NRA:

$App(n,t)=\frac{g\left(0\right)}{{\left(1+\frac{a}{n}{t}^2\right)}^n},$ (2.1)

where $a=-\frac{g\text{'}\text{'}\left(0\right)}{2g\left(0\right)}$. This particular choice of the parameter a makes App equal to g(t) up to the second derivative at t = 0. The parameter n is determined following the procedure described in Subsection 2.1. The existence of such n is guaranteed from the properties of the approximant App, as it is proved in the Appendix A. Then, the second term in (1.3), T ₂, can be expressed by

$$\begin{gathered} T_2=C\frac{2^{5/2}{\pi }^{1/2}}{D^5}{\int }_0^{\infty }dt\cos \left(Dt\right)g\left(t\right) \\ \simeq C\frac{2^{5/2}{\pi }^{1/2}}{D^5}{\int }_0^{\infty }dt\cos \left(Dt\right)App(n,t) \\ =C\frac{2^{5/2}{\pi }^{1/2}}{D^5}g\left(0\right){\left(\frac{n}{a}\right)}^n\frac{\Gamma \left(\frac{1}{2}\right)}{\Gamma \left(n\right)}{\left(\frac{D{a}^{1/2}}{2n^{1/2}}\right)}^{n-\frac{1}{2}}{K}_{n-\frac{1}{2}}\left(D{\left(\frac{n}{a}\right)}^{1/2}\right) \end{gathered}$$ (2.2)

2.1 Calculation of the parameter n

Two different strategies are used to estimate the parameter n. At first, the calculation of n can be carried out by means of the least squares approximation. Since this technique is computationally expensive, the computation must be done from the interpolation at a point, setting a criterion for the point selection that takes into account the behaviour of the function to be approximated. This last approach naturally leads to calculate the estimate of n using the fixed point method, as it is detailed below. To determine an estimate of n it is used the fixed point method for a previously selected value of t. For this purpose, some transformations are performed to make the problem more manageable. The problem of approximating g by App from expression (2.1) can be reformulated considering the functions $\frac{g\left(0\right)}{g\left(t\right)}$ and ${\left(1+\frac{a}{x}{t}^2\right)}^x$ or, equivalently, $\ln \left(\frac{g\left(0\right)}{g\left(t\right)}\right)$ and $x\ln \left(1+\frac{a}{x}{t}^2\right)$. Thus, the fixed point iteration function is given by

$h\left(x\right)=\frac{\ln \left(\frac{g\left(0\right)}{g\left(t\right)}\right)}{\ln \left(1+\frac{a}{x}{t}^2\right)}$ (2.3)

Taking into account the behaviour of the function g(t), it is decided to choose t as $t=2\sqrt{a^{-1}}$, which is the value of the abscissa where $\frac{g\left(0\right)}{1+a{t}^2}=\frac{1}{2}g\left(0\right)$. In order to reduce the computational effort required to achieve the estimating calculus from the least squares approximation technique, the parameter n is estimated by n _pf iterating the function h(x) six times. It is important to note that for each test example it is verified that the fixed point iteration converges on the interval [1, 10]. This convergence interval is dependent of the studied examples.

3 RESULTS

The performance of the proposed method is analyzed within the examples from Table 1. The corresponding results are shown in Tables 2 and 3. Table 2 presents the reference values that correspond to the direct numerical evaluation of equation (1.1). We include the approximated value, the relative errors, and n estimation by using the least squares method and the fixed point techniques. Table 3 presents the second term reference’s values, including the approximated value and the relative errors by using both techniques.

Table 1:
Definition of the examples: $\overrightarrow{A}$, $\overrightarrow{B}$, $\overrightarrow{R}$, α and β .

Ex	A x	A y	A z	B x	B y	B z	R x	R y	R z	α	β
1	0	0	0	1	0	0	0.5	2	0	8	1
2	0	0	0	1	0	0	0.5	0.5	0	8	1
3	0	0	0	1	0	0	0.1	0.1	0	8	1
4	0	0	0	2	0	0	0.5	1	0	2	1
5	-2.286613	-1.767509	0	-2.286613	1.767509	0	2.286613	1.767509	0	1	3
6	-2.286613	-1.767509	0	-2.286613	1.767509	0	2.286613	1.767509	0	1	8

Table 2:
RV: Reference value, corresponds to the direct numerical evaluation of equation (1.1). AV_ms: Approximated value by using the least squares method. AV_fp: Approximated value by using the fixed point technique. RE_ms: Relative error by least squares. RE_fp: Relative error by fixed point. n_ms: n estimation by using the least squares method. n _fp: n estimation by fixed point.

Ex	RV	AV ms	AV fp	RE ms	RE fp	n ms	n fp
1	8.515504843e-3	8.515504856e-3	8.515504857e-3	2e-9	2e-09	3.500	3.607
2	2.565282072e-2	2.566088705e-2	2.56639191e-2	3e-4	4e-04	1.549	1.477
3	6.388687665e-2	6.524496021e-2	6.559809098e-2	2e-2	3e-02	1.606	1.555
4	2.959858143e-1	3.231339878e-1	3.229071644e-1	7e-4	5e-04	4.148	4.224
5	6.405626287e-3	6.40562629e-3	6.405626288e-3	5e-10	6e-11	9.831	9.997
6	3.171312569e-4	3.171312569e-4	3.171312569e-4	<1e-13	<1e-13	7.207	7.461

Table 3:
T₂RV: Second term reference’s value. T ₂ AV _ms: Approximated value’s second term by using the least squares method. T ₂ AV _fp: Approximated value’s second term by using the fixed point technique. RE _ms: Relative error by least squares. RE _fp: Relative error by fixed point.

Ex	T 2 RV	T 2 AV ms	T 2 AV fp	RE ms	RE fp
1	1.537019371e-9	1.524412743e-9	1.5234828e-9	8e-3	9e-3
2	4.255363114e-4	4.1746998e-4	4.144379321e-4	2e-2	3e-2
3	7.95831452e-2	7.822506164e-2	7.787193087e-2	2e-2	2e-2
4	5.408322576e-2	5.387010451e-2	5.394448295e-2	4e-3	3e-3
5	3.177659062e-9	3.174721566e-9	3.177279844e-9	9e-4	1e-4
6	0	0	0	0	0

4 CONCLUDING REMARKS

On the one hand, one of the most important contribution of this work is the definition of the first term in the partition of equation (1.3) which results a one dimensional integral with a numerical evaluation quite simple, i.e. contrasting with the equation (1.1). Also, it can be interpreted as being proportional to the weighted average of |p(u)|⁻¹. These facts, together with the geometrical meaning of p(u), i.e. it defines a triangular region with vertices given by $\overrightarrow{A}$, $\overrightarrow{B}$ and $\overrightarrow{R}$ when all of them are not in the same line, make T ₁ a friendly quantity.

On the other hand, the second term in equation (1.7) also has a very simple evaluation after the determination of the parameters g(0), g''(0) and n. Although this makes a reduction in the computing effort, the relative errors may have a broad range of magnitudes, and the order of these moves between values less than 10⁻¹³ and 10⁻². However, it can be seen that the weight of T ₂ relative to T ₁ is important in the result. This is shown by examples 1, 2 and 3, where there is a systematic deterioration of the App when p is diminishing. The fact that T ₁ >> T ₂, gives better precision, see examples 1, 5 and 6. We consider that this last behaviour is a good result achieved from the approach presented here, which can accomplish an improved overall final error for a given set of examples, i.e. to that presented in reference ⁴.

A perspective of this work is the possibility of considering other approximants, such as those kinds of Rational Approximants presented in reference ¹, i.e. perhaps with lesser deterioration than the corresponding to NRA. It is important to note that both approaches with the integration in the complex plane are related.

Acknowledgments

Grants from SECYT (UNRC) and CONICET are gratefully acknowledged.

REFERENCES

• ¹
  J.C. Cesco, J.E. Pérez, C.C. Denner, G.O. Giubergia & A.E. Rosso. Rational Approximants to Evaluate four-center Electron Repulsion Integrals for 1s Hydrogen Slater Type Functions. Applied Numerical Mathematics, 55 (2005), 179-190.
• ²
  J.D. Jackson. “Classical Electrodynamics”. John Wiley and Sons, Inc. (1999).
• ³
  C.J.A. Lanzardo, J. Pérez & J. Cesco. Aproximación de integrales impropias cuyo integrando tiene una parte oscilante. In “Matemática Aplicada Computacional e Industrial”, volume 7. ASAMACI (2019), p. 181-184.
• ⁴
  C.J.A. Lanzardo, J. Pérez , M.L. Tardivo, G. Fraschetti &J.C. Cesco . Potencial de Coulomb para funciones 1s de Slater expresado en términos de funciones de Bessel. In “Matemática Aplicada Computacional e Industrial”, volume 8. ASAMACI (2021), p. 81-84.
• ⁵
  J.E. Pérez , O.E. Taurian, A. Bouferguene & P.E. Hoggan. On a transformation for the electrostatic potential, generated by the product of two 1s Slater type orbitals, giving an efficient expression. Advances in Quantum Chemistry, 67 (2013), 65-71.
• ⁶
  I. Stegun & M. Abramowitz. “Handbook of mathematical functions and tables, Bessel function of fractional order”. National Bureau of Standards, App Math. Series, Dover (1972).

APPENDIX A

Theorem A.1.For t > 0 fixed and large enough, exists n ∈ R, n > 0 such that the function$ϵ\left(x\right)=\frac{g\left(t\right)}{g\left(0\right)}-{\left(1+\frac{a}{x}{t}^2\right)}^{-x}$has a zero in x = n.

Proof. To prove this result, just see that for t > 0 big enough, the function ε(x) changes sign at least once in the interval (0, ∞).

Given that $\frac{g\left(t\right)}{g\left(0\right)}<1,\forall t>0$, and $\underset{x\to 0^{+}}{\lim }{\left(1+\frac{a}{x}{t}^2\right)}^{-x}=1$, the function $ϵ\left(x\right)<0$ as $x⟶0^{+}$.

On the other hand, given that $\underset{x\to \infty }{\lim }{\left(1+\frac{a}{x}{t}^2\right)}^{-x}=e^{-a{t}^2}$ and that asymptotically, i.e. for all t > 0 large enough, we have that $\frac{g\left(t\right)}{g\left(0\right)}>e^{-a{t}^2}$, then $ϵ\left(x\right)>0$ as $x⟶\infty$.

Therefore, exists $n\in (0,\infty )$ such that ε(n) = 0. □

This theorem ensures that, choosing t conveniently, it is possible to construct an approximant of the form (2.1) for g(t).

It is important to mention that the fact that $\underset{x\to \infty }{\lim }{\left(1+\frac{a}{x}{t}^2\right)}^{-x}=e^{-a{t}^2}$ shows the positive function $H\left(x\right)={\left(1+\frac{a}{x}{t}^2\right)}^{-x}$ is decreasing for t > 0 fixed. In fact, $H\text{'}\left(x\right)=H\left(x\right)\left(\frac{a{t}^2}{n\left(1+\frac{a}{x}{t}^2\right)}-\ln \left(1+\frac{a}{x}{t}^2\right)\right)$ becomes less than zero as $x⟶\infty$.

APPENDIX B

Theorem B.2 (The asymptotic behaviour of g (t )). The function

$g\left(t\right)={\int }_0^{1}du\frac{\left[u\right(1-u){]}^{-2}}{p{\gamma }^6}\left(\exp (-p\gamma )\left(\frac{(p\gamma {)}^2}{8}+\frac{5p\gamma }{8}+1\right)\right),$

where $\gamma =\frac{[t^2+u{\alpha }^2+(1-u){\beta }^2{]}^{1/2}}{\left[u\right(1-u){]}^{1/2}}$ , satisfy $g\left(t\right)=O\left(\frac{1}{t^2}\right)$ .

Proof. We have that $\gamma (0,u)<\gamma (t,u)$, for all t > 0, and p(u) > 0 for all u ∈ (0, 1). Then, $\exp (-\gamma (t,u\left)\right)<\exp (-\gamma (0,u\left)\right)$.

On the other hand,

$\frac{1+\frac{5p\gamma }{8}+\frac{(p\gamma {)}^2}{8}}{p{\gamma }^6}=p\frac{\frac{1}{(p\gamma {)}^2}+\frac{5}{8p\gamma }+\frac{1}{8}}{{\gamma }^4}$

Therefore,

$\frac{{\displaystyle 1}}{{\displaystyle (p\gamma {)}^2}}+\frac{{\displaystyle 5}}{{\displaystyle 8p\gamma }}+\frac{{\displaystyle 1}}{{\displaystyle 8}}<\frac{{\displaystyle 1}}{{\displaystyle \left(p\gamma \right(0,u){)}^2}}+\frac{{\displaystyle 5}}{{\displaystyle \left(8p\gamma \right(0,u\left)\right)}}+\frac{{\displaystyle 1}}{{\displaystyle 8}}.$

Calling $F\left(u\right):=\frac{{\displaystyle 1}}{{\displaystyle \left(p\gamma \right(0,u){)}^2}}+\frac{{\displaystyle 5}}{{\displaystyle \left(8p\gamma \right(0,u\left)\right)}}+\frac{{\displaystyle 1}}{{\displaystyle 8}}$, we have

$g\left(t\right)<\frac{1}{t^2}{\int }_0^{1}du\exp (-p\gamma (0,u\left)\right)p\left(u\right)F\left(u\right)$

Then, $g\left(t\right)=O\left(\frac{{\displaystyle 1}}{{\displaystyle t^2}}\right)$. □

Publication Dates

• Publication in this collection
  27 Mar 2023
• Date of issue
  Jan-Mar 2023

History

• Received
  29 Sept 2021
• Accepted
  30 June 2022