THE JOURNAL OF CHEMICAL PHYSICS 141, 224104 (2014)
Modeling light-induced charge transfer dynamics across a metal-molecule-metal junction: Bridging classical electrodynamics and quantum dynamics Zixuan Hu, Mark A. Ratner, and Tamar Seidemana) Department of Chemistry, Northwestern University, Evanston, Illinois 60208-3113, USA
(Received 21 May 2014; accepted 19 November 2014; published online 8 December 2014) We develop a numerical approach for simulating light-induced charge transport dynamics across a metal-molecule-metal conductance junction. The finite-difference time-domain method is used to simulate the plasmonic response of the metal structures. The Huygens subgridding technique, as adapted to Lorentz media, is used to bridge the vastly disparate length scales of the plasmonic metal electrodes and the molecular system, maintaining accuracy. The charge and current densities calculated with classical electrodynamics are transformed to an electronic wavefunction, which is then propagated through the molecular linker via the Heisenberg equations of motion. We focus mainly on development of the theory and exemplify our approach by a numerical illustration of a simple system consisting of two silver cylinders bridged by a three-site molecular linker. The electronic subsystem exhibits fascinating light driven dynamics, wherein the charge density oscillates at the driving optical frequency, exhibiting also the natural system timescales, and a resonance phenomenon leads to strong conductance enhancement. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4903046] I. INTRODUCTION
Nanoplasmonics has grown into an exciting frontier of both theoretical and experimental research.1, 2 Rapid advances in both nanofabrication techniques3, 4 and numerical modeling5–7 allow precise synthesis and handling of nanoscaled structures, as well as reliable prediction and interpretation of plasmonic effects, leading to a variety of applications. One example is surface-enhanced Raman spectroscopy (SERS)8, 9 that utilizes the light enhancement effect of plasmonic resonances to greatly increase the ability to detect and distinguish single or a few molecules.10 The light enhancement effect has also been applied to other spectroscopic techniques such as surface-enhanced fluorescence11 and surface-enhanced infrared,12 as well as to the development of plasmon-enhanced solar cells,13 sensors and medical diagnostics,14 to mention but a few of the established examples. Another important area of applications utilizes the light confinement property of various metallic nanoconstructs15–17 for the development of photonic circuits, which have an important role in the further miniaturization of electronic devices.18 Many studies in nanoplasmonics are done with nanostructures that are in close vicinity of each other, constituting a nano-scaled junction.19–23 The majority of applications and theories developed so far have focused on the enhancement and/or guiding of electromagnetic energy,22, 24, 25 but other opportunities are equally inviting. In particular, in the well-established field of molecular electronics,26 charge transfer dynamics across a metal-molecule-metal nano-junction has been studied extensively.27–29 The vast majority of experimental and a) Author to whom correspondence should be addressed. Electronic mail:
[email protected] 0021-9606/2014/141(22)/224104/8/$30.00
numerical studies on nano-junction charge transfer have been carried out in the energy domain, under constant potential biases. The combination of ultrafast (short-pulse induced) nanoplasmonics and charge transfer through a nano-junction could introduce both interesting phenomena and new applications. Furthermore, experimentally short-pulse driven, plasmon-enhanced molecular junctions are currently becoming a reality,30, 31 and may be expected to grow in scope and application. In particular, light-excited plasmons have been demonstrated to alter the photoconductivity of conjugated organic molecules.32 Photon-coupled electron transfer from an STM tip to a molecule has also been observed.33 Although theorists have proposed several imaginative concepts for light-controlled conductance junctions as well as formal approaches for their descriptions,34, 35 the numerical simulation of a combination of time-domain nanoplasmonics with transport via molecular junctions remains a challenge. One difficulty is the proper numerical description of two vastly disparate length scales: the metallic nanoparticles with a typical size-scale of several tens of nanometers and the molecularsized bridge. A second, related challenge is the combination of the classical description used in the vast majority of studies for the plasmonic subsystem and the quantum description that is typically essential for the molecular subsystem. In the present work we develop a numerical method for studying light-triggered, plasmon-enhanced transport via a nanoparticle-molecule-nanoparticle junction, and apply it to directly calculate the coupling between the plasmons and the molecule and thus explore the ensuing electronic dynamics in a simple model system consisting of a linear molecular chain in contact with two silver nanoparticles. The former challenge, namely, bridging the drastic size-scale disparity between the nanojunction and the molecule, is addressed using a recently developed subgridding scheme36 within the
141, 224104-1
© 2014 AIP Publishing LLC
224104-2
Hu, Ratner, and Seideman
J. Chem. Phys. 141, 224104 (2014)
finite-difference time-domain (FDTD)37 method. The subgridding approach36, 38 has been shown to greatly increase the resolution in the critical junction region, allowing accurate and seamless description of both the molecule and the nanoparticles. The electric current density and charge density calculated within classical electrodynamics are transformed into a quantum mechanical electronic wavefunction, which is coupled with the molecular component of the junction within Fermi’s golden rule. Finally the charge dynamics and transport properties are studied with the Heisenberg equations of motion. The details of the theory are discussed in Sec. II. The specific structural details of the model and parameters used in simulations are presented in Sec. III. Section IV presents and discusses the results and the final section summarizes our conclusions and suggests avenues for future research.
gridding provides greater main-sub ratio (greater subgrid resolution), little numerical reflection at the main grid-sub grid boundary, good late time numerical stability, and easy implementation. For details of the Huygens subgridding implementation, see Ref. 36. Note that in our model only the classical contribution to the plasmonic effect is considered, and the nonlinear feedback effect from the molecule to the plasmonic fields as suggested by Ref. 23 is not considered. The strongest excitation field to be used in our model is 0.030 V/nm, significantly lower than the threshold for nonlinear molecular feedback.23 Under field strength of this magnitude, high molecular density is required for the linear molecular feedback to affect the plasmonic fields, as supported by numerical calculations.41, 42 B. Wavefunction generation and electron transfer from metal to molecule
II. THEORY A. FDTD Huygens subgridding
The plasmonic response on the nanojunction can be efficiently and accurately simulated by classical numerical methods.6 The finite-difference time-domain (FDTD) method37, 39 has been widely used due to its accuracy, versatility, and ease of implementation. The time-domain nature of FDTD makes it especially suitable for studying lightinduced dynamics. As do many numerical modeling techniques, FDTD discretizes space and time to propagate the electromagnetic fields on numerical grids. Often in plasmonic studies on nanostructures, a critical region requires a much finer grid than the other regions. In the current study, the region between the nanoparticles requires a fine grid. This calls for an effective subgridding technique that can increase the resolution in a specific region by embedding a subgrid into the main computational grid. One such subgridding algorithm, called the Huygens subgridding, has been developed for Lorentz media in a previous work.36 Huygens subgridding is a subgridding technique for FDTD38, 40 that utilizes the Love’s equivalence principle in electrodynamics to pass fields information across the main grid and the subgrid. Love’s equivalence principle states that the fields inside a closed surface (also called a Huygens surface) are uniquely determined by the boundary fields through the relations ⎧ ∂B ⎪ ⎪ ⎨ ∇ × E = − ∂t − Ks , (1) ⎪ ⎪ ⎩ ∇ × H = ∂D + J s ∂t Js = n × Hi ,
Ks = −n × Ei .
(2)
Equations (1) are the curl Maxwell equations with additional terms Js and Ks as boundary electric and magnetic current densities. Equations (2) show how Js and Ks are calculated from the boundary fields Hi and Ei as produced by electromagnetic sources outside, n being a unit vector normal to the boundary surface. With (1) and (2) we can pass fields information across the boundaries between the main grid and the subgrid, as required by the subgridding algorithm. Compared to other FDTD subgridding algorithms, Huygens sub-
The fields calculated within the sub-gridding approach of Subsection II A are subsequently used to calculate the charge density ρ cl (cl stands for classical since it is calculated by solving the classical Maxwell equations) and electric current density J at the junction interface. The calculation of charge and electric current densities from selected Drude and Lorentz poles (in this study the Drude pole is used since we are interested in free electrons) has been described in a previous work.43 So far we have been using a classical algorithm (FDTD) to calculate classical physical quantities. To propagate electron density through the molecular linker requires that the classical quantities be converted into forms that can be used for quantum mechanical treatments. Here we introduce a single-electron approximation that ignores electron-electron interaction in the metal. The probability density ρ pr associated with the single electron wavefunction ψ is then given within the standard quantum mechanical postulates as the squared modulus of ψ. To relate ρ pr to the classical charge density ρ cl , we introduce the renormalization ρcl
ρpr ≡ V
,
(3)
ρcl dx
where V is the total volume where the wavefunction ψ resides, which is highly localized on the edge of the metal lead. The validity of Eq. (3) is substantiated by first noting that ρ dx ρ pr ≥ 0 and V ρpr dx = V ρcl dx = 1, hence the probability V cl density is well-defined and normalized. Second, we observe that V ρcl dx is the total charge Q of the volume V , hence Eq. (3) can be rewritten as Q · (ρpr dV ) = (ρcl dV ),
(4)
where dV is a volume element. Equation (4) has a clear physical interpretation: the total charge multiplied by the probability of finding an electron in a differential volume element gives the expected charge in that volume. Here we have assumed that all electrons in the volume element have the same ρ pr . Thus, Eq. (3) effectively relates the classical quantity ρ cl to the quantum quantity ρ pr by a probabilistic interpretation.
224104-3
Hu, Ratner, and Seideman
J. Chem. Phys. 141, 224104 (2014)
Similarly we determine the single electron probability density current jpr by using (3) to relate the classical electric current J to jpr , jpr = V
J
.
(5) T =
ρcl dx
Given the single electron charge density and current as determined from the plasmon-mediated optical excitation event, we proceed to determine the amplitude and phase of the wavefunction within quantum mechanics using the relations ⎧ ⎪ ⎨ jpr = ¯ (ψ ∗ ∇ψ − ψ∇ψ ∗ ) 2mi , (6) ⎪ ⎩ ∗ ρpr = ψ ψ where ¯ is Planck’s √ constant divided by 2π , m is the electron mass, and i = −1. Adding the gradient of ρ pr to jpr , we obtain ∇ρpr = ψ ∗ ∇ψ + ψ∇ψ ∗ ,
(7)
2mi (8) j + ∇ρpr = 2ψ ∗ ∇ψ. ¯ pr When ψ is zero, both ρ pr and jpr are trivially zero. When ρ
ψ is not zero, we can write ψpr = ψ ∗ , which, upon substitution into Eq. (8) and several rearrangements yields ∇ρpr mi jpr + · ψ = ∇ψ ¯ ρpr 2ρpr ⎧ ⎪ ∂ψ mi ∂ 1 ⎪ ⎪ = ρ · j + ·ψ ⎪ ⎪ ⎪ ∂x ¯ρpr x ∂x pr 2ρpr ⎪ ⎪ ⎪ ⎪ ⎨ ∂ψ mi ∂ 1 (9) → ·ψ . = j + ρ · ⎪ ∂y ¯ρ y ∂y pr 2ρpr ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ mi ∂ψ ∂ 1 ⎪ ⎪ ·ψ ⎪ ⎩ ∂z = ¯ρ jz + ∂z ρpr · 2ρ pr Equation (9) is a system of partial differential equations that are easy to solve due to their separability. If we write the solution in the form ψ = Aei(κx x+κy y+κz z) , then ψ = Aei(κx x+κy y+κz z) , ρcl
, A = ρpr =
ρcl dx
(10)
V
κν =
mjν mJ = ν ¯ρpr ¯ρcl
With the single electron wave function ψdetermined, we proceed by using the Fermi golden rule to calculate the transition rate from the metal to the molecule via a single atomic orbital,
(ν = x, y, z).
Note that although the form of the solution resembles a planewave, it is actually not a planewave because the κ ν depend on Jν and ρ cl , which in turn are functions of space and time – although the κ ν ‘s give information about the space variation of the wavefunction, they are not the wavevector components. With Eqs. (3)–(10), we have converted ρ cl and J, both available from the classical calculation, into ψ, a quantum mechanical wavefunction.
2π | ψ| W |φ |2 δ( E), ¯
(11)
where T is the probability of transition per unit time, ψ is calculated in (10), φ is the atomic orbital nearest the metal contact, W is the interaction potential, δ is the Dirac delta function, and E is the energy difference between the metal and the adjacent atomic state eigenenergies. The transition rate Tin Eq. (11) applies to a single state in the metal. If we make the conventional assumption that all the states in the metal interact in the same way with the molecule,44 then the total transition rate summed over all available states is R=
2π | ψ| W |φ |2 ρDoS . ¯
(12)
where ρ DoS is the density of states that are available for the transition – here they are the excess charge states excited by the light in the volume V adjacent to the molecule, whose number of states n is given by 1 Q = · ρ dx, (13) n= −e −e V cl where e is the elementary charge and V is a small volume near the edge of the metal, defined by the region where the metal wavefunction overlaps with the nearest atomic orbital n and substituting (13) of the molecule. Using ρDoS = unit energy into (12) we write the desired expression for the rate of electron transfer from the metal to the nearest atomic site on the molecule as ρcl dx 2π V | ψ| W |φ |2 . · (14) R= −e · (unit energy) ¯ From Eq. (14) we can see that the plasmonic effect on the charge transfer rate enters through the integral term V ρcl dx and the wavefunction ψ. Since ψ is normalized, the mag nitude of the rate is mostly determined by V ρcl dx, which grows linearly with the excitation light source. Above we discussed the situation when the metal edge is negatively charged. We now introduce the assumption that when the metal edge is positively charged, a total charge −Q in a small region near the edge of the metal has exactly the opposite effect to a total charge Q. That is, if Q results in electron deposition (from metal to molecule) at the rate R, then −Q results in hole injection (from metal to molecule) at the same rate R.
C. Electron propagation through molecular sites
The transition rate R calculated in the last section serves as an electron source for the electron transport via the molecule. To illustrate how electron transport is calculated through the molecule, we consider a model system consisting of N identical atomic sites as an example. In second quantiza-
224104-4
Hu, Ratner, and Seideman
J. Chem. Phys. 141, 224104 (2014)
tion the Hamiltonian of this system can be written as H =α
N
†
ak ak + β
k=1
N−1
†
†
(ak ak+1 + ak+1 ak ),
(15)
k=1
where α is the site energy, β is the hopping integral between † adjacent sites, ak is the creation operator, and ak is the annihilation operator. Using the Heisenberg equations of motion for operators, we have
Initial values of the various operators in Eqs. (16a) and (16b) at t = 0 are all set to zero. As expected, it is readily shown that α in Eqs. (16a) and (16b) cancels out (as the zero of energy is immaterial), hence its actual value is irrelevant to the results. The β parameter is taken, reasonable for a qualitative analysis, to be −3 eV (0.11 hartree).45 In Sec. IV we solve Eqs. (16a) and (16b) for a simple 3-site model (N = 3) to illustrate the charge transport dynamics on the atomic sites. III. COMPUTATIONAL DETAILS
i d † † a a = [H, a1 a1 ] + S, dt 1 1 ¯ d † i † a a = [H, aj aj ], dt j j ¯
j = 2, . . . , N − 1 ,
(16a)
d † i † aN aN = [H, aN aN ] − D, dt ¯ where S, the electron source, is added to the first site to account for electrons entering the molecular sites and D, the electron drain, is subtracted from the last site to account for electrons leaving the molecular sites. We note that the nomenclature “source” and “drain” refer to the overall flow of electrons from the excited metal cylinder to Site 1, through the molecular chain to Site N, and then to the metal cylinder on the other end. At an arbitrary time during the fluctuation electrons might leave through the electron source or enter through the electron drain. When S is positive, the electrons are added from the metal lead to the 1st molecular site at the rate R as calculated in the last section; when S is negative, the electrons are drained from the 1st molecular site to the nearest metal lead, at a rate proportional to the electron population on † the 1st site, i.e., R · a1 a1 . Similar dynamics apply to the other side of the molecular chain: when D is negative, the electrons are added from the metal to the Nth molecular site, at the rate D = R; when D is positive, the electrons are drained from the Nth site to the nearest metal lead, at a rate proportional to the † electron population on the Nth site, i.e., D = R · aN aN . Equation (16a) is supplemented by analogous relations for the off-diagonal elements i d † † a a = [H, aj aj ], dt j j ¯
j = j .
(16b)
The FDTD-subgridding computational setup used in this study is shown in Figure 1. The two silver cylinders are equally sized at 7 nm radius and the distance between them is 1 nm. IS and OS stand for the inner and outer surfaces, respectively, used in the Huygens subgridding approach to pass fields information across the main-grid/subgrid boundary.36 The main grid spatial size is 0.1 nm and the main-grid/subgrid ratio is 9, making the subgrid step size approximately 0.011 nm. An infinite line light source (a point source in the 2D simulation setup) is placed 5 nm below the lower cylinder. For all simulations in this work, the light source is a Gaussian pulse with 10 fs FWHM and 400 nm wavelength. The electric field of the source is polarized along the long axis of the silver dimer. The material parameters for silver are obtained from Ref. 46. The charge and electric current densities calculated by FDTD-subgridding are then used to calculate the transition rate R as shown in Eq. (14). In this study we use a simple molecular model of three silver atoms equally spaced between the silver cylinders, as shown in the focus view in Figure 1. The silver valence 5s orbitals on the first and third atoms are used as the atomic wavefunctions φ in (14). The transition rate is subsequently used as the electron source and drain in (16a). The Heisenberg equations of motion (16a) and (16b) are then solved numerically with the commercial software Maple 12.
IV. RESULTS AND DISCUSSIONS A. Qualitative study of the charge transfer dynamics
The qualitative behavior of the charge transfer dynamics is determined by the electron source S and the electron drain D
FIG. 1. The computational setup, where IS and OS denote the inner and outer surfaces used in Huygens subgridding. On the right side is the focus view of the gap region, with the spacing of the atoms shown in between the gap.
224104-5
Hu, Ratner, and Seideman
FIG. 2. Electron dynamics on three atomic sites with R = 0.004 · sin (0.11t). † † Red curve corresponds to a1 a1 , green curve corresponds to a2 a2 , and blue † curve corresponds to a3 a3 .
in (16a). To gain a basic understanding of the dynamic system, we first study several simple examples for the electron source and electron drain. In the case where the two nanoparticles are of the same material and size and the incident light is placed so as to excite both symmetrically, the electron source and electron drain are symmetric, that is, the coupling strength R is the same for both S and D. If the nanoparticles are excited by a single frequency harmonic field, we have R = Asin (ω · t), where A is an amplitude determined by the plasmonic effect through Eq. (14), and ω coincides with the angular frequency of the incident light. We begin our survey with a weak source, A = 0.004 atomic units, and an incident frequency equal to the plasmon resonance – approximately 0.11 atomic units, corresponding to a 400 nm wavelength. With these parameters substituted in Eqs. (16a) and (16b) for N = 3 (see Figure 1), solution of the differential equations system yields the results shown in Figure 2. Here we observe electron population starting to build up on all atomic sites upon excitation while also exhibiting distinct oscillations. A Fourier transform of the dynamics in † † Figure 2 shows that the oscillations of a1 a1 and a3 a3 contain two frequencies, one at the source frequency as expected, and another one at ω = 0.157 a.u., which suggests the existence of a natural frequency for the dynamic system. This is confirmed by setting R = 0.004 · sin (0.157t), with which the dynamics is shown in Figure 3. Here we observe a striking resonance effect, where the oscillations at all three atomic sites are greatly enhanced, † and the unity population cap is quickly reached. Both a1 a1
† and a3 a3 exhibit perfect harmonic behavior at ω1 = ω3 † = 0.157 a.u., while a2 a2 oscillates at double that value at ω2 † † = 0.314 a.u. This is because a1 a1 and a3 a3 are perfectly † out of phase, hence their collective effect on a2 a2 is an os† cillation of frequency 2ω = ω2 , by which a2 a2 is readily driven. Additional calculations show that these natural system frequencies are proportional to β, the hopping integral term in the tight-binding Hamiltonian in Eq. (15). This is expected since β determines the tendency for the molecular sites to exchange electrons. The above results point to an interest-
J. Chem. Phys. 141, 224104 (2014)
FIG. 3. Electron dynamics on three atomic sites with R = 0.004 · † sin (0.157t). Red curve corresponds to a1 a1 , green curve corresponds to † † a2 a2 , and blue curve corresponds to a3 a3 .
ing control opportunity with potentially useful applications. By setting the excitation frequency ω equal to the natural system frequency of 0.157 atomic units, we generate a resonance phenomenon wherein the charge oscillation amplitudes are enhanced by close to an order of magnitude as compared to their values in Fig. 2. In addition, the electron populations grow much faster. It is important to note that the effect of the excitation frequency on the amplitude of charge oscillation is twofold: the plasmonic resonance frequency is determined by the dielectric function and geometry of the nanostructure and the molecular transfer resonance frequency is determined by the Hamiltonian of the system described in Eq. (15). Hence in experiments the ideal amplification in Figure 4 may not be realizable, requiring that a balance between the two resonance frequencies would be maintained for the charge transfer amplitude to be maximized. This result will be explored in a future study, where the focus will be on the optical control of the charge transfer dynamics. For the purpose of this discussion we simply use ω = 0.11 atomic units for maximized plasmonic resonance. For ω = 0.11, we observe the steady-state behavior without extended time by increasing the amplitude to A = 0.04.
FIG. 4. Electron dynamics on three atomic sites with R = 0.04 · sin (0.11t). † † Red curve corresponds to a1 a1 , green curve corresponds to a2 a2 , and blue † curve corresponds to a3 a3 .
224104-6
Hu, Ratner, and Seideman
FIG. 5. Electron dynamics on three atomic sites with S = R = 0.02 and D † = 2R = 0.04. Red curve corresponds to a1 a1 , green curve corresponds to † † a2 a2 , and blue curve corresponds to a3 a3 .
The dynamics is shown in Figure 4, where we see that the † † populations a1 a1 and a3 a3 quickly grow to saturation at unity and continue to oscillate in a symmetric but out of phase † fashion. The bridging site population a2 a2 oscillates at a different frequency and with a smaller amplitude. There are two major differences between the dynamics in Figure 4 and that in Figure 3: (1) all oscillation amplitudes are significantly smaller in Figure 4, which is reasonable because the excitation is off-resonance with the natural system time scale; (2) both † † a1 a1 and a3 a3 oscillate at the source frequency, not their † natural frequencies, whereas a2 a2 oscillates at double the source frequency. This is also expected because the steadystate of a driven oscillator should have the frequency of the driving force – in this case the excitation source. To conclude our qualitative study, we test the reliability of the proposed method by examining the model results in the well studied DC limit. To avoid overwhelming the molecular chain with electrons we intentionally set the source weaker than the drain. (This problem does not apply to the AC source as used in the rest of Secs. IV A and IV B, because in that case both the source and the drain rapidly oscillate and thus do not cause electrons to overcrowd on the molecular chain.) Here we set S = R = 0.02 and D = 2R = 0.04, leading to the DC dynamics shown in Figure 5, where we find that under a DC source, equivalent to a constant potential bias, each of the three sites quickly reaches a constant steady-state. The steady-state populations are all low because the drain is efficient and the source is not sufficiently strong to overwhelm the electron transfer between sites. This dynamics is expected for a constant potential bias, to which our model reduces in the limit of a DC source.
B. Quantitative rate of charge transfer
Subsection IV A has examined in detail the qualitative charge transfer dynamics governed by the Heisenberg equations of motion. In this subsection we combine the Heisenberg equations with the FDTD subgridding method of Ref. 36 to gain a quantitative appreciation of the magnitude of
J. Chem. Phys. 141, 224104 (2014)
FIG. 6. Electron dynamics on three atomic sites with R = 0.04 · sin (0.11t) † and D = 0.76S. Red curve corresponds to a1 a1 , green curve corresponds to † † a2 a2 , and blue curve corresponds to a3 a3 .
the charge transfer dynamics. As in the calculations discussed above, we take the incident frequency to equal the plasmon resonance of the silver double cylinder structure considered, ω = 0.11 atomic units, corresponding to 400 nm wavelength. At this wavelength we carry out the FDTD-subgridding simulation according to the parameters described in Sec. III. For an incident field strength of 0.030 V/nm (or translated into laser intensity ≈2.4 × 108 W/cm2 ), Eqs. (10) and (14) give R = 0.04 · sin (0.11t)atomic units. In the discussions of Figures 2–5 we set the electron drain D to have exactly the same amplitude as R so as to begin with a simple symmetric model that merely extracts the basic dynamics. In the actual simulation setup we place the light source behind one of the silver cylinders (Figure 1) to model a realistic lightdriven charge transport experiment. The result is a smaller amplitude (by about 24%) for the drain D compared to the source S, creating a similar but different dynamics from Figure 4. The electron population evolution obtained is shown in Figure 6. The charge transfer dynamics in Figure 6 is similar to that † † in Figure 4, with two major differences: (1) a1 a1 and a3 a3
now have slightly different amplitudes due to the asymmetry † between the source and the drain; (2) the a2 a2 curve can be seen to oscillate according to several frequencies, rather than the single frequency exhibited in Figure 4. A Fourier † transform of a2 a2 confirms this observation, as illustrated in Figure 7, which shows that 2 additional frequencies in † the Fourier transform of a2 a2 arise with the asymmetric simulation setup. The peak at 0.22 a.u. corresponds to the natural system time-scale as in Figure 4. The peak at 0.11 a.u. – the driving frequency, which here is taken to be the plasmon resonance frequency – illustrates that the asymmetry of the source and the drain gives rise to effect of the incident frequency on the population at all sites. Finally, the peak at 0.33 a.u. is an overtone corresponding to the sum of the incident frequency at 0.11 a.u. and the natural system frequency at 0.22 a.u. The dominating peak is still at 0.22 a.u., hence Figure 6 is not significantly different from Figure 4, but it shows that the dynamics is richer with a realistic simulation setup.
224104-7
Hu, Ratner, and Seideman
J. Chem. Phys. 141, 224104 (2014)
†
FIG. 7. Fourier transform for a2 a2 data in Figure 6. The DC component has been truncated for better scaling.
V. CONCLUSION
In this study we have proposed a viable approach to simulate light-induced charge transfer dynamics across a metalmolecule-metal junction, with a novel scheme to calculate the coupling between the metal and the molecule. To that end we combined several classical and quantum mechanical methods: (1) use of a newly developed Huygens subgridding approach36 to efficiently handle the vastly disparate size-scales of the gap and the two nanoparticles; (2) transformation of the current and charge densities computed by solution of the Maxwell equations into a quantum mechanical wavefunction; (3) use of the Fermi golden rule to calculate the rate of electron transfer from the plasmonic excitation to the molecular orbital; and (4) use of the Heisenberg equations of motion to propagate the electron population through the molecular bridge. For the purpose of this paper we have focused on clarifying the concepts and the methods, using to that end a highly simplified (although very common) model molecule. The approach is general, however, and we expect it to prove capable of simulating a wide variety of light-driven and light-controlled plasmon-mediated phenomena. These include junctions with realistic (ab initio) molecules, unidirectional current in asymmetric junctions, transport in STM tip-molecule-substrate environments, transport between semiconductor electrodes, and transport controlled by the incident excitation, e.g., pulse trains, shaped pulses, and two pathway excitations.
ACKNOWLEDGMENTS
T.S. is grateful to the National Science Foundation (Award No. CHEM-1012207) and the Department of Energy (Award No. DE-FG02-09ER16109) for support. M.A.R. thanks the National Science Foundation for support (CHE 1058896). 1 N.
J. Halas, Nano Lett. 10, 3816–3822 (2010). I. Stockman, Opt. Express 19, 22029 (2011). 3 J. Henzie, J. Lee, M. H. Lee, W. Hasan, and T. W. Odom, Annu. Rev. Phys. Chem. 60, 147–165 (2009). 4 H. Duan, A. I. Fernandez-Dominguez, M. Bosman, S. A. Maier, and J. K. W. Yang, Nano Lett. 12, 1683–1689 (2012). 5 F. L. Teixeira, IEEE Trans. Antennas Propag. 56(8), 2150 (2008). 2 M.
6 J.
Zhao, A. O. Pinchuk, J. M. McMahon, S. Li, L. K. Ausman, A. L. Atkinson, and G. C. Schatz, Acc. Chem. Res. 41(12), 1710 (2008). 7 R. Esteban, A. G. Borisov, P. Nordlander, and J. Aizpurua, Nat. Commun. 3, 825 (2012). 8 K. A. Willets and R. P. V. Duyne, Annu. Rev. Phys. Chem. 58, 267–297 (2007). 9 J. M. McMahon, A.-I. Henry, K. L. Wustholz, M. J. Natan, R. G. Freeman, R. P. V. Duyne, and G. C. Schatz, Annl. Bioanal. Chem. 394, 1819–1825 (2009). 10 S. Nie and S. R. Emory, Science 275, 1102–1106 (1997). 11 E. Fort and S. Gresillon, J. Phys. D: Appl. Phys. 41, 013001 (2008). 12 A. Pucci, F. Neubrech, D. Weber, S. Hong, T. Toury, and M. L. d. l. Chapelle, Phys. Status Solidi B 247, 2071–2074 (2010). 13 S. Pillai, K. R. Catchpole, T. Trupke, and M. A. Green, J. Appl. Phys. 101, 093105 (2007). 14 T. Vo-Dinh, A. M. Fales, G. D. Griffin, C. G. Khoury, Y. Liu, H. Ngo, S. J. Norton, J. K. Register, H.-N. Wang, and H. Yuan, Nanoscale 5, 10127 (2013). 15 S. A. Maier and H. A. Atwater, J. Appl. Phys. 98, 011101 (2005). 16 P. Berini, R. Charbonneau, N. Lahoud, and G. Mattiussi, J. Appl. Phys. 98, 043109 (2005). 17 S. A. Maier, P. E. Barclay, T. J. Johnson, M. D. Friedman, and O. Painter, Appl. Phys. Lett. 86, 071103 (2005). 18 E. Ozbay, Science 311(5758), 189–193 (2006). 19 S. K. Gray and T. Kupka, Phys. Rev. B 68, 045415 (2003). 20 J. P. Kottmann and O. J. F. Martin, Opt. Express 8(12), 655–663 (2001). 21 D. C. Marinica, A. K. Kazansky, P. Nordlander, J. Aizpurua, and A. G. Borisov, Nano Lett. 12(3), 1333–1339 (2012). 22 E. Hao and G. C. Schatz, J. Chem. Phys. 120, 357 (2004). 23 K. Lopata and D. Neuhauser, J. Chem. Phys. 130, 104707 (2009). 24 A. J. Haes, C. L. Haynes, A. D. McFarland, G. C. Schatz, R. P. V. Duyne, and S. Zou, MRS Bull. 30, 368 (2005). 25 M. Sukharev and T. Seideman, J. Chem. Phys. 124(14), 144707 (2006). 26 J. C. Cuevas and E. Scheer, Molecular Electronics (World Scientific Publishing Company, 2010). 27 S. Datta, Quantum Transport: Atom to Transistor (Cambridge University Press, 2005). 28 M. D. Ventra, Electrical Transport in Nanoscale Systems (Cambridge University Press, 2008). 29 A. Nitzan and M. A. Ratner, Science 300, 1384 (2003). 30 J. M. Klingsporn, M. D. Sonntag, T. Seideman, and R. P. V. Duyne, J. Phys. Chem. Lett. 5, 106–110 (2014). 31 P. Dombi, A. Horl, P. Racz, I. Marton, A. Trugler, J. R. Krenn, and U. Hohenester, Nano Lett. 13, 674–678 (2013). 32 P. Banerjee, D. Conklin, S. Nanayakkara, T.-H. Park, M. J. Therien, and D. A. Bonnell, ACS Nano 4, 1019–1025 (2010). 33 S. W. Wu, N. Ogawa, and W. Ho, Science 312, 1362–1365 (2006). 34 B. D. Fainberg, M. Jouravlev, and A. Nitzan, Phys. Rev. B 76, 245329 (2007). 35 B. D. Fainberg, M. Sukharev, T.-H. Park, and M. Galperin, Phys. Rev. B 83, 205425 (2011). 36 Z. Hu, M. A. Ratner, and T. Seideman, J. Chem. Phys. 137(20), 204111 (2012).
224104-8
Hu, Ratner, and Seideman
37 A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-
Difference Time-Domain Method, 3rd ed. (Artech House, Boston, 2005). 38 J.-P. Berenger, J. Comp. Phys. 230, 5635 (2011). 39 K. S. Yee, IEEE Trans. Antennas Propag. 14, 302–307 (1966). 40 J.-P. Berenger, IEEE Trans. Antennas Propag. 57, 3860–3867 (2009). 41 A. Salomon, R. J. Gordon, Y. Prior, T. Seideman, and M. Sukharev, Phys. Rev. Lett. 109, 073002 (2012).
J. Chem. Phys. 141, 224104 (2014) 42 M.
Sukharev, T. Seideman, R. J. Gordon, A. Salomon, and Y. Prior, ACS Nano 8(1), 807–817 (2014). 43 Z. Hu, M. A. Ratner, and T. Seideman, J. Chem. Phys. 136, 054504 (2012). 44 V. Mujica, M. Kemp, and M. A. Ratner, J. Chem. Phys. 101, 6849 (1994). 45 Y. Xie and J. A. Blackman, Phys. Rev. B 64, 195115 (2001). 46 A. D. Rakic, A. B. Djurisic, J. M. Elazar, and M. L. Majewski, Appl. Opt. 37(22), 5271–5283 (1998).
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
Modeling light-induced charge transfer dynamics across a metal-molecule-metal junction: Bridging classical electrodynamics and quantum dynamics Zixuan Hu, Mark A. Ratner, and Tamar Seidemana) Department of Chemistry, Northwestern University, Evanston, Illinois 60208-3113, USA
(Received 21 May 2014; accepted 19 November 2014; published online 8 December 2014) We develop a numerical approach for simulating light-induced charge transport dynamics across a metal-molecule-metal conductance junction. The finite-difference time-domain method is used to simulate the plasmonic response of the metal structures. The Huygens subgridding technique, as adapted to Lorentz media, is used to bridge the vastly disparate length scales of the plasmonic metal electrodes and the molecular system, maintaining accuracy. The charge and current densities calculated with classical electrodynamics are transformed to an electronic wavefunction, which is then propagated through the molecular linker via the Heisenberg equations of motion. We focus mainly on development of the theory and exemplify our approach by a numerical illustration of a simple system consisting of two silver cylinders bridged by a three-site molecular linker. The electronic subsystem exhibits fascinating light driven dynamics, wherein the charge density oscillates at the driving optical frequency, exhibiting also the natural system timescales, and a resonance phenomenon leads to strong conductance enhancement. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4903046] I. INTRODUCTION
Nanoplasmonics has grown into an exciting frontier of both theoretical and experimental research.1, 2 Rapid advances in both nanofabrication techniques3, 4 and numerical modeling5–7 allow precise synthesis and handling of nanoscaled structures, as well as reliable prediction and interpretation of plasmonic effects, leading to a variety of applications. One example is surface-enhanced Raman spectroscopy (SERS)8, 9 that utilizes the light enhancement effect of plasmonic resonances to greatly increase the ability to detect and distinguish single or a few molecules.10 The light enhancement effect has also been applied to other spectroscopic techniques such as surface-enhanced fluorescence11 and surface-enhanced infrared,12 as well as to the development of plasmon-enhanced solar cells,13 sensors and medical diagnostics,14 to mention but a few of the established examples. Another important area of applications utilizes the light confinement property of various metallic nanoconstructs15–17 for the development of photonic circuits, which have an important role in the further miniaturization of electronic devices.18 Many studies in nanoplasmonics are done with nanostructures that are in close vicinity of each other, constituting a nano-scaled junction.19–23 The majority of applications and theories developed so far have focused on the enhancement and/or guiding of electromagnetic energy,22, 24, 25 but other opportunities are equally inviting. In particular, in the well-established field of molecular electronics,26 charge transfer dynamics across a metal-molecule-metal nano-junction has been studied extensively.27–29 The vast majority of experimental and a) Author to whom correspondence should be addressed. Electronic mail:
[email protected] 0021-9606/2014/141(22)/224104/8/$30.00
numerical studies on nano-junction charge transfer have been carried out in the energy domain, under constant potential biases. The combination of ultrafast (short-pulse induced) nanoplasmonics and charge transfer through a nano-junction could introduce both interesting phenomena and new applications. Furthermore, experimentally short-pulse driven, plasmon-enhanced molecular junctions are currently becoming a reality,30, 31 and may be expected to grow in scope and application. In particular, light-excited plasmons have been demonstrated to alter the photoconductivity of conjugated organic molecules.32 Photon-coupled electron transfer from an STM tip to a molecule has also been observed.33 Although theorists have proposed several imaginative concepts for light-controlled conductance junctions as well as formal approaches for their descriptions,34, 35 the numerical simulation of a combination of time-domain nanoplasmonics with transport via molecular junctions remains a challenge. One difficulty is the proper numerical description of two vastly disparate length scales: the metallic nanoparticles with a typical size-scale of several tens of nanometers and the molecularsized bridge. A second, related challenge is the combination of the classical description used in the vast majority of studies for the plasmonic subsystem and the quantum description that is typically essential for the molecular subsystem. In the present work we develop a numerical method for studying light-triggered, plasmon-enhanced transport via a nanoparticle-molecule-nanoparticle junction, and apply it to directly calculate the coupling between the plasmons and the molecule and thus explore the ensuing electronic dynamics in a simple model system consisting of a linear molecular chain in contact with two silver nanoparticles. The former challenge, namely, bridging the drastic size-scale disparity between the nanojunction and the molecule, is addressed using a recently developed subgridding scheme36 within the
141, 224104-1
© 2014 AIP Publishing LLC
224104-2
Hu, Ratner, and Seideman
J. Chem. Phys. 141, 224104 (2014)
finite-difference time-domain (FDTD)37 method. The subgridding approach36, 38 has been shown to greatly increase the resolution in the critical junction region, allowing accurate and seamless description of both the molecule and the nanoparticles. The electric current density and charge density calculated within classical electrodynamics are transformed into a quantum mechanical electronic wavefunction, which is coupled with the molecular component of the junction within Fermi’s golden rule. Finally the charge dynamics and transport properties are studied with the Heisenberg equations of motion. The details of the theory are discussed in Sec. II. The specific structural details of the model and parameters used in simulations are presented in Sec. III. Section IV presents and discusses the results and the final section summarizes our conclusions and suggests avenues for future research.
gridding provides greater main-sub ratio (greater subgrid resolution), little numerical reflection at the main grid-sub grid boundary, good late time numerical stability, and easy implementation. For details of the Huygens subgridding implementation, see Ref. 36. Note that in our model only the classical contribution to the plasmonic effect is considered, and the nonlinear feedback effect from the molecule to the plasmonic fields as suggested by Ref. 23 is not considered. The strongest excitation field to be used in our model is 0.030 V/nm, significantly lower than the threshold for nonlinear molecular feedback.23 Under field strength of this magnitude, high molecular density is required for the linear molecular feedback to affect the plasmonic fields, as supported by numerical calculations.41, 42 B. Wavefunction generation and electron transfer from metal to molecule
II. THEORY A. FDTD Huygens subgridding
The plasmonic response on the nanojunction can be efficiently and accurately simulated by classical numerical methods.6 The finite-difference time-domain (FDTD) method37, 39 has been widely used due to its accuracy, versatility, and ease of implementation. The time-domain nature of FDTD makes it especially suitable for studying lightinduced dynamics. As do many numerical modeling techniques, FDTD discretizes space and time to propagate the electromagnetic fields on numerical grids. Often in plasmonic studies on nanostructures, a critical region requires a much finer grid than the other regions. In the current study, the region between the nanoparticles requires a fine grid. This calls for an effective subgridding technique that can increase the resolution in a specific region by embedding a subgrid into the main computational grid. One such subgridding algorithm, called the Huygens subgridding, has been developed for Lorentz media in a previous work.36 Huygens subgridding is a subgridding technique for FDTD38, 40 that utilizes the Love’s equivalence principle in electrodynamics to pass fields information across the main grid and the subgrid. Love’s equivalence principle states that the fields inside a closed surface (also called a Huygens surface) are uniquely determined by the boundary fields through the relations ⎧ ∂B ⎪ ⎪ ⎨ ∇ × E = − ∂t − Ks , (1) ⎪ ⎪ ⎩ ∇ × H = ∂D + J s ∂t Js = n × Hi ,
Ks = −n × Ei .
(2)
Equations (1) are the curl Maxwell equations with additional terms Js and Ks as boundary electric and magnetic current densities. Equations (2) show how Js and Ks are calculated from the boundary fields Hi and Ei as produced by electromagnetic sources outside, n being a unit vector normal to the boundary surface. With (1) and (2) we can pass fields information across the boundaries between the main grid and the subgrid, as required by the subgridding algorithm. Compared to other FDTD subgridding algorithms, Huygens sub-
The fields calculated within the sub-gridding approach of Subsection II A are subsequently used to calculate the charge density ρ cl (cl stands for classical since it is calculated by solving the classical Maxwell equations) and electric current density J at the junction interface. The calculation of charge and electric current densities from selected Drude and Lorentz poles (in this study the Drude pole is used since we are interested in free electrons) has been described in a previous work.43 So far we have been using a classical algorithm (FDTD) to calculate classical physical quantities. To propagate electron density through the molecular linker requires that the classical quantities be converted into forms that can be used for quantum mechanical treatments. Here we introduce a single-electron approximation that ignores electron-electron interaction in the metal. The probability density ρ pr associated with the single electron wavefunction ψ is then given within the standard quantum mechanical postulates as the squared modulus of ψ. To relate ρ pr to the classical charge density ρ cl , we introduce the renormalization ρcl
ρpr ≡ V
,
(3)
ρcl dx
where V is the total volume where the wavefunction ψ resides, which is highly localized on the edge of the metal lead. The validity of Eq. (3) is substantiated by first noting that ρ dx ρ pr ≥ 0 and V ρpr dx = V ρcl dx = 1, hence the probability V cl density is well-defined and normalized. Second, we observe that V ρcl dx is the total charge Q of the volume V , hence Eq. (3) can be rewritten as Q · (ρpr dV ) = (ρcl dV ),
(4)
where dV is a volume element. Equation (4) has a clear physical interpretation: the total charge multiplied by the probability of finding an electron in a differential volume element gives the expected charge in that volume. Here we have assumed that all electrons in the volume element have the same ρ pr . Thus, Eq. (3) effectively relates the classical quantity ρ cl to the quantum quantity ρ pr by a probabilistic interpretation.
224104-3
Hu, Ratner, and Seideman
J. Chem. Phys. 141, 224104 (2014)
Similarly we determine the single electron probability density current jpr by using (3) to relate the classical electric current J to jpr , jpr = V
J
.
(5) T =
ρcl dx
Given the single electron charge density and current as determined from the plasmon-mediated optical excitation event, we proceed to determine the amplitude and phase of the wavefunction within quantum mechanics using the relations ⎧ ⎪ ⎨ jpr = ¯ (ψ ∗ ∇ψ − ψ∇ψ ∗ ) 2mi , (6) ⎪ ⎩ ∗ ρpr = ψ ψ where ¯ is Planck’s √ constant divided by 2π , m is the electron mass, and i = −1. Adding the gradient of ρ pr to jpr , we obtain ∇ρpr = ψ ∗ ∇ψ + ψ∇ψ ∗ ,
(7)
2mi (8) j + ∇ρpr = 2ψ ∗ ∇ψ. ¯ pr When ψ is zero, both ρ pr and jpr are trivially zero. When ρ
ψ is not zero, we can write ψpr = ψ ∗ , which, upon substitution into Eq. (8) and several rearrangements yields ∇ρpr mi jpr + · ψ = ∇ψ ¯ ρpr 2ρpr ⎧ ⎪ ∂ψ mi ∂ 1 ⎪ ⎪ = ρ · j + ·ψ ⎪ ⎪ ⎪ ∂x ¯ρpr x ∂x pr 2ρpr ⎪ ⎪ ⎪ ⎪ ⎨ ∂ψ mi ∂ 1 (9) → ·ψ . = j + ρ · ⎪ ∂y ¯ρ y ∂y pr 2ρpr ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ mi ∂ψ ∂ 1 ⎪ ⎪ ·ψ ⎪ ⎩ ∂z = ¯ρ jz + ∂z ρpr · 2ρ pr Equation (9) is a system of partial differential equations that are easy to solve due to their separability. If we write the solution in the form ψ = Aei(κx x+κy y+κz z) , then ψ = Aei(κx x+κy y+κz z) , ρcl
, A = ρpr =
ρcl dx
(10)
V
κν =
mjν mJ = ν ¯ρpr ¯ρcl
With the single electron wave function ψdetermined, we proceed by using the Fermi golden rule to calculate the transition rate from the metal to the molecule via a single atomic orbital,
(ν = x, y, z).
Note that although the form of the solution resembles a planewave, it is actually not a planewave because the κ ν depend on Jν and ρ cl , which in turn are functions of space and time – although the κ ν ‘s give information about the space variation of the wavefunction, they are not the wavevector components. With Eqs. (3)–(10), we have converted ρ cl and J, both available from the classical calculation, into ψ, a quantum mechanical wavefunction.
2π | ψ| W |φ |2 δ( E), ¯
(11)
where T is the probability of transition per unit time, ψ is calculated in (10), φ is the atomic orbital nearest the metal contact, W is the interaction potential, δ is the Dirac delta function, and E is the energy difference between the metal and the adjacent atomic state eigenenergies. The transition rate Tin Eq. (11) applies to a single state in the metal. If we make the conventional assumption that all the states in the metal interact in the same way with the molecule,44 then the total transition rate summed over all available states is R=
2π | ψ| W |φ |2 ρDoS . ¯
(12)
where ρ DoS is the density of states that are available for the transition – here they are the excess charge states excited by the light in the volume V adjacent to the molecule, whose number of states n is given by 1 Q = · ρ dx, (13) n= −e −e V cl where e is the elementary charge and V is a small volume near the edge of the metal, defined by the region where the metal wavefunction overlaps with the nearest atomic orbital n and substituting (13) of the molecule. Using ρDoS = unit energy into (12) we write the desired expression for the rate of electron transfer from the metal to the nearest atomic site on the molecule as ρcl dx 2π V | ψ| W |φ |2 . · (14) R= −e · (unit energy) ¯ From Eq. (14) we can see that the plasmonic effect on the charge transfer rate enters through the integral term V ρcl dx and the wavefunction ψ. Since ψ is normalized, the mag nitude of the rate is mostly determined by V ρcl dx, which grows linearly with the excitation light source. Above we discussed the situation when the metal edge is negatively charged. We now introduce the assumption that when the metal edge is positively charged, a total charge −Q in a small region near the edge of the metal has exactly the opposite effect to a total charge Q. That is, if Q results in electron deposition (from metal to molecule) at the rate R, then −Q results in hole injection (from metal to molecule) at the same rate R.
C. Electron propagation through molecular sites
The transition rate R calculated in the last section serves as an electron source for the electron transport via the molecule. To illustrate how electron transport is calculated through the molecule, we consider a model system consisting of N identical atomic sites as an example. In second quantiza-
224104-4
Hu, Ratner, and Seideman
J. Chem. Phys. 141, 224104 (2014)
tion the Hamiltonian of this system can be written as H =α
N
†
ak ak + β
k=1
N−1
†
†
(ak ak+1 + ak+1 ak ),
(15)
k=1
where α is the site energy, β is the hopping integral between † adjacent sites, ak is the creation operator, and ak is the annihilation operator. Using the Heisenberg equations of motion for operators, we have
Initial values of the various operators in Eqs. (16a) and (16b) at t = 0 are all set to zero. As expected, it is readily shown that α in Eqs. (16a) and (16b) cancels out (as the zero of energy is immaterial), hence its actual value is irrelevant to the results. The β parameter is taken, reasonable for a qualitative analysis, to be −3 eV (0.11 hartree).45 In Sec. IV we solve Eqs. (16a) and (16b) for a simple 3-site model (N = 3) to illustrate the charge transport dynamics on the atomic sites. III. COMPUTATIONAL DETAILS
i d † † a a = [H, a1 a1 ] + S, dt 1 1 ¯ d † i † a a = [H, aj aj ], dt j j ¯
j = 2, . . . , N − 1 ,
(16a)
d † i † aN aN = [H, aN aN ] − D, dt ¯ where S, the electron source, is added to the first site to account for electrons entering the molecular sites and D, the electron drain, is subtracted from the last site to account for electrons leaving the molecular sites. We note that the nomenclature “source” and “drain” refer to the overall flow of electrons from the excited metal cylinder to Site 1, through the molecular chain to Site N, and then to the metal cylinder on the other end. At an arbitrary time during the fluctuation electrons might leave through the electron source or enter through the electron drain. When S is positive, the electrons are added from the metal lead to the 1st molecular site at the rate R as calculated in the last section; when S is negative, the electrons are drained from the 1st molecular site to the nearest metal lead, at a rate proportional to the electron population on † the 1st site, i.e., R · a1 a1 . Similar dynamics apply to the other side of the molecular chain: when D is negative, the electrons are added from the metal to the Nth molecular site, at the rate D = R; when D is positive, the electrons are drained from the Nth site to the nearest metal lead, at a rate proportional to the † electron population on the Nth site, i.e., D = R · aN aN . Equation (16a) is supplemented by analogous relations for the off-diagonal elements i d † † a a = [H, aj aj ], dt j j ¯
j = j .
(16b)
The FDTD-subgridding computational setup used in this study is shown in Figure 1. The two silver cylinders are equally sized at 7 nm radius and the distance between them is 1 nm. IS and OS stand for the inner and outer surfaces, respectively, used in the Huygens subgridding approach to pass fields information across the main-grid/subgrid boundary.36 The main grid spatial size is 0.1 nm and the main-grid/subgrid ratio is 9, making the subgrid step size approximately 0.011 nm. An infinite line light source (a point source in the 2D simulation setup) is placed 5 nm below the lower cylinder. For all simulations in this work, the light source is a Gaussian pulse with 10 fs FWHM and 400 nm wavelength. The electric field of the source is polarized along the long axis of the silver dimer. The material parameters for silver are obtained from Ref. 46. The charge and electric current densities calculated by FDTD-subgridding are then used to calculate the transition rate R as shown in Eq. (14). In this study we use a simple molecular model of three silver atoms equally spaced between the silver cylinders, as shown in the focus view in Figure 1. The silver valence 5s orbitals on the first and third atoms are used as the atomic wavefunctions φ in (14). The transition rate is subsequently used as the electron source and drain in (16a). The Heisenberg equations of motion (16a) and (16b) are then solved numerically with the commercial software Maple 12.
IV. RESULTS AND DISCUSSIONS A. Qualitative study of the charge transfer dynamics
The qualitative behavior of the charge transfer dynamics is determined by the electron source S and the electron drain D
FIG. 1. The computational setup, where IS and OS denote the inner and outer surfaces used in Huygens subgridding. On the right side is the focus view of the gap region, with the spacing of the atoms shown in between the gap.
224104-5
Hu, Ratner, and Seideman
FIG. 2. Electron dynamics on three atomic sites with R = 0.004 · sin (0.11t). † † Red curve corresponds to a1 a1 , green curve corresponds to a2 a2 , and blue † curve corresponds to a3 a3 .
in (16a). To gain a basic understanding of the dynamic system, we first study several simple examples for the electron source and electron drain. In the case where the two nanoparticles are of the same material and size and the incident light is placed so as to excite both symmetrically, the electron source and electron drain are symmetric, that is, the coupling strength R is the same for both S and D. If the nanoparticles are excited by a single frequency harmonic field, we have R = Asin (ω · t), where A is an amplitude determined by the plasmonic effect through Eq. (14), and ω coincides with the angular frequency of the incident light. We begin our survey with a weak source, A = 0.004 atomic units, and an incident frequency equal to the plasmon resonance – approximately 0.11 atomic units, corresponding to a 400 nm wavelength. With these parameters substituted in Eqs. (16a) and (16b) for N = 3 (see Figure 1), solution of the differential equations system yields the results shown in Figure 2. Here we observe electron population starting to build up on all atomic sites upon excitation while also exhibiting distinct oscillations. A Fourier transform of the dynamics in † † Figure 2 shows that the oscillations of a1 a1 and a3 a3 contain two frequencies, one at the source frequency as expected, and another one at ω = 0.157 a.u., which suggests the existence of a natural frequency for the dynamic system. This is confirmed by setting R = 0.004 · sin (0.157t), with which the dynamics is shown in Figure 3. Here we observe a striking resonance effect, where the oscillations at all three atomic sites are greatly enhanced, † and the unity population cap is quickly reached. Both a1 a1
† and a3 a3 exhibit perfect harmonic behavior at ω1 = ω3 † = 0.157 a.u., while a2 a2 oscillates at double that value at ω2 † † = 0.314 a.u. This is because a1 a1 and a3 a3 are perfectly † out of phase, hence their collective effect on a2 a2 is an os† cillation of frequency 2ω = ω2 , by which a2 a2 is readily driven. Additional calculations show that these natural system frequencies are proportional to β, the hopping integral term in the tight-binding Hamiltonian in Eq. (15). This is expected since β determines the tendency for the molecular sites to exchange electrons. The above results point to an interest-
J. Chem. Phys. 141, 224104 (2014)
FIG. 3. Electron dynamics on three atomic sites with R = 0.004 · † sin (0.157t). Red curve corresponds to a1 a1 , green curve corresponds to † † a2 a2 , and blue curve corresponds to a3 a3 .
ing control opportunity with potentially useful applications. By setting the excitation frequency ω equal to the natural system frequency of 0.157 atomic units, we generate a resonance phenomenon wherein the charge oscillation amplitudes are enhanced by close to an order of magnitude as compared to their values in Fig. 2. In addition, the electron populations grow much faster. It is important to note that the effect of the excitation frequency on the amplitude of charge oscillation is twofold: the plasmonic resonance frequency is determined by the dielectric function and geometry of the nanostructure and the molecular transfer resonance frequency is determined by the Hamiltonian of the system described in Eq. (15). Hence in experiments the ideal amplification in Figure 4 may not be realizable, requiring that a balance between the two resonance frequencies would be maintained for the charge transfer amplitude to be maximized. This result will be explored in a future study, where the focus will be on the optical control of the charge transfer dynamics. For the purpose of this discussion we simply use ω = 0.11 atomic units for maximized plasmonic resonance. For ω = 0.11, we observe the steady-state behavior without extended time by increasing the amplitude to A = 0.04.
FIG. 4. Electron dynamics on three atomic sites with R = 0.04 · sin (0.11t). † † Red curve corresponds to a1 a1 , green curve corresponds to a2 a2 , and blue † curve corresponds to a3 a3 .
224104-6
Hu, Ratner, and Seideman
FIG. 5. Electron dynamics on three atomic sites with S = R = 0.02 and D † = 2R = 0.04. Red curve corresponds to a1 a1 , green curve corresponds to † † a2 a2 , and blue curve corresponds to a3 a3 .
The dynamics is shown in Figure 4, where we see that the † † populations a1 a1 and a3 a3 quickly grow to saturation at unity and continue to oscillate in a symmetric but out of phase † fashion. The bridging site population a2 a2 oscillates at a different frequency and with a smaller amplitude. There are two major differences between the dynamics in Figure 4 and that in Figure 3: (1) all oscillation amplitudes are significantly smaller in Figure 4, which is reasonable because the excitation is off-resonance with the natural system time scale; (2) both † † a1 a1 and a3 a3 oscillate at the source frequency, not their † natural frequencies, whereas a2 a2 oscillates at double the source frequency. This is also expected because the steadystate of a driven oscillator should have the frequency of the driving force – in this case the excitation source. To conclude our qualitative study, we test the reliability of the proposed method by examining the model results in the well studied DC limit. To avoid overwhelming the molecular chain with electrons we intentionally set the source weaker than the drain. (This problem does not apply to the AC source as used in the rest of Secs. IV A and IV B, because in that case both the source and the drain rapidly oscillate and thus do not cause electrons to overcrowd on the molecular chain.) Here we set S = R = 0.02 and D = 2R = 0.04, leading to the DC dynamics shown in Figure 5, where we find that under a DC source, equivalent to a constant potential bias, each of the three sites quickly reaches a constant steady-state. The steady-state populations are all low because the drain is efficient and the source is not sufficiently strong to overwhelm the electron transfer between sites. This dynamics is expected for a constant potential bias, to which our model reduces in the limit of a DC source.
B. Quantitative rate of charge transfer
Subsection IV A has examined in detail the qualitative charge transfer dynamics governed by the Heisenberg equations of motion. In this subsection we combine the Heisenberg equations with the FDTD subgridding method of Ref. 36 to gain a quantitative appreciation of the magnitude of
J. Chem. Phys. 141, 224104 (2014)
FIG. 6. Electron dynamics on three atomic sites with R = 0.04 · sin (0.11t) † and D = 0.76S. Red curve corresponds to a1 a1 , green curve corresponds to † † a2 a2 , and blue curve corresponds to a3 a3 .
the charge transfer dynamics. As in the calculations discussed above, we take the incident frequency to equal the plasmon resonance of the silver double cylinder structure considered, ω = 0.11 atomic units, corresponding to 400 nm wavelength. At this wavelength we carry out the FDTD-subgridding simulation according to the parameters described in Sec. III. For an incident field strength of 0.030 V/nm (or translated into laser intensity ≈2.4 × 108 W/cm2 ), Eqs. (10) and (14) give R = 0.04 · sin (0.11t)atomic units. In the discussions of Figures 2–5 we set the electron drain D to have exactly the same amplitude as R so as to begin with a simple symmetric model that merely extracts the basic dynamics. In the actual simulation setup we place the light source behind one of the silver cylinders (Figure 1) to model a realistic lightdriven charge transport experiment. The result is a smaller amplitude (by about 24%) for the drain D compared to the source S, creating a similar but different dynamics from Figure 4. The electron population evolution obtained is shown in Figure 6. The charge transfer dynamics in Figure 6 is similar to that † † in Figure 4, with two major differences: (1) a1 a1 and a3 a3
now have slightly different amplitudes due to the asymmetry † between the source and the drain; (2) the a2 a2 curve can be seen to oscillate according to several frequencies, rather than the single frequency exhibited in Figure 4. A Fourier † transform of a2 a2 confirms this observation, as illustrated in Figure 7, which shows that 2 additional frequencies in † the Fourier transform of a2 a2 arise with the asymmetric simulation setup. The peak at 0.22 a.u. corresponds to the natural system time-scale as in Figure 4. The peak at 0.11 a.u. – the driving frequency, which here is taken to be the plasmon resonance frequency – illustrates that the asymmetry of the source and the drain gives rise to effect of the incident frequency on the population at all sites. Finally, the peak at 0.33 a.u. is an overtone corresponding to the sum of the incident frequency at 0.11 a.u. and the natural system frequency at 0.22 a.u. The dominating peak is still at 0.22 a.u., hence Figure 6 is not significantly different from Figure 4, but it shows that the dynamics is richer with a realistic simulation setup.
224104-7
Hu, Ratner, and Seideman
J. Chem. Phys. 141, 224104 (2014)
†
FIG. 7. Fourier transform for a2 a2 data in Figure 6. The DC component has been truncated for better scaling.
V. CONCLUSION
In this study we have proposed a viable approach to simulate light-induced charge transfer dynamics across a metalmolecule-metal junction, with a novel scheme to calculate the coupling between the metal and the molecule. To that end we combined several classical and quantum mechanical methods: (1) use of a newly developed Huygens subgridding approach36 to efficiently handle the vastly disparate size-scales of the gap and the two nanoparticles; (2) transformation of the current and charge densities computed by solution of the Maxwell equations into a quantum mechanical wavefunction; (3) use of the Fermi golden rule to calculate the rate of electron transfer from the plasmonic excitation to the molecular orbital; and (4) use of the Heisenberg equations of motion to propagate the electron population through the molecular bridge. For the purpose of this paper we have focused on clarifying the concepts and the methods, using to that end a highly simplified (although very common) model molecule. The approach is general, however, and we expect it to prove capable of simulating a wide variety of light-driven and light-controlled plasmon-mediated phenomena. These include junctions with realistic (ab initio) molecules, unidirectional current in asymmetric junctions, transport in STM tip-molecule-substrate environments, transport between semiconductor electrodes, and transport controlled by the incident excitation, e.g., pulse trains, shaped pulses, and two pathway excitations.
ACKNOWLEDGMENTS
T.S. is grateful to the National Science Foundation (Award No. CHEM-1012207) and the Department of Energy (Award No. DE-FG02-09ER16109) for support. M.A.R. thanks the National Science Foundation for support (CHE 1058896). 1 N.
J. Halas, Nano Lett. 10, 3816–3822 (2010). I. Stockman, Opt. Express 19, 22029 (2011). 3 J. Henzie, J. Lee, M. H. Lee, W. Hasan, and T. W. Odom, Annu. Rev. Phys. Chem. 60, 147–165 (2009). 4 H. Duan, A. I. Fernandez-Dominguez, M. Bosman, S. A. Maier, and J. K. W. Yang, Nano Lett. 12, 1683–1689 (2012). 5 F. L. Teixeira, IEEE Trans. Antennas Propag. 56(8), 2150 (2008). 2 M.
6 J.
Zhao, A. O. Pinchuk, J. M. McMahon, S. Li, L. K. Ausman, A. L. Atkinson, and G. C. Schatz, Acc. Chem. Res. 41(12), 1710 (2008). 7 R. Esteban, A. G. Borisov, P. Nordlander, and J. Aizpurua, Nat. Commun. 3, 825 (2012). 8 K. A. Willets and R. P. V. Duyne, Annu. Rev. Phys. Chem. 58, 267–297 (2007). 9 J. M. McMahon, A.-I. Henry, K. L. Wustholz, M. J. Natan, R. G. Freeman, R. P. V. Duyne, and G. C. Schatz, Annl. Bioanal. Chem. 394, 1819–1825 (2009). 10 S. Nie and S. R. Emory, Science 275, 1102–1106 (1997). 11 E. Fort and S. Gresillon, J. Phys. D: Appl. Phys. 41, 013001 (2008). 12 A. Pucci, F. Neubrech, D. Weber, S. Hong, T. Toury, and M. L. d. l. Chapelle, Phys. Status Solidi B 247, 2071–2074 (2010). 13 S. Pillai, K. R. Catchpole, T. Trupke, and M. A. Green, J. Appl. Phys. 101, 093105 (2007). 14 T. Vo-Dinh, A. M. Fales, G. D. Griffin, C. G. Khoury, Y. Liu, H. Ngo, S. J. Norton, J. K. Register, H.-N. Wang, and H. Yuan, Nanoscale 5, 10127 (2013). 15 S. A. Maier and H. A. Atwater, J. Appl. Phys. 98, 011101 (2005). 16 P. Berini, R. Charbonneau, N. Lahoud, and G. Mattiussi, J. Appl. Phys. 98, 043109 (2005). 17 S. A. Maier, P. E. Barclay, T. J. Johnson, M. D. Friedman, and O. Painter, Appl. Phys. Lett. 86, 071103 (2005). 18 E. Ozbay, Science 311(5758), 189–193 (2006). 19 S. K. Gray and T. Kupka, Phys. Rev. B 68, 045415 (2003). 20 J. P. Kottmann and O. J. F. Martin, Opt. Express 8(12), 655–663 (2001). 21 D. C. Marinica, A. K. Kazansky, P. Nordlander, J. Aizpurua, and A. G. Borisov, Nano Lett. 12(3), 1333–1339 (2012). 22 E. Hao and G. C. Schatz, J. Chem. Phys. 120, 357 (2004). 23 K. Lopata and D. Neuhauser, J. Chem. Phys. 130, 104707 (2009). 24 A. J. Haes, C. L. Haynes, A. D. McFarland, G. C. Schatz, R. P. V. Duyne, and S. Zou, MRS Bull. 30, 368 (2005). 25 M. Sukharev and T. Seideman, J. Chem. Phys. 124(14), 144707 (2006). 26 J. C. Cuevas and E. Scheer, Molecular Electronics (World Scientific Publishing Company, 2010). 27 S. Datta, Quantum Transport: Atom to Transistor (Cambridge University Press, 2005). 28 M. D. Ventra, Electrical Transport in Nanoscale Systems (Cambridge University Press, 2008). 29 A. Nitzan and M. A. Ratner, Science 300, 1384 (2003). 30 J. M. Klingsporn, M. D. Sonntag, T. Seideman, and R. P. V. Duyne, J. Phys. Chem. Lett. 5, 106–110 (2014). 31 P. Dombi, A. Horl, P. Racz, I. Marton, A. Trugler, J. R. Krenn, and U. Hohenester, Nano Lett. 13, 674–678 (2013). 32 P. Banerjee, D. Conklin, S. Nanayakkara, T.-H. Park, M. J. Therien, and D. A. Bonnell, ACS Nano 4, 1019–1025 (2010). 33 S. W. Wu, N. Ogawa, and W. Ho, Science 312, 1362–1365 (2006). 34 B. D. Fainberg, M. Jouravlev, and A. Nitzan, Phys. Rev. B 76, 245329 (2007). 35 B. D. Fainberg, M. Sukharev, T.-H. Park, and M. Galperin, Phys. Rev. B 83, 205425 (2011). 36 Z. Hu, M. A. Ratner, and T. Seideman, J. Chem. Phys. 137(20), 204111 (2012).
224104-8
Hu, Ratner, and Seideman
37 A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-
Difference Time-Domain Method, 3rd ed. (Artech House, Boston, 2005). 38 J.-P. Berenger, J. Comp. Phys. 230, 5635 (2011). 39 K. S. Yee, IEEE Trans. Antennas Propag. 14, 302–307 (1966). 40 J.-P. Berenger, IEEE Trans. Antennas Propag. 57, 3860–3867 (2009). 41 A. Salomon, R. J. Gordon, Y. Prior, T. Seideman, and M. Sukharev, Phys. Rev. Lett. 109, 073002 (2012).
J. Chem. Phys. 141, 224104 (2014) 42 M.
Sukharev, T. Seideman, R. J. Gordon, A. Salomon, and Y. Prior, ACS Nano 8(1), 807–817 (2014). 43 Z. Hu, M. A. Ratner, and T. Seideman, J. Chem. Phys. 136, 054504 (2012). 44 V. Mujica, M. Kemp, and M. A. Ratner, J. Chem. Phys. 101, 6849 (1994). 45 Y. Xie and J. A. Blackman, Phys. Rev. B 64, 195115 (2001). 46 A. D. Rakic, A. B. Djurisic, J. M. Elazar, and M. L. Majewski, Appl. Opt. 37(22), 5271–5283 (1998).
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
