+ Theoretical study on collision dynamics of H + CH4 at low energies Cong-Zhang Gao, Jing Wang, Feng Wang, and Feng-Shou Zhang

Citation: The Journal of Chemical Physics 140, 054308 (2014); doi: 10.1063/1.4863635 View online: http://dx.doi.org/10.1063/1.4863635 View Table of Contents: http://aip.scitation.org/toc/jcp/140/5 Published by the American Institute of Physics

Articles you may be interested in Collision dynamics of proton with formaldehyde: Fragmentation and ionization The Journal of Chemical Physics 140, 124306 (2014); 10.1063/1.4868985 Perspective: Fundamental aspects of time-dependent density functional theory The Journal of Chemical Physics 144, 220901 (2016); 10.1063/1.4953039 Accuracy and computational efficiency of real-time subspace propagation schemes for the time-dependent density functional theory The Journal of Chemical Physics 144, 204125 (2016); 10.1063/1.4952646 Time-dependent density functional theory Ehrenfest dynamics: Collisions between atomic oxygen and graphite clusters The Journal of Chemical Physics 126, 134307 (2007); 10.1063/1.2713391 Time-dependent density functional theory: Past, present, and future The Journal of Chemical Physics 123, 062206 (2005); 10.1063/1.1904586 + Observation of vibrationally resolved charge transfer in H +H2 at ECM=20 eV The Journal of Chemical Physics 87, 2685 (1998); 10.1063/1.453107

THE JOURNAL OF CHEMICAL PHYSICS 140, 054308 (2014)

Theoretical study on collision dynamics of H+ + CH4 at low energies Cong-Zhang Gao,1,2 Jing Wang,1,3 Feng Wang,4 and Feng-Shou Zhang1,3,5,a) 1

The Key Laboratory of Beam Technology and Material Modification of Ministry of Education, College of Nuclear Science and Technology, Beijing Normal University, Beijing 100875, China 2 Laboratoire de Physique Théorique-IRSAMC, Université Paul Sabatier, F-31062 Toulouse Cedex, France and CNRS, UMR5152, F-31062 Toulouse Cedex, France 3 Beijing Radiation Center, Beijing 100875, China 4 Laser Micro/Nano Fabrication Laboratory, School of Mechanical Engineering, Beijing Institute of Technology, Beijing 100081, China 5 Center of Theoretical Nuclear Physics, National Laboratory of Heavy Ion Accelerator of Lanzhou, Lanzhou 730000, China

(Received 31 July 2013; accepted 17 January 2014; published online 4 February 2014) In this work we make an investigation on collision dynamics of H+ + CH4 at 30 eV by using time-dependent density functional theory coupled with molecular dynamics approach. All possible reactions are presented based on 9 incident orientations. The calculated fragment intensity is in nice agreement with experimental results. The mechanism of reaction transition for dissociation and proton exchange processes is explained by the intra-molecule energy transfer. However, the energy loss of the proton is in poor agreement with experimental results. The discrepancy is attributed to the mean-field treatment of potential surface. We also studied the dependence on initial velocity of both proton and methane. In addition, we find that for dynamical evolution a different self-interaction correction (SIC) may lead to different results, but with respect to the position of rainbow angle, average-density SIC seems to have reasonable correction. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4863635] I. INTRODUCTION

Ion-molecule collisions are considered to play an important role in many fields of science, such as chemistry, plasma physics, material science, astrophysics, and radiation therapy. Accordingly, a large number of ion-molecule collision experiments1 have been conducted, especially at ∼keV energies. Most of their concerned interests are the measurement of cross sections. These experiments promoted the development of early theoretical models relevant to ions’ collisions with molecules. In view of the limitation of experimental techniques (time and energy resolution), there have been few experiments on ion-molecule collisions in the energy range from a few eV to tens of eV. When a projectile with ∼eV collision energy hits a target, some phenomena would occur, such as vibrational excitation, dissociation, charge transfer, and nuclear exchange. The complexity and diversity of collision processes in this energy range presented challenges in experimental measurement. In many experimental ion-molecule collision investigations, the proton as an ion-projectile was favored because of its missing electronic structure and capability of capturing electrons. For this reason, a large amount of experiments have been performed on proton collisions with various molecules,2–8 particularly with hydrocarbon molecules (Cm Hn ). Proton-hydrocarbon collisions have become increasingly important in some applications. For instance, in fusion a) Author to whom correspondence should be addressed. Electronic mail:

[email protected]. Tel.: +86-10-62205602. Fax: +86-10-62231765.

0021-9606/2014/140(5)/054308/12/$30.00

device, the low temperature plasma would lead to the intense chemical erosion of carbon materials covered upon the divertor,9 and then Cm Hn molecules were produced. The collisions between proton and Cm Hn would become extensive, and these have been the subject of many studies. Recently, the hydrocarbon molecules have been detected in the atmospheres of the outer planets,10, 11 which is regarded as the primary evidence of the origin of life. The unknown collision mechanisms, such as energetic particles (i.e., proton) irradiating hydrocarbons, are the most urgent problem. In the family of hydrocarbon, methane (CH4 ) is the simplest, and it could be converted to chain hydrocarbons under the photochemical reactions. Therefore, the study on proton-methane collision is the most fundamental. Several research groups have conducted proton collisions with methane at different collision energies. Koopman12 measured the charge exchange cross sections for the process H+ + CH4 over the energy range 0.1–1.5 keV. McNeal13 reported the total ionization cross section and total electron capture cross section by proton beam scattering with CH4 gas in the range 1–25 keV, and the results had been in agreement with other measurements within the suitable experimental errors. Rudd et al.14 obtained the ionization and electroncapture cross sections for 5–4000 keV proton impacting on CH4 , and the data were fitted to semi-empirical equations. Gao et al.15 provided absolute differential cross sections for charge-transfer scattering of proton from CH4 molecule over the laboratory angular range 0.02◦ –1.0◦ at 1.5 keV, and the data exhibited an oscillatory structure which arose from transitions between the initial and final states of system.

140, 054308-1

© 2014 AIP Publishing LLC

054308-2

Gao et al.

