Non-Adiabatic Molecular Dynamics Simulations of Non-Charge-Transfer and Charge-Transfer Scattering in H++CO2at ELab=30 eV

Yun-An Yn,Jorge A.Morles

a.Guizhou Provincial Key Laboratory of Computational Nano-material Science Guizhou Education University,Guiyang 550018,China

b.Department of Chemistry and Biochemistry,Texas Tech University,PO Box 41061,Lubbock,TX 79409-1061,USA

I.INTRODUCTION

Research on high-energy proton-molecule reactions is of primary importance in experimental and theoretical chemistry.In the present context,high energy means hyper-thermal collision energies ranging from about 10 eV up to hundreds of keV.The relevance of protonmolecule reactions arises from their preponderance in the earth’s upper atmosphere,other planets’and satellites’atmospheres,astrophysical media,plasmas,and particle accelerators.Furthermore,in recent years,this relevance has extended to the area of proton cancer therapy(PCT),a relatively new medical treatment that employs high-energy proton projectiles to obliterate cancerous cells[1−4].Aside from their relevance in the above systems,proton-molecule reactions also offer a unique opportunity to probe collision-induced rovibrational,energy-and charge-transfer processes.This results from the lack of structure of proton projectiles that simplifies the interpretation of proton-molecule spectra.Consequently,numerous high-energy proton-molecule reactions have been investigated in cross-beam scattering experiments such as H++M,with M=H2[5,6],N2[6−9],O2[6,8,10],CO[7−9],HF[7,11],HCl[7],H2O[12],CO2[6,11,13],N2O[13],CH4[14],C2H2[15],CF4[16],SF6[16],etc.Those experiments provided accurate dynamical properties of the investigated reactions and prompted a few theoretical chemistry studies(cf.Refs.[4,5,17−20]).

Out of various proton-molecule reactions investigated experimentally,H++CO2at the collision energy ELab=30 eV has attracted considerable attention due to its importance for upper atmosphere chemistry[6,11,13].Specifically,H+ions from the solar wind collide with CO2molecules in the upper atmosphere prompting various chemical reactions.Moreover,in a PCT context,the system can be included among the computationally feasible prototypes to simulate PCT reactions,such as H++(H2O)1−6to model water radiolysis[4,21−24]and H++DNA bases to model DNA damage[4,22,25−27];H++CO2can be used to model proton-induced cleavages on carboxylated biomolecules.

Among the experimental studies on H++CO2at ELab=30 eV[6,11,13],the one conducted by the Toennies group[13]has provided a detailed account of the various collision-induced processes occurring in this system:non-charge-transfer scattering,charge-transfer scattering,vibrational excitations,etc.However,unlike other similar experiments[4,5,17−22,25,28−35],no theoretical studies have been conducted on this experiment[11]to examine and further expand its findings.This may be due to the fact that it is challenging to have a single theoretical method that can describe all the simultaneous physical and chemical processes occurring in this system.Fortunately,our featured method does have such a capability and is therefore applied to H++CO2for the first time.

The traditional approach to describing ion-molecule reactions involves solving the time-independent Schrödinger equation of the corresponding scattering problem[36],a procedure requiring predetermined potential energy surfaces(PESs).Since the solution of full close-coupling scattering equations is computationally prohibitive,different decoupling approximations among degrees of freedom are usually applied.A common approach is the infinite-order sudden approximation(IOSA)[17,36−38]that prescribes the decoupling of orbital and rotational angular momenta at high-energy(“sudden”)collisions.Different types of IOSA schemes have been applied to proton-molecule reactions such as H++H2[17],N2[18],CO[19,37],and NO[20].While those investigations are highly valuable,they remain computationally onerous due to the cost of constructing complete PESs and of solving the close-coupling scattering equations(the latter task remains arduous even with the IOSA).Furthermore,the IOSA schemes are time-independent and cannot therefore provide a detailed account of a reaction at any time of its evolution from reactants to products.

To overcome the limitations of the above methods,we have successfully applied the electron nuclear dynamics(END)theory[4,39−43]to various high-energy proton-molecules reactions.END is a time-dependent,variational,on-the- fly,and non-adiabatic approach to chemical dynamics that assumes different versions according to the level of accuracy of the nuclear and electronic descriptions[4,39−43].In the simplest-level END(SLEND)version used herein,the nuclei are described classically while the electrons receive a quantum description via a Thouless single-determinantal wavefunction[44].Thus,SLEND presents a perfect balance between accuracy and computational cost that makes it competitive in comparison to traditional scattering methods.SLEND has accurately reproduced a large array of experimental properties(scattering patterns,energy transfers,non-charge-transfer,and charge-transfer cross sections,etc.)of various protonatom and-molecule reactions including H++H[45],He[45],H2[28,29,45],CH4[46],H2O[47],C2H2[30,48],HF[31],CF4[32],N2[33],CO[34],NO[35],(H2O)2−6[4,21,22,49],DNA bases[4,22,25]and other molecular systems[50].The success of SLEND in simulating those systems can be attributed to two distinctive features:(i)its on-the- fly dynamics nature that makes possible its application to diverse systems without the cumbersome predetermination of complete PESs,and(ii)its versatility that permits the description of the different types of chemical processes occurring simultaneously in the above systems(e.g.projectile scattering,charge transfers,rearrangements,dissociations,etc.).

In continuation of our SLEND investigations,we present here in the first theoretical study of the H++CO2system at ELab=30 eV experimentally investigated by the Toennies group[13].This reactive system displays a wealth of scattering processes and chemical reactions that certainly calls for a detailed theoretical investigation.Toward that end,the first part of this investigation(Section IV.A)concentrates on providing a detailed account of all the chemical processes occurring during H++CO2collisions(e.g.projectile inelastic scattering,bond dissociations(cleavages),charge-transfer reactions,etc.).Notice that this important chemical information revealed herein was not obtained by the Toennies experiment[13].Subsequently,this investigation concentrates on the scattering patterns of the outgoing projectiles(Section IV.B);this information was only partially captured by the Toennies experiment[13].The scattering patterns are relevant for the dynamics of the H+irradiation on CO2in the upper atmosphere and in PCT processes.Finally,this investigation proceeds to the calculation of differential cross sections for non-charge-transfer and charge-transfer processes and their comparison with experimental data(Section IV.C).

