Vibron-vibron coupling from ab initio molecular dynamics simulations of a silicon cluster

We study the temperature dependent dynamical processes of a Si10H16 cluster and obtain a blue shift of the Si-Si vibrational modes with transverse acoustic character and a red shift of the other vibrational modes with increasing temperature. We link this behavior to the bond length expansion and the varying sign of the Grueneisen parameter. We further present a computational approach able to extract the vibron-vibron coupling strength in clusters or molecules. Our approach is based on ab initio Born-Oppenheimer molecular dynamics and a projection formalism able to deliver the individual vibron occupation numbers. From the Fourier transform of the vibron energy autocorrelation function we obtain the coupling strength of each vibron to the most strongly coupled vibronic states. We find vibron-vibron coupling strength up to 2.5 THz with a moderate increase of about 5 % when increasing the temperature from 50 to 150 K.


Introduction
Colloidal semiconductor nanoclusters (NCs) have undergone an enormous development in the fields of optoelectronics, spintronics, and bio-labeling over the past two decades. The good control of the NC size makes it possible to tailor their electronic and optical properties. Many successful applications in these fields were already reported [1,2,3,4,5,6,7]. The experiments are usually performed at room temperature, making the treatment of vibrational properties, especial the anharmonic effects, crucial [8,9,10]. A solid understanding of the anharmonic effects of vibrations, is therefore a decisive step for the real world application of nanostructures, where the physical properties such as thermal conductivity in nanowires [11], non-radiative transition via phonon in NCs [12], and Raman spectra broadening in nanostructures [13] are dominated by the phonon lifetime. There are basically two approaches to address this problem theoretically.
The first one is to calculate the electron-phonon coupling [14,15,16,17] and phonon-phonon coupling [18,19,20,21,22] terms via perturbative approaches. The temperature effects are subsequently included using the Bose-Einstein distribution of the lattice vibrations (phonons). The physical properties, such as spectral broadening, spinflip, loss of quantum coherence, relaxation of charge carriers, and thermal conductivity have been studied theoretically [19,20,16,17,24], mostly considering harmonic or quasiharmonic (harmonic approximation performed at different volumes) approximations. Some have considered the third-order and fourth-order derivatives of the total energy [21,22], but the large computational demand for these higher-order derivatives limit the study to very small nanostructures [25] or bulk materials [21,22]. The 2n + 1 theorem, on the other hand limits the studies based on density functional perturbation theory (DFPT) to third order process [26].
The second one is to use molecular dynamics simulations [27,28,29,30], where the temperature effects including all anharmonic effects are intrinsically contained in the atomic trajectories of the simulations. The temperature-dependent vibrational density of states (DOS) and the thermal conductivity of nanostructures have already been studied using classical molecular dynamics simulations [48]. The required empirical force field potentials are limited by the lack of transferability to different systems and by the inability to correctly predict chemical bonding processes. These difficulties are overcome in the case of ab initio molecular dynamics (AIMD) simulations [28,29,30,31,32,20,33], which does not build upon the harmonic approximation but allows for a self-consistent rearrangement of the changes, including all anharmonic effects. Based on the accurately calculated forces from the electronic structure calculations, AIMD enables us to monitor the changes in the electronic and vibrational properties on-the-fly, with thermal effects included directly. Although AIMD simulations have been successfully applied to study a variety of problems [34], and the phonon lifetime in bulk has been studied using classical molecular dynamics simulations [35,36,48], a vibron-vibron interaction extracted from AIMD simulation is still lacking.
In this work, we perform AIMD simulations within the Born-Oppenheimer (BO) approximation of a Si 10 H 16 cluster to study the geometry and the vibrational spectra from Fourier transformed velocity auto-correlation functions. The converged vibrational modes are obtained from a trajectory of about 46 ps (corresponding to 96,000 steps). We find a blue-shift of the Si-Si vibrational modes with transverse acoustic (TA) characters and a redshift of the other vibrational modes with increasing temperature. We also see a broadening of all the vibrational modes with increasing temperature. The former can be linked to the negative (blue-shift) and positive (red-shift) Grüneisen parameters along with the extended bond lengths. The latter is attributed to the low symmetry (proximity to the surface) enhancing anharmonic effects at high temperatures. We further present a computational approach that enables the extraction of inter-vibron coupling strength, taking all the anharmonic effects into account. We find, for our cluster, couplings in the range of 0.15 to 2.5 THz which correspond to "Rabi" oscillations of the vibron occupation number in the range of 0.4 to 6 ps.

