Electronic quantum fluxes in vibrating symmetric and polar single, double and triple bonds

ABSTRACT We investigate electronic quantum fluxes during large amplitude nuclear vibrations of hydrocarbon (HnCCHn), hydrosilicon (HnSiSiHn) and organosilicon (HnSiCHn) compounds with n = 1, 2, 3. The total electronic fluxes are analysed in terms of contributions from localised molecular orbital densities. Furthermore, the vibrationally induced charge transfer in the polar compounds is investigated in terms of the underlying fluxes.


Introduction
The chemical bond is one of the most important concepts in chemistry. It is used in all fields of chemistry and on all levels, from school education to research, to describe and understand molecules and their reactions. Several approaches exist to identify and classify chemical bonds in a molecule, once a static electronic structure has been obtained. Some rely on localised orbitals, like the natural bond analysis [1]. Others analyse the electron density distribution, for example, Bader's theory of atoms in molecules [2] and the electron localisation function (ELF) [3]. In view of the progress in experimental techniques that allow to monitor electronic motions in molecules, it is desirable to extend bonding analyses to the time-dependent domain. For example, a time-dependent ELF for excited electronic states has already been developed [4][5][6]. Here, we analyse the electron dynamics during large amplitude nuclear vibrations of different types of chemical bonds in the electronic ground state in terms of electronic fluxes.
The real-time observation, control and modelling of chemical reactions and vibrations at the molecular level has made significant advances. The field of femtosecond chemistry is well established and provides insights on the nuclear dynamics of molecules [7]. In quantum reaction dynamics, chemical reactions are studied via the evolution of nuclear wave packets on potential energy surfaces. When the nuclei move, so do of course the electrons. In the organic chemists language, arrows indicate the shift of electrons and electron pairs in the respective Lewis structures during a chemical reaction. However, in general, the actual electron dynamics are different. Recent advances in the generation of ultrashort pulses allow the study of electron dynamics even on the attosecond time scale [8,9]. For example, the following techniques are used to image temporal changes in the electronic structure of atoms and molecules: electron diffraction [10][11][12], X-ray scattering [10,[13][14][15][16] and high harmonic generation spectroscopy [17][18][19][20][21][22][23]. Such experiments call for theoretical investigations regarding the dynamics of the electronic structure associated with different kinds of chemical bonds during nuclear motion.
In this study, we analyse changes in the electronic density that occur during large amplitude nuclear vibrations of symmetric and polar single, double and triple bonds. To monitor and quantify the nuclear driven time-dependent changes in the electron density distribution, we employ the concept of quantum fluxes through observer planes as introduced in Ref. [24]. This approach involves a quantum treatment of the electron-nuclear motion on the Born-Oppenheimer (BO) [25,26] electronic ground state potential energy surface. Previously, this method was applied to the investigation of coupled electronic and nuclear fluxes in the vibrating molecules H + 2 [24,27,28] and H 2 [29], to pericyclic reactions [30][31][32] and to the CC bond stretching vibrations of H n CCH n molecules, namely ethane, ethene and ethyne [33]. In the latter study, it was found that, contrary to chemical intuition based on the respective Lewis structures, the electronic rearrangement during large amplitude vibrations is highest for ethane followed by ethene and lowest for ethyne. Here, we further analyse these findings and extent the study to homonuclear H n SiSiH n compounds and to heteronuclear polar single, double and triple bonds in H n SiCH n molecules with n = 1, 2, 3 and investigate how the polarity of the SiC bond is reflected in the electronic flux.
Specifically, we partition the total electronic flux into contributions from different localised molecular orbitals [34], thereby allowing to distinguish the respective contributions from individual chemical bonds. Such a partitioning was applied also in the study on bond-to-bond fluxes during the Cope rearrangement of semibullvalene [30,31] and during tunnelling in cyclooctatetraene [32]. Furthermore, these different contributions to the total electronic flux were recently predicted to become observable experimentally via time-resolved X-ray scattering [35].

Theory and methods
The coupled dynamics of electrons and nuclei is obtained as solution of the time-dependent molecular Schrödinger equation Here, mol into an electronic wave function ψ (q; Q) and a timedependent nuclear wave function χ (Q, t ). The electronic wave function ψ (q; Q) is a solution of the timeindependent electronic Schrödinger equation with the electronic HamiltonianĤ el (r el ; R) and depends parametrically on the nuclear position coordinates. In our model, the nuclear motion is restricted to the reduced nuclear coordinate ζ , which is equal to the distance between the centres of mass of rigid XH n and YH n groups (X = C, Si, Y = C, Si, n = 1, 2, 3). The origin of the coordinate system is chosen as the centre of mass (COM) and the XY-bonds are oriented along the z-axis, c.f. Figure 1.
The corresponding time-dependent nuclear Schrödinger where the SiC bond is elongated and at t = t tp , where the bond is shortened with respect to the equilibrium bond length during the investigated model vibration of the SiC bond. Also shown are three sets of electronic observer planes: the black plane, located at the midpoint of the bond, monitors vibrationally induced charge transfer and separates the domains associated with silicon and carbon groups, S SiH 2 and S CH 2 , respectively. The green and red observer planes, located at symmetric positions above and below the midpoint, monitor the electronic rearrangement close to the carbon and silicon groups, respectively. The direction of electronic flux associated with positive values is indicated by the arrows attached to the observer planes and the centre of mass (COM) is marked by a cross.
where the reduced mass of the systems is μ = M XH n · M YH n /(M XH n + M YH n ), and M ZH n is the mass of the ZH n group. The nuclear spin degrees of freedom are neglected in our model. For the electronic flux calculation, we follow the approach given in Ref. [24] and first calculate the timedependent electron density ρ el (r, t ) from the timedependent BO nuclear density ρ nu (ζ , t) and the timeindependent BO electron densities ρ el (r; ζ ) according to Next, the electronic flux through observer planes perpendicular to the z-axis ( Figure 1) centred at variable positions z obs is calculated according to whereρ el (z, t ) = ρ el (r, t )dxdy denotes the timedependent reduced electron density along the z-axis and N el (z obs (t ), t ) is the mean number of electrons in the given volume. The time-integration of the electronic flux gives the electronic yield (7) which quantifies the number of electrons that flow through the given electronic observer plane within the time interval [0,t].

Computational details
We model the organosilicon and hydrosilicon compounds as structurally equivalent to the respective hydrocarbons, thus constraining HXYH to linear and H 2 XYH 2 to planar structures. For studies on the structure of hydrosilicon and organosilicon compounds and multiple bonding between group 14 elements, see, for example, Ref. [36][37][38][39][40][41][42][43][44][45] and references therein. The structures and potential energy curves (PEC) are obtained using the quantum chemistry package MOLPRO [46]. Structure optimisations are carried out at the CCSD(T)/cc-pVTZ [47] level of theory. The one-dimensional PEC for the electronic ground states were calculated at the MRCI/cc-pVTZ [48][49][50] level of theory. The underlying active space in the CASSCF calculations [51,52] is chosen to allow a description of the homolytic dissociation of the respective bonds. Hence, for molecules with single bonds, the [2,2] active space including the σ XY -bonding and the corresponding anti-bonding orbital is used. Active spaces [4,4] and [6,6] containing σ XY and π XY bonding and anti-bonding orbitals are used for molecules with double and triple bonds, respectively. The time-independent electron densities are calculated on grids along ζ by means of DFT-B3LYP/cc-pVTZ [53][54][55] calculations, yielding excellent agreement with the corresponding densities from Ref. [33] obtained from MRCI/cc-pVTZ calculations. Pipek-Mezey [34] localisation is used to localise the Kohn-Sham molecular orbitals. The numerical calculation of the densities from the electronic wave functions is carried out with the software Orbkit [56]. For each molecule, 64 electron densities are sampled at constant step sizes within specific ranges of ζ given in the Supplemental data. The electron densities are calculated with a spatial resolution of 1.85 · 10 −5Å3 in boxes of dimension 601(z) × 301(y) × 301(x) centred at the COM.
The total electronic flux is partitioned into contributions from different localised orbital densities. Hence, for molecules with single bonds, the total electron density is partitioned as ρ el = ρ core + ρ σ XY + ρ σ X(Y)H , and for double and triple bonds as ρ el = ρ core + ρ σ XY + ρ π XY + ρ σ X(Y)H . The core densities are comprised of 2 electrons per C and 10 electrons per Si nucleus. The valence density, ρ val = ρ el − ρ core , holds 10 (n = 1), 12 (n = 2) and 14 (n = 3) electrons. The X(Y)H n density represents 2n electrons. The density we associate with a XY bond is separated in a σ -type part with 2 electrons and a π-type part with 2 (n = 2), and 4 (n = 1) electrons. These partitioned valence densities are depicted in Figure 2. Visualisations of densities were created using the software ZIBAmira [57].
The nuclear Schrödinger equation is solved numerically with the split operator and the Fourier grid hamiltonian methods, as implemented in the program WavePacket [58]. In each case, temporal steps t = 0.1 fs and grids of 512 equally spaced points within specific ranges along ζ , given in the Supplemental data, are used.
Throughout, we compare the electronic yield for bonds vibrating with equal initial amplitudes of A = 0.5Å. Specifically, the initial nuclear wave packet is the vibrational ground state of the respective PEC, which is displaced in the direction of increasing bond distance ζ , such that all molecules exhibit the same initial amplitude of vibration A = ζ (t = 0) − ζ (t = t p ). Here, ζ (t) denotes the centre of the nuclear wave packet at time t, and t p is the time when the respective nuclear wave packet reaches the inner turning point the first time. This corresponds to exciting the different molecules by a different amount of energy, which was previously shown to be an adequate strategy for the comparison of different chemical bonds in terms of electronic fluxes [33]. The specific parameters for these equal amplitude nuclear vibrations, as well as vibrational periods are given in the Supplemental data. Also given in the Supplemental data are depictions of the PEC and superimposed nuclear wave packets. Moreover, the reader is referred to, respectively, Refs [59,60] and [61,62] for theoretical and experimental work on the initialisation of small and large amplitude molecular vibrations.

Results and discussion
In order to compare and analyse the electronic rearrangement of the symmetric and asymmetric compounds in the bonding region, we determine the electronic flux and yield through two observer planes at symmetric positions above and below the midpoint of the respective bonds. These observers are depicted as green and red lines in Figure 1. The displacement from the midpoint, z obs (Figure 1), is chosen such that the observer planes coincide with the mean positions of the respective C-or Sinuclei at their respective first turning points. Such molecular specific observer planes were shown in Ref. [33] to be well suited to compare the electronic rearrangement of single, double and triple bonds in the bonding region during large amplitude vibrations in the symmetric carbon compounds. Moreover, it was shown in Ref. [33] that the drawn conclusions are stable over wide ranges of the displacement z obs .
The time evolution of the electronic yield for the nine molecules under investigation with respect to molecular specific observer planes near the C-and Si-atoms are shown in Figure 3 as solid and dashed lines, respectively. Here, |Y el | corresponds to electron accumulation in the bonding region during the large amplitude vibration, c.f. Equation (7). The period of the timedependent electronic yield is equal to the vibrational period of the respective system, i.e. the electrons follow the nuclei. An increase in the electronic yield happens as the SiC bond shortens. The maximal electronic yield is reached as the bonds are most compressed. During bond stretching, the net electronic flux changes direction, such that a decrease in the electronic yield occurs until the outer turning point is reached. Those time-dependent Figure . Time evolution of electronic yields for symmetric molecules H n CCH n (left) and H n SiSiH n (right) and for asymmetric molecules H n SiCH n (middle) with n =  (black), n =  (blue) and n =  (red). Results for molecular specific observer planes ( Figure  and text) located near the C-atoms (Si-atoms) are shown as solid (dashed) lines. Moreover, the total electronic yields (|Y el |) are decomposed into contributions from core densities (core |Y el |), valence densities (val |Y el |), densities associated with σ -bonds (σ |Y el |), with π-bonds (π |Y el |) and from CH-and SiH-bonds (σ CH Y el and σ SiH |Y el |, respectively). In all cases, |Y el | corresponds to electron accumulation in the bonding region during the large amplitude vibration, c.f. Equation ().  (Figure ), partial charges δ + (t) in sector S SiHn and the maximal electronic yields Y el, max , as monitored by the electronic observer plane located in the midpoint of the bond ( Figure ). All values are given at t =  fs and t = t tp and are calculated with respect to the valence, σ SiC and π SiC densities of the H n SiCH n molecules.  [24,[27][28][29]33]. The total (Y el ) and valence (val Y el ) electronic yields of the hydrocarbons (left panel) show the trend, that more electrons participate in the electronic rearrangement in vibrating ethane compared to ethene and ethyne [33]. The same is true for the hydrosilicons (right panel) in that disilane has a higher valence electronic yield than disilene and disilyne, with the latter two having nearly identical maximal valence electronic yields. Insight into these trends may be obtained from the separation of the total yields into contributions from different localised orbital densities [30,31,35]. The bottommost panel in Figure 3 shows that the high contribution of the CH-and the SiH-densities are the reason for this dominating electronic rearrangement in the bonding region for the single bonded systems. Subtracting the corresponding fluxes from the total electronic flux gives indeed the intuitive picture expected from the respective Lewis-structures, i.e. largest electronic rearrangement in the bonding region for the compounds with triple bonds, then for double bonds, and last for single bonds.
Moreover, the valence electronic rearrangement (val Y el ) in the bonding region is higher for the hydrocarbons than for the hydrosilicons. This is due to the more diffuse valence electron density in the hydrosilicons ( Figure 2). Also, due to the higher electronegativity of carbon compared to silicon, this effect is even more pronounced in the polar organosilicon compounds (middle panel of Figure 3).
For the analysis of the electronic rearrangement in the polar H n SiCH n -compounds, we monitor the electronic yields through an electronic observer plane placed in the midpoint of the respective bonds (black line in Figure 1). Initially, the electron densities are polarised in the direction of the more electronegative carbon-atom with the polarity increasing with increasing bond order (Table 1 and Figure 2). The corresponding time-evolution of the electronic yields monitoring the dynamical bond polarity is shown in Figure 4. Here a positive electronic yield corresponds to vibrationally induced charge transfer from the C-side to the Si-side of the molecule upon bond compression. We observe an increase of the charge transfer from the single, via the double to the triple-bonded molecules. While the contribution of the σ SiC densities to the overall flux is comparable in all cases, 0.09 electrons are added to the electronic flux through the bond midpoint due to the π SiC density in the case of silene and 0.15 electrons in the case of silyne. This vibrationally induced charge transfer from the silicon to the carbonatom is supported by the fluxes of the SiH-bond densities, while it is diminished by the oppositely directed fluxes of the CH-bond densities (bottommost panel in Figure 4). The flux associated with the CH-bond densities is stronger than for the SiH-bond densities, which is mostly due to shorter CH-bond distances and less diffuse densities. This leads to an overall decrease of vibrationally induced charge transfer, which is strongest for silane.

Summary
We have presented an analysis of electronic quantum fluxes during large amplitude vibrations of polar and non-polar single, double and triple bonds, including the separation of the total electronic fluxes into contributions from individual chemical bonds. The dominant electron rearrangement in the bonding region of ethane compared to ethene and ethyne [33] is confirmed for the corresponding symmetric hydrosilicon and asymmetric organosilicon compounds. Here, we attribute this counter-intuitive result to the large contribution of the CH (SiH) bond densities in the single-bonded systems, surpassing by far the flux contribution of σ and also of π bond densities. Furthermore, the less diffuse electron Figure . Time evolution of the electronic yield through one observer plane located at the SiC bond center of the H n SiCH n molecules. Positive yields are due to electronic flux from the SiH n to the CH n side, negative yields are due to fluxes in the opposite direction. The first sub-figure shows the total electronic yield due to the changes in the total electron density. The second to fifth sub-figures show the electronic yields due to changes in the core electron, σ SiC bond , π SiC bond and σ SiH and CH electron densities. The last figure shows the results for the σ SiH and σ CH bond densities separately. density of the hydrocarbons compared to the hydrosilicons is reflected in stronger electronic fluxes of valence densities in the hydrocarbons. Finally, we showed that the vibrationally induced charge transfer in the polar organosilicon compounds is nicely reflected in the total electronic fluxes, showing strongest associated flux for the triple-bonded system. In summary, electronic quantum fluxes can be regarded as a valuable tool for timedependent bonding analyses.