Charge transfer in ultrafast isomerization of acetylene ions

First-principle calculations are employed to investigate the ultrafast isomerization of the acetylene cation and dication. We use the time-dependent density functional theory together with the Ehrenfest dynamics to track the coupled electron-nuclear dynamics. For both the acetylene cation and the dication, we observe nonadiabatic behaviors during the isomerization. We ﬁnd that the charge transfer not only alters the electronic structure through nonadiabatic transitions, but also plays a key role in the subsequent hydrogen migration. We show that nonadiabatic transitions affect the structural modiﬁcation of the excited potential energy surface, and also facilitate the ultrafast isomerization through the creation of a channel of increased negative charge that facilitates the proton movement. For the acetylene cation, we ﬁnd a timescale for hydrogen isomerization of 66 ± 4 fs, which is consistent with previous pump-probe experiments and on-the-ﬂy calculations. For the dication, we ﬁnd nonadiabatic transitions occur before the isomerization and identify a similar channel for the proton. Moreover, we ﬁnd the formation of vinylidene-like structures is always accompanied by a characteristic charge separation on the carbon skeleton. These heuristics will be useful in identifying tautomers and motivating the ultrafast charge-transfer detection methods for future experiments.

First-principle calculations are employed to investigate the ultrafast isomerization of the acetylene cation and dication. We use the time-dependent density functional theory together with the Ehrenfest dynamics to track the coupled electron-nuclear dynamics. For both the acetylene cation and the dication, we observe nonadiabatic behaviors during the isomerization. We find that the charge transfer not only alters the electronic structure through nonadiabatic transitions, but also plays a key role in the subsequent hydrogen migration. We show that nonadiabatic transitions affect the structural modification of the excited potential energy surface, and also facilitate the ultrafast isomerization through the creation of a channel of increased negative charge that facilitates the proton movement. For the acetylene cation, we find a timescale for hydrogen isomerization of 66 ± 4 fs, which is consistent with previous pump-probe experiments and on-the-fly calculations. For the dication, we find nonadiabatic transitions occur before the isomerization and identify a similar channel for the proton. Moreover, we find the formation of vinylidene-like structures is always accompanied by a characteristic charge separation on the carbon skeleton. These heuristics will be useful in identifying tautomers and motivating the ultrafast charge-transfer detection methods for future experiments. DOI: 10.1103/PhysRevA.106.033111

