Open Access
Add to BibSonomyAbstract
Understanding the formation of electrodynamically interacting assemblies of metal nanoparticles requires accurate computational methods for determining the forces and propagating trajectories. However, since computation of electromagnetic forces occurs on attosecond to femtosecond timescales, simulating the motion of colloidal nanoparticles on milliseconds to seconds timescales is a challenging multi-scale computational problem. Here, we present a computational technique for performing accurate simulations of laser-illuminated metal nanoparticles. In the simulation, we self-consistently combine the finite-difference time-domain method for electrodynamics (ED) with Langevin dynamics (LD) for the particle motions. We demonstrate the ED-LD method by calculating the 3D trajectories of a single 100-nm-diameter Ag nanoparticle and optical trapping and optical binding of two and three 150-nm-diameter Ag nanoparticles in simulated optical tweezers. We show that surface charge on the colloidal metal nanoparticles plays an important role in their optically driven self-organization. In fact, these simulations provide a more complete understanding of the assembly of different structures of two and three Ag nanoparticles that have been observed experimentally, demonstrating that the ED-LD method will be a very useful tool for understanding the self-organization of optical matter.
© 2015 Optical Society of America
1. Introduction
A focus of nanophotonics is to control light on length scales below the diffraction limit, which can be achieved by plasmonic interactions among metal nanoparticles [1–4]. Ordered assemblies of metal nanoparticles generate collective optical properties that can be tuned via the spatial arrangements of the particles for various nanophotonics applications [5, 6]. Light itself can be used to create self-organized systems of colloidal metal nanoparticles [7–10] via optical trapping and optical binding [11–13]. The creation of light-controlled particle arrangements, commonly known as optical matter [14], is an active area of research [7, 15, 16]. The particles in optical matter are confined in a light field by the forces arising from the gradient of the intensity [7, 12] and the phase [13, 17] of the light beam, while the specific spatial arrangements among the particles result from the interparticle electrodynamic interactions, also known as the optical binding forces, as well as the overall boundary conditions of the incident optical field and other material boundaries (e.g., glass cover slides) [12]. On the other hand, Brownian motion of particles opposes the order and assembly of optical matter thereby making the self-organization of an extended optical matter system challenging. To understand the assembly and the stability of optical matter it is essential to study the coupled mechanics and electrodynamics for experimentally realistic conditions.
The dynamics of optical trapping have been studied computationally by Ng et al. [18], Gordon et al. [19], and more recently by Volpe and Volpe [20], and Wang et al. [21]. However, a complete Langevin dynamics description of colloidal metal nanoparticles with a self-consistent full electrodynamic description of the optical trapping beam is still lacking. The electrodynamic description of the trapping beam and its interaction with the trapped particles have commonly been treated with approximations such as extended Mie theory [18, 21], generalized multipole techniques [19], or by simple harmonic potentials [20]. However, these approximations would not be suitable for simulations of complex focused beams and irregularly shaped particles.
In this paper, we describe a new computational scheme to perform realistic simulations of colloidal metal nanoparticles illuminated by a focused laser beam. It combines the finite-difference time-domain (FDTD) method [22] to solve Maxwell’s equations and simulate the electrodynamics (ED) of metal nanoparticles with full Langevin dynamics (LD) [23], to simulate particle motions. We illustrate the ED-LD method by simulating optical trapping of a single 100-nm-diameter Ag nanoparticle in a focused Gaussian beam (i.e., optical tweezers) and by simulating optical trapping and optical binding of two and three 150-nm-diameter Ag nanoparticles. In agreement with experiments [12, 24, 25], our simulations demonstrate stronger confinement of a Ag nanoparticle in a focused Gaussian beam when it is trapped closer to the focal region due to larger intensity gradients. Our simulations also reproduce the different experimentally observed self-organized structures formed by two and three Ag nanoparticles [12]. We find that the optical binding separation between two metal nanoparticles depends on the properties of the illuminating beam, with focused beams resulting in shorter optical binding distances. Furthermore, our simulations indicate that the Coulomb forces, due to surface charges on the metal nanoparticles, can play an important role in the optical self organization. These results demonstrate that our ED-LD method provides realistic simulations and insights into the nature of optically induced motion and interactions of metallic nanoparticles.
2. The coupled ED-LD method
2.1. The finite-difference time-domain (FDTD) technique
FDTD is a commonly used numerical method [22] in which the time-dependent Maxwell’s curl equations are discretized using centered differences in both space and time and solved using leapfrog time integration yielding a fully explicit second-order accurate scheme. The E⃗, and H⃗ field components are spatially staggered on a grid cell such that each E⃗ (H⃗) component is surrounded by circulating H⃗ (E⃗) components, thereby implicitly enforcing Gauss’ law. Moreover, efficient solutions of absorbing boundaries and plane-wave sourcing have been implemented in FDTD via the convolutional perfectly matched layer (CPML) absorbing boundary condition and the total-field/scattered field (TF/SF) technique, respectively [22]. We model the dispersive and lossy metal nanoparticles in FDTD using an auxiliary differential equation approach with a single-pole Drude model for the conductivity of the metal. Although, a simple Drude model does not yield a good fit to experimentally measured permittivities of metals like Au and Ag [26] over a wide frequency range, one can achieve very good agreement for a small range (provided that the fitting parameters are not interpreted with any physical significance) and therefore achieve accurate modeling of the metal nanoparticle [27].
We calculate the optical forces using the Maxwell stress tensor (MST) formalism following the description in Nome et al. [28]. The time-averaged optical force, where the average is taken over several optical cycles, is given by the surface integral of the MST,
where T̄ is the MST (Strictly speaking, Eq. (1)b is the Minkowski tensor [29] since it contains the dielectric constant of the surrounding medium), n̂ is a unit vector normal to the surface S, Ec (Hc) represents the complex-valued phasor electric (magnetic) field, ⊗ denotes the outer product between the phasor quantities, and Ī is the identity tensor. We validated our FDTD simulation and optical force calculation by computing the force between two 50-nm-diameter spherical Au nanoparticles separated by 16 nm in water illuminated by a linearly polarized plane wave (polarization direction along the interparticle axis; 550 nm wavelength in vacuum) with an intensity of 5 MW/cm2. Using a 1 nm grid-cell, the same as Ref. [28], our calculation yielded a force of −1.76 pN, which is in very good agreement with the −1.78 pN obtained by Nome et al. [28]. [The Au Drude parameters used for this validation simulation were ε∞ = 13.26, ωp = 103.7 eV, and Γp = 0.3104 eV.]Optical trapping and optical binding experiments have been carried out with variously shaped, focused beams [11, 24, 30]. While there are approaches to encode such beams within standard FDTD procedures [31], we have chosen a conceptually and algoritmatically far simpler approach. We introduce a lossless dielectric plano-convex lens in the path of the plane wave injected by the TF/SF method. The radius of curvature and the permittivity of the lens are chosen to minimize reflections of the plane wave from the lens and to achieve the required focal length. The numerical aperture (NA) of the lens in our simulations is typically in the range of 1 to 1.25. The lens used in the simulation does not represent a physically viable focusing scheme; a realistic lens is much larger than our simulation domain. However, this lens-in-simulation-volume approach allows a simple and relatively computationally inexpensive approximation of a realistic lens due to the freedom of choosing its refractive index and thus obtaining experimentally appropritate NAs. The simplicity of this approach over that in Ref. [31] is a major asset due to the ease with which one could parallelize the FDTD simulation. We do not use a more rigorous scheme for describing strongly focused beams such as the vector Debye method [32] because they are not easily compatible with the TF/SF method in FDTD. A schematic of the simulation geometry is shown in Fig. 1(a). We assume the metal nanoparticle to be in water, bounded on top and bottom by the lens and a glass slide, respectively. The incident wave is propagated from top to bottom in Fig. 1(b)–1(d). Figures 1(b) and 1(c) show a time snapshot of the electric field intensity distribution in the simulation domain, demonstrating the focusing effect of the lens, and with a Ag nanoparticle, showing the near-field enhancement, respectively.
Fig. 1 Illustration of the simulation domain and electric field intensity from a focused incident wave. Schematic representing the geometry of the simulation showing (a) the xy and (b) the yz plane. The total-field scattered-field (TF/SF) and convolutional perfectly matched layer (CPML) boundaries are also shown. An example of the electric field intensity (electric field polarized along y) distribution in the yz plane (c) showing the focusing effect of the lens and (d) with a Ag nanoparticle showing the near-field intensity enhancement. The axes in panels (c) and (d) represent grid cells numbers. A grid cell is a cubic box with each side of length 5 nm.
2.2. The Langevin dynamics (LD) simulation
The LD equation of motion for each nanoparticle is [23]
Here, r⃗ is the (center-of-mass) position of a nanoparticle, F⃗ is the force, λv is the friction coefficient given by Stokes’ law (λv = 6πνR, where R is the radius of the sphere and ν is the dynamic viscosity of the medium), and η is the stochastic thermal noise term with a Gaussian probability distribution. Note that in the case of two or more nanoparticles, F⃗ may include not only optical forces but electrostatic interactions between particles. The fluctuation-dissipation theorem connects λv and η such that 〈ηi(t)ηj(t′)〉 = 2λkBTδijδ(t − t′), where the angle brackets denote an average over time, δij denotes a Kronecker delta function over the coordinate indices i, j, and k, and δ(t − t′) is a delta function of time [33]. Equation (2) describes the Brownian motion of the colloidal metal nanoparticles.We solve the LD equation [Eq. (2)] using the scheme (the BAOAB method) described in Refs. [34] and [35]. In the BAOAB scheme, the LD equation of motion is separated into two first order differential equations:
where p⃗ is the momentum of the nanoparticle, γ = λv/m, and . Rn is a 3D vector containing normally distributed random numbers. We use the Box-Muller [36] method to generate Rn from a set of random numbers distributed uniformly between 0 and 1. In BAOAB, Eqs. (3)a and (3)b are solved by using a splitting-method time-integration scheme: Eqs. (3) are split into three parts labelled A, B, and O corresponding to the kinetic term, potential term, and friction-thermal fluctuation term, respectively. BAOAB denotes the sequence of solving the constituent A, B, and O equations, where the repeated parts (A, B) are integrated over a half time step, resulting in the following update equations for r⃗ and p⃗ [35]:We find the BAOAB integration scheme to be more accurate and stable than a simple centered-difference solution of Eq. (2). The LD time step is chosen such that the displacement of the nanoparticle during a time step is not more than 1/10th of its diameter. A typical LD time step in our simulations is between 0.1 μs and 0.5 μs. In simulations containing two or more nanoparticles, the collisions among them are resolved by assuming the nanoparticles to be hard spheres [37], although weaker (e.g., Lennard-Jones) interactions could easily be introduced.
2.3. Computational details
The FDTD and LD methods are coupled in time as shown in the flow chart in Fig. 2 and the algorithm is implemented in FORTAN 90. We start the LD time-stepping loop following the initialization of the simulation domain and the assignment of the material properties (relative permittivity, conductivity, etc.) to the FDTD grid. In the LD loop, the positions and velocities of the nanoparticles are updated according to the total force acting on them, i.e. using Eqs. (4)a–d. The new positions are then used to reassign the grid-based material properties. The FDTD time-stepping loop is started subsequently. The continuous-wave incident fields are propagated in the FDTD loop, allowing the fields to evolve to a steady-state and then electrodynamic forces are recalculated [Eq. (1)]. The updated forces are then used in the LD loop to drive the nanoparticles [Eq. (4)e followed by Eqs. (4)a–d] and this process is repeated until the desired duration of the simulation.
Fig. 2 Flow chart of the coupled ED-LD method. The outer loop contains the LD update equations. The inner loop containing the FDTD update equations is run every third time step of the outer loop.
Since, FDTD requires the most computer time to update the fields and calculate the optical forces, we calculate the updated optical forces in FDTD every third time step of the LD loop. Having separate loops for the LD and ED, and by calculating the ED forces every third time step amounts to a quasi-static-force approximation for particle motion, i.e. the forces on the particles are assumed to be constant during the LD time step. For the source intensities and the LD time steps used in the simulations here, the displacement of a metal nanoparticle every LD time step is less than 1/10th of its diameter (we consider nanoparticles of 100–150 nm diameter), whereas the wavelength of the trapping beam is roughly 600 nm in water (800 nm in vacuum). As a result, the error in the ED forces and therefore the nanoparticle trajectories due to this quasi-static-force approximation is small, thus yielding accurate mean positions and distributions in steady-state. (Note that the FDTD method solves the full time-dependent form of Maxwell’s equations taking into account retardation effects and the quasi-static-force approximation is not to be confused with the quasi-static approximation in electrodynamics.) It should also be noted that the over-damped or Brownian approximation [i.e. m → 0 in Eq. (2)] is often used as a valid description of nanoparticle motion in solutions. However, we use the full LD equation since the main computational burden in our simulations is due to the FDTD part and the over-damped dynamics approximation would not provide any noticeable increase in computational speed. Moreover, because of the large intensity gradients in the optical near fields, ignoring the inertial terms could lead to errors in the trajectories.
For the FDTD method in our simulation, we use a 3D computational domain with roughly 200 × 200 × 200 grid cells including 15-cell-thick CPML absorbing boundaries on all sides. We use a grid cell size of 1/25th of the diameter of the metal nanoparticle. We assume a linearly polarized source wave (i.e. trapping laser) with a wavelength of 800 nm in vacuum. A continuous-wave sinusoidal function with an initial Gaussian increase is used for the incident wave to eliminate high frequency components due to any abruptness at the start of the TF/SF source injection [i.e. we use E⃗s(t) = ŷE0 exp[(n − 3n0)/n0]2 sin[ω(n − n0)Δt − kz] for t < n0Δt, and E⃗s(t) = ŷE0 sin[ω(n − n0)Δt − kz] for t ← n0Δt. n0 is the rise time of the Gaussian, which is set to one period of the optical cycle]. Depending on the size of the simulation domain in the z direction, we calculate the optical forces by averaging over 10–17 time periods of the incident wave. The surface integral of the MST is performed over a box that is three grid cells larger than the diameter of the metal nanoparticle. The motion of the nanoparticle is restricted within a box that is 5–10 grid cells smaller than the TF/SF region in the x and y dimensions (denoted by white dashed line in Figs. 3 and 6) and 5–10 grid cells below and above the lens and glass interface, respectively, in the z dimension. When a nanoparticle hits the boundaries of this box, it undergoes an elastic reflection. In our simulations here, we assume Ag nanoparticles. The single-pole Drude model parameters for the Ag nanoparticles are found by fitting the experimental dielectric function [26] for wavelength in the range of 750–850 nm. We obtain the following Drude parameters: ε∞ = 3.045, ωp = 81.08 eV, and Γp = 0.063 eV. The simulations presented here were performed on a single node with one CPU (2.53–2.60 GHz) and took up to 300 hours of CPU time while requiring less than 2 GB of memory each. The ED-LD method can be made significantly faster by parallelizing the FDTD computations [38] while keeping the LD computations serial. Alternative approaches to calculate the fields surrounding the nanoparticle such as the transition matrix method [39] could be used in lieu of FDTD, possibly making the simulation faster. Nevertheless, the FDTD method is more generally applicable and therefore could be used in future simulations involving particles with sharper features, such as bipyramids, nanotriangles, and nanoplates.
Fig. 3 Results of an ED-LD simulation of a 100-nm-diameter Ag nanoparticle in a focused Gaussian optical beam. (a) Trajectory of the nanoparticle with (maroon trace) and without a source (black trace) (see Visualization 1). The electric field [polarized along y] intensity of the focused Gaussian beam in an xy plane 80 nm from the water-glass interface is shown as the colored background. (b) Optical forces mapped in the small region of the xy plane for the Ag nanoparticle when illuminated with the focused Gaussian beam.
Laser-illuminated metal nanoparticles generate localized heat due to their enhanced absorption near the plasmon resonance. The temperature increase due to localized heating for typical source intensities was estimated to be about 20 K [13, 40]. To account for the temperature increase, we use the value of dynamic viscosity of water corresponding to 320 K, i.e. ν = 0.6 mPa-s, which affects the nanoparticle motion via the viscous drag force coefficient γ in Eq. (3)b. Hydrodynamic interactions also affect the dynamics of colloidal nanoparticles; for example in Ref. [19] an optohydrodynamic theory was used to treat the interactions of the nanoparticles in optical traps and the hydrodynamic interactions were shown to affect the rate of approach to the traps. However, the hydrodynamic interactions do not affect the final stable configurations [19]. Since the focus of our computations is to simulate the self-organization of metal nanoparticles into different structures that have been observed in experiments [12, 41–43] and to understand the effects of different focused beams and surface charges on the final configurations with reasonably accurate particle trajectories, in this work we ignore the corrections to the particle dynamics that can be obtained by including hydrodynamic interactions.
3. Results
We illustrate the ED-LD method by simulating the motion of Ag nanoparticles illuminated by a focused laser beam. We first examine optical trapping of a single Ag nanoparticle in a focused Gaussian beam (optical tweezers) and then investigate the optical trapping and optical binding of two and three Ag nanoparticles. For simulations involving multiple nanoparticles, we consider different incident beams: plane wave, focussed Gaussian, and line trap. Moreover, in multi-nanoparticle simulations, we also consider the effects of the surface charge on the colloidal metal nanoparticles.
3.1. Optical trapping of a single Ag nanoparticle
The results of an optical trapping simulation of a single 100-nm-diameter spherical Ag nanoparticle in a focused Gaussian beam are shown in Fig. 3 [also see Visualization 1]. The focusing lens has a focal length that is 107 nm farther from the water-glass interface inside the glass region [see Fig. 1 for the simulation geometry]. Figure 3(a) shows the trajectory of a Ag nanoparticle in the xy plane without any incident wave (black trace) and with the focused Gaussian beam (maroon trace). The electric field intensity in an xy plane 80 nm from the water-glass interface is shown as the colored background. The total duration of the simulation for both the trapped and the free particle is 1.5 ms (corresponding to 3000 LD time steps for a LD time step of 0.5 μs). It is clearly seen that when illuminated by the focused incident beam, the Ag nanoparticle is confined in the xy plane, whereas, without the source, it undergoes free Brownian motion. Figure 3(b) shows a zoomed-in view of the positions of the trapped Ag nanoparticle(i.e. with focused Gaussian source) as circles. The optical forces at each of these positions are shown by arrows, with the length of a given arrow representing its magnitude and the direction being the direction of the force at that time step. As expected, the optical forces point towards the center of the xy plane (i.e., the center of the beam axis at x = 0 and y = 0). In the z direction, radiation pressure pushes the Ag nanoparticle towards the water-glass interface (in the direction of propagation of the incident wave). When the particle collides with the bounding surface, it undergoes elastic reflection. In steady-state along the z direction, the nanoparticle bounces back and forth from the glass-water interface, with its mean position being 80 nm from the interface. The center-of-mass positions of the nanoparticle in the xz and yz planes are shown in Figs. 4(a) and 4(b). The glass surface in this simulation is located at the z = 600 nm plane. Figure 4(c) shows the potentials of mean force [Eq. (5a)] as a function of the distance of the center of the nanoparticle from the water-glass interface.
Fig. 4 Trapping along the z direction at the water-glass interface. In this simulation, the glass surface is located at the z = 600 nm plane. Center-of-mass positions of a 100-nm-diameter Ag particle in a focused Gaussian beam in (a) the xz plane and (b) the yz plane. (c) Potentials of mean force (pmf) in units of kBT as a function of the distance of the center of the nanoparticle from the water-glass interface. The solid curve is a parabolic fit to the pmf.
Consistent with the experimental results [12,24,25], we find that the confinement of the particle is stronger when the focal region of the source is farther from the water-glass interface. We vary the focal length of the lens by changing its refractive index. We obtain 2-dimensional histograms of positions in the xy plane of a 100-nm-diameter Ag nanoparticle, where the distance between the focal region and the water-glass interface in decreased. Trapping in the xy plane is more efficient when the focal region is closer to the trapping plane (i.e. the water-glass interface) due to the increase in the intensity gradient forces. Using the probability density function from the position histograms, one can calculate the potentials of mean force (pmf), which represents an effective trapping potential, and the trap force constant (κ), which represents the positional stiffness of the trap [44]:
Here, the parameter ξ is the position along x, y, or z, P(ξ) is the probability density of ξ, and the angled brackets denote an average. The pmf is shown in Figs. 5(a)–5(c) in units of kBT. The stronger confinement correlates with deeper trapping potentials. The trap force constants κ calculated for the three different focused Gaussian beams in Figs. 5(a)–5(c) are 3.8 pN/μm, 5.6 pN/μm, and 6.3 pN/μm, respectively.Fig. 5 Two-dimensional potentials of mean force determined from histograms of the positions calculated for a 100-nm-diameter Ag nanoparticle illuminated by a focused Gaussian beam for varying focal lengths. The distance between the water-glass interface and focal length of the lens are (a) 202 nm, (b) 107 nm, and (c) 32 nm. Refractive index of the lens (nl) and the corresponding focal lengths (f) are (a) nl = 2.0, f = 1110 nm, (b) nl = 2.062, f = 1015 nm, and (c) nl = 2.121, f = 940 nm. The corresponding peak intensity of the electric field at the water-glass interface is 2.5 MW/cm2, 5 MW/cm2, and 8 MW/cm2, respectively.
Fig. 6 Optical trapping and binding of two 150-nm-diameter Ag nanoparticles. (a) The trajectories of the two nanoparticles in a line trap along x (maroon and pink traces) and without any incident field (black and grey traces) (see Visualization 2). The electric field (polarized along y) intensity of the line trap in the xy plane at the glass-water interface is shown as the colored background. (b) Histogram (bin width 5 nm) of the separation between the two Ag nanoparticles in the line trap (maroon) and without any incident wave (black) calculated from a simulation of 2 ms total duration.
3.2. Optical trapping and binding of two Ag nanoparticles
Next, we simulate optical trapping and optical binding of two Ag nanoparticles in a line trap created by an anisotropically focused plane wave. A line trap is simulated by using a cylindrical lens for focusing the incident TF/SF plane wave. The axis of the cylindrical lens, and therefore that of the line trap, is along the x direction and two 150-nm-diameter Ag nanoparticles are initialized with a separation of 550 nm between their centers, close to the left boundary. Figure 6(a) shows the trajectories of the two Ag nanoparticles calculated from a total simulation duration of 2 ms [also see Visualization 2] using a LD time step of 0.5 μs. Without any incident electromagnetic field, the two nanoparticles undergo free Brownian motion, as seen from the black and gray traces representing their trajectories. On the other hand, when the TF/SF source is turned on and is focused by the cylindrical lens, the two nanoparticles (maroon and pink traces) are trapped along the y dimension due to the strong gradient force created by focusing. Moreover, the two nanoparticles in the line trap drift together towards the center of the trap in the x direction demonstrating optical binding between them. The separation between the two nanoparticles without any incident field (black) and with the line trap (maroon) is plotted in Fig. 6(b) as a histogram showing a strong peak for the line trap at the optical binding separation of 550 nm.
Figure 7(a) shows the paths of the two Ag nanoparticles illuminated by a plane wave (green and light green traces) and the line trap (maroon and pink traces). We find that in both cases the two nanoparticles move together as a rigid body demonstrating optical binding. However, the steady-state average optical binding distance between the nanoparticles depends on the shape of the incident beam. In Fig. 7(b), we plot the separation between the two nanoparticle centers as a histogram while Fig. 7(c) shows the pmf as a function of the separation from which we infer that the plane-wave illuminated nanoparticles have an average optical binding separation of 595 nm compared to the 550 nm for the line trap. In Fig. 7(c), the solid lines are parabolic fits to the potentials of mean force. As expected due to focusing, the line trap shows a deeper trapping potential compared to the plane wave.
Fig. 7 Comparing optical-binding separations of two 150-nm-diameter Ag nanoparticles. (a) Trajectories of the nanoparticles under plane-wave illumination (green and light green traces) and in a line trap (maroon and pink traces) starting from the same pair of initial positions. (b) Histogram of the separation between the two Ag nanoparticles and (c) Potential of mean force calculated as a function of the separation for the different incident waves (Solid lines are parabolic fits). (d) Average separation between two Ag nanoparticles in a line trap (maroon squares) and under plane wave illumination (green circle) as a function of the ζ-potential. Error bars denote standard deviation from the mean.
3.3. Effect of surface charge on metal nanoparticle assembly
To mimic the colloidal Ag nanoparticles used in experiments [12], we simulated the surface charge on the nanoparticles by assigning a total charge, Qnp, to each nanoparticle based on the following parameters; Qnp = 4πR2σ where, σ = ε0εrζ/λD, is the surface charge density [45] and λD = [ε0εrkBT/2e2NAI]1/2 is the Debye length [46]. Here, ζ is the zeta-potential, I is the ionic concentration, e is the elementary charge, and NA is Avogadro’s number. The measured ζ-potential that is typical in our experiments is in the range of −20 mV to −50 mV [12,47] and we assume negatively charged monovalent ions with ionic concentration in the range of 0.01 mmol/L to 0.2 mmol/L. In Fig. 7(d), we show the average separation between two nanoparticles with different ζ-potentials (maroon squares) in a line trap and with plane wave illumination (green circle). The error bars denote the standard deviation from the mean. We find that the average optical-binding separation between two 150-nm-diameter Ag nanoparticles in a line trap depends on the total charge on their surface for the ζ-potential values found in actual samples. The optical-binding separation for these nanoparticles computed with no surface charge is 550 nm, in agreement with previous calculations [12]. However, the experimentally measured separation is close to 600 nm [12]. By taking into account the center-to-center Coulomb repulsion between the nanoparticles (by varying the ζ-potential in our simulations), we obtain an optical-binding separation of 600 nm for two particles.
We also examined the effect of surface charge on the structures formed by two Ag nanoparticles in a focussed Gaussian beam [Note that surface charge is simulated by a center-to-center electrostatic Coulomb interaction among the nanoparticles]. We find that the surface charge plays an important role when the optical forces on the nanoparticles are attractive towards the center of the beam, while the Coulomb interaction between them is repulsive. Figures 8(a)–8(b) show the trajectories of two 150-nm-diameter Ag nanoparticles in a focused Gaussian beam, where the total charge Qnp on the nanoparticles in panel (b) is three times that in (a). We assume the same ζ-potential of −50 mV, whereas we use different ionic concentrations [0.02 mmol/L for (a) and 0.18 mmol/L for (b)]. In both simulations, the Ag nanoparticles are initialized at the same locations (shown by the X marks). The black circles (diameter 150 nm) represent the mean positions of the nanoparticles in steady-state after the particles are trapped (i.e. after the mean position becomes localized). Other than the surface charge, we used exactly the same initial conditions (we also used the same seed for the random number generator). Clearly, the nanoparticles form completely different final arrangements due to the effect of the repulsive Coulomb interaction. In Fig. 8(a), the optical gradient forces win over the Coulomb repulsion between the nanoparticles, attracting them both towards the center of the focused beam. At small interparticle separations (< 300 nm), the near-field interactions start to dominate, which results in the nanoparticles aligning themselves along the polarization direction (y-axis) of the electric field. On the other hand, in Fig. 8(b), Coulomb repulsion dominates, preventing the nanoparticles from getting closer; ultimately, the nanoparticles are trapped at approximately the first optical binding separation [i.e. at d ≈ λ/n (n=1.33)]. Analogous simulations can be used to understand the experimental work of Bao et al. [41] and from Refs. [42, 43] in which the surface charge of the nanoparticles and the substrate was shown to have an important effect on their positions upon optical printing.
Fig. 8 ED-LD simulation of two 150-nm-diameter Ag nanoparticles in a focused Gaussian beam. Trajectories and final assembly of the two nanoparticles (blue and maroon lines) assuming ionic concentrations of (a) 0.02 mmol/L and (b) 0.18 mmol/L. The mean trapped center-of-mass position of (a) 200 nm and (b) 525 nm is shown by the black circles (diameter 150 nm) and initial position (same in both panels) is shown by the X-marks. The direction of polarization of the incident beam is along the y axis. The LD time step used for these simulations was 0.25 μs and the total duration was 1 ms.
3.4. Optical self-assembly of three Ag nanoparticles
Finally, we show that the structures that form with multiple particles depend on a balance of all forces present; intensity gradient from the trap size and shape, optical binding, and Coulomb. We simulated the dynamics of three 150-nm-diameter Ag nanoparticles in a focussed Gaussian beam to examine their steady-state self-organization. We assume the ζ-potential to be −50 mV and vary the ionic concentration I from 0.01 mmol/L to 0.08 mmol/L in steps of 0.01 resulting in Qnp to be in the range of −1.75 × 10−18 C to −4.95 × 10−18 C. Figures 9(a)–9(h) show the trajectories of the three Ag nanoparticles with the different Qnp but with the same initial positions. The axis of the focused Gaussian beam is centered at the x = 0, y = 0 point. The X marks show the initial positions of the nanoparticles, while the circles (diameter 150 nm) represent the average positions after trapping (i.e. after the average positions become persistent). We use a LD time step of 0.1 μs and the trajectories represent a total simulated duration of 0.5 ms. Again, we find that the self-organized structure depends on the strength of the Coulomb interaction between the nanoparticles. Moreover, these results show a crossover from a regime where the Coulomb interaction is weak and the self-assembled structure is dominated by the optical near-field effect [I < 0.02 mmol/L, Figs. 9(a) and 9(b)] to a regime where the repulsive Coulomb forces are stronger, leading to a structure created by the optical binding interactions [I > 0.07 mmol/L, Figs. 9(g) and 9(h)]. In the intermediate-Coulomb-interaction regime, [I ≈ 3–4 mmol/L, Figs. 9(c) and 9(d)], the final arrangement revealed from our simulations exactly mirrors the experimental work reported in [12] [see Fig. 2a – 3] that involved 150 nm Ag nanoparticles and a focused Gaussian beam. In our simulation, the trapping plane along the z direction is only 200 nm from the focal region. As a result, we find that the intensity gradient forces are stronger than the optical binding interaction, and therefore, the self-organized structure is not centered with the beam axis (x = 0, y = 0 in Fig. 9) [In the experimental observation in [12] the focal plane is much farther (few microns) from the trapping plane, presumably resulting in the centered arrangement]. These results further demonstrate the importance of surface charge in optical matter structures. In fact, these subtleties of electrostatic versus electrodynamic forces were not fully appreciated in our earlier experimental work [12, 13, 47].
Fig. 9 ED-LD simulation of three 150-nm-diameter Ag nanoparticles in a focused Gaussian optical beam. The X marks represent the initial positions, which are same in all panels. The circle of radius 75 nm, represents the average position of the nanoparticle (averaged over the last 1000 LD time steps). The ζ-potential is fixed at −50 mV, while the ionic concentration (I, in units of mmol/L) is varied to simulate different surface charge conditions: (a) I = 0.01, (b) I = 0.02, (c) I = 0.03, (d) I = 0.04, (e) I = 0.05, (f) I = 0.06, (g) I = 0.07, and (h) I = 0.08. The electric field of the incident beam is polarized along the y axis. Note that the initial positions are identical in all the panels.
4. Summary and conclusions
We have presented the ED-LD method, a rigorous computational technique for performing accurate simulations of metal nanoparticles illuminated by different focused laser beams. The simulation combines fully electrodynamic solutions of Maxwell’s equations using the FDTD method with the full Langevin dynamics description of the motion of colloidal metal nanoparticles. We demonstrate the ED-LD method by simulating optical trapping of a single 100-nm-diameter Ag nanoparticle. The position density of a Ag nanoparticle trapped at different distances from the focal region calculated from our simulation is consistent with experimental results [12, 24, 25]. We also simulated the motion of two 150-nm-diameter Ag nanoparticles in a line trap and in plane-wave illumination. We demonstrate optical binding of the nanoparticles and show that the optical binding separation depends on the shape of the incident beam, with shorter optical binding separation resulting from focused line traps as compared to a plane wave. We also examine the effect of inter-particle Coulomb interactions on the self-organization of two and three 150-nm-diameter Ag nanoparticles, showing that in a focused Gaussian beam, the surface charge and the resulting Coulomb interactions play an important role in determining the steady-state positions and arrangements supporting various experimental observations [12, 41–43]. From these results, we can better appreciate the surface charge effects on optical matter assembly, which are often ignored in simulations [18–21].
The simulations presented here demonstrate the potential of the ED-LD method for understanding the self-organization of metal nanoparticles in optical matter systems. In particular, simulations of realistic systems are possible by taking into account the shape of the optical beam and surface charge on the nanoparticles. The advantage of the FDTD method lies in its ability to compute the optical forces acting on irregularly shaped, lossy, and dispersive metal nanoparticles and to treat complex shaped beams. Combined with LD, one can calculate accurate trajectories of these metal nanoparticles. Although, accurate simulations of large optical matter assemblies of metal nanoparticles are possible using the ED-LD method, the computational requirements scale rapidly with increasing size of the simulation domain making simulations of more than three nanoparticles prohibitive on a single processor. Efficient parallelization of the FDTD method could in principle allow simulations with more nanoparticles and is the focus of our future work. We also plan to add hydrodynamic interactions in future studies that threat dynamically driven nanoparticle systems. Understanding and characterizing the dynamical behavior of optical matter assembly could lead to design of more stable and novel structures of optical matter with potential applications in photonics and nano-optics.
Acknowledgments
We thank Patrick Figliozzi for the information about zeta-potentials relevant to the experimental conditions. This work was supported by a grant from the Department of Defense (National Security Science and Engineering Faculty Fellowship). Computer time at the Center for Nanoscale Materials, an Office of Science user facility, was supported by the U. S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357.
References and links
1. V. Giannini, A. I. Fernández-Domínguez, Y. Sonnefraud, T. Roschuk, R. Fernández-García, and S. A. Maier, “Controlling light localization and light–matter interactions with nanoplasmonics,” Small 6, 2498–2507 (2010). [CrossRef] [PubMed]
2. P. Wang, B. Huang, Y. Dai, and M.-H. Whangbo, “Plasmonic photocatalysts: Harvesting visible light with noble metal nanoparticles,” Phys. Chem. Chem. Phys. 14, 9813–9825 (2012). [CrossRef] [PubMed]
3. T. Jensen, L. Kelly, A. Lazarides, and G. Schatz, “Electrodynamics of noble metal nanoparticles and nanoparticle clusters,” J. Clust. Sci. 10, 295–317 (1999). [CrossRef]
4. P. Nordlander, C. Oubre, E. Prodan, K. Li, and M. I. Stockman, “Plasmon hybridization in nanoparticle dimers,” Nano Lett. 4, 899–903 (2004). [CrossRef]
5. O. M. Marago, P. H. Jones, P. G. Gucciardi, G. Volpe, and A. C. Ferrari, “Optical trapping and manipulation of nanostructures,” Nat Nano 8, 807–819 (2013). [CrossRef]
6. A. Wei, Metal Nanoparticle Ensembles: Collective Optical Properties (Taylor and Francis, 2009).
7. T. Čižmár, L. C. D. Romero, K. Dholakia, and D. L. Andrews, “Multiple optical trapping and binding: new routes to self-assembly,” J. Phys. B: At. Mol. Opt. Phys. 43, 102001 (2010). [CrossRef]
8. T. M. Grzegorczyk, J. Rohner, and J.-M. Fournier, “Optical mirror from laser-trapped mesoscopic particles,” Phys. Rev. Lett. 112, 023902 (2014). [CrossRef] [PubMed]
9. V. Karásek, T. Čižmár, O. Brzobohatý, P. Zemánek, V. Garcés-Chávez, and K. Dholakia, “Long-range one-dimensional longitudinal optical binding,” Phys. Rev. Lett. 101, 143601 (2008). [CrossRef] [PubMed]
10. C. D. Mellor, T. A. Fennerty, and C. D. Bain, “Polarization effects in optically bound particle arrays,” Opt. Express 14, 10079–10088 (2006). [CrossRef] [PubMed]
11. V. Demergis and E.-L. Florin, “Ultrastrong optical binding of metallic nanoparticles,” Nano Lett. 12, 5756–5760 (2012). [CrossRef] [PubMed]
12. Z. Yan, S. Gray, and N. F. Scherer, “Potential energy surfaces and reaction pathways for light-mediated self-organization of metal nanoparticle clusters,” Nat. Commun. 5, 3751 (2014). [CrossRef] [PubMed]
13. Z. Yan, M. Sajjan, and N. F. Scherer, “Fabrication of a material assembly of silver nanoparticles using the phase gradients of optical tweezers,” Phys. Rev. Lett. 114, 143901 (2015). [CrossRef] [PubMed]
14. M. M. Burns, J.-M. Fournier, and J. A. Golovchenko, “Optical matter: Crystallization and binding in intense optical fields,” Science 249, 749–754 (1990). [CrossRef] [PubMed]
15. K. Dholakia and P. Zemánek, “Colloquium: Gripped by light: Optical binding,” Rev. Mod. Phys. 82, 1767–1791 (2010). [CrossRef]
16. A. Jonáš and P. Zemánek, “Light at work: The use of optical forces for particle manipulation, sorting, and analysis,” Electrophoresis 29, 4813–4851 (2008). [CrossRef]
17. Y. Roichman, B. Sun, Y. Roichman, J. Amato-Grill, and D. G. Grier, “Optical forces arising from phase gradients,” Phys. Rev. Lett. 100, 013602 (2008). [CrossRef] [PubMed]
18. J. Ng, Z. F. Lin, C. T. Chan, and P. Sheng, “Photonic clusters formed by dielectric microspheres: numerical simulations,” Phys. Rev. B 72, 085130 (2005). [CrossRef]
19. R. Gordon, M. Kawano, J. T. Blakely, and D. Sinton, “Optohydrodynamic theory of particles in a dual-beam optical trap,” Phys. Rev. B 77, 245125 (2008). [CrossRef]
20. G. Volpe and G. Volpe, “Simulation of a brownian particle in an optical trap,” Am. J. Phys. 81, 224–230 (2013). [CrossRef]
21. N. Wang, J. Chen, S. Liu, and Z. Lin, “Dynamical and phase-diagram study on stable optical pulling force in bessel beams,” Phys. Rev. A 87, 063812 (2013). [CrossRef]
22. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech house, 2005), 3rd ed.
23. R. Biswas and D. R. Hamann, “Simulated annealing of silicon atom clusters in langevin molecular dynamics,” Phys. Rev. B 34, 895–901 (1986). [CrossRef]
24. P. M. Hansen, V. K. Bhatia, N. Harrit, and L. Oddershede, “Expanding the optical trapping range of gold nanoparticles,” Nano Lett. 5, 1937–1942 (2005). [CrossRef] [PubMed]
25. Z. Li, M. Käll, and H. Xu, “Optical forces on interacting plasmonic nanoparticles in a focused gaussian beam,” Phys. Rev. B 77, 085412 (2008). [CrossRef]
26. P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6, 4370–4379 (1972). [CrossRef]
27. S. Gray and T. Kupka, “Propagation of light in metallic nanowire arrays: Finite-difference time-domain studies of silver cylinders,” Phys. Rev. B 68, 045415 (2003). [CrossRef]
28. R. A. Nome, M. J. Guffey, N. F. Scherer, and S. K. Gray, “Plasmonic interactions and optical forces between au bipyramidal nanoparticle dimers,” J. Phys. Chem. A 113, 4408–4415 (2009). [CrossRef] [PubMed]
29. R. Pfeifer, T. Nieminen, N. Heckenberg, and H. Rubinsztein-Dunlop, “Colloquium: Momentum of an electromagnetic wave in dielectric media,” Rev. Mod. Phys. 79, 1197–1216 (2007). [CrossRef]
30. G. Milne, K. Dholakia, D. McGloin, K. Volke-Sepulveda, and P. Zemánek, “Transverse particle dynamics in a bessel beam,” Opt. Express 15, 13972–13987 (2007). [CrossRef] [PubMed]
31. I. Capoglu, A. Taflove, and V. Backman, “Computation of tightly-focused laser beams in the fdtd method,” Opt. Express 21, 87–101 (2013). [CrossRef] [PubMed]
32. J. Ng, Z. Lin, and C. T. Chan, “Theory of optical trapping by an optical vortex beam,” Phys. Rev. Lett. 104, 103601 (2010). [CrossRef] [PubMed]
33. R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University, 2001).
34. B. Leimkuhler and C. Matthews, “Robust and efficient configurational molecular sampling via Langevin dynamics,” J. Chem. Phys. 138, 174102 (2013). [CrossRef] [PubMed]
35. B. Leimkuhler and C. Matthews, “Rational construction of stochastic numerical methods for molecular sampling,” Appl. Math. Res. Express 2013, 34–56 (2013).
36. G. E. P. Box and M. E. Muller, “A note on the generation of random normal deviates,” Ann. Math. Statist. 29, 610–611 (1958). [CrossRef]
37. P. Allen, D. Frenkel, and J. Talbot, Molecular Dynamics Simulation Using Hard Particles (North-Holland, 1989).
38. W. Yu, R. Mittra, T. Su, Y. Liu, and X. Yang, Parallel Finite-Difference Time-Domain Method (Artech House, 2006).
39. F. Borghese, P. Denti, R. Saija, and M. A. Iatì, “Optical trapping of nonspherical particles in the t-matrix formalism,” Opt. Express 15, 11984–11998 (2007). [CrossRef] [PubMed]
40. G. Baffou, R. Quidant, and F. J. G. de Abajo, “Nanoscale control of optical heating in complex plasmonic systems,” ACS Nano 4, 709–716 (2010). PMID: [PubMed] . [CrossRef]
41. Y. Bao, Z. Yan, and N. F. Scherer, “Optical printing of electrodynamically coupled metallic nanoparticle arrays,” J. Phys. Chem. C 118, 19315–19321 (2014). [CrossRef]
42. A. S. Urban, A. A. Lutich, F. D. Stefani, and J. Feldmann, “Laser printing single gold nanoparticles,” Nano Lett. 10, 4794–4798 (2010). [CrossRef] [PubMed]
43. M. J. Guffey and N. F. Scherer, “All-optical patterning of au nanoparticles on surfaces using optical traps,” Nano Lett. 10, 4302–4308 (2010). [CrossRef] [PubMed]
44. Z. Yan, J. E. Jureller, J. Sweet, M. J. Guffey, M. Pelton, and N. F. Scherer, “Three-dimensional optical trapping and manipulation of single silver nanowires,” Nano Lett. 12, 5155–5161 (2012). [CrossRef] [PubMed]
45. H.-J. Butt, K. Graf, and M. Kappl, Physics and Chemistry of Interfaces (Wiley-VCH, 2003). [CrossRef]
46. W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal Dispersions (Cambridge University, 1989). [CrossRef]
47. Z. Yan, R. A. Shah, G. Chado, S. Gray, M. Pelton, and N. F. Scherer, “Guiding spatial arrangements of silver nanoparticles by optical binding interactions in shaped light fields,” ACS Nano 7, 1790–1802 (2013). [CrossRef] [PubMed]