II.THEORY

The END theory and its SLEND version have been discussed in a series of review articles[4,39−42]and a book[43].Therefore,in the following,we present a succinct discussion of END and SLEND.END is a time dependent,variational,on-the- fly and non-adiabatic framework to simulate a large variety of chemical reactions and scattering processes.For a molecular system with NNnuclei and Neelectrons,SLEND,prescribes a total trial wavefunction Ψtotalthat is a product of nuclear ΨNand electronic Ψeparts:Ψtotal=ΨNΨe.ΨNis a product of 3NNnarrow-width,frozen Gaussian wave packets:

with average positions Ri,average momenta Pi,and widths∆Ri.In order to save computational cost,the zero-width limit∆Ri→0,∀i,is applied to ΨNafter constructing the quantum Lagrangian for the SLEND dynamical equations(cf.next paragraph).This originates a nuclear classical dynamics that still retains all the nucleus-electron non-adiabatic coupling terms.Adoption of a nuclear classical dynamics does not harm accuracy in the present type of high-energy reaction as demonstrated in previous SLEND investigations[4,39].

where the non-orthogonal Thouless dynamical spinorbitalsetc.,are singly-,doubly-excited,etc.,determinants out of|0⟩.The{ψh,ψp}MSOs are constructed with conventional atomic basis set functions centered on the nuclear wave packets;therefore,the MSOs depend upon{Ri}and so do|z;⟩.The less conventional Thouless wavefunction provides a minimum and non-redundant representation of the most general single-determinantal state|z;⟩out of and nonorthogonal to the reference[44].This reduces the number of electronic parameters{zph}and saves computational cost and also provides an invertible mapping between independent single-determinantal states and electronic parameters that avoids singularities during time evolution(cf.Refs.[4,39]for full details).Both the MSOs and DSOs are unrestricted with respect to spin so that|z;⟩is flexible to describe bond breaks as discussed in Section IV.

Having defined Ψtotal,the quantum Lagrangian L[51]for such a trial function is:

whereis the total(nuclear and electronic)Hamiltonian of the system,andandtake the time derivative to the right and to the left,respectively[51];use of those two time derivatives assures a symmetric form of the Lagrangian that is more apt for mathematical manipulations[51].After applying the zero-width limit to the wave packets in L,the time-dependent variational principle(TDVP)is applied to the quantum action S associated with L to obtain the SLEND equations.That step involves enforcing the stationarity on S(δS=δL(t)dt=0).The resulted SLEND equations for the time-dependent variational parameters{Ri,Pi,zph,}are[4,39]:

where the electron-electron C and nucleus-electron CRand CRRmetric matrix coupling terms are[4,39]:

where Etotalis the total(nuclear and electronic)energy and S=⟨z′;R′|z;is the overlap between two Thouless states.Eq.(4)expresses the coupled classical nuclear and quantum electronic dynamics in a generalized quantum Hamiltonian format[51].

In simulations,the reactants are prepared with initial positions and momentaand in the electronic stateto define the initial wavefunction:At the initial time,all is the ground UHF state corresponding to the incoming H++CO2reactants.A chemical reaction simulation is computed by integrating Eq.(4).As time elapses,the system will access the non-adiabatic regime where the Thouless parameters{zph(t)}cease to be zero and the excited-state determinants in Eq.(2)arise.At final time,all the relevant properties of a chemical reaction(e.g.electron-transfer probabilities,differential cross sections,etc.)can be calculated from the final wavefunction.For instance,in the present system,an outgoing projectile can in practice capture up to twothe probability of forming the unstable H1−nions with n≥3 is negligible.Under these conditions,the probabilities of transferring n electrons from the CO2target to the H+projectile Pn-ET,0≤n≤2,are[52]

where Nαand Nβare the Mulliken population numbers of alpha and beta-spin electrons,respectively,in the outgoing projectile.

At this point,it is opportune to explain some details of SLEND regarding electron correlation effects.SLEND utilizes an electronic unrestricted singledeterminantal wavefunction[4,39−42]that lacks dynamical correlation and has a limited amount of non-dynamical correlation;one may wonder whether this situation may affect the accuracy of the presented results.Regarding dynamical correlation,we developed the SLEND/Kohn-Sham-density-functionaltheory(SLEND/KSDFT)[4,41,42]that incorporates such effect at the level of time-dependent KSDFT[53].SLEND/KSDFT performs better for molecular geometries and vibrational frequencies than SLEND but inherits from time-dependent KSDFT a tendency to overestimate electron-transfer processes[4,41,42].For that reason,we applied SLEND to H++CO2,a system where electron-transfer processes are important.While lacking dynamical correlation,SLEND has provided good predictions of electron-transfer properties(e.g.integral and differential cross sections)in comparison with experimental data as documented in Refs.[4,39].Regarding non-dynamical correlation(correct bond breaking),a multi-configuration END theory having high levels of this correlation has been derived[40];however,that challenging method has not been coded yet.Therefore,we cannot provide numerical appraisals about non-dynamical correlation effects on SLEND results.However,from previous results[4,39],we can judge that the level of non-dynamical correlation in the unrestricted single-determinantal SLEND wavefunction provides reasonable bond-breaking descriptions,at least qualitatively.