I. INTRODUCTION
Isomerization, as an elementary bond rearrangement process, plays a key role in chemistry and biology in which the conformational change can be induced by light absorption or ionizing radiation [1][2][3][4][5]. A simple isomerization occurs in acetylene ions following excitation or ionization, i.e., [ shift reactions widely existing in organic chemistry [6]. In recent years, both experiments and theoretical calculations confirmed an extremely short timescale (usually below 100 fs) of the ionic isomerization [7][8][9]. This ultrafast isomerization of acetylene ions involves the rapid hydrogen migration driven by complex many-body interactions, which can potentially reveal the correlations between electrons and nuclei that are otherwise hidden in adiabatic properties of a system. These correlations, and the coupled electron-nuclear dynamics they evoke, impact fundamental processes, such as electronic transport, catalytic efficiency, or photoelectric energy conversion.
In both the cation and the dication exists a very interesting feature that the isomerization usually goes hand in hand with the ultrafast electronic decay. Jiang et al. reported a cationic isomerization time of 52 ± 12 fs and confirmed the isomerization is related to the nonradiative decay of the A 2 + g state [8]. Madjet et al. substantiated that this cationic isomerization is associated with nonadiabatic transitions from the A 2 + g to X 2 u state [25]. For the dication, Li et al. pointed out that a nonadiabatic pathway involving the decay of the 1π −1 u 3σ −1 g state can possibly account for the sub-100-fs isomerization [26]. Moreover, the ultrafast decay of the A 2 + g state is responsible for the isomerization recurrence behavior between A and V structures [22], and as such to and fro isomerization can be interpreted as the result of the internal conversion from excited states to the ground X 3 − g state [19,20]. Although a large number of experiments or electronic structure calculations [27][28][29] were carried out for the acetylene isomerization, a time-resolved electron-nuclear dynamics from first principles, in particular, involving the ultrafast electronic decay, is still sparse.
In this paper, we focus our attention on the charge transfer in the ultrafast isomerization of both the acetylene cation and the dication by employing the time-dependent density functional theory together with the Ehrenfest dynamics (Ehrenfest-TDDFT) [30,31]. The time-resolved relaxation dynamics of acetylene isomerization is presented by the Ehrenfest-TDDFT method. Different from charge migration purely driven by electronic correlations [32], charge transfer in this work is found to consist of two processes coupled with ionic motion, i.e., nonadiabatic electronic transitions between different potential energy surfaces (PESs) and the intramolecular charge transfer together with the ionic motion. We find that the charge density difference (CDD) is a signature of the nonadiabatic charge transfer, quantitatively featuring the structural modification of the excited PES [30,33]. As an inherited effect from the CDD, we find the hydrogen isomerization as well as the to-and-fro isomerization can be facilitated by the creation of an attractive surplus of negative charge that acts as a channel for the proton movement. In this work, this proton channel is a dynamical consequence of the intramolecular charge transfer. Moreover, we find the isomerization is constantly accompanied by asymmetric charge distributions on the carbon skeleton, which may be used as an indicator to identify tautomers in future experiments. This paper is organized as follows. In Sec. II, we briefly introduce the Ehrenfest-TDDFT method and numerical details. In Sec. III, we discuss the isomerization mechanism, especially the intramolecular charge transfer. Finally, the conclusion is drawn in Sec. IV.

A. Ehrenfest-TDDFT method
The electron-nuclear dynamics is modeled by the Ehrenfest-TDDFT: the electrons are described by the Kohn-Sham (KS) TDDFT, nonadiabatically coupled to the nuclei that are described as classical point charges. Specifi-cally, we solve the following set of coupled equations (atomic units): where ψ j (rσ, t ) are the KS orbitals, r is the electronic coordinates, and σ is the electronic spin. R I (t ) are the nuclear coordinates and N ion is the total number of nuclei. M I and Z I are, respectively, the mass and charge of the Ith nuclear. The KS effective potential consists of the Hartree potential V H , the exchange-correlation (xc) potential V xc , and the electronnucleus potential V en . The total electronic density ρ(r, t ) is computed by summing up all the spin-up and -down densities To see the impact of charge transfer on the ultrafast isomerization, we introduce the charge density difference ρ q+ CDD for an instantaneous ionic configuration, defined as the time-dependent electronic density subtracted from the instantaneous DFT electronic density , t] is computed using the standard DFT method, ensuring the holes describe the same as that in TDDFT.

B. Numerical details
Our simulations begin with the preparation of initial states. The initial ionic configurations of [HCCH] q+ are generated using a similar technique as that in Ref. [25]: the two hydrogens are sampled from a Wigner distribution of positions and momenta on the ground-state PES of the neutral acetylene, while the two carbons are initially fixed at the equilibrium positions. Second, for the specified ionic positions, the electronic states of the ionized acetylene are initialized using DFT: we use the standard energy-functional minimization with a hole localized in 3σ g for the [HCCH] + and two holes localized in 1π u 3σ g (with one electron in each of π and σ orbitals) for the [HCCH] 2+ , which leads to the A 2 + g and 1 u states [26,34], respectively. This is an improvement over the sudden removal of electrons from Kohn-Sham orbitals (KSOs), which could cause an incorrect charge dynamics [35,36]. Throughout this paper, the state label is referring to the electronic state in the linear acetylene structure regardless of the molecular geometry.
The coupled electron-nuclear dynamics is tracked using the real-time and real-space Ehrenfest-TDDFT in a development version of OCTOPUS [37], in which the electronic dynamics is described quantum-mechanically [30], nonadiabatically coupled to the ions through the Ehrenfest molecular dynamics [31]. For the cation, 563 trajectories are used to extract the ensemble-averaged results and the isomerization is identified in 225 trajectories. For the dication, we find the isomerization is a rare channel. All the trajectories are propagated with a time step of 8.45×10 −4 fs, and a local density approximation (LDA) [38] is employed for the exchange-correlation functional. The vinylidene-like (V-like) structure is defined as that in which the two hydrogens belong to the same side with respect to the imaginary carbon-carbon (C-C) midplane.

