Classical and quantum evaluation in the low temperature limit of the Keesom integral for the interaction between permanent dipoles

Michele Battezzati; Valerio Magnasco

doi:10.1016/j.crci.2008.09.026

Résumé

In a previous paper on the Keesom effect, we presented a procedure for the evaluation of the average potential energy of interaction between two dipoles of strength μ a distance R apart in the limit of small temperature T. The resulting expansion in terms of the dimensionless T-dependent parameter a = −μ²/R³kT was then proved by us to be an asymptotic one. Here we reformulate our method in a simplified form, extending it to include the quantum effects which must be taken into account in the low-temperature regime. The resulting formulae show explicitly the transition from the quantum regime to the classical one as a function of temperature.

1 Introduction

In 1921 Keesom [1] pointed out that two molecules having permanent dipoles μ_A and μ_B and whose separation is R attract each other at temperature T according to the well-known formula [2]

E_{6} = - \frac{2 μ_{A}^{2} μ_{B}^{2}}{3 k T} \frac{1}{R^{6}} = - \frac{C_{6} (Keesom)}{R^{6}}

(1)

where k is the Boltzmann constant. This formula, which was derived well before London's work on the attraction between induced dipoles [3–7], shows a dependence of the interaction on R⁻⁶ with a C₆ coefficient which depends on T, decreasing in importance with increasing temperature. However, the Keesom formula looses its validity at low temperatures and for molecules having large values of μ, i.e. for large values of the modulus of the dimensionless T-dependent factor

a = - \frac{μ_{A} μ_{B}}{R^{3} k T} .

(2)

Corrections to the Keesom C₆ coefficients for several molecular systems at room temperature were derived by us [8] using the complete series expansion for the angle-dependent exponential function occurring in the Keesom integral K(a)

K (a) = \int_{Ω} ⅆ Ω \exp [a F (Ω)]

(3)

with

F (Ω) = \sin θ_{A} \sin θ_{B} \cos ϕ - 2 \cos θ_{A} \cos θ_{B} = - \frac{V (Ω)}{a k T}

where Ω = (θ_A, θ_B, ϕ) is the angle describing the mutual orientations of dipoles μ_A and μ_B, dΩ the (internal) configurational weight, and V(Ω) the potential of interaction of the two dipoles.

In the same paper [8], we extended the calculation to the evaluation of Keesom coefficients up to the R⁻¹⁰ term for small values of the parameter a.

For large values of a, which include the low-temperature regime, an asymptotic series expansion in (−a)⁻ⁿ was also derived by us [9], the leading terms giving for the Keesom integral the simple formula

K_{∞} (a) \sim \frac{4 π}{3} \frac{\exp [- 2 a]}{a^{2}} (1 - \frac{2}{3 a}) .

(4)

The resulting expansion in terms of the inverse of the interaction parameter a was then proved by us [10] to be an asymptotic representation, being subsequently rederived by Abbott [11] using standard integral tables [12] and the Laplace expansion theorem for asymptotic integrals.

We propose here to extend our method of calculation, reformulated in a considerably simplified form by using a modified coordinate transformation, to include the quantum effects, which must be taken into account in the low-temperature regime. The resulting formulae show the transition from the quantum regime to the classical one as a function of temperature. Our new approach takes as a starting point the oscillations of the dipoles around their equilibrium positions [9] by expanding the interaction potential in powers of the deviation angles from the axis joining the centres of attraction. The measure for the calculation of averages is taken that of the harmonic oscillator, corrected for the curvature, all the resulting integrals being over gaussian functions. The quantum density matrix for this model [13a] can be constructed so that the averages in the quantum regimes are straightforward. Taking into account the leading non-linear terms in the potential energy necessitates evaluation to first order of the perturbed density matrix, which could be obtained by many-body theoretical techniques [13a,14].

2 The coordinate system

We introduce the coordinate system ξ, η, ζ defined as follows: let ζ denote the axis originating in O, which is the location of the point-like dipole A, and pointing from A to B. Dipole A is situated in the plane (Aζ), containing the ζ axis, and similarly it is defined the plane (Bζ). Let ϕ denote the angle between the planes (Aζ) and (Bζ), rotating from A toward B. Further we consider the plane orthogonal to the ζ-axis in the point O, which is the location of dipole A, and introduce a pair of mutually orthogonal axis ξ and η inside this plane, so as to make the coordinate system (ξ, η, ζ) right-handed. Then θ_A and θ_B are the polar angles of the corresponding dipoles measured with respect to the ζ-axis. Then are defined, in analogy with Ref. [9], the components of the vectors

