text:
# Three-body interactions involving clusters and films

## I Introduction

The growing interest in problems involving interfaces has stimulated a large literature involving computations of the properties of such systems book ; davis . Most of these studies have assumed that the potential energy of the system is well approximated by a sum of two-body interactions between pairs of the constituent particles. In some of these instances, the pair potentials are not well known, so the implicit or explicit neglect of many-body interactions is appropriate, in view of the manifest quantitative limitations of these particular calculations. In some other cases, however, the pair potentials employed in the calculations have been constructed from empirical fits to properties of the relevant bulk systems; some authors have suggested that such potentials inherently include higher order, many-body contributions milen ; barker . Implicit in the latter point of view is the assumption that the many-body interactions in the given problem are of a similar form and magnitude to those in the bulk system. That assumption is evidently worth assessing in those cases where it can be tested. In this paper, we explore the importance of three-body van der Waals (VDW) interactions between a single atom and clusters of varying sizes and shapes, from nearly spherical to flat. In all cases, the atoms in the cluster are located at sites of a lattice of lattice constant $`a`$, like that shown in Fig. 1. Our study leads us to conclude that one should not expect three-body interactions to be properly reckoned when an empirical pair potential is employed.

Our analysis of the three-body energy is based on the Axilrod-Teller-Muto (ATM) expression for that energy AT ; M . For specificity, we consider the case when all three atoms are identical; then the interaction may be written:

$$V_{ATM}(𝐫_1,𝐫_2,𝐫_3)=\frac{9}{16}\frac{E_a\alpha _0^3}{r_{12}^3r_{23}^3r_{13}^3}(1+3\mathrm{cos}\theta _1\mathrm{cos}\theta _2\mathrm{cos}\theta _3)$$
(1)
In this expression, the three angles represent the inner angles of the triangle connecting the atoms (which has sides of length $`r_{12}`$, $`r_{13}`$ and $`r_{23}`$), the quantity $`\alpha _0`$ represents the static polarizability of the atom and $`E_a`$ is an energy (of order the ionization energy) characterizing the assumed frequency dependence of the polarizability at imaginary frequency:

$$\alpha (i\omega )=\frac{\alpha _0}{1+(\frac{\mathrm{}\omega }{E_a})^2}$$
(2)
Representative values of these parameters for some atoms appear in Table I.

The preceding ATM expression is a dispersion interaction, derived from third-order perturbation theory, ignoring exchange and overlap of the three atoms’ wave functions; thus, it is known to be valid at large separation, where “large” means relative to the atom’s hard-core interaction diameter. Empirical studies have often employed the interaction at closer approach but there is some disagreement about the validity of such applications say ; lucas ; bar1 ; bar2 ; proposals have been presented that modify the behavior at small separation damp . That important question is beyond the scope of the present work. The origin of the ATM interaction is not particularly intuitive. For that reason, among others, we discuss in Appendix A its origin in a particularly simple limit: the case when one atom (A) is far from the other two (B and C). In the limit when the characteristic frequency of A is much less than that of B and C, the three body (ATM) energy coincides with the mean electrostatic interaction between dipoles induced on B and C by the dipolar fluctuations of A.

The outline of this paper is the following. Section II presents a study of the VDW interactions between an atom and a monolayer film, from which we make a straightforward extension to include the case of an entire 3D crystal. Section III describes the three-body interactions between an atom and a cluster of atoms. It will be shown in both of these sections that the three-body energy $`V_3`$, relative to the two-body energy $`V_2`$, is proportional to the product $`n_s\alpha _0`$ of the number density of the solid and the polarizability (Table 1). Because of this dependence, we focus our attention on a dimensionless ratio

$$R=\frac{V_3}{V_2}\frac{1}{n_s\alpha _0}$$
(3)
$`R`$ depends very much on the specific geometry and can be either positive or negative. In particular, $`R`$ has a drastically different dependence on separation in the two problems investigated here. In the cluster case, the ratio decreases with distance (vanishing at large distance) while in the atom-solid case, the magnitude of $`R`$ increases with distance, asymptoting to a constant value that may be derived from the continuum theory of VDW interactions. For reference, we state the corresponding result in the case of the cohesive energy of inert gas solids; the ratio of three-body to two-body energies then is typically less than 10 per cent bar1 ; bar2 . This small value reflects in part a significant cancellation of ATM terms of opposite sign in a bulk material. The present study finds similar cancellation in some geometries, but not others. In the latter case, the three-body energy represents an important term that should be included in calculations of system properties. Section IV discusses the implications of this work.

## II Three-body atom-monolayer and atom-solid interactions

It was shown by Dzyaloshinskii, Lifshitz and Pitaevski DLP  that the long range interaction between a nonpolar atom dip  at position r=(0,0,z) and a semi-infinite solid, that occupies the half-space $`z<0`$, assumes the form

$$V(𝐫)\frac{C_3}{z^3}$$
(4)
$$C_3=\frac{\mathrm{}}{4\pi }𝑑\omega G(\omega )\alpha (i\omega )$$
(5)
Here $`G(\omega )`$, the surface-response function at imaginary frequency, is determined from the dielectric function $`\epsilon (i\omega )`$

$$G(\omega )=\frac{\epsilon (i\omega )1}{\epsilon (i\omega )+1}$$
(6)
Among the questions addressed in this section is the relative importance of contributions to $`C_3`$ from pair interactions and three-body interactions; we call these contributions $`C_3^{(2)}`$ and $`C_3^{(3)}`$, respectively. We first consider the problem of an atom situated at position r=(0,0,z) interacting with a square lattice ($`s\times s`$) of atoms that comprise a monolayer film, in the $`z=0`$ plane; its crystalline lattice directions are oriented along the $`x`$ and $`y`$ axes. In this case, the two-body van der Waals interaction may be written:

$$V_2=\underset{i}{}\frac{C_6}{r_i^6}$$
(7)
Here, the sum is over all atoms ($`i`$) of the film, each at distance $`r_i`$ from the external atom, and the coefficient of the asymptotic VDW pair interaction is

$$C_6=\frac{3}{4}E_a\alpha _0^2$$
(8)
This expression (first obtained by London london ) is derived from the frequency-dependent polarizability, Eq. (2). For $`s=60`$, Fig. 2 presents the dimensionless two-body potential

$$U_2=\frac{V_2}{ϵ_2};ϵ_2=E_a(\alpha _0n_s)^2$$
(9)
Similarly, we define

$$U_3=\frac{V_3}{ϵ_3};ϵ_3=ϵ_2\alpha _0n_s$$
(10)
Values of $`ϵ_2`$ and $`ϵ_3`$ are presented in Table I. One observes a significant dependence of the energy on lateral position only at small $`z/a`$, since at $`z>a`$ the film looks continuous. At large distance, one finds for $`V_2`$ a result that can be derived from an integration over the monolayer, which has a two-dimensional (2D) density $`1/a^2`$:

$$V_2(𝐫)\frac{C_{mono}^{(2)}}{z^4};C_{mono}^{(2)}=\frac{\pi }{2}\frac{C_6}{a^2}$$
(11)
The superscript (2) in the coefficient indicates that this energy originates from the pair interactions, i.e. Eq. (7). Also shown in Fig. 2 is the total three-body energy, obtained by summing ATM interactions between all pairs of atoms in the monolayer and the external atom (at r):

$$V_3(𝐫)=\underset{i<j}{}V_{ATM}(𝐫_i,𝐫_j,𝐫)$$
(12)
This energy is positive because the contribution of the repulsive terms in the sum exceed those of attractive terms. In Fig. 3 we show the ratio $`R`$ (Eq. (3)) for this case. One observes in the inset of this figure that at a distance $`z`$ much larger than the lattice constant $`a`$ (but smaller than the monolayer’s lateral extent) the ratio of these interactions is $`R_{monolayer}1.6`$. We will derive this result below. For distances larger than the size of the monolayer (equal to $`60a`$ in the case of the Fig. 3), the value of $`(V_3/V_2)`$ approaches a finite value that depends on the size of the monolayer. This value is discussed in the next section.

Now, let us address the problem of the three-body interaction in the case of an atom above a half-space solid, consisting of layers located at $`z=na`$, where $`n=0,1,2,\mathrm{}`$ First, we evaluate the two-body contribution by summing the monolayer energy over all layers. We note, from the Euler-Maclaurin theorem and Eq. (11), that the asymptotic dependence is

