Convergence properties of the multipole expansion of the exchange contribution to the interaction energy\dagger {}^{\dagger} Dedicated to Professor Andreas Savin on the occasion of his 65 birthday.

Convergence properties of the multipole expansion of the exchange contribution to the interaction energy thanks: Dedicated to Professor Andreas Savin on the occasion of his 65 birthday.

Piotr Gniewek and Bogumił Jeziorski
Faculty of Chemistry, University of Warsaw, Pasteura 1, 02-093 Warsaw, Poland
Corresponding author, e-mail:
July 12, 2019

The conventional surface integral formula and an alternative volume integral formula are used to compute the asymptotic exchange splitting of the interaction energy of the hydrogen atom and a proton employing the primitive function in the form of its truncated multipole expansion. Closed-form formulas are obtained for the asymptotics of and , where is the multipole expansion of truncated after the term, being the internuclear separation. It is shown that the obtained sequences of approximations converge to the exact results with the rate corresponding to the convergence radius equal to 2 and 4 when the surface and the volume integral formulas are used, respectively. When the multipole expansion of a truncated, th order polarization function is used to approximate the primitive function the convergence radius becomes equal to unity in the case of . At low order the observed convergence of is, however, geometric and switches to harmonic only at certain value of dependent on . An equation for is derived which very well reproduces the observed -dependent convergence pattern. The results shed new light on the convergence properties of the conventional SAPT expansion used in applications to many-electron diatomics.


10.1080/0026897YYxxxxxxxx \issn1362–3028 \issnp0026–8976 \jvol00 \jnum00 \jyear2015


Research Article

Figure 1: Graphical abstract: The leading term of exchange energy can be calculated from the multipole expansion of the wave function using variational principle, surface integral, or perturbation theory.

1 Introduction

Energetic effects of noncovalent molecular interactions are very small, compared to the energies of noninteracting subsystems, nevertheless they determine the properties of condensed matter [1, 2] and influence chemical reactivity [3]. The theoretical framework for the understanding of the nature and origin of noncovalent interactions is provided by Symmetry Adapted Perturbation Theory (SAPT) [4, 5], which introduces rigorous decomposition of the interaction energy into electrostatic, induction, dispersion and exchange contributions. At large interatomic distances the first three of these constituents can be expanded in powers of , leading to the so-called multipole expansion of interaction energy


with the ’s coefficients referred to as the van der Waals constants. While the expansion of Eq. (1) is divergent for any [6, 7], it is known to provide the correct asymptotic representation of for large [8, 9]. From the perspective of the large asymptotics of , the exchange effects might seem to be insignificant. However, at intermediate and short distances the importance of the exchange part of is hard to overestimate, as it is responsible for the splitting of potential energy curves and for the strong repulsion in the case of interaction of closed-shell atoms when it provides the necessary quenching of the attractive induction and dispersion components of the interaction energy.

The exchange effects are caused by the resonant tunneling of one or more electrons between the interacting systems and are considered to be non-perturbative, that is, an application of the straightforward Rayleigh-Schrödinger perturbation treatment is unable to provide even qualitative description of the exchange contribution to [10, 11]. This problem can be circumvented by the symmetry forcing procedure [12], in which the wave function is correctly symmetrized in each order of perturbation theory. Different variants of symmetry forcing lead to different SAPT expansions [12, 13], of which the most important is the Symmetrized Rayleigh-Schrödinger (SRS) perturbation approach [10, 14]. SRS is nowadays routinely applied in studies of noncovalent interactions of many-electron systems [15, 16, 17, 18].

The success of SRS approach is somewhat bewildering, as it is based on the so-called polarization approximation (PA), which is the straightforward application of the Rayleigh-Schrödinger (RS) perturbation theory to molecular interactions. PA is known to diverge when at least one of the interacting subsystems has more than two electrons [19, 20, 21, 22]. On the other hand, PA converges for H and H to the ground state wave function and energy [10, 11, 23]. For large this convergence is very slow, as the convergence radius of PA exceeds one by a quantity exponentially decaying with [24, 11, 23, 25]. Irrespectively of the ultimate convergence or divergence, PA provides the asymptotic approximation to the primitive function [19] of the interacting system, in the sense that [8, 26]


where is the wave function of the state , is the symmetry operator associated with this state, is the PA wave function computed through the th order, and equals two if at least one of the interacting subsystems is charged and three if both are neutral.

From the practical point of view it is important to know the limits of validity and convergence properties of various approaches to calculation of the exchange contribution to . In this paper we present exact results concerning the rate of convergence of the exchange energy of the hydrogen molecular ion H when computed using the multipole expansion of the primitive wave function. The exchange energy is defined in this context as


