Convergence of the Laplace and the alternative multipole expansion approximation series for the Coulomb potential

Multipole expansion is a powerful technique used in many-body physics to solve dynamical problems involving correlated interactions between constituent particles. The Laplace multipole expansion series of the Coulomb potential is well established in literature. We compare its convergence with our recently developed perturbative and analytical alternative multipole expansion series of the Coulomb potential. In our working, we confirm that the Laplace and the analytical alternative multipole expansion series are equivalent as expected. In terms of performance, the perturbative alternative multipole expansion series underapproximate the expected results to some extent while the Laplace and the analytical alternative multipole expansion series yield results which are relatively accurate but oscillatory in nature even with a higher number of angular momentum terms employed. As a practical example, we have evaluated the Slater double integrals for two-electron systems using the multipole expansion techniques and a mean field approximation. The estimated results show that only spherical interactions are dominant while the higher-order interactions are negligible. To highlight the discrepancy in the application of each of the formulations of the multipole expansion series for the electron-electron interaction potential, an estimation of the non-relativistic groundstate energies of some helium-like systems, evaluated using the spherical approximation of the multipole potential, is provided. Our findings are likely to be useful in the treatment of the Coulomb potential in electronic structure calculations as well as in celestial mechanics.

The Laplace multipole expansion series is established in the works of Laplace and Legendre in their search for solutions to the problem of attractions.The historical developments that led to the derivation of the expansion series and the introduction of the Legendre polynomials, for the first time, as the coefficients used in the Laplace expansion are captured in Laden's thesis 1 .The Laplace multipole expansion has become conventional knowledge in physics textbooks 2 and it is quite useful in solving the many-body physics problems in celestial mechanics, quantum physics and chemistry, nuclear physics, and condensed matter physics.
Naturally, the multipole expansion becomes convenient to use in solving physical problems in 3D if expressed in the spherical polar coordinates.This decomposes the problem as a product of both radial and angular parts.The radial part can be treated as a 1D case while the well defined angular momentum algebra 3 can be used to simplify the angular parts.Several studies have employed multipole expansion techniques in the recent past in solving physical problems of interest [4][5][6][7][8][9] .
In our alternative multipole expansion of the Coulomb potential 10,11 , we stated that the Laplace multipole expansion series of the Coulomb repulsion term is incomplete, and therefore inaccurate.Vaman clarified that both the Laplace and the alternative multipole expansion are indeed equivalent 12 .Since the Laplace expansion series is a single-index summation while the alternative method is a double-index summation series, it becomes necessary to test the conditions for convergence of the two methods.We also compare the accuracy of the Laplace expansion method, relative to our perturbative and analytical alternative multipole expansion methods, in estimating the expected results.We have seen in literature that such a comparison, not exactly similar to the current study, is reported in ref. 13,14 .Comparison of different methods allows characterization of relative accuracy and capabilities, which is quite instrumental in guiding application 13 .This is particularly important given the fact that the alternative multipole expansion has already been successfully employed in determining the electronic structure for neutral atoms 15,16 .