Sanders et al.16 employed 60–120 keV protons in collisions with methane molecules to obtain the electron-capture cross sections, and the cross sections were observed to decrease with increasing energy. Among the numerous experiments on proton collisions with methane, studies on collision energies below 100 eV have been rarely performed. To our knowledge, Toennies et al.17–19 conducted the most comprehensive experimental research for H+ + CH4 system in the energy range of 9.8 eV ≤ ELab ≤ 30 eV. They measured angular distributions and time-of-flight spectra for the protons (inelastic collision) and H atoms (charge transfer collision). They found that in the case of inelastic collisions the average energy transfer (E) showed a nearly linear increase with collision energy within the experimental errors. However, the relative energy transfer did not scale with the collision strength, and this deviation was explained by a strong charge transfer which attended the energy transfer. The relative fragment ion intensity was analyzed by mass spectrum, for example, at 30 eV, the + + products CH+ 4 , CH3 , and CH2 accounted for 74%, 22%, and 4%, respectively. Not only have Toennies et al. performed the abundant experiments, but they also proposed a mechanism based on vibronic symmetry correlation theory to interpret their experimental results. Their experiments have enriched our knowledge of collision dynamics of the studied system at low energies, and stimulated theoretical investigations on this system. For theory aspect, a complete description for protonmethane collisions requires us to solve the full timedependent Schrödinger equations (TDSE), which is forbidden by immense computational costs. Consequently, some approximations are needed to be made to obtain the reasonable results. Kimura et al.20 have studied electron capture and direct elastic scattering in H+ + CH4 system using a molecular model combined with a molecular-orbital expansion method below 1.5 keV. They employed three incident orientations: the proton approaches along a C–H direction and its opposite direction, and a bisector of a H–C–H bond angle direction. The results from the three orientations were different to some extent, and showed the dynamics of the incident orientation effect. Jacquemin et al.21 have investigated dynamics of proton collisions with methane at 30 eV with the electron nuclear dynamics (END) approach. They described nuclei with classical mechanics and treated electrons with a quantum mechanical method. By employing the six basic orientations, they predicted the possible reaction processes. Total differential cross sections and the fragmentation processes were reported. The average energy transfer has been made in comparison with experimental results. Overall, both simulations displayed valuable mechanistic information which was not accessible by experimental measurements, and helped to deepen our views on proton-molecule collision dynamics. At present, we study the H+ + CH4 collision at 30 eV by employing 9 incident orientations. Our simulations are based on a combination of time-dependent density functional theory (TDDFT)22 and molecular dynamics (MD),23 using the timedependent local density approximation (LDA).24 TDDFTMD model is a non-adiabatic approach which uses classi-

J. Chem. Phys. 140, 054308 (2014)

cal mechanics to describe the motion of the nuclei and employs a quantum representation to characterize electron dynamics. The electron and ion dynamics are calculated simultaneously in real time and real space. It has been applied to investigate the stopping power of proton and antiproton in metals,25, 26 the excitation and ionization of molecules by atom or proton scattering,27–32 and electronic excitation induced by ion collision with graphene.33–35 These studies provided us reasonable qualitative (or semi-quantitative) results. To the end, the well established approach is demonstrated to be feasibility to deal with diverse collision issues. In this work, we are devoted to investigate the mechanism of fragmentation as well as orientation effect in H+ + CH4 collision at low energies. The dependence of scattering properties on initial velocity of both projectile and target is studied. Moreover, we also study the influence of a different self-interaction correction (SIC) scheme on dynamical evolution. The article is organized as follows. In Sec. II, we explain the theoretical framework and the numerical details. In Sec. III, we present and discuss our results compared with other theoretical and experimental results. Finally, the conclusions are drawn in Sec. IV. Atomic units are employed throughout this work, except where specifically mentioned.

II. THEORY AND COMPUTATIONAL DETAILS A. Ehrenfest dynamics within TDDFT

The detailed reviews of TDDFT method can be found in Refs. 36–39, only a brief summary is presented here. In TDDFT, the degrees of freedom are the single-particle wave functions, {ϕ j (r), j = 1, . . . , Nel }, and the positions of the ionic cores, {RI , I = 1, . . . , Nion }. The time-dependent density can be obtained from the single-particle wave function ϕ j (r, t) which satisfies the time-dependent Kohn-Sham (TDKS) equation   1 2 ∂ (1) i ϕj (r, t) = − ∇ + UKS [ρ](r, t) ϕj (r, t), ∂t 2 where the electronic density is ρ(r, t) =

Nel 

|ϕj (r, t)|2 .

(2)

j =1

The KS effective potential UKS [ρ](r, t) is decomposed into three terms:  ρ(r , t) + υxc [ρ](r, t), UKS [ρ](r, t) = υext (r, t) + d 3 r  |r − r | (3) the external potential υ ext (r, t) consists of the ionic core potential and a time-dependent field (i.e., a by-passing ion). The middle term is the Hartree potential. υ xc [ρ](r, t) is the exchange-correlation (xc) potential. Although xc potential has been studied for years, its exact formulation is still unknown and must be approximated for practical application.

054308-3

Gao et al.

J. Chem. Phys. 140, 054308 (2014)

The most widely used one is time-dependent adiabatic LDA (TDALDA), which assumes that the functional only depends on ρ(r, t). In other words, no memory effect40 is included. This functional is local both in space and time, and it should be only available in the case where time-dependent density changes slowly. However, it usually provides reasonable results for various systems even far from this limitation.37, 38, 41 To go beyond the adiabatic approximation, the non-adiabatic functionals are constructed, see Refs. 42–46. The equations of motion for ions are determined by variation with respect to the ionic variables RI (position) and PI (momentum)  ∂2 MI 2 RI (t) = − drρ(r, t)∇RI υne (r − RI ) ∂t  ZI ZJ − υext,ion (RI , t), − ∇RI |RI (t) − RJ (t)| I =J (4) where MI and ZI represent the mass and the charge of ion number I. υ ne denotes the electron-nucleus potential. υ ext, ion is the interaction between ions and the external field. From Eq. (4), it shows that the ionic movement could be obtained by the time-dependent electronic density ρ(r, t). This fact is the foundation of the framework of TDDFT based on Ehrenfest MD (E-TDDFT). It should be noticeable that the non-adiabatic effect in current strategy arises from the non-adiabatic coupling between ions and electrons. Generally there are two methods to describe the non-adiabaticity of ion-electron coupling: surface-hopping47, 48 and Ehrenfest dynamics.49 Surface-hopping method and its derivatives50 treat electronic transitions with a hopping probability between different excited-states. Because it requires an explicit computation of adiabatic time-dependent excited-state, it is prohibitive for large system. While, present Ehrenfest approach is a mean-field treatment where ionic dynamics evolves on an averaged potential surface, and it represents the wavefunction as a superposition state instead of exact calculations for excited-states, thus potential, force, and electronic density are calculated as mean expectation values. Even though Ehrenfest method may lose its physical meaning in asymptotic limit, it could still give acceptable results for describing non-adiabatic collisions.

B. Self-interaction correction

As shown in Eq. (3), the Hartree potential includes the self-interaction (SI) of electrons, which will induce a wrong asymptotics of the electronic potential. The first selfinteraction correction scheme is proposed by Perdew and Zunger, which subtracts SI contributions from the total energy:51 UαSI C = UKS − [UC (ρα ) + Uxc (ρα )].

(5)

The net result is that the obtained Hamiltonian depends on the orbital on which it acts. The non-hermiticity of Hamiltonian

can therefore induce a violation of orthonormality in a time propagation. Many other SIC ideas have been proposed. For a recent review on this topic, see the literature.52 Among these schemes, we consider in this paper only two53 of them, that is, Slater SIC and the average-density SIC (ADSIC):    ρi Slater SI C U (r) , = (6) Uα ρ i i ρ   ρ

+ Uxc UαADSI C = UKS − UC ,0 . N N

(7)

It is clear that Slater SIC in Eq. (6) first constructs potentials state by state, then averages them by taking the weight of density ratio. This state-dependence potential computation will increase the computational cost. While ADSIC treated the SI contribution of each electron on the same footing by averaging the SIC term in Eq. (7) over all valence electrons. It can thus be interpreted as an approximation to the PerdewZunger one. Moreover, ADSIC presents numerous formal and numerical advantages, one of them rendering the Hamiltonian state-independent. C. Numerical details