{\begin{cases} θ_{A} = θ_{A}^{(ξ)}, θ_{A}^{(η)} \\ θ_{B} = θ_{B}^{(ξ)}, θ_{B}^{(η)} \end{cases}

(5)

from the angles θ_A and θ_B, so as to make Eq. (33) of Ref. [9] satisfied, for all values of the variables in the interval (0, π), viz.

{\begin{cases} θ_{A} = \sqrt{θ_{A}^{{(ξ)}^{2}} + θ_{A}^{{(η)}^{2}}} \\ θ_{B} = \sqrt{θ_{B}^{{(ξ)}^{2}} + θ_{B}^{{(η)}^{2}}} \end{cases}

(6)

and moreover

\cos ϕ = \frac{θ_{A}^{(ξ)} θ_{B}^{(ξ)} + θ_{A}^{(η)} θ_{B}^{(η)}}{θ_{A} θ_{B}}

(7)

therefore obviously

{\begin{cases} θ_{A}^{(ξ)} = θ_{A} \cos φ_{A} \\ θ_{B}^{(ξ)} = θ_{B} \cos φ_{B} \end{cases}, {\begin{cases} θ_{A}^{(η)} = θ_{A} \sin φ_{A} \\ θ_{B}^{(η)} = θ_{B} \sin φ_{B} \end{cases}

(8)

ϕ = φ_{B} - φ_{A} .

The postulated variable transformation is canonical; therefore it conserves the measure in phase space. It would therefore be possible to compute the phase space measure in the transformed variables, in order to evaluate the Keesom integral. However, we compute the configurational measure in the three dimensional Euclidean space of each dipole, which is proportional to the solid angle spanned by it, after factorizing out the integral over Euclidean momenta. In the linear approximation, for small θ_A and θ_B, there follows

ⅆ φ_{A} ⅆ Ω = ⅆ θ_{A}^{(ξ)} ⅆ θ_{A}^{(η)} ⅆ θ_{B}^{(ξ)} ⅆ θ_{B}^{(η)} = θ_{A} θ_{B} ⅆ θ_{A} ⅆ θ_{B} ⅆ ϕ ⅆ φ_{A}

(9)

where dφ_A accounts for the rotational symmetry as a whole around the ζ-axis.

By comparison of Eq. (9) with Eq. (3) it results that

ⅆ φ_{A} ⅆ Ω = \sin θ_{A} \sin θ_{B} ⅆ θ_{A} ⅆ θ_{B} ⅆ ϕ ⅆ φ_{A} = \frac{\sin θ_{A} \sin θ_{B}}{θ_{A} θ_{B}} ⅆ θ_{A}^{(ξ)} ⅆ θ_{A}^{(η)} ⅆ θ_{B}^{(ξ)} ⅆ θ_{B}^{(η)}

(10)

where, on account of symmetry, the reduced domain of variation may be considered:

{\begin{cases} 0 \leq θ_{A} \leq \frac{π}{2} \\ 0 \leq θ_{B} \leq π . \end{cases}

(10

3 The average potential energy in the classical regime

In the classical regime the kinetic energy can be factorized in the evaluation of the thermodynamic potential, yielding a constant factor independent of the configuration. The potential energy function may be expanded yielding to second order

F (Ω) = θ_{A}^{{(ξ)}^{2}} + θ_{A}^{{(η)}^{2}} + θ_{B}^{{(ξ)}^{2}} + θ_{B}^{{(η)}^{2}} + θ_{A}^{(ξ)} θ_{B}^{(ξ)} + θ_{A}^{(η)} θ_{B}^{(η)} - 2 .

(11)

This potential function may be diagonalized by the linear orthogonal coordinate transformation

{\begin{cases} α = θ_{A} + θ_{B} \\ β = - θ_{A} + θ_{B} \end{cases}

(12)

The resulting potential function expanded to fourth order leads to

F (Ω) = \frac{3}{4} α^{2} + \frac{1}{4} β^{2} - \frac{1}{16} α^{4} - \frac{1}{48} β^{4} - \frac{1}{12} {| α \times β |}^{2} - 2 + h.o.t .

(13)

The Boltzmann average potential energy is therefore obtained through Eqs. (9) and (10) as

〈 V (Ω) 〉 = \frac{\int_{Ω} ⅆ Ω V (Ω) \exp {a F (Ω)}}{K (a)} = \frac{2 \int_{0}^{\frac{π}{2}} ⅆ θ_{A} \sin θ_{A} \int_{0}^{π} ⅆ θ_{B} \sin θ_{B} \int_{0}^{2 π} ⅆ φ V (Ω) \exp {a F (Ω)}}{K (a)}

(14)

where the factor 2 accounts for the equivalent configurations with reversed dipoles.

Expanding in inverse powers of a, and by keeping terms up to O(1/a³) while neglecting the exponentially damped terms arising from the upper limits of integration, Eq. (14) yields, through the change of variables (12)

\begin{array}{l} 〈 V (Ω) 〉 \sim - \frac{a k T \exp {- 2 a}}{4 π K (a)} \int_{- ∞}^{+ ∞} ⅆ α^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ α^{(η)} \int_{- ∞}^{+ ∞} ⅆ β^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ β^{(η)} \times [1 - \frac{1}{12} (α^{2} + β^{2})] \exp {\frac{3}{4} a α^{2} + \frac{1}{4} a β^{2}} \times [- 2 + \frac{3}{4} α^{2} + \frac{1}{4} β^{2}] - \frac{a k T \exp {- 2 a}}{4 π K (a)} \int_{- ∞}^{+ ∞} ⅆ α^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ α^{(η)} \int_{- ∞}^{+ ∞} ⅆ β^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ β^{(η)} \exp {\frac{3}{4} a α^{2} + \frac{1}{4} a β^{2}} \\ \times [- \frac{1}{16} α^{4} - \frac{1}{48} β^{4} - \frac{1}{12} {| α \times β |}^{2}] - \frac{a^{2} k T \exp {- 2 a}}{4 π K (a)} \int_{- ∞}^{+ ∞} ⅆ α^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ α^{(η)} \int_{- ∞}^{+ ∞} ⅆ β^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ β^{(η)} \exp {\frac{3}{4} a α^{2} + \frac{1}{4} a β^{2}} \times [- \frac{1}{16} α^{4} - \frac{1}{48} β^{4} - \frac{1}{12} {| α \times β |}^{2}] \times [- 2 + \frac{3}{4} α^{2} + \frac{1}{4} β^{2}] + h.o.t = - a k T (- 2 - \frac{2}{a} + \frac{2}{3 a^{2}}) + O (\frac{1}{a^{2}}) \end{array}

(15)

where, up to O(1/a³):

\begin{array}{l} K (a) \sim - \frac{a k T \exp {- 2 a}}{4 π} \int_{- ∞}^{+ ∞} ⅆ α^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ α^{(η)} \int_{- ∞}^{+ ∞} ⅆ β^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ β^{(η)} [1 - \frac{1}{12} (α^{2} + β^{2})] \\ \times \exp {\frac{3}{4} a α^{2} + \frac{1}{4} a β^{2} - \frac{3}{48} a α^{4} - \frac{1}{48} a β^{4} - \frac{1}{12} a {| α \times β |}^{2} - 2 a} + O (\frac{1}{a^{4}}) \\ = - \frac{a k T}{4 π} \times \exp {- 2 a} \int_{- ∞}^{+ ∞} ⅆ α^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ α^{(η)} \int_{- ∞}^{+ ∞} ⅆ β^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ β^{(η)} [1 - \frac{1}{12} (α^{2} + β^{2}) - \frac{1}{16} a α^{4} - \frac{1}{48} a β^{4} - \frac{1}{12} a {| α \times β |}^{2}] \times \exp {\frac{3}{4} a α^{2} + \frac{1}{4} a β^{2}} + O (\frac{1}{a^{4}}) \\ = \frac{\exp {- 2 a}}{a^{2}} \frac{4 π}{3} (1 - \frac{2}{3 a}) + O (\frac{1}{a^{4}}) . \end{array}

(15

This result coincides with that obtained by us in Ref. [9] by a more laborious procedure.

4 Quantum evaluation of the Keesom integral

In the limit of a very small temperature the average in Eq. (14) should be evaluated using quantum statistics [13]

〈 \hat{V} (Ω) 〉 = \frac{tr 〈 \hat{V} (Ω) \hat{ρ} (Ω) 〉}{tr 〈 \hat{ρ} (Ω) 〉}

(16)

where

\hat{ρ} (Ω)

is the quantum statistical operator, or density matrix, at the temperature T. Eq. (16) must be used in the place of Eq. (14) when the temperature is comparable or smaller than the energy differences between the lowest levels of the oscillators. In this case, only a few of the lowest levels are effectively occupied, so that the displacements of the coordinates from the equilibrium positions are small in the mean. The velocity v of a bead placed in the extremity of a dipole (A or B) oscillating around its equilibrium position along the ζ-axis is

v^{2} = σ^{2} ({\dot{θ}}^{2} + \sin^{2} θ {\dot{φ}}^{2})

(17)

\dot{θ} = \frac{ⅆ}{ⅆ t} | θ |

(17

σ being the distance from the centre of oscillation of the dipole. For small θ the kinetic energy of that extremity is written

\frac{1}{2} m v^{2} ≅ \frac{1}{2} m σ^{2} ({\dot{θ}}^{2} + θ^{2} {\dot{φ}}^{2}) = \frac{1}{2} I ({\dot{θ}}^{{(ξ)}^{2}} + {\dot{θ}}^{{(η)}^{2}})

(18)

where I = mσ² is a moment of inertia. The rhs of Eq. (18) is a notation of more general validity because it includes an arbitrary linear distribution of masses in the system [15]. Consequently, the system of coupled dipoles behaves like an assembly of four coupled harmonic oscillators, whose average potential energy would amount classically to 2kT, being independent of the frequencies of oscillation (see Eq. (15)).

It is then possible to write the total energy operator of the system, putting I_A = I_B = I, μ_A = μ_B = μ

{\hat{H}}_{0} (Ω, R) = \frac{1}{2} I ({\dot{θ}}_{A}^{{(ξ)}^{2}} + {\dot{θ}}_{A}^{{(η)}^{2}} + {\dot{θ}}_{B}^{{(ξ)}^{2}} + {\dot{θ}}_{B}^{{(η)}^{2}}) + \frac{μ^{2}}{R^{3}} (θ_{A}^{{(ξ)}^{2}} + θ_{A}^{{(η)}^{2}} + θ_{B}^{{(ξ)}^{2}} + θ_{B}^{{(η)}^{2}} + θ_{A}^{(ξ)} θ_{B}^{(ξ)} + θ_{A}^{(η)} θ_{B}^{(η)} - 2)

(19)

where Ω denotes here the complete configuration space of the system, including rotational symmetry as a whole. This Hamiltonian can be diagonalized by the variable transformation (12), so as to give the separable Hamiltonian

{\hat{H}}_{0} (Ω, R) = \frac{1}{4} I ({\dot{α}}^{{(ξ)}^{2}} + {\dot{α}}^{{(η)}^{2}} + {\dot{β}}^{{(ξ)}^{2}} + {\dot{β}}^{{(η)}^{2}}) + \frac{3}{4} \frac{μ^{2}}{R^{3}} (α^{{(ξ)}^{2}} + α^{{(η)}^{2}}) + \frac{1}{4} \frac{μ^{2}}{R^{3}} (β^{{(ξ)}^{2}} + β^{{(η)}^{2}}) - 2 \frac{μ^{2}}{R^{3}}

(20)

where the oscillation frequencies are

{\begin{cases} ω_{α}^{2} = \frac{3 μ^{2}}{R^{3} I} \\ ω_{β}^{2} = \frac{μ^{2}}{R^{3} I} . \end{cases}

(21)

According to the Gibbs canonical distribution, at equilibrium with temperature T, the density matrix in the coordinate representation has elements

ρ (Ω, Ω^{'}) = \sum_{n} \sum_{n^{'}} 〈 n | \exp {- \frac{\hat{H}}{k T}} n^{'} 〉 ψ_{n} (Ω) ψ_{n^{'}}^{*} (Ω^{'})

(22)

where, if

\hat{H}

is diagonal in the {n} basis,

〈 n | \hat{H} n^{'} 〉 = ɛ_{n} δ_{n, n^{'}}

(23)

then

ρ (Ω, Ω^{'}) = \sum_{n} \exp {- \frac{ɛ_{n}}{k T}} ψ_{n} (Ω) ψ_{n}^{*} (Ω^{'}) .

(24)

Expression (24) may be evaluated for a linear harmonic oscillator whose frequency is ω₀, and the energy levels

ɛ_{n} = ℏ ω_{0} (n + \frac{1}{2}) n = 0, 1, 2, \dots

(25)

thus obtaining [13a] the normalized diagonal density elements

ρ_{0} (α^{(ξ)}, α^{(ξ)}) = {(\frac{I ω_{α}}{2 π ℏ} \tan h \frac{ℏ ω_{α}}{2 k T})}^{\frac{1}{2}} \times \exp {- α^{{(ξ)}^{2}} \frac{I ω_{α}}{2 ℏ} \tanh \frac{ℏ ω_{α}}{2 k T}}

(26)

and the similar expressions for the remaining variables appearing in Eq. (20). There follows that

\begin{array}{l} ρ_{0} (Ω, Ω) = ρ_{0} (α^{(ξ)}, α^{(ξ)}) ρ_{0} (α^{(η)}, α^{(η)}) ρ_{0} (β^{(ξ)}, β^{(ξ)}) ρ_{0} (β^{(η)}, β^{(η)}) \\ = (\frac{I ω_{α}}{2 π ℏ} \tan h \frac{ℏ ω_{α}}{2 k T}) \exp {- α^{2} \frac{I ω_{α}}{2 ℏ} \tan h \frac{ℏ ω_{α}}{2 k T}} \times (\frac{I ω_{β}}{2 π ℏ} \tan h \frac{ℏ ω_{β}}{2 k T}) \exp {- β^{2} \frac{I ω_{β}}{2 ℏ} \tan h \frac{ℏ ω_{β}}{2 k T}} . \end{array}

(27)

It is therefore possible to replace Eq. (15) by

{〈 \hat{V} (Ω) 〉}_{0} \sim \exp {- 2 a} \int_{- ∞}^{+ ∞} ⅆ α^{(ξ)} \int_{- ∞}^{+ ∞} ⅆ α^{(η)} \int_{- ∞}^{+ ∞} ⅆ β^{(ξ)} \times \int_{- ∞}^{+ ∞} ⅆ β^{(η)} \times ρ_{0} (Ω, Ω) \times (- a k T) (- 2 + \frac{3}{4} α^{2} + \frac{1}{4} β^{2}) .

(28)

Therefore

{〈 \hat{V} (Ω) 〉}_{0} \sim \frac{1}{2} ℏ ω_{α} \coth \frac{ℏ ω_{α}}{2 k T} + \frac{1}{2} ℏ ω_{β} \coth \frac{ℏ ω_{β}}{2 k T} + 2 a k T .

(29)

In the limiting case T → ∞, or equivalently

ℏ \to 0

, the coth terms in Eq. (29) may be expanded into powers of

ℏ

and yield, when combined with Eq. (15), the leading terms of the low temperature expansion with quantum corrections for

ℏ^{2} | a | / 3 I

small

〈 \hat{V} (Ω) 〉 \sim 2 a k T + 2 k T + \frac{ℏ^{2} (ω_{α}^{2} + ω_{β}^{2})}{12 k T} - \frac{2 k T}{3 a} + O (\frac{1}{a^{2}})

(30)

where the average non-linear terms were evaluated classically because they become important at the higher temperatures [13b].

The opposite limiting case T → 0 can be handled easily, though the resulting formula is not an analytic function of T in this limit, and yields obviously

{〈 \hat{V} (Ω) 〉}_{0} \sim 2 a k T + \frac{1}{2} ℏ (ω_{α} + ω_{β}) + ℏ ω_{α} \exp {- \frac{ℏ ω_{α}}{k T}} + ℏ ω_{β} \exp {- \frac{ℏ ω_{β}}{k T}} + \dots

(31)

In order to proceed further with quantum effects, it is required to evaluate the non-linear terms of the Hamiltonian function, which are both kinetic and potential energy terms, and the quantum corrections to the diagonal elements of the density matrix in the coordinate representation. In order to evaluate these to first order in the perturbation

\hat{W}

, it is sufficient to calculate

ρ^{(1)} (Ω, Ω) = \sum_{n} 〈 n | \exp {- \frac{{\hat{H}}_{0} + \hat{W}}{k T}} | n 〉 ψ_{n} (Ω) ψ_{n}^{*} (Ω) ≅ \sum_{n} 〈 n | \exp {- \frac{{\hat{H}}_{0}}{k T}} (1 - \frac{\hat{W}}{k T}) | n 〉 ψ_{n} (Ω) ψ_{n}^{*} (Ω)

(32)

because the changes in the eigenvectors

| n >

and the corresponding eigenfunctions

ψ_{n}^{*} (Ω)

would contribute to second order only. The necessary first order corrections to the density matrix can be evaluated from Eq. (32). Calculation of the quantum corrections to the third term of the asymptotic expansion (Eq. (15)) is in progress.

5 Conclusions

In this paper we studied the small oscillations of two dipoles which are coupled through electromagnetic interaction, both classically and quantum mechanically. In the low-temperature regime in which the dipoles perform small oscillations around their equilibrium positions it is convenient to introduce a coordinate system which reduces, in the limit of small oscillations, to a system of four coupled harmonic oscillators, whose density matrix can be evaluated exactly in the quantum as well as in the classical regime. The position of each dipole is represented by the distance of the vertex of the arrow from a pole of the sphere, measured over the spherical surface, and then referred to an orthogonal coordinate system inside the plane tangent to the sphere in the same pole. The average of successive terms of this series expansion yields, as it was proved in Ref. [10], an asymptotic expansion of the averaged interacting potential in inverse powers of the parameter a.

The model allows one to evaluate non-linear effects which are due to the curvature of the configuration space and to higher potential energy terms. Since these affect mainly the higher quantum states, they can be evaluated classically [13b].

Eq. (30) shows the effect of the lowering of the temperature so as to make the quantum separation of levels, though small, not negligible compared to thermal fluctuations. The effect is a decrease of the static interaction in the points of equilibrium by a term which is written $ℏ^{2} | a | / 3 I$ (third term of Eq. (30), the quantum correction to the second term of the equation) and is proportional to the coefficient of the Laplace operator in the Schroedinger equation in the angle variables. This term can be interpreted as the effect of quantum fluctuations around equilibrium, averaged at the temperature T.

The condition for localization of quantum fluctuations [3,16,17] is that the distance between rotatory and oscillatory levels should be small toward the height of the barrier between the two symmetrical wells of potential energy, which is of the order of magnitude of ∣a∣kT [9], which gives

\frac{ℏ^{2}}{I} < < | a | k T .

(33)

Quantum corrections in Eqs. (29) and (30) are becoming important when

ℏ ω_{α}

is comparable to or smaller than kT [18], which means that

\frac{ℏ^{2}}{I} \approx \frac{k T}{| a |} .

(34)

Consequently, if this condition is satisfied inside some temperature interval and

| a | > > 1

, the condition (33) is automatically satisfied, though inequality (33) is actually independent of temperature.

Eq. (31) is expected to hold at very small temperatures and does not modify essentially the R-dependence of the interaction energy, since the leading term has the same R⁻³ dependence upon distance, like in the high-temperature classical regime. It is, however, noticeable that, quantum mechanically, also the operator kinetic energy of oscillation has a dependence upon angular coordinates.

Acknowledgments

Support by the Italian Ministry for Education University and Research (MIUR) under Grant No. 2006 03 0944 003, and by the University of Genoa is gratefully acknowledged.

Appendix Extension to infinity of the integration limits in Eqs. (15) and (15′)

It is noticed first that the integration domain Ω may be subdivided, owing to symmetries of the integrand, in such a way (see Fig. 1 in Ref. [9]) that the integral in Eq. (14) is performed over the region E defined by the inequalities

θ_{A} + θ_{B} \leq π, θ_{A} \geq 0, θ_{B} \geq 0 .

After expansion of the integrand in powers of θ_A, θ_B as shown in the main text, the domain of integration is extended by adding the region E of the plane (θ_A, θ_B) satisfying

θ_{A} + θ_{B} \geq π, θ_{A} \geq 0, θ_{B} \geq 0 .

The integral extended to the region E must therefore be shown to be negligible. From Eq. (11) it is obtained

θ_{A}^{2} + θ_{B}^{2} + θ_{A} \cdot θ_{B} - 2 \geq \frac{1}{2} (θ_{A}^{2} + θ_{B}^{2}) - 2 .

Therefore, if a < 0,

\exp {a (θ_{A}^{2} + θ_{B}^{2} + θ_{A} \cdot θ_{B} - 2)} \leq \exp {- 2 a + \frac{1}{2} a (θ_{A}^{2} + θ_{B}^{2})}

There follows the bound, for integer m, n,

\iint_{E} ⅆ θ_{A} ⅆ θ_{B} θ_{A}^{2 m} θ_{B}^{2 n} \exp {a (θ_{A}^{2} + θ_{B}^{2} + θ_{A} \cdot θ_{B} - 2)} \sin θ_{A} \sin θ_{B} \leq 2^{m + n} \exp {- 2 a} \iint_{E} ⅆ (\frac{1}{2} θ_{A}^{2}) ⅆ (\frac{1}{2} θ_{B}^{2}) {(\frac{1}{2} θ_{A}^{2})}^{m} {(\frac{1}{2} θ_{B}^{2})}^{n} \exp {\frac{1}{2} a (θ_{A}^{2} + θ_{B}^{2})} = 2^{m + n} \frac{\exp {- 2 a}}{{| a |}^{m + n + 2}} \iint_{E} ⅆ x ⅆ y x^{m} y^{n} \exp {- (x + y)} .

From the inequalities stated above follows that in E

θ_{A}^{2} + θ_{B}^{2} \geq \frac{1}{2} π^{2},

and, therefore,

x + y \geq 1 / 4 | a | π^{2}

. Hence

\begin{array}{l} 2^{m + n} \frac{\exp {- 2 a}}{{| a |}^{m + n + 2}} \iint_{E} ⅆ x ⅆ y x^{m} y^{n} \exp {- (x + y)} < 2^{m + n} \frac{\exp {- 2 a}}{{| a |}^{m + n + 2}} \int_{\frac{1}{4} | a | π^{2}}^{∞} ⅆ y y^{n} e^{- y} \int_{0}^{∞} ⅆ x x^{m} e^{- x} + 2^{m + n} \frac{\exp {- 2 a}}{{| a |}^{m + n + 2}} \int_{0}^{\frac{1}{4} | a | π^{2}} ⅆ y y^{n} e^{- y} \int_{\frac{1}{4} | a | π^{2} - y}^{∞} ⅆ x x^{m} e^{- x} \\ < C_{1} 2^{m + n} \frac{\exp {- 2 a}}{{| a |}^{m + n + 2}} {(\frac{1}{4} | a | π^{2})}^{n} \exp {- (\frac{1}{4} | a | π^{2})} + C_{2} 2^{m + n} \frac{\exp {- 2 a}}{{| a |}^{m + n + 2}} \exp {- (\frac{1}{4} | a | π^{2})} {(\frac{1}{4} | a | π^{2})}^{n + 1} {(\frac{1}{4} | a | π^{2})}^{m} = C_{1} 2^{m + n} {(\frac{1}{4} π^{2})}^{n} \frac{\exp {a (\frac{1}{4} π^{2} - 2)}}{{| a |}^{m + 2}} + C_{2} 2^{m + n} {(\frac{1}{4} π^{2})}^{m + n + 1} \frac{\exp {a (\frac{1}{4} π^{2} - 2)}}{| a |} \end{array}

where C₁, C₂ are suitable numerical constants, which can be estimated by evaluating the integrals.

^☆ Presented at the 4^èmes Journées Franco-Italiennes de Chimie (JFIC4), Nice, Sophia Antipolis, 17–18 April 2008.