Theory
The Coulomb repulsion potential term can be expressed as which reduces to the Laplace multipole expansion series, where t = r < /r > , r > = max{r i , r j } , r < = min{r i , r j } , x = cos θ , with θ being the relative angle between the posi- tion vectors r i and r j , l are non-negative integers, and P l (x) are the lth order Laplace coefficients of tl , also known as the Legendre polynomials.It is important to note that the form given by Eq. ( 1) is considered as the generating function for the Legendre polynomials 2,17 .
In the alternative approach 10,11 , the multipole expansion of the Coulomb potential can also be expressed in the basis of Legendre polynomials, where the coefficients are a function of the spherical Bessel-like functions, jl ( t) , which can be expanded in the perturbative polynomial form as 10,11 or analytically as a differential equation 11 with and defined in terms of t in this case.The equivalence of Eqs. ( 2) and (3) shows that is an identity.
From the identity relation in Eq. ( 10), we can further infer that: (1) Using the relations given by Eq. ( 9), we have analytically tested and confirmed the inferences given by Eqs. ( 11)- (12) for the first two orders of the spherical Bessel-like functions herebelow.The zeroth-order spherical Bessel-like function simplifies to: Likewise, the first-order spherical Bessel-like function simplifies to: The use of the recurrence relations given by Eqs. ( 13) and ( 14) can be useful in eliminating singularities associated with the analytical expression of the spherical Bessel-like functions, jl ( t) 11 , as t → 0.
As a practical example, we use the Laplace and the alternative multipole expansion series, within a meanfield approximation, to estimate the Slater double integrals F l18 for helium-like systems, where the higher-order terms involve the exchange of angular momentum quantum number between the s and the l th orbital.The optimization is based on the root mean value of t, per solid angle, obtained by determining the root mean square value of t scaled by 2π as given by Eq. (32) of ref. 16 .This, consequently, yields The unscreened hydrogenic radial orbitals can be employed as the trial wavefunctions in evaluating Eq. (17).
In solving the Schrödinger equation for the two-electron atomic systems, with single-electron wavefunction φ(r i ) and the orbital energy eigenvalues ǫ α i , the single-electron Hamiltonian operator (in atomic units) can be expressed as where γ l i = 1 2 is the orbital angular momentum-dependent partitioning fraction 15,16 for the systems in their groundstate.The first, second, and third terms in Eq. ( 21) are the kinetic energy term, the electron-nuclear interaction, with Z as the nuclear charge, and the electron-electron interaction respectively.The second and the third terms form the potential energy function for the systems.
Within the spherical approximation, the potential energy functions for the system reduce to www.nature.com/scientificreports/using the Laplace multipole expansion series, or to using the alternative multipole expansion series.Making use of the mean field approximation given by Eq. ( 18), and the minimization of the lowest-order term of the multipole potential energy function, as derived in refs. 15,16, Eq. ( 23) can be re-written as a singleelectron potential with an effective nuclear charge The second term of Eq. ( 25) can be seen as a nuclear charge screening parameter which is dependent on the magnitude of nuclear charge and the partitioning fraction.
Likewise, Eq. ( 22) can be expressed as with expressed in terms of the spatial coordinate, r i , using the Slater integral as defined in Eq. ( 17) or probabilisti- cally as assuming that r i and r j have an equal probability of being greater.Other possible assumed scenarios can include r i ≤ r j or r j ≤ r i .
For the alternative multipole expansion method, the energy eigenvalue for an orbital with the principal quantum number n i , can be obtained directly by using the effective potential defined in Eq. (25).However, for the Laplace multipole expansion method, the Schrödinger equation must be solved numerically using the potential defined by Eqs. ( 26) and (27) or perturbatively using as the adjusted effective charge corresponding to the scenario assumed in Eq. ( 28).With the energy eigenvalues, the non-relativistic groundstate energy for two-electron atomic systems, with 1s 2 configuration, can be determined by 16 employing the anti-symmetric indistinguishable Hartree-Fock wavefunctions with the spatial and spin coordinates � q i = (� r i , � s i ) , in the solution of the Schrödinger equation for the separable Hamiltonian derived using the multipole expansion series.