The explicit numerical method is introduced elsewhere.54, 55 In calculations, we use ALDA functional with the parameterization by Perdew and Wang56 complemented by ADSIC. Pseudopotentials57 are employed to describe ionic core potential, which consists of a local part and a non-local one. The ground state is optimized by the simulated annealing method and the damped gradient method. For propagating electronic wave-function, the time-splitting method58, 59 is used. The ionic motion is solved by using Verlet algorithm.60 In order to describe electronic excitation, the total number of escaped electrons is expressed as Nesc (t) = N (t = 0) − V drρ(r, t). To avoid escaped electrons being reflected back to the numerical box when they reach the boundary, absorbing boundary conditions (ABS)61 are used. The ionization probabilities P(k+) (t) are evaluated by an approximate method,62 and because no sizable contributions from k > 1, the high charge states are omitted here. All above numerical techniques are implemented in code package TELEMAN.54, 55 In order to describe the KohnSham orbitals, densities, and potentials, the three-dimensional Cartesian coordinate system is represented by 72 × 72 × 72 grids with a spacing of 0.412 a0 . The equilibrium geometry of methane is calculated. The present equilibrium C–H bond length is 1.04 Å and the angle H–C–H is 109.5◦ , which are in good agreement with previous theoretical and experimental results.63 The calculated orbital energies, ε(1t2 )=14.58 eV and (2a1 )=22.34 eV, compare well with photoelectron spectra experiments64 within a 3% discrepancy. The initial collision configuration is displayed in Fig. 1, where carbon, hydrogen, and proton are labelled by C, Hi (i = 1, 2, 3, 4), and p+ , respectively. For convenience of description, some special points are introduced (see Fig. 1(a)):

054308-4

Gao et al.

J. Chem. Phys. 140, 054308 (2014)

FIG. 1. The initial geometry of H+ + CH4 system. Black, blue, and red spheres represent carbon, hydrogen, and proton, respectively. The three main configurations, (a) Corner, (b) Face, and (c) Edge, generate 9 initial directions of attack for the projectile labelled by I, II, and III. O denotes the center of triangle H2H3H4 and A represents the midpoint of H2H3. θ denotes scattering angle. The center of mass of targets is initially in the origin.

O is the center of triangle H2H3H4 and A denotes the midpoint of H2H3. Obviously, there are six identical triangular pyramids H1OH2A, which is supposed to be the minimal unit cell. Therefore, 9 initial orientations are taken as follows. For Corner collisions, see Fig. 1(a), the proton attacks the target parallel to H1O, the orientation where impact parameter b increases along OA, OH2 and their bisector corresponds to corners I, II, and III. For Face cases in Fig. 1(b), the proton moves to the target along OH1, the orientation where b increases along OH2, OA and their bisector is faces I, II, and III. For Edge cases, see Fig. 1(c), the proton approaches the target along AC, the orientation where b increases along H4H1, AH2 and their bisector is edges I, II, and III. It is worthwhile to note that the orientations used in Ref. 20 correspond to corner II, face I, and edge I, while the configurations considered in Ref. 21 are represented by corner I, corner II, face I, face II, edge I, and edge II in the present study. In experiments, the methane molecule is probably in vibrational and rotational states at the initial moment due to a non-zero temperature, which might produce at least 9 vibrational modes.65 To calculate collision dynamics with these initial states, a large number of events are required. However, this is beyond the topic of present study. Methane is initially at rest with the center of mass in the origin at T = 0 K, since nuclei is treated as classical particles with non-zero point vibrations in our theoretical model. The proton is placed at 25 a0 distant from the methane, and it is given initial velocity corresponding to the collision energy. b varies from 0 a0 to 5 a0 with b = 0.1 a0 and from 5 a0 to 10 a0 with b = 0.5 a0 . The total simulation time is 40 fs with t = 0.6 as. The small time step could ensure the stability of time propagation in the calculations. For a few representative trajectories, longer evolution time (i.e., 100 fs) is required. All the dynamical collision processes are considered in the laboratory system in present simulation.

III. RESULTS AND DISCUSSION

For the dynamical process we first performed the numerical check. As shown in Fig. 2, the time evolution of en-

FIG. 2. The kinetic, potential, and total energies in eV unit (left axis), and the relative distance between atoms (right axis) as a function of time. The initial direction corner I with b = 1.9 a0 for H+ + CH4 collision at 30 eV. The vertical dashed line represents the time of energy change.

ergies and the relative distance between atoms are plotted from the direction corner I with b = 1.9 a0 . Within 15 fs, the distance between C and proton decreases monotonically as time increases, and bond lengths of CH4 have a slight oscillation by the weak proton-methane interaction. When the proton approaches close to CH4 at 15 fs, a portion of potential energy (Epot ) is converted into the kinetic energy (Ekin ). At this moment, the bond length starts to increase obviously, and an inter-ion oscillation is motivated. The closest point appears at 17.5 fs. After this point, the proton leaves the system, the Coulomb attraction is dominant and decelerates the proton, thus it transforms a fraction of kinetic energy (Ekin ) into the potential energy (Epot ) until 20 fs. After collision, bond lengths exhibit harmonic oscillations. Although there occurs energy transfer between Ekin and Epot during the entire collision process, one can observe the invariance of the total energy (Etot ) except that the proton reached the boundary at 33.8 fs. The reason that Etot increases at 33.8 fs is that the absorbing boundary absorbed the electrons carried off by the proton. In conclusion, the law of conservation of energy is satisfied in collision process.

A. Reaction channel

At low energies, the collision of H+ + CH4 exhibits a variety of reaction channels, for example, dissociation, proton exchange, and molecular formation. The information on fragmentation of CH4 is included in these reactive processes. Hence it is beneficial to understand the reactive dynamics of collision in low impact energy regions. For the sake of comparing with experiments and other calculations, we mainly focus on collision energy at 30 eV in present simulations. The predicted reaction channels at 30 eV are shown below:

054308-5

Gao et al.

p+ + CH4 ⎧ (a) Non-charged transfer (NCT), ⎪ ⎪ ⎪ p+ + CH4 , ⎪ ⎪ ⎪ ⎪ ⎪ (b) Charge transfer (CT), ⎪ ⎪ ⎪ ⎪ H0 + [CH4 ]+ , ⎪ ⎪ ⎪ ⎪ (c) Protonation (Pro), ⎪ ⎪ ⎪ ⎪ [CH4 p]+ , ⎪ ⎪ ⎪ ⎪ (d) Proton exchange (Ex), ⎪ ⎪ ⎪ Hq + [CH p]q  , q + q  = +1, ⎪ 3 ⎨ dissociation (Dis), (8) → (e) Collision-induced ⎪ q  q q   ⎪ p + H + [CH ] , q + q + q = +1, ⎪ 3 ⎪ ⎪ ⎪ (f) Hydrogen molecular formation (For), ⎪ ⎪  ⎪ ⎪ [Hp]q + [CH3 ]q , q + q  = +1, ⎪ ⎪ ⎪ ⎪ (g) Proton exchange with dissociation (Ex + Dis), ⎪ ⎪ ⎪   ⎪ ⎪ Hq + Hq + [CH2 p]q , q + q  + q  = +1, ⎪ ⎪ ⎪ ⎪ (h) Hydrogen molecular formation with ⎪ ⎪ ⎪ ⎪ dissociation (For + Dis), ⎪ ⎩   [Hp]q + Hq + [CH2 ]q , q + q  + q  = +1, where q, q , and q represent the charge state of ions. For present collision system, we have demonstrated that the pure ionization, i.e., direct electron emission, can hardly happen at such a low collision energy.66 The non-charged transfer (NCT) process is an inelastic scattering process where the proton is scattered away while exciting the vibrational degree of freedom of methane (see Fig. 2). Reactions (d) to (h) also contribute to charge transfer process since the outgoing parti cle Hq (or Hq ) becomes q(q )= 0 with a finite probability. One should keep in mind that the non-reactive process is mostly a hybrid collision process between NCT and CT with a respective probability for each process.28

