Molecular Dynamics Study of Plasmon-Mediated Chemical Transformations

Heterogeneous catalysis of adsorbates on metallic surfaces mediated by plasmon has potential high photoelectric conversion efficiency and controllable reaction selectivity. Theoretical modeling of dynamical reaction processes provides in-depth analyses complementing experimental investigations. Especially for plasmon-mediated chemical transformations, light absorption, photoelectric conversion, electron-electron scattering, and electron-phonon coupling occur simultaneously at different timescales, rendering it very challenging to delineate the complex interplay of different factors. In this work, a trajectory surface hopping non-adiabatic molecular dynamics method is used to investigate the dynamics of plasmon excitation in an Au$_{20}$-CO system, including hot carrier generation, plasmon energy relaxation, and CO activation induced by electron-vibration coupling. The electronic properties indicate that when Au$_{20}$-CO is excited, a partial charge transfer takes place from Au$_{20}$ to CO. On the other hand, the dynamical simulations show that hot carriers generated after plasmon excitation transfer back and forth between Au$_{20}$ and CO. Meanwhile, the C-O stretching mode is activated due to the non-adiabatic couplings. The efficiency of plasmon-mediated transformation ($\sim$40\%) is obtained based on the ensemble average of these quantities. Our simulations provide important dynamical and atomistic insights into plasmon-mediated chemical transformation from the perspective of non-adiabatic simulations.


Introduction
Plasmon-mediated chemical transformations, which utilize plasmonic nanostructures as catalysts, is an emerging technology that has attracted extensive attentions. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] Here, energy of light is collected by the collective oscillations of electrons, resulting in a localized plasmon resonance (LSP) which has scattering cross section much larger than relevant geometric sizes. [16][17][18] After being excited, LSP decays simultaneously through radiative and nonradiative pathways. Usually, the majority of the energy stored in the plasmonic field is dissipated through nonradiative decay, resulting in the formation of energetic (hot) electrons (HEs) and holes. 19 Surface scattering, electron-phonon coupling, and electron-electron scattering are the major mechanisms that determine the distribution of these carriers. [20][21][22][23] The energies of HEs (and their corresponding holes) are then redistributed via electron-electron scattering, released to phonon modes via electron-phonon interactions, and ultimately dissipated into environment via thermal conduction. 24,25 Thanks to their tunable optical properties, plasmons have been extensively explored in the context of driving chemical transformations via different mechanisms. So far, four major mechanisms have been identified in these reactions, 8 including 1) plasmon-enhanced intramolecular excitations (or near-field effect), where the plasmon resonance overlaps with the electronic transition energies, [26][27][28][29] 2) indirect HE transfer from metal nanoparticles (NPs) to the adsorbed molecules, where HE generation within the nanostructures follows by a charge transfer process to the adsorbed molecules, [30][31][32] 3) direct charge transfer mechanism, where the electrons are directly excited from the valence band of plasmonic materials to the virtual molecular orbital (MO) of the adsorbates, [33][34][35][36] and 4) thermal activation (or local heating effect) due to the HE relaxation. [37][38][39] In particular, HE-mediated reactions have higher tunability compared to conventional temperature-driven catalysis because the energy can be selectively deposited into particular reaction coordinates in the former process.
Hence, HE-mediated chemical transformation of adsorbates on the surface of nanoparticles induced by the unique light-matter interactions of plasmonic NPs has attracted significant attention in recent decades. For instance, plasmonic excitations of Ag and Au NPs have shown to activate carbonyl hydrogenation, dissociation of H 2 and reduction of nitroaromatics, etc. 32,[40][41][42][43][44] In addition, for HE-mediated processes, recent reports have suggested that interactions with semiconductors or adsorbates enhance the reaction efficiency through interfacial states. 45,46 These interfacial states provide an additional energy dissipation pathway for the plasmons or serve as transient reservoirs of these hot carriers. 2,3,45 Other studies have investigated the tunability of energy distribution of hot carriers to enhance quantum yield and also to realize controllability of reaction pathways. For example, Manjavacas and co-workers found that small-size NPs and longer carrier lifetimes result in higher hot carrier energies but lower hot-carrier production rates, and vice versa. 19 Numerous studies have reported ways of improving HE production in plasmonic nanomaterials, including the use of small NPs with large surface-to-volume ratio, designing plasmonic materials with a longer carrier mean free path, constructing hybrid nanostructures with plasmon resonances in the red and infrared regions, and designing NPs with extended plasmonic hot spots. [22][23][24]47 These findings have inspired further interests in the field of plasmonic catalysis as it offers opportunities to develop new and selective catalytic processes that were previously inaccessible.
In addition to HE-mediated reactions, heating effect due to the HE relaxation may activate chemical reactions as well. Recent debates on the thermal impact underpin the necessity of atomistic and dynamical insights via simulations of the HE generation, transfer, and relaxation processes on equal footing. This demands for special theoretical tools such as our recent development of NEXMD-DFTB method. 48 Dynamic simulations of plasmon-enhanced catalysis were previously performed using different theoretical methods. For example, Meng and co-workers have investigated H 2 O splitting catalyzed by Au 20 cluster using real-time timedependent density functional theory (RT-TDDFT) and Ehrenfest dynamics. 49 Their results show significant energy and spatial overlap between oscillating electron density within the Au cluster and MOs of H 2 O. Our previous work demonstrated that HEs transfer to the antibonding state of the H 2 molecule from the metal NP upon photoexcitation, a process that introduces a repulsive potential and drives the cleavage of the H-H bond. 50 These calculations assume high-intensity photoexcitation, which depends on the strength of the external field.
Recently, we also proposed a mechanism for H 2 dissociation, which suggests that the charge transition from HE states within the metal NP to charge transfer (CT) excited states is the main channel to trigger H 2 reaction. 30 Importantly, this study only involved one electron (one-hole) or mono-electronic excitation process, which facilitates the understanding of the basic principle of such reactions. Furthermore, this investigation explained the energy and electron transfer in terms of reaction coordinates, which elucidates the details of reaction pathway. However, the simple model system in this study does not support plasmon excitation, and therefore precludes the description of processes including plasmon decay into hot carriers, electron transfer to adsorbed H 2 , and hot carriers decay through electron-phonon interaction.
In this work, an Au 20 cluster is chosen, which supports plasmon-like excitation with a symmetric geometry and high stability. 51 A CO molecule is adsorbed on the Au 20 which gives rise to CT states. 52 For this system, the coupling strength between Au 20 and CO is in the intermediate regime, so both direct and indirect mechanisms can be studied simultaneously, together with their competition with HE relaxation. Through molecular dynamics simulations, it is find that the C-O stretching vibration is excited by HE transfer via CT excited states. These states are excited via either direct or indirect HE transfer processes.
Moreover, our results reveal that both direct and indirect processes are significantly faster than the energy relaxation process in Au 20 , which demonstrates that the HE transfer is the dominating process that governs the chemical transformation of the CO molecule.

Methods
Density Functional Tight-Binding Theory (DFTB) Calculations. Geometry optimizations and ground-state properties are carried out with the DFTB+ code. 53 The auorg-1-1 parameter set, designed to describe optical excitations of organic molecules on gold nanoclusters, are employed for all computations in this work. 54 To obtain a correct desorption reaction coordinate, the DFTB parameters associated with the repulsive potential between . .} and {a, b, · · · } denote the occupied and virtual KS orbitals, respectively. K is the exchange-correlation kernel. Within the Mulliken approximation and LR-TDDFTB formalism, the kernel is simplified by em- , S is the overlap matrix of atomic orbitals, {µ, ν, · · · } denote atomic orbitals, and γ includes the Coulomb and gradient-corrected exchange-correlation kernels. The element F ia denotes the contribution coefficient for KS orbital transition from an occupied orbital i to an virtual orbital a. Hence, for the Au 20 −CO system, 8395 possible transitions are included for each excited state. Because plasmon excitation usually corresponds to higher energy excited state, hundreds of excited states are generally required. To reduce the computational cost, only occupied-virtual MO pairs with energy differences smaller than 3.85 eV are included in the excited state calculations.

Trajectory Surface Hopping (TSH) non-adiabatic Molecular Dynamics (NAMD)
Simulations. The TSH driver within the LR-TDDFTB framework is used to simulate the NAMD of Au 20 −CO initiated by plasmon excitation. This is recently implemented by combining the open-source software packages DFTB+ and NEXMD. 48 The adiabatic electronic structures, including ground, excited state energies, gradients, excited states transition densities, and non-adiabatic couplings (NACs) between excited states, are described at the TDDFTB level. The non-adiabatic effects originating from coupled motions of electrons and nuclei are treated by the TSH algorithm. 48,57,58 First, the system is optimized on its singlet ground state at the DFTB level. This is followed by a 90 ps Born-Oppenheimer ground state MD simulation with a time step ∆t = 0.1 fs. Here, a Langevin thermostat with a damping rate of 20 ps −1 is used to keep the temperature fluctuating around 300 K. 59 After a 10 ps equilibration period, 300 snapshots of initial geometries together with nuclear velocities are sampled from the ground state trajectory with 166 fs intervals, which are used to calculate the average absorption spectrum and serve as starting points for subsequent excited state dynamics. To obtain the average absorption spectrum, vertical excitation energies and oscillator strengths of all snapshots are computed using the LR-TDDFTB. 54 For NAMD simulations, all excited state trajectories are prepared by plasmon excitation (i.e., 2.7 eV) according to the absorption spectrum and are propagated for 1 ps. The timesteps for nuclear (classical) and electronic (quantum) equations of motions are set to 0.4 fs and 0.1 fs, respectively. 160 excited states are taken into account (N st = 160) in our simulations to ensure the inclusion of plasmon excitation and allow for possible transitions to higher excited states. As described in previous work, the instantaneous decoherence corrections (IDC) and trivial unavoided crossings methods are adopted to improve the accuracy of calculations. 59 The quantum timestep is further refined by a factor of 10 when potential trivial crossings are detected.

Results and Discussions
Geometry and Electronic Properties of Au 20 −CO Different conformations of adsorbates on the NP surface strongly influence the electronic coupling, ultimately determine the mechanism of hot carrier transfer and plasmon decay. 2,60,61 The inset of Figure 1 shows the ground state equilibrium geometry of Au 20 −CO optimized by DFTB method. 54 This structure with CO molecule adsorbed at the apex site of the Au 20 cluster is consistent with previous findings. 62 To further confirm this adsorption configuration is favorable during the dynamic process, ground state potential energy surface (PES) along the adsorption reaction coordinate (Figure 1) is calculated. As compared with DFT (BP86 63 ) result, it can be clearly seen that the auorg-1-1 parameter set 54 significantly underestimates the energy barrier of CO desorption. Therefore, in order to obtain a correct topology of PES, we reoptimize the repulsive potential between atom C (O) and Au in the auorg-1-1 parameter set. 54,64 As shown in Figure 1, the potential profile obtained using the optimized parameter is consistent with that of DFT (BP86 63 and LC-ωPBE 65  permits an analysis of the collective character of plasmon excitation. As shown in Table   S1, the excitation energy of 2.71 eV composes of multiple KS transitions with comparatively similar weights. The energy distributions of holes and electrons resulting from these multiconfigurational transitions contain only a few peaks and have different shapes owing to the quantum confinement effect and discretized energy levels ( Figure S1). For Au 20 -CO, the adsorption of CO induces a broadening and slightly red shift of the major absorption peak (black curve in Figure 2(a)), which are consistent with previous studies . 60,66 New peaks also emerge in the lower energy range of the absorption spectrum of Au 20 -CO, which are absent in the absorption spectra of bare Au 20 and CO molecule. This signifies that the presence of CT excited states upon CO adsorption. Note that CT excited states defined here refer to excited states for which at least 0.1e is transferred from the Au 20 cluster to the CO molecule upon excitation (Table S4). In combination with the analysis of transition component of the first excited state (S 1 , Table S3) and the spatial distributions of MOs (Figure 2 It is well-known that TDDFT calculations with local and semilocal functionals lead to spurious low-lying CT states. The parameterization in DFTB is however generally based on PBE semilocal functional which inherits this error in the description of CT states. Figure   S2 shows the absorption spectra of Au 20 calculated with TDDFTB, TDDFT(PBE), 68 and TDDFT(LC-ωPBE). 65 As expected, the TDDFTB and TDDFT(PBE) spectra are densely populated with low-energy excited states as compared with that of long-range corrected LC-ωPBE calculation. However, the sparse manifold of excited states calculated with LC-ωPBE functional does not fit the characteristics of plasmon excitations. On the other hand, DFTB/TDDFTB with GGA functional are more suitable for qualitative description of the dense manifold of excited states in plasmon nanoparticles and the competition between HE relaxation and HE transfer-induced chemical reactions. Therefore, in this work, we study the plasmon-like excitation dynamics of the Au 20 systems at the level of TDDFTB. Nevertheless, for a more accurate description of excited state properties, long-range corrected DFTB can be applied which exhibits similar accuracy as range-separated DFT methods at significantly reduced computational cost. 69,70 To summarize, the excited state properties of the Au 20 -CO system indicate the existence of two qualitatively different manifolds of states. One manifold of excited states is dominated by the electronic excitation confined in the Au 20 cluster alone, which can be likened to the HE state. The other manifold involves the excitations to the delocalized hybrid Au 20 −CO MOs, which is likened to CT states. 30 Since HE and CT states can cross, indirect HE transfer from Au 2 0 cluster to CO molecule can occur via the non-adiabatic transitions between HE and CT states. 30 Hence, both direct and indirect HEs transfer may contribute to the relaxation of the plasmon excitation and activation of CO vibrations. At the heart of both direct and indirect HEs transfer mechanisms is the coupling between CT states and vibrational states of Au 20 -CO complex, which can lead HEs to dispose a portion of their energy into the vibrational motion of CO during the relaxation process. To verify this scenario and to clarify the dynamical details, including the direct/indirect HE transfer, the non-adiabatic relaxation, and the competition between the HE transfer and energy relaxation, we next perform direct NAMD simulations of plasmon excitation.

Real-time Simulations of Plasmon Excitation of Au 20 -CO
The processes of plasmon relaxation in the Au 20 -CO system are summarized in Figure 3.
The NAMD simulations allow us to delineate different processes during plasmon relaxation, as shown schematically in Figure 3(a). To analyze the energy evolution during the relaxation process, a time-energy 2D map is constructed by calculating the probability of finding the system at a given energy level and time. Specifically, the probability is defined as where N E (t) and N tot stand for the number of trajectories with excitation energy E and the total number of trajectories, respectively. Such representation is similar to that often used in experiments, 71 and thus provides a direct comparison between simulation and experimental data. Figure 3 at each time step, the lifetime of energy relaxation can be obtained by fitting the data to an exponential decay function where τ is the energy relaxation time. A and B are determined according to the conditions: , and E 1 stands for the energy of the lowest excited state (S 1 ).
From Figure 3(c), the lifetime of energy relaxation process in Au 20 -CO is estimated to be ∼ 1.0 ps, which is significantly faster than that of bare Au 20 (∼2.7 ps). Here, non-adiabatic coupling between ground and excited states is not considered. Thus, the system is allowed to relax to S 1 only. This result indicates that hybridization of the MOs of CO molecule and Au 20 accelerates the energy relaxation process by offering an additional pathway for plasmon energy decay. It is postulated that the CO stretching mode participates in the relaxation process via the CT excited states, which accelerates the energy relaxation.
Apart from ensemble average of all the trajectories in the dynamics, each trajectory in our simulations give insights into different possible reaction pathways. Here, we analyze a typical trajectory to shed light on the details of plasmon relaxation. Figure 4(a-d) shows the evolution of the potential energy, excited state occupations, number of HEs on the CO molecule, and the bond length of CO. The situation when the potential energy curve shows vibrations with high frequency and small amplitude is highlighted with green shadow. The appearance of this vibration is found to be associated with non-adiabatic hops between excited states, as shown in Figure 4(b). The hops between excited states is accompanied with electron transfer to CO molecule (Figure 4(c)) and also an elongation of C-O bond ( Figure 4(d)). Such a strong concurrent behavior across all these quantities, especially the that the bond length of CO oscillates around 1.12 Å, and its oscillation amplitude gradually increases. We then estimate the ionic temperature of CO and Au 20 according to the equipartition theorem (E kin = 3 2 N k B T , k B is the Boltzmann constant, E kin is the kinetic energy). As shown in Figure 5(c), the temperature increase of CO is significantly faster than that of Au 20 , and the final temperature of CO is 4 times higher than that of Au 20 after the 1 ps simulation. This steep temperature gradient strongly indicates that activation of CO originates from HE transfer-induced vibrational excitations rather than thermal effect. According to Figure 5(b), we define CO to be activated if the amplitude of CO vibrational motion is larger than 1.2 Å, and the corresponding probability of CO activation is shown in Figure 5(d). The evolution of HEs on CO shows a similar trend with the increase of CO bond length, and the probability of activation reaches ∼40% at 1 ps, demonstrating that the HEs transfer can effectively activate CO by exciting the C-O stretching mode. Finally, we estimate the timescale of HEs transfer from Au 20 to CO by calculating the accumulated probability that HEs on CO is greater than 0.1e at each time step. However, the HEs on CO oscillate with time, indicating back-and-forth CT between Au 20 and CO.
In order to analyze the accumulated probability, we think that for each trajectory, once the HEs on CO is greater than 0.1e at a certain time, the HEs on CO at all subsequent time steps is set to 0.1 e in the statistics. Here, the reason we choose 0.1e as the target value for determining whether HEs transfer occurs is that the initial plasmon excitation has ∼ 0.1e charge transfer from Au 20 to CO. As shown in Figure 5(d) and Figure S4, initially, a majority of HEs are distributed on Au 20 , and the probability of direct HEs transfer is only ∼7%. However, HEs rapidly transfer to CO (0.1e) within 100 fs. Such a fast process is also consistent with previous experimental investigations. 72 Importantly, the energy relaxation time (∼1 ps) is significantly longer than the HE transfer timescale (∼100 fs), indicating that HE transfer occurs before electronic energy relaxation completes. Therefore, it is expected that the plasmon excitation in Au 20 can effectively activate CO adsorbate and mediate subsequent chemical reactions before the nonradiative decay of HEs.

Conclusion
The plasmon-induced bond activation of CO adsorbed on Au 20 cluster triggered by the HEs transfer has been investigated by a non-adiabatic molecular dynamics simulation combining LR-TDDFTB theory and TSH method. The simulations naturally treat the entire dynamical processes and multiple interactions in plasmon-mediated chemistry on an equal footing, including the plasmon excitation, HEs relaxation, direct/indirect HE transfer, and the activation of CO vibration mode induced by HE transfers. These dynamical processes compete with the vibrations scattering process inside Au 20 that results in local heating of the system.
Our simulations provide a comprehensive picture of plasmon-induced chemistry in atomic details.
The simulations demonstrate that CT states are crucial in plasmon-induced CO activation. The direct and indirect HEs transfers are generated by excitation of the CT states, leading to the activation of CO stretching mode, which is also an essential feature of the plasmon energy relaxation process. This activation of CO vibrational mode has been observed with high efficiency (∼ 40%). Importantly, the HEs transfer is faster than the conventional scattering process of the Au 20 . The simulations show that HEs transfer from Au 20 to CO completes within ∼100 fs, while the energy relaxation occurs on a ∼1 ps timescale.
The present simulation employs the LR-TDDFTB-based TSH method to simulate non-adiabatic molecular dynamics following the plasmon excitation directly. At each time step, 160 excited states are calculated based on the LR-TDDFTB method, pushing the limits of currently feasible theoretical methods. Nevertheless, our direct atomistic simulations outline the dynamics of the plasmon-mediated chemical transformations in terms of the evolution of the potential energy with the reaction coordinates. This enables a straightforward demonstration of the energy relaxation and HE transfer during the reaction and elucidates the details of the reaction pathways. Our numerical simulations give clear dynamical insights into CO bond activation on Au clusters. We believe the study reported in this manuscript paves the way for simulating and understanding plasmon-mediated photochemistry.