Finally,it is opportune to elucidate a few details of the SLEND theory in comparison with other dynamical methods.Like other approaches,the SLEND total wavefunction is a product of nuclear and electronic wavefunctions,Eqs.(1)and(2).The last equation shows that during dynamics,the electronic SLEND wavefunction becomes a superposition of the ground state⟩and all the singly-,doubly-excited,etc.,states out of⟩.This superposition creates a mean- field potential under which the nuclei evolve.In this regard,SLEND is similar to other mean- field methods such as Ehrenfest dynamics[54].However,SLEND differs from those methods in a number of ways(e.g.by its use of a Thouless electronic wavefunction,Eq.(2)and by its full inclusion of non-adiabatic coupling terms,Eq.(5),absent in Ehrenfest dynamics[54]).In some methods,the nuclear wavefunction appears in the denominator of some expressions and may therefore cause singularity problems(e.g.in Bohmian mechanics[55]).The SLEND nuclear wavefunction,Eq.(1),appears in both the numerator and denominator of the quantum Lagrangian,Eq.(3).However,due to its mathematical form,the SLEND nuclear wavefunction never becomes zero,thus avoiding singularities(this can be corroborated by inspecting the values of the Gaussian wave packets in Eq.(1)for all the values of the nuclear positions and momenta).Ultimately,the zero-width limit is applied to the SLEND nuclear wavefunction in the numerator and denominator of Eq.(3)to bring about nuclear classical mechanics;that limit is finite and the resulting SLEND equations,Eq.(4),are singularity-free.More details about SLEND in comparison with other dynamical methods are discussed in Refs.[4,39]

FIG.1 H++CO2reactants initial conditions.Balls represent nuclear wave packets for the nuclei with projectile impact parameter b and target orientation[α,β].

III.COMPUTATIONAL DETAILS

The nuclear parameters defining the reactants’initial conditions for the present calculations are depicted in FIG.1.Nuclei in the CO2target are labeled:C,O(1),and O(2),respectively,whereas the nucleus making the incoming projectile H+is labeled H.The 6-31G basis set[56]is placed on the four moving nuclei to construct the electronic MSO{ψh,ψp}.Medium-size basis sets like 6-31G have provided accurate results with SLEND without posing a high computational cost[33].The CO2target is initially at rest(Pc=PO(1)=PO(2)=0),with its center of mass at(0,0,0)(cf.FIG.1)and in its electronic UHF-SCF/6-31G[56]ground state.The angular orientation[α,β]of the CO2target is determined by the polar α(0◦≤α≤180◦)and azimuthal β(0◦≤β≤360◦)angles of its inter-nuclear C−O(1)semiaxis(cf.FIG.1).The H+projectile is initially placed at(b,0,−15 a.u.),where b(0≤b≤+∞)is the impact parameter.The incoming projectile initially travels along the positive z-axis with momentumcorresponding to a kinetic energy ELab=30 eV.Bene-tions not replicating symmetry-equivalent trajectories can be generated by constraining α and β in the ranges:0◦≤α≤180◦and 0◦≤β≤90◦,respectively(e.g.the[α,β]and[α,180◦−β]orientations generate equivalent trajectories,cf.FIG.1).Two sets of trajectories starting from relevant target orientations[α,β]are calculated.Set 1 contains target initial orientations[α,0◦]with 0◦≤α≤180◦in steps of∆α=15◦and with b varying:0.0≤b≤9.0 a.u.,in steps of∆b=0.1 a.u.in the 0.0≤b≤7.0 a.u.subinterval,and in steps of∆b=0.2 a.u.in the 7.2≤b≤9.0 a.u.one.Set 2 contains collisions of H+projectile impacting orthogonally on the CO2axis and arising from targets at initial orientations[90◦,β]with β in two ranges:0◦≤β≤15◦in steps of∆β=1◦,and 20◦≤β≤90◦in steps of∆β=5◦,respectively,and with b varying as in Set 1.These initial conditions produce a total of 42 symmetry-independent target orientations and 3402 independent simulation trajectories.Set 1 and Set 2 were selected herein because they can capture some of the most reactive situations(e.g.an incoming H+projectile directly encountering one C=O bond on its way).All present simulations were run for a total time of 1,000 a.u.;that time allows the final fragments to have separations equal to or larger than those of the initial reactants.

IV.RESULTS AND DISCUSSION

A.Reactive processes and trajectory analysis

Theoretically predicted reactions during the present simulations include:

(i)Non-charge-transfer scattering(NCTS),where the H+projectile scatters off after exciting the rotational,vibrational and electronic degrees of freedom of CO2:

(ii)Charge-transfer scattering(CTS),where essentially one electron is transferred to the incoming projectile H+from the target CO2:

NCTS and CTS processes happen together on a given classical nuclear trajectory,with their respective probabilities given by P0-ETand P1-ETin Eq.(6).

(iii)Collision-induced C=O bond dissociation(C=O BD),where one C=O bond of the CO2molecule is cleaved by the incoming projectile H+:

FIG.2 C=O bond dissociations(C=O BD)from the[90◦,β]orientation in a 2-D polar representation as a function of the impact parameter b(the radius in a.u.)and the angle β(in degree).The scale of β is not linear to facilitate visualizations at low β values.There is no C=O BD outside the plotted data.

Above,the quantum number sets:(li,vi,ni)and(lf,vf,nf),label the rotational,vibrational and electronic eigenstates of the initial and final chemical species,respectively,because all channels exhibit rovibrational inelasticity and electronic excitations.Additional hypothetical reactions,such as double C=O BD:H++CO2→Hq+Cq′+Oq′′+Oq′′′(q+q′+q′′+q′′′=+1);CO/OH formation:H++CO2→COq+OHq′(q+q′=+1);and H+capture:H++CO2→HCO2+,were not observed during these simulations.The majority of the present simulations predict NCTS/CTS reactions.The simulations predicting C=O BD reactions are mostly from[90◦,β]orientations in Set 2 and from relatively small impact parameters because those initial conditions compel the incoming projectile H+to directly encounter one C=O bond.Representative[90◦,β]C=O BD reactions are listed in Table I,entered by orientations and impact parameter intervals within which the C=O BD reactions occur.The reactive channels listed in Table I are also shown in the 2-D polar plot in FIG.2.Additional C=O BD simulations from Set 1 not listed in Table I include:[75◦,0◦]/b=1.0−1.4,and[105◦,0◦]/b=0.3.The extension of the C=O BD reaction in CO2is smaller than those of analogous collision-induced reactions in previously investigated targets(e.g.H++H2[28,29],+CH4[46]and+C2H2[30]at ELab=30 eV)likely due to the higher C=O bond strength.