J. Chem. Phys. 140, 054308 (2014)

In Fig. 3, the reactive probability is shown by comparing with other theoretical calculation results.21 Apparently, we obtained similar reaction type as Ref. 21, as one can see from Figs. 3(a)–3(e). Moreover, we observed another two new reactions with a small probability (see in Fig. 3(f)): protonation (Pro) and proton exchange with dissociation (Ex+Dis). Pro reaction yields CH+ 5 complex, while Ex+Dis reaction is responsible for Ex reaction transition from H1 to H2 site as analyzed below (see Sec. III B). For NCT/CT in Fig. 3(a), we predicted larger probability for impact parameter 0.8 < b < 2.8 a0 , while for 0 < b < 0.8 a0 the case is reversed. Nevertheless, both agree well for larger impact parameter (b > 3.4 a0 ), in this impact parameter region no reactions happen. For Dis reaction (see Fig. 3(b)), we obtained smaller probability except some points in considered impact parameter region. For Ex reaction in Fig. 3(c), the results are completely reversed comparing with Dis case. In Fig. 3(d), Jacquemin’s For reaction results are dominant over our ones. Note that both predicted a small probability for For+Dis reaction shown in Fig. 3(e). The above discrepancy is ascribed to the different treatment of nuclear-electron coupling. In Ref. 21, Jacquemin et al. used the extended Lagrangian molecular dynamics to describe collision process, where a single electronic potential energy surface (PES) determines the dynamics. This treatment is obviously an adiabatic description. However, the collision problems typically involve several PESs, which requires a non-adiabatic approach. Our present study used Ehrenfest dynamics which considers the nuclear-electron coupling by an averaged multipotential surface. Although it is still a single PES, the nuclear-electron non-adiabatic coupling is indeed involved, especially for small impact parameters

FIG. 3. ((a)–(f)) The reaction probability of all channels is compared with others’ calculation.21 The obtained reactions are non-charged transfer (NCT), charge transfer (CT), proton exchange (Ex), collision-induced dissociation (Dis), proton exchange with dissociation (Ex+Dis), protonation (Pro), hydrogen molecular formation (For), and hydrogen molecular formation with dissociation (For+Dis).

054308-6

Gao et al.

FIG. 4. The fragment yield is calculated to compare with experimental and previous theoretical results.

(b < 3.0 a0 ). The applicability of the present method can be also found when we calculated the intensity of reaction product displayed in Fig. 4. On the whole, we obtain the same fragment species, + + 19 and othCH+ 4 , CH3 , and CH2 , with experimental results 21 ers’ calculations. The estimated fragment intensities are + + 79.6%, 20.1%, and 0.2% for CH+ 4 , CH3 , and CH2 , respectively, which is in agreement with experimental values + + (CH+ 4 :CH3 :CH2 = 74%:22%:4%), while results from Ref. + + 21 (CH4 :CH3 :CH+ 2 = 49%:50%:1%) drastically deviated from ones from experiments. This indicates that the intensity of CH+ 4 was underestimated by employing the adiabatic PES, and in this case that of CH+ 3 was overestimated. The results are consistent with what has been shown in Fig. 3. Moreover, there is a small intensity peak at 17 amu in experimental measurement, which is also found by present study instead of in calculations of Ref. 21. However, the origin of this small peak is different. It was considered as an isotope effect (13 CH4 /12 CH4 ) in experiment,19 while we conclude that CH+ 5 fragment may also contribute to this small intensity peak. B. Fragmentation process

In this section, we shall elucidate the mechanism of fragmentation at 30 eV. We have studied the proton exchange and dissociation processes by means of ionic trajectory and the energy loss of proton in detail in Ref. 31, for the case of H+ + CH4 collision, we will use similar analysis method to explain the patterns of fragmentation. In order to distinguish the same reaction type occurs at different hydrogen sites (H1 or H2), we have to define a special label for channels, that is “abbreviation + number,” such as Dis1 means dissociation results from C–H1 bond breaking and Ex1+Dis2 indicates that proton exchange occurs at H1 site and meanwhile the cleavage of C–H2 bond happens. The refined channels ensure us a better understanding of reaction dynamics. In Fig. 3, it is remarkable that there are two main fragmentation patterns: collision-induced dissociation (Dis) and proton exchange (Ex) processes. We find that the reaction transition from H1 to H2 for each pattern occurs in some directions, such as Dis1 to Dis2 in corner II and Ex1 to Ex2 in

J. Chem. Phys. 140, 054308 (2014)

face I. It is instructive to understand the reactive mechanism of these two reaction transitions. With regard to dissociative processes, we take the cases in direction edge I at b = 0.7, 1.1, and 1.5 a0 . The bond lengths and the corresponding kinetic energies are given in Fig. 5. For b = 0.7 a0 , the cleavage of the C–H1 bond occurs when the proton moves to the closest distance, while as seen from Fig. 5(b), the energy loss of proton is 17.4 eV and H1 hydrogen obtains an energy of 9.6 eV. After the closest point, as time increases, the distance of C–p and C–H1 increases monotonically and other atoms oscillate weakly around the equilibrium position. Although H1 and proton are scattered, both are decelerated under the influence of Coulomb attraction with the asymptotic values 2.84 eV and 2.26 eV for H1 and proton, respectively. In terms of the definition for reaction type, this belongs to the Dis1 reaction channel. In Fig. 5(e), there is no difference with Dis1 process except that the bond breaking occurs at C–H2 instead of C–H1 bond, thus it is a Dis2 reaction channel. Note that in this case H2 and proton escape with 14.84 eV and 3.5 eV, respectively, which is much larger than in the Dis1 case. In edge I, the dissociation from H1 to H2 (Dis1 to Dis2), an intermediate NCT/CT process occurs, which is depicted in Figs. 5(c) and 5(d). In Fig. 5(c), there is no bond breaking, although the bond vibration is much more violent than that in Dis1 and Dis2 cases. This is due to the amount of energy transferred from the proton to methane. In Fig. 5(d), at the closest distance, H1 and H2 obtain an energy of 7.2 eV and 9.4 eV, respectively, which do not induce the bond breaking of C–H1 and C–H2, albeit the obtained energies are much larger than the chemical dissociation energy of methane:19 CH4 → CH3 + H, Ed = 4.4 eV. In fact, we observe an intra-molecule energy transfer from H1 and H2 to H3 and H4, which permits the stabilization of methane with causing an anti-phase oscillation between C–H1(H2) and C–H3(H4) in Fig. 5(c). From Dis1 to NCT/CT to Dis2 channel, energy transferred mode from proton to H1 (Dis1 case) first becomes that from proton to H1 and H2, and H2 gets more energies than H1 (NCT/CT case), finally becomes that from proton to H2 (Dis2 case). Similar Dis1 → NCT/CT → Dis2 reaction transition is also found in corners II and III. As a result, NCT/CT channels appear in small impact parameter region is supposed to be a necessary reaction path for dissociation from H1 to H2. Another significant reaction transition is the proton exchange (Ex) from H1 to H2, which is reported in Fig. 6 at b = 0.3, 0.6, and 0.9 a0 in direction face I. In Fig. 6(a), the projectile proton is caught by methane at closest collision point with losing considerable energy, and then H1 hydrogen obtains 26.6 eV energy (see Fig. 6(b)) to escape the system. There is a time delay between proton capture and H1 emission, which is important for the energy transfer process from proton to H1. During H1 scattering process, it loses another 4.9 eV energy influenced by the Coulomb attraction interaction and finally gains 21.7 eV kinetic energy. In the meantime, other atoms are excited in vibration modes around the equilibrium position. Apparently, it belongs to Ex1 process. For the case of the Ex2 process in Fig. 6(e), proton capture and H2 hydrogen emission almost happen at the same time, while H2 hydrogen obtains 4.5 eV energy in Fig. 6(f).