where and are the Born-Oppenheimer energies of the lowest gerade and ungerade states, respectively. The H ion is particularly suitable for studies of the convergence of exchange effects for two reasons: () the analysis is simplified by the fact that there are no electron correlation effects involved and () the asymptotic expansion of is known exactly for this system [7]. H is one of the most important models of the theory of molecular interactions and there is a wealth of analytical results available for this species [27, 28, 13, 29, 7]. H is also a prototypical double-well system [30, 31]. The asymptotic expansion of exchange energy of H can be conveniently written in the form


where , and . The higher coefficients can be found in Ref. [7].

The most widely used method of calculating employs the surface-integral formula given by [32, 33, 34, 35]


where the primitive function is assumed to be localized in the “left” half of space and is the median plane dividing the space into “left” and “right” halves. Atomic units () are used in Eq. (5) and throughout the paper.

Alternatively, the exchange energy can be calculated from the SAPT formula [25]


where is the interaction operator, is the unperturbed wave function, and is the operator inverting the electronic coordinates with respect to the midpoint of the internuclear axis. When is approximated by the polarization function and the r.h.s. of Eq. (6) is expanded in powers of , one obtains the SRS expansion of the exchange energy through the th order.

In our recent work [36] we showed that the exchange energy can be very accurately calculated from the variational principle based formula


where is the full Hamiltonian of the interacting system. Numerical results illustrating the convergence of , , and with respect to the order of multipole expansion used to approximate the primitive wave function were presented in our previous work [36]. By analyzing these numerical results we have found that the convergence rate of the leading term, , of the asymptotic expansion of Eq. (4) corresponds to the convergence radius equal to 2, 1, and 4 for the surface-integral, SAPT, and variational formula, respectively.

In this communication we shall present a rigorous mathematical derivation of the discovered convergence rates and the corresponding convergence radii (with respect to the order of multipole expansion of ) for calculated from the surface and variational expression (in the case of such analysis has already been presented in Ref. [36]). Furthermore we derive closed-form formulas for partial sums of the series approximating and show that these partial sums converge to the correct value . We shall also explain why in the calculations of using the primitive function approximated by a finite polarization expansion the convergence is initially geometric but switches to harmonic when the multipole expansion of is carried out to sufficiently high order.

2 Asymptotic expansion of the polarization wave function

The non-relativistic Born-Oppenheimer Hamiltonian for an interacting system consisting of two molecules A and B decomposes naturally as , where is the unperturbed part ( and being the Hamiltonians of isolated monomers A and B, respectively) and is the perturbation accounting for the Coulomb interactions between charged particles belonging to different monomers. In the case of H the pertinent partitioning of is


where and are the distances of the electron from the nuclei and , respectively.

Polarization approximation is obtained when one applies the usual RS perturbation theory to the above definition of and , with the the unperturbed wave function taken as the product of the ground state wave functions of isolated subsystems (in case of H the unperturbed wave function is the ground state wave function of a hydrogen atom , ). The polarization corrections to energy are given by


whereas the wave function corrections are defined recursively by equations


and the intermediate normalization condition holding for . For we have . The wave function corrections can be written in a product form, separating the zeroth-order wave function, , where and for the factors satisfy the equation [37, 38]


When the perturbation in Eq. (11) is represented by its multipole expansion


one obtains the asymptotic expansion for


where is the inclination angle of the electron in a coordinate system centered at the nucleus and is the Legendre polynomial. Using Eq. (11) one can verify that the factors are finite order polynomials in and


such that and . Eq. (15) may be viewed as an expansion of Eq. (28) of Ref. [36] in powers of the perturbation [the function in Eq. (28) of Ref. [36] is defined via ; this definition has been erroneously omitted in the line preceding Eq. (28)].

Since and the asymptotically leading term of , defined by the coefficient in Eq. (4), can be obtained by considering a simplified formula


Evaluation of Eq. (16) with the multipole expansion of the polarization function truncated after the th term,


leads to a linear combination of integrals


The asymptotic formula (18) shows that only the terms with in Eq. (15) contribute to of the expansion (4). For this reason we will further consider only the dominant, part of , denoted , given by


where has been denoted by . Inserting Eqs. (14) and (15) into Eq. (11) and comparing coefficients at terms proportional to one obtains the relation


In deriving Eq. (20) we took advantage of the fact that and the last term in Eq. (11) do not contain terms proportional to .