III. RESULTS AND DISCUSSIONS
For the acetylene cation, we first show in Fig. 1(a) the ensemble-averaged vinylidene fraction. It is seen that the vinylidene emerges at around 38 fs, followed by a rapid increase until 77 fs. The half of the maximum fraction is reached at 66 ± 4 fs in our simulations, which is consistent with the measured timescale of 52 ± 12 fs found in Ref. [8]. These results are also consistent with previous surface-hopping calculations (∼60 fs) [25].
For one of the representative isomerization trajectories, the electron-nuclear dynamics is presented in Fig. 1(b). Within the first 30 fs, the CDD is negligible, which means that the time-dependent propagation during this period almost occurs on the A 2 + g PES. After 32.1 fs, the CDD is clearly seen and it later coherently evolves with the isomerization, i.e., from that moment on the time-dependent potential is remarkably different from the A 2 + g PES. A movie of the CDD for this cation trajectory is available in Gitlab Digital Commons [39]. Notably, the conformation transforms from the A-to the V-like structure at 106.5 fs and returns again to the A-like structure at 124.4 fs, which may be facilitated by the appearance of a region of negative CDD that acts as a channel of attraction for the proton. As the isomerized hydrogen reaches the C-C imaginary midplane to form the V-like structure, the C-H bond region is consistently dominated by the negative CDD (golden isosurfaces). One can analytically derive from the Ehrenfest-TDDFT that the negative CDD generates an attractive potential in addition to the excited A 2 + g potential (see more details in Appendix A), and such a negative CDD region opens a barrierless route for the isomerized hydrogen. It should be mentioned that all isomerization trajectories show clearly that the nonadiabatic transition occurs before the A-to-V conformational change, which is consistent with on-the-fly calculations. For example, based on the complete active space self-consistent field (CASSCF) method with surface-hopping, Madjet et al. concluded that the electronic decay happens first and then isomerization follows [40].
To further investigate nonadiabatic effects during the isomerization dynamics, we analyze the time evolution of its energies computed by the Ehrenfest-TDDFT. For comparison, the Car-Parrinello (CP) dynamics is also run with identical initial conditions, in which the ions propagate rigorously on the adiabatic A 2 + g PES. One can see in Fig. 1(b) how Ehrenfest energies are hardly different from those obtained with the CP dynamics before 30 fs, confirming again that the Ehrenfest-TDDFT dynamics evolves on the A 2 + g PES for a certain period. Afterwards, in order for an isomerization pathway to open, nonadiabatic effects are essential since a sufficient amount of kinetic energy is required to overcome the isomerization barrier, whereas a kinetic energy is missing in the adiabatic trajectories. It is emphasized that the total energy of the cationic trajectory is conserved because the entire system is isolated. In Fig. 1(b), the total energy can be divided into three parts: the potential energy (E pot ), the ionic kinetic energy of the isomerized H b hydrogen (E iso kin ), and the total ionic kinetic energy of the three other atoms (i.e., H a , C a , C b ). For the Ehrenfest-TDDFT trajectory in Fig. 1(b), the initial the ions is 0.27 eV, which is less than the A-to-V (more than 1 eV [40]) and the V-to-A barrier (∼0.5 eV [29]). In the final state at t = 150 fs, the total ionic kinetic energy is 2.10 eV. In Ehrenfest-TDDFT simulations, the E pot drops quickly after 30 fs, and is partially converted to E iso kin exceeding 1 eV. In contrast, the E pot in the CP dynamics always vibrates around the equilibrium and the E iso kin obtains less than 0.5 eV. Therefore, the CDD modifies the potential felt by the isomerized hydrogen along the isomerization route, allowing more energy to be redistributed to E iso kin . For the cation trajectory, there is a barrier region along the isomerization path (see more details in Appendix B). When the isomerized hydrogen (H b ) approaches the barrier region, more kinetic energy transfers from the E iso kin to the E pot . It was shown that the fast intramolecular energy redistribution is very important to the isomerization. For example, Madjet et al. also demonstrated a vibrational energy redistribution from isomerization coordinates to slower modes [25]. Based on the multiconfiguration time-dependent Hartree (MCTDH) method, Richter proposed that a fast energy redistribution exists in cis-trans delocalized modes [41].
In Fig. 2(a), we show the potential energy curves (PECs) and projection probabilities of the instantaneous Hamiltonian along the isomerization path to evaluate the nonadiabaticity and the rearrangement of the electronic structure. Immediately, we see a crossing between the excited A 2 + g and the ground X 2 u (with a hole localized in 1π u ) curves around 30 fs, indicating the possibility of state mixtures. This mixing results in the nonadiabatic electronic transition from the 1π u of A 2 + g to the 3σ g of X 2 u , as is quantified by projection probabilities of the transition that ∼81% single-electron population transfer after the mixing. One can find that this nonadiabatic transition causes the pronounced CDD during the isomerization, which can be understood straightforwardly: The nonadiabatic charge transfer makes the potential depart from the A 2 + g PES, resulting in the hole localized on the 3σ g of A 2 + g being filled (see Appendix C for more details). Thus, at the original 1π u location, a positive CDD is identified, whereas a negative CDD is obtained at the original 3σ g location. Since the 3σ g orbital usually distributes along the internuclear axis, the negative CDD shows an intensive distribution on the C-H bond region.
To assess the charge transfer inside the acetylene cation, we compute the number of electrons localized on each atom, i.e., the partial charge [42]. As shown in Fig. 2(b), it fluctuates around the equilibrium charges. However, a distinct feature, depending on the ionic conformation, is immediately noticeable when the cation conforms the V-like structure (shadow region), i.e., an asymmetric charge distribution on the carbons. The maximum of charge separation is as large as 0.45 electrons at 114.6 fs, showing how the hole density [ρ H (r, R, t ) = ρ 0 (r, R, t ) − ρ + (r, R, t ), i.e., the difference between the instantaneous neutral ground-state electronic density and the time-dependent electronic density of the cation] localizes on the C b site [see the inset of Fig. 2(b)]. We explicitly examined this feature in other isomerization trajectories and found that the V-like acetylene structure persistently displays a pronounced charge separation on the carbon skeleton. For nonisomerization trajectories or A-like structures, the partial charges of carbons are almost equal.
Furthermore, we quantify the carbon-carbon charge separation by defining a dimensionless asymmetry parameter as η(τ ) = (P C a − P C b )/(P C a + P C b ) (P C i denotes the partial charge of the carbon i) as a function of the isomerization time delay τ . Figure 3 depicts the asymmetry parameter distributions of the acetylene cation. In this plot, τ = 0 means the time instant of the first A-to-V conformational change for each isomerization trajectory. For τ < 0 (i.e., before isomerization), the cation keeps A-like structure with negligible η values, as expected. After encountering the first A-to-V structural transformation (τ > 0), η splits into two remarkable bridge-like patterns when the V-like conformation exists, and they converge again at a late τ (∼38.7 fs) due to the conformational recurrence: some trajectories return to the A-like structure, in reasonable agreement with the surface-hopping calculations [40]. For the V-like structures, the maximum absolute value |η| reaches 0.05, which may be used in optical experiments for identifying the tautomers. Since the carbons' charge separation is related to the hole distribution, we suggest that future experiments with more advanced orientation control technology may be possible to link such a charge separation with experimental probes, such as positioning holes using the high-order harmonic generation (HHG) spectroscopy [43] and some other time-resolved ultrafast techniques for characterizing the hole migration dynamics [44] or the transition hole dynamics [45].
For the dication [HCCH] 2+ , the nonadiabatic isomerization pathway was theoretically predicted [26], but it is not clear whether or not it follows a similar dynamics as the cation. For this reason, we simulate the dynamics with two types of initial vacancies that are suggested by experiments [23,46] and calculations [26]. For 1π −2 u vacancies, it was found that all simulated trajectories remain in A-like structures, reflecting that the hydrogen isomerization initialized by 1π −2 u vacancies is highly infeasible, which is consistent with a high isomerization barrier of ∼2 eV with 1π −2 u vacancies [27,28]. In contrast, for the dynamics started with the 1π −1 u 3σ −1 g vacancies, all the trajectories dissociate via the C-H bond breakup when the C-C bond length is started at the value in the neutral acetylene. However, if the C-C bonding is initially relaxed to the dication equilibrium length, we identify a nonadiabatic isomerization pathway. As shown in Fig. 4, the dication dynamics starts on an adiabatic propagation with the 1π −1 u 3σ −1 g double hole states, and then undergoes a nonadiabatic transition to 1π −2 u , resulting in the hole on the 3σ g being filled. For the typical trajectory in Fig. 4, initial 1π −1 u 3σ −1 g double hole states correspond to the total energy of −302.11 eV, and at t = 0, an ionic kinetic energy of 0.20 eV is set to initiate the isomerization. In the final state at t = 234.4 fs, the total ionic kinetic energy is 0.94 eV. Note that the nonadiabatic transition happens at around 120 fs, before the A-to-V change at 153.1 fs. Such a nonadiabatic transition produces a pronounced CDD, which facilitates the hydrogen isomerization. In addition, when the dication transforms to vinylidene during the 153.1-189.6 fs time interval, the carbon-carbon charge separation is about 0.4 electrons.
According to the above analysis of nonadiabatic pathways in the cation and dication, our results show that nonadiabatic transitions occur before the A-to-V conformational change. Therefore, under current conditions, the lower time limit of the A-to-V isomerization can be predicted. For the cation, the ensemble-averaged lower onset for nonadiabatic transitions is about 58.5 fs, close to the lower limit of the isomerization time 66 ± 4 fs. For the dication, nonadiabatic transitions happen between 26.3 and 133.2 fs, suggesting that the dication isomerization may be not strictly forbidden within the first 100 fs.