054308-7

Gao et al.

J. Chem. Phys. 140, 054308 (2014)

FIG. 5. Dissociation process (Dis) transition from H1 to H2 in direction edge I at 30 eV. Panels (a) and (b) illustrate a dissociation process that happens at H1 site (Dis1) at impact parameter b = 0.7 a0 , panels (c) and (d) indicate a NCT/CT process at impact parameter b = 1.1 a0 , and panels (e) and (f) show a dissociation process that happens at H2 (Dis2) at impact parameter b = 1.5 a0 . Panels (a), (c), and (e) illustrate the relative atomic distance vs. time, while panels (b), (d), and (f) show the time evolution of kinetic energies.

Therefore, the longer time delay between proton capture and H1 (or H2) emission promotes outgoing particle (H1 or H2) acquiring more energy. After the closest point, other atoms strongly oscillate with a relatively large amplitude, influenced by the scattering particle at lower velocity (3.0 eV). In face I, proton exchange from H1 to H2 (Ex1 to Ex2), there is an intermediate Ex1+Dis2 reaction process shown in Figs. 6(c) and 6(d). In Fig. 6(c), when the proton moves to the closest distance with methane, proton exchange with H1 and C–H2 bond dissociation occurs. During the H1 and H2 scattering process, they are influenced by each other. This is the reason that there is no asymptotic energy for outgoing particles in Fig. 6(d). From Ex1 to Ex1+Dis2 reactions, the way of energy transferred is similar with Dis1 to NCT/CT, but C–H2 dissociated due to no intra-molecule energy transfer. Similar Ex1 → Ex1+Dis2 → Ex2 reaction transition is also found in corner II. Consequently, Ex+Dis reaction is associated with the transition for proton exchange from H1 to H2. From the results of others’ work,67 CH+ 5 was believed to be the intermediate product of reaction H+ + CH4 + → CH+ 5 → CH3 + H2 . In Fig. 3, the protonated methane appears only at a couple of impact parameters in edges I and II, which is consistent with theoretical calculations.68, 69 Nevertheless, we do not obtain the final separate products (or secondary products) from protonation, because the CH+ 5 ion

reached the absorbing boundary with a small translational velocity. In fact, the protonation process is exactly a part of another reaction process For1+Dis2 (such as 0.4 a0 in corner II and 0.3 a0 in face III) as discussed below. Fig. 7(a) shows the relative atomic distance with kinetic energies vs. time from the initial direction corner II/b = 0.4 a0 . When the proton approaches the methane molecule at the closest distance, it mostly transfers its kinetic energy completely to H1, and embeds itself into the methane molecule 19 to form a strongly distorted CH+ 5 quasimolecular complex. Note that in a very short time, the kinetic energy of H1 is converted into the potential energy of system. Then the motion of system is similar to the protonation process. However, the complex is labile, the proton makes bonding with H1 hydrogen at 61 fs and then escapes the system with a small translational energy. The remaining methyl group (CH+ 3 ) has an energy fluctuation, which results in the dissociation of C–H2 bond. This dynamical process holds for For2+Dis1 channel which occurs at 0.2 a0 in face III. As for For2 channels from the initial direction edge II/b = 1.2 a0 displayed in Fig. 7(b), the approaching proton loses about 18.7 eV at the closest distance, and directly runs away from system together with H2 (hydrogen molecule) at some certain velocity. This reactive dynamics is the same as the manifested F channels in the calculation performed by Jacquemin et al.21 In their study, they

054308-8

Gao et al.

J. Chem. Phys. 140, 054308 (2014)

FIG. 6. Proton exchange process (Ex) transition from H1 to H2 in direction face I at 30 eV. Panels (a) and (b) illustrate a proton exchange process that happens at H1 site (Ex1) at impact parameter b = 0.3 a0 , panels (c) and (d) indicate a proton exchange with dissociation (Ex1+Dis2) process at impact parameter b = 0.6 a0 , and panels (e) and (f) show a proton exchange process that happens at H2 (Ex2) at impact parameter b = 0.9 a0 . Panels (a), (c), and (e) illustrate the relative atomic distance vs. time, while panels (b), (d), and (f) show the time evolution of kinetic energies.

observed subsequent dissociation of For channels, which was the principal way to the formation of For+Dis in Ref. 21. However, this subsequent reactions are not observed in the

present study. In our simulation the For+Dis channel results from the dissociation of CH+ 5 which plays the role of an intermediate reaction product. The For channel derives from the direct bond making while leaving a stable methyl group. C. The scattering properties

The scattering angle θ is defined as the projectile deflected from initial momentum direction (see Fig. 1). The expression70 reads 1/2  2 Px + Py 2 sin θ = , (9) |P|

FIG. 7. The relative atomic distance and kinetic energies vs. time of simulation of H+ + CH4 at 30 eV. (a) Simulation of hydrogen molecular formation with dissociation (For+Dis) reaction from the initial direction corner II and impact parameter b = 0.4 a0 . This is For1+Dis2 reaction channel. (b) Simulation of hydrogen molecular formation (For) reaction from the initial direction edge II and impact parameter b = 1.2 a0 . This indicates For2 reaction channel. The insets show the evolution of kinetic energies for the system corresponding to each case.

where Px and Py are x and y components of the final momentum P of the outgoing particle. The scattering angle θ as a function of b for 9 investigated incident directions at 30 eV is plotted in Fig. 8. For these 9 orientations, there is at least one extremum (maximum) called primary rainbow angle (θr1 ) which occurs at relatively large impact parameter br1 where dθ | 1 = 0. The primary rainbow angle is the feature of iondb b=br molecule scattering and it refers to the most obvious attractive interaction. For b < br1 , θ decreases smoothly to the glory angle (θ g ) where θg (b)|b=bg = 0◦ , which signifies the projectile has no deflection. However, corner III and face III orientations exhibit the secondary rainbow angle (θr2 ) instead of the

054308-9

Gao et al.

J. Chem. Phys. 140, 054308 (2014)

