Combined Linear Response Quantum Mechanics and Classical

Article pubs.acs.org/JPCA

Combined Linear Response Quantum Mechanics and Classical Electrodynamics (QM/ED) Method for the Calculation of SurfaceEnhanced Raman Spectra Jonathan Mullin and George C. Schatz* Department of Chemistry, Northwestern University, 2145 Sheridan Road, Evanston, Illinois 60208-3113, United States ABSTRACT: A multiscale method is presented that allows for evaluation of plasmon-enhanced optical properties of nanoparticle/molecule complexes with no additional cost compared to standard electrodynamics (ED) and linear response quantum mechanics (QM) calculations for the particle and molecule, respectively, but with polarization and orientation effects automatically described. The approach first calculates the total field of the nanoparticle by ED using the finite difference time domain (FDTD) method. The field intensity in the frequency domain as a function of distance from the nanoparticle is calculated via a Fourier transform. The molecular optical properties are then calculated with QM in the frequency domain in the presence of the total field of the nanoparticle. Back-coupling due to dipolar reradiation effects is included in the single-molecule plane wave approximation. The effects of polarization and partial orientation averaging are considered. The QM/ED method is evaluated for the wellcharacterized test case of surface-enhanced Raman scattering (SERS) of pyridine bound to silver, as well as for the resonant Raman chromophore rhodamine 6G. The electromagnetic contribution to the enhancement factor is 104 for pyridine and 102 for rhodamine 6G.

1. INTRODUCTION The ability to generate a well-defined and unique molecular signature is key for any spectroscopic technique that will be of interest for the identification and detection of chemical systems. Two common approaches for this purpose are absorption and Raman spectroscopy1 as these provide unique signals in response to an external radiation field. Signal to noise is an issue with Raman spectroscopy, but the development of surface-enhanced Raman scattering (SERS) circumvents this problem for many applications.2,3 The enhancement comes from both chemical and electromagnetic effects, with the latter arising from plasmon-enhanced local fields on the surface of the nanoparticles. The challenge to build a theory for calculating SERS optical response that includes both electrodynamics (ED) and quantum mechanics (QM) stems from the order of magnitude differences in length scales needed for the two calculations. Purely chemical models of SERS based on QM methods have thus been limited to ∼1 nm in size, including the silver nanoparticle, while electromagnetic field evaluations are usually based on grids or finite elements that have 1 nm dimensions at the minimum.4 Theoretical treatments of SERS have thus usually taken one of two paths. One approach neglects the chemical enhancement (CHEM) and focuses on the dominant electromagnetic enhancement (EM).5 Other work determines the CHEM enhancement using a small cluster model of the nanoparticle.4 Although small clusters can have plasmon-like optical spectra, there are usually substantial sizedependent shifts compared to what is found for larger © 2012 American Chemical Society

nanoparticles, and the fields associated with these clusters are only vaguely related to those of larger nanoparticles. In the EM model, the EM enhancement is approximately | Eloc(ω)|2|Eloc(ω′)|2, where Eloc is the local field enhancement and ω and ω′ are the incident and Stokes-shifted frequencies. Further, this formula is often approximated in the zero Stokes shift limit, where the enhancement factor simplifies to | Eloc(ω)|4.6−8 This assumes that dipole reradiation (DR) effects can be evaluated using the plane wave (PW) response of the nanoparticle, instead of that arising from the induced dipole in the adsorbed molecule at the Stokes frequency. This issue has been known since the work of Kerker et al.,9 and recently, the comparison of DR and PW results was explicitly explored.10 It was found that the PW approximation is adequate in the quasistatic (small particle limit) approximation when observing at certain locations such as in the back-scattered direction,11−15 but for large particles, there can be important differences. Corni and Tomasi explicitly coupled electronic structure calculations in the molecule to an ED calculation of the particle in the frequency domain through effective charges that were included in the molecular Hamiltonian in the quasi-static approximation.16 Neuhauser developed a local two-level random phase approximation (RPA) model for density matrix evolution to evaluate the molecule’s population transfer rate, Received: September 12, 2011 Revised: January 26, 2012 Published: January 27, 2012 1931

dx.doi.org/10.1021/jp2087829 | J. Phys. Chem. A 2012, 116, 1931−1938

The Journal of Physical Chemistry A

Article