IV. CONCLUSION
In conclusion, we studied the ultrafast isomerization dynamics of both the acetylene cation and the dication within the Ehrenfest-TDDFT method. We presented a time-resolved picture for the ultrafast charge transfer dynamics, featuring the nonadiabatic electronic transition from 1π u to 3σ g along the isomerization route. We find the charge transfer is an important element to understand the isomerization because it not only alters the electronic structure instantaneously by nonadiabatic transitions, but also plays a key role in the subsequent hydrogen migration by modifying certain parts of the excited PES. For both the cation and the dication, we find the CDD facilitates the hydrogen isomerization, as it creates a channel of surplus negative charge, and the carbon-carbon charge asymmetry might be useful for identifying isomers in future experiments. Since the Ehrenfest-TDDFT is a commonly used scheme for describing the coupling between classical and quantum subsystems, our results, which provide insights into how charge transfer impacts the ultrafast isomerization of acetylene ions, have implications for other polyatomic ions.

033111-5 APPENDIX A: CREATION OF A CHANNEL FOR THE PROTON MOVEMENT BY THE CHARGE DENSITY DIFFERENCE
The creation of a channel of negative charge that attracts the proton movement and facilitates the isomerization, thanks to the nonadiabaticity, can be understood within the Ehrenfest-TDDFT method. In this work, the equations for the ionic movement that we used are where r and R i denote the electronic and the ith ionic positions, respectively. ρ TD is the time-dependent electronic density. Z i is the charge of the ith ion. On the right-hand side (RHS) of Eq. (A1), the first potential term is the electron-ion potential (V en ), and the second potential term is the ion-ion potential (V nn ). Then, let us consider the following electronion potential difference (the V nn term is trivial): where the term ρ DFT (r, t ) − ρ TD (r, t ) is the charge density difference (ρ CDD ). Hence, Eq. (A2) can be reformulated as From Eq. (A3), it is clear that the negative parts of the ρ CDD play a role as negative charges, generating an attractive potential in addition to the instantaneous DFT potential. On the contrary, positive parts of the ρ CDD generate a repulsive potential.