Results
Our goal in this work is to test the convergence of the Laplace and the alternative multipole expansion series and also to compare the performance of both methods in estimating the exact function given by Eq. ( 1).The spherical Bessel-like functions, jl ( t) , used in the alternative multipole expansion can be evaluated perturbatively as given by Eqs. ( 5) and ( 6) or analytically as given by Eq. ( 7).Our calculations for convergence and performance are computed both perturbatively and analytically.As an example to test the convergence of the two approaches in real physical applications, we have calculated the Slater integrals and the non-relativistic groundstate energy of two-electron helium-like systems within the spherical approximation of the electron-electron interaction.In Fig. 1, we plot the convergence of the first two orders of the Laplace functions, tl , relative to the alternative multipole expansion functions, h l ( t) , as given by Eq. (10).The domain 0 ≤ t ≤ 1 has been chosen to coincide with the regime of convergence of the Laplace multipole expansion series.The convergence tests should confirm the validity of the identity relations given by the stated equation.Since h l ( t) is an infinite series function, it can be seen that only three terms (with k max = 2 ) of the summation series already yield reasonable trend of conver- gence, albeit slowly.In subsequent figures, we use h k max =2 l ( t) as our best converged perturbative results.From Eqs. ( 4)-( 6), it can be seen that the divergence between the Laplace functions and the perturbative alternative multipole expansions stems from the approximation of the factor 1 + t2 → 1 as k max → 0.
In Fig. 2a, we compare the convergence of the perturbative results with the corresponding analytical, h l ( t) , functions and the Laplace basis functions as given by Eqs.(10) and (12) for the first six orders of l.As already shown in Fig. 1, except at lower values of t , the perturbative basis functions do not agree fully with the cor- responding Laplace basis functions in all the cases considered.As expected, the analytical basis functions, on the other hand, show an excellent agreement with the corresponding Laplace basis functions.In the figure, the analytical h l ( t) functions and the Laplace basis functions are overlapping in the entire domain.In Fig. 2b, we show the relative deviation between the analytical and the Laplace basis functions.The relative deviations are calculated as the absolute difference between the analytical h l ( t) and the Laplace f l ( t) = tl functions divided by the Laplace functions.The observed relative deviations can be attributed to numerical noise as well as to the divergences due to singularities in the analytical function as t → 0.
Because of the slow convergence of the perturbative functions, it became of importance to test the performance of the expansions in approximating the value of the analytic function given by Eq. (1) for various values of t across the angular spectrum.The performance results are summarized in Fig. 3 for all values of x = cos θ .For lower values of t , fewer angular momentum values are necessary for convergence.For t = 0.75 , reasonable convergence is obtained with l max = 10 .The perturbative expansion on the other hand converges faster with fewer values of l max and k max , although the expected results are underapproximated to some extent using this approximation.In particular, complete convergence for the perturbative expansion is obtained using l max = 5 and k max = 2 only.As t → 1 , a higher number of angular momenta are necessary for convergence if the Laplace or the analytical multipole functions are used.In Fig. 4, we show that for t = 1 convergence of the expected func- tion is not yet achieved even with l max = 30 for the Laplace and the analytical multipole expansion.It can also be observed in Fig. 4 that as the angular momenta increases, the period and the amplitude of oscillation of the Laplace and the analytical multipole expansion results reduces.The perturbative expansion, on the other hand, is converged with less angular momenta and shows remarkable stability in the approximation of the expected function.The perturbative results offer the possibility to isolate features that are dependent on the lower order terms of the multipole expansion of the Coulomb potential.
The equivalence between the Laplace and the analytical alternative multipole expansion methods provides a wider choice of techniques to use when dealing with the Coulomb repulsion term.The Laplace basis functions , summed up to the maximum value ( k max ), plotted using left and right hand side of Eq. ( 10) respectively.The black solid line corresponds to the Laplace basis functions, tl .appear simpler, in comparison with the analytical functions, but the underlying difficulty lies in the uncertainty of identifying the r > and r < variables.In the analytical alternative multipole expansion, on the other hand, it is not necessary to identify the r > and r < variables in the mathematical manipulation of the electron-electron interaction.Additionally, as already shown in references 15,16 , the correlated term becomes separable in the alternative multipole expansion within some meanfield approximation, yielding a definite nuclear charge screening parameter, making it quite favourable to use for computations.
As a practical example, we have compared the Slater integrals 18 for the 1s − nl interacting states estimated using the Laplace and the lowest-order perturbative alternative multipole expansion series for two electron systems as expressed in Eq. (17).The results are presented in Table 1.The calculations have been done using the unscreened hydrogenic radial wavefunctions and the mean-field approximation in Eq. (19).As expected the analytical alternative multipole expansion yields results equivalent to the Laplace multipole expansion results.The lowest-order perturbative alternative multipole expansion results, on the other hand, are slightly less by a factor.From the results presented in the table, it is evident that the higher order multipole interactions are negligible and only become important when the lower -order interactions vanish.
with the value k max = 2 , and the analytical h l ( t) = fl ( t) , plotted using left and right hand side of Eq. ( 12) respectively.The solid and the dash-dot lines represent the perturbative and the analytical h l ( t) functions, as given by Eqs. ( 5)- ( 7), while the dashed lines represent the Laplace basis functions, f l ( t) = tl , respectively.The analytical h l ( t) functions and the Laplace basis functions are overlapping in the figure.(b) The relative deviation given as the absolute difference between the analytical h l ( t) and the Laplace f l ( t) = tl functions divided by the Laplace functions.As a test of convergence in the application of both multipole series expansion methods, we have estimated the non-relativistic groundstate energies of helium-like systems, with 1 ≤ Z ≤ 10 , using the spherical approximation of the series.For ease of comparison, we have generated separable interaction potentials corresponding to both expansion series.Making the Laplace potential separable is quite a challenging task though.To achieve this goal, we have created three different scenarios representing all possible configurations in the electronic distributions.Firstly, we have assumed that the spatial coordinate r i is greater that the other coordinate r j in the entire domain.Secondly, we assumed that the reverse is true ( r j > r i ), and lastly we have assumed that either r i or r j has an equal probability of being greater.With separable potentials, it is possible to solve the Schrödinger equation analytically for the two-electron systems.To enhance clarity in our comparison, we use the Hartree-Fock expansion of the wavefunction for all the cases considered.We have also avoided the use of self consistent calculations, which usually require several iterations for convergence, by choosing to use the separable potentials.The results of our estimation are presented in Table 2 in comparison with the literature values 19 .
From the results in Table 2, it can be seen that, except for the negatively charged hydrogen ion which has varied results, the Laplace multipole expansion approximation yields tight binding groundstate energies as presented in the three scenarios.The scenario with r i ≤ r j have the highest binding energies while the scenario with r j ≤ r i have the least among the three scenarios presented.The correlation between the scenarios is also evident, that is, the corresponding results have a common difference among the positively charged ions.Using the literature values 19 as our benchmark, we can say that Lap 1 results are the best case scenario for the Laplace multipole expansion for two-electron systems since the groundstate energy of the negatively charged hydrogen ion is quite close to the reference value.The effect of the electron-electron interaction for the Laplace scenarios given by Eq. (1) , at t = 1.00 , as a function of the angular momenta ( L max = 10 and L max = 30 ) with k max = 2 .The black solid line is the expected curve.The Laplace functions are denoted by g L L while the perturbative function by g k L .The logarithmic scale has been used for clarity.
Table 1.Comparison of Slater integrals for the 1s − nl interacting states evaluated using Eqs.(17) and (19) for the Laplace and the lowest-order perturbative alternative multipole expansion of the Coulomb repulsion term.The calculations have been done using the unscreened hydrogenic radial wavefunctions.www.nature.com/scientificreports/can be seen to be a weak perturbation at higher Z values.This is because the nuclear charge screening term appears to be independent of the nuclear charge.
The alternative multipole approximation results in Table 2, on the other hand, presents the least binding energies for the positive charged ions.This is because of the optimization used in generating the separable potential, which associates the nuclear charge screening parameter with the nuclear charge itself and in the process, maximising the electron-electron interaction.However, it can be observed that, except for the negatively charged hydrogen ion, the binding energies for the alternative multipole expansion approximation are still higher than the literature values 19 .The disparity with literature values widens with increase in charge polarisation, Z p = Z − Z e , with Z as the nuclear charge and Z e as the total electronic charge.This could be a pointer to a possible inadequacy of the Hamiltonian used in accounting for some of the interactions.Future investigations can probe the disparity further, for example, by including the higher-order multipole interactions as well as the relativistic effects like spin-orbit coupling and spin-spin interactions.