with FDTD used to model the surface plasmons.17 Masiello and Schatz applied a many body Greens function and the quasistatic approximation to evaluate plasmon-enhanced molecular absorption with electrodynamics and electronic structure explicitly coupled.18 Morton and Jensen have described a discrete interaction model/quantum mechanics method to explicitly model nanoparticle interaction.19,20 With an atomistic representation of the nanoparticle, the influence of the local environment is an explicit model for the optical properties of a molecule interacting with the nanoparticle surface. Chen and Schatz described a hybrid QM/ED method in the time domain using RT-TDDFT and FDTD.21 Although time domain calculations were used for both calculations, it was necessary to transform the field to the frequency domain when passing between FDTD and RT-TDDFT due to the different boundary conditions used in the two calculations (i.e., the FDTD calculations take the initial field to be an asymptotically defined plane wave, while the RT-TDDFT uses a field at the molecule that is turned on with a femtosecond pulse at the beginning of the calculation). These methods and other methods have recently been reviewed.22,23 In this work, we present a formalism that couples linear response TDDFT in the frequency domain with the particle modeled by FDTD in the time domain. This avoids the back Fourier transform used in the Chen et al. work, which is a minor issue, but of more interest is that it provides access to alternative electronic structure methods, allowing for explicitly resolving the QM eigenstates if desired. Our theory also includes for explicit treatment of the “back coupling” using the PW approximation,21 and this can be extended to explicit DR calculations if desired. Also, we show how polarization averaging and partial orientation averaging are implemented to mimic the SERS experiment. The formalism is general for a full range of optical response properties that can be studied with frequency domain TDDFT. The formalism is briefly outlined in section 2, while section 3 outlines the computational methodology. The numerical results are discussed in section 4, while section 5 discusses the implications of our QM/ED methods for linear optical response and its relation to other methods.

difference. Here, we have added an energy broadening parameter, γ, to include for the finite lifetime of the excited state, using the same broadening for each transition. This approximation has been described previously,4 along with theories (not implemented here) that provide a rigorous way to calculate γ.18 With the broadening included, we need to consider both the real and imaginary responses of the system; this is described below. The change in the effective potential is given as pert

V st′ eff (ω) = Vst (ω) +

(2.1.3)

uv pert

where V (ω) = E(ω)rα, with the E(ω) field coming from the FDTD calculation, as described in section 2.2, and the coupling term K has a Coulombic and exchange correlation term. By substituting eq 2.1.3 into eq 2.1.2, a set of linear equations can be used to describe the first-order density matrix as nt − ns pert ′ (ω)] [V (ω) + ∑ Kst , uvPuv Pst′ (ω) = ω − iγ − ωst st uv

(2.1.4)

which can be solved iteratively without explicitly constructing the coupling matrix.28 Equation 2.1.4 is an alternative form of the RPA equations,29 presented here with the additional damping term. Expressing the first-order density matrix as Pst′ (ω) = Pst′ R (ω) + iPst′ I(ω) R

(2.1.5) I

where P′ is the real contribution and P′ is the imaginary contribution, the real polarizability is R (ω) = − Tr[H αR (ω)P′βR (ω) − H αI(ω)P′βI (ω)] α αβ (2.1.6)

while the imaginary part is α Iαβ(ω) = −Tr[H αR (ω)P′βI (ω) + H αI(ω)P′βR (ω)] (2.1.7) α

Here, H is the dipole operator associated with the direction α, PβR(ω) is the real density response in the direction β, and P′βI is the corresponding imaginary part. The imaginary polarizability is proportional to the absorption cross section via

2. THEORY 2.1. Linear Response Time-Dependent Density Functional Theory. The time-dependent Kohn−Sham (TD-KS) equations24−27 allow for an external time-dependent perturbation in order to study the dynamic response of a system. Response theory is used for small perturbations of a system to yield perturbative solutions to the TD-KS equations. This work focuses on linear response properties; therefore, we consider the first-order change in the density

σ(ω) =

4πω I α̅ αβ(ω) c

(2.1.8)

where c is the speed of light. α̅ is the isotropic polarizability, given by α̅ =

ρ′(r , ω) = ∑ Pst′ (ω)φs(r )φ*t (r ) s,t

′ (ω) ∑ Kst , uvPuv

1 (αxx + α yy + αzz) 3

(2.1.9)

For some of the expressions that we use later, we will also use the anisotropic polarizability, which is given by

(2.1.1)

where P′ is the first-order density matrix and φs(r) are the KS orbitals. Expansion of the KS equations to first order with a frequency-dependent perturbation, Vper(ω), yields the firstorder density matrix nt − ns V st′ eff (ω) Pst′ (ω) = ω + iγ − ωst (2.1.2)

β(α)2 =

1 [(αxx − α yy)2 + (α yy − αzz)2 2 + (αzz − αxx)2 + 6(α2xy + α2xz + α2yz)] (2.1.10)

2.2. Finite-Difference Time-Domain (FDTD) method. A nanostructure in FDTD is discretized into small cubic blocks that are each characterized by a dielectric permittivity, ε(r), and

where the change in the effective potential is V′eff, n is the orbital occupation number, and ωst is the orbital energy 1932



Article

a magnetic permeability, μ(r), which is unity in our case. In this formalism, the Maxwell equations30 ε(r )

∂⇀ E (r , t ) = ∇ × ⇀ H (r , t ) − ⇀ J (r , t ) ∂t

(2.2.1)

μ(r )

∂⇀ H (r , t ) = −∇ × ⇀ E (r , t ) ∂t

(2.2.2)

system to the particle-fixed coordinate system. To simplify things, we assume that the polarizability tensors of both the particle and molecule are diagonal and that the coordinate frames are aligned. In this case, when the direction n̂ between the molecule and particle is aligned with the applied field direction, then the dipole field of the particle, evaluated at the molecule, is given by

are solved in the time domain, and the evolution of the electric field, ⇀ E (r,t), magnetic field, ⇀ H (r,t), and electric current density, ⇀ J (r,t), are obtained. Fourier transforms are then used to determine frequency domain electromagnetic properties. The incident field is taken to be a time-shifted Gaussian wave, which allows a broad spectral range to be studied.31 For a given observation point, the combination of the scattered field, ⇀ E p(r,ω), and the incident field, ⇀ E 0(r,ω), yields the total electric field ⇀ E (r , ω) = ⇀ E (r , ω) = ⇀ E (r , ω) + ⇀ E (r , ω) total

p

0

Ep =

R

μp =

3 cos(0) − 1 3

R

μp =

2 R3

μp

⎛ 2 ⎞ μm = α m⎜E0 + 3 μp⎟ ⎠ ⎝ R

(2.3.3)

(2.3.4)

Similarly, the dipole of the particle is related to the dipole of the molecule via

(2.2.3)

⎞ ⎛ 2 μp = α p⎜E0 + 3 μm⎟ ⎠ ⎝ R

(2.3.5)

By substituting eq 2.3.4 into eq 2.3.5, the total dipole of the molecule is

μm =

⎛ ⎞ 2 α m⎜1 + 3 α p⎟E0 ⎝ ⎠ R 2

2

1 − 3 αm 3 αp R R

(2.3.6)

Because, in general, αm ≪ R3, the denominator may be approximated as unity.32 This shows that the molecular dipole is enhanced by a factor of g = (1 + (2/R3)αp) compared to what it would be in the absence of the particle. To evaluate the Raman intensity, we must evaluate the Stokes side band associated with this induced dipole. This is obtained from

dμ p μm = μ 0 + Q dQ

(2.3.7)

therefore, the Raman emitting dipole is μ′m = Q

⎛ 2 ⎞ = Q α′m⎜1 + 3 α p⎟E0 ⎝ ⎠ dQ R

dμ m

(2.3.8)

The dipole field associated with this now interacts with the particle at the Stokes shifted frequency ω′, leading the an induced dipole in the particle equal to μ′p =

(2.3.1)

and similarly for the particle μp = α p(E0 + Em)

3

where the distance between the particle and the molecule is R. The induced dipole of the molecule, μm, is then expressed in terms of the particle’s dipole μp via

The local field enhancement is measured by ⇀ E (r,ω), which is a complex tensor that is dependent on the polarization and propagation direction of the incident light. This field is contained in the perturbation, Vpert, in section 2.1. 2.3. Toward a Combined Quantum Mechanics/ Classical Electrodynamics (QM/ED) Method. To understand how the QM and ED calculations can be combined, we now present a quasistatic description of the QM/ED method using point dipole approximations for the fields induced in the nanoparticle and in the molecule. Later, we show how this can be generalized beyond the point dipole approximation and beyond the quasistatic approximation, but the present derivation, which is different from an earlier derivation presented by Masiello and Schatz,18 provides important insights. Given the quasistatic/point dipole assumptions, the Raman scattering associated with a molecule near a silver sphere can be expressed in two steps. First, the applied electric field induces a dipole in the particle, and the combination of the applied field and induced field act on the molecule to induce an oscillating dipole that has a Raman side band. Second, this Raman dipole acts as an emitter that interacts with the particle to induce a dipole in the particle, and emission by molecular and particle dipoles defines the signal that is detected. Note that this description assumes that there is only one molecule adsorbed on the particle; therefore, the back-coupling is only treated in the zero coverage limit.9,10 Considering the first step, the induced dipole in the molecule can be expressed as the product of the molecular polarizability, αm, and the electric field E from eq 2.2.3 as follows

μm = α m(E0 + E p)

3nn̂ ̂ − 1

2 R3

α p(ω′)μ′m

To include for DR effects, the input for the electrodynamics calculation would be eq 2.3.8, and the intensity would be extracted from the outgoing scattered wave. Implicitly, this includes for effects arising from the induced dipole μp′ as well as any higher multipoles induced in the particle. Returning to the dipole treatment, the total emitting dipole is

(2.3.2)

where Ep is the field due to the induced dipole in the particle evaluated at the position of the molecule and Em is the field of the molecule evaluated at the position of the particle. Note that there is a time dependence exp(−iωt) to these dipoles and fields that we have omitted. These equations are coupled as Ep depends on μp and Em depends on μm. To express these relations, we need to relate the molecule-fixed coordinate

⎛ ⎞⎛ ⎞ 2 2 μ′ = Q α′m⎜1 + 3 α p(ω)⎟⎜1 + 3 α p(ω′)⎟E0 ⎝ ⎠⎝ ⎠ R R 1933

(2.3.9)



Article

at the origin, and the molecule is along the x axis at a distance R from the surface of the particle. The molecular orientation is assumed to be fixed in space so that no orientation averaging is needed. In this case, the Raman differential cross section, dσ/ dΩ, involves an average over the x and y polarization directions as follows

This expression, when converted to a Raman intensity, leads to the SERS intensity expression I = |μ′|2 2 ⎞⎛ ⎞ 2 2 2 ′ 2 ⎛ = Q |α m| ⎜1 + 3 α p(ω)⎟⎜1 + 3 α p(ω′)⎟ I0 ⎝ ⎠⎝ ⎠ R R

dσ (180°) dΩ x x 1 = K [|⇀ E (ω′)·α′m·⇀ E (ω)|2 2 p 2E0 y y + |⇀ E (ω′)·α′ ·⇀ E (ω)|2 ]

I(180°) =

(2.3.10)

This expression involves the well-known SERS enhancement factor Enh = |g (ω)g (ω′)|2

(2.3.11)

m

and it is easily generalized beyond the dipole approximation by taking g(ω) = E(ω)/E0, where E(ω) is the field from eq 2.2.3. Note that DR effects10 are included in this formula in the plane wave approximation. Although this expression was developed for the case where the polarization vector and axis between the particle and molecule are collinear, it is not hard to write a completely general expression for the intensity I = Q 2 |⇀ E (ω′)·α′ ·⇀ E (ω)|2 I /E 4 m

0

0

Here, ⇀ E x and ⇀ E y are the local fields at the position of the molecule that are respectively associated with x- and y-polarized applied fields. A further step is to consider that the molecule may rotate about the axis along which the molecule is bound to the particle, the x axis, as described here. Considering this effect on one of the terms in eq 2.3.17 is done by applying a rotation matrix about the x axis

(2.3.12)

⎛1 0 0 ⎞ ⎛1 0 0 ⎞ ⎜ ⎟ ⎜ ⎟ y ⇀ E (ω′)·⎜ 0 cos(ϕ) sin(ϕ) ⎟α′m⎜ 0 cos(ϕ) −sin(ϕ)⎟ ⎜ ⎟ ⎜ ⎟ ⎝ 0 −sin(ϕ) cos(ϕ)⎠ ⎝ 0 sin(ϕ) cos(ϕ) ⎠ y ·⇀ E (ω)

Polarization averages can also be determined from this expression, as we describe below. Also note that this same derivation leads to an enhanced polarization α = α ·⇀ E (ω)/E Enh

m

0

(2.3.18)

(2.3.13)

Here, the azimuthal angle ϕ is the angle between the lab-fixed y axis and the rotating molecular y′ axis. Averaging eq 2.3.17 over ϕ using eq 2.3.18 leads to the following expression, which we will call the azimuthal polarization averaged intensity

from which the enhanced absorption cross section can be derived. 2.3.2.1. Evaluation of Optical Properties. Absorption Spectra. The gas-phase absorption cross section is obtained from the imaginary part of the polarizability 4πω I σ(ω) = α Enh c

dσ (180°) dΩ ⎡ 1 = K α′ 2 ([E x]4 + [Exy]4 ) 2 p⎢⎣ m, xx x 2E

I(180°) = (2.3.14)

A similar expression works for the enhanced cross section, except that the enhanced polarization in eq 2.3.13 is used. This cross section includes for both molecule and particle contributions to absorption. 2.3.2.2. Raman Spectra. The gas-phase Raman intensity, assuming a backward or forward scattering geometry, is given by33,34 I(180°) =

(2.3.17)

0

1 2 + (2α′m2 + β′̃ m)({[Eyx]2 + [Ezx]2 }2 2 2 2 + {[Eyy]2 + [Ezy]2 }2) + 2(α′m, xy + α′m, xz + α′m, xxα̃ ′){[Exx]2 ([Eyx]2 + [Ezx]2 ) ⎤ + [Exy]2 ([Eyy]2 + [Ezy]2 )}⎥ ⎦

⎡ 90α′2 + 14β(α′)2 ⎤ dσ ⎥ (180°) = K p⎢ ⎢⎣ ⎥⎦ 90 dΩ

(2.3.19)

where the superscript on the symbol E denotes the polarization direction while the subscript indicates the space-fixed laboratory component of the field. Also, for the azimuthal average, the trace and anisotropy from eqs 2.1.9 and 2.1.10 are replaced by partial traces and anisotropies as follows

(2.3.15)

where the isotropic eq 2.1.8 and anisotropic eq 2.1.9 polarizability tensors are used. K p is independent of experimental setup; however, it depends on both the incident and scattered frequencies as follows

α̃′m =

2

1 π h K p = 2 (υ̃ in − υ̃ p)4 2 1 exp[ − − hc υ̃ p/kBT ] 8π c υ̃ p ε0

1 (α̃′m, yy + α̃′m, zz) 2

1 2 2 2 β′̃ m = (4α̃′m, yz + (α̃′m, yy − α̃′m, zz) ) 4

(2.3.16)

where for the pth vibrational mode, the υ̃in and υ̃p are the frequencies of the incident and scattered light. The surface Raman intensities will be determined with the assumption that the wavevector of the incident light is along the z axis, the center of the particle (assumed spherical) is taken

(2.3.20)

(2.3.21)

3. COMPUTATIONAL METHODS 3.1. Electrodynamics: Finite Difference Time Domain (FDTD). For the calculation of the total electric field, the JFDTD3D package35 was used. A Cartesian coordinate system 1934



Article

local field for x polarization is about twice that for y polarization. The total field strength is roughly halved over the 1−5 nm distance. However, the field at 1 nm separation is over an order of magnitude stronger than the asymptotic laser field; therefore, the 20 nm silver sphere is a reasonable model for surface-enhanced optical properties as used in section 4.2 and 4.3. The results also show that the total field is 2 orders of magnitude stronger along the polarization direction than that along the orthogonal directions; therefore, for all practical purposes, the off-diagonal terms may be neglected for ⇀ E (r,ω). 4.2. Raman Spectra of Pyridine. Pyridine was chosen because it has been extensively studied and has a typical SERS enhancement factor.46,47 The gas-phase normal Raman spectrum of pyridine is presented in Figure 1, here assuming

with a unit length of 1 nm was used for the silver sphere calculations. The grid size was 0.25 nm with a box length of 40 nm. A silver sphere of 20 nm diameter was placed at the origin, with a box side length of 40 nm yielding 4096000 grid cells. Each box in the simulation is characterized by ε and μ values that were determined experimentally.36 The incident light pulse originated from the z direction. Parameters for the pulse were chosen to cover the visible spectrum (300−800 nm). The incident light pulse was injected from the plane at z = −16 using the functional form 2 2 E0(t ) = e−(t − to) / σ sin(ωot )

(3.1.1)

where τ0 = 10.0 fs, σ = 0.7 fs, and ω0 = 600.0 nm. In order to absorb the light, perfectly matched layers were placed above and below boundary planes at z = ±20. Periodic boundary conditions were applied in the x and y directions. A simulation time of 100 fs, with a time step of 4.57 × 10−4 fs, was used. A separate calculation was done for a pulse polarization in the x and y directions. The electric field vectors were output every 10 steps at 1, 2, 3, 4, and 5 nm distances from the nanoparticle surface along the x axis. The electric field was Fouriertransformed into a frequency-dependent electric field. 3.2. Quantum Mechanics: Time-Dependent Density Functional Theory (TDDFT). All calculations were done using a locally modified version of NWChem 5.1.37 The molecular structures were determined with DFT functionals. The generalized gradient approximation (GGA) Perdew− Burke (PBE) exchange−correlation functional38 was used to verify previous work on pyridine with DFT-specific DZVP basis sets.39 Rhodamine 6G (R-6G) was calculated with the BP86 GGA functional40,41 along with DFT-specific TZVP basis sets39 to be consistent with previous studies.42 All structures were optimized via the xfine grid, with an energy convergence of 10−8 h. A gradient convergence criterion of 10−5 h and an energy convergence criterion of 5 × 10−6 h were used in order to obtain converged geometries. Hessians were calculated to verify the positive definite nature of the minima. The diagonalized Hessian was used for normal mode and frequency analysis. The real and imaginary polarizabilities were computed using the AORESPONSE28,43,44 module. This assumes a short-time approximation to evaluate the perturbed density matrix (and, from this, the frequency-dependent polarizability). A damping parameter of Γ = 0.0037 au was used (consistent with previous work3). It was necessary to set nolevelshifting for the polarization calculation. The Raman intensity was determined through finite differencing of normal mode displacements. The infinitely narrow intensity from the TDDFT calculation was broadened by a Lorentzian line shape with a 20 cm−1 width. All visualizations were done in MacMolPlt.45

Figure 1. Numerical calculation of a normal Raman spectrum at 514.5 nm for gas-phase pyridine. The spectrum was broadened with a Lorentzian line shape with a 10 cm−1 width.

a 514.5 nm excitation wavelength. The two major peaks at 976 and 1020 cm−1 are the well-known ring breathing modes. The secondary peak at 589 cm−1 is an example of a ring deformation mode. Ring stretching modes appear at 1207 and 1588 cm−1. These modes are illustrated in Figure 2. The isolated pyridine Hessian is used for the SERS calculation as we are imaging that the pyridine is removed from the surface by 1−5 nm. Also, the C2v axis of the molecule is taken along the x axis, with the nitrogen atom pointing toward the nanoparticle. The frequencies of the modes are about 10 cm−1 lower than the experimental SERS results,48 but we are neglecting solvent effects,23 and our SERS model will not have chemical effects as currently implemented. Figure 3 shows the calculated spatially fixed orientation SERS spectrum at a laser excitation of 355 nm. The plane of the pyridine is assumed to be fixed in the xy plane for this calculation. The ring breathing modes are enhanced by a factor on the order of 105 when the pyridine is 1 nm from the surface. The 1588 cm−1 mode peak is enhanced but not to same degree as the ring breathing modes. These other modes show changes that reflect whether the modes involve significant motion in the x direction or not.46 Figure 4 shows the calculated azimuthal average SERS spectrum at a laser excitation of 355 nm. The pyridine C2v axis is taken to be fixed along the x axis; therefore, it corresponds to the R = ∞ limit. As with the spectrum in Figure 3, the two main peaks are enhanced by a factor on the order of 104. The overall spectrum is similar to what has been seen

4. RESULTS 4.1. The 20 nm Silver Sphere Calculation. The total field ⇀ ( E (r,ω)) was calculated via FDTD simulations for x and y polarization, as described in section 3.1. This is the same FDTD calculation as that used by Chen and Schatz. We used this in order to eliminate the electrodynamics calculation as a source E (r,ω)| is a of differences seen between the two methods. |⇀ function of distance for 1−5 nm separations of the center of mass of the molecule from the nanoparticle surface for the 355 nm wavelength. Because the molecule is along the x axis, the 1935



Article

Figure 2. Important vibrational normal modes of pyridine.

Schatz using RT-TDDFT and by Yang and Schatz using static DFT.49 The azimuthal average spectrum is a new result to this work, and we see that it is mostly similar to the spatially fixed spectrum but with stronger intensity at 1588 cm−1, and both components of the doublet near 600 cm−1 are visible. These changes compared to the spatial spectrum are in better agreement with experiment, although the bluer component of the doublet near 600 cm−1 is less intense in the experiment. Strong variation of the intensities with R is apparent in Figures 3 and 4, where it is seen that the peak heights decrease rapidly with the molecule−particle separation. At a 5 nm separation, the electromagnetic enhancements are down by a factor of approximately 50 compared to the results at 1 nm separation. It should be noted that this neglects all chemical effects and is purely an electromagnetic enhancement. As a comparison, a generalized Mie theory evaluation of the maximum |E|4 enhancement at 1 nm separation using the Johnson and Christy values36 yields 1.15 × 105 at 357 nm. This shows that the QM/ED approach gives results that are within a factor of 10 of the purely electromagnetic results; however, polarization and orientation effects that are not included in the purely electrodynamic calculation do lead to the factor of 10 for quantitative calculations QM/ED to be preferred. In addition, there is no additional computational effort required in the frequency domain version of the theory compared to separately doing the electrodynamics and electronic structure calculations. This is in contrast to the time domain approach, where each polarization and each orientation require new calculations. Also note that chemical enhancements may be included in QM/ED by including explicit silver atoms in the QM region; we plan to do this in future work. 4.3. Raman Spectra of Rhodamine 6G (R-6G). R-6G has the largest known Raman cross section.42,50 To calculate the SERS spectrum, we assume that the xanthene ring of R-6G is in the xy plane with the long axis of the ring parallel to the surface. The gas-phase normal Raman spectrum at 1064 nm is shown in Figure 5. The spectrum shows approximately a 20 cm−1 shift in vibrational frequencies compared to the experimental (SERS) values, which may come from the lack of chemical interactions with the particle or from solvent effects. The mode at 1638 cm−1 is a xanthene ring stretch with in-plane C−H bends. There are modes at 1604 and 1545 cm−1 that are a xanthene ring stretch with in-plane N−H bends, but they are not significant in the gas-phase spectrum. The 1492 and 1515 cm−1

Figure 3. Spatially fixed surface-enhanced resonance Raman spectrum of pyridine at 355 nm with pyridine for a series of separations (1−5 nm) from the Ag sphere (20 nm diameter). The spectrum was broadened with a Lorentzian line shape and a 10 cm−1 width.

Figure 4. Azimuthal averaged surface-enhanced resonance Raman spectrum of pyridine at 355 nm for a series of separations (1−5 nm) from the Ag sphere (20 nm diameter). The spectrum was broadened with a Lorentzian line shape and a 10 cm−1 width.

experimentally, though the relative heights of the two main peaks are reversed. The appearance of the spatially fixed spectrum is quite similar to the spectra obtained by Chen and 1936



Article

Figure 5. Gas-phase normal Raman spectrum at 1064 nm for isolated R-6G. The spectrum was broadened with a Lorentzian line shape and a 10 cm−1 width.

Figure 7. Azimuthal averaged surface-enhanced resonance Raman spectrum at 476 nm of R-6G for a series of separations (1−5 nm) from the Ag sphere (20 nm diameter). The spectrum was broadened with a Lorentzian line shape and a 10 cm−1 width. The inset was taken from ref 50, used with permission of the American Chemical Society.

modes are xanthene ring stretches with C−H and N−H bends. The 1342 cm−1 mode is a xanthene ring stretch. The mode at 1181 cm−1 is an in-plane xanthene ring deformation with C−H and N−H bends. The results in Figure 5 are in good agreement with those from Jensen and Schatz based on comparable DFT calculations but using the ADF code rather than NWChem.37,42 Figure 6 shows the spatially fixed SERS spectrum of R-6G at 476 nm. This wavelength corresponds to the wavelength of

effects. Despite these effects, there is generally good agreement between the calculated spectrum and experiment, with the pattern of the most intense peaks being similar. The integrated SERS enhancements compared to the gas-phase spectrum are on the order of 107. This consists of effects of ∼104 from molecular resonance enhancement along with a roughly ∼101 effect from the ω4 prefactor in the cross section that arises because the nonresonant spectra are at 1064 nm while the resonant spectra are at 476 nm. This leaves ∼102 of the enhancement stemming from electromagnetic effects at 476 nm. The enhancement falls off by a factor of 40 as the R-6G is pulled back to a 5 nm separation from the nanoparticle. The 105 nonelectromagnetic enhancement factor for R-6G (i.e., 104 × 101) is similar to a value estimated at the same wavelength by Jensen and Schatz using ADF.42 However, the ∼5 × 10−24 cm2 value of the total nonelectromagnetic cross section is well below the ∼2.3 × 10−22 cm2 value estimated from FSRS (femtosecond stimulated Raman scattering) measurements by Shim et al.51 This suggests that the 0.1 eV width factor should be smaller by at least an order of magnitude in order to properly describe the resonance Raman R-6G cross section. Of course, the 0.1 eV width factor was estimated for the case where R-6G interacts with a plasmonic particle; therefore, the overestimate for R-6G in the absence of a plasmonic particle is understandable.

Figure 6. Spatially fixed surface-enhanced resonance Raman spectrum of R-6G at 476 nm for a series of separations (1−5 nm) from the Ag sphere (20 nm diameter). The spectrum was broadened with a Lorentzian line shape and a 10 cm−1 width.

5. DISCUSSION AND CONCLUSIONS The combined frequency domain QM and time domain ED method (QM/ED) that we have presented has been demonstrated to reproduce SERS spectra previously generated in the time domain within the limitations of the different electronic structure methods that are used for the frequency and time domain calculations. This method is extremely sensitive to the orientation of the molecule and the relative distance from the molecule to the nanoparticle. Polarization effects are also important, and slight changes in the propagation and polarization direction of the incident light have a large spectroscopic effect. One major benefit of this method compared to previous work with RT-TDDFT is the ability to select experimentally relevant wavelengths for study. This can save significant

maximum absorption of isolated R-6G. The plasmon resonance of the 20 nm silver sphere is at a wavelength of 355 nm. Because of this mismatch, the effects of the resonant molecular spectra are emphasized in this presentation. The most intense peak is at 1604 cm−1, followed by peaks at 1638, 1545, and 1342 cm−1. Figure 7 shows the surface-enhanced resonance Raman (azimuthal average) spectrum of R-6G at 476 nm. The major difference compared to the spatially fixed spectrum is the larger intensity of the 1492 and 1465 cm−1 peaks. The experimental results (in the inset) refer to the 514 nm wavelength, as noted above, and we note that there is a 50 cm−1 width in the experimental spectrum that probably arises from heterogeneous 1937



Article

(12) Shegai, T.; Li, Z.; Dadosh, T.; Zhang, Z.; Xu, H.; Haran, G. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 16448. (13) Suh, Y.; Schenter, G.; Zhu, L.; Lu, H. Ultramicroscopy 2003, 97, 89. (14) Pristinski, D.; Le Ru, E. C.; Tan, S.; Sukhishvili, S.; Du, H. Opt. Express 2008, 16, 20117. (15) Zhu, L.; Schenter, G. K.; Micic, M.; Suh, Y. D.; Klymyshyn, N.; Lu, H. P. Proc. SPIE 2003, 4962, 70. (16) Corni, S.; Tomasi, J. J. Chem. Phys. 2002, 116, 1156. (17) Lopata, K.; Neuhauser, D. J. Chem. Phys. 2009, 130, 104707. (18) Masiello, D. J.; Schatz, G. C. Phys. Rev. A 2008, 78, 042505. (19) Morton, S. M.; Jensen, L. J. Chem. Phys. 2011, 135, 134103. (20) Morton, S. M.; Jensen, L. J. Chem. Phys. 2010, 133, 074103. (21) Chen, H.; McMahon, J.; Ratner, M.; Schatz, G. J. Phys. Chem. C 2010, 161. (22) Tong, L.; Zhu, T.; Liu, Z. Chem. Soc. Rev. 2011, 40, 1296. (23) Morton, S. M.; Silverstein, D. W.; Jensen, L. Chem. Rev. 2011, 111, 3962. (24) Gross, E.; Kohn, W. Density functional theory of Many-Fermion systems; Advances in Quantum Chemistry; Elsevier: New York, 1990; Vol. 21, p 255. (25) Casida, M. E. Time-Dependent Density Functional Response Theory of Molecular Systems: Theory, Computational Methods, and Functionals; Elsevier: Amsterdam, 1996; Vol. Vol. 4. (26) Van Leeuwen, R. Int. J. Mod. Phys. B 2001, 15, 1969. (27) Runge, E.; Gross, E. Phys. Rev. Lett. 1984, 52, 997. (28) Jensen, L.; Autschbach, J.; Schatz, G. C. J. Chem. Phys. 2005, 122, 224115. (29) McLachlan, A.; Ball, M. Rev. Mod. Phys. 1964, 36, 844. (30) Maxwell, J. C. Philos. Trans. R. Soc. London 1865, 155, 459. (31) Furse, C.; Mathur, S.; Gandhi, O. IEEE Trans. Microwave Theory 1990, 38, 919. (32) Gersten, J.; Nitzan, A. J. Chem. Phys. 1980, 73, 3023. (33) Long, D. The Raman effect: a unified treatment of the theory of Raman scattering by molecules; John Wiley and Sons: Chichester, U.K., 2002. (34) Long, D. Raman Spectroscopy; McGraw-Hill: New York, 1977. (35) McMahon, J. M.; Wang, Y.; Sherry, L. J.; Van Duyne, R. P.; Marks, L. D.; Gray, S. K.; Schatz, G. C. J. Phys. Chem. C 2009, 113, 2731. (36) Johnson, P.; Christy, R. Phys. Rev. B 1972, 6, 4370. (37) Mullin, J. M.; Autschbach, J.; Schatz, G. C. Comput. Theor. Chem. 2011, DOI: 10.1016/j.comptc.2011.11.001. (38) Perdew, J.P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (39) Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Can. J. Chem. 1992, 70, 560. (40) Becke, A. D. Phys. Rev. A 1988, 38, 3098. (41) Perdew, J. P. Phys. Rev. B 1986, 33, 8822. (42) Jensen, L.; Schatz, G. C. J. Phys. Chem. A 2006, 110, 5973. (43) Krykunov, M.; Banerjee, A.; Ziegler, T.; Autschbach, J. J. Chem. Phys. 2005, 122, 074105. (44) Krykunov, M.; Autschbach, J. J. Chem. Phys. 2005, 123, 114103. (45) Bode, B. M.; Gordon, M. S. J. Mol. Graphics Modell. 1998, 16, 133. (46) Zhao, L. L.; Jensen, L.; Schatz, G. C. J. Am. Chem. Soc. 2006, 128, 2911. (47) Aikens, C. M.; Schatz, G. C. J. Phys. Chem. A 2006, 110, 13317. (48) Golab, J.; Sprague, J.; Carron, K.; Schatz, G.; Van Duyne, R. J. Chem. Phys. 1988, 88, 7942. (49) Yang, W.; Schatz, G. J. Chem. Phys. 1992, 97, 3831. (50) Michaels, A.; Nirmal, M.; Brus, L. J. Am. Chem. Soc. 1999, 121, 9932. (51) Shim, S.; Stuart, C. M.; Mathies, R. A. ChemPhysChem 2008, 9, 697.

computational effort over propagation in time when only one wavelength is of interest. Further, the highly scalable nature of the NWChem package allows TDDFT calculations with several hundred atoms or more if necessary. Although the applications that we considered used the zero Stokes shift approximation, the method is not limited to this approximation. Chen et al. explored the coupling of the dipole−quadrupole polarizability to the field gradient associated with the plasmonenhanced local field. However, this changed the SERS differential cross section by less than 1%; therefore, similar tests were not done here. The FDTD method was used in this work so that the differences in the current approach compared to that from Chen et al. would be isolated to the QM methods and not the underlying ED calculations. A better complement to the current work would involve using a frequency domain method such as the finite element method (FEM) or Mie theory for the electrodynamics calculations as then, the entire calculation can be done for a single wavelength. This would also allow a more accurate description of the sphere electrodynamics as FDTD must use a square grid to form a curved surface. Another important capability for future work is to include silver atoms explicitly in the QM calculations. There is no fundamental constraint on the number of incident light pulses or the direction of the propagation or polarization. Further work to extend this method to nonlinear optical (NLO) materials is limited by the underlying methods for the QM and ED calculations and not the approach outlined for the combination framework.

■

AUTHOR INFORMATION

Corresponding Author

* E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by DOE SISGR Grant DESC0001785. The authors would like to thank Prof. Jochen Autschbach, Dr. Hanning Chen, and Dr. Martin Blaber for useful discussions. This research was supported in part through the computational resources and staff contributions provided by Information Technology at Northwestern University as part of its shared cluster program, Quest.

■

REFERENCES

(1) Raman, C.; Krishnan, K. Nature 1928, 121, 501. (2) Jeanmaire, D. L.; Van Duyne, R. P. J. Electroanal. Chem. 1977, 84, 1. (3) Albrecht, M. G.; Creighton, J. A. J. Am. Chem. Soc. 1977, 99, 5215. (4) Jensen, L.; Aikens, C. M.; Schatz, G. C. Chem. Soc. Rev. 2008, 37, 1061. (5) Zhao, J.; Pinchuk, A. O.; Mcmahon, J. M.; Li, S.; Alisman, L. K.; Atkinson, A. L.; Schatz, G. C. Acc. Chem. Res. 2008, 41, 1710. (6) Camden, J. P.; Dieringer, J. A.; Wang, Y.; Masiello, D. J.; Marks, L. D.; Schatz, G. C.; Van Duyne, R. P. J. Am. Chem. Soc. 2008, 130, 12616. (7) Zeman, E.; Schatz, G. J. Phys. Chem. 1987, 91, 634. (8) Franzen, S. J. Phys. Chem. C 2009, 113, 5912. (9) Kerker, M.; Wang, D.; Chew, H. Appl. Opt. 1980, 19, 4159. (10) Ausman, L. K.; Schatz, G. C. J. Chem. Phys. 2009, 131, 084708. (11) Geshev, P.; Dickmann, K. J. Opt. A: Pure Appl. Opt. 2006, 8, S161. 1938


Combined Linear Response Quantum Mechanics and Classical

Recommend Documents