The two predicted reactions in the H++CO2system:C=O BD and NCTS/CTS are depicted in FIGs.3−11 for simulations from the initial conditions:[90◦,6◦]/b=0.8 a.u.(C=O BD,FIGs.3−5),[90◦,15◦]/b=2.0 a.u.(NCTS/CTS with negligible charge-transfer probability,FIGs.6−8),and[90◦,85◦]/b=1.6 a.u.(NCTS/CTS with high charge-transfer probability,FIGs.9−11),respectively.FIGs.3−11 share common features:FIGs.3,6,and 9 depict the four nuclei trajectories projected on appropriate planes(x-z plane in FIGs.3 and 6 and y-z plane in FIG.9)to obtain their clearest visualizations;FIGs.4,7,and 10 depict the relative distances between the target nuclear pairs:C−O(1),C−O(2)and O(1)−O(2)vs.time;FIGs.5,8,and 11 depict the Mulliken population analysis of the alpha-spin electrons on the nuclei:C,O(1),O(2)and H vs.time.

TABLE I C=O bond dissociation(C=O BD)reactions for representative[90◦,β◦]orientations vs.impact parameter b.

FIGs.3−5 show a C=O reaction from the[90◦,6◦]/b=0.8 a.u.initial conditions.FIG.3 portrays the incoming projectile H+(H)initially approaching the C−O(1)bond of CO2and eventually transferring a considerable amount of momentum to the O(1)atom;after that closest approach,the incoming projectile H+(H)scatters backward tracing a knot-like trajectory while O(1)apparently dissociates from CO2.The C=O BD process suspected in FIG.3 is verified in FIG.4.Therein,the three relative distances between pairs of target nuclei:C−O(1),C−O(2),and O(1)−O(2),initially remain constant at the values of the CO2equilibrium bond lengths;however,at the time of the projectiles closest approach(about 400 a.u.),both C−O(1)and O(1)−O(2)monotonically increase and remain so until the final time,whereas C−O(2)becomes oscillatory.Those simultaneous motions indicate the dissociation of CO2into O(1)and CO(2)fragments,the latter arising vibrationally excited.FIG.5 is also consistent with C=O BD and illustrates concomitant charge-transfer processes:the incoming H+(H)projectile starts with no electrons but when closest to CO2,it transitorily acquires a Mulliken population number of alpha-spin electrons from the CO2nuclei that is mostly lost after the H nucleus departs from the dissociating CO2.On the other hand,the alpha-Mulliken populations on C,O(1),and O(2)initially remain constant at their CO2equilibrium values;however,at and after the time of the projectile’s closest approach,the Mulliken populations on C and O(2)become synchronically oscillatory,whereas that on O(1)initially undergoes a brief and slight oscillation to become constant at the final time;these two phenomena are consistent with the formation of CO(2)and O(1)fragments:the former exhibiting intramolecular oscillation-driven electron transfers between C and O(2),and the latter dissociating from those two nuclei as indicated by its nonoscillatory alpha-Mulliken population that reflects no electron transfers through bonding.

FIG.3 SLEND simulation of a C=O bond dissociation from the[90◦,6◦]/b=0.8 a.u.initial conditions,nuclei trajectories on the x-z plane.

FIG.4 SLEND simulation of a C=O bond dissociation from the[90◦,6◦]/b=0.8 a.u.initial conditions,target nuclear relative distances vs.time.

FIG.5 SLEND simulation of a C=O bond dissociation from the[90◦,6◦]/b=0.8 a.u.initial conditions,alpha-Mulliken populations vs.time.

FIG.6 SLEND simulation of an almost pure non-chargetransfer scattering from the[90◦,15◦]/b=2.0 a.u.initial conditions,nuclei trajectories on the x-z plane.

FIG.7 SLEND simulation of an almost pure non-chargetransfer scattering from the[90◦,15◦]/b=2.0 a.u.initial conditions,target nuclear relative distances vs.time.

FIG.8 SLEND simulation of an almost pure non-chargetransfer scattering from the[90◦,15◦]/b=2.0 a.u.initial conditions,alpha-Mulliken populations vs.time.

FIG.9 SLEND simulation of a non-charge transfer and charge-transfer scattering with a high charge-transfer probability from the[90◦,85◦]/b=1.6 a.u.initial conditions,nuclei trajectories on the y-z plane.

FIGs.6−8 show a NCTS/CTS reaction with negligible charge-transfer probability from the[90◦,15◦]/b=2.0 a.u.initial conditions.FIG.6 portrays the incoming projectile H+(H)initially approaching the O(1)terminus of the CO2molecule and eventually transferring momentum to O(1),C,and O(2);after that closest approach,the H nucleus repulsively scatters in the right-backward direction,whereas the CO2molecule recoils,rotates,and vibrates without visible fragmentation.FIG.7 confirms a NCTS/CTS reaction without target dissociation because the relative distances:C−O(1),C−O(2),and O(1)−O(2)initially remain constant at the values of the CO2equilibrium bond lengths,but eventually all become oscillatory around the time of the projectile’s closest approach(about 400 a.u.)and remain so until the final time;those simultaneous oscillations indicate that CO2remains bonded and vibrationally excited after collision.FIG.8 is also consistent with NCTS/CTS and illustrates concomitant chargetransfer processes:the incoming H+(H)projectile starts again with no electrons but when closest to CO2,it transitorily acquires a partial Mulliken population number of alpha-spin electrons from the CO2nuclei that is mostly lost after the H nucleus departs from CO2.FIGs.9−11 show another NCTS/CTS reaction from the[90◦,85◦]/b=1.6 a.u.initial conditions but with high charge-transfer probability.The interpretation of FIGs.9−11 is similar to that of FIGs.6−8.Worth noticing,FIG.11 indicates that this NCTS/CTS reaction has a high charge-transfer probability.Therein,the incoming H+(H)projectile starts again with no electrons but when closest to CO2,it acquires a Mulliken population number of alpha-spin electrons through a process exhibiting a two-jump pattern.Unlike the previous cases,that electron transfer becomes permanent after the H nucleus and CO2separate.