$$V_2(𝐫)C_{mono}^{(2)}\underset{n}{}\frac{1}{(z+na)^4}\frac{\pi }{6}n_s\frac{C_6}{(z\frac{a}{2})^3}$$
(13)
This leads to the result that the two-body contribution to the asymptotic potential of Eq. (4) involves a coefficient

$$C_3^{(2)}=\frac{\pi }{6}n_sC_6$$
(14)
We now evaluate the analogous three-body asymptotic coefficient, $`C_3^{(3)}`$. The expression is obtained from Eq. (12), with the sum of the ATM interactions taken over pairs of atoms throughout the half-space $`z0`$. To determine this quantity requires consideration of two-kinds of terms in the sum, leading to:

$$V_3=V_3^{inter}+V_3^{intra}$$
(15)
The interlayer term $`V_3^{inter}`$ represents the three-body contribution arising from pairs of solid atoms coming from different planes, while the intralayer term $`V_3^{intra}`$ involves pairs of atoms from the same plane. Fig. 4 shows the resulting $`R^{inter}`$ and $`R^{intra}`$ terms.

In Appendix B, we prove that the first term in Eq. (15) is identically zero at large separation, leaving just the second term. We note that the ratio of this surviving intralayer term to the analogous pair interaction term was calculated above for a single layer, as shown in Fig. 3; the ratio asymptotes to a constant value mentioned above. This means that, in the semi-infinite solid case, the same ratio applies to contributions from all planes, leading to exactly the same asymptotic ratio for the solid and the monolayer. Thus, asymptotically,

$$\left(\frac{V_3}{V_2}\right)_{solid}\frac{C_3^{(3)}}{C_3^{(2)}}1.6n_s\alpha _0$$
(16)
We now explain why this numerical result arises (and could have been anticipated without any numerical calculations). To do so, we employ the Clausius-Mossotti relation:

$$\frac{\epsilon (i\omega )1}{\epsilon (i\omega )+2}=\frac{4\pi }{3}n_s\alpha (i\omega )\xi (i\omega )$$
(17)
Assuming that the quantity $`\xi `$ is small (the case of small $`V_3/V_2`$), we may expand the relation for the asymptotic VDW coefficient, Eq. (5):

$$C_3=\frac{\mathrm{}}{4\pi }d\omega \alpha (i\omega )\left[\frac{3}{2}\xi (i\omega )(1\frac{\xi (i\omega )}{2}\mathrm{}.)\right]$$
(18)
$$C_3=C_3^{(2)}+C_3^{(3)}\mathrm{}.$$
(19)
The identification of successive terms with increasing number of interacting particles is physically plausible and is demonstrated here for the first two terms. If one inserts into these equations the Drude expression for the polarizability, Eq. (2), assuming the same species for external atom and solid atoms, one obtains

$$C_3^{(2)}=\frac{\mathrm{}}{2}n_s\alpha _0^2𝑑\omega \frac{1}{\left[1+(\frac{\mathrm{}\omega }{E_a})^2\right]^2}$$
(20)
$$C_3^{(3)}=\frac{\mathrm{}\pi }{3}n_s^2\alpha _0^3𝑑\omega \frac{1}{\left[1+(\frac{\mathrm{}\omega }{E_a})^2\right]^3}$$
(21)
The integrals may be evaluated, leading to values of the two relevant contributions to the dispersion coefficient:

$$C_3^{(2)}=\frac{\pi }{8}n_s\alpha _0^2E_a=\frac{\pi }{6}n_sC_6$$
(22)
$$C_3^{(3)}=\frac{\pi ^2}{16}n_s^2\alpha _0^3E_a$$
(23)
The ratio of coefficients results:

$$\frac{C_3^{(3)}}{C_3^{(2)}}=\frac{\pi }{2}n_s\alpha _0$$
(24)
This result is fully consistent with the behavior found numerically, with coefficient $`R=\pi /2`$ -1.6, Eq. (16). This confirms our attribution of the physical origin of the first two terms in the expansion of $`C_3`$.

We also explore the dependence of $`V_3`$ on the orientation of the monolayer. This is shown in Fig. 5. The result is that $`V_3`$ depends very strongly on the orientation, becoming attractive when the single atom and the monolayer are coplanar. The origin of this behavior will be explained in the next section.

## III Three-body interaction between a single atom and a cluster

Now, we address the problem of a single atom at a distance $`z`$ from a quasi-spherical cluster. By “quasi-spherical” we mean a cluster of $`N`$ atoms constructed from the closest $`N`$ points to the center of a simple cubic lattice of sites. The resulting ratio $`R`$ in this configuration is shown in Fig. 6 for different sizes and orientations of the cluster. There are at least three particularly interesting dependences seen in this figure. The first is that $`R`$ goes asymptotically to zero, in contrast to the nonzero constant value in the case of a monolayer at large distance. The second is the presence of a maximum in the magnitude of $`R`$. A third is that $`R`$ depends strongly on the orientation of the cluster. For example, in the case of the clusters with 7 and 33 atoms, the three-body interaction is attractive for $`\theta =0`$ but repulsive for large values of $`\theta `$. In the other two cases shown in this figure (clusters with 19 and 93 atoms), the trend is reversed.

We now explain why $`R`$ vanishes at large $`z`$ for quasi-spherical clusters. The ATM expression Eq. (1) at large $`z`$ is simplified because two of the three interatomic distances are essentially equal to $`z`$ and two of the interior angles of the three-atom triangle sum to $`\pi `$, while the third angle nearly vanishes. Then, the energy may be written

$$V_{ATM}(z)=\frac{C}{z^6}\underset{i<j}{}\frac{(13\mathrm{cos}^2\varphi _{ij})}{r_{ij}^3};C=\frac{9}{16}E_a\alpha _0^3$$
(25)
Here, the sum is over pairs of atoms in the cluster, separated by distance $`r_{ij}`$ and $`\varphi _{ij}`$ is the angle between their separation vector $`𝐫_{ij}`$ and the vector connecting the atoms to the external atom. Since the term in parentheses can vary between (-2) and 1, the sign of the resultant sum is not obvious. However, in a case where the separation vectors are randomly oriented, the average value of cos$`{}_{}{}^{2}\varphi `$ is 1/3. Similarly, for an infinite simple cubic lattice, the average is 1/3, as may easily be established (for any crystalline orientation). In either of these cases, therefore, the sum in Eq. (25) vanishes identically. In the case of a finite quasi-spherical cluster, the ATM sum is nonzero, due to residual surface contributions, but it is small compared to the two-body energy, which is proportional to the number of atoms in this large separation limit. Thus, it is expected that $`R`$ should be very small for symmetric clusters at large distances.

To understand the intriguing angular dependence in Fig. 6, we undertake additional calculations, as follows. We consider the long range interaction for two geometries shown in Fig. 7. In Fig. 7(a) appears the interaction for anisotropic clusters consisting of ($`s`$=1,2 or 3) rows of atoms, with $`N`$ atoms per row, oriented parallel to the $`z`$ axis, on which the atom is located. One observes for $`s=1`$ row that $`R`$ is positive ($`V_3<0`$) for all values of $`N`$. This behavior can be understood readily from Eq. (25) because $`\varphi _{ij}=0`$ for all pairs within the cluster. Hence, every term is negative and the sum increases essentially linearly with the number of nearest neighbor pairs (and hence $`N`$) for $`z/a>>N`$. Cases involving $`s>1`$ yield somewhat less attractive ATM sums, relative to two-body sums, as indicated in the $`s=2`$ and 3 curves in Fig. 7(a). To summarize, needle-like clusters, pointed toward the external atom, have negative ATM interactions, enhancing the two-body attraction.

What about more symmetric structures? In Fig. 7(b), we show results for systems consisting of atoms located at the sites of $`N`$ parallel squares of size $`s`$. The distance from the external atom to the first square is $`z=10s`$, and $`R`$ is plotted as a function of $`N/s`$. We see that for $`N/s<1`$ the 3-body contribution is repulsive, while for $`N/s>1`$ it is attractive. In the case $`N/s=1`$, that corresponds to a cubic configuration, the ratio $`R`$ is almost zero at this distance. Thus, the asymptotic behavior of the cube is analogous to that of the small clusters. We note from Fig. 7 that this vanishing is an “accident”, in that nonzero values of opposite signs occur for $`N>s`$ or $`N<s`$.

