{'fraction_non_alphanumeric': 0.03352194787379972, 'fraction_numerical': 0.011445473251028807, 'mean_word_length': 4.7025177218284036, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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][3][4] For these models, the nearby nanostructure is typically too large to be routinely treated by Quantum Chem-istry (QM) methods as well. Presently, metal NPs of up to around 1000 atoms, i.e., a few nm 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][8][9][10][11] or as a collection of polarizable (and chargeable) atoms. [12][13][14][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 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][17][18] It is also a convenient approach for non-linear optical properties 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, 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) requires 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 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 is 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 been already 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 quasistatic limit).
The choice of a wavefunction 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][31][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 to 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 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 multi-electronic dynamics for molecules in gas-phase. [36][37][38][39] Here, we shall limit ourselves to the simplest CI method, CI Single (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 been already 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 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 Polarizable Continuum Model (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 quasistatic 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 home-made 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 has been already described elsewhere. 7,22,23 Here we summarize them shortly. The NP is described as a continuous dielectric with a complex shape and characterized, in the present work, by a Drude-Lorentz dielec- The metal NP is polarized by both 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, such 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 quasistatic limit (i.e., in the limit of electric field wavelength much larger than the NP) 50  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 interac-tions are instead written exploiting the charges q ref (t) and q pol (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 andV is a vector operator representing the electrostatic potential of the solute at the representative points on the NP surface where the apparent charges q ref (t) and q pol (t) are also located. Notably, q pol (t) depend on the electrostatic potential originated by the molecule V (t ) at all the previous instants t < t, giving rise to a non-linear and non-localin-time evolution problem. 42 In the next subsection, we shall focus on the calculations of q ref (t) and q pol (t) for a NP described by the Drude-Lorentz dielectric constant.

NP
Following previous works, 7,42 in the frequency domain q pol (ω) and q ref (ω) are calculated from the electrostatic potential produced by the molecular charge density on the NP discretized surface (V(ω)) and the potential associated to the incident field V inc ( r) = − r · E inc ( r M , ω), respectively, as follow: 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 with charges and potentials defined in the time domain and Q(t −t ) being the kernel (or memory) matrix that determines the non-equilibrium 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 build by the eigenvectors and the eigenvalues, respectively, of the symmetric matrix 1/2 S −1/2 DAS 1/2 + S 1/2 AD † S −1/2 ≈ S −1/2 DAS 1/2 (the identity is exact for the corresponding integral operators).
Taking into account the specific form of the Drude-Lorentz dielectric function Eq.(1), K ii (ω) in Eq. (8) can be written as: We compute then the Fourier transform to obtain the following form of the diagonal kernel function: Θ(t − t ) is the Heaviside step function. The equations for the apparent charges are then obtained by substituting in Eq.(7) the expression of the diagonal kernel function K ii (t − t ) given in Eq. (12).
It is now possible to write an EOM for the charges q pol (t). The procedure consists in taking the first and the second derivatives w.r.t. time of Eq. (7), and use them to eliminate the convolution integral from the definition of q pol (t), similarly to what done in ref. 42. This leads to the EOM: and analogously for the reflected field charges: with K f being the diagonal matrix defined before in Eq.(10) and K ω still diagonal and with elements defined by: The form of Eqs. (15) and (18) is clearly that of (coupled) forced damped harmonic oscillators i.e., the well-known frequencies of the multipolar plasmons of a spherical particle.
From a practical point of view, numerically integrating the EOM Eq.(15) 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 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 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 non-linear Schrödinger Approximated solutions can be obtained by writing |Ψ(t) as a linear combination with time- where 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 |Φ I (with I > 0) those obtained from a configuration interaction expansion where ∆q pol (t), µ IJ and V IJ are defined by q pol (t) in Eq. (27)

Computational Details
The

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, E inc 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 1 we compare the profile of µ N (t) obtained from the EOM BEM simulation µ sim with the analytical one, Eq. (30), in the case of a sinusoidal external field with maximum amplitude of 10 −5 au and a 2.5nm-radius NP. s i 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 = 0fs and t = 120fs. The frequency of the external field is ω inc = 0.26au=7.1eV. 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.5nm. The labels "sim" and "anl" refer to the results of the simulation and to the analytical profile, Eq. (30), respectively. The frequency of the incident field is ω inc = 0.26au.
(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. (30) and considering, as an incident field, the electric field generated by the non polarizable dipole moment at the center of the NP. In a simulation we model such a limit 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 25nm 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.425D, the difference between the calculated and the analytical value 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 sphere, each with radius 5nm and having centers displaced of 4nm 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 Fig.2). We compute the frequency-dependent the NP manifests as a blue shift on the NPs excitation energies in 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   3). The new peaks at a frequency lower than 7eV 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 4 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. Finally in Figure 5 we compare the profiles of the imaginary part of α N,xx for the nanoparticle in vacuum calculated (i) using Eq. (32) 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.

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 wavefunction 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 method. Moreover, here the molecular nuclei have been frozen, that is suitable if only the electronic dynamics is of interest. The study of several other phenomena ranging from vibrational spectroscopy to photochemistry, requires such 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. Non-local aspects of the metal-NP response (that are included for example in proper non local dielectric functions 57 ) may be also 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 on-going in our group.