FIG. 8. Scattering angle θ vs. b at 30 eV corresponding to the 9 orientations of Fig. 1. The position of glory and rainbow angle is observed in the range of 1.8 ≤ bg ≤ 3.8 a0 and 2.9 ≤ br ≤ 4.6 a0 , respectively.

glory angle. The reason is that the scattering angle arises from the net effect of the Coulomb repulsion and attraction during the collision process.31 Only when the net effect of the force is zero the glory angle can be observed. Although, there are different systems, where the secondary rainbow angles are detected in experiments for H+ + (H2 ,71 CF4 19 ) collisions and found in H+ + (H2 ,72, 73 C2 H2 ,74 HF75 ) theoretical calculations. No secondary rainbow angle is discerned in six directions’ calculations21 of H+ + CH4 collision. Comparing with present results, this indicates that adding more orientations helps to discover more detailed microscopic scattering patterns. The position of glory and rainbow angle is 1.8 ≤ bg ≤ 3.8 a0 and 2.9 ≤ br ≤ 4.6 a0 , which is in accordance with 1.8 ≤ bg ≤ 3.5 a0 and 3.0 ≤ br ≤ 4.3 a0 for the glory and rainbow angle results in Ref. 21. Furthermore, the mean primary rainbow angle in present calculations is 9.6◦ , which is in nice agreement with experimental results of 10◦19 and calculation values of 12◦ .21 We have calculated the total charge transfer cross section76, 77 by the equation  bmax CT σ = 2π bP (b)db, (10) 0

where P(b) is the charge transfer probability. The averaged result of this integral equation is 10.87 Å2 , which is much larger than the calculation result of END 3.3 Å2 and phenomenological formula result 5.8 Å278 based on experimental data at keV energies. However, the result from available semi-empirical expression79 recommended by fusion plasma simulations is 28.55 Å2 , which is roughly a factor of 3 larger than present result. In order to further compare with experiments19 as well as previous theoretical calculations,21 the averaged energy loss of the proton vs. θ is displayed in Fig. 9. Comparing with END calculations, we obtain less reasonable results, since the present results underestimate the results at larger angles. The discrepancy between experiments and our calculations may result from the theoretical defect. Proton-molecule collisions are full quantum problems at a microscopic level and

FIG. 9. The averaged energy loss of proton vs. θ at 30 eV along with experimental and END results.

should include all freedoms of nuclei and electrons. In that sense, E-TDDFT is by no means the accurate method, and it may produce some artifacts due to the mean-field treatment, especially in asymptotic range. On the other hand, we also made some approximations, such as the initial vibration of methane is frozen and the density changes slowly. In addition, for experimental measurement the signal noise increases with scattering angle decreasing due to energy resolution, anyway, there should be error bar for the results. Although present results agree poorly with experiments, it still reveals a general trend that for a given energy, larger scattering angle, originating from a strong projectile-molecule interaction, would lead to the more energy loss of proton. For this angle range, the largest energy loss is less than 2.0 eV. D. The dependence on initial velocity

1. The initial velocity of projectile

First, we investigated the energy dependence of scattering properties for proton energies 20, 30, and 40 eV. Note that in this case the target is initially at rest. The scattering angle and energy loss as a function of the impact parameter are shown in Fig. 10 for the orientation edge III. In Fig. 10(a), for a given impact parameter, the scattering angle decreases with increasing energy, but the profiles of the three curves are very similar. The large angle scattering occurs for small impact parameters, in particular at b = 0 a0 the proton directly bounces back with θ = 180◦ , since the core repulsive force is dominant in this impact parameter range. With increasing the impact parameter, the scattering angle decreases smoothly to the minimum θ min which results from the net effect of the Coulomb repulsion and attraction. The extreme case is that the net effect of the force is zero. In this case, θ min equals zero which becomes the glory angle θ g . In fact, θ min is the transition point of the interaction force. The minimum angle divides the scattering area into two different parts: the repulsion and the attraction. Through the minimum, the scattering angle increases gradually to the primary rainbow angle (θ r ). The position of the rainbow angle,

054308-10

Gao et al.

J. Chem. Phys. 140, 054308 (2014)

and is almost equivalent at 30 and 40 eV. The transition point of the energy-loss happens to be at b = 1.1 a0 , which is supported by the fact as follows. If the C–H1 bond (∼1.96 a0 ) is projected into the orientation along which the impact parameter increases (y axis), as one can see in the inset, the value is exactly 1.1 a0 . On the basis of the amount of energy loss of the proton, the small b collisions are more violent than larger ones.

2. The initial velocity of target

Here, we considered a “special” initial state of methane, where all atoms move along the same direction with a small initial velocity, see the collision configuration shown in the inset of Fig. 11. The proton is given an initial kinetic energy of 30 eV (υ 0pro = 1.43 a0 /fs), and the initial velocity of methane varies as 0, 0.01υ 0pro , and 0.02υ 0pro . The results are plotted from the orientation edge II. In Fig. 11(a), for a given b in small impact parameter range (b < 3 a0 ), with increasing methane velocity υ 0 , the scattering angle first increases and then decreases. It indicates that for smaller υ 0 the proton-methane interaction will slightly become stronger due to the motion of methane. For

FIG. 10. The orientation edge III (see Fig. 1(c)) at impact energy 20, 30, and 40 eV. Methane molecule is initially at rest. (a) Scattering angle θ vs. b. The position of the rainbow angle, 20 eV (3.2 a0 , 17.3◦ ), 30 eV (3.1 a0 , 11.8◦ ), and 40 eV (3.0 a0 , 8.8◦ ). (b) The energy loss of proton vs. b. The geometry projection of C–H1 bond is depicted in the inset to interpret the appearance of transition point at b = 1.1 a0 .

20 eV (3.2 a0 , 17.3◦ ), 30 eV (3.1 a0 , 11.8◦ ), and 40 eV (3.0 a0 , 8.8◦ ), has a shift toward the smaller impact parameter with energy increasing. This rainbow angle shift was also found in the H+ + CF4 collision.80 The larger energy indicates a shorter interaction time. In order to obtain the same attractive effect, the new rainbow angle could be observed only when the proton impacts closer to the target. For b > br , the scattering angle becomes smaller and smaller. During the collision, although the proton is accelerated in a certain period, it has ultimately been observed as an energyloss behavior (E) as shown in Fig. 10(b). Basically, the energy loss decreases with increasing impact parameter regardless of the collision energy. For small impact parameters (b ≤ 1.1 a0 ), the energy loss of proton is directly proportional to the collision energy. The largest energy loss occurs at b = 0 a0 , where a head-on collision happens (see the inset). Emax is 7.0, 9.89, and 12.76 eV for 20, 30, and 40 eV, respectively. With increasing the impact parameter, proton energy loss decreases steeply. Note that for b = 1.1 a0 the energy loss becomes the same (E = 3.4 eV) for the three impact energies. For b > 1.1 a0 , the proton energy loss is the largest at 20 eV

FIG. 11. Collisions for edge II at impact energy of 30 eV, and methane is given to an initial velocity. (a) Scattering angle as a function of b. (b) The energy loss of proton vs. b. The collision configuration is shown in the inset.

054308-11

Gao et al.