Since the asymptotic value of the integral of Eq. (18) is independent of and , the angular dependence of does not affect the asymptotic estimate of Eq. (16) and, therefore, has no influence on . In other words, the value of the asymptotically leading term of depends only on the values of on the line joining the nuclei and , where , similarly as it was observed for in Ref. [36]. Thus, in the analysis of the convergence properties of approximations to it is sufficient to calculate the values of the functions on the line joining the nuclei,


where is defined as the sum


Substituting in Eq. (20) we obtain a simple recurrence relation for the coefficients


with the initial, values given by and for . In setting the lower summation limit in Eq. (23) we used the fact that for . In Eq. (23) and throughout the paper we use the convention that when the lower summation limit is greater than the upper one the sum is zero. Using Eq. (23) it is straightforward to show that


which furnishes a fast algorithm for the calculation of .

The values of the multipole expansion of , Eq. (17), truncated after term and computed on the line joining the nuclei with replaced by its dominant part , Eq. (19), are given by




The multipole expansion of the primitive function truncated after , denoted by , is


where denotes the integral part of . The finite upper limit of the summation in Eq. (27) is a consequence of the fact that for . The dominant part of , denoted by , is defined in the same way as except that is replaced by . On the line joining the nuclei the function is given by Eq. (25) with superscripts omitted and defined by


We shall now derive a closed-form formula for . This can be done via a recurrence for differences of consecutive ’s with , which using Eqs. (23) and (28) takes the form


where the lower limit of the summation was set equal to 1 since for the term does not contribute to the summation in Eq. (28). The summation over is actually finite and limited by the condition for .

Extracting the term form the left inner sum and combining the remaining part with the other sum gives


Comparison with Eqs. (23) and (26) allows to rewrite Eq. (30) as


which together with the fact that and gives


Eq. (32) agrees with Eq. (36) of Ref. [36] obtained using a different derivation, not related to the polarization expansion of the wave function.

One may note that after multiplication by the coefficients and have a combinatorial interpretation: is equal to the so-called derangement number (see Ref. [39] p. 59, Ref. [40] p. 180, and Ref. [41] entry A000166) defined as the number of permutations with no fixed points (ie. permutations in which all elements change places), while is equal to the associated Stirling number , defined as the number of derengements of elements having exactly permutation cycles (Ref. [40] p. 256, Ref. [39] p. 73, and Ref. [41] entry A008306).

3 Closed-from expressions for the partial sums and the convergence rates of the multipole expansion of the exchange energy

Approximations to obtained from and , where is any approximation to the primitive function , will be denoted by and , respectively. The convergence of these approximations resulting when is substituted by the multipole expanded polarization function can be characterized with the help of the increment ratios defined as


where is the multipole expansion of truncated after the term, cf. Eq. (17). Increment ratios are defined analogously. The convergence of results obtained when is substituted by will be characterized by the increment rations and , defined as in Eq. (33), but with superscripts omitted. Similar convention will be applied to other quantities as well — symbols without the superscript will denote the limits of the corresponding expressions carrying this superscript. It should be noted that the r.h.s. of Eq. (33) remains unchanged if the functions or are replaced by or , respectively.

3.1 Surface-integral formula

It has been shown in Ref. [38] that the large asymptotics of is determined by the value of at the mid-point of the line joining the nuclei. Specifically, if then


Thus, in view of Eq. (25) the asymptotics of is determined by the constant




Using mathematical induction one can easily prove that the limit of , denoted by , is given by


where denotes the exponential sum function


Combining Eqs. (35) and (37) one obtains a simple closed-form formula for , that is for the exchange energy determined by partial sums of the multipole expansion of the primitive function:


Taking the limit we see that the leading, term in the expansion of Eq. (4) is correctly obtained in this way, i.e,


in agreement with the results of Ref. [36].

Equation (39) can be compared with the formula for obtained by Tang et al. [38] using the surface-integral expression and the polarization function . Their Eq. (7.25) can be recast in our notation as follows:


where the terms in the square brackets involving tuple summation originate from the th order polarization function . The last term (not explicitly written) involves the tuple summation and originates from . Eq. (41) is a sum of terms of increasingly high order in , each of which is an infinite sum corresponding to the expansion in powers of . On the other hand, Eq. (39) contains only single sums over terms of different order in , each of which is a closed-form sum of contributions of different orders in .

Knowing the expression for , the increment ratio , defined by the limit of Eq. (33), can be written in the form


When the ratio is equal to 1 with an error of the order of Estimating the remaining factor using Eq. (37) one finds that