B.Scattering angles,rainbow and glory scattering angles

The scattering angle θ of the outgoing projectile Hq(0≤q≤+1)with respect to the incoming z-direction is

where(i=x,y,z)are the outgoing projectile momentum components at the final time.Scatteringangle functions θ(b)vs.impact parameters b from selected[90◦,β]target orientations in Set 2 are shown in FIG.12 to illustrate the effect of the orientation angle β on θ(b).In addition,values and locations of critical angles θ(b)(rainbow and glory scattering angles[57])are listed in Table II for selected[α,β]orientations in Sets 1 and 2.

TABLE II Position and values of primary(classical and uniform Airy-treated)and secondary(classical)rainbow scattering angles,and positions of glory scattering angles from selected target orientations[α,β].

FIG.10 SLEND simulation of a non-charge transfer and charge-transfer scattering with a high charge-transfer probability from the[90◦,85◦]/b=1.6 a.u.initial conditions,target nuclear relative distances vs.time.

FIG.11 SLEND simulation of a non-charge transfer and charge-transfer scattering with a high charge-transfer probability from the[90◦,85◦]/b=1.6 a.u.initial conditions,alpha-Mulliken populations vs.time.

FIG.12 Scattering angle θ(b)vs.impact parameter[90◦,β]for selected orientations.

In FIG.12,θ(b)from the[90◦,0◦]orientation exhibits a primary rainbow angleat the high impact parameter=4.69 a.u.(cf.Table II;henceforth,the superscript CL will denote angles calculated from the classical nuclear trajectories via Eq.(10)).Rainbow angle features are related to parameters defining model projectile-target potential energy functions[57].Therefore,accurate predictions of rainbow angles features will be a measure of the correctness of the projectile-target interactions with SLEND.Away from the rainbow area,θ(b)becomes a decreasing function at impact parameters b somewhat lower thanand reaches a zero value θG(bG)=0◦known as the glory scattering angle[57]at the intermediate impact parameter bG(cf.Table II).Other zero and nearly zero scattering angles θS→0◦,known as small angles[57],also occur from all orientations at very large b(cf.FIG.12)due to weak projectiletarget interactions at long separations.Finally,where b

ondary rainbow angle approximately corresponds to a maximum off-plane scattering[57],a phenomenon also inferred experimentally and predicted with SLEND and other theories in other proton-molecule reactions(e.g.H++H2[5,28,29]and H++C2H2[15,30]).

C.Differential cross sections

The most important result from the Toennies group experiment on H++CO2at ELab=30 eV is the total NCTS and CTS differential cross sections(DCS)[13].In the present context,quantum inelastic DCS[dσi→f(θ)/dΩ]from initial(i)to final(f)states(i→f)can be calculated as[57]:

where fi→f(θ)is the scattering amplitude,ki(kf)is the projectile initial( final)wave vector,l is orbital angular momentum quantum numbers,(l)is the T-matrix,and Pl(cosθ)is Legrende polynomials.Since SLEND employs classical nuclear dynamics,the quantum DCS in Eq.(11)should be evaluated using semi-classical techniques[4,57,58]as follows.When θ(b)is not close to rainbow angles,the stationary phase approximation[57]gives the semi-classical DCS

wherebels the N classical projectile trajectories(branches)that from impact parameters bmlead to the same scattering direction θ;in the present simulations,the maximum values of N are:3 from[α,0◦]and[90◦,0◦<β≤50◦],1 from[90◦,55◦<β≤85◦],and 5 from[90◦,90◦],respectively(cf.FIG.12).In addition,

pattern;that classical DCS is the lowest-level approximation to the experimental DCS.Both Eqs.(12)and(14)exhibit unphysical singularities at the primary and secondary rainbow angle directions because their common denominator terms,[dθ(b)/db]b=bm,become zero in those directions.Those singularities result from the breakdown of the stationary phase approximation[57],implicit in Eqs.(12)and(14),when dθ(b)/db=0◦;more sophisticated semi-classical techniques are required in those situations.The uniform Airy[57,59]approximation is the most rigorous semi-classical treatment of an isolated rainbow angle that occurs in all[α,0◦]simulations.There,the scattering amplitudes f2(θ)and f3(θ)from two attractive trajectories near the primary rainbow angle direction(Eq.(12))are replaced with a combined single amplitude(θ):the Airy function and its first derivative[60].The uniform-Airy DCS from Eq.(15)evolves uniformly to the DCS in Eq.(12)when θ is away from the rainbow angle and to the transitional Airy approximation[59]DCS when θ is close to the rainbow angle.The uniform Pearcey[59]approximation is the most rigorous semi-classical treatment of the two adjacent rainbow angles that occur in all[90◦,0◦<β≤50◦]simulations but it is too cumbersome to apply.However,that approximation approaches the uniform-Airy expression near a rainbow angle that is somewhat distant from the other[59].Therefore,a simpler but still accurate treatment successfully used in previous SLEND studies[30−32]is adopted herein.That approach involves approximating the Pearcey scat-near the rainbow anglesrespectively,and switching between those amplitudes at the inflexion point bIPofThe above uniform approximations are valid for rainbow angles not close to 0◦;therefore,treatments of a few low-lying rainbow angle singularities are not possible with those techniques and are left uncorrected(e.g.the two rainbow angles:simulation).Contrastingly,the[90◦,55◦<β≤85◦]simulations do not exhibit rainbow angles and their DCS can be directly calculated with Eq.(12)with only one branch(N=1).