AIMD simulations
We construct a sila-adamantane (Si 10 H 16 ) cluster that has the T d point group symmetry [38,39], as shown in Figure 1. In this cage-shaped cluster, four silicon atoms are bonded to three other silicon atoms and terminated by one hydrogen atom, while the remaining six silicon atoms are bonded to two silicon atoms and saturated with two hydrogen atoms. The simulations are performed using density functional theory (DFT) within the local density approximation (LDA) and Trouiller-Martin norm-conserving pseudopotentials with an energy cutoff of 20 Ry [40]. The supercell is simple cubic with an extent of 15Å in each direction. The initial geometry of the cluster is optimized until the minimum forces acting on the Si and H atoms are less than 3×10 −6 a.u. under constrained symmetry.
The AIMD simulations are initially performed at a constant temperature (NVTensemble) using the Nosé-Hoover chain thermostat [41] with time step of 20 a.u. (about 0.48 fs) in order to equilibrating the system. Following the 2 ps equilibration in the NVTensemble, we start a constant energy (NVE-ensemble) simulation and the trajectories are recorded at each time step. All the AIMD simulations are performed with the CPMD code [40]. In order to improve the statistics, we chop the NVE-ensemble simulation into   Schematic diagram of the AIMD simulation process. We start by an NVT simulation of 5.75 ps. From this simulation, we extract 15 initial configurations for the NVE runs. These are taken after an initial 2 ps of simulation time, with a time interval ∆t of 0.25 ps. We perform 15 NVE simulations ("slices") lasting 3 ps (12 ps) for the calculation of vibrational DOS (vibrational cooling). For each slice we perform a simple moving average, as sketched on the lower right. several slices. Each slice starts from a different NVT-ensemble equilibration time and ends with the same number of NVE-ensemble simulation time steps.

Vibrational DOS
In order to obtain the temperature-dependent vibrational DOS from the AIMD simulation, we calculate the velocity auto-correlation function C(t) [42], where, denotes the time averaged value calculated along the entire trajectory, v i (t k 0j ) is the velocity vector of atom i in slice k at time origin point j. The number of atoms in the cluster, the number of time origin points, and the number of slices are labeled as N, n t , and n s , respectively. The time origin points labeled as t 01 , t 02 , t 0j and t 0nt are given in Figure 2. In the present work, the number of slices n s is taken to be 15 with 6400 time steps and 3200 time origin points in each slice. The NVE-ensemble simulations are performed for a total of 96000 time steps corresponding to about 46 ps simulation time. One slice corresponds to a simulation time of approximately 3 ps.
The vibrational DOS are calculated using the Fourier transform [42] V DOS(ω) and for demonstration of the convergence with the number of slices we use the change of the vibrational DOS showing the deviation of the vibrational DOS between a simulation with k and a simulation with k − 1 time slices. Figure 3 (a) and (b) shows ∆V DOS k (ω) for the low (a) and high (b) frequency range at a temperature of 800 K. In this figure, the maximum deviation is 30%. We observe that the vibrational modes start to converge after 35 ps (12 slices).
Since the limited statistic we obtained from AIMD at low temperature even for a small cluster, we use the DFPT results based on the harmonic approximation of lattice dynamics as a low temperature limit. In this case, the DFPT represents a very good classical (neglecting zero-point motion) approximation. The vibrational frequencies ω and the vibrational eigenvectors U i are obtained from the eigenvalue equation: [43,44] where, M i(j) is the mass of atom i (j), V is the potential, R i(j) denotes the atomic positions. The dynamical matrix elements ∂R i ∂R j are calculated using linear response as implemented in the CPMD code [40]. In this approach, the second derivative of the potential is computed from the linear response of the system to an infinitesimal displacement based on perturbation theory [26].

Energy relaxation
The extraction of the potential energy of a certain vibrational mode, E ν p (t), was first proposed by Ladd et al. [35] and subsequently modified by McGaughey et al. [36] to include the total energy of each mode instead of the potential energy only [36,37]. In the quasi-harmonic approximation, where the changes in bond length due to thermal expansion are included but further anharmonic effects are excluded [20], the total energy of each vibrational mode E ν (t) -proportional to the occupation number of the vibrational mode-is expressed in terms of the time-dependent normal coordinates Q ν (t) where The three-component vectors U ν i in Eq. (6) represent the three components belonging to atom i of the vibrational eigenvectors obtained from Eq. (4). R 0 i is the equilibrium position of atom i obtained from the trajectory of the AIMD simulation as, Based on the quasi-harmonic approximation, the vibrational vectors U ν i used in Eq. (6) are calculated using DFPT with the atomic positions obtained from Eq. (7), i.e. from an AIMD simulation at a certain temperature.
The attenuation of the vibrational amplitude reflects the mode relaxation processes and can be described by the auto-correlation function of the energy fluctuation, written as where, denotes the time averaged value calculated along the entire trajectory and δE ν (t) = E ν (t) − E ν is the deviation from the average vibrational energy E ν .