Now suppose, instead, that our system consists of atoms on sites of three parallel squares, as in Fig. 8. As expected, if the side of the squares is exactly 3, $`R`$ goes to zero, while the asymptotic value of $`R`$ becomes finite (and negative) for $`s>3`$. Note that this trend is consistent with the asymptotic value $`R3.1`$ for one layer, with $`s=60`$, shown in Fig. 3, and can be understood by the following argument. At large $`z`$ ($`sa`$), the film atoms are esentially equidistant form the external atom. Hence, the two-body energy is $`V_2=NC_6/z^6`$. As for the ATM energy, in the case of normal incidence at large $`z`$, the term involving 3 cosines in Eq. (1) vanishes, leaving:

$$V_{ATM}=\frac{9}{16}E_a\alpha _0^3\mathrm{\Sigma }$$
(26)
The sum $`\mathrm{\Sigma }`$ is over all pairs of atoms; for a film of very large extent, this may be written

$$\mathrm{\Sigma }=\underset{i<j}{}\frac{1}{r_{ij}^3}N_a\frac{k}{a^3};(N_a1)$$
(27)
Here, $`N_a`$ is the total number of atoms and $`k=4.7`$ is a numerical factor evaluated from the sum. Since $`C_6=(3/4)E_a\alpha _0^2`$, the ratio $`R`$ becomes

$$R=\frac{3}{4}k3.5$$
(28)
The large $`z`$, large $`s`$, limit of Fig. 3 approaches this value. For large $`z`$, but $`s`$ not very large compared to 1, the contribution of sites on the boundary of the film causes $`\mathrm{\Sigma }`$ to be less than $`Nk/a^3`$, so the ratio $`R>3.5`$ (as seen in Fig. 8). This is a significant reduction in magnitud for finite $`s`$ because the fraction of sites on the border of an $`s\times s`$ square is $`4(s1)/s^2`$.

These results lead to a relatively simple and comprehensive picture of the role of three-body interactions. Lines of atoms pointing toward the external atom lead to attraction, while planes of atoms perpendicular to the cluster-atom separation vector lead to repulsive interactions overall (see Appendix A). Note, in particular, that the explanation clarifies the peculiar angle-dependence in Fig. 6. In the process of making the clusters nearly spherical, those with $`N=7`$ and 33 end up with single atom prominences on their external ”faces”. The $`N=7`$ cluster consists of an atom at the origin and one along each of the +/- axes. The $`N=33`$ cluster is composed by a cube of side 3 plus an additional atom at the center of each face. In particular, when $`\theta =0`$ these configurations mimic the effect of a line of atoms pointing toward the external atom originating an attractive interaction as indicated by the dashed lines in the left panels of Fig. 6. Instead, the other two clusters with $`N=19`$ or 93 face the external atom with an almost complete plane of atoms perpendicular to the cluster-atom separation vector, leading to repulsive interactions at small angle $`\theta `$ (dashed lines in right panels of Fig. 6).

## IV Discussion and final remarks

The calculations in this paper use the Axilrod-Teller-Muto three-body interaction to supplement the contribution from two-body calculations of van der Waals dispersion forces. Specific calculations reveal that material geometry - both the morphology and orientation of clusters - is important, and that changing the orientation of a cluster can cause the overall VDW potential to change from attractive to repulsive. This is especially true for clusters with high polarizability.

Our calculations suggest three heuristics for predicting van der Waals dispersion interactions between an atom A and a material B (e.g., surfaces, clusters, or atoms) interacting by a force along a unit vector $`\widehat{e}`$: 1) If B has a large dimension perpendicular to $`\widehat{e}`$, this tends to produce a repulsive 3-body force with A, thus reducing the total dispersion attractions; 2) if B has a large dimension parallel to $`\widehat{e}`$ (i.e., depth), this tends to produce attractive 3-body forces with A, thus increasing the total dispersion attractions; 3) if B is quasi-spherical, this tends to decrease the magnitude of 3-body forces with A, meaning that the pairwise interaction sum gives a good approximation to the total interaction.

This work has produced other important conclusions as well. For example, in assessing 3-body interactions between an atom and a semi-infinite solid, only the intra-plane contributions to the 3-body effect are important; that is, distinct layers of the solid do not interact significantly with the external atom to produce 3-body ATM effects. In addition, analytical calculations show that the 3-body effect asymptotes to zero for quasi spherical clusters, as opposed to semi-infinite materials.

At least two interesting questions arise from this work. First, since we did not assess 4-body interactions, it remains to be seen whether 4-body interactions are of the same order of magnitude as 3-body interactions. Presumably such interactions scale as one higher power of the expansion parameter $`n_s\alpha _0`$. Second, the effect of coating layers on particles can be examined by a similar methodology. In 2-body effects only, the coating layer might not alter the VDW forces significantly; however, based on the present work, the 3-body forces could produce even repulsive net interactions.

###### Acknowledgements.

We are grateful to L.W. Bruch for stimulating and helpful discussions. This research has been supported by the National Science Foundation, the Ben Franklin Technology Center and the Environmental Protection Agency.

## Appendix A Nature of the ATM interaction at large distance

In this section we explore the microscopic nature of the ATM interaction by deriving its form (in a particular case) using a distinct alternative to the conventional approach to deriving the three-body interaction. The method described here provides a simple interpretation to this energy that is consistent with the general notion that such dispersion interactions arise from interactions between fluctuating dipole moments.

To be specific, we consider the situation when two atoms (B and C) are close together while the third atom (A) is located at large distance from the pair. Furthermore, for convenience, we make a specific assumption about the values of their characteristic energies that is different from that used to derive Eq. (1). The latter assumes that all three energies are the same; here, instead, we assume that $`E_A`$, the characteristic energy of atom A, is small compared to the energies (both $`E^{}`$) of B and C, meaning that A has a slowly fluctuating dipole moment. In this limit, we evaluate the triple dipole interaction by setting it equal to the electrostatic interaction between the two instantaneous dipoles (on B and C) that are induced by dipolar fluctuations of the third atom, A.

Assume that the external atom (A) is located at the origin and the B, C pair of atoms is located at distance $`z`$ on the z axis. A dipole moment fluctuation $`𝐩`$ on A will produce an electric field $`𝐄=[𝐩+3p_z\widehat{𝐳}]/z^3`$ at the neighboring positions of the B, C pair, where $`\widehat{𝐳}`$ is the unit vector in the z direction . Atoms B and C polarize instantaneously (in comparison with A), so that they develop an identical induced dipole moment $`𝐩_{induced}=\alpha _0𝐄`$. The electrostatic interaction energy between two such dipoles is then obtained in the usual way, resulting in:

$$E_{BC}=\frac{\alpha _0^2}{z^6}𝐩^2_A\frac{(13\mathrm{cos}^2\varphi _{BC})}{r_{BC}^3}$$
(29)
Here, we have have averaged over the time-dependent fluctuations of $`𝐩`$ on atom A, with the average denoted by brackets. The angle $`\varphi _{BC}`$ is that between the BC separation vector and the AB separation vector. Note that this energy is usually nonzero, varies as the inverse cube of the B-C separation and the inverse sixth power of the distance to A and has an interesting angle dependence. To explore the comparison with Eq. (1), we evaluate $`𝐩^2_A`$ within the Drude approximation fano , Eq. (2). The result is $`𝐩^2_A=(3/2)E_A\alpha _0`$. Hence, we arrive at the expression for the three-body interaction in this approximation:

$$V(z)=\frac{3}{2}E_A\frac{\alpha _0^3}{z^6}\frac{(13\mathrm{cos}^2\varphi _{BC})}{r_{BC}^3}$$
(30)
With one exception, this equation coincides with the result obtained from Eq. (1) in the present geometry (two angles summing to $`\pi `$ and the third vanishing). The exception is the value of the numerical coefficient, which is here (3/2) instead of (9/16) in Eq. (1). This difference arises because Eq. (1) is obtained from the ATM model for a case when the characteristic energies coincide, while Eq. (A1) assumes that the energies are very different, $`E_AE^{}`$. If this latter set of energies were inserted in the general ATM expression, it would yield the same (3/2) coefficient found here. Thus, we have identified a plausible and intuitive explanation (induced dipole interactions) of the origin of this ATM interaction and shown that the energy computed in the present limit agrees with ATM formula in the case of the same interacting atoms.

## Appendix B Interlayer three-body interaction

The total interlayer ATM interaction for an atom above a half-space crystal may be written as a sum of terms in which the external atom ($`0`$) is fixed in position and one sums over all pairs of internal atoms. We evaluate the sum by considering a particular individual atom ($`i`$) within the solid and summing over atoms in all layers below it, which we represent as a continuum. McLachlan evaluated the many-body interaction for precisely this geometry mc1 ; mc2 , depicted in Fig. 9. Atom $`i`$ lies at distance $`b`$ (of order a lattice constant) above the continuum while the external atom lies at distance $`h=z+na+b`$ above this continuum; here $`n`$ is the layer number for atom $`i`$. The atoms labelled $`0`$ and $`i`$ are separated laterally by a distance $`R`$. McLachlan’s resulting expression for this many-body interaction energy involves two coefficients:

$$C_{s1}=\frac{3\mathrm{}}{\pi }\xi (i\omega )\alpha ^2(i\omega )𝑑\omega $$
(31)
$$C_{s2}=\frac{3\mathrm{}}{\pi }\xi ^2(i\omega )\alpha ^2(i\omega )𝑑\omega $$
(32)
Since the first of these coefficients is proportional to the third power of the polarizability, it represents a contribution to the three-body energy explored in this paper. The $`C_{s2}`$ coefficient is of fourth order in the polarizability, as discussed explicitly below. Corresponding to these coefficients, the expression of McLachlan yields a limiting behavior, for $`h>>b`$, as follows:

$$V_M=C_{s1}U_1(R,z)+C_{s2}\frac{1}{(R^2+h^2)^3}$$
(33)
$$U_1(R,z)=\frac{1}{(R^2+h^2)^3}\left[\frac{4}{3}\frac{2h^2}{R^2+h^2}\right]$$
(34)
Note that the $`U_1`$ term is negative for $`R/h<2^{1/2}`$ and positive for larger $`R/h`$, with a maximum at $`R=h`$. Hence, there tends to be cancellation between these regimes of distance. To determine the total contribution from all of the atoms in layer $`n`$, we need to sum this expression for $`V_M`$ over a monolayer of atoms; for $`h>>b`$, we can achieve this by integrating this term over all positions $`\stackrel{}{R}`$ in the plane of atom $`i`$. Note that

$$U_1(R,z)d^2\stackrel{}{R}=0$$
(35)
We conclude that the three-body contribution from interlayer interactions vanishes at large distance. The remaining interlayer contribution is fourth order, which is not considered in this paper. For completeness, however, we evaluate this contribution from layer $`n`$, of two-dimensional density $`1/a^2`$:

$$V_{4n}=\frac{1}{a^2}\frac{C_{s2}}{(R^2+h^2)^3}d^2\stackrel{}{R}=\frac{\pi }{2}\frac{C_{s2}}{a^3}\frac{1}{h^4}$$
(36)
We may then integrate this contribution over all of the layers ($`n`$), yielding

$$V_{4inter}=\frac{\pi }{6}\frac{C_{s2}}{a^2}\frac{1}{z^3}$$
(37)
The ratio of this term to the two-body term is found to be

$$\frac{V_{4inter}}{V_2}=\frac{C_{s2}}{C_6}=\frac{5}{8}\left(\frac{2\pi \alpha _0}{a^3}\right)^2$$
(38)
This is smaller by a factor $`\alpha _0/a^3`$ than the corresponding ratio in Eq. (14), for the three-body intralayer term. We note, however, that we have not considered fourth order intralayer terms, which might be larger and/or of opposite sign.


format: text

text:
# On the inner structure of confined positronium

## I Introduction

Shortly after the discovery of the positron (e+) Anderson (1933), it was recognized that e+ may form an ‘atom’ (e+e-) which, from the point of view of spectroscopy, may be looked upon as an ultralight isotope of hydrogen Mohorovičić (1934). With the same right, (e+e-) may be considered as an ultralight isotope of antihydrogen. Together with the photon, the Majorana neutrino, and the electrically neutral $`\pi `$-meson ($`\pi ^0`$), (e+e-) thus belongs to those ‘particles’ that are their own ‘antiparticles’. In 1945 (e+e-), at that time still hypothetical, was named positronium Ruark (1945), later to be equipped with the ‘chemical’ symbol Ps McGervey and DeBenedetti (1959). At about the same time its main properties were established Wheeler (1946). Among these are that Ps is rather short-lived due to the particle–antiparticle annihilation and that its ground state splits into two substates, a singlet state l1S0 and a triplet state l3S1.

In the l1S0 state, Ps is called parapositronium (p-Ps), and in the l3S1 state, orthopositronium (o-Ps). Under most circumstances, the energy excess of o-Ps over p-Ps, $`\mathrm{\Delta }E_{\mathrm{o}\mathrm{p}}=0.8412\times 10^4`$ eV, is negligibly small, hence the ratio of the production rates of o-Ps and p-Ps is usually equal to that of the statistical weights of the substates, namely 3:1. The most important differences between o-Ps and p-Ps result from the different annihilation modes of these two substates, which, in turn, are consequences of the conservation of energy and angular momentum. The self-annihilation of p-Ps results in the generation of two $`\gamma `$-quanta of equal energies $`E_\gamma =m_\mathrm{e}c^2=0.511`$ MeV, where $`m_\mathrm{e}`$ is the electron mass and $`c`$ the speed of light in vacuum, whereas the self-annihilation of o-Ps requires the generation of, at least, three $`\gamma `$-quanta. Owing to the smallness of the coupling constant of quantum electrodynamics ($`\alpha 1/137`$) and to various statistical factors, the ratio between the 2$`\gamma `$ annihilation rate, $`\mathrm{\Gamma }_{2\gamma }=7.99\times 10^9`$ s-1, and the 3$`\gamma `$ annihilation rate, $`\mathrm{\Gamma }_{3\gamma }=7.04\times 10^6`$ s-1, exceeds $`10^3`$ and is therefore highly significant.

The preceding statements refer to Ps in vacuum or in extremely dilute gases. In dense gases and in condensed matter the fate of o-Ps is radically changed, since now the e+–e- annihilation may (and in general does) involve electrons other than the 1s electrons of the Ps atoms. This leads to 2$`\gamma `$ annihilation by the ‘pick-off’ mechanism. This mechanism may be visualized as follows. An electron in the neighbourhood with spin direction opposite to that of the e+ ‘picks off’ the positron from a Ps atom and annihilates it by a 2$`\gamma `$ reaction. The pick-off reaction rate $`\mathrm{\Gamma }_{\mathrm{po}}`$ is determined by the probability density of the electrons with opposite spin at the e+ location. This density is clearly less than that of electrons in the 1s orbit of Ps, in condensed matter typically by about a factor of ten. Hence $`\mathrm{\Gamma }_{\mathrm{po}}`$ is smaller than the self-annihilation rate $`\mathrm{\Gamma }_{2\gamma }`$ but still large compared to $`\mathrm{\Gamma }_{3\gamma }`$. This has several important consequences:

* Positron annihilation in condensed matter takes place virtually exclusively by the generation of two $`\gamma `$ quanta.
* The annihilation rate of o-Ps depends significantly on the environment in which the positronium atom is located.
* Because of (ii), the o-Ps annihilation rate may be used as indicator of the preferred sites of Ps in condensed matter.
* In Ps-forming solids containing imperfections that attract Ps, the site-sensitivity of $`\mathrm{\Gamma }_{\mathrm{po}}`$ may be used to detect and monitor such imperfections.

The pick-off annihilation rates of o-Ps fall into a range that allows accurate measurements of their inverses, the o-Ps lifetimes Hautojärvi (1979). Often it is possible to perform e+ lifetime spectrometry, i.e., to distinguish different positron annihilation sites by the different e+ lifetimes associated with them. In order to gain a theoretical understanding of the site dependence of the e+ lifetimes, one has to calculate the wavefunctions of Ps atoms that are confined to the various imperfections that are capable of trapping Ps. This has turned out to be a very formidable task.

The electrostatic forces experienced by electrons and positrons are equal and opposite. In most situations the variation of the electrostatic potentials over the diameter of a Ps atom (about two Bohr radii, see below) is negligibly small compared with the positronium binding energy $`E_{\mathrm{free}\mathrm{PS}}`$ (in vacuum about 6.8 eV). Unless a Ps atoms enters a kind of hydrogen bond, the electrostatic forces on its constituents cancel to an excellent approximation. Then the interaction of Ps with its condensed-matter environment is entirely due to the exchange interaction between the Ps-electrons and the host electrons since for e+ implanted into matter the Pauli principle is irrelevant. The fact that the ensuing repulsion of Ps by the host atoms is short-ranged suggests that the confinement of Ps to interstices or locations where atoms are missing may be represented by a rectangular potential well at these sites.