In order to illustrate the above procedure to evaluate DCS,calculated total NCTS DCS for H++CO2at ELab=30 eV from the individual target orientation[90◦,0◦]are shown in FIG.13,both at the classical(Eq.(14))and the uniform-Airy(Eqs.(12)and(15))levels,and along with the experimental result[13].The latter was reported un-normalized but is normalized in FIG.13 to allow comparisons.The normalization criterion involves attaining a minimum deviation between experimental and uniform-Airy results over all the investigated scattering angles θ.The classical DCS shows an unphysical rainbow angle singularly atwhile the uniform-Airy CS DCS shows a bounded rainbow angle maximum near θ=10◦in agreement with the experimental result.Except for that singularity,the classical DCS agrees fairly with the measured data.The uniform-Airy DCS agrees better and even displays an oscillatory pattern produced by the coherent sums of the action terms inEq.(13),and the Airy terms,Eq.(15),[57].Those genuine quantum oscillations are artificially ab-

FIG.13 Individual non-charge-transfer differential cross section for the H++CO2system at ELab=30 eV from the[90◦,0◦]orientation;quasi-classical,uniform Airy-treated,and experimental results.

FIG.14 Differential cross sections for the charge-transfer and non-charge-transfer channels in H++CO2system at ELab=30 eV;calculated and experimental results.

where A23(θ)and ξ23(θ)are functions of the action S along those trajectories[59],and Ai and Ai′are sent in the classical DCS(due to incoherent sum and ages of the DCSs from the simulated target orientations;they are shown in FIG.14 at the uniform-Airy level along with their corresponding experimental results[13].The un-normalized experimental NCTS and CTS DCS are normalized in FIG.13 as has been done in FIG.13 to allow comparisons.The theoretical DCS in FIG.14 are averages of individual DCSs from the simulated target orientations.Those averages considerably level out the oscillatory patterns displayed by the individual DCS(cf.FIG.13)and produce results closely resembling the slightly oscillatory experimental data.The calculated NCTS DCS in FIG.14 agrees well with its measured counterpart over all the investigated scattering angles θ.Conversely,the calculated CTS DCS in FIG.14 agrees well with its measured counterpart at large scattering angles θ≥8◦but it is lower in value than that counterpart at small scattering angles θ<8◦.Notably,both calculated NCTS and CTS DCS exhibit bounded rainbow angle features at θ=10◦in very good agreement with those in their experimental counterparts[13].As stated previously,rainbow angle features are related to parameters defining model projectile-target potential energy functions[57].Therefore,the current accurate prediction of rainbow angles features is a measure of the correctness of the projectile-target interactions with SLEND.The positions of the primary rainbow angle features in the callack of Airy terms,cf.Eq.(14))and are scarcely perceptible in the experimental one.The experiment collects DCS contributions from all target orientations;that“experimental”average considerably levels out the oscillatory pattern in the final DCS as well as the calculated average over uniform-Airy DCS as shown below.The comparison between classical(Eq.(14))and uniform-Airy(Eqs.(12)and(15))DCS in FIG.14 demonstrates the better performance of the latter;therefore,only uniform-Airy DCS free of singularities will be used in the final DCS calculations shown below.With the uniform Airy approximation,a classical rainbow singularity is substituted by the first bounded maximum of the Airy function.For a primary rainbow angle,that bounded peak always occurs at a scatter-culated DCS are averages of thevalues per target orientation weighted with the individual DCS intensities;those positions distribute closely around θ=10◦(cf.Table II)and thus produce clear primary rainbow features in the calculated DCS.The experimental NCTS DCS in FIG.14 seems to display a faint secondary rainbow feature at around θ=7◦[13]that is not clearly reproduced by the calculated NCTS DCS.Individual DCS per target orientation do exhibit secondary rainbow angle features but they get blurred on average because their positions distribute broadly over target orientations(cf.Table II);that theoretical finding also explains the faintness of the secondary rainbow features in the experimental results.Further simulations from additional orientations may eventually sharpen the secondary rainbow angles’distribution and enhance the secondary rainbow features in the calculated DCS.Additional aspects of the calculated total DCS are discussed in Section V.values provide a more accurate prediction of the experimental rainbow angle

V.CONCLUSION

The SLEND method[4,39]has been applied to the H++CO2reactive system at ELab=30 eV[13]for the first time to assess and further elucidate previous experimental results from this system[13].The studied system is relevant to atmospheric chemistry,astrophysics and proton cancer therapy(PCT)research.SLEND[4,39]describes nuclei in terms of classical mechanics and electrons in terms of a single-determinantal Thouless wavefunction[44].The present investigation from selected orientations are listed in Table II.