Temperature dependent vibrational properties
We now describe the temperature dependence of the vibrational properties. In Figure 4, we plot the converged vibrational DOS calculated using Eq. (1) and (2) at (a) 1600 K, (b) 1400 K, (c) 1200 K, (d) 1000 K, and (e) 800 K along with (f) the zero temperature linear-response results from Eq. (4) (all ignore the zero point vibrations). In this work, we use the linearresponse results as a low temperature limit because of the limited statistic we obtained from AIMD at low temperature even for the small cluster. The vibrational eigenmodes obtained from the linear-response calculations are analyzed in terms of their displacement pattern by a projection onto bulk phonon modes (see Ref. [45]) and are divided into acoustic-like Si-Si modes and optical-like Si-Si modes. The vibrations between the silicon and hydrogen atoms are divided into bending Si-H and rotation H-Si-H modes, shear H-Si-H modes and stretching Si-H modes according to the displacement of the atoms. These vibrational modes are listed in Table 1 and labeled as M 1 -M 5 in Figure 4 (f). In the AIMD simulations, the character of the vibrational eigenmodes cannot be obtained directly [46].
From Figure 4, we observe that the vibrational DOS of all the modes, especially for the Si-H modes, show a broadening with increasing temperature and that the vibrational density of the shear H-Si-H modes and the stretching Si-H modes decreases with increasing temperature.   There is also a red shift with increasing temperature of all the vibrations except for acousticlike Si-Si modes, which correspond to the TA phonon modes of bulk Si.
To explain the reason for this temperature dependence, we plot in Figure 5 the average bond lengths of (a) Si-Si bonds and (b) Si-H bonds as a function of the distance of the respective bond center to the cluster center at T = 800, 1200, and 1600 K. In the cageshaped Si 10 H 16 cluster, the bond lengths of each Si-Si bond (Si-H bond) are uniform after geometry optimization at zero temperature. We plot these optimized bond lengths as dashed lines in Figure 5. We see from Figure 5  temperatures. The increase of bond length and its large distribution along with the positive Grüneisen parameters (describing the change in phonon frequencies with volume) for the optical branches and longitudinal acoustic (LA) branch result in the softening (red shifting) and broadening of the optical-like and the LA-like Si-Si modes with increasing temperature. In contrast to the relatively small extension of Si-Si bonds (≈ 3%), the Si-H bonds show a large extension and distribution with increasing temperature. A bond length extension of as much as 9% at T = 1600 K is seen in Figure 5 (b). We attribute the large extension of the Si-H bonds to the geometry of the Si 10 H 16 cluster. In contrast to the silicon atoms, which are localized at the center of the tetrahedral structure (see Figure 1), the hydrogen atoms at the cluster surface are free to move outwards. This introduces a potential asymmetry towards the vacuum side and increases the anharmonicity. A relatively large extension of Si-H bonds is therefore obtained at high temperature, which explains the red shift and the broadening of the Si-H vibrational modes, considering the positive Grüneisen parameters. Especially, the vibrations of the high frequency H-Si-H shear and Si-H stretching modes mainly consist of the displacements of hydrogen atoms. This leads to a significant reduction and broadening of the vibrational peaks corresponding to the H-Si-H shear modes and Si-H stretching modes in Figure 4 (a)-(e). Finally, the blue shift of the TA-like Si-Si modes observed in Figure 4 is the result of the bond length extension along with the negative Grüneisen parameters of the TA modes.