larger υ 0 the scattering angle decreases because of the short interaction time, which is similar with that of Fig. 10(a). It is noticeable that the glory angle disappears when increasing υ 0 , and the minimum angle is shifted. For rainbow angle, it does not change the position when υ 0 increases to 0.01υ 0pro , and it finally has a shift toward the smaller b when increases to 0.02υ 0pro , which is also consistent with Fig. 10(a). For b > 4 a0 , the scattering angle does not depend on initial methane velocity strongly. The energy loss in b < 3 a0 range is in agreement with scattering angle variation as shown in Fig. 11(b). The transition point occurs around the glory angle, where E = 2.12 eV. For b > 4 a0 , the proton energy loss is almost the same for three initial υ 0 . We found that the energy loss behavior in this case is completely different with that in Fig. 10(b). In general, it implies that initial methane velocity has an obvious effect for scattering angle and energy loss for small impact parameter, while for larger b (b > br ) no difference is observed. E. Comments on SIC in dynamical process

In order to study the sensitivity of xc functionals in dynamical collision, we performed the calculation on LDA level, and compared it with results from some choices of SIC: Slater and ADSIC, see Fig. 12. The collision orientation is the same as Fig. 11. For b < 3.1 a0 , we found that two SIC results are totally different comparing with LDA. Slater SIC exhibits larger scattering angle which corresponds to larger energy loss (see the inset), while ADSIC means smaller scattering angle due to weak interaction. Note that the glory angle occurs in ADSIC case. For b > 3.5 a0 , both SIC have a stronger scattering than LDA, and ADSIC has the largest rainbow angle. Concerning the position of rainbow angle, for ADSIC it occurs the same b as LDA, but for Slater it is off the LDA case. This is also observed in energy loss. In this sense, ADSIC seems to work better than Slater, and has reasonable correction for LDA. For further study, other SIC choices should be also considered, and SIC effect in dynamical collision needs to be checked in other collision systems.

J. Chem. Phys. 140, 054308 (2014)

IV. CONCLUSIONS

By employing TDDFT based on Ehrenfest dynamics approach, we investigated the mechanism of fragmentation in H+ + CH4 collision at low energies, and studied the initial velocity dependence of scattering properties. In this paper, all possible reactions are presented from 9 incident orientations at 30 eV. The present fragment intensity is found to be closer to experimental results than others’ calculations, the discrepancy is explained by nuclear-electron coupling. The mechanism of reaction transition is attributed to the intra-molecule energy transfer. We also found that the protonation product CH+ 5 plays the role of an intermediate reaction product for the For+Dis channel, and the For channels derive from the direct bond making while leaving a stable methyl group. The energy loss of the proton is in poor agreement with experimental results. The dependence of scattering angle and energy loss on initial velocity of proton is different from that of methane. Finally, different SIC for dynamical evolution may induce different results, but with respect to the position of rainbow angle, ADSIC seems to have reasonable correction for LDA. In conclusion, present E-TDDFT works better than the adiabatic method on most aspects, and it also gives more fruitful results. Nevertheless, E-TDDFT is a mean-field approach, and it averaged excited-states potential to a single potential surface, thus it may obtain some unphysical results for the case where different states would induce different reactions, such as fragmentation. On the other hand, the non-adiabatic effect from xc functionals is absent, and the temperature effect of collision is not clarified. In this sense present results give us a qualitative picture of collision dynamics. Besides, more advanced treatment for non-adiabaticity of xc functional and nuclear-electron coupling should be welcome in the future work, but for the time being, to the best of our knowledge, it still constitutes a great numerical challenge.

ACKNOWLEDGMENTS

We thank E. Suraud, P. M. Dinh, and F. Calvayrac for many constructive discussions. C.-Z.G. thanks the financial support from China Scholarship Council (CSC) (No. [2013]3009). This work was supported by the National Natural Science Foundation of China under Grant Nos. 11025524 and 11161130520, the National Basic Research Program of China under Grant No. 2010CB832903, and the European Commissions 7th Framework Programme (FP7-PEOPLE2010-IRSES) under Grant Agreement Project No. 269131. 1 C.

FIG. 12. The sensitivity of xc functionals for dynamical collision. The comparison is made between LDA, Slater, and ADSIC.

F. Barnett, J. A. Ray, E. Ricci, M. I. Wilker, E. W. McDaniel, E. W. Thomas, and H. B. Gilbody, Atomic Data for Controlled Fusion Research (National Laboratory, Oak Ridge, Tennessee, 1977), Vol. 1. 2 J. G. Collins and P. Kebarle, J. Chem. Phys. 46(3), 1082 (1967). 3 L. H. Toburen, M. Y. Nakai, and R. A. Langley, Phys. Rev. 171(1), 114 (1968). 4 W. E. Wilson and L. H. Toburen, Phys. Rev. A 11(4), 1303 (1975). 5 W. R. Gentry, H. Udseth, and C. F. Giese, Chem. Phys. Lett. 36, 671 (1975). 6 H. Tawara, At. Data Nucl. Data Tables 22, 491 (1978). 7 H. Knudsen, U. Mikkelsen, K. Paludan, K. Kirsebom, S. P. Moller, E. Uggerhoj, J. Slevin, M. Charlton, and E. Morenzoni, J. Phys. B 28, 3569 (1995). 8 F. Gobet, B. Farizon, M. Farizon, M. J. Gaillard, M. Carré, M. Lezius, P. Scheier, and T. D. Märk, Phys. Rev. Lett. 86, 3751 (2001).

054308-12

Gao et al.

J. Chem. Phys. 140, 054308 (2014)

9 R. K. Janev, Atomic and Molecular Process in Fusion Edge Plasmas, edited

43 G.