The final predicted total DCS for both NCTS and CTS processes in H++CO2at ELab=30 eV are aver-comprises a total of 3402 reaction simulations from 42 independent CO2orientations.The first part of this investigation(Section IV.A)focused on predicting the chemical reactions occurring in the H++CO2system;this important chemical information could not be obtained by the experiment[13]and is therefore revealed by the present simulations.SLEND calculations only predicted three types of reactive processes:non-charge-transfer scattering(NCTS),charge-transfer scattering(CTS),and C=O bond breaking(C=O BD),cf.Eqs.(7)−(9).The relevant features of these reactions’mechanisms were displayed and discussed in full detail in terms of the nuclei motions and charge-transfer processes(cf.FIGs.3−11).H++CO2at ELab=30 eV turned out to be less reactive in terms of both extension and variety of reactions than similar proton-molecule systems simulated with SLEND(e.g.H++H2[28,29],+CH4[46]and+C2H2[30]at ELab=30 eV);this may be due to the higher strength of the C=O bonds.In the second part of this investigation(Section IV.B),a detailed account of the scattering patterns of the outgoing H+projectile was provided.This scattering information could only be inferred in a partial form from the experiment[13].This information is very important to understand the scattering of H+projectiles irradiated on a CO2mass.The calculated scattering functions revealed complex scattering patterns in terms of rainbow and glory angles.Simulations from the[α,0◦]orientations(i.e.those with the initial H+and CO2on the same plane)and the[90◦,90◦]one(i.e.that with the initial projectile hitting orthogonally to a plane containing CO2)keep the H+projectile on the same plane during scattering due to the symmetry of the projectile-target interaction forces.As a result,scattering functions from those orientations display a primary rainbow angle and a glory scattering angle.In addition,simulations from orientations[90◦,0◦<β<90◦](i.e.those where the above planar symmetries between H+and CO2at initial time are missing)do not keep the H+projectile on the same plane during scattering.As a result,simulations from the[90◦,0◦<β≤50◦]orientation produce scattering functions displaying a primary and a secondary[61]rainbow angle and no glory angle.This phenomenon involving secondary rainbows was also inferred experimentally in H+and CO2and predicted with SLEND and other theories in other proton-molecule reactions(e.g.H++H2[5,29]and H++C2H2[15,30]).Moreover,in simulations from the[90◦,55◦≤β<90◦]orientations,the primary and secondary rainbows coalesce into an inflexion point so that all rainbow angle features disappear in the scattering functions from those orientations.Finally,in the last part of this investigation(Section IV.C),NCTS and CTS differential cross sections(DCSs)were calculated and compared with their experimental counterparts[13].The reported DCSs are obtained with sophisticated semiclassical techniques and compared with less accurate classical DCSs.SLEND NCTS DCS agrees very well with its experimental counterpart[13]over all the experimentally measured scattering angles 0◦<θ≤18◦.On the other hand,SLEND CTS DCS agrees well with its experimental counterpart[13]at large scattering angles θ≥8◦but is lower in value than its experimental counterpart at small scattering angles θ<8◦.Given all the available information,it is not possible to determine at this moment whether the experiment overestimates or SLEND underestimates the CTS DCS at low scattering angles.Additional experiments and additional calculations with higher END versions having higher accuracy will shed light on this issue.Remarkably,both the NCTS and CTS SLEND DCSs display a primary rainbow angle peak at the scattering angle θ=10◦as do their experimental counterparts[13].Rainbow angle features are related to the projectile-target interactions;therefore,this good rainbow angle prediction indicates that SLEND reproduce proper projectile-target forces.Like the experimental data[13],SLEND NCTS and CTS DCSs do not clearly show a secondary rainbow angle peak.However,the SLEND calculations demonstrate that such a feature cannot manifest clearly in the DCSs due to the wide distribution of the secondary rainbow angle values(cf.Table II).

VI.ACKNOWLEDGMENTS

Present calculations were performed at the Texas Tech University High Performance Computer Center and the Texas Advanced Computing Center at the University of Texas at Austin.Prof.Morales acknowledges financial support from the Cancer Prevention and Research Institute of Texas(CPRIT)grant RP140478.Prof.Yan acknowledges the financial support from the National Natural Science Foundation of China(No.21373064)and the Program for Innovative Research Team of Guizhou Province(No.QKTD[2014]4021).

[1]A.V.Solov’yov,E.Surdutovich,E.Scifoni,I.Mishustin,and W.Greiner,Phys.Rev.E 79,011909(2009).

[2]E.Surdutovich and A.V.Solov’yov,Eur.Phys.J.D 68,353(2014).

[3]W.P.Levin,H.Kooy,J.S.Loeffler,and T.F.DeLaney,Br.J.Cancer 93,849(2005).

[4]C.Stopera,T.V.Grimes,P.M.McLaurin,A.Privett,and J.A.Morales,Adv.Quantum Chem.66,113(2013).

[5]G.Niedner,M.Noll,J.P.Toennies,and C.Schlier,J.Chem.Phys.87,2685(1987).

[6]U.Hege and F.Linder,Z.Phys.A 320,95(1985).

[7]H.Udseth,C.F.Giese,and W.R.Gentry,J.Chem.Phys.60,3051(1974).

[8]F.A.Gianturco,U.Gierz,and J.P.Toennies,J.Phys.B 14,667(1981).

[9]J.Krutein and F.Linder,J.Chem.Phys.71,599(1979).

[10]M.Noll and J.P.Toennies,J.Chem.Phys.85,3313(1986).

[11]B.Friedrich,F.A.Gianturco,G.Niedner,M.Noll,and J.P.Toennies,J.Phys.B 20,3725(1987).

[12]B.Friedrich,G.Niedner,M.Noll,and J.P.Toennies,J.Chem.Phys.87,5256(1987).

[13]G.Niedner,M.Noll,and J.P.Toennies,J.Chem.Phys.87,2067(1987).

[14]Y.N.Chiu,B.Friedrich,W.Maring,G.Niedner,M.Noll,and J.P.Toennies,J.Chem.Phys.88,6814(1988).

[15]N.Aristov,G.Niedner-Schatteburg,J.P.Toennies,and Y.N.Chiu,J.Chem.Phys.95,7969(1991).

[16]U.Gierz,M.Noll,and J.P.Toennies,J.Chem.Phys.83,2259(1985).

[17]M.Baer,G.Niedner-Schatteburg,and J.P.Toennies,J.Chem.Phys.91,4169(1989).

[18]F.A.Gianturco,S.Kumar,T.Ritschel,R.Vetter,and L.Zlicke,J.Chem.Phys.107,6634(1997).

[19]T.J.D.Kumar and S.Kumar,J.Chem.Phys.121,191(2004).

[20]S.Amaran and S.Kumar,J.Chem.Phys.128,124306(2008).

[21]A.J.Privett,E.S.Teixeira,C.Stopera,and J.A.Morales,PLoS One 12,e0174456(2017).

[22]P.M.McLaurin,A.J.Privett,C.Stopera,T.V.Grimes,A.Perera,and J.A.Morales,Mol.Phys.113,297(2015).

[23]M.Murakami,T.Kirchner,M.Horbatsch,and H.J.Ldde,Phys.Rev.A 85,052704(2012).

[24]C.Champion,M.A.Quinto,J.M.Monti,M.E.Galassi,P.F.Weck,O.A.Fojón,J.Hanssen,and R.D.Rivarola,Phys.Med.Biol.60,7805(2015).

[25]A.J.Privett and J.A.Morales,Chem.Phys.Lett.603,82(2014).

[26]C.Champion,P.F.Weck,H.Lekadir,M.E.Galassi,O.A.Fojn,P.Abufager,R.D.Rivarola,and J.Hanssen,Phys.Med.Biol.57,3039(2012).

[27]H.Lekadir,I.Abbas,C.Champion,and J.Hanssen,Nucl.Instrum.Methods Phys.Res.Sect.B 267,1011(2009).

[28]J.A.Morales,A.C.Diz,E.Deumens,and Y.Öhrn,Chem.Phys.Lett.233,392(1995).

[29]J.Morales,A.Diz,E.Deumens,and Y.Öhrn,J.Chem.Phys.103,9968(1995).

[30]J.A.Morales,B.Maiti,Y.A.Yan,K.Tsereteli,J.Laraque,S.Addepalli,and C.Myers,Chem.Phys.Lett.414,405(2005).

[31]B.Maiti,R.Sadeghi,A.Austin,and J.A.Morales,Chem.Phys.340,105(2007).

[32]B.Maiti,P.M.McLaurin,R.Sadeghi,S.Ajith Perera,and J.A.Morales,Int.J.Quantum Chem.109,3026(2009).

[33]C.Stopera,B.Maiti,T.V.Grimes,P.M.McLaurin,and J.A.Morales,J.Chem.Phys.134,224308(2011).

[34]C.Stopera,B.Maiti,T.V.Grimes,P.M.McLaurin,and J.A.Morales,J.Chem.Phys.136,054304(2012).

[35]C.Stopera,B.Maiti,and J.A.Morales,Chem.Phys.Lett.551,42(2012).

[36]T.P.Tsien,G.A.Parker,and R.T.Pack,J.Chem.Phys.59,5373(1973).

[37]F.A.Gianturco,U.T.Lamanna,and M.Attimonelli,Chem.Phys.48,399(1980).

[38]S.V.K.Kumar,V.S.Ashoka,and E.Krishnakumar,Phys.Rev.A 70,052715(2004).

[39]E.Deumens,A.Diz,R.Longo,and Y.Öhrn,Rev.Mod.Phys.66,917(1994).

[40]E.Deumens and Y.Öhrn,J.Phys.Chem.A 105,2660(2001).

[41]K.Tsereteli,Y.A.Yan,and J.A.Morales,Chem.Phys.Lett.420,54(2006).

[42]S.A.Perera,P.M.McLaurin,T.V.Grimes,and J.A.Morales,Chem.Phys.Lett.496,188(2010).

[43]F.Hagelberg,Electron Dynamics in Molecular Interactions:Principles and Applications,Singapore:World Scientific,(2014).

[44]D.J.Thouless,Nucl.Phys.21,225(1960).

[45]R.Longo,E.Deumens,and Y.Öhrn,J.Chem.Phys.99,4554(1993).

[46]D.Jacquemin,J.A.Morales,E.Deumens,and Y.Öhrn,J.Chem.Phys.107,6146(1997).

[47]M.Hedström,J.A.Morales,E.Deumens,and Y.Öhrn,Chem.Phys.Lett.279,241(1997).

[48]S.A.Malinovskaya,R.Cabrera-Trujillo,J.R.Sabin,E.Deumens,and Y.Öhrn,J.Chem.Phys.117,1103(2002).

[49]O.Quinet,E.Deumens,and Y.Öhrn,Int.J.Quantum Chem.109,259(2009).

[50]F.Hagelberg and E.Deumens,Phys.Rev.A 65,052505(2002).

[51]P.Kramer and M.Saraceno,Geometry of the Time-Dependent Variational Principle in Quantum Mechanics,Berlin:Springler-Verlag,(1981).

[52]P.M.McLaurin,PhD Dissertation,Texas:Texas Tech University,(2011).

[53]R.G.Parr and W.Yang,Density-Functional Theory of Atoms and Molecules,New York:Oxford University Press,(1989).

[54]X.S.Li,J.C.Tully,H.B.Schlegel,and M.J.Frisch,J.Chem.Phys.123,084106(2005).

[55]R.E.Wyatt,Quantum Dynamics with Trajectories,New York:Springer Velag,(2005).

[56]A.Szabo and N.S.Ostlund,Modern Quantum Chemistry:Introduction to Advanced Electronic Structure Theory,New York:Dover Publications,Inc.,(1989).

[57]M.S.Child,Molecular Collision Theory,Mineola:Dover Publications,Inc.,(1984).

[58]J.A.Morales,Mol.Phys.108,3199(2010).

[59]J.N.L.Connor and D.Farrelly,J.Chem.Phys.75,2831(1981).

[60]M.Abramowitz and I.A.Stegun,Handbook of Mathematical Functions:with Formulas,Graphs,and Mathematical Tables,New York:Dover Publications,Inc.,(1972).

[61]C.F.Giese and W.R.Gentry,Phys.Rev.A 10,2156(1974).