Since even this model is still a nontrivial quantum mechanical problem, Tao and Eldrup initiated a further simplified model in which the positronium is considered as a point particle in a spherical infinite potential well, being in the ground state Tao (1972); Chuang and Tao (1973); Eldrup et al. (1981), extending the ideas of Ferrell Ferrell (1957) and Brandt et al. Brandt et al. (1960). In this model, pick-off annihilation happens in a layer $`r_0rR_0`$ near the surface of the spherical volume of radius $`R_0`$, the layer representing the medium outside the hole. The width of this layer is a parameter fitted experimentally. Later, this model was refined by considering a potential wall of finite height Nakanishi and Jean (1988) and of unsharp shape Mukherjee et al. (1998). In these modifications, the pick-off annihilation of the positron occurs in the region $`R_0r<\mathrm{}`$, which seems to be a better approximation from the physical point of view and provides a better fit to the experimental data. Nevertheless, as stressed by Mukherjee et al. Mukherjee et al. (1998), this approach is still an oversimplification since it treats the positronium as a single point particle.

In this paper, we provide an approach that avoids this simplification. We determine the quantum mechanical ground state of the positron-electron two-particle system with an infinite potential well that acts only on the electron, in addition to the Coulomb attraction between the positron and the electron. Hence, we give account of the arising nontrivial internal structure of a confined positronium.

Due to the complicated nature of the problem at hand, we calculate the ground state numerically, in terms of an expansion with respect to an appropriate set of basis functions (truncated to contain the first 216 basis functions) and via a variational method, for various pore sizes. The probability of finding the positron outside the void and the overlap of the positron and the electron inside the hole are also determined. We analyse and interpret the results, and discuss how further refinements of the model and the explanation of the observed temperature dependence of the pick-off rate can be implemented.

It is to be noted that a related investigation was done by Sommerfeld and Welker Sommerfeld and Welker (1938), where the ground-state energy and wave function of a hydrogen atom in a hole is calculated. In that study, the proton is considered to be infinitely heavy and to reside in the center of the spherical hole, and the electron is confined to the cavity by an infinite potential well. We make a comparison between those results and our ones, and provide a simple common explanation of the behavior of both systems.

## II The model

Using the nonrelativistic framework and observing that, under the present conditions, spatial and spin degrees of freedom are decoupled, we formulate our model expounded above with a normalized electron–positron two-particle scalar wave function $`\mathrm{\Psi }(𝐫_\mathrm{e},𝐫_\mathrm{p})`$, the Hamiltonian

$$=\frac{\mathrm{}^2}{2m_\mathrm{e}}(_{𝐫_\mathrm{e}}^2+_{𝐫_\mathrm{p}}^2)\frac{e^2}{4\pi \epsilon _0}\frac{1}{\left|𝐫_\mathrm{e}𝐫_\mathrm{p}\right|},$$
(1)
and the boundary condition

$$\mathrm{\Psi }(𝐫_\mathrm{e},𝐫_\mathrm{p})|_{r_\mathrm{e}=R_0}=0,𝐫_\mathrm{p}.$$
(2)
Here, $`𝐫_\mathrm{e}`$ denotes the electron co-ordinates and $`𝐫_\mathrm{p}`$ the positron co-ordinates, $`R_0`$ is the radius of the potential well to which the electron is confined, $`e`$ is the elementary charge, $`\epsilon _0=8.85418710^{12}`$ As/Vm is the dielectric constant of the vacuum, and $`2\pi \mathrm{}=h`$ is Planck’s constant.

We wish to determine the ground-state wave function — from now on, $`\mathrm{\Psi }(𝐫_\mathrm{e},𝐫_\mathrm{p})`$ will refer to the ground state — together with the ground-state energy, and to calculate the following two ground-state related quantities. First, the probability

$$p_{\mathrm{out}}=\underset{r_pR_0}{}\text{d}^3𝐫_\mathrm{p}\underset{r_\mathrm{e}R_0}{}\text{d}^3𝐫_\mathrm{e}|\mathrm{\Psi }(𝐫_\mathrm{e},𝐫_\mathrm{p})|^2$$
(3)
that the positron will be found outside the hole. This quantity is a global measure of the pick-off annihilation rate of o-Ps. Second, the electron–positron contact parameter

$$\kappa =\frac{\underset{r_\mathrm{e}R_0}{}\text{d}^3𝐫_\mathrm{e}|\mathrm{\Psi }(𝐫_\mathrm{e},𝐫_\mathrm{p}=𝐫_\mathrm{e})|^2}{\text{d}^3𝐫_\mathrm{e}|\mathrm{\Psi }_{\mathrm{free}\mathrm{PS}}(𝐫_\mathrm{e},𝐫_\mathrm{p}=𝐫_\mathrm{e})|^2},$$
(4)
where $`\mathrm{\Psi }_{\mathrm{free}\mathrm{PS}}`$ is the ground-state wave function of the free positronium, having the form

$$\frac{1}{\sqrt{8\pi a_\mathrm{B}^3}}\mathrm{exp}\left(\frac{\left|𝐫_\mathrm{p}𝐫_\mathrm{e}\right|}{2a_\mathrm{B}}\right)$$
(5)
after separating the centre-of-mass motion, with $`a_\mathrm{B}=4\pi ϵ_0\mathrm{}^2/\left(m_\mathrm{e}e^2\right)=0.529\times 10^{10}`$ m denoting the Bohr radius. This quantity measures the overlap of the elecron and the positron, and is equal to the ratio of the self-annihilation rate of confined Ps to that of free Ps.

Because of the spherical symmetry of the system, the ground-state wave function depends on three independent scalar variables only. For these we choose the distance of the electron from the hole centre, $`r_\mathrm{e}=|𝐫_\mathrm{e}|`$, the electron–positron separation $`R=|𝐫_\mathrm{p}𝐫_\mathrm{e}|`$, and the angle $`\chi `$ between $`𝐫_\mathrm{e}`$ and $`𝐫_\mathrm{p}𝐫_\mathrm{e}`$, satisfying $`\mathrm{cos}\chi =𝐫_\mathrm{e}(𝐫_\mathrm{p}𝐫_\mathrm{e})/(r_\mathrm{e}R)`$.

A lengthy yet straightforward calculation shows that, in terms of these variables, the Hamiltonian reads