by R. K. Janev (Plenum, New York, 1995), pp. 1–13. 10 M. Burgdorf, G. Orton, J. van Cleve, V. Meadows, and J. Houck, Icarus 184, 634 (2006). 11 M. R. Swain, G. Vasisht, and G. Tinetti, Nature 452, 329 (2008). 12 D. W. Koopman, J. Chem. Phys. 49, 5203 (1968). 13 R. J. McNeal, J. Chem. Phys. 53(11), 4308 (1970). 14 M. E. Rudd, R. D. DuBois, L. H. Toburen, C. A. Ratcliffe, and T. V. Goffe, Phys. Rev. A 28, 3244 (1983). 15 R. S. Gao, L. K. Johnson, C. L. Hakes, K. A. Smith, and R. F. Stebbings, Phys. Rev. A 41, 5929 (1990). 16 J. M. Sanders, S. L. Varghese, C. H. Fleming, and G. A. Soosai, J. Phys. B 36, 3835 (2003). 17 T. Ellenbroek and J. P. Toennies, Chem. Phys. 71, 309 (1982). 18 T. Ellenbroek, U. Gierz, M. Noll, and J. P. Toennies, J. Phys. Chem. 86, 1153 (1982). 19 Y. N. Chiu, B. Friedrich, W. Maring, G. Niedner, M. Noll, and J. P. Toennies, J. Chem. Phys. 88, 6814 (1988). 20 M. Kimura, Y. Li, G. Hirsch, and R. J. Buenker, Phys. Rev. A 52, 1196 (1995). 21 D. Jacquemin, J. A. Morales, E. Deumens, and Y. Öhrn, J. Chem. Phys. 107, 6146 (1997). 22 E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984). 23 O. Sugino and Y. Miyamoto, Phys. Rev. B 59, 2579 (1999). 24 A. Zangwill and P. Soven, Phys. Rev. Lett. 45, 204 (1980). 25 M. Quijada, A. G. Borisov, I. Nagy, R. Díez Muiño, and P. M. Echenique, Phys. Rev. A 75, 042902 (2007). 26 J. M. Pruneda, D. Sánchez-Portal, A. Arnau, J. I. Juaristi, and E. Artacho, Phys. Rev. Lett. 99, 235501 (2007). 27 Z. P. Wang, P. M. Dinh, P.-G. Reinhard, E. Suraud, G. Bruny, C. Montano, S. Feil, S. Eden, H. Abdoul-Carime, B. Farizon, M. Farizon, S. Ouaskit, and T. D. Märk, Int. J. Mass Spectrom. 285, 143 (2009). 28 Z. P. Wang, P. M. Dinh, P.-G. Reinhard, E. Suraud, and F. S. Zhang, Int. J. Quantum Chem. 111, 480 (2011). 29 Z. P. Wang, J. Wang, and F. S. Zhang, Chin. Phys. Lett. 27, 053401 (2010). 30 A. Castro, M. Isla, J. I. Mart´lnez, and J. A. Alonso, Chem. Phys. 399, 130 (2012). 31 C. Z. Gao, J. Wang, and F. S. Zhang, Chem. Phys. 410, 9 (2013). 32 J. Wang, C. Z. Gao, and F. S. Zhang, Chem. Phys. Lett. 556, 256 (2013). 33 C. M. Isborn, X. Li, and J. C. Tully, J. Chem. Phys. 126, 134307 (2007). 34 Y. Miyamoto and H. Zhang, Phys. Rev. B 77, 161402 (2008). 35 S. Bubin, B. Wang, S. Pantelides, and K. Varga, Phys. Rev. B 85, 235435 (2012). 36 C. Fiolhais, F. Nogueira, and M. A. L. Marques, A Primer in Density Functional Theory (Springer, Berlin, 2003). 37 M. A. L. Marques, C. A. Ullrich, F. Nogueira, A. Rubio, K. Burke, and E. K. U. Gross, Time Dependent Density Functional Theory, Lecture Notes in Physics Vol. 706 (Springer, Berlin, 2006). 38 M. A. L. Marques, N. T. Maitra, F. Nogueira, E. K. U. Gross, and A. Rubio, Fundamentals of Time-Dependent Density Functional Theory, Lecture Notes in Physics Vol. 837 (Springer, Berlin, 2012). 39 C. A. Ullrich, Time-Dependent Density-Functional Theory: Concepts and Applications (Oxford University Press, Oxford, 2012). 40 N. T. Maitra, K. Burke, and C. Woodward, Phys. Rev. Lett. 89, 023002 (2002). 41 R. Baer, J. Mol. Struct.: THEOCHEM 914, 19 (2009). 42 E. K. U. Gross and W. Kohn, Phys. Rev. Lett. 55, 2850 (1985).

44 I.

Vignale and W. Kohn, Phys. Rev. Lett. 77, 2037 (1996). V. Tokatly, Phys. Rev. B 71, 165105 (2005). 45 C. A. Ullrich and I. V. Tokatly, Phys. Rev. B 73, 235102 (2006). 46 M. Ruggenthaler, S. E. B. Nielsen, and R. van Leeuwen, Phys. Rev. A 88, 022512 (2013). 47 J. C. Tully and R. K. Preston, J. Chem. Phys. 55, 562 (1971). 48 J. C. Tully, J. Chem. Phys. 93, 1061 (1990). 49 D. A. Micha, J. Chem. Phys. 78, 7138 (1983). 50 B. F. E. Curchod, U. Rothlisberger, and I. Tavernelli, Chem. Phys. Chem. 14, 1314 (2013). 51 J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 52 S. Kümmel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008). 53 C. Legrand, E. Surand, and P.-G. Reinhard, J. Phys. B 35, 1115 (2002). 54 F. Calvayrac, P.-G. Reinhard, E. Suraud, and C. A. Ullrich, Phys. Rep. 337, 493 (2000). 55 P.-G. Reinhard and E. Suraud, Introduction to Cluster Dynamics (Wiley, New York, 2003). 56 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 57 S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B 54, 1703 (1996). 58 M. D. Feit, J. A. Fleck, and A. Steiger, J. Comput. Phys. 47, 412 (1982). 59 A. Castro, M. A. L. Marques, and A. Rubio, J. Chem. Phys. 121, 3425 (2004). 60 L. Verlet, Phys. Rev. 159, 98 (1967). 61 J. L. Krause, K. J. Schafer, and K. C. Kulander, Phys. Rev. A 45, 4998 (1992). 62 C. A. Ullrich, J. Mol. Struct.: THEOCHEM 501, 315 (2000). 63 P. Pulay, G. Fogarasi, F. Pang, and J. E. Boggs, J. Am. Chem. Soc. 101, 2550 (1979). 64 G. Bieri and L. Asbrink, J. Electron. Spectrosc. Relat. Phenom. 20, 149 (1980). 65 D. L. Gray and A. G. Robiette, Mol. Phys. 37, 1901 (1979). 66 C. Z. Gao, J. Wang, and F. S. Zhang, Nucl. Instrum. Methods Phys. Res., Sect. B 307, 225 (2013). 67 H. Müller, W. Kutzelnigg, J. Noga, and W. Klopper, J. Chem. Phys. 106, 1863 (1997). 68 D. Marx and M. Parrinello, Nature 375, 216 (1995). 69 D. Marx and M. Parrinello, Science 284, 59 (1999). 70 E. Deumens, A. Diz, R. Longo, and Y. Öhrn, Rev. Mod. Phys. 66, 917 (1994). 71 G. Niedner, M. Noll, J. P. Toennies, and Ch. Schlier, J. Chem. Phys. 87, 2685 (1987). 72 C. F. Giese and W. R. Gentry, Phys. Rev. A 10, 2156 (1974). 73 J. A. Morales, A. Diz, E. Deumens, and Y. Öhrn, J. Chem. Phys. 103, 9968 (1995). 74 J. A. Morales, B. Maiti, Y. Yan, K. Tsereteli, J. Laraque, S. Addepalli, and C. Myers, Chem. Phys. Lett. 414, 405 (2005). 75 B. Maiti, R. Sadeghi, A. Austion, and J. A. Morales, Chem. Phys. 340, 105 (2007). 76 O. Knospe, J. Jellinek, U. Saalmann, and R. Schmidt, Phys. Rev. A 61, 022715 (2000). 77 R. Cabrera-Trujillo, I. Alvarez, C. Cisneros, A. Guerrero, and J. C. Poveda, Phys. Rev. A 82, 022718 (2010). 78 Y. Nakai and T. Shirai, At. Data Nucl. Data Tables 37, 69 (1987). 79 S. Fukata, M. Kimura, Y. Li, H.-P. Liebermann, R. J. Buenker, and L. Pichl, Plasma Fusion Res. 3, S1040 (2008). 80 B. Maiti, P. M. McLaurin, R. Sadeghi, S. A. Perera, and J. A. Morales, Int. J. Quantum Chem. 109, 3026 (2009).