APPENDIX B: ENERGY VARIATION ALONG REACTION COORDINATES
We extract reaction coordinates from the typical cationic trajectory to clarify the variation of E pot and E kin using the Ehrenfest-TDDFT. Analogous to the reaction coordinates in Ref. [25], Figs. 5(a) and 5(b) show the energy variation along the intersection angle (θ ) and C b -H b bond length (B C b H b ). The cationic trajectory starts with the point S 0 , and then reaches the nonadiabatic transition point at ∼30 fs. After the transition point, it is notable that a barrier region locates at intermediate values of the angle θ (≈100 • ) before the A-to-V conformational change. Accordingly, as shown in Fig. 5(b), when the cationic trajectory reaches the barrier region before the A-to-V change, more kinetic energy transfers from the isomerized hydrogen to the potential energy. Moreover, for comparison, we plot Figs. 5(c) and 5(d) to clarify the energy variation along another important reaction coordinate, i.e., the C a -C b bond length (B C a C b ). The barrier region located at intermediate θ values (≈100 • ) is also identified.
The existence of the isomerization barrier on the X 2 u state before the A-to-V change is also confirmed by previous calculations. For instance, using the highly correlated CCSD(T) and MRCI methods in Ref. [29], Boyé-Péronne et al. found that the planar isomerization route exists at a barrier hight of ∼2 eV from the X 2 u state of acetylene to V-like structures. More importantly, in Ref. [22], Ibrahim et al. pointed out that employing 800-nm pump schemes is not capable of launching hydrogen migration in cationic states, mainly due to the fact that overcoming the potential barrier (∼2 eV) separating the A-like structure in the ground state X 2 u from the V-like structure is quite challenging. In additon, for the dication trajectory in Fig. 4, we also identify a barrier region at intermediate θ values (≈95 • ) after nonadiabatic transitions, which was predicted by static calculations (Refs. [27,28]) using an isomerization barrier of ∼2 eV on 1π −2 u double hole states.

APPENDIX C: PROJECTION PROBABILITY
The nonadiabatic transition is quantified by the timedependent projection probability, defined as the squared modulus of the overlap between the time-dependent Kohn-Sham orbital ( TD 1π u ,↑ ) and the eigenstate of the instantaneous Hamiltonian ( S i,↑ ), where S i,↑ is either the 1π u eigenstate of A 2 + g ( A 2 + g 1π u ,↑ ), or the 3σ g eigenstate of X 2 u ( X 2 u 3σ g ,↑ ). Figure 6 shows the time evolution of projection probability for the typical cationic trajectory. After the crossing at ∼30 fs, the projection probability on the A 2 + g 1π u ,↑ state goes down to ∼0.09, and the projection probability on X 2 u 3σ g ,↑ jumps to ∼0.81, demonstrating that the nonadiabatic transition occurs before the A-to-V conformational change and the population transfer from the excited A 2 + g to the ground X 2 u is mainly attributed to the single-particle component from the to the X 2 u 3σ g ,↑ eigenstate. Moreover, the projection probability P( TD 1π u ,↑ , g ,↑ ) jumps to ∼82% after nonadiabatic transitions, demonstrating that the 3σ −1 g hole state of A 2 + g is filled.