$``$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m_\mathrm{e}}}({\displaystyle \frac{^2}{r_\mathrm{e}^2}}+2{\displaystyle \frac{^2}{R^2}}+{\displaystyle \frac{2r_\mathrm{e}^2+R^2+2r_\mathrm{e}R\mathrm{cos}\chi }{r_\mathrm{e}^2R^2}}{\displaystyle \frac{^2}{\chi ^2}}2\mathrm{cos}\chi {\displaystyle \frac{^2}{r_\mathrm{e}R}}+2{\displaystyle \frac{\mathrm{sin}\chi }{R}}{\displaystyle \frac{^2}{r_\mathrm{e}\chi }}`$ (6)
$`+2{\displaystyle \frac{\mathrm{sin}\chi }{r_\mathrm{e}}}{\displaystyle \frac{^2}{R\chi }}+{\displaystyle \frac{2}{r_\mathrm{e}}}{\displaystyle \frac{}{r_\mathrm{e}}}+{\displaystyle \frac{4}{R}}{\displaystyle \frac{}{R}}+{\displaystyle \frac{2r_\mathrm{e}^2\mathrm{cos}\chi +R^2\mathrm{cos}\chi +2r_\mathrm{e}R}{r_\mathrm{e}^2R^2\mathrm{sin}\chi }}{\displaystyle \frac{}{\chi }}){\displaystyle \frac{e^2}{4\pi \epsilon _0}}{\displaystyle \frac{1}{R}}`$
The Hamiltonian (6) is not separable and does not contain a small parameter on which a perturbation expansion might be based, therefore, we perform the subsequent calculations numerically.

## III The method of computation

We use a variational method to compute the approximate ground-state wave function, minimizing the energy $`\mathrm{\Psi },\mathrm{\Psi }`$ in the function space

$$\mathrm{\Psi }=\underset{n=1}{\overset{N}{}}C_n\mathrm{\Psi }_n,$$
(7)
where the functions $`\mathrm{\Psi }_n(n=1,\mathrm{},N)`$ form a set of appropriately chosen orthonormal base functions. The minimization is carried out by finding the eigenvector of the matrix of the elements $`_{mn}=\mathrm{\Psi }_m,\mathrm{\Psi }_n`$ with the lowest eigenvalue. The components of this eigenvector are then identified with the coefficients $`C_n`$ of the approximate ground state wave function (7). The base functions are chosen in such a way that
(i) the $`r_\mathrm{e}`$-dependence of the wave function is similar to the ground state of a single electron in an infinitely deep potential well,

$$\psi _\mathrm{e}(r_\mathrm{e})=\text{const. }\frac{\mathrm{sin}(\pi r_\mathrm{e}/R_0)}{r_\mathrm{e}}$$
(8)
(see, e.g., the work of Galindo and Pascual Galindo and Pascual (1990)),
(ii) the $`R`$-dependence resembles that of the free positronium (cf. Eq. 5), and
(iii) that the matrix elements can be computed relatively fast.
To this end, we introduce the functions

$`\phi _i^{(1)}(r_\mathrm{e})=`$ $`(r_\mathrm{e}^2R_0^2)^i`$ $`(i=1,\mathrm{},N_1),`$ (9)
$`\phi _j^{(2)}(R)=`$ $`R^{j1}\mathrm{e}^{R/\xi }`$ $`(j=1,\mathrm{},N_2),`$ (10)
$`\phi _k^{(3)}(\chi )=`$ $`\mathrm{cos}^{k1}\chi `$ $`(k=1,\mathrm{},N_3),`$ (11)
from which, by means of the Gram–Schmidt orthonormalization procedure Schmidt (1907), we obtain the orthonormalized functions $`\psi _i^{(1)}(r_\mathrm{e})`$, $`\psi _j^{(2)}(R)`$ and $`\psi _k^{(3)}(\chi )`$, respectively. The products

$$\mathrm{\Psi }_{ijk}(r_\mathrm{e},R,\chi )=\psi _i^{(1)}(r_\mathrm{e})\psi _j^{(2)}(R)\psi _k^{(3)}(\chi )$$
(12)
serve as the base functions $`\mathrm{\Psi }_n`$, i.e., the expansion (7) is implemented as

$$\mathrm{\Psi }=\underset{i=1}{\overset{N_1}{}}\underset{j=1}{\overset{N_2}{}}\underset{k=1}{\overset{N_3}{}}C_{ijk}\mathrm{\Psi }_{ijk}.$$
(13)
The variational method is improved by considering $`\xi `$ as an additional, nonlinear adjustable parameter, with respect to which the energy is also minimized (for a free positronium, $`\xi =2a_\mathrm{B}=1.058\times 10^{10}`$ m).

For the first functions $`\psi _1^{(1)}(r_\mathrm{e})`$, $`\psi _1^{(2)}(R)`$ and $`\psi _1^{(3)}(\chi )`$, the orthonormalization procedure means only a normalization. It is plausible to expect that the product

$$\mathrm{\Psi }_{111}=\text{const. }\phi _1^{(1)}(r_\mathrm{e})\phi _1^{(2)}(R)\phi _1^{(3)}(\chi )$$
(14)
proves to be the dominant basis function in the ground state (13).

The matrix components $`_{mn}(R_0,\xi )`$ were calculated analytically by the software Maple V map  and subsequently exported to the format of programming language C, and the numerical solution of the eigenproblem, as a function of $`R_0`$ and $`\xi `$, was performed by a C program using routines from Numerical Recipes Press et al. (1992). For a fixed $`R_0`$, the lowest eigenvalue was minimized as a function of $`\xi `$, and the eigenvector corresponding to this eigenvalue was determined. The lowest eigenvalue gave the approximate energy of the ground state, while the components of the eigenvector provided the coefficients $`C_{ijk}`$ of the approximate ground-state function (cf. Eq. 13). Finally, using the resulting $`C_{ijk}`$ values, $`p_{\mathrm{out}}`$ and $`\kappa `$ were calculated again by Maple V. The computations took approximately two months on a workstation.

## IV Results

It was found that a base-function set with $`N_1=N_2=N_3=2`$ sufficed for determining the energy, $`N_1=N_2=N_3=3`$ for $`p_{\mathrm{out}}`$, and $`N_1=N_2=N_3=4`$ for $`\kappa `$ — in the sense that any further increase of the number of base functions produced only a negligible change in the values of $`E`$, $`p_{\mathrm{out}}`$, or $`\kappa `$. The results presented below were calculated using the base-function set $`N_1=N_2=N_3=6`$. The calculations were performed for 33 different values of the hole size in the range $`a_\mathrm{B}R_098a_\mathrm{B}`$.

The resulting ground-state energy $`E`$ is shown in Fig. 1 as the function of the hole radius. For large $`R_0`$ values, $`E`$ approaches the ground-state energy $`E_{\mathrm{free}\mathrm{PS}}=6.80\text{eV}`$ of the free positronium. The dependence of the energy on $`R_0`$ is well described by the function

$$E=|E_{\mathrm{free}\mathrm{PS}}|\left[8.75\left(\frac{R_0}{a_\mathrm{B}}\right)^2+11.78\left(\frac{R_0}{a_\mathrm{B}}\right)^{2.53}1\right].$$
(15)
The error of the numerical parameters in this fitted function is less then 3.7% . For comparison, Fig. 1 displays the ground-state energy of a single electron confined to a hole, $`E_\mathrm{e}=\pi ^2\mathrm{}^2/(2m_\mathrm{e}R_0^2)`$ (see, e.g., the work of Galindo and Pascual Galindo and Pascual (1990)), and the ground-state energy of a hydrogen atom in a hole with its proton in the centre (reproduced from the work of Sommerfeld and Welker Sommerfeld and Welker (1938)) as well.

The determined ground state $`\mathrm{\Psi }`$ of the positronium in a pore is characterized by the coefficients $`C_{ijk}`$ and $`\xi `$, Fig. 2 presents the most dominant coefficients and $`\xi `$, as the function of the hole radius. The coefficients are displayed in their squared form, $`|C_{ijk}|^2`$, which give the physical weights of the corresponding base functions $`\mathrm{\Psi }_{ijk}`$ in the ground state $`\mathrm{\Psi }`$. All coefficients $`C_{ijk}`$ were found real, as was to be expected.

The calculated $`R_0`$-dependence of the probability $`p_{\mathrm{out}}`$ that the positron can be found outside the void is shown in Fig. 3. The data may be fitted by the function

$$p_{\mathrm{out}}=\frac{1+0.018\left(\frac{R_0}{a_\mathrm{B}}\right)^2}{1+0.34\left(\frac{R_0}{a_\mathrm{B}}\right)^2+0.004\left(\frac{R_0}{a_\mathrm{B}}\right)^5}$$
(16)
with an uncertainty of the coefficients of less than 4% .

Similarly, the results for the density parameter $`\kappa `$ (cf. Fig. 4) may be described by the function

$$\kappa =1+\frac{1}{0.41\frac{R_0}{a_\mathrm{B}}+0.024\left(\frac{R_0}{a_\mathrm{B}}\right)^{3.22}}$$
(17)
with uncertaities of the fit parameters of less then 2% .

## V Analysis of the results

The first observation to make is that, as expected, for large hole radii, the system tends to the free positronium. This can be seen on each the quantities $`E`$, $`\xi `$, $`p_{\mathrm{out}}`$ and $`\kappa `$. Also, as expected, it is $`\mathrm{\Psi }_{111}`$ among the base functions $`\mathrm{\Psi }_{ijk}`$ whose presence is the strongest in the determined wave function, and this dominance increases with the size of the hole to almost $`100\%`$.

To investigate the ground state more closely, let us study its dependence on the variables $`r_\mathrm{e}`$, $`R`$ and $`\chi `$. This can be done by calculating the weights of the base function components $`\psi _i^{(1)}(r_\mathrm{e})`$, $`\psi _j^{(2)}(R)`$ and $`\psi _k^{(3)}(\chi )`$ in the wave function. These weights are defined as

$$W_i^{(1)}=\underset{j=1}{\overset{N_2}{}}\underset{k=1}{\overset{N_3}{}}|C_{ijk}|^2$$
(18)
for a $`\psi _i^{(1)}(r_\mathrm{e})`$, and analogously for a $`\psi _j^{(2)}(R)`$ and a $`\psi _k^{(3)}(\chi )`$. The calculated weights are shown on Fig. 5, as functions of the hole radius.

Concerning the $`r_\mathrm{e}`$-dependence, we can observe that the weights $`W_i^{(1)}`$ are almost independent of $`R_0`$, and that they are very close to the weights of $`\psi _i^{(1)}(r_\mathrm{e})`$s in the ground state of a single electron in a hole \[cf. Eq. (8)\],

$$|c_1|^2=\frac{945}{\pi ^3}0.983,|c_2|^2=\frac{93555(10\pi ^2)^2}{\pi ^{10}}0.017$$
(19)
(the higher weights $`|c_3|^2,|c_4|^2,\mathrm{}`$ are negligible so it is enough to consider the first two ones). Therefore, the $`r_\mathrm{e}`$-dependent part of $`\mathrm{\Psi }`$ is close to the wave function of the confined electron, even for small values of $`R_0`$, and with a better and better accuracy as $`R_0`$ increases.

The $`R`$-dependence is more influenced by $`R_0`$. In a small void, the components higher than the first one contribute up to $`20\%`$ and this amount decreases to zero as the hole size is increased. For large $`R_0`$, only the first component remains, and in the meantime its parameter $`\xi `$ approaches the value of the free positronium. The $`\chi `$-dependence contains a $`10\%`$ amount of the higher base functions for small $`R_0`$, which decreases to zero for increasing hole radius. Since $`\psi _1^{(3)}(\chi )=\text{const. }`$, this means that, for larger pore sizes, the ground state becomes $`\chi `$-independent.

Turning to the ground-state energy, a simple interpretation of its found $`R_0`$-dependence can be that the presence of the finite hole increases the energy of the electron by the amount of $`E_\mathrm{e}=\pi ^2\mathrm{}^2/(2m_\mathrm{e}R_0^2)`$, similarly as it does with a single electron. Fig. 6 shows that this interpretation is fairly good. In Appendix A, we show that, with an only somewhat more involved physical argument, it is possible to give an even more precise approximation to the energy as the function of $`R_0`$.

The probability that the positron can be found outside the hole, $`p_{\mathrm{out}}`$, is one for holes much smaller than $`a_\mathrm{B}`$, since the characteristic space region needed for the positron is of the order of $`a_\mathrm{B}`$, like in the free positronium (see Appendix A as well on this subject). On the other side, for large hole sizes, $`p_{\mathrm{out}}`$ falls down to zero as $`R_0^3`$ (cf. Eq. 16). This latter can also be understood, an explanation for this is presented in Appendix B.

The density parameter $`\kappa `$ also behaves as expected. For $`R_0`$ large, it is 1 since the positronium is almost free. By gradually squeezing the hole, the allowed space for the electron becomes smaller, it can no more go so far away from the positron as in a free positronium. This results in some increasing in the overlap-measuring $`\kappa `$. We can guess what happens when $`R_0`$ becomes much smaller than the Bohr radius: The electron is practically confined at the origin, the wave function of the positron becomes similar to what it is like in the antihydrogen atom, and, therefore, the $`R`$-dependent part of $`\mathrm{\Psi }`$ becomes

$$\frac{1}{\sqrt{\pi a_\mathrm{B}^3}}\mathrm{exp}\left(\frac{R}{a_\mathrm{B}}\right).$$
(20)
The $`\kappa `$ value corresponding to this wave function is 8 \[cf. Eqs. (4) and (5)\] so the increase that can be observed in Fig. 4 is expected to slow down for $`R_0<a_\mathrm{B}`$ and to converge to the value 8 reached at $`R_0=0`$.

We mention that a similar behavior of $`\kappa `$ is expected if one replaces the infinite potential well by some other confining potential. Limiting the electron to a region is plausible to cause a kind of squeezing the two particles onto each other, irrespective of the concrete form of the confining mechanism.

## VI Discussion

In the framework of our model, the lifetime of an orthopositronium is determined by the pick-off annihilation rate, proportional to $`p_{\mathrm{out}}`$, and the intrinsic $`3\gamma `$ decay rate of the orthopositronium, which is considered to be its vacuum value multiplied by $`\kappa `$. The angular correlation of the two emitted photons of the pick-off process can be obtained by generalizing the Fourier approach of Mukherjee and coworkers Mukherjee et al. (1998) to the total momentum of the present two-particle system.

The infinite potential well, chosen in our calculations presented here, can be a good approximation for the confining mechanism for the electron in the pore for some materials and less good for others. The same variational method can be applied for other potential shapes, modifying only the electronic base functions, in the appropriate way. Naturally, finding and using the corresponding electronic base functions for various potential shapes may make the calculations hard to carry out in practice. In this respect, it is worth mentioning a feasible alternative approach for this purpose. It is known from quantum mechanics that the boundary condition of requiring the vanishing of the wave function $`\mathrm{\Psi }`$, Eq. 2, is not the only possibility for confining a particle to a finite region Fülöp and Tsutsui (2000). There exists actually a one-parameter family of quantum mechanically allowed boundary conditions, of the form $`\mathrm{\Psi }+L\mathrm{\Psi }/r_\mathrm{e}=0`$, that is, the vanishing of a combination of the wave function and its derivative, where $`L`$ is an arbitrary (real) length scale parameter. The well-known and most frequently used case, the vanishing of the wave function itself is simply the special case $`L=0`$. Instead of choosing different potential shapes it is simpler to choose only a different value of $`L`$ in the boundary condition. The corresponding electronic base functions can still be chosen as polynomials, fulfilling now the new boundary condition. One can expect that the precise shape of the confining potential is not so important in practice. Bearing in mind that these boundary conditions express nothing but that the quantum probability cannot flow outside from the holeFülöp and Tsutsui (2000), and that any confining potential shape wishes to ensure the same — although not completely strictly —, we can guess that all confining potentials will, in their effect, resemble one or another boundary condition (potential well) with an appropriate $`L`$. Therefore, the possibilities provided by the general potential shapes may be fairly well represented with the one-parameter family of sharp potential wells.

In this paper, we have calculated the ground state only. To give account of the observed temperature dependence of the pick-off decay rate, it would be useful to determine the excited states, too. Indeed, as demonstrated by Goworek and coworkers Goworek et al. (1998) and Gidley and coworkers Gidley et al. (1999) for some variants of the Tao-Eldrup model, this temperature dependence can be attributed to the fact that it is also possible for the positronium to be in one of its excited states, approximately with the thermal equilibrial probability $`e^{E/kT}`$. Explaining the effect in our framework, the excited states possess different $`p_{\mathrm{out}}`$ and, consequently, a different annihilation rate. Calculating the excited states is also possible in the variational method, after suitable modifications.

In the end we mention that it would be interesting to perform a calculation similar to the one presented here for systems when the two particles are of different mass, e.g., when the positron is replaced with a proton or a $`\mu ^+`$ particle, as an interpolation between the positron and the infinitely heavy limit discussed by Sommerfeld and Welker. Unfortunately, our results cannot be simply ’renormalized’ to be applicable to such situations. The reason for this is that, for different masses, the Hamiltonian has a more general form than Eq. (6) has, and these differences are not only some easily rescalable numerical factors but mean a more general and complicated structure of the Hamiltonian. However, the considerations of Section 4 and the Appendices may more easily be adjusted to the case of different masses, and in this way at least some approximate information could be obtained for those systems.

###### Acknowledgements.

We wish to thank Nikolay Djourelov for useful references.

## Appendix A Estimating the ground state energy

A closer, yet not involved interpretation of the found $`R_0`$-dependence of the energy of the positronium in a hole is motivated by the simple derivation of the ground-state energy of the free hydrogen atom. There, based on the uncertainty relation, one approximates the average momentum of the electron by $`p=\mathrm{}/r`$, where $`r`$ is the size of the average space region “run” by the electron. Thus the energy of the electron is

$$\frac{p^2}{2m_\mathrm{e}}\frac{e^2}{4\pi \epsilon _0}\frac{1}{r}=\frac{\mathrm{}^2}{2m_\mathrm{e}}\frac{1}{r^2}\frac{e^2}{4\pi \epsilon _0}\frac{1}{r}.$$
(21)
Minimizing this expression in $`r`$ one reaches just the correct ground-state energy, and the corresponding $`r`$ is also nothing else than the Bohr radius $`a_\mathrm{B}`$.

As a next step, let us consider the free positronium. There, the average distance between the electron and the positron is $`2a_\mathrm{B}`$ (cf. Eq. 5). Neglecting the motion of the center-of-mass, the average region run by the effective particle of reduced mass $`m_\mathrm{e}/2`$ belonging to the relative motion is then $`2a_\mathrm{B}`$. The energy is, therefore, considered as

$$\frac{\mathrm{}^2}{2(m_\mathrm{e}/2)}\frac{1}{(2a_\mathrm{B})^2}\frac{e^2}{4\pi \epsilon _0}\frac{1}{2a_\mathrm{B}},$$
(22)
which proves to be just the correct energy value $`E_{\mathrm{free}\mathrm{PS}}`$. If we return to the ’two particles’-picture, the kinetic energy in Eq. (22) is shared by the two particles so each possesses the half of it,

$$\frac{1}{2}\frac{\mathrm{}^2}{2(m_\mathrm{e}/2)}\frac{1}{(2a_\mathrm{B})^2}=\frac{\mathrm{}^2}{2m_\mathrm{e}}\frac{1}{(2a_\mathrm{B})^2}.$$
(23)
Using this in the reverse way, we find that the characteristic space region for both particles is $`2a_\mathrm{B}`$.

Now, when we confine the electron in a hole, its characteristic region will be influenced by the hole as well. To this end, let us determine the characteristic space size of a single electron in a hole with the same logic as before, from its known energy:

$$E_\mathrm{e}=\frac{\pi ^2\mathrm{}^2}{2m_\mathrm{e}R_0^2}=\frac{\mathrm{}^2}{2m_\mathrm{e}}\frac{1}{(R_0/\pi )^2}$$
(24)
“implies” that the corresponding characteristic distance is $`R_0/\pi `$. foo

When the electron is both confined to a hole and attracted by a positron, its allowed space region will be determined by the smaller of the two corresponding length scales — except for some narrow intermediate region when the two length scales are equal or similar. A simple as well as reasonable formula for estimating such a joint length scale from two ones is

$$\frac{1}{l_{\mathrm{joint}}^2}=\frac{1}{l_1^2}+\frac{1}{l_2^2},$$
(25)
which in our case reads

$$\frac{1}{r_\mathrm{e}^2}=\frac{1}{(2a_\mathrm{B})^2}+\frac{1}{(R_0/\pi )^2}.$$
(26)
Concerning the space region of the positron, we make the simplest choice to consider it unaltered with respect to the case of the free positronium, i.e., to be $`2a_\mathrm{B}`$. The average distance between the electron and the positron, $`r_{\mathrm{ep}}`$, also has to be estimated: For this we can use the formula

$$r=\sqrt{r_1^2+r_2^2}$$
(27)
as being a “half-way” between the minimal distance $`|r_1r_2|`$ and the maximal one $`r_1+r_2`$ :

$`|r_1r_2|`$ $`=`$ $`\sqrt{r_1^22r_1r_2+r_2^2}\sqrt{r_1^2+r_2^2}`$ (28)
$``$ $`\sqrt{r_1^2+2r_1r_2+r_2^2}=r_1+r_2`$
($`r_1,r_20`$). Putting all these together, we estimate the energy as

$`E^{(\mathrm{approx})}`$ $`=`$ $`K_\mathrm{e}+K_\mathrm{p}+V(r_{\mathrm{ep}})={\displaystyle \frac{\mathrm{}^2}{2m_\mathrm{e}}}{\displaystyle \frac{1}{r_\mathrm{e}^2}}+{\displaystyle \frac{\mathrm{}^2}{2m_\mathrm{e}}}{\displaystyle \frac{1}{(2a_\mathrm{B})^2}}`$ (29)
$`{\displaystyle \frac{e^2}{4\pi \epsilon _0}}{\displaystyle \frac{1}{\sqrt{(r_\mathrm{e})^2+(2a_\mathrm{B})^2}}}.`$
As can be seen on Fig. 6, this formula provides a pretty good approximation for the ground-state energy of the confined positronium.

We can apply the same style of approach for the case of the hydrogen atom in a hole as well. There, the joint length scale of the length scale $`a_\mathrm{B}`$ of the free hydrogen and the one of the hole, $`R_0/\pi `$, is given by

$$\frac{1}{r_\mathrm{H}^2}=\frac{1}{a_\mathrm{B}^2}+\frac{1}{(R_0/\pi )^2}.$$
(30)
In the potential energy this same length scale will appear, since the proton stands in the center of the hole. Consequently, our estimate for the energy of a hydrogen atom in a hole is

$$E_\mathrm{H}^{(\mathrm{approx})}=\frac{\mathrm{}^2}{2m_\mathrm{e}}\frac{1}{r_\mathrm{H}^2}\frac{e^2}{4\pi \epsilon _0}\frac{1}{r_\mathrm{H}}.$$
(31)
Fig. 6 shows that our reasoning gives a fairly good approximation for this system as well.

## Appendix B The probability $`p_{\mathrm{out}}`$ for large holes

An explanation for the found large-$`R_0`$ asymptotic behavior $`p_{\mathrm{out}}R_0^3`$ can be given by the following argument: Let us rewrite Eq. (3) as

$$p_{\mathrm{out}}=8\pi ^2\underset{0}{\overset{\mathrm{}}{}}R^2\text{d}R\underset{0}{\overset{R_0}{}}r_\mathrm{e}^2\text{d}r_\mathrm{e}\underset{\chi _1}{\overset{\chi _2}{}}\mathrm{sin}\chi \text{d}\chi |\mathrm{\Psi }(r_\mathrm{e},R,\chi )|^2,$$
(32)
where the integration boundaries $`\chi _10`$ and $`\chi _2\pi `$ for $`\chi `$ are determined from the condition $`r_\mathrm{p}=\sqrt{r_\mathrm{e}^2+R^2+2r_\mathrm{e}R\mathrm{cos}\chi }R_0`$. Expanding $`|\mathrm{\Psi }(r_\mathrm{e},R,\chi )|^2`$ corresponding to Eq. (13) and also expanding all $`\psi `$s in it in terms of $`\phi `$s, let us examine the integral of one term,

$$\phi _i^{(1)}(r_\mathrm{e})\phi _i^{}^{(1)}(r_\mathrm{e})\phi _j^{(2)}(R)\phi _j^{}^{(2)}(R)\phi _k^{(3)}(\chi )\phi _k^{}^{(3)}(\chi ).$$
(33)
Concerning the variable $`r_\mathrm{e}`$, the main contribution to the integral must come from $`r_\mathrm{e}R_0`$, within the range $`|r_\mathrm{e}R_0|\xi `$, since otherwise $`R`$ has to be around $`R_0`$ to ensure $`r_\mathrm{p}R_0`$, but that is exponentially suppressed by $`e^{2R/\xi }e^{2R_0/\xi }`$ in the integrand. Now, for $`r_\mathrm{e}R_0`$ and $`RR_0`$, the condition $`r_\mathrm{p}R_0`$ is fulfilled by $`\chi `$ in the region $`\mathrm{\hspace{0.17em}0}\chi \pi /2`$. The integral of $`\mathrm{sin}\chi \phi _k^{(3)}(\chi )\phi _k^{}^{(3)}(\chi )`$ over this region gives an $`R_0`$-independent numerical factor. Turning to the $`r_\mathrm{e}`$-integration, for $`r_\mathrm{e}R_0`$, the factor $`r_\mathrm{e}^2\phi _i^{(1)}(r_\mathrm{e})\phi _i^{}^{(1)}(r_\mathrm{e})`$ behaves as

$`r_\mathrm{e}^2(r_\mathrm{e}^2R_0^2)^{i+i^{}}`$ $`=`$ $`r_\mathrm{e}^2(r_\mathrm{e}+R_0)^{i+i^{}}(r_\mathrm{e}R_0)^{i+i^{}}`$ (34)
$``$ $`R_0^2(2R_0)^{i+i^{}}(r_\mathrm{e}R_0)^{i+i^{}}`$
so its integral within the range $`\mathrm{\hspace{0.17em}0}|r_\mathrm{e}R_0|\xi `$ is proportional to $`R_0^{i+i^{}+2}\xi ^{i+i^{}+1}`$. The normalizing factors standing before $`\phi _i^{(1)}(r_\mathrm{e})`$ and $`\phi _i^{}^{(1)}(r_\mathrm{e})`$ bring in a factor $`R_0^{(2i+2i^{}+3)}`$, resulting in a behavior $`R_0^{(i+i^{}+1)}`$. The last integration, with respect to $`R`$, is an independent one,

$$\underset{0}{\overset{\mathrm{}}{}}R^2\text{d}R\phi _j^{(2)}(R)\phi _j^{}^{(2)}(R),$$
(35)
with a result depending only on $`\xi `$ so it does not modify the $`R_0`$-asymptotics (note that, for large $`R_0`$, $`\xi `$ is independent of the hole size, it tends to the constant $`2a_\mathrm{B}`$). Hence, we can see that the strongest $`R_0`$-asymptotics among the terms in $`p_{\mathrm{out}}`$ is caused by the $`i=i^{}=1`$ term, and is found to be $`R_0^3`$.