J. Phys. Chem. C 2008, 112, 3283-3293
3283
Simulations of Laser-Induced Glass Formation in Ag-Cu Nanoparticles Charles F. Vardeman II and J. Daniel Gezelter* Department of Chemistry and Biochemistry, UniVersity of Notre Dame, Notre Dame, Indiana 46556 ReceiVed: October 16, 2007; In Final Form: NoVember 30, 2007
Using molecular dynamics simulations, we have simulated the rapid cooling experienced by bimetallic nanoparticles following laser excitation at the plasmon resonance and find evidence that glassy beads, specifically Ag-Cu bimetallic particles at the eutectic composition (60% Ag, 40% Cu), can be formed during these experiments. The bimetallic nanoparticles are embedded in an implicit solvent with a viscosity tuned to yield cooling curves that match the experimental cooling behavior as closely as possible. Because the nanoparticles have a large surface-to-volume ratio, experimentally realistic cooling rates are accessible via relatively short simulations. The presence of glassy structural features was verified using bond orientational order parameters that are sensitive to the formation of local icosahedral ordering in condensed phases. As the particles cool from the liquid droplet state into glassy beads, a silver-rich monolayer develops on the outer surface and local icosahedra can develop around the silver atoms in this monolayer. However, we observe a strong preference for the local icosahedral ordering around the copper atoms in the particles. As the particles cool, these local icosahedral structures grow to include a larger fraction of the atoms in the nanoparticle, eventually leading to a glassy nanosphere.
1. Introduction Excitation of the plasmon resonance in metallic nanoparticles has attracted enormous interest in the past several years. This is partly due to the location of the plasmon band in the near IR for particles in a wide range of sizes and geometries. Living tissue is nearly transparent in the near IR, and for this reason there is an unrealized potential for metallic nanoparticles to be used in both diagnostic and therapeutic settings.1,2 One of the side effects of absorption of laser radiation at these frequencies is the rapid (subpicosecond) heating of the electronic degrees of freedom in the metal. This hot electron gas transfers heat to the phonon modes of the particle quickly, resulting in a rapid heating of the lattice of the metal particles. Because metallic nanoparticles have a large surface area to volume ratio, many of the metal atoms are at surface locations and experience relatively weak bonding. This is observable in a lowering of the melting temperatures of these particles when compared with bulk metallic samples.3,4 One of the side effects of the excitation of small metallic nanoparticles at the plasmon resonance is the facile creation of liquid metal droplets.5-10 Much of the experimental work on this subject has been carried out in the Hartland, El-Sayed, and Plech groups.6,7,9-13 These experiments mostly use the technique of time-resolved optical pump-probe spectroscopy, where a pump laser pulse serves to excite conduction band electrons in the nanoparticle and a following probe laser pulse allows observation of the time evolution of the electron-phonon coupling. Hu and Hartland have observed a direct relation between the size of the nanoparticle and the observed cooling rate using such pumpprobe techniques.14 Plech et al. have use pulsed X-ray scattering as a probe to directly access changes to atomic structure following pump excitation.9 They further determined that heat transfer in nanoparticles to the surrounding solvent is goverened by interfacial dynamics and not the thermal transport properties * Corresponding author. E-mail:
[email protected].
of the solvent. This is in agreement with Cahill,15 but opposite to the conclusions in ref 14. Because these experiments are carried out in condensed-phase surroundings, the large surface area to volume ratio makes the heat transfer to the surrounding solvent a relatively rapid process. In our recent simulation study of the laser excitation of gold nanoparticles,16 we observed that the cooling rate for these particles (1011-1013 K/s) is in excess of the cooling rate required for glass formation in bulk metallic alloys.17 Given this fact, it may be possible to use laser excitation to melt, alloy, and quench metallic nanoparticles in order to form glassy nanobeads. To study whether or not glass nanobead formation is feasible, 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.18 The Hume-Rothery rules suggest that alloys composed of copper and silver should be miscible in the solid state because their lattice constants are within 15% of each another.19 Experimentally, however, Ag-Cu alloys are a wellknown exception to this rule and are only miscible in the liquid state given equilibrium conditions.20 Below the eutectic temperature of 779 °C and composition (60.1% Ag, 39.9% Cu), the solid alloys of Ag and Cu will phase separate into Ag and Cu rich R and β phases, respectively.21,22 This behavior is due to a positive heat of mixing in both the solid and liquid phases. For the one-to-one composition fcc solid solution, ∆Hmix is on the order of +6 kJ/mol.22 Nonequilibrium solid solutions may be formed by undercooling, and under these conditions, a compositionally disordered γ fcc phase can be formed.23 Metastable alloys composed of Ag-Cu were first reported by Duwez in 1960 and were created by using a “splat quenching’’ technique in which a liquid droplet is propelled by a shock wave against a cooled metallic target.24 Because of the small positive ∆Hmix, supersaturated crystalline solutions are typically obtained rather than an amorphous phase. Higher ∆Hmix systems, such as Ag-Ni, are immiscible even in liquid
10.1021/jp710063g CCC: $40.75 © 2008 American Chemical Society Published on Web 02/01/2008
3284 J. Phys. Chem. C, Vol. 112, No. 9, 2008 states, but they tend to form metastable alloys much more readily than Ag-Cu. If present, then the amorphous Ag-Cu phase is usually seen as the minority phase in most experiments. Because of this unique crystalline-amorphous behavior, the Ag-Cu system has been widely studied. Methods for creating such bulkphase structures include splat quenching, vapor deposition, ion beam mixing, and mechanical alloying. Both structural25 and dynamic18 computational studies have also been performed on this system. Although bulk Ag-Cu alloys have been studied widely, this alloy has been mostly overlooked in nanoscale materials. The literature on alloyed metallic nanoparticles has dealt with the Ag-Au system, which has the useful property of being miscible in both solid and liquid phases. Nanoparticles of another miscible system, Au-Cu, have been constructed successfully using techniques such as laser ablation26 and the synthetic reduction of metal ions in solution.27 Laser-induced alloying has been used as a technique for creating Au-Ag alloy particles from core-shell particles.12 To date, attempts at creating AgCu nanoparticles have used ion implantation to embed nanoparticles in a glass matrix.28,29 These attempts have been largely unsuccessful in producing mixed alloy nanoparticles and instead produce phase-segregated or core-shell structures. One of the more successful attempts at creating intermixed Ag-Cu nanoparticles used alternate pulsed laser ablation and deposition in an amorphous Al2O3 matrix.30 Surface plasmon resonance (SPR) of bimetallic core-shell structures typically show two distinct resonance peaks where mixed particles show a single shifted and broadened resonance.11 The SPR for pure silver occurs at 400 nm and for copper at 570 nm.31 On Al2O3 films, these resonances move to 424 and 572 nm for the pure metals. For bimetallic nanoparticles with 40% Ag, an absorption peak is seen between 400 and 550 nm. With increasing Ag content, the SPR shifts toward the blue, with the peaks nearly coincident at a composition of 57% Ag. Gonzalo et al. cited the existence of a single broad resonance peak as evidence of an alloyed particle rather than a phase-segregated system. However, spectroscopy may not be able to tell the difference between alloyed particles and mixtures of segregated particles. High-resolution electron microscopy has so far been unable to determine whether the mixed nanoparticles were an amorphous phase or a supersaturated crystalline phase. Characterization of glassy behavior by molecular dynamics simulations is typically done using dynamic measurements such as the mean squared displacement, 〈r 2(t)〉. Liquids exhibit a mean squared displacement that is linear in time (at long times). Glassy materials deviate significantly from this linear behavior at intermediate times, entering a sublinear regime with a return to linear behavior in the infinite time limit.32 However, 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 (R). In these confined geometries, 〈r 2(t)〉 approaches a limiting value of 3R2/40.33 This limits the utility of dynamical measures of glass formation when studying nanoparticles. However, glassy materials exhibit strong icosahedral ordering among nearest neghbors (in contrast with crystalline and liquidlike configurations). Local icosahedral structures are the threedimensional equivalent of covering a two-dimensional plane with five-sided tiles; they cannot be used to tile space in a periodic fashion and are therefore an indicator of nonperiodic packing in amorphous solids. Steinhart et al. defined an orientational bond order parameter that is sensitive to icosahedral
Vardeman and Gezelter
Figure 1. Methodology used to mimic the experimental cooling conditions of a hot nanoparticle surrounded by a solvent. Atoms in the core of the particle evolved under Newtonian dynamics, whereas atoms that were in the outer skin of the particle evolved under Langevin dynamics. The radius of the spherical region operating under Newtonian dynamics, rNewton, was set to be 4 Å smaller than the original radius (R) of the liquid droplet.
ordering.34 This bond order parameter can therefore be used to characterize glass formation in liquid and solid solutions.35 Theoretical molecular dynamics studies have been performed on the formation of amorphous single-component nanoclusters of either gold,36-38 or nickel,39,40 by rapid cooling(∼1012-1013K/ s) from a liquid state. All of these studies found icosahedral ordering in the resulting structures produced by this rapid cooling, which can be evidence of the formation of an amorphous structure.41 The nearest-neighbor information was obtained from pair correlation functions, common neighbor analysis, and bond order parameters.34 It should be noted that these studies used single-component systems with cooling rates that are only obtainable in computer simulations and particle sizes less than 20 Å. Single-component systems are known to form amorphous states in small clusters42 but do not generally form amorphous structures in bulk materials. Because the nanoscale Ag-Cu alloy has been largely unexplored, many interesting questions remain about the formation and properties of such a system. Does the large surface area to volume ratio aid Ag-Cu nanoparticles in rapid cooling and formation of an amorphous state? Nanoparticles have been shown to have a size-dependent melting transition (Tm),3,4 so we might expect a similar trend to follow for the glass transition temperature (Tg). By analogy, bulk metallic glasses exhibit a correlation between Tm and Tg although such a dependence is difficult to establish because of the dependence of Tg on cooling rate and the process by which the glass is formed.43 It has also been demonstrated that there is a finite size effect depressing Tg in polymer glasses in confined geometries.44 In the sections below, we describe our modeling of the laser excitation and subsequent cooling of the particles in silico to mimic real experimental conditions. The simulation parameters have been tuned to the degree possible to match experimental conditions, and we discusss both the icosahedral ordering in the system as well as the clustering of icosahedral centers that we observed.
Simulations of Laser-Induced Glass Formation
J. Phys. Chem. C, Vol. 112, No. 9, 2008 3285
2. Computational Methodology 2.1. Initial Geometries and Heating. Cu-core/Ag-shell and random alloy structures were constructed on an underlying FCC lattice (4.09 Å) at the bulk eutectic composition Ag6Cu4. Both initial geometries were considered although experimental results suggest that the random structure is the most likely structure to be found following synthesis.30,45 Three different sizes of nanoparticles corresponding to a 20 Å radius (2382 atoms), 30 Å radius (6603 atoms), and 40 Å radius (15 683 atoms) were constructed. These initial structures were relaxed to their equilibrium structures at 20 K for 20 ps and again at 300 K for 100 ps sampling from a Maxwell-Boltzmann distribution at each temperature. All simulations were conducted using the OOPSE molecular dynamics package.46 To mimic the effects of the heating due to laser irradiation, we allowed the particles to melt by sampling velocities from the Maxwell Boltzmann distribution at a temperature of 900 K. The particles were run under microcanonical simulation conditions for 1 ns of simualtion time prior to studying the effects of heat transfer to the solvent. In all cases, the center of mass translational and rotational motion of the particles were set to zero before any data collection was undertaken. Structural features (pair distribution functions) were used to verify that the particles were indeed liquid droplets before cooling simulations took place. 2.2. Modeling Random Alloy and Core-Shell Particles in Solution-Phase Environments. To approximate the effects of rapid heat transfer to the solvent following a heating at the plasmon resonance, we utilized a methodology in which atoms contained in the outer 4 Å radius of the nanoparticle evolved under Langevin dynamics
m
x ∂ 2b ) Fsys(x b(t)) - 6π aηV b(t) + Fran ∂t2
(1)
with a solvent friction (η) approximating the contribution from the solvent and capping agent. Atoms located in the interior of the nanoparticle evolved under Newtonian dynamics. The setup of our simulations is nearly identical to the “stochastic boundary molecular dynamics’’ (SBMD) method that has seen wide use in the protein simulation community.47-49 A sketch of this setup can be found in Figure 1. In eq 1, the frictional forces of a spherical atom of radius a depend on the solvent viscosity. The random forces are usually taken as Gaussian random variables with zero mean and a variance tied to the solvent viscosity and temperature
〈Fran(t) Fran(t′)〉 ) 2kBT(6πη a)δ(t - t′)
(2)
Because of the presence of the capping agent and the lack of details about the atomic-scale interactions between the metallic atoms and the solvent, the effective viscosity is a essentially a free parameter that must be tuned to give experimentally relevant simulations. The viscosity (η) can be tuned by comparing the cooling rate that a set of nanoparticles experience with the known cooling rates for similar particles obtained via the laser heating experiments. Essentially, we tune the solvent viscosity until the thermal decay profile matches a heat-transfer model using reasonable values for the interfacial conductance and the thermal conductivity of the solvent. Cooling rates for the experimentally observed nanoparticles were calculated from the heat transfer equations for a spherical particle embedded in a ambient medium that allows for diffusive
heat transport. Following Plech et al.,9 we use a heat transfer model that consists of two coupled differential equations in the Laplace domain
McP‚(s‚Tp(s) - T0) + 4π R2G‚(Tp(s) - Tf (r ) R,s)) ) 0 (3)
(∂r∂ T (r,s)) f
r)R
+
G (T (s) - Tf (r,s)) ) 0 K p
(4)
where s is the time-conjugate variable in Laplace space. The variables in these equations describe a nanoparticle of radius R, mass M, and specific heat cp at an initial temperature T0. The surrounding solvent has a thermal profile Tf (r,t), thermal conductivity K, density F, and specific heat c. G is the interfacial conductance between the nanoparticle and the surrounding solvent and contains information about heat transfer to the capping agent as well as the direct metal-to-solvent heat loss. The temperature of the nanoparticle as a function of time can then obtained by the inverse Laplace transform
Tp(t) ) 2kR2g2T0 π
exp(-κu2t/R2)u2
∫0∞ (u2(1 + Rg) - kRg)2 + (u3 - kRgu)2 du
(5)
For simplicity, we have introduced the thermal diffusivity κ ) K/(Fc), and defined k ) 4πR3Fc/(Mcp) and g ) G/K in eq 5. Equation 5 was solved numerically for the Ag-Cu system using mole-fraction weighted values for cp and Fp of 0.295 (Jg-1K-1) and 9.826 × 106 (gm-3) respectively. Because most of the laser excitation experiments have been done in aqueous solutions, the parameters used for the fluid are K of 0.6(Wm-1K-1), F of 1.0 × 106 (gm-3), and c of 4.184(Jg-1K-1). Values for the interfacial conductance have been determined by a number of groups for similar nanoparticles and range from a low 87.5 × 106(Wm-2K-1) to 130 × 106(Wm-2K-1).15 Wilson et al. worked with Au, Pt, and AuPd nanoparticles and obtained an estimate for the interfacial conductance of G ) 130(Wm-2K-1).15 Similarly, Plech et al. reported a value for the interfacial conductance of G ) 105 ( 15(Wm-2K-1) for Au nanoparticles.9 We conducted our simulations at both ends of the range of experimentally determined values for the interfacial conductance. This allows us to observe both the slowest and fastest heat transfers from the nanoparticle to the solvent that are consistent with experimental observations. For the slowest heat transfer a value for G of 87.5 × 106 (Wm-2K-1) was used, and for the fastest heat transfer a value of 117 × 106 (Wm-2K-1) was used. On the basis of calculations we have done using raw data from the Hartland group’s thermal half-time experiments on Au nanospheres,50 the true G values are probably in the faster regime: 117 × 106(Wm-2K-1). The rate of cooling for the nanoparticles in a molecular dynamics simulation can then be tuned by changing the effective solvent viscosity (η) until the nanoparticle cooling rate matches the cooling rate described by the heat transfer, eq 4. The effective solvent viscosity (in Pa s) for a G of 87.5 × 106 (Wm-2K-1) is 4.2 × 10-6, 5.0 × 10-6, and 5.5 × 10-6 for 20, 30, and 40 Å particles, respectively. The effective solvent viscosity (again in Pa s) for an interfacial conductance of 117 × 106 (Wm-2K-1) is 5.7 × 10-6, 7.2 × 10-6, and 7.5 × 10-6 for 20, 30, and 40 Å particles. These viscosities are essentially gas-phase values, a fact that is consistent with the initial temperatures of the particles being well into the supercritical region for the aqueous environment. Gas bubble generation has
3286 J. Phys. Chem. C, Vol. 112, No. 9, 2008
Vardeman and Gezelter
Figure 2. Thermal cooling curves obtained from the inverse Laplace transform heat model in eq 4 (solid line) as well as from molecular dynamics simulations (circles). Effective solvent viscosities of 4.2-7.5 × 10-6 Pa s (depending on the radius of the particle) give the best fit to the experimental cooling curves. This viscosity suggests that the nanoparticles in these experiments are surrounded by a vapor layer (which is a reasonable assumption given the initial temperatures of the particles).
Figure 3. Illustrative cooling profile for the 40 Å nanoparticle evolving under stochastic boundary conditions corresponding to G ) 87.5 × 106(Wm-2K-1). At temperatures along the cooling trajectory, configurations were sampled and allowed to evolve in the NVE ensemble. These subsequent trajectories were analyzed for structural features associated with bulk glass formation.
also been seen experimentally around gold nanoparticles in water.51 Instead of a single value for the effective viscosity, a time-dependent parameter might be a better mimic of the cooling vapor layer that surrounds the hot particles. This may also be a contributing factor to the size dependence of the effective viscosities in our simulations. Cooling traces for each particle size are presented in Figure 2. It should be noted that the Langevin thermostat produces cooling curves that are consistent with Newtonian (singleexponential) cooling, which cannot match the cooling profiles from eq 5 exactly. Fitting the Langevin cooling profiles to a single-exponential produces τ ) 25.576 ps, τ ) 43.786ps, and τ ) 56.621 ps for the 20, 30, and 40 Å nanoparticles and a G of 87.5 × 106(Wm-2K-1). For comparison’s sake, similar singleexponential fits with an interfacial conductance of G of 117 ×
106 (Wm-2K-1) produced τ ) 13.391 ps, τ ) 30.426ps, and τ ) 43.857 ps for the 20, 30, and 40 Å nanoparticles. 2.3. Potentials for Classical Simulations of Bimetallic Nanoparticles. Several different potential models have been developed that reasonably describe interactions in transition metals. In particular, the Embedded Atom Model (EAM)52 and Sutton-Chen (SC)53 potential have been used to study a wide range of phenomena in both bulk materials and nanoparticles.18,33,54-57 Both potentials are based on a model of a metal that treats the nuclei and core electrons as pseudoatoms embedded in the electron density due to the valence electrons on all of the other atoms in the system. The sc potential has a simple form that closely resembles that of the ubiquitous Lennard-Jones potential
Simulations of Laser-Induced Glass Formation
Utot )
J. Phys. Chem. C, Vol. 112, No. 9, 2008 3287
[
]
1
DijV pair ∑i 2 ∑ ij (rij) - ci DiixFi j*i
(6)
where V pair ij and Fi are given by
V pair ij (r) )
() Rij rij
nij
,
Fi )
∑ j*i
() Rij rij
mij
(7)
V pair ij is a repulsive pairwise potential that accounts for interactions between the pseudoatom cores. The xFi term in eq 6 is an attractive many-body potential that models the interactions between the valence electrons and the cores of the pseudo-atoms. Dij and Dii set the appropriate overall energy scale, ci scales the attractive portion of the potential relative to the repulsive interaction, and Rij is a length parameter that assures a dimensionless form for F. These parameters are tuned to various experimental properties such as the density, cohesive energy, and elastic moduli for FCC transition metals. The quantum Sutton-Chen (Q-SC) formulation matches these properties while including zero-point quantum corrections for different transition metals.58 This particular parametarization has been shown to reproduce the experimentally available heat of mixing data for both FCC solid solutions of Ag-Cu and the hightemperature liquid.25 In contrast, the eam potential does not reproduce the experimentally observed heat of mixing for the liquid alloy.59 In this work, we have utilized the Q-SC formulation for our potential energies and forces. Combination rules for the alloy were taken to be the arithmetic average of the atomic parameters with the exception of ci because its value is only dependent on the identity of the atom where the density is evaluated. For the Q-SC potential, cutoff distances are traditionally taken to be 2Rij and include up to the sixth coordination shell in FCC metals. To better understand the structural changes occurring in the nanoparticles throughout the cooling trajectory, configurations were sampled at regular intervals during the cooling trajectory. These configurations were then allowed to evolve under NVE dynamics to sample from the proper distribution in phase space. Figure 3 illustrates this sampling. 3. Analysis Frank first proposed local icosahedral ordering of atoms as an explanation for supercooled atomic (specifically metallic) liquids and further showed that a 13-atom icosahedral cluster has a 8.4% higher binding energy either face-centered cubic (FCC) or hexagonal close-packed (HCP) crystal structures.60 Icosahedra also have six 5-fold symmetry axes that cannot be extended indefinitely in three dimensions; long-range translational order is therefore incommensurate with local icosahedral ordering. This does not preclude icosahedral clusters from possessing long-range orientational order. The “frustrated’’ packing of these icosahedral structures into dense clusters has been proposed as a model for glass formation.61 The size of the icosahedral clusters is thought to increase until frustration prevents any further growth.62 Molecular dynamics simulations of a two-component Lennard-Jones glass showed that clusters of face-sharing icosahedra are distributed throughout the material.63 Simulations of freezing of single-component metalic nanoclusters have shown a tendency for icosohedral structure formation particularly at the surfaces of these clusters.36,39,64 Experimentally, the splitting (or shoulder) on the second peak of the X-ray structure factor in binary metallic glasses has been
attributed to the formation of tetrahedra that share faces of adjoining icosahedra.65 Various structural probes have been used to characterize structural order in molecular systems including common neighbor analysis, Voronoi tesselations, and orientational bond-order parameters.66-69 The method that has been used most extensively for determining local and extended orientational symmetry in condensed phases is the bond-orientational analysis formulated by Steinhart et al.34 In this model, a set of spherical harmonics is associated with each of the near neighbors of a central atom. Neighbors (or “bonds’’) are defined as having a distance from the central atom that is within the first peak in the radial distribution function. The spherical harmonic between a central atom i and a neighboring atom j is
Ylm(b r ij) ) Ylm(θ(b r ij),φ(b r ij))
(8)
where, {Ylm(θ,φ)} are spherical harmonics and θ(r b) and φ(r b) are the polar and azimuthal angles made by the bond vector b r with respect to a reference coordinate system. We chose for simplicity the origin as defined by the coordinates for our nanoparticle. (Only even-l spherical harmonics are considered because permutation of a pair of identical particles should not affect the bond-order parameter.) The local environment surrounding atom i can be defined by the average over all neighbors, Nb(i), surrounding that atom
qjlm(i) )
Nb(i)
1 Nb(i)
Ylm(b r ij) ∑ j)1
(9)
We can further define a global average orientational-bond order over all qjlm by calculating the average of qjlm(i) over all N particles N
Q h lm )
Nb(i) qjlm(i) ∑ i)1 (10)
N
Nb(i) ∑ i)1 The Q h lm contained in eq 10 is not necessarily invariant under rotations of the arbitrary reference coordinate system. Secondand third-order rotationally invariant combinations, Ql and Wl, can be taken by summing over m values of Q h lm
Ql )
(
l
4π
∑
2l + 1 m)-l
| |)
1/2
Q h lm
2
(11)
and
W ˆl)
Wl
(∑| |) l
Q h lm
2
(12)
3/2
m)-l
where
Wl )
∑ m ,m ,m 1
2
3
(
)
l l l m1 m2 m3 × Qlm1Qlm2Qlm3
(13)
m1+m2+m3)0
The factor in parentheses in eq 13 is the Wigner-3j symbol, and the sum is over all valid (|m| el) values of m1, m2, and m3 which sum to zero.
3288 J. Phys. Chem. C, Vol. 112, No. 9, 2008
Vardeman and Gezelter
Figure 4. Cutaway views of 30 Å Ag-Cu nanoparticle structures for random alloy (top) and Cu (core)/Ag (shell) initial conditions (bottom). Shown from left to right are the crystalline, liquid droplet, and final glassy bead configurations.
Figure 5. Distributions of the bond orientational order parameter (W ˆ 6) at different temperatures. The upper, middle, and lower panels are for 20, 30, and 40 Å particles, respectively. The left-hand column used cooling rates commensurate with a low interfacial conductance (87.5 × 106 Wm-2K-1), whereas the right-hand column used a more physically reasonable value of 117 × 106 Wm-2K-1. The peak at W ˆ 6 ≈ -0.17 is due to local icosahedral structures. The different curves in each of the panels indicate the distribution of W ˆ 6 values for samples taken at different times along the cooling trajectory. The initial and final temperatures (in K) are indicated on the plots adjacent to their respective distributions.
For ideal face-centered cubic (FCC), body-centered cubic (BCC), and simple cubic (SC) as well as hexagonally closepacked (hcp) structures, these rotationally invariant bond order parameters have fixed values independent of the choice of coordinate reference frames. For ideal icosahedral structures, the l ) 6 invariants, Q6 and W ˆ 6 are also independent of the coordinate system. Q4 and W ˆ 4 will have nonvanishing values for indiVidual icosahedral clusters, but these values are not invariant under rotations of the reference coordinate systems. Similar behavior is observed in the bond-orientational order parameters for individual liquid-like structures. Additionally, both Q6 and W ˆ 6 are thought to have extreme values for the icosahedral clusters.34 This makes the l ) 6 bond-orientational order parameters particularly useful in identifying the extent of local icosahedral ordering in condensed phases. For example, a local structure that exhibits W ˆ 6 values near -0.17 is easily
identified as an icosahedral cluster and cannot be mistaken for distorted cubic or liquid-like structures. One may use these bond orientational order parameters as an averaged property to obtain the extent of icosahedral ordering in a supercooled liquid or cluster. It is also possible to accumulate information about the distributions of local bond orientational order parameters, p(W ˆ 6) and p(Q6), which provide information about individual atomic sites that are central to local icosahedral structures. The distributions of atomic Q6 and W ˆ 6 values are plotted as a function of temperature for our nanoparticles in Figures 5 and 6. At high temperatures, the distributions are unstructured and are broadly distributed across the entire range of values. As the particles are cooled, however, there is a dramatic increase in the fraction of atomic sites that have local icosahedral ordering around them. (This corresponds to the sharp peak appearing in
Simulations of Laser-Induced Glass Formation
J. Phys. Chem. C, Vol. 112, No. 9, 2008 3289
Figure 6. Distributions of the bond orientational order parameter (Q6) at different temperatures. The curves in the six panels in this figure were computed at conditions identical to those of the same panels in Figure 5.
local minima once the temperature has fallen below 528 K. As the temperature falls, it is possible for substructures to become trapped in the local icosahedral well, and if the cooling is rapid enough this trapping leads to vitrification. A similar analysis of the free energy surface for orientational order in bulk glass formers can be found in the work of van Duijneveldt and Frenkel.70 We have also calculated the fraction of atomic centers that have strong local icosahedral order
ficos )
Figure 7. Free energy as a function of the orientational order parameter (W ˆ 6) for 40 Å bimetallic nanoparticles as they are cooled from 902 to 310 K. As the particles cool below 528 K, a local minimum in the free energy surface appears near the perfect icosahedral ordering (W ˆ6 ) -0.17). At all temperatures, liquid-like structures are a global minimum on the free energy surface; but if the cooling rate is fast enough, substructures may become trapped with local icosahedral order, leading to the formation of a glass.
Figure 5 at W ˆ 6 ) -0.17 and to the broad shoulder appearing in Figure 6 at Q6 ) 0.663.) The probability distributions of local order can be used to generate free energy surfaces using the local orientational ordering as a reaction coordinate. By making the simple statistical equivalence between the free energy and the probabilities of occupying certain states
g(W ˆ 6) ) -kBT ln p(W ˆ 6)
(14)
we can obtain a sequence of free energy surfaces (as a function of temperature) for the local ordering around central atoms within our particles. Free energy surfaces for the 40 Å particle at a range of temperatures are shown in Figure 7. Note that at all temperatures the liquid-like structures are global minima on the free energy surface, whereas the local icosahedra appear as
∫-∞w p(Wˆ 6)dWˆ 6 i
(15)
where wi is a cutoff value in W ˆ 6 for atomic centers that are displaying icosahedral environments. We have chosen a (somewhat arbitrary) value of wi ) -0.15 for the purposes of this work. A plot of ficos(T) as a function of temperature of the particles is given in Figure 8. As the particles cool, the fraction of local icosahedral ordering rises smoothly to a plateau value. The smaller particles (particularly the ones that were cooled in a higher viscosity solvent) show a slightly larger tendency toward icosahedral ordering. Because we have atomic-level resolution of the local bondorientational ordering information, we can also look at the local ordering as a function of the identities of the central atoms. In Figure 9, we display the distributions of W ˆ 6 values for both the silver and copper atoms and we note a strong predilection for the copper atoms to be central to icosahedra. This is probably due to local packing competition of the larger silver atoms around the copper, which would tend to favor icosahedral structures over the more densely packed cubic structures. The locations of these icosahedral centers are not distrubted uniformly throughout the particles. In Figure 10 we show snapshots of the centers of the local icosahedra (i.e., any atom that exhibits a local bond orientational order parameter W ˆ6 < -0.15). At high temperatures, the icosahedral centers are transitory, existing only for a few femtoseconds before being reabsorbed into the liquid droplet. As the particle cools, these centers become fixed at certain locations and additional icosahedra develop throughout the particle, clustering around the sites where the structures originated. There is a strong preference for icosahedral ordering near the surface of the particles.
3290 J. Phys. Chem. C, Vol. 112, No. 9, 2008
Vardeman and Gezelter
Figure 8. Temperautre dependence of the fraction of atoms with local icosahedral ordering, ficos(T) for 20, 30, and 40 Å particles cooled at two different values of the interfacial conductance.
Figure 9. Distributions of the bond orientational order parameter (W ˆ 6) for the two different elements present in the nanoparticles. This distribution was taken from the fully cooled 40 Å nanoparticle. Local icosahedral ordering around copper atoms is much more prevalent than that around silver atoms.
Figure 10. Centers of local icosahedral order (W ˆ 6 < 0.15) at 900, 471, and 315 K for the 30 Å nanoparticle cooled with an interfacial conductance G ) 87.5 × 106 Wm-2K-1. Silver atoms (blue) exhibit icosahedral order at the surface of the nanoparticle, whereas copper icosahedral centers (green) are distributed throughout the nanoparticle. The icosahedral centers appear to cluster together, and these clusters increase in size with decreasing temperature.
Identification of these structures by the type of atom shows that the silver-centered icosahedra are evident only at the surface of the particles.
In contrast with the silver ordering behavior, the copper atoms that have local icosahedral ordering are distributed more evenly throughout the nanoparticles. Figure 11 shows this tendency as
Simulations of Laser-Induced Glass Formation
J. Phys. Chem. C, Vol. 112, No. 9, 2008 3291
Figure 11. Appearance of icosahedral clusters around central silver atoms is largely due to the presence of these silver atoms at or near the surface of the nanoparticle. The upper panel shows the fraction of icosahedral atoms ( ficos(r) for each of the two metallic atoms as a function of distance from the center of the nanoparticle (r). The lower panel shows the radial density of the two constituent metals (relative to the overall density of the nanoparticle). Icosahedral clustering around copper atoms are more evenly distributed throughout the particle, whereas icosahedral clustering around silver is largely confined to the silver atoms at the surface.
TABLE 1: Values of Bond Orientational Order Parameters for Simple Structures Cooresponding to l ) 4 and l ) 6 Spherical Harmonic Functions35 a fcc hcp bcc sc icosahedral (liquid)
Q4
Q6
W ˆ4
W ˆ6
0.191 0.097 0.036 0.764
0.575 0.485 0.511 0.354 0.663
-0.159 0.134 0.159 0.159
-0.013 -0.012 0.013 0.013 -0.170
a Q4 and W ˆ 4 have values for indiVidual icosahedral clusters, but these values are not invariant under rotations of the reference coordinate systems. Similar behavior is observed in the bond-orientational order parameters for individual liquid-like structures.
a function of distance from the center of the nanoparticle. Silver, because it has a lower surface free energy than copper, tends to coat the skins of the mixed particles.71 This is true even for bimetallic particles that have been prepared in the Ag (core)/ Cu (shell) configuration. Upon forming a liquid droplet, approximately 1 monolayer of Ag atoms will rise to the surface of the particles. This can be seen visually in Figure 4 as well as in the density plots in the bottom panel of Figure 11. This observation is consistent with previous experimental and theoretical studies on bimetallic alloys composed of noble metals.72-74 Bond order parameters for surface atoms are averaged only over the neighboring atoms, so packing constraints that may prevent icosahedral ordering around silver in the bulk are removed near the surface. It would certainly be interesting to see if the relative tendency of silver and copper to form local icosahedral structures in a bulk glass differs from our observations on nanoparticles. Similar behavior has been observed by Luo et al. in their work on amorphous Ni-Ag alloys. They used a common neighbor analysis (CNA) technique that identified icosahedral ordering from simulated structures that match experimental EXAFS spectra. Their simulated structures exhibited icosahedral structures that were nearly always centered on the smaller Ni
TABLE 2: Estimates of the Glass Transition Temperatures (Tg) for Three Different Sizes of Bimetallic Ag6Cu4 Nanoparticles Cooled under Two Different Values of the Interfacial Conductance, G radius (Å)
interfacial conductance
effective cooling rate (× 1013K/s)
Tg (K)
20 20 30 30 40 40
87.5 117 87.5 117 87.5 117
2.4 4.5 1.3 1.9 1.0 1.3
477 502 491 493 476 487
atoms in the sample.75 Details of the common neighbor analysis technique can be found in Sheng et al.’s work on the glass transition in bulk Ag-Cu alloys.25 In the bulk Ag-Cu alloys, high quench rates do lead to an increase in icosahedral ordering, although the onset is much more gradual than what we have observed in the bimetallic nanoparticles. Sheng et al. also estimated the glass transition temperature (Tg) in bulk Ag-Cu alloys by locating a discontinuity in the slope of the average atomic volume, 〈V〉/N, or enthalpy when plotted against the temperature of the alloy. They obtained a bulk glass transition temperature of Tg ) 510 K for a quenching rate of 2.5 × 1013 K/s. For simulations of nanoparticles, there is no periodic box, and therefore no facile way of exactly computing the volume. Instead, we estimate the volume of our nanoparticles using Barber et al.’s very fast quickhull algorithm to obtain the convex hull for the collection of 3-d coordinates of all of atoms at each point in time.76,77 The convex hull is the smallest convex polyhedron that includes all of the atoms, so the volume of this polyhedron is an excellent estimate of the volume of the nanoparticle. This method of estimating the volume will be problematic if the nanoparticle breaks into pieces (i.e., if the bounding surface becomes concave); but for the relatively short trajectories used in this study, it provides an excellent measure of particle volume as a function of time (and temperature).
3292 J. Phys. Chem. C, Vol. 112, No. 9, 2008 Using the discontinuity in the slope of the average atomic volume versus temperature, we arrive at an estimate of Tg that is approximately 488 K. We note that this temperature is somewhat below the onset of icosahedral ordering exhibited in the bond orientational order parameters. It appears that icosahedral ordering is initiated while the system is still somewhat fluid and is locked in place once the temperature falls below Tg. We did not observe any dependence of our estimates for Tg on either the nanoparticle size or the value of the interfacial conductance. However, the cooling rates and size ranges we utilized are all sampled from a relatively narrow range, and it is possible that much larger particles would have substantially different values for Tg. Our estimates for the glass transition temperatures for all three particle sizes and both interfacial conductance values are shown in Table 2. 4. Conclusions Our heat-transfer calculations have utilized the best current estimates of the interfacial heat transfer coefficient (G) from recent experiments. Using reasonable values for thermal conductivity in both the metallic particle and the surrounding solvent, we have obtained cooling rates for laser-heated nanoparticles that are in excess of 1013 K/s. To test whether or not this cooling rate can form glassy nanoparticles, we have performed a mixed molecular dynamics simulation in which the atoms in contact with the solvent or capping agent are evolved under Langevin dynamics while the remaining atoms are evolved under Newtonian dynamics. The effective solvent viscosity (η) is a free parameter that we have tuned so the particles in the simulation follow the same cooling curve as their experimental counterparts. From the local icosahedral ordering around the atoms in the nanoparticles (particularly Copper atoms), we deduce that it is likely that glassy nanobeads are created via laser heating of bimetallic nanoparticles, particularly when the initial temperature of the particles approaches the melting temperature of the bulk metal alloy. Improvements to our calculations would require (1) explicit treatment of the capping agent and solvent, (2) another radial region to handle the heat transfer to the solvent vapor layer that is likely to form immediately surrounding the hot particle,14,51 and (3) larger particles in the size range most easily studied via laser heating experiments. The local icosahedral ordering we observed in these bimetallic particles is centered almost completely around the copper atoms, and this is likely due to the size mismatch leading to a more efficient packing of five-membered rings of silver around a central copper atom. This size mismatch should be reflected in bulk calculations, and work is ongoing in our lab to confirm this observation in bulk glass-formers. The physical properties (bulk modulus, frequency of the breathing mode, and density) of glassy nanobeads should be somewhat different from their crystalline counterparts. However, observation of these differences may require single-particle resolution of the ultrafast vibrational spectrum of one particle both before and after the crystallite has been converted into a glassy bead. Acknowledgment. We thank Greg Hartland, Dani Meisel and Kristina Furse for reading through early versions of this paper. Greg Hartland was also kind enough to supply us with raw data on the thermal relaxation half-times of gold nanospheres in aqueous environments. Support for this project was provided by the National Science Foundation under grant CHE0134881.
Vardeman and Gezelter References and Notes (1) West, J.; Halas, N. Annu. ReV. Biomed. Eng. 2003. (2) Hu, M.; Chen, J.; Li, Z.-Y.; Au, L.; Hartland, G. V.; Li, X.; Marquez, M.; Xia, Y. Chem. Soc. ReV. 2006. (3) Buffat, P.; Borel, J.-P. Phys. ReV. A 1976, 13, 2287-2298. (4) Dick, K.; Dhanasekaran, T.; Zhang, Z.; Meisel, D. J. Am. Chem. Soc. 2002, 124, 2312-2317. (5) Mafune, F.; Kohno, J.; Takeda, Y.; Kondow, T. J. Phys. Chem. B 2001, 105, 9050-9056. (6) Hartland, G.; Hu, M.; Sader, J. J. Phys. Chem. B 2003, 107, 74727478. (7) Link, S.; El-Sayed, M. A. Int. ReV. Phys. Chem. 2000, 19, 409453. (8) Plech, A.; Kurbitz, S.; Berg, K.; Graener, H.; Berg, G.; Gresillon, S.; Kaempfe, M.; Feldmann, J.; Wulff, M.; von Plessen, G. Europhys. Lett. 2003, 61, 762-768. (9) Plech, A.; Kotaidis, V.; Gresillon, S.; Dahmen, C.; von Plessen, G. Phys. ReV. B 2004, 70, 195423. (10) Plech, A.; Cerna, R.; Kotaidis, V.; Hudert, F.; Bartels, A.; Dekorsy, T. Nano Lett. 2007, 7, 1026-1031. (11) Hodak, J. H.; Henglein, A.; Giersig, M.; Hartland, G. V. J. Phys. Chem. B 2000, 104, 11708-11718. (12) Hartland, G.; Guillaudeu, S.; Hodak, J. Mol. Compon. Electron. DeVices 2003. (13) Petrova, H.; Hu, M.; Hartland, G. V. Z. Phys. Chem. 2007, 221, 361-376. (14) Hu, M.; Petrova, H.; Hartland, G. V. Chem. Phys. Lett. 2004, 391, 220-225. (15) Wilson, O.; Hu, X.; Cahill, D.; Braun, P. Phys. ReV. B 2002, 66. (16) Vardeman, C.; Conforti, P.; Sprague, M.; Gezelter, J. J. Phys. Chem. B 2005, 109, 16695-16699. (17) Greer, A. L. Science 1995, 267, 1947-1953. (18) Vardeman, II. C. F.; Gezelter, J. D. J. Phys. Chem. A 2001, 105, 2568. (19) Kittel, C. Introduction to Solid State Physics, 7th ed.; Wiley: New York, 1996. (20) Massalski, T. B.; Murray, J. L.; Bennett, L. H.; Baker, H. Binary Alloy Phase Diagrams; American Society for Metals: Metals Park, OH, 1986. (21) Banhart, J.; Ebert, H.; Kuentzler, R.; Voitla¨nder, J. Phys. ReV. B 1992, 46, 9968-9975. (22) Ma, E. Prog. Mater Sci. 2005, 50, 413-509. (23) Najafabadi, R.; Srolovitz, D. J.; Ma, E.; Atzmon, M. J. Appl. Phys. 1993, 74, 3144-3149. (24) Duwez, P.; Willens, R. H.; Klement, W. Jr. J. Appl. Phys. 1960, 31, 1136-1137. (25) Sheng, H. W.; He, J. H.; Ma, E. Phys. ReV. B 2002, 65, 184203. (26) Malyavantham, G.; O’Brien, D. T.; Becker, M. F.; Keto, J. W.; Kovar, D. J. Nanopart. Res. 2004, 6, 661-664. (27) Kim, M.; Na, H.; Lee, K. C.; Yoo, E. A.; Lee, M. J. Mater. Chem 2003, 13, 1789-1792. (28) De, G.; Gusso, M.; Tapfer, L.; Catalano, M.; Gonella, F.; Mattei, G.; Mazzoldi, P.; Battaglin, G. J. Appl. Phys. 1996, 80, 6734-6739. (29) Magruder, III. R. H.; Osborne, Jr. D. H.; Zuhr, R. A. J. Non-Cryst. Solids 1994, 176, 299-303. (30) Gonzalo, J.; Babonneau, D.; Afonso, C. N.; Barnes, J.-P. J. Appl. Phys. 2004, 96, 5163-5168. (31) Henglein, A. J. Phys. Chem. B 2000, 104, 1206-1211. (32) Kob, W. J. Phys. Condens. Matter 1999, 11, R85-R115. (33) Shibata, T.; Bunker, B.; Zhang, Z.; Meisel, D.; Vardeman, C.; Gezelter, J. J. Am. Chem. Soc. 2002, 124, 11989-11996. (34) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Phys. ReV. B 1983, 28, 784-804. (35) ten Wolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. J. Chem. Phys. 1996, 104, 9932-9947. (36) Chen, Y.; Bian, X.; Zhang, J.; Zhang, Y.; Wang, L. Modell. Simul. Mater. Sci. Eng. 2004, 12, 373-379. (37) Cleveland, C. L.; Landman, U.; Schaaff, T. G.; Shafigullin, M. N.; Stephens, P. W.; Whetten, R. L. Phys. ReV. Lett. 1997, 79, 1873-1876. (38) Cleveland, C. L.; Landman, U.; Shafigullin, M. N.; Stephens, P. W.; Whetten, R. L. Z. Phys. D 1997, 40, 503-508. (39) Gafner, Y. Y.; Gafner, S. L.; Entel, P. Phys. Solid State 2004, 46, 1327-1330. (40) Qi, Y.; Cagin, T.; Johnson, W. L. III. W. A. G. J. Chem. Phys. 2001, 115, 385-394. (41) Strandburg, K. J. Bond-Orientational Order in Condensed Matter Systems; Springer-Verlag: New York, 1992. (42) Breaux, G. A.; Cao, B.; Jarrold, M. F. J. Phys. Chem. B 2005, 109, 16575-16578. (43) Wang, W.; Wen, P.; Zhao, D.; Pan, M.; Wang, R. J. Mater. Res. 2003, 18, 2747-2751.
Simulations of Laser-Induced Glass Formation (44) Alcoutlabi, M.; McKenna, G. J. Phys.: Condens. Matter 2005, 17, R461-R524. (45) Jiang, H.; sik Moon, K.; Wong, C. P. Proc. Int. Symp. AdV. Packag. Mater. 2005, 173-177. (46) Meineke, M. A.; Vardeman, II. C. F.; Lin, T.; Fennell, C. J.; Gezelter, J. D. J. Comput. Chem. 2005, 26, 252-271. (47) Brooks, C.; Brunger, A.; Karplus, M. Biopolymers 1985, 24, 843865. (48) Brooks, C.; Karplus, M. J. Chem. Phys. 1983, 79, 6312-6325. (49) Brunger, A.; Brooks, C.; Karplus, M. Chem. Phys. Lett. 1984, 105, 495-500. (50) Hu, M.; Hartland, G. J. Phys. Chem. B 2002, 106, 7029-7033. (51) Kotaidis, V.; Dahmen, C.; von Plessen, G.; Springer, F.; Plech, A. J. Chem. Phys. 2006, 124, 184702. (52) Foiles, S. M.; Baskes, M. I.; Daw, M. S. Phys. ReV. B 1986, 33, 7983-7991. (53) Sutton, A. P.; Chen, J. Philos. Mag. Lett. 1990, 61, 139-146. (54) Sankaranarayanan, S.; Bhethanabotla, V.; Joseph, B. Phys. ReV. B 2005, 71. (55) Chui, Y.; Chan, K. Phys. Chem. Chem. Phys. 2003, 5, 2869-2874. (56) Wang, G.; Van Hove, M.; Ross, P.; Baskes, M. J. Phys. Chem. B 2005, 109, 11683-11692. (57) Medasani, B.; Park, Y. H.; Vasiliev, I. Phys. ReV. B 2007, 75. (58) Qi, Y.; C¸ ain, T.; Kimura, Y.; Goddard, III. W. A. Phys. ReV. B 1999, 59, 3527-3533. (59) Murray, J. L. Metall. Trans. 1984, 15, 261-268. (60) Frank, F. C. Proc. R. Soc. London, Ser. A 1952, 215, 43-46. (61) Steinhardt, P. J. Science 1987, 238, 1242-1247.
J. Phys. Chem. C, Vol. 112, No. 9, 2008 3293 (62) Hoare, M. Ann. N.Y. Acad. Sci. 1976, 279, 186-207. (63) Jo´nsson, H.; Andersen, H. C. Phys. ReV. Lett. 1988, 60, 22952298. (64) Nam, H.-S.; Hwang, N. M.; Yu, B. D.; Yoon, J.-K. Phys. ReV. Lett. 2002, 89, 275502. (65) van de Waal, B. W. J. Non-Cryst. Solids 1995, 189, 118-128. (66) Honeycutt, J. D.; Andersen, H. C. J. Phys. Chem. 1987, 91, 49504963. (67) Iwamatsu, M. Mater. Sci. Eng., A 2007, 449-451, 975-978. (68) Hsu, C. S.; Rahman, A. J. Chem. Phys. 1979, 71, 4974-4986. (69) Nose, S.; Yonezawa, F. J. Chem. Phys. 1986, 84, 1803-1814. (70) van Duijneveldt, J. S.; Frenkel, D. J. Chem. Phys. 1992, 96, 46554668. (71) Zhu, L.; DePristo, A. E. J. Catal. 1997, 167, 400-407. (72) Mainardi, D.; Balbuena, P. Langmuir 2001, 17, 2047-2050. (73) Huang, S.-P.; Balbuena, P. J. Phys. Chem. B 2002, 106, 72257236. (74) Ramirez Caballero, G. E.; Balbuena, P. B. Mol. Simul. 2006, 32, 297-303. (75) Luo, W. K.; Sheng, H. W.; Alamgir, F. M.; Bai, J. M.; He, J. H.; Ma, E. Phys. ReV. Lett. 2004, 92, 145502. (76) Barber, C. B.; Dobkin, D. P.; Huhdanpaa, H. T. ACM Trans. Math. Software 1996, 22, 469-483. (77) “Qhull’’, 1993 software library is available from the National Science and Technology Research Center for Computation and Visualization of Geometric Structures (The Geometry Center), University of Minnesota. http://www.geom.umn.edu/software/qhull/.