Real-Time Description of the Electronic Dynamics for a Molecule Close to a Plasmonic Nanoparticle.

The optical properties of molecules close to plasmonic nanostructures greatly differ from their isolated molecule counterparts. To theoretically investigate such systems from a quantum-chemistry perspective, one has to take into account that the plasmonic nanostructure (e.g., a metal nanoparticle-NP) is often too large to be treated atomistically. Therefore, a multiscale description, where the molecule is treated by an ab initio approach and the metal NP by a lower level description, is needed. Here we present an extension of one such multiscale model [Corni, S.; Tomasi, J. J. Chem. Phys.2001, 114, 3739], originally inspired by the polarizable continuum model, to a real-time description of the electronic dynamics of the molecule and of the NP. In particular, we adopt a time-dependent configuration interaction (TD CI) approach for the molecule, the metal NP is described as a continuous dielectric of complex shape characterized by a Drude-Lorentz dielectric function, and the molecule-NP electromagnetic coupling is treated by an equation-of-motion (EOM) extension of the quasi-static boundary element method (BEM). The model includes the effects of both the mutual molecule-NP time-dependent polarization and the modification of the probing electromagnetic field due to the plasmonic resonances of the NP. Finally, such an approach is applied to the investigation of the light absorption of a model chromophore, LiCN, in the presence of a metal-NP of complex shape.

Introduction

The response of a molecule to an applied electromagnetic field is strongly modified by the nearby presence of metal nanoparticles (NPs) or nanostructures. The resulting effects define the field of molecular plasmonics(1) and include many striking phenomena such as surface enhanced Raman scattering (SERS), surface enhanced infrared absorption (SEIRA), and metal enhanced fluorescence. Several theoretical models have been proposed through the years to understand such phenomena, some of which exploit a quantum-chemical description of the molecule.[2−4] For these models, the nearby nanostructure is typically too large to be routinely treated by quantum-chemistry (QM) methods as well. Presently, metal NPs of up to around 1000 atoms, i.e., a few nanometers in diameter, can be described by QM methods,[5,6] but they are much smaller than many NPs of practical interest. To tackle this limitation, and taking into account the good quality that a classical electromagnetic (EM) description gives of the optical properties of the NP, in the past different groups devised multiscale methods that keep the quantum-mechanical description of the molecule but include the NP as a classical object, either as a continuous dielectric[7−11] or as a collection of polarizable (and chargeable) atoms.[12−15]

Most of these models treat the applied electric field, and the corresponding molecule + NP response, in the frequency domain, i.e., assuming a monochromatic sinusoidal EM perturbation acting on the system. While this is, in principle, a very general approach (any time-dependent EM field can be Fourier-decomposed in a superposition of sinusoidal monochromatic fields), it is not the only one. A fully real-time description of optical phenomena involving molecules is also possible, i.e., a numerical approach that solves the time-dependent Schrödinger equation (or the time-dependent Kohn–Sham equation for time-dependent density functional theory, TDDFT, approaches) for the molecule in the presence of a time-dependent field. Such a real-time description is quite suitable when one wants to simulate the interaction of molecules with short light pulses (that span many frequencies) such as those employed in ultrafast spectroscopies, also in the presence of plasmonic nanostructures.[16−18] It is also a convenient approach for nonlinear optical property simulations.[19] Finally, a real-time description also offers an alternative way to interpret optical phenomena involving molecules, that may disclose more clearly the physics taking place in the system.[20] On the down side, a real-time description typically requires large computational efforts if fine features in a spectra should be solved, and many of the analysis that are straightforward in the frequency domain (e.g., orbital excitations contributing to a given transition; inspection of the transition density) require extra steps.

In this work we shall present an evolution of the polarizable continuum model (PCM)[21]-inspired approach that has been developed in the past to treat optical properties of molecules close to a metal NP[7] to the real-time domain. In that model, the molecule was treated initially with Hartree–Fock/time-dependent Hartree–Fock (HF/TDHF)[21] and then by DFT/TDDFT[22,23] and ZINDO,[24] and was coupled to a classical dielectric description of the metal NP employing an empirical frequency-dependent dielectric constant corrected for quantum-size effects.[7] The boundary element method (BEM) was used to solve the quasi-static electromagnetic problem needed to write the effective Hamiltonian for the molecule close to the NP. Many phenomena were treated with this model, such as SERS,[25] SEIRA,[26] metal-enhanced fluorescence,[23,24,27,28] and excitation energy transfer (EET) close to a NP.[28,29] Here the same model will be extended to simulations in real-time. The optical phenomena that will be considered here are the light absorption of the molecule–NP system, with the goal to extend the method to a range of phenomena in the near future.

In the field of quantum-chemical models for molecules close to NPs, a real-time approach has already been used in the past.[8,10,11] The approach that we are proposing here differs from previous ones in two main aspects. The first is the choice of the quantum-mechanical method for the molecule. All the cited works exploit real-time TDDFT (RT-TDDFT). Here, for the molecule we shall instead solve the time-dependent Schrödinger equation projected on a finite basis of multielectronic states, obtained by configuration interaction (TD CI approach). The second difference is that previous works used the finite difference time domain (FDTD) method for the NP polarization, while here we shall exploit BEM (in the quasi-static limit).

The choice of a wave function based method (CI) instead of RT-TDDFT is motivated by the qualitative artifacts that plague the current versions of RT-TDDFT when the molecular state departs substantially from the ground state. For example, the behavior in response to a light pulse that should take most of the population to a given excited state is unsound with respect to the pulse parameter,[30−32] and Rabi oscillations are not properly described.[33] The cause of these artifacts is believed to be the adiabatic approximation, i.e., the assumption that the exchange-correlation (xc) potential at a given time depends only on the density at that time rather than on all the previous instants.[34,35] Therefore, although RT-TDDFT will likely become in the future as accurate in the study of electronic dynamics as DFT and TDDFT are now for ground and excited state properties, at present, such artifacts make it difficult to correctly assign the physical origin of features observed in the electronic dynamics of a molecule close to a metal NP. TD CI offers the simplest theoretical framework to bypass these problems and has been already applied to the investigation of multielectronic dynamics for molecules in the gas phase.[36−39] Here, we shall limit ourselves to the simplest CI method, CI Singles (CIS). Although CIS is known to have a limited accuracy in terms of excitation energies,[40] it provides qualitatively correct behaviors also for charge transfer states, and can be systematically improved (e.g., the CIS(D) correction has already been used in a TD context).[36,37]

The classical electromagnetic description of the NP via BEM is an alternative approach to FDTD. BEM and FDTD are both widespread methods to study plasmonic nanosystems.[41] Here, BEM has some computational advantages when coupled to a molecule described with a localized basis set, that will be discussed in the theoretical section. BEM is most naturally a technique applied in the frequency domain, and it needs to be translated to the time domain to be directly coupled to TD CI. We have recently investigated the electronic dynamics of molecules in solution by a novel extension to the time domain of the PCM implicit solvation approach.[42,43] To this aim, we have developed equations of motion (EOMs) for the apparent charges that describe the evolution of time-dependent solvent polarization. Such EOMs are able to describe the delayed solvent response coded in the Debye dielectric constant. Formally, the electrostatic problem involved in PCM and in a quasi-static BEM description of a molecule close to a continuous metal NP is the same.[44] As such, the EOM method developed in ref (42) can be extended to NPs as well, as anticipated in that work. However, the Debye dielectric constant used there is not appropriate for metal NPs. Therefore, here we shall develop EOM suitable for the more complex Drude–Lorentz dielectric function,[45] and adopt a numerical approach consistent with such EOMs. The real-time TD CI coupled with EOM BEM evolution is implemented in the stand-alone homemade software WaveT that is interfaced with a local version of the widespread quantum-chemistry code GAMESS.[46,47]

The article is organized as follows: in section 2, we shall present the theoretical basis of the current approach, including the derivation of the EOM for the NP apparent charges that originates from the Drude–Lorentz dielectric function. Then, in section 3, we shall apply this theory to study the absorption spectrum of a molecule close to a complex-shaped metal NP. Finally, we shall draw the conclusions and discuss the development perspectives of this modeling approach.

Theory

The basic features of the model we are using for a molecule close to a NP have already been described elsewhere.[7,22,23] Here we summarize them briefly. The NP is described as a continuous dielectric with a complex shape and characterized, in the present work, by a Drude–Lorentz dielectric function ϵ(ω):[45]The metal NP is polarized both by the externally applied, incident time-dependent electric field E⃗inc(r⃗,t) and by the time-dependent EM field originated by the time-varying electronic charge density of the molecule. In turn, these two polarizations create time-dependent electric fields (sometimes dubbed reflected field E⃗ref(r⃗,t) and image or polarization field E⃗pol(r⃗,t), respectively)[48,49] that act on the molecule. In the quasi-static limit (i.e., in the limit of electric field wavelength much larger than the NP),[50] such electric fields are locally associated with time-dependent electrostatic potentials (Vref(r⃗,t) and Vpol(r⃗,t)) that satisfy proper quasi-static Poisson equations and the electrostatic boundary conditions at the NP surface.[7,51] Following standard integral equation theory and the boundary elements method (IEF-BEM), each of these electrostatic potentials (reflected and image) can be written as originated by a proper set of point charges (qref(t) and qpol(t)) located in representative points on the NP surface.[7]

The molecule is described by an effective Hamiltonian Ĥ(t) that includes the electromagnetic interaction with the incident field E⃗inc(r⃗,t), the reflected field E⃗ref(r⃗,t), and the polarization field E⃗pol(r⃗,t). The first is included in the usual dipole approximation; the second and the third interactions are instead written exploiting the charges qref(t) and qpol(t):Here, Ĥ0 is the Hamiltonian of the isolated molecule, is the dipole operator, E⃗inc(r⃗M,t) is evaluated in the molecular center of charge r⃗M, and V̂ is a vector operator representing the electrostatic potential of the solute at the representative points on the NP surface where the apparent charges qref(t) and qpol(t) are also located. Notably, qpol(t) depends on the electrostatic potential originated by the molecule V(t′) at all the previous instants t′ < t, giving rise to a nonlinear and nonlocal-in-time evolution problem.[42]

In the next subsection, we shall focus on the calculations of qref(t) and qpol(t) for a NP described by the Drude–Lorentz dielectric constant.

Equations of Motion for the Apparent Surface Charges for a Drude–Lorentz NP

Following previous works,[7,42] in the frequency domain, qpol(ω) and qref(ω) are calculated from the electrostatic potential produced by the molecular charge density on the NP discretized surface (V(ω)) and the potential associated with the incident field Vinc(r⃗) = −r⃗·E⃗inc(r⃗,ω), respectively, as follows:where Q(ω) is the response matrix in the frequency domain:A is a diagonal matrix with elements equal to the tessera areas; S and D are the representative matrices of the Calderon’s projectors:[21]s⃗ are the representative points on the nanoparticle surface and n⃗ are unit vectors normal to the j-th tessera and pointing outward. The generalization to the time domain of eq is straightforward (for the sake of brevity, we focus here on qpol and V only; the equations are identical for the pair qref and Vinc):with charges and potentials defined in the time domain and Q(t – t′) being the kernel (or memory) matrix that determines the nonequilibrium response of the NP in the time domain, obtained by Fourier-transforming Q(ω). Following the route presented in ref (42), we rewrite the Q(t – t′) matrix in terms of a diagonal representation:where the matrices T and Λ are built by the eigenvectors and the eigenvalues, respectively, of the symmetric matrix 1/2(S–1/2DAS1/2 + S1/2AD†S–1/2) ≈ S–1/2DAS1/2 (the identity is exact for the corresponding integral operators).

Taking into account the specific form of the Drude–Lorentz dielectric function eq , K(ω) in eq can be written asWe compute then the Fourier transformto obtain the following form of the diagonal kernel function:withΘ(t – t′) is the Heaviside step function. The equations for the apparent charges are then obtained by substituting in eq the expression of the diagonal kernel function K(t – t′) given in eq .

It is now possible to write an EOM for the charges qpol(t). The procedure consists in taking the first and the second derivatives w.r.t. time of eq , and use them to eliminate the convolution integral from the definition of qpol(t), similarly to what was done in ref (42). This leads to the EOM:and analogously for the reflected field charges:with K being the diagonal matrix defined before in eq and Kω still diagonal and with elements defined byThe form of eqs and 18 is clearly that of (coupled) forced damped harmonic oscillators (ẍ = −γẋ – ω2x + F). The squared frequencies of the oscillators are those given in eq and depend on two quantities: one is the Lorentz resonance frequency ω0 characteristic of the material making up the NP, the other (K) is a purely geometric term that is related to the shape of the NP and the various excitations it can sustain. We have already discussed the possible values of Λ for a spherical NP in ref (42). Based on those results, we find that for a spherical metal NP described by a Drude dielectric constant (i.e., ω0 = 0 and A = Ω2 in eq , where Ω is the plasma frequency), the squared frequencies of the plasmonic modes (i.e., Kω,) arei.e., the well-known frequencies of the multipolar plasmons of a spherical particle.

From a practical point of view, numerically integrating the EOM eq requires some care.[52] Here we exploit the velocity–Verlet scheme proposed in ref (53), second-order accurate in the numerical integration step dt.

Extension to the dielectric constant made of sums of Drude–Lorentz terms is straightforward, as they create a system of coupled oscillators.

Coupling to TD CI

The theory of TD CI in the presence of a solvent with a delayed response has been recently developed.[43] The equivalent TD CI theory for a molecule coupled to a metal NP is very similar, and will be only briefly summarized here.

The state vector |Ψ(t)⟩ of the molecule satisfies the time-dependent nonlinear Schrödinger equationApproximated solutions can be obtained by writing |Ψ(t)⟩ as a linear combination with time-dependent coefficients C(t) of a reference ground state |Φ0⟩ and a finite set of excited states |Φ⟩By using eq , the time-dependent Schrödinger equation, eq , becomeswhere the Hamiltonian matrix H has elements:We choose a reference state given by the Hartee–Fock determinant |Φ0⟩ = |HF⟩ for the molecule equilibrated with the NP described as a perfect (and charge neutral) conductor, and we use as the excited states |Φ⟩ (with I > 0) those obtained from a configuration interaction expansion limited to single excitations (CIS). Such CIS excited states |Φ⟩ and the corresponding energies are obtained here by solving, in the subspace of the single excited determinants |ψ⟩ = â†î|HF⟩ (↠and î are electron creation and annihilation operators, respectively), the time-independent Schrödinger equation for the molecular solute in the presence of the fixed Hartee–Fock polarization charges:With this choice, the elements of the Hamiltonian matrix in eq becomewhere Δqpol(t), μ⃗, and V are defined byqpol(t) in eq and qref(t) in eq are obtained by the EOM approach described in the previous section. Numerically, eq is solved by a simple Euler method at the first order,[37] using the same dt as for the apparent charge EOMs.

We close this section by remarking that, within the present BEM approach, the electrostatic potentials V in eq have to be calculated only for points on the NP surface. In a FDTD approach, they would be needed on a volume including the NP, making the calculation rather cumbersome.

Computational Details

The geometry of the LiCN molecule used in the calculations is the one of ref (36) (LiC and CN bond lengths of 1.949 and 1.147 Å, respectively). We have used a Drude dielectric function (i.e., ω0 = 0) with parameters for silver:[7] plasma frequency Ω = 0.332 au = 9.03 eV (A = Ω2) and γ = 1.515 × 10–3 au = 4.123 × 10–2 eV. No interband transition is included in the dielectric constant; therefore, the NP absorption spectra are not directly comparable to silver NP experimental ones. The calculations of the HF and CIS states in the presence of the NP, their energies (at frozen NP dielectric response), the expectation values and transition integrals of the dipole moment, and the electrostatic potential operators are performed using the 6-31G(d) basis set[54] exploiting a locally modified version of the GAMESS code.[46,47] Only the first 15 excited states are kept in the Hamiltonian. The coupled propagation of the wave function of the molecule, the polarization, and reflected field charges is performed using Wave-T, a homemade code. The time step of the propagation simulations is dt = 4.838 as (0.2 au). For the calculations with the spherical NP, 240 tesseras were used; in those with the elongated NP, the number of tesseras was 472 (tests with 676 and 850 tesseras were also done). We could perform 10.9 (7.4) ps/day of simulation on a Intel Xeon E5-2650 core and 472 (676) tesseras (no parallelization of the code was attempted).

Results and Discussion

Tests for the Polarization of a Spherical NP

To test our approach, we present the results of the simulations for the induced dipole of the NP in two cases where an analytical result is available, i.e. (i) a spherical NP subject to a monochromatic time-dependent electric field and (ii) a spherical NP interacting with a static dipole placed at a large distance.

(i) The analytical expression for the induced dipole moment of the NP (μ⃗N) generated by a monochromatic time-dependent electric field is the following:[55]where a is the radius of the NP, Einc is the complex amplitude of the incident field, ε⃗ is a unit vector defining its direction, and ϵ(ωinc) is the dielectric function evaluated at the frequency ωinc. In Figure we compare the profile of μ⃗N(t) obtained from the EOM BEM simulationwith the analytical one, eq , in the case of a sinusoidal external field with maximum amplitude of 10–5 au and a 2.5-nm-radius NP. s⃗ are the positions of the tessera representative points. During the simulation, the incident electric field is slowly turned on by modulating the sinusoidal form with a linear function between t = 0 fs and t = 120 fs. The frequency of the external field is ωinc = 0.26 au = 7.1 eV. The simulated profile and the analytical reference perfectly agree.

Figure 1

Simulation profiles of the x-component of μ⃗N for a spherical NP with a radius of 2.5 nm. The labels “sim” and “anl” refer to the results of the simulation and to the analytical profile, eq , respectively. The frequency of the incident field is ωinc = 0.26 au.

(ii) When a spherical NP is polarized by a static dipole at large distance, μ⃗N may be evaluated by taking the t → 0 limit in eq and considering, as an incident field, the electric field generated by the nonpolarizable dipole moment at the center of the NP. In a simulation we model such a limiting case by placing the molecule far away from the NP and by freezing its charge density as it was in vacuum. The value of |μ⃗N| computed for a system’s configuration where the center of the NP is placed at a distance of 25 nm from the center of charge of the molecule is 9.4184 × 10–3 D. The direction of the molecular dipole is perpendicular to the line joining the center of the molecule and the NP (the NP is the same used in the previous test). Taking into account that the magnitude of the dipole moment of LiCN in vacuum is 9.425 D, the difference between the calculated and the analytical values is less than 0.1%.

These simple tests back up the numerical accuracy of the model. In the next section we shall apply it to a physical problem.

LiCN Physisorbed on an Elongated NP: Absorption Spectra

In this section we discuss the results obtained from EOM BEM TD CIS simulations of a LiCN–NP system where the metal NP is build as the union of two intersecting spheres, each with radius 5 nm and having centers displaced 4 nm along the x direction. The LiCN molecule is also placed along x, with the N atom closest to the NP and at 4 Å from its surface. The surface of the NP is divided in tesseras that are smaller closer to the molecule (see Figure ). We compute the frequency-dependent linear polarizability tensors (αX) from the Fourier transforms of the induced dipole signals for the molecule (X = M), the NP (X = N), and the entire system (X = S). The elements of αX may be written as follows[56]here j and k refer to two of the three directions in Cartesian space, Einc,(ω) is the component of the Fourier transform of the perturbing field along k, μX,(t) = μX,(t) − μX,(0) is the component of the induced dipole profile along j, and τM is a damping constant that allows modeling a finite lifetime for the excited states (τM = 2000 au = 48.38 fs).[56] A narrow Gaussian pulse at t = 0, mimicking a δ-pulse, directed along the x axis is used to perturb the system within the linear regime. The perturbing field has a width of 48.38 as and a maximum intensity of 10–6 au. The imaginary part of αX,(ω) is proportional to the absorption cross section.[50] In Figure we compare such imaginary parts for the molecular polarizabilities in vacuum, and close to the NP. Two simulations have been performed to address the effect of the NP polarization on the molecular properties, one with the NP polarization frozen at its equilibrium value before the action of the perturbation (Δqpol(t) = qref(t) = 0; “NP frozen” in Figure ), and one with the full propagation of both the NP and molecule electronic dynamics as detailed in section 2 (“full” in Figure ). The effect of the equilibrium polarization of the NP manifests as a blue-shift on the molecule excitation energies in the gas phase. This is reasonable, since the dipole moments in the excited states corresponding to these excitations are lower than in the ground state, and therefore, they are less stabilized by the NP presence. In addition to this, accounting for the full time-dependent NP polarization results in a further, but smaller, blue-shift of the molecular excitations and a strong enhancement of the peak intensities (see the inset in Figure ). The new peaks at a frequency lower than 7 eV are reminiscent of the NP absorption. They can be negative, as the physical requirement of a positive absorption (starting from the ground state) refers to the overall (molecule + NP) system. In Figure we report the imaginary parts of the NP and system polarizabilities for the simulation with the full propagation of the system degrees of freedom. Because of the difference in size between the NP and the molecule, the responses have different orders of magnitude and the effect of the molecule on the NP response is not detectable.

Figure 2

Side (left) and front (right) views of the NP used in section 3.2. The long axis of the NP is directed along x. The NP surface is divided in 472 tesseras, refined in the region closer to the molecule.

Figure 3

Imaginary part of αM, plotted as a function of frequency for (i) a simulation of the molecule in vacuum (vacuum), (ii) a simulation with the NP polarization frozen at its equilibrium value when no perturbation is acting on the system (NP frozen), and (iii) a simulation with the full propagation of the molecule and NP degrees of freedom (full—see section 2). Inset: a zoom in the 8.1–8.6 eV region.

Figure 4

Imaginary parts of αS, (S) and αN, (N) plotted as a function of frequency for a simulation where the full propagation of a molecule and NP degrees of freedom is performed (see section 2).

Finally in Figure we compare the profiles of the imaginary part of αN, for the nanoparticle in vacuo calculated (i) using eq and (ii) directly in the frequency domain. A good agreement is found between the two methods. We close this section by remarking that the results, in general, will depend on the accuracy of the nanoparticle tessellation. This is discussed in the Supporting Information, where it is shown that increasing the number of tesseras from 472 to 676 can shift the nanoparticle peaks by 0.1–0.2 eV and that peak intensities can also change. In practical applications where specific experiments are targeted, it is therefore necessary to check the convergence of the numerical results with the tesselation.

Figure 5

Imaginary part of αN, as a function of the frequency calculated using eq (time) and directly in the frequency domain (freq).

Conclusions and Perspectives

In this work, we have presented a real-time methodology to study the electronic dynamics of a molecule close to a metal NP, based on a polyelectronic (TD CI) description of the molecule and a EOM BEM treatment of the molecule–metal NP electromagnetic coupling (in the quasi-static limit). We have then applied this approach to study the light absorption by a Drude metal NP plus molecule system. The present work represents an initial step to provide a comprehensive modeling of the real-time dynamics of molecules close to metal NPs, that will require developments in various directions. For what regards the molecule, the quality of the wave function approach used here can be improved by obtaining the necessary quantities (excitation energies, dipole and electrostatic potential matrix elements) via methods more accurate than CIS, that can be coded in richly featured and flexible quantum-chemistry software such as GAMESS. The coupling with RT-TDDFT, notwithstanding its current drawbacks related to the adiabatic approximation, is also important to benefit from the advantageous accuracy vs computational cost balance of such a method. Moreover, here, the molecular nuclei have been frozen, which is suitable if only the electronic dynamics is of interest. The study of several other phenomena, ranging from vibrational spectroscopy to photochemistry, requires such an assumption to be relaxed.

On the side of the NP description, there are some extensions that are relatively straightforward, such as to consider more complex dielectric functions that involve a sum of Drude–Lorentz terms. With this kind of ϵ(ω), one can reproduce quite faithfully empirical permittivities that also account for interband transitions (disregarded here), as done in FDTD approaches. Corrections for the limited mean free-path of the electrons in the confining NP can be included. Nonlocal aspects of the metal–NP response (that are included, for example, in proper nonlocal dielectric functions[57]) may also be included at different levels in a IEF BEM approach.[22,58] In perspective, a supermolecular description of the molecule + NP system may also be envisaged, where part of the NP is treated quantum-mechanically together with the molecule.[10]

Here we have exploited the quasi-static approximation, that becomes less and less accurate for increasing NP size.[59] BEM can be extended to a full electromagnetic description as well,[60,61] that is still compatible with coupling with a quantum-chemistry molecule.[51] Finally, we have neglected here the possible presence of a solvent/dielectric matrix, that can be included as an additional dielectric (with its own time-dependent response) hosting the molecule in a cavity, as done previously in the frequency domain.[7]

In conclusion, the development presented here represents a promising starting point to investigate the optical properties of molecules close to metal NPs, that is amenable to various developments currently ongoing in our group.