which is the result discovered numerically in Ref. [36]. We see that the sequence of approximations to obtained using the surface integral formula and the multipole expansion of the wave function converges like a series with the convergence radius equal to 2. In Sec. 4.1 we show that for finite the convergence radius corresponding to the sequence remains also equal to 2, although the rate of convergence becomes somewhat slower in this case.

Concluding this Subsection it may be remarked that the authors of Ref. [38] were only concerned with the limit of Eq. (41) and did not consider the convergence rate of the obtained series expansions. In fact, this rate cannot be simply inferred from Eq. (41). On the other hand, the convergence rate corresponding to Eq. (39) can be easily obtained and is given by Eq. (43).

3.2 Variational formula

In the case of the variational volume-integral formula it is easier to analyse the convergence rate than to obtain a closed-form expression for the partial sum of the corresponding expansion. We therefore start by considering the increment ratio of Eq. (33) first.

3.2.1 Convergence rate of the variational formula

Employing Eqs. (18) and (25) one obtains




The numerator in Eq. (33) multiplied by can now be written as




A closed-form formula for the the limits of , denoted by , can be obtained using the summation by parts formula,




We shall use this formula with and . In this case and


The latter identity can be proved with mathematical induction. The evaluation of the summation on the r.h.s. of Eq. (48) reduces now to sums of the form


where and are non-negative integers. For while for one easily obtains the recurrence relation


The latter relation allows us to express in the closed form as


Therefore the limit of , denoted by , is


and consequently


in agreement with the numerical results of Ref. [36]. Thus, the series of approximations to the exchange energy obtained from the multipole expansion of the primitive function and the variational formula converges at large like a geometric series with the term ratio equal to 4. In Section 4 we shall show that, surprisingly enough, this convergence changes to harmonic when is finite, i.e., when is replaced with its approximation by a finite-order polarization function .

3.2.2 Partial sums for the variational formula

Recalling the definition of as the numerator of Eq. (33) times one can easily see that the constant calculated from the variational formula with the multipole expansion of the primitive function through the terms can be obtained by summing the increments up to . Using Eq. (54) one finds


where the summation in Eq. (56) was separated into the contribution containing terms quadratic in and the part linear in . Both these sums can be calculated using the summation by parts formula of Eq. (48). The sum is obtained from this formula by setting and The differences of ’s are then given by , where


The sums can be proved with the mathematical induction to be equal to , where


The summation formula of Eq. (48) gives now


Using the definitions of and the term in Eq. (56) can be written in the form


where now we have set and


Note that the sums containing cancel out when the sum is taken. To execute the first sum on the r.h.s. of Eq. (60) using the summation by parts formula we need an expression for the sums . With the help of the mathematical induction we find that , where


Application of the summation by parts formula to the first term in Eq (60) gives


Now the only contributions to that still have to summed are the sum of in Eq. (59) and the sum on the r.h.s. of (63). These sums can be combined to yield


Using this result we obtain the desired closed form expression


It is readily seen that the limit of this formula, resulting from the first term on the r.h.s., is equal to in agreement with the exact value . After some manipulations one can also verify that the corresponding increment ratio of Eq. (33) agrees with the estimate of Eq. (55).

4 Convergence rate in the case of truncated polarization expansion

4.1 The surface-integral formula

When is finite the increment ratio of Eq. (33) takes the form


where and


Equation (67) may be regarded a finite precursor of Eq. (42). In view of Eq. (101) the asymptotics of is given by


After expanding Eq. (67) one obtains the following asymptotic expression for




Equation (36) shows that is a sum of rapidly decaying terms, and therefore converges quickly to its limit. Specifically it can be showed that for




For this reason asymptotics of is determined by the large behavior of . Eq. (101) implies then that at large and fixed we have


The limit of (with fixed ) is therefore


On the other hand the limit of for a fixed is


Ultimately the limit of is


in agreement with Eq. (43). Note that the combined limit of does not exist, as the result depends on the order of limits.

The asymptotics of the increment ratio of is given by


Eq. (77) agrees with the results presented for in Ref. [36]. This equation shows that for finite the convergence radius corresponding to and equal 2 is the same as for . It may seem disturbing that the limit of Eq. (77) does not coincide with Eq. (43). This apparent inconsistency results from the facts that the and limits of do not commute.

4.2 The variational formula

In this subsection we shall show that for finite the convergence of of corresponds to the convergence radius of 1, i.e. becomes harmonic, in contrast to the fast geometric convergence found for .

The inverse d’Alembert ratio can be factored as


similarly to Eq. (66) for the surface-integral formula. is defined as