Conclusion
The convergence as well as the performance of the Laplace multipole expansion of the Coulomb potential, in comparison with our recently developed alternative multipole expansion series, is investigated in this study.We have confirmed that the Laplace and the analytical alternative multipole expansion series are indeed equivalent and offer a higher degree of accuracy if a larger l max is used in the approximation.The perturbative alternative multipole expansion, on the other hand, converges with a much lower value of l max and k max and is stable against oscillations in results as t → 1 but the converged results underapproximate the expected results to some extent at all angles.The stability of the perturbative results may be useful in isolating physically meaningful features even with less angular momenta in converged results.
As a practical example, we have shown that, despite the equivalence of the two formulations of the multipole expansion series, the results generated by each form can be very different depending on the optimization procedures used.In the present case, in the solution of the non-relativistic groundstate energy for helium-like systems with 1 ≤ Z ≤ 10 , the alternative multipole expansion leads to an effective nuclear charge with a charge-dependent nuclear screening parameter.The Laplace multipole expansion, on the other hand, leads to a nuclear screening parameter which is charge independent.The disparity in the results obtained in the two approaches is large.Comparatively, the alternative multipole expansion performs better relative to the literature values, with the disparity in results varying as a function of the charge polarisation between the nucleus and electrons.
Table 2. Comparison of the non-relativistic groundstate (1s 2 ) energies of helium-like systems with 1 ≤ Z ≤ 10 estimated using Eqs.( 22) and (23) for the Laplace and the alternative forms of multipole expansion of the Coulomb repulsion term respectively.Lap. 1 and Lap. 2 results have been evaluated using Eq. ( 22) assuming that r > = r i and r > = r j respectively, while Lap. 3 results have been evaluated using Eq. ( 28), considering both possibilities with equal distribution.The calculations have been benchmarked with ref. 19 data extracted from the provided ionization energies and electron affinities, that is, E 1s−1s = − 1 2 Z 2 − I p , where I p is the ionization energy (or electron affinity). https://doi.org/10.1038/s41598-023-42724-8

Figure 2 .
Figure 2. (Color online) (a) Comparison of the six functions of f l( t) = t l , h l ( t) = f k max

Figure 3 .
Figure 3. (Color online) Convergence of the Laplace and the perturbative alternative multipole expansion series in comparison to the expected function g(x, t) = (1 − 2x t + t2 ) − 1 2 given by Eq. (1) , at t = 0.75 , as a function of: (a) the angular momenta L with k max = 2 and, (b) k max with L max = 10 .The black solid line is the expected curve.The blue and red solid lines in (a) are overlapping.The Laplace functions are denoted by dashed lines in (a) and g L in (b).

Figure 4 .
Figure 4. (Color online) Convergence of the Laplace and the perturbative alternative multipole expansion series in comparison to the expected function g(x, t) = (1 − 2x t + t2 ) − 1 2given by Eq. (1) , at t = 1.00 , as a function of the angular momenta ( L max = 10 and L max = 30 ) with k max = 2 .The black solid line is the expected curve.The Laplace functions are denoted by g L L while the perturbative function by g k L .The logarithmic scale has been used for clarity.