Does Vibrational Delocalization Stabilize Multiply-Charged Neon

pubs.acs.org/JPCL

Does Vibrational Delocalization Stabilize Multiply-Charged Neon Clusters? F. Calvo* LASIM, Universite de Lyon and CNRS UMR 5579, 43 Bd du 11 Novembre 1918, F69622 Villeurbanne Cedex, France

clusters is theoretically investigated by ABSTRACT The dissociation of Ne2þ n extending the liquid drop model to include vibrational delocalization and finite temperature effects as well as by using explicit path-integral molecular dynamics simulations with a polarizable many-body potential. Both methods show that quantum nuclear effects are not sufficient to explain the particularly low appearance size of nc = 288 reported experimentally, hereby suggesting that other mechanisms, possibly electronic, play an important role. SECTION Dynamics, Clusters, Excited States

kþ (q-k)þ that the cluster Xqþ n dissociates into Xm and Xn-m and aims at determining the energy barrier E†(n,q;m,k) associated with this reaction. The appearance size nc(q) is then the lowest number n for which E†(n,q;m,k) > 0 for all m,k. According to the conventional model, the energy E(m,k) of Xkþ m is the sum of the neutral (van der Waals, or vdW) contribution E(m,0) for the m-atom cluster and of the purely Coulombic repulsion Ec(m,k) of the cluster, the k charges being distributed homogeneously in the spherical volume of the cluster. The vdW term is usually expanded in powers of n1/3, the successive terms accounting for volume, surface, edge, and vertices contributions9

M

ultiple ionization of an isolated system can destabilize it through various processes ranging from fission into few, comparably large fragments up to Coulomb explosion if the repulsion between charges vastly exceeds chemical bonding. The so-called Rayleigh instability1 is the engine of the electrospray device and is manifested on macroscopic,2 mesoscopic,3 molecular,4 and nuclear5 scales. It has been specifically investigated in atomic and molecular clusters as different as carbon,6 metals7,8 and rare gases.9-11 One fundamental property of multiply charged clusters is their appearance size, or the minimum number nc(q) of constituents X needed to stabilize Xqþ n against spontaneous dissociation on time scales relevant in mass spectrometry. Recently, the appearance sizes of cationic neon clusters have been measured10 and found to be significantly lower than the predictions of the liquid drop model (LDM), yet known to be accurate for most other weakly bound species.9 With respect to argon or xenon, neon clusters exhibit a significant quantum behavior, which suffices to change the stable structures12,13 or the melting point12,14,15 because of zero-point effects. Experimentally, the special stabilities associated with magic numbers have been found to differ partially from those in heavier rare gases.16 Quantum nuclear effects, although mentioned by the authors of ref 10, have not been characterized yet for charged clusters. Our purpose in this Letter is precisely to address this issue, by revisiting the LDM and incorporating vibrational delocalization in the calculation of energies and diameters. Beyond this phenomenological approach, we have performed quantum molecular dynamics simulations of the dissociation in selected Ne2þ n clusters using an explicit polarizable potential with accurate pair interactions. These simulations provide further insight into the dissociation mechanisms, including the barrier at finite temperature and the distribution of fragments. We start along the tracks of Echt and coworkers,9 who developed the original LDM by assimilating the clusters as homogeneous dielectric media. In brief, the LDM assumes

r 2010 American Chemical Society

Eðm, 0Þ ¼ ab m þ as m2=3 þ ae m1=3 þ av

ð1Þ

The expression for the Coulomb energy of a uniformly dense sphere of radius, r, and dielectric constant, ε, is9 Ec ðm, kÞ ¼

ε - 1 k2 3kðk - 1Þ þ 2ε rðmÞ 5εrðmÞ

ð2Þ

where atomic units have been used. The radius r(m) derives from the size m through the simple relation r(m) = (3m/4π)1/3σ, σ being the van der Waals radius. The two above expressions can be corrected for vibrational delocalization, but before doing so, we consider the contribution of a finite temperature. Replacing the energy E(m,0) by a Helmholtz free energy F(m,0;T)and assuming an harmonic approximation to estimate the vibrational partition function, F reads17 Fðm, 0; TÞ = Eðm, 0Þ þ Ezpe ðm, 0Þ þ kB Tð3m - 6Þ lnf1 - exp½ - βpωðmÞg

ð3Þ

Received Date: July 15, 2010 Accepted Date: August 19, 2010 Published on Web Date: August 26, 2010

2637

DOI: 10.1021/jz100959q |J. Phys. Chem. Lett. 2010, 1, 2637–2641

pubs.acs.org/JPCL

where kB and p are the Boltzmann and reduced Planck constants, respectively, β = 1/kBT, and Ezpe is the zero-point vibrational correction. ω h(m) is the geometric average of the vibrational frequencies at size m, and the arithmetic average Æωæ(m) enters the expression of the zero-point energy correction as 1 ð4Þ Ezpe ðm, 0Þ ¼ ð3m - 6ÞpÆωæðmÞ 2

2þ Table 1. Appearance Size nc(q = 2) Obtained for Ne2þ n and Arn Clusters Using the Liquid Drop Model in Its Classical and Quantum Versions at Zero Temperature and Close to the Melting Point

argon

T=0 T ≈ Tmelt

Expansions in powers of m are also used for both ω h(m) and Æωæ(m) in the range 1 e m e 1000.17 In addition to energies, temperature and quantum delocalization alter the cluster radius, an effect that can be captured in an harmonic oscillator picture. Using the results of ref 18, the average cluster radius Æræ(m,T) can be estimated from the mean square radius Ærc2æ(m,T) = r2(m) þ δr2(m,T) as 1/3

Æræðm, TÞ = rðmÞ þ with18 δr2 ðm, TÞ ¼

δr2 ðm, TÞ 2rðmÞ

1 X p 2mμ R ωR th βpωR =2

ð5Þ

ð6Þ

The final component of the LDM is the barrier energy E†, obtained from the contact sphere model9 as E† ðn, q; m, kÞ ¼ Eðm, kÞ þ Eðn - m, q - kÞ - Eðn, qÞ kðq - kÞ rðmÞ þ rðn - mÞ

ð8Þ

in which the last term accounts for the recoil kinetic energy after fission, equivalent to the work needed to bring the two fragments from infinity into contact. The main ingredients of the model are the expansion coefficients of E(m,0), ω h(m), and Æωæ(m), which were obtained by least-squares fitting on the entire set of Lennard-Jones clusters up to size n = 1000.19 These parameters are given as Supporting Information. The calculations were performed for both neon and argon, at temperatures of 0 K and near the cluster melting point Tmelt = 12 (neon) or 40 K (argon), assuming classical or quantum statistics. The appearance sizes nc(q = 2) resulting from these calculations are summarized in Table 1. The results for the classical case at T = 0 match those reported in the literature,9,20 within small deviations originating from the different coefficients in the expansion of E(m,0). At finite temperature, the appearance size increases as a consequence of the decreasing free energy in smaller clusters.17 This is in agreement with the study by Nakamura,20


quantum

classical

quantum

93 104

94 94

632 703

656 626

who further found that geometrical shells have a minor influence. Quantum vibrational effects at T = 0 also lead to slightly higher appearance sizes, which is a manifestation of the greater relative destabilization of small clusters, in an analogous way to pure thermal effects.21 Finally, once both temperature and quantum effects are included, a minor compensation occurs in the quantum free energies, which yields appearance sizes that are close to the original LDM data. In all cases, including delocalization and temperature effects never produces significant changes in nc. A closer inspection at direct (thermodynamical) and indirect (cluster radius) contributions reveals that each of them accounts for at most a few percent variation, very far from the factor of three experimentally observed.10 One possible flaw in the LDM applied to rather large neon clusters is the assumption that the charge is uniformly distributed inside the cluster. Singly charged rare-gas clusters consist of a small, two-to-four-atom flexible core carrying the charge, solvated by the remaining neutral atoms bound by vdW and polarization forces.22,23 In a large, multiply charged cluster, the charge is expected to be split into as many covalently bound cores; the core themselves repel each other because of the Coulomb interaction. Unfortunately, a proper modeling at the atomistic level is challenging because it involves a quantum description of both charge delocalization and correlation, along lines similar to the empirical valencebond model for multiprotonated water developed by Voth and coworkers.24 For the purpose of simulating the dissociation at the level of full atomistic details with some quantum treatment of the nuclei, it is not unrealistic to assume that the charged cores are can be only monomers and that the total cluster Neqþ n described as two Neþ ions solvated by (n - 2) Ne atoms. In this approach, which was used for instance by Gay and Berne 25 the cluster is bound in their related study of Xe2þ n clusters, by pairwise vdW interactions V2(r), the Coulomb repulsion Vc(r) = 1/r between the ions, plus some many-body polarization energy Vpol({ri}) with a standard form and no shielding.26 The explicit form of V2 was taken as an extended LennardJones expression, which is quantum chemically accurate against coupled-cluster calculations.27 The full expression and parameters of the potential energy are given as Supporting Information. The dissociation of two selected Ne2þ n clusters was simulated using path-integral molecular dynamics (PIMD) trajectories, with n chosen below and above the experimentally measured appearance size, namely, n = 267 and 309. The first task consisted of locating the most stable structures of these clusters, particularly the positions of the two ions in the cluster. These stable structures were obtained by basin-hopping global

where μ is the mass of a neon atom and the sum runs over all (3m - 6) modes. To a good approximation, this sum can be replaced by (3m - 6)/(ω h /2), which was checked to be h th βpω accurate within 15% for clusters containing up to 1000 atoms. In the classical limit f 0, the above expression leads to a thermal expansion that is linear in temperature. However, zero-point effects have a finite expanding contribution already at T = 0 ð3m - 6Þp ð7Þ δr2 ðm, T ¼ 0Þ ¼ 2mωμ

þ

neon

classical

2638


pubs.acs.org/JPCL

optimization, starting from the icosahedral global minima of the corresponding LJ clusters.19 Direct (classical) MD simulations of the dissociation process were attempted by constant-energy trajectories initially thermalized at 12 K; however, only a few of these trajectories turned out to dissociate before 100 ps, and mostly for the smaller cluster. To overcome the time scale gap in the experiment, we have instead performed driven dissociation trajectories by pulling the two charges away from each other with a time-dependent harmonic potential, similar to pulling experiments on polymers.28 This procedure is convenient to extract the free energy along the dissociation pathway or the potential of mean force (PMF). We follow here the approach of Hummer and Szabo28 by biasing the potential energy V(R) of the coordinates R with an umbrella contribution Vbias as 1 ~ VðR; d, tÞ ¼ VðRÞ þ ks ½rq - dðtÞ2 ¼ VðRÞ þ Vbias ðR; d, tÞ 2 ð9Þ where rq is the distance between the two charges, d(t) = d0 þ λt is a linearly increasing distance, and ks is a spring constant. The PMF is obtained by28 F ðrq Þ ¼ F 0 - kB T lnÆδðrq0 - rq Þ expð - β½WðtÞ - Vbias ðrq0 ÞÞæ ð10Þ where F 0 is an additive constant and W(t) is the work accumulated at time t28 Z t Z t D V~ Dd 0 WðtÞ ¼ dt ¼ k λ ½dðt0 Þ - rq dt0 ð11Þ s 0 0 Dd Dt 0 Figure 1. Potential of mean force for the separation of the two 2þ charges in Ne2þ 267 (upper panel) and Ne309 (lower panel) clusters at T=12 K, as obtained from classical (P = 1) and quantum pathintegral (P=24) molecular dynamics simulations. The most stable structures of the cluster are depicted, with bonds drawn only between close-by neutral and ionic atoms.

In eq 10, the brackets indicate a canonical average under initial conditions. Quantum nuclear effects are included through the path-integral formalism,29 for which the above nonequilibrium method holds.30,31 In practice, we use a Trotter discretization number of P = 24 and solve the equations of motion in normal mode coordinates with all masses equal to μ. An initial sampling stage at fixed temperature of T=12 K and fixed reference distance d(t) = d0 is carried out using Nose-Hoover thermostats on each normal mode. The distance d0 corresponds to the equilibrium value in each global minimum structure, d0=16.5 and d0 = 17.25 Å for Ne2þ 267 and Ne2þ 309, respectively. These thermalized trajectories provide a canonical sampling of 200 initial conditions that are saved every picosecond. Each initial condition is then propagated with the time-dependent biasing potential, without thermostatting. The simulations are stopped when the distance d(t) of the bias reaches 37 Å after 100 ps. It should be emphasized here that the biasing potential directly acts on the two neon cations but only indirectly influences the neutral atoms through their interactions with the ions. In particular, this procedure may result in a variety of dissociation channels, as expected from the actual experiments. Other technical details of the simulations, such as the PIMD equations of motion, the propagation algorithms, and the values of the various parameters used are given as Supporting Information. 2þ The free-energy profiles obtained for Ne2þ 267 and Ne309 are represented in Figure 1 in both the classical (P = 1) and


quantum (P=24) descriptions. The additive constant F 0 in eq 10 was chosen such that the long distance behavior recovers the (dominant) Coulomb repulsion, F (r) ≈ 1/r. For the two cluster sizes, the PMF exhibits comparable variations with a fission barrier of 0.2 to 0.5 eV in the classical case and about half due to vibrational delocalization. The softer barriers obtained in the quantum case are consistent with the smearing out of thermodynamical properties reported in several cluster studies.12,14,15,31 However, this result disqualifies quantum nuclear effects as the main cause for the enhanced stability of small Ne2þ n clusters because dissociation is exponentially enhanced with decreasing barrier height. Yet, the significant increase in the barrier when the cluster enlarges from 267 to 309 atoms is compatible with the appearance size being close to this range. The barrier shifts further away from the undissociated free-energy minimum with increasing size in the quantum system, but an opposite situation is found in the classical case. This effect is difficult to quantify because the minimum itself is not accurately defined. Additionally, we find that the dissociation barriers are also broad and extend over 10 Å, thereby ruling out tunneling as a likely mechanism. The simulations still provide useful insight into the dissociation products, whose size is not constrained by the

2639


pubs.acs.org/JPCL

implicitly incorporate quantum vibrational effects, because they were fitted to reproduce experimental properties of the neon fluid. The more accurate potential V2(r) used in the molecular dynamics simulations reflects a proper electronic, Born-Oppenheimer ground-state energy curve, and the lowering of the dissociation barrier is at variance with the possible explanation for the higher stability of charged neon clusters as being due to quantum nuclear effects. The treatment of vibrational delocalization through the path-integral formalism is currently one of the only few that is practical for such large clusters, but it would be useful to extend the present nonequilibrium simulations to estimate proper dissociation rates, for instance using transition path sampling methods. At a deeper level, a quantum chemical approach to the electronic structure in Ne2þ n would be highly desirable to address the issue of charge delocalization on more realistic grounds. A fully quantum modeling of both electrons and nuclei may then, hopefully, shed more light onto the surprisingly high stability of multiply charged neon clusters.

SUPPORTING INFORMATION AVAILABLE All parameters used in the liquid drop expansions for the neutral energies E(m,0) and average vibrational frequencies ω h(m) and Æωæ(m). Path-integral MD equations of motion and numerical parameters used for the simulations and biasing potential. This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION Corresponding Author: *To whom correspondence should be addressed. E-mail: fcalvo@ lasim.univ-lyon1.fr.

REFERENCES (1)

Figure 2. Distribution of fragments obtained from (a) Ne2þ 267 and (b) Ne2þ 309 at the end of the dissociation process. The results of classical (upper panels) and quantum (lower panels) path-integral simulations are shown with the corresponding most abundant small cluster products.

(2)

(3)

computational procedure. We have represented in Figure 2 the distributions of fragment sizes obtained at the end of the pulling process for the two cluster sizes and in both classical and quantum modelings. In all cases, and in qualitative agreement with the LDM, fission is highly asymmetric. The most abundant small fragment contains between 12 and 22 atoms, and the distributions become more symmetric but also broader upon including quantum effects. The 13-atom cluster, which is a magic size for the present classical potential, is not particularly abundant as a dissociation product, therefore suggesting that fission is not entirely determined by energetics but that kinetic details should not be neglected. Even though quantum vibrational effects have a qualitative impact on the dissociation product, they are unsufficient to explain the much higher stability of doubly charged neon clusters, and this conclusion has no reason to be altered for higher charges. At this stage, it is useful to remember that the LJ parameters used in the original LDM analyses for neon9,10


(4) (5) (6)

(7)

(8)

(9)

2640

Lord Rayleigh, L. On the Equilibrium of Liquid Conducting Masses Charged with Electricity. Philos. Mag. 1882, 14, 184–186. Pruvost, L.; Serre, I.; Duong, H. T.; Jortner, J. Expansion and Cooling of a Bright Rubidium Three-Dimensional Optical Molasses. Phys. Rev. A 2000, 61, 053408-1–053408-9. Duft, D.; Lebius, H.; Huber, B. A.; Guet, C.; Leisner, T. Shape Oscillations and Stability of Charged Microdroplets. Phys. Rev. Lett. 2002, 89, 084503-1–084503-4. Echt, O.; Sattler, K.; Recknagel, E. Critical Sizes for Stable Multiply Charged CO2 Clusters. Phys. Lett. A 1982, 90, 185–189. Bohr, N.; Wheeler, J. A. The Mechanism of Nuclear Fission. Phys. Rev. 1939, 56, 426–450. Chabot, M.; Mezdari, F.; B eroff, K.; Martinet, G.; Hervieux, P.-A. Scaling Law for the Partitioning of Energy in Fragmenting Multicharged Carbon Clusters. Phys. Rev. Lett. 2010, 104, 043401-1–043401-4. Br echignac, C.; Cahuzac, Ph.; Carlier, F.; de Frutos, M.; Barnett, R. N.; Landman, U. Dynamics and Energy Release in Fission of Small Doubly Charged Clusters. Phys. Rev. Lett. 1994, 72, 1636–1639. Calvo, F. Role of Charge Localization on the Coulomb Fragmentation of Large Metal Clusters: A Model Study. Phys. Rev. A 2006, 74, 043202-1–043202-14. Echt, O.; Kreisle, D.; Recknagel, E.; Saenz, J. J.; Casero, R.; Soler, J. M. Dissociation Channels of Multiply Charged van der Waals clusters. Phys. Rev. A 1988, 38, 3236–3248.


pubs.acs.org/JPCL

(10)

(11)

(12)

(13)

(14)

(15)

(16) (17) (18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

M€ahr, I.; Zappa, F.; Denifl, S.; Kubala, D.; Echt, O.; M€ ark, T. D.; Scheier, P. Multiply Charged Neon Clusters: Failure of the Liquid Drop Model? Phys. Rev. Lett. 2007, 98, 023401-1– 023401-4. Hoener, M.; Bostedt, C.; Schorb, S.; Thomas, H.; Foucar, L.; Jagutzki, O.; Schmidt-B€ ocking, H.; D€ orner, R.; M€ oller, T. From Fission to Explosion: Momentum-Resolved Survey over the Rayleigh Instability Barrier. Phys. Rev. A 2008, 78, 021201(R)1–021201(R)-4. Calvo, F.; Doye, J. P. K.; Wales, D. J. Quantum Partition Functions from Classical Distributions: Application to RareGas Clusters. J. Chem. Phys. 2001, 114, 7312–7329. Deckman, J.; Frantsuzov, P. A.; Mandelshtam, V. A. Quantum Transitions in Lennard-Jones Clusters. Phys. Rev. E 2008, 77, 052102-1–052102-4. Neirotti, J. P.; Freeman, D. L.; Doll, J. D. A Heat Capacity Estimator for Fourier Path Integral Simulations. J. Chem. Phys. 2000, 112, 3990–3996. Predescu, C.; Sabo, D.; Doll, J. D.; Freeman, D. L. Heat Capacity Estimators for Random Series Path-Integral Methods by Finite-Difference Schemes. J. Chem. Phys. 2003, 119, 12119-1–12119-10. M€ark, T. D.; Scheier, P. Production and Stability of Neon Cluster Ions up to Neþ 90. Chem. Phys. Lett. 1987, 137, 245–249. Doye, J. P. K.; Calvo, F. Entropic Effects on the Structure of Lennard-Jones Clusters. J. Chem. Phys. 2002, 116, 8307–8317. Calvo, F.; Doye, J. P. K.; Wales, D. J. Equilibrium Properties of Clusters in the Harmonic Superposition Approximation. Chem. Phys. Lett. 2002, 366, 176–183. Wales, D. J.; Doye, J. P. K.; Dullweber, A.; Hodges, M. P.; Naumkin, F. Y.; Calvo, F.; Hern andez-Rojas, J.; Middleton, T. F. The Cambridge Cluster Database. http://www-wales.ch.cam. ac.uk/CCD.html. Nakamura, M. Shell Effects on Stability and Fragmentation of Multiply Charged Rare Gas Clusters. Chem. Phys. Lett. 2007, 449, 1–5. Deckman, J.; Mandelshtam, V. A. Quantum Disordering versus Melting in Lennard-Jones Clusters. Phys. Rev. E 2009, 79, 022101-1–022101-4. Naumkin, F. Y.; Wales, D. J. Structure and Properties of Neþ n Clusters from a Diatomics-in-Molecules Approach. Mol. Phys. 1998, 93, 633–648. Bonhommeau, D.; Viel, A.; Halberstadt, N. Dissociative Ionization of Neon Clusters Nen, n = 3 to 14: A Realistic Multisurface Dynamical Study. J. Chem. Phys. 2005, 123, 0543161–054316-13. Wang, F.; Voth, G. A. A Linear-Scaling Self-Consistent Generalization of the Multistate Empirical Valence Bond Method for Multiple Excess Protons in Aqueous Systems. J. Chem. Phys. 2005, 122, 144105-1–144105-9. Gay, J. G.; Berne, B. J. Computer Simulation of Coulomb Explosions in Doubly Charged Xe Microclusters. Phys. Rev. Lett. 1982, 49, 194–198. Stillinger, F. H. Dynamics and Ensemble Averages for the Polarization Models of Molecular Interactions. J. Chem. Phys. 1979, 71, 1647–1651. Schwerdtfeger, P.; Gaston, N.; Krawczyk, R. P.; Tonner, R. E.; Moyano, G. E. Extension of the Lennard-Jones Potential: Theoretical Investigations into Rare-Gas Clusters and Crystal Lattices of He, Ne, Ar, and Kr Using Many-Body Interaction Expansions. Phys. Rev. B 2006, 73, 064112-1–064112-19. Hummer, G.; Szabo, A. Free Energy Reconstruction from Nonequilibrium Single-Molecule Pulling Experiments. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 3658–3661.


(29) (30)

(31)

2641

Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; MacGraw-Hill: New York, 1965. van Zon, R.; Hern andez de la Pe~ na, L.; Peslherbe, G. H.; Schofield, J. Quantum Free-Energy Differences from Nonequilibrium Path Integrals. I. Methods and Numerical Application. Phys. Rev. E 2008, 78, 041103-1–041103-11. Hern andez de la Pe~ na, L.; Peslherbe, G. H. Quantum Effects on the Free Energy of Ionic Aqueous Clusters Evaluated by Nonequilibrium Computational Methods. J. Phys. Chem. B 2010, 114, 5404–5411.


Does Vibrational Delocalization Stabilize Multiply-Charged Neon

Recommend Documents