Magic Numbers in Boson ⁴He Clusters: The Auger Evaporation Mechanism

Abstract

The absence of magic numbers in bosonic ⁴He clusters predicted by all theories since 1984 has been challenged by high-resolution matter-wave diffraction experiments. The observed magic numbers were explained in terms of enhanced growth rates of specific cluster sizes for which an additional excitation level calculated by diffusion Monte Carlo is stabilized. The present theoretical study provides an alternative explanation based on a simple independent particle model of the He clusters. Collisions between cluster atoms in excited states within the cluster lead to selective evaporation via an Auger process. The calculated magic numbers as well as the shape of the number distributions are in quite reasonable agreement with the experiments.

1. Introduction

Van der Waals interatomic bonds are the weakest in nature and, unlike covalent bonds, are isotropic, Thus the atomic clusters of heavy rare gases tend to grow in close-packed icosahedral shapes, with magic numbers determined by pure geometrical constraints [1]. Not so for ⁴He clusters: in view of their liquid and generally superfluid state, magic numbers would seem to be unlikely. Although some evidence of magic numbers in ⁴He clusters not corresponding to special geometries has been reported as early as 1983 [2], theoretical studies of cluster energies have indeed failed to provide any evidence for magic numbers. Several calculations based on different Monte Carlo methods [3,4,5,6,7] that the binding energy per atom $\epsilon (N)$ of small ⁴He clusters, and in general for boson clusters [8], is a smooth monotonic function of the atom number N. Nevertheless, high-resolution matter-wave diffraction experiments by Brühl et al. [7,9] have shown that the abundance ρ(N) of ⁴He clusters formed during a free jet supersonic expansion into vacuum from a cryogenic source is not monotonously decreasing for increasing N, as expected from the behavior of $\epsilon (N)$, but shows a series of maxima at certain magic numbers. The numbers N = 10, 14, 26, 44, ... at which maxima occur [9] do not correspond to any of the geometric magic numbers for the heavier rare gas clusters nor of the possible magic numbers of the ions formed in the detection process. It is however interesting that, in case of stronger polarization forces exerted by an ion on the surrounding He atoms, close-packed arrangements and therefore geometrical magic numbers are favored. Note that ionization occurs after the free jet clusters have reached thermodynamic equilibrium and cannot alter the number distribution. see, e.g., [10]. Rather they could be explained on the basis of a thermodynamic equilibrium model for the growth of the clusters [7]. Each time a new collective state with quantum numbers (n,l), obtained from a dedicated diffusion Monte Carlo calculation, is bound, a sudden jump in the ⁴He cluster partition function Z(N) occurs. As a consequence the equilibrium constant K(N) for the aggregation process ⁴He_N−1 + ⁴He ↔ ⁴He_N, which is given by the function Z(N)/Z(N − 1), exhibits sharp peaks at N = N_nl in good agreement with the experimental magic numbers. Thus, the enhanced growth of the magic number clusters reconciles the apparent disagreement with the predicted monotonous behavior of the binding energy per atom.

2. Theory

The present theory provides an alternative explanation by accounting for the role of the elementary processes of evaporation occurring in a supersonic cluster beam with a very high speed-ratio as in the matter wave experiments [9]. It is assumed that the final asymptotic cluster size distribution is a result of two sequential processes. First the clusters grow rapidly and, as a result of the energy of recombination, are internally hot; then, in a next step they cool down to their final temperature by evaporation of atoms [11,12,13].

Note that all these theories are based on thermodynamics and predict only a gradual dependence of the rate of evaporation on the cluster number size.

Since the densities and the relative velocities of the clusters and single atoms and consequently their collision probabilities in the final stages of the expansion are comparatively small, the kinetics are dominated by unimolecular evaporation of single atoms. As a result, the hot large clusters are transformed into smaller and cooler clusters until a final distribution ρ(N) stabilizes at some average temperature (0.37 K for ⁴He droplets [14]). The rate equation for the densities ρ(N) can be written as

$$\begin{array}{c}\frac{d\rho (N)}{dt}=-P(N) \rho (N)+P(N+1) \rho (N+1)- {\displaystyle {\sum }_{N\prime }\sigma (N,N\prime )} \rho (N) \rho (N\prime )\\ +{\displaystyle {\sum }_{N\prime <N}\sigma (N-N\prime ,N\prime )} \rho (N-N\prime ) \rho (N\prime ) \end{array}$$

(1)

where P(N) is the rate of evaporation of a cluster consisting of N atoms and $\sigma (N,N\prime )$ is the probability of coalescence of two clusters of size N and N′. The fission of large clusters to smaller ones is neglected in Equation (1). The asymptotic steady state conditions far downstream from the nozzle are obtained by setting $d\rho (N)/dt=0$. Moreover, for clusters moving all at about the same speed the relative velocities corresponding to beam temperatures in the mK range [15] are so low that coalescence processes may be neglected. Under these conditions ρ(N) is, in a first approximation, inversely proportional to the ratio of the evaporation P(N).

Next it is assumed that an Auger process, involving two-body collisions, governs the evaporation of single atoms. As argued in the discussion below, phonon-induced evaporation is under the present physical conditions much less probable and is not considered. In the Auger process one bound He atom gains enough energy to leave the cluster in a collision with another atom, which drops from an initial excited state to a lower energy state (Figure 1). In order to investigate the basic physics of this process in boson clusters, a simplified independent particle model is considered in which a small number (N < 125) of ⁴He atoms is trapped in a spherical box potential of radius R = R(N) and depth $-V_0\cong -\epsilon (N)$. The model shows that the Auger evaporation rate strongly depends on the presence of bound one-particle states near the vacuum threshold, namely around those atom numbers at which a bunch of new bound states is created which are depleted compared to the adjacent cluster sizes. The role of bound states in determining the size distributions is similar to the analysis by Guardiola et al. [7] and suggests that Auger evaporation is actually the microscopic mechanism leading to the magic numbers.

In order calculate the rate of Auger collisions the wave functions of the atoms inside the clusters are required and for this reason an independent particle description is invoked. Although most current models of helium clusters are based on the diffusion Monte Carlo method (T = 0) and Path Integral Monte Carlo methods (T ≠ 0), or related algorithms, these models suffer from providing little information on the behavior of the individual constituents. The independent particle method, first introduced to ⁴He clusters and implemented via a variational Monte Carlo method by Schmid et al. in 1965 [16] and further refined by Lewart et al. [17] and Ramakrishna and Whaley [18], has been revived in connection with theories of Bose-Einstein condensed ultra-cold gases [19]. For calculating the collision dynamics between two bosons an exchange-symmetric two-atom wavefunction is required

$${\psi }_{ij}({\mathbf{r}}_1,{\mathbf{r}}_2)=\frac{1}{\sqrt{2(1+{\delta }_{ij})}}[{\varphi }_i({\mathbf{r}}_1){\varphi }_j({\mathbf{r}}_2)+{\varphi }_j({\mathbf{r}}_1){\varphi }_i({\mathbf{r}}_2)]$$

(2)

where ${\delta }_{ij}$ is the Kronecker symbol. Since the single particles move freely inside the cluster they are approximated by a plane-wave representation, ${\phi }_j(r)\propto \exp (i{k}_j\cdot r)$. Then the two-atom wavefunction can be expressed as a function of the internal coordinate $r=r_2-r_1$ and center-of-mass coordinate $r_{CM}=(r_2+r_1)/2$:

$$\begin{array}{ll}{\psi }_{ij}({\mathbf{r}}_1,{\mathbf{r}}_2)=& \frac{1}{\mathcal{N}}\exp (i {\mathbf{k}}_{CM}\cdot {\mathbf{r}}_{CM}) \cos (\mathbf{k}\cdot \mathbf{r})\\ & =\frac{\sqrt{4\pi }}{\mathcal{N}}\exp (i{\mathbf{k}}_{CM}\cdot {\mathbf{r}}_{CM}){\displaystyle \sum _{n=0}^{\infty }{(-1)}^n{(4n+1)}^{1/2}{j}_{2n}(kr) Y_{2n,0}(\theta )}\end{array}$$

(3)

where $k=k_j-k_i$, $k_{CM}=(k_j+k_i)/2$, m* is the atomic effective mass, $N$ a normalization constant, θ is the angle of r with respect to the conventional z axis. The atom-pair angular quantum number $L\equiv 2n$, labelling the m = 0 spherical harmonics $Y_{L,0}(\theta )$ and the spherical Bessel functions $j_L(kr)$ in Equation (3), is an even integer originating from the composition of the individual atom angular momenta l₁ and l₂. If the eigenfunctions of Equation (2) are expressed via free-particle spherical waves, ${\phi }_{klm}(\mathbf{r})\propto j_l(kr)Y_{l m}(\theta , \phi )$ and exchange is neglected, the matrix element $\left.〈{\psi }_{ij}\right|{\mathbf{k}}_j\cdot {\mathbf{k}}_i\left|{\psi }_{ij}\right.〉=0$, which ensures that the expectation value of the atom-pair internal kinetic energy ${\hslash }^2 k^2/4m^{\ast }$ is just $(E_i+E_j)/2$.

The coupling leading to Auger processes is essentially given by the He-He interatomic potential, which is here chosen in the Tang-Toennies (TT) form [20,21] restricted to the first dipole-dipole dispersion term in the London expansion,

$$V(r)=\left[A-\frac{C_6}{r_{}^6}{\displaystyle \sum _{k=7}^{\infty }\frac{{(\beta r)}^k}{k!}}\right] e^{-\beta r}$$

(4)

with $r=| r_2-r_1|$ the distance between atoms at positions r₁ and r₂, A = 22.16 a.u., β = 2.388 a.u., and C₆ = 1.461 a.u. [22]. The transition rate encompassing the two processes of Figure 1 is expressed as

$$\hspace{1em}W(E_i,E_j\to E_{i\prime },E_{j\prime })=\frac{2\pi }{\hslash }{\left| \left.〈{\psi }_{i\prime j\prime }({\mathbf{r}}_1,{\mathbf{r}}_2)\right|V(r) \left|{\psi }_{ij}({\mathbf{r}}_1,{\mathbf{r}}_2)\right.〉 \right|}^{2}g(E_i+E_j-E_{j\prime })$$

(5)

where g(E) is the free-atom density of states. Due to the short range of $V(r)$ as compared to the cluster diameter, surface effects are neglected.

Since the total energy is conserved ($E_{i\prime }+E_{j\prime }=E_i+E_j$) and also L is conserved due to the spherical symmetry of $V(r)$, a diagonal matrix element occurs in Equation (5). Although the TT potential $V(r)$ has no divergence, the Born approximation adopted in the Fermi golden rule of Equation (5) may not be sufficient due to the size of the repulsive potential, and distorted waves decaying exponentially for r smaller than the classical turning-point distance d₀ should be used. This complication is avoided by restricting the integrations in Equation (5) to the range d₀ ≤ r ≤ R₀, with the minimum distance $d_0$ corresponding to the atomic co-volume $b=4\pi d_0^3/3$ = 40 ų (excluded volume approximation) and R₀~R(N) the cluster radius, and using wavefunctions constructed on the basis of free-particle spherical waves.

In the present approximation the confinement due to the finite cluster radius effects acts through the discretization of energy levels $E_j\equiv E(n_j,l_j)$ appropriate to the cluster size, with n_j and l_j the radial and angular quantum numbers of the j-th level, respectively. An acceptable approximation for the energy eigenvalues of a spherical box of radius R and potential bottom −V₀ is [23]

$$E(n,l)\cong -V_0+\frac{{\pi }^2{\hslash }^2}{8m^{\ast }{R}^2}{(2n+l-1)}^2, (n=1,2, 3\dots ; l=0, 1, 2,\dots )$$

(6)

It is important to remark that this equation becomes exact for $E(n,l)\to 0^{-}$(threshold conditions) [23]. Moreover, by assuming a liquid drop model for the clusters, so that $R=R(N)=2.22N^{1/3}$Ả, it is found that Equation (7) reproduces the threshold (vanishing eigenvalue) condition for (³He)_N clusters (N = 40, 70, 112, 168) as calculated with a density functional [6] and the same value of $c\equiv b{(8m^{\ast }{V}_0)}^{1/2}/\pi \hslash =1.76\pm 0.03$. Since the evaporation probability, as shown below, depends essentially on the threshold quantum numbers, i.e., on the number of bound states, and much less on their detailed energy distribution, the following calculations based on the spherical-box model may be compared with the experimental data. For (⁴He)_N clusters no similar check can be made with density-functional results. However, Equation (6) provides a good fit of the single-particle levels obtained from the variational Monte Carlo calculations by Lewart et al. for (⁴He)₇₀ clusters [18] with c = 2.47 and a threshold number $2n+l-1$ between 10 and 11. Moreover a convenient expression for R(N) can be fitted to Pandharipande et al. Green’s function Monte Carlo results [3]:

$$R(N)=b_1{N}^{1/3}+\frac{b_2}{N^{1/3}-1}-\frac{3\pi \hslash }{{(8m^{\ast }{V}_0)}^{1/2}}$$

(7)

with b₁ = 1.88Å and b₂ = 1.13Å. With V₀ corresponding to the binding energy per atom for the bulk liquid (7.2 K [7]), an effective mass m* = 3.2 × 4 a.u. is obtained. Here only values of N above the minimum of R(N) at about N = 7 shall be considered. It should be noted that V₀ is actually a smooth function of N [7]. Assuming a constant V₀, however, is affecting the amplitude of magic number maxima but essentially not their position, which is the scope of the present analysis. The effects of a V₀(N) fitted to Guardiola et al. result [7] are discussed in Ref. [24].

The total Auger evaporation rate of one atom from an N-atom cluster is given by

$$P(N)={\displaystyle \int dE {\displaystyle \sum _{E_1,E_2,E_2\prime }\hspace{1em}W(E_1,E_2\to E, E_2\prime )}} n(E_1)n(E_2)[1+n(E_2\prime )]$$

(8)

where $n(E)={[e^{(E-{\mu }_B)/k_{B}T}-1]}^{-1}$ is the Bose-Einstein occupation number and the chemical potential ${\mu }_B\equiv {\mu }_B(T,N)$ is determined at each temperature from normalization to the atom number N. In the kinetic regime discussed above, controlled by single-atom evaporation, $\rho (N)\propto 1/P(N)$.

It appears from various numerical tests that the integrated evaporation rate, Equation (8), is rather insensitive to the choice of V₀ and m* as long as the number of one-particle bound states remains the same, thus confirming that the threshold behavior, i.e., the insertion of new states, rather than the bound state spectrum, determines the stepwise dependence on N of the evaporation rate. In general terms this is consistent with the thermodynamic interpretation by Guardiola et al. [7] based on a diffusion Monte Carlo analysis of the N-dependent partition function, though their calculation refers to a larger temperature (1.7 K) and to collective excitations with smaller threshold numbers.

The calculations show that the distribution $\rho (N)\propto 1/P(N)$ remains within the same order of magnitude for N varying from 10 to 100, whereas the size distribution function ${\rho }_{expt}(N)$ of the incident cluster beams, reported by Brühl et al. for different source pressures P₀, decreases for increasing N over a few orders of magnitude [9]. The latter were obtained from the diffraction experiments via multiplication by a Jacobian factor and a factor correcting for the cluster ionization probability, altogether giving a factor $\propto N^{-3}$. To isolate and identify the magic-number oscillations, the actual ${\rho }_{expt}(N)$ is divided by its smoothed version, ${\rho }_{fit}(N)$, obtained from averaging over the oscillations, actually a power law, ${\rho }_{fit}(N)$ ∝ $N^{-\alpha }$ with α ranging from 1.74, for a cluster beam source pressure P₀ = 1.33 bar, to 3.43 for P₀ = 1.10 bar. The dependence of the exponent α on P₀ clearly concerns the cluster formation kinetics close to the source, whereas the experiment indicates that the magic number oscillations are practically independent of P₀ and are therefore linked with the kinetics inside the travelling beam.

In Figure 3 the experimental distributions ${\rho }_{expt}(N)/{\rho }_{fit}(N)$ as reported by Brühl et al. [9] are compared with the corresponding theoretical distributions for two different temperatures and α = 3. For a closer comparison at large N the calculated ρ (N) has been replaced by the Gaussian convolution

$${\rho }^{\ast }(N)=\frac{1}{N^{2}s\sqrt{\pi }}{\displaystyle {\sum }_{N\prime }\rho (N\prime )\exp \left[-\frac{{(N\prime -N)}^2}{s^2{N}^4}\right]}$$

(9)

in order to reproduce the finite instrumental angular resolution in the cluster number $\mathrm{\Delta }N=s{N}^2$, where s ≅ 0.002.

The calculated distributions, plotted in Figure 2b for three different temperatures and N > 7, appear to reproduce very well all the salient features of the experiments. Some of the sharp peaks occurring in the calculated distributions are quite dependent on temperature: those at N = 11 and 21 are prominent at T = 0.37 K, but barely visible at 0.8 K, whereas the features for N > 40 increase with temperature. This can be interpreted as because clusters get smaller via more evaporation processes and are therefore favored at smaller temperatures. Experiments show a gradual switch of intensity from the peak at N = 21 to that at N = 26 for a decreasing source pressure P₀, whereas theory shows a similar effect for increasing cluster temperature.

The experimental oscillations in cluster abundance at larger N are better seen when plotted as a function of 1/N (Figure 3a). The calculated distribution for 0.8 K shows four peaks of comparable intensities: those at smaller N correspond to the experimental peaks for smaller source pressures, whereas those at larger N correspond to the experimental peaks at larger source pressures.

3. Conclusions

The good qualitative agreement between experiment and the calculated ⁴He cluster size distribution, although based on a simple free-particle model, provides a convincing argument in favor of the Auger evaporation mechanism for magic numbers (actually stability regions) of boson clusters. It is important to note that the Auger evaporation approach is conceptually similar to Guardiola et al. [6,7] quantum Monte Carlo thermo-dynamic approach in terms of threshold states, but they are alternative in terms of who is doing the job: single particle (actually pair collisions) or collective excitations?

Another issue concerning Brühl et al. experiments [9] is cluster formation kinetics. The observed similarity of peak distributions obtained at different source pressures after dividing out the respective smoothed distribution speaks in favor of in-beam kinetics. According to the present analysis Auger evaporation appears to be the basic kinetic mechanism for magic numbers in ⁴He clusters as generated in supersonic beams.