Lithium ions solvated in helium.

We report on a combined experimental and theoretical study of Li+ ions solvated by up to 50 He atoms. The experiments show clear enhanced abundances associated with HenLi+ clusters where n = 2, 6, 8, and 14. We find that classical methods, e.g. basin-hopping (BH), give results that qualitatively agree with quantum mechanical methods such as path integral Monte Carlo, diffusion Monte Carlo and quantum free energy, regarding both energies and the solvation structures that are formed. The theory identifies particularly stable structures for n = 4, 6 and 8 which line up with some of the most abundant features in the experiments.

Introduction

Ultra-cold helium nanodroplets have on numerous occasions proven to be a powerful tool for studying systems ranging from individual atoms to complex nanoparticles.1,2 Despite their low temperature (0.37 K) and weak interactions, they are remarkably good at solvating dopants, which typically reside near the core of the droplets and can, for example, be utilized for assembling clusters of one or more species. An interesting exception to this are the alkali metals that, at least in the case of individual atoms and small clusters, are “heliophobic” and reside in dimples on the surface of the droplets.3–5 This is due to a balance between the Pauli repulsion between the unpaired valence electron in the metal atoms and the closed shell He electrons on the one hand, and the surface tension of the droplet on the other. In contrast to neutral alkali metal atoms, alkali cations interact very strongly with He and are readily solvated, forming several layers of solvation shells around the ions. For the innermost layers, the He density can surpass the density of solid He, giving structures known as Atkins snowballs.6

Several theoretical and experimental studies have focused on the solvation structures of He snowballs around alkali metal cations. Reatto et al. employed variational Monte Carlo (MC) simulations with shadow wave functions to probe alkali ion impurities (Li+, Na+, K+, and Cs+) in liquid helium for equilibrium densities at 0 K.7–9 The chemical potential, local order, single particle excitation, and effective ionic mass were determined from the simulation model. A substantial difference in the snowball for corresponding ions could be observed. It was predicted that only Na+ and K+ have a tendency to form a solid snowball whereas the localization is not as prominent for Li+ and Cs+ species. Gianturco et al. later conducted a dedicated study based on the solvation of Li and other alkali metals in helium matrices employing a combination of classical energy minimization techniques and of exact quantum Diffusion Monte Carlo (DMC) methods.10–13 Small HeLi+ clusters with n ≤ 30 were considered for their investigations and they treated the full cluster interaction as a sum of pairwise potentials for Li+–He and He–He. It could be deduced that three particularly stable structures exist at n = 6, 8, and 10 with the most stable structure being found for n = 6. Additionally, evaluation of single particle evaporation energies, employing classical and quantum techniques, shed light on the formation of a rigid layer of helium with approximately 8 atoms being more tightly bound to the central ion. After this first shell, the evaporation energy was mainly governed by He–He interaction and not by interactions with the ionic core. The behavior was found to be similar to Na+ and K+ doped helium clusters, where the initial rigid layer was comprised of 9 and 12 He atoms, respectively. This rigid behavior of a fully developed first solvation shell for n = 8 was also reported in the ground state path integral calculation performed by Paolini et al.14 who found a stable structure of He atoms forming two parallel squares rotated by π/4 with respect to each other repeated in successively larger clusters (n ∼ 70, for example).

Previous investigations found that three-body (3B) contributions are rather insignificant in the stability of these helium clusters doped with alkali ions. Marinetti et al.13 observed some shifts of the radial distributions to slightly larger distances in their study on HeLi+ when the coupling between induced dipoles on the He atoms were taken into account. The ab initio calculations of the potential energy curve of HeLi+ and optimal structures for HeLi+ with n = 1–6 performed by Sebastianelli et al.12 concluded that the overall interactions were in fact governed mainly by diatom-like interactions between the ion and He atoms and that the pairwise approximation turns to be an acceptable description for these systems. The theoretical investigations performed by Issaoui et al.15 to study HeNa+ clusters added a self consistent many-body contribution between induced dipoles to the pairwise diatomic energy curves. The analysis performed on these studies revealed a notable overestimation of the energies predicted by the 2B approximation in larger clusters, which was found to delay the onset of delocalization and snowball features.15