Vibron-Vibron Coupling
The vibrational energy of a certain mode is proportional to the mode occupation number E ν (t) ∝ n ν (t) (see Eq. 5), and thus the time evolution of E ν (t) carries the information about the coupling between different vibrational modes. In the following we will refer to E ν (t) as to the "occupation auto-correlation function", to avoid confusions with the many variants of the energy auto-correlation function. The total vibrational enegy, ν E ν , is conserved in our NVE simulation and energy is allowed to flow between vibrational modes. To extract a meaningful statistical average out of E ν (t) we use the correlation function C ν E (t) (Eq. 8), which decays to zero for large times t, when the occupation number is unrelated to the initial situation (at time t 0 ). We have calculated the occupation auto-correlation functions for all the vibrational modes of the Si 10 H 16 cluster at T = 50, 100, 125, and 150 K. In figure 6 (a)-(f) we plot the vibrational occupation autocorrelation function of six selected vibrational modes of the Si 10 H 16 cluster at T = 100 K, which will be discussed in more detail.
It should be pointed out that the AIMD simulations in this work are performed only at relatively low temperatures in order not to violate the quasi-harmonic approximation which is used to estimate the vibrational eigenmode energy [35,36]. Indeed, the vibrational vectors U ν i are calculated using DFPT with the average finite-temperature atomic positions obtained from Eq. (7). Imaginary vibrational frequencies are introduced when the atoms deviate far from their equilibrium positions at high temperatures.
First, we have to prominently assert that a lifetime, as can be extracted from a damped Rabi oscillation or from transitions treated at the level of Fermi's golden rule does not exist in our cluster. The NVE simulation does not allow for energy dissipation, that would possibly lead to exponential decay of high energy phonon modes, from which a lifetime could be extracted. On the other hand, the discrete and energetically sparse nature of the vibrons does not allow for a treatment following Fermi's golden rule, where transitions into a perfectly dissipative continuum are assumed. Our occupation numbers describe the time evolution of a many vibronic level system without dissipative term. The change of occupation of the individual levels reflects the interlevel coupling. The strongest coupling will dominate the short time evolution.
In figure 6 (a)-(f) we observe a rather different behavior of the different vibrational modes. In all cases, however, there is an ultrafast oscillation superimposed on a slower variation. To analyze these results in a quantitative way, we replot in figure 7 (a) the occupation autocorrelation function for the vibrational mode with frequency 387.9 cm −1 (Fig. 6 (c)) along with its Fourier transform (panel b). The Fourier analysis shows that the high frequency oscillation (Period II in 7 (a)), in this case at 11.86 THz, corresponds to the vibron frequency. The origin of this oscillation resides in errors in the projection (Eq. (6)). Indeed, we project our AIMD anharmonic results onto the harmonic vibrations of the NC. This represents an approximation with an error that is largest when the atomic displacements have the largest amplitude. Accordingly, the error increases and decreases periodically with a frequency equal to the vibrational mode's frequency. The slower variation (Period I in 7 (a)) is very clear from Fig. 7 (b) and has a frequency of 1.89 THz. We identify this feature as the frequency with which the vibron state undergoes periodic oscillations with strongly coupled vibronic states. It describes how the specific vibron mode can transfer energy to other modes, i.e., decay into lower energy vibrations or how two, or more, low energy vibrations can create a high energy vibron. Both processes are present in our AIMD calculations. The fact that we are dealing with a many-level (we have 78 vibrational eigenmodes) system explain that we do not see clear Rabi oscillations between two distinct levels but a fast oscillation (1.89 THz) modulated in time by the interaction with more weakly coupled levels. We cannot extract the weak coupling components from our results (that are all hidden in the low frequency Fourier transform Figure 7 (b), and would require much longer simulation times to be properly captured), but can clearly obtain the strongest coupling.
In Figure 7 (c) we summarize our Fourier analysis of all the vibron modes and plot the oscillation period as a function of the vibron frequency. We find strongly coupled (short oscillation periods) modes at all available frequencies, starting from the bulk acoustic-like modes to the passivant hydrogen vibrations. The coupling strength varies between 0.15 and 2.5 THz.
We now study the temperature effects on the coupling strength of a certain vibrational mode. In Figure 8 (a)-(d), we plot the occupation autocorrelation functions of a certain acoustic-like mode with varying temperatures and in Fig. 8(e) the extracted coupling strengths. We find that the vibrational frequencies show a red shift with increasing temperature which can be related to the thermal expansion and behaves as expected, and that the coupling strength increases slightly, from 2.06 to 2.16 THz. This increase can be attributed to higher-order anharmonic effects [47,48]. With the increase in temperature, the effects of higher order anharmonicity, which includes four-vibron, five-vibron, and even higher order interactions, become stronger [32].

Summary
We performed first-principles AIMD simulations, within the BO approximation, to study the temperature dependent vibrational properties of a Si 10 H 16 cluster. We first calculate the vibrational DOS using the Fourier transformed velocity autocorrelation functions and compare the results to DFPT, obtaining good agreement. We quantify the softening and broadening, with increasing temperature, of the Si-Si LA-like modes, the Si-Si optical modes and the Si-H modes (especially for the H-Si-H shear modes and Si-H stretching modes) and analyze the results in terms of bond length variations. Subsequently, we suggest a method to calculate the vibron-vibron coupling strength, including all the anharmonic effects, using a Fourier transformation of the occupation autocorrelation function. For our Si cluster, we find a coupling strength between 0.15 and 2.5 THz, which corresponds to oscillations in the occupation of vibron states (akin Rabi oscillations) in the range of 0.4 to 6 ps. The obtained coupling parameters could be used in further study of the dynamical processes in nanosystem based on rate equations and potentially including dissipation through the coupling to an external phonon bath [49,50]. We conclude that, although our approach enables the extraction of anharmonic vibron coupling parameters otherwise not available, it still requires trajectories of around 24 ps, which presently limits it's applicability to rather small structures. Further challenges involve the non-adiabatic coupling to electronic states and the coupling to a dissipative environment, which we have ignored.