In collaboration with other research groups at the University of Notre Dame, one of our primary research topics is the structure and dynamics of nanoparticles containing coinage (Group IB) metals. We use molecular dynamics to study structural and dynamic properties of nanoparticles that would be difficult to obtain experimentally. All of our calculations were performed using our group simulation code, OOPSE [1]. Any image on this page can be enlarged by clicking on it.

Dynamics of Alloying in Bimetallic Nanoparticles

Using X-ray fine structure absorption spectroscopy (XAFS), the Bunker group in the department of Physics observed that gold-core / silver-shell nanoparticles in aqueous solution at ambient temperature form an alloy near the interface between the two metals.[2] The amount of inter-diffusion is highly dependent on the size of the nanoparticle core with larger cores alloying on a substantially slower time scale. This alloying occurs with a timescale of approximately three days (which is 7-9 orders of magnitude too fast to be explained by normal solid-state diffusion). Diffusion of atoms in a liquid droplet would cause the entire particle to alloy in a few microseconds. Our group along with the Meisel group of the Notre Dame Radiation Laboratory and the Bunker group, we have helped to explain this phenomenon. Our calculations show that the fast alloying times are due to vacancies present at the interface layer between the two metals. To gain insight into the alloying mechanism, molecular dynamics simulations were carried out on nanoparticles composed of Ag and Au atoms interacting under the Embedded Atom Method (EAM) potential. [3] Gold and silver have nearly identical lattice constants, so the nanoparticles were constructed on a perfect FCC lattice using a lattice constant of 4.085 Å. Gold atoms were placed inside of the core radius (rcore) and silver atoms were filled in out to the particle radius (rshell ). We designated a 2 Å "interface" region centered on rcore where either 5% or 10% of the interfacial atoms were removed at random for the simulations that involved vacancies. For these simulations, rcore = 1.25 nm and rshell = 1.998 nm matching the smallest of the nanoparticles studied by our collaborators. Before starting the molecular dynamics simulations, a relatively short steepest-descent minimization was performed to relax the lattice in the initial configuration. During the initial 30 ps of each trajectory, velocities were repeatedly sampled from a Maxwell-Boltzmann distribution matching the target temperature for the run. Following this initialization procedure, trajectories were run in the microcanonical (constant-NVE) ensemble with zero initial total angular momentum. In order to compare the structural features obtained from the NVE-ensemble molecular dynamics, additional trajectories were run using a modified Nosé-Hoover thermostat to maintain constant temperature (NVT) and zero total angular momentum. Data collection for all of the simulations started after the 30 ps equilibration period had been completed. We simulated particles with the above-mentioned interfacial vacancy density at 100 K intervals from 500 to 1200 K. Given the masses of the constituent atoms, we were able to use time steps of 5 fs while maintaining excellent energy conservation. Typical simulation times were 100 ns for nanoparticles simulated at 500-600 K and 12-24 ns for particles at 600-1200 K. From simulation results, we have constructed a radial density profile of the two constituents as a function of distance from the center of the nanoparticle. This density profile was obtained from the last 2 ns of the 24 ns simulations at 800 K. Remarkably, the interfacial vacancies result in a substantial smoothing of the peaks in the density profile, indicating that those nanoparticles are closer to the melting transition. The region near 1.25 nm in the radial density profile shows the interfacial vacancies result in significantly enhanced radial diffusion of the Ag into the core region. Much of the displacement occurs along the interface (i.e., at constant r) which can be seen from the broadening of the density peaks.

An Au(core)-Ag(shell) bimetallic nanoparticle with vacancy defects (red) at the interface

Radial Density Profile

Interdiffusion in Nanoparticles

Arrhenius Plot of Relaxation Rates

Atoms confined to a spherical volume of radius R have a mean square displacement the approaches 6R2/5 in the infinite time limit. (In bulk systems, the mean squared displacement is unbounded and can grow linearly with time.) At short times, a useful approach might be to connect the observable displacements of atoms in the simulation to solutions of the appropriate diffusion equation. For a spherically symmetric volume and reflecting boundary conditions at R, the solutions to the diffusion equation are given by where the coefficients an are determined by the initial displacements of the atoms. Proper mapping of the displacements of the individual atoms would require projecting initial and final positions of each atom onto the diffusional modes expressed in the summation. An alternative approach is to compute the relaxation time of the mean squared displacement using the correct infinite time limit and the Kohlraush-William-Watts (KWW) relaxation law, [4]

The mixing relaxation times are displayed in a standard Arrhenius plot to calculate the energy barrier associated with the diffusion mechanism. At the high end of the temperature scale the nanoparticles are liquid droplets, and since diffusion in the melt occurs via a different mechanism than in the solid, these temperatures have been excluded from the plot. Temperatures lower than ~500 K result in atomic displacements too small to provide a accurate estimate of the relaxation time. With no vacancies we obtain a activation energy of 1.85 eV (consistent with experimental determinations). [5] Nanoparticles with 5% and 10% vacancy levels have activation energies of 0.92 and 0.11 eV respectively. Projecting to 300 K, the room temperature alloying times are in excess of 10⁷ years with no vacancies, a few days for 5% vacancies and minutes to seconds for particles containing 10% vacancies. The alloying time associated with 5% interfacial vacancies is closest to the experimentally observed result.

Vibrational Dynamics in Gold Nanoparticles

Nanoparticles exhibit behavior not displayed in the bulk because of a high surface area to volume ratio. Energy due to photons absorbed during ultrafast laser excitation is transfered into thermal excitation of atomic degrees of freedom, coherently exiting the breathing mode of the nanoparticle. Dephasing of this breathing mode has been seen in laser experiments by Greg Hartland's group of the Department of Chemistry and Biochemsitry. [6] Dephasing following the coherent excitation of the breathing mode may be due to inhomogeneous size distributions in the sample, but may also be due to softening of the breathing mode vibrational frequency following a melting transition. We have performed molecular dynamics simulations that mimic the laser excitation experiments of Hartland's lab to gain insight into the possible mechanism for this breathing mode dephasing. [7] Spherical Au nanoparticles were created from the bulk in a standard FCC latice at four different radii [20 Å (1926 atoms), 20 Å (6602 atoms) and 35 Å (10606 atoms)] using a lattice spacing of 4.08 Å. Potentials for the Au interactions were once again provided by the EAM model. The ultrafast laser excitation event was modeled by sampling from a Boltzmann Distribution at twice the target temperature. Following the excitation event, nanoparticles evolved in a canonical ensemble using Nosé-Hoover NVT at the target temperature for 40 ps. Five independent samples with target temperature spanning a range from 300 K to 1350 K were studied. Volume was determined as a function of time by using the convex hull geometric construction which defines the smallest convex polyhedron which include all of the atoms as calculated by the Quickhull algorithm. [8] From these simulations, we find that there is both a size and temperature-dependent softening of these nanoparticles even at temperatures below the melting temperature. In a plot of the Bulk Modulus (K) for the various size nanoparticles studied as a function of final temperature after excitation, we can see a dramatic (and size-dependent) drop in K at temperature well below the melting temperature of bulk polycrystalline gold. (Polycrystalline gold has a bulk modulus of 220 GPa at a temperature of 300 K.) Surface melting occurs at even lower temperatures as indicated by our calculations of the radial-density profile which shows a merging of the crystalline peaks in the outer layer of the nanoparticle much like that seen in the Ag-Au nanoparticle with vacancies. The time-dependent estimate for the bulk modulus indicates a time scale for softening of about 10 ps independent of particle size. Correspondingly, the smaller particles exhibit coherent breathing vibrations for a few periods before melting and larger particles melt within the first vibrational period.

35 Å Gold Nanoparticles and their Convex Hulls Before and After Excitation

Bulk Modulus of Gold Nanoparticles

Glass Nanobeads: rapid cooling of silver-copper nanoparticles

Cutaway views of 30 Å Silver-Copper Nanoparticle Structures for Random (top) and Copper (core) - Silver (shell) (bottom) in Solid, Liquid, and Glass

In our study of the laser excitation of gold nanoparticles, we observed that the cooling rate for these particles (1011-10 12 K/s) which is in excess of the cooling rate required for glass formation in bulk metal alloys. We believe that it may be feasible to use laser excitation to melt, alloy and quench a metallic nanoparticle in order to form a metallic glass nanobead. We have chosen the bimetallic alloy of Silver (60%) and Copper (40%) as a model system because it is an experimentally known glass former and has been used previously as a theoretical model for glassy dynamics. [9] We have modeled the laser excitation and cooling by constructing and relaxing the eutectic composition (Ag6Cu4) on a FCC lattice with a lattice constant of 4.09 Å for 20, 30 and 40 Å radius nanoparticles. The nanoparticles are melted at 900 K and allowed to mix for 1 ns. Resulting structures are then quenched using a implicit solvent model where Langevin dynamics is applied to the outer 4 Å radius of the nanoparticle and normal Newtonian dynamics are applied to the rest of the atoms. By fitting to experimentally-determined cooling rates, we find that collision frequencies of 3.58 fs-1 for Ag and 5.00 fs-1 for Cu lead to nearly exact agreement with the Temperature vs. time data. The cooling rates are therefore 2.37 x 1013 K/s, 1.37 x 1013 K/s and 1.06 x 1013 K/s for the 20, 30 and 40 Å radius nanoparticles respectively. Structural Measures for Glass Formation Characterization of glassy behavior by molecular dynamics simulations is typically done using dynamic measurements such as the mean squared displacement, <r2(t)>. Liquids exhibit a mean squared displacement that is linear in time. Glassy materials deviate significantly from this linear behavior at intermediate times, entering a sub-linear regime with a return to linear behavior in the infinite time limit. Diffusion in nanoparticles differs significantly from the bulk in that atoms are confined to a roughly spherical volume and cannot explore any region larger than the particle radius. In these confined geometries, <r2(t)> in the radial direction approaches a limiting value of 6R2/40. However, glassy materials exhibit strong icosahedral ordering among nearest-neghbors in contrast to crystalline or liquid structures. Steinhart, et al., defined an orientational bond order parameter that is sensitive to the nearest-neighbor environment by using invariant combinations of spherical harmonics Yl,m(θ,φ).[10] Spherical harmonics involving the Y6,m(θ,φ) are particularly sensitive to icosohedral order among nearest neighbors as can be seen in the cartoon to the left. The second and third-order invariants, Q6 and W6 are used to determine the level of icosahedral order present in a quenched nanoparticle. Perfect icosahedral structures have a maximal value of 0.663 for Q6 and -0.170 for W6. A plot of the distributions of Q6 and W6 with cooling temperature indicates increasing icosahedral order with decreasing temperature. This is a clear indication that glassy structures are forming as the nanoparticles are quenched.

Combinations of Y6,m spherical harmonics indicating icosohedral ordering of the red atoms around a central (green) atom

Bond order distributions for 30 Å nanoparticle sampled at various temperatures

References [1] Meineke M. A.; Vardeman II, C. F.; Lin T.; Fennell C. J.; Gezelter J. D. J. Comp. Chem., 2005 26, 252. [2] Shibata T.; Bunker B.A.; Zhang Z.; Meisel D.; Vardeman II C. F.; Gezelter J. D. J. Am. Chem. Soc., 2002 124, 11989. [3] Foiles S. M.; Baskes M. I.; Daw M. S. Phys. Rev. B, 1986, 33, 7983 [4] Williams G.; Watts D. C. Trans. Faraday Soc., 1970, 66, 80. [5] Tu, K.-N.; Mayer, J. W. Electronic Thin Film Science; Macmillan: New York, 1992. [6] Hartland, G. V.; Hu, M.; Sader, J. E. J. Phys. Chem. B, 2003, 107, 7472 [7] Vardeman II, C. F.; Conforti, P. F.; Sprague, M. M.; Gezelter, J. D. J. Phys. Chem. B, 2005, 109, 16695 [8] Qhull, 1993. Software library is available from the National Science and Technology Research Center for Computation and Visualization of Geometric Structures, University of Minnesota. http://www.qhull.org [9] Vardeman II, C. F.; Gezelter, J. D. J. Phys. Chem. A, 2001, 105, 2568 [10] Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Phys. Rev. B, 1982, 28, 784