Müller et al. carried out a systematic investigation on the formation and stability of helium snowballs created by employing femtosecond photoionization (PI) and electron impact ionization (EII) of alkali clusters (Na, K, Rb and Cs).16 From PI spectra, it could be deduced that alkali metal ions that originate from fragmented alkali clusters are more likely to constitute snowball complexes than their ionized monomer counterparts. This could be attributed to the fragmentation of clusters into singly charged ions, due to multiple ionization. Additionally, it was concluded that the size of a snowball with respect to the mass of alkali metals is a function of the kinematics of photofragmentation. For Na+ and K+ ions, they only observed the formation of small snowball sizes (up to 3 and 10 He atoms, respectively), which prohibited a direct comparison with predicted first shell closures from theory. However, with the heavier Rb+ and Cs+ ions they observed the formation of snowballs with up to 41 He atoms and identified the closures of the first solvation layers at He14Rb+ and He16Cs+, somewhat smaller than the shell sizes predicted by theory.9,16 Later, An der Lan et al. studied He droplets containing Na and K monomer and dimer cations.17 They reported on snowballs containing up 30 He atoms, with the first shell closures identified after He9Na+ and He12K+. The lightest alkali ion, Li+, was excluded from both of these experimental studies (and others like them) as the small mass and isotopic composition of the Li ion could potentially obstruct the evaluation of mass spectrometric data, corresponding to alkali-helium snowball complexes.

In this present work we report on the solvation of Li+ ions in helium, evaluated with high-resolution mass spectrometry measurements and different theoretical methods. In our experiments, HeLi+ complexes containing several tens of He atoms are identified and anomalies in specific cluster size yields let us probe the ion stabilities of these systems. These results are compared with both classical and quantum mechanical (QM) simulations of a Li+ ion solvated with He atoms. In particular, and in a similar fashion as previous investigations of clusters formed doping coronene molecules with rare gas atoms and H2,18–20 we have carried out basin-hopping (BH), DMC and path integral Monte Carlo (PIMC) calculations. In addition to this, estimations of the quantum free energy (QFE) have been calculated, leading to very similar results to those obtained with QM corrections of the BH results including zero-point energy (ZPE) effects. Geometries and energies of the stable configurations observed for the different HeLi+ clusters have been investigated and, in particular, the behavior as a function of the size of each cluster has been analyzed in an attempt to understand the abundances observed in the experiment for each n.

The structure of the paper is as follows: in Section 2 we present the essential details of the experimental setup, in Section 3 we present the theoretical approaches employed in this work, and in Section 4 results are shown and discussed. Finally in Section 5 the conclusions are listed.

Experimental details

Pure, pre-cooled helium (purity 99.9999%) with a stagnation pressure of 25 bar, was expanded through a 5 μm nozzle cooled to 6.5 K, leading to the formation of droplets with a broad size distribution (averaging about 107–108 He atoms). The resulting supersonic beam passed through a conical skimmer (diameter 0.8 mm) which is located 8 mm downstream from the nozzle. The skimmed beam travels across a 20 cm long differentially pumped pick up region where it is doped with high purity lithium (99% trace metals basis from Sigma Aldrich). The lithium sample was introduced into a cylindrical pickup cell under an inert atmosphere and covered with hexane to prevent oxidation during the transfer to the vacuum chamber. After vaporization of the hexane by evacuation at room temperature, the pickup-cell was resistively heated to a temperature of 750 K. The doped helium droplets were ionized by electron impact with kinetic energies of 70 eV. Resulting cations were then extracted from the ion source and guided into the extraction region of a commercial reflectron time-of-flight mass spectrometer (Tofwerk AG, model HTOF, mass resolution Δm/m = 1 : 5000). Additional experimental details have been described elsewhere.21,22

Theoretical methods

Potential energy surface

The employed force field is based on the sum of two-body (2B) He–Li+ and He–He non-covalent interaction contributions. For the He–He interaction we use the potential reported in ref. 23 while for the He–Li+ contribution we have developed a new potential based on accurate CCSD(T) results obtained in the complete basis set limit. In both cases the adopted analytical representation exploits the improved Lennard Jones (ILJ) formulation given by:24where ε is the potential depth, rm the position of the minimum and n(r) is defined as follows:24

The corresponding parameters for both the He–He and He–Li+ potentials using the ILJ analytical expression are given in Table 1. The effects of 3B terms are investigated by introducing an induced dipole-induced dipole interaction as that employed in previous studies13,25 with damping functions in our PES:where α = 1.31a03 for the He polarizability, r and r are He–Li+ distances, r is the He–He distance, and g(r) = f(r)/rni, where f(r) are the damping functions expressed as:26with b being equal to 2.9a0–1 or 3.2a0–1 for the He–Li+ and He–He interaction, respectively. The 3B calculation has been performed with a value of β = 9 in the ILJ description of the He–He interaction.

Table 1

Parameters for the ILJ potentials for the He–Li+ and He–He interactions. rm is given in Å, ε in meV; m and β are dimensionless

	m	r m	ε	β
He–Li+	4	1.90	81.3	4.2
He–He	6	2.97	0.947	8

Basin-hopping

The BH27 is a stochastic method to obtain the global minima of a potential energy surface (PES). This technique transforms the surface into a collection of basins which are explored by hopping between the local minima. Both local and global minima are preserved under this transformation. A Metropolis criterion using the energies of the initial and final minima in each step at a fictitious temperature determines if the attempted steps are accepted or rejected. The algorithm is particularly efficient since downhill barriers between different basins are removed and trapping is usually avoided. Moreover, size steps are typically larger than those employed for thermal sampling in MC simulations. The calculation was performed with a constant fictitious temperature such that kBT = 1.5 meV. The BH approach has been successfully employed for a large series of molecular clusters,27–36 and in particular, it has been an extremely useful tool in recent investigations of coronene doped with rare gas atoms and molecular hydrogen.18–20 Quantum effects can be included by means of the ZPE in the harmonic approximation,18,20 a calculation which requires to construct a database of local minima close to the global minimum for each cluster size.

Quantum free energy

The QFE of a specific minimum α at a temperature T of a cluster with n He atoms is given by:where kB is the Boltzmann constant and Z(T) is the partition function of the minimum α at constant temperature T. Under a harmonic approximation, this partition function can be written as:37where 𝒪 is the order of the point group corresponding to the minimum α, ℏ is the Planck constant, β = (kBT)–1, and ωαi is the i-th vibrational frequency associated with the minimum with energy E. In the present calculations we consider T = 2 K.

There is an alternative version of this method in which the partition function given in eqn (6) is replaced by its classical expression. This classical free energy tends asymptotically to the BH results when the temperature is decreased to zero, whereas the QFE tends to the QM corrected BH + ZPE values.

Diffusion Monte Carlo

QM calculations were carried out by means of the DMC method. In this algorithm, the time-dependent Schrödinger equation is transformed into a diffusion equation after substitution of the real time t by the imaginary time τ = it. The ground state wave function can be then obtained from the longest lasting term (τ → ∞) in the solution of the diffusion equation. Details of the method can be found elsewhere.18,38,39 Ground state energies and probability densities were computed using a code developed by Sandler and Buch,40,41 assisted with the descendant weighting method. For a given cluster size, six simulations were typically performed, each of them involving nine generations for the descendant weighting procedure. About 12 000 replicas were propagated with time steps ranging from 40 to 80 a.u. and for around 6000 steps. The initial set of replicas consisted in Gaussian spatial distributions (widths between 0.3–0.4 Å) around the classical equilibrium cluster geometry. It was found that the calculations can run optimally if the initial distribution is obtained by scaling the equilibrium geometries by a factor of 1.1–1.5.

Path integral Monte Carlo

The PIMC method has been described in detail before42,43 and therefore here we will restrict ourselves to give the most relevant aspects for the present calculation. In essence this approach is based on the expression of the density matrix at a temperature T as the product of M density matrices at a higher-temperature T′ = T × M. The density matrix is therefore evaluated in a collection of quantities ℛ ≡ {rα1,…,rαN}, where α runs over the M quantum beads, containing the position vectors rαi of the N particles which form the cluster: the Li+ and the n He atoms. In particular, the energy of each cluster can be obtained by means of the thermodynamic estimator developed by Baker:44with λm = ℏ2/2m, m being the mass of either He or Li+ and τ = β/M. The first term in eqn (7) corresponds to the classical kinetic energy multiplied by M and the average of the energy due to the spring-like interaction assumed between consecutive beads in the same ring describing a specific particle and the potential energy V is performed over the MC steps. The PIMC calculation has been performed at 2 K using M = 200 beads, which are moved in groups of 10 following a staging method.45,46

Results and discussion

Experimental results

A cationic mass spectrum from helium nanodroplets doped with lithium is shown in Fig. 1. The raw mass spectrum (in black) is dominated by the pure He+ cluster series and by subtracting a background measurement we can largely isolate the products containing lithium (in blue). The background measurement is performed under identical experimental conditions as the main experiment, but with the lithium-containing oven operating at a lower temperature (<600 K). This lower temperature is insufficient for vaporizing Li, but is sufficient for vaporizing Na, which is the main pollutant in our lithium sample and forms He snowballs of similar masses. From the reduced spectrum (blue in Fig. 1) we deduce the abundances of HeLi+ snowballs by fitting each peak with a Gaussian profile, two examples with particularly large overlaps with other peaks are shown in the insets. Beyond the HeLi+ series we observe pure Li+ clusters for n = 2 and 3, but no larger clusters irrespective of the pickup conditions (even with an oven temperature as high as 1100 K).

Fig. 1

Mass spectrum from helium nanodroplets doped with lithium. The black spectrum is the raw spectrum as measured and the blue spectrum is what remains following background subtraction (mainly ions containing Li). The background spectrum was measured with identical experimental conditions except that the Li-containing oven was operated at a lower temperature. This lower temperature is sufficient to vaporize the main pollutants in the sample (e.g. Na), but not Li. Gaussian fits (orange dashed lines in insets) are performed for each HeLi+ peak.

The integrated counts from each HeLi+ complex are shown in Fig. 2. The standout features in this spectrum are the local maxima observed for n = 2, 6, 8, and 14, indicating that these are particularly stable systems (compared to their neighbors). At higher masses there are a few dips in the spectrum at n = 21, 24, and 27–28 on top of an underlying distribution that smoothly tapers off towards large cluster sizes.

Fig. 2

Measured abundances of HeLi+ as a function of number of He atoms, n. Maxima are observed for n = 2, 6, 8, and 14, while distinct minima over the underlying size distribution are seen for n = 21, 24, 27, and 28. The statistical errors of the integrated yields are smaller than the data markers.

Theoretical results

From a theoretical point of view, we have investigated the behavior of some quantities as a function of the number of He atoms n in the cluster, searching for features which may point out the stability of specific sizes of HeLi+ consistent with the measured abundances shown in Fig. 2. Thus, in Fig. 3 we show the evaporation energies, ΔE = E – E, obtained by means of those methods discussed in Section 3. Despite the expected quantitative differences, the classical approach we show here, the BH method (see Section 3.2) exhibits qualitatively the same trend as the corresponding QM counterparts employed, that is, QFE, DMC and PIMC (discussed in Sections 3.3, 3.4 and 3.5, respectively). Results obtained with the QM corrected version of the BH approach, the BH + ZPE method, are almost identical to the QFE energies and are not included in the figure. The main abrupt changes observed in the corresponding energy curves as a function of the size of the cluster occur at n = 4, 6 and 8. Jumps in the evaporation energies for these same sizes were also reported in previous investigations of the system10,11,14 and were interpreted as indications of the presence of stable structures for He6Li+ and He8Li+.11

Fig. 3

Evaporation energies obtained by means of the BH (green circles), QFE (full red circles), DMC (blue triangles) and PIMC (black squares) methods.

Second energy differences, defined as Δ2E = E + E – 2E, are also a useful magnitude to search for stable HeLi+ clusters at specific numbers of He atoms. The results obtained by means of the PIMC, BH + ZPE and QFE approaches are shown in Fig. 4. As expected the features observed in the curve of the evaporation energy of Fig. 3 also manifest as peaks when we plot these Δ2E differences. Thus, noticeable maxima are observed also at n = 4, 6, and 8, which, in view of the BH result also included in Fig. 4, seem to have their origin in the minima of the PES. The comparison between the second energy differences obtained by means of the BH and those calculated with the other methods reveals however noticeable discrepancies between classical and QM approaches. The classical result suggests a similar feature at n = 10 as well but the QM calculations do not entirely confirm this regard, in apparent agreement with the experiment. The integrated yield shown in Fig. 2 also exhibits maxima at n = 6 and 8, whereas for n = 4, only a suggested shoulder is seen.

Fig. 4

Second difference energies (see text for details) obtained with the PIMC (black), QFE (red), BH + ZPE (blue) and BH (green) theoretical methods.

Two of the most prominent maximum peaks observed in Fig. 4, those for He6Li+ and He8Li+, are also clearly visible in the measured abundances shown in Fig. 2. The equilibrium geometries associated to the minimum energy configurations have been investigated before in previous works: He6Li+ exhibits a symmetrical octahedral configuration with the He atoms coordinating the Li+ impurity, located at the center of the cluster.10,13 He8Li+, on the other hand, has a stable configuration formed by two parallel squares rotated by π/4 to each other surrounding the Li+ ion,14 which was found also as the inner core of larger clusters, thus suggesting that it corresponds to the geometry of the first solvation shell. This hypothesis was confirmed by the integration of the shell performed by Paolini et al.14 with ground state path integral calculations which yielded a value of 8.24 He atoms.

In this work we have performed DMC and PIMC calculations of the probability density functions corresponding to the HeLi+ clusters with n = 4, 6 and 8, those which correspond to special features in the curves as a function of the number of He atoms shown in Fig. 3 and 4. In the top panels of Fig. 5 we show the PIMC distributions obtained using a representation on the Eckart frames for specific snapshots of the quantum beads for each atom and their corresponding average represented as a cloud surrounding the expected location of both the He and Li+ atoms. The choice of a system satisfying Eckart conditions47 to guarantee an optimal separation between rotation and vibration, is made here only for pictorial purposes. Analogously, geometries obtained by averaging the positions of the DMC replicas (after rotation to a common body-fixed frame) are included in the figure. Both methods yield distributions which are not far from the equilibrium structures predicted by classical energy minimization algorithms.10 Thus, the QM approaches find a structure for He4Li+ which contains the ionic impurity caged inside a tetrahedron formed by the four He atoms and, for n = 6 and 8, the above mentioned octahedral and parallel squares structures found in previous investigations are reproduced here by means of the DMC and PIMC calculations. Although the probability density functions (not shown here) for the inter-particle distances, He–He and He–Li+, and the corresponding angles obtained with the DMC approach, certainly exhibit an intrinsic broadening, the maxima are only slightly deviated with respect to the stick values predicted in classical energy minimization studies.

Fig. 5

Probability density functions for He4Li+ (left), He6Li+ (middle) and He8Li+ (right) obtained by means of the PIMC (top) and DMC (bottom) methods. Green color in the PIMC results are for the average location of the He atoms and red for the Li+ atom. The average position of He (Li+) atoms in the replica of the DMC distributions are in white (purple). See text for details.

Our calculations also reveal the stability of the structure found for He8Li+. Thus, Fig. 6 shows BH and PIMC results for He10Li+ indicating that both the classical optimized geometry and the QM probability density function consists, in essence, of the core observed at n = 8 with the two extra He atoms located over the center of each parallel square. This result is consistent with previous findings for this particular cluster.10 This trend is maintained even for larger cluster sizes, and the analysis of radial and angular distributions reveals that the inner shell is quite similar to the structure seen for He8Li+ (see ESI†).

Fig. 6

BH optimized geometry (left) and PIMC probability density function (right) for He10Li+.

The comparison between these theoretical results and the experimental abundances reveals agreement for peaks at n = 6 and 8. However, the prominent maximum seen for n = 2 in the experimental data in Fig. 2 does not have a definitive direct explanation from the theory. This abundance anomaly from n = 2 is also observed in experiments with Na+ and K+ ions17 which suggests that this is a product of the ionization mechanism itself, which is not covered by the simulations. One possible explanation is that a He2* is formed by the initial electron impact which then through associative Penning ionization forms a He2Li+ complex. Furthermore, the high abundance of the He14Li+ complex in the experiments is not reproduced by calculations. This structure could be explained by the nesting of a parallel square structure like He8Li+ in an octahedron like He6Li+, or vice versa, similar to the nested solvation shells observed for the HeAr+ 48 and H– 49 clusters. However, a particularly stable structure with such a geometry is not observed in the present simulations. In fact, further DMC calculations were carried out starting with a geometry where a He6Li+ octahedron is nested inside a cube formed with eight He atoms (a higher energy classical local minimum) but, after the simulations, the cluster rearranged to a structure with a core formed by eight atoms. This final geometry was not particularly stable as compared with their closest neighbors n = 13 and 15.

In an attempt to test the effect of 3B terms on our present results we introduce a conveniently damped induced dipole-induced dipole interaction contributions as in ref. 25 and calculate the corresponding evaporation energies. The comparison of this magnitude as a function of the number of He atoms obtained by means of the present BH and DMC methods is shown in Fig. 7. Some differences are certainly observed between those energies calculated only with the 2B pairwise description and those with the 3B terms included, especially for the smallest clusters n < 8. Beyond that size evaporation energies are practically the same regardlessly the potential interaction employed. However the qualitative trend is the same for both the classical and the QM results in the figure. In addition, the second difference energies (not shown here) calculated with the 3B effects do not exhibit substantially different features in comparison with Fig. 4, and in particular, no new peaks are seen. This suggests that contributions from terms beyond a mere 2B description, being significant in terms of the absolute energies of the clusters, do not improve in essence the comparison between theoretical and experimental results shown in this work.

Fig. 7

Evaporation energy as a function of the number of He atoms n for HeLi+ clusters calculated by means of the BH (black squares) and DMC (red circles) approaches using merely a pairwise description (empty symbols) and including 3B effects (solid symbols). See text for further details.

Conclusions

We have studied the formation and structures of Atkins snowballs around Li+ ions using high resolution mass spectrometry and a number of different theoretical methods. The experiments show a series of particularly abundant HeLi+ complexes at n = 2, 6, 8, and 14, as well as some weaker features such as minima at n = 21, 24, 27, and 28.

The theoretical results show that classical approaches such as basin-hopping (BH) predicts qualitatively similar cluster properties as the quantum mechanical approaches of Quantum Free Energy (QFE), Diffusion Monte Carlo (DMC), and Path Integral Monte Carlo (PIMC). The simulations identify three particularly stable He shells surrounding the Li+ ions for n = 4, 6, 8, and a slightly weaker structure at 10. The sizes of n = 6 and n = 8 line up well with the experimental findings, suggesting that these features in the experimental mass spectrum are the results of these clusters stable octahedral and parallel square geometries, respectively. A larger magic structure observed for He14Li+ is observed in the experiments, but not in the simulations, could indicate the formation of multiple rigid and nested solvation shells.

Conflicts of interest

There are no conflicts to declare.

Click here for additional data file.