NonAdiabatic Molecular Dynamics Simulations of NonChargeTransfer and ChargeTransfer Scattering in H^{+}+CO_{2} at E_{Lab}=30 eV
Ⅰ. INTRODUCTION
Research on highenergy protonmolecule reactions is of primary importance in experimental and theoretical chemistry. In the present context, high energy means hyperthermal 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 highenergy proton projectiles to obliterate cancerous cells [14]. Aside from their relevance in the above systems, protonmolecule reactions also offer a unique opportunity to probe collisioninduced rovibrational, energy and chargetransfer processes. This results from the lack of structure of proton projectiles that simplifies the interpretation of protonmolecule spectra. Consequently, numerous highenergy protonmolecule reactions have been investigated in crossbeam scattering experiments such as H$^+$+M, with M=H$_2$ [5, 6], N$_2$ [69], O$_2$ [6, 8, 10], CO [79], HF [7, 11], HCl [7], H$_2$O [12], CO$_2$ [6, 11, 13], N$_2$O [13], CH$_4$ [14], C$_2$H$_2$ [15], CF$_4$ [16], SF$_6$ [16], etc. Those experiments provided accurate dynamical properties of the investigated reactions and prompted a few theoretical chemistry studies (cf. Refs.[4, 5, 1720]).
Out of various protonmolecule reactions investigated experimentally, H$^+$+CO$_2$ at the collision energy $E_{\textrm{Lab}}$=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 CO$_2$ molecules 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$^+$+(H$_2$O)$_{16}$ to model water radiolysis [4, 2124] and H$^+$+DNA bases to model DNA damage [4, 22, 2527]; H$^+$+CO$_2$ can be used to model protoninduced cleavages on carboxylated biomolecules.
Among the experimental studies on H$^+$+CO$_2$ at $E_{\textrm{Lab}}$=30 eV [6, 11, 13], the one conducted by the Toennies group [13] has provided a detailed account of the various collisioninduced processes occurring in this system: nonchargetransfer scattering, chargetransfer scattering, vibrational excitations, etc. However, unlike other similar experiments [4, 5, 1722, 25, 2835], 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$^+$+CO$_2$ for the first time.
The traditional approach to describing ionmolecule reactions involves solving the timeindependent Schrödinger equation of the corresponding scattering problem [36], a procedure requiring predetermined potential energy surfaces (PESs). Since the solution of full closecoupling scattering equations is computationally prohibitive, different decoupling approximations among degrees of freedom are usually applied. A common approach is the infiniteorder sudden approximation (IOSA) [17, 3638] that prescribes the decoupling of orbital and rotational angular momenta at highenergy ("sudden") collisions. Different types of IOSA schemes have been applied to protonmolecule reactions such as H$^+$+H$_2$ [17], N$_2$ [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 closecoupling scattering equations (the latter task remains arduous even with the IOSA). Furthermore, the IOSA schemes are timeindependent 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, 3943] to various highenergy protonmolecules reactions. END is a timedependent, variational, onthefly, and nonadiabatic approach to chemical dynamics that assumes different versions according to the level of accuracy of the nuclear and electronic descriptions [4, 3943]. In the simplestlevel END (SLEND) version used herein, the nuclei are described classically while the electrons receive a quantum description via a Thouless singledeterminantal 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, nonchargetransfer, and chargetransfer cross sections, etc.) of various protonatom and molecule reactions including H$^+$+H [45], He [45], H$_2$ [28, 29, 45], CH$_4$ [46], H$_2$O [47], C$_2$H$_2$ [30, 48], HF [31], CF$_4$ [32], N$_2$ [33], CO [34], NO [35], (H$_2$O)$_{26}$ [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: (ⅰ) its onthefly dynamics nature that makes possible its application to diverse systems without the cumbersome predetermination of complete PESs, and (ⅱ) 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 herein the first theoretical study of the H$^+$+CO$_2$ system at $E_{\textrm{Lab}}$=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 Ⅳ. A) concentrates on providing a detailed account of all the chemical processes occurring during H$^+$+CO$_2$ collisions (e.g. projectile inelastic scattering, bond dissociations (cleavages), chargetransfer 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 Ⅳ. 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 CO$_2$ in the upper atmosphere and in PCT processes. Finally, this investigation proceeds to the calculation of differential cross sections for nonchargetransfer and chargetransfer processes and their comparison with experimental data (Section Ⅳ. C).
Ⅱ. THEORY
The END theory and its SLEND version have been discussed in a series of review articles [4, 3942] and a book [43]. Therefore, in the following, we present a succinct discussion of END and SLEND. END is a timedependent, variational, onthefly and nonadiabatic framework to simulate a large variety of chemical reactions and scattering processes. For a molecular system with $N_\textrm{N}$ nuclei and $N_\textrm{e}$ electrons, SLEND, prescribes a total trial wavefunction $\Psi_{\textrm{total}}$ that is a product of nuclear $\Psi_\textrm{N}$ and electronic $\Psi_\textrm{e}$ parts: $\Psi_{\textrm{total}}$=$\Psi_\textrm{N}$$\Psi_\textrm{e}$. $\Psi_\textrm{N}$ is a product of 3$N_\textrm{N}$ narrowwidth, frozen Gaussian wave packets:
$
\begin{align}
& {{\Psi }_{\text{N}}}(\mathbf{X};\mathbf{R},\mathbf{P})= \\
& \underset{i=1}{\overset{3{{N}_{\text{N}}}}{\mathop{\prod }}}\,\exp \left[ {{\left( \frac{{{X}_{i}}{{R}_{i}}}{2\Delta {{R}_{i}}} \right)}^{2}}+i{{P}_{i}}({{X}_{i}}{{R}_{i}}) \right] \\
\end{align}
$

(1) 
with average positions $R_i$, average momenta $P_i$, and widths $\Delta R_i$. In order to save computational cost, the zerowidth limit $\Delta R_i$$\rightarrow$0, $\forall i$, is applied to $\Psi_\textrm{N}$ after constructing the quantum Lagrangian for the SLEND dynamical equations (cf. next paragraph). This originates a nuclear classical dynamics that still retains all the nucleuselectron nonadiabatic coupling terms. Adoption of a nuclear classical dynamics does not harm accuracy in the present type of highenergy reaction as demonstrated in previous SLEND investigations [4, 39]. For $\Psi_\textrm{e}$, $K$ fermion creation (annihilation) operators $b_i^\dagger $ ($b_i$) of $N_\textrm{e}$ occupied $\{\Psi_{{h}}\}$ and $K$$N_\textrm{e}$ virtual $\{\Psi_\textrm{p}\}$ orthogonal molecular spinorbitals (MSO) are considered to generate particlehole pair operators $b_{p}^\dagger b_{h}$. Taking as a reference state the unrestricted HartreeFock (UHF) state $$0$\rangle$=$b_{N_\textrm{e}}^+$$\cdots$$b_{l}^+$$$vac$\rangle$, $\Psi_\textrm{e}$ is a singledeterminantal wavefunction in the Thouless representation $\bf{z};\bf{R}\rangle$ [44]
$
\begin{eqnarray}
\Psi _\textrm{e } \hspace{0.15cm}&=&\hspace{0.15cm} \left {{\bf{z}};{\rm{ }}{\bf{R}}} \right\rangle \\
&=& \hspace{0.15cm}\det \left[ {\chi _{h} \left( {{\bf{z}};{\rm{ }}{\bf{R}}} \right)} \right] \\
&=&\hspace{0.15cm} \exp \left( {\sum\limits_{p = N_\textrm{e} + 1}^K {\sum\limits_{h = 1}^{N_\textrm{e} } {z_{ph} b_p^\dagger b_h } } } \right)\left 0 \right\rangle \\
{\rm{ }} &=&\hspace{0.15cm} \left {\rm{0}} \right\rangle + \sum\limits_{p, h} {z_{{ph}} } \left( t \right)b_{p}^\dagger b_{h} \left {\rm{0}} \right\rangle + \\
&\hspace{0.15cm}\sum\limits_{p, h, p', h'} {z_{{ph}} } \left( t \right)z_{{p'h}'} \left( t \right)b_p^\dagger b_{p'}^\dagger b_{h'} b_h \left {\rm{0}} \right\rangle {\rm{ }} + {\rm{ }}...\\
\chi _h\hspace{0.15cm}&=&\hspace{0.15cm}\psi _h+\sum\nolimits_{p = N_\textrm{e} + 1}^K {{\rm{ }}\psi _p } z_{{ph}} (1\leq h\leq\left. {N_\textrm{e} } \right)
\end{eqnarray}
$

(2) 
where the nonorthogonal Thouless dynamical spinorbitals (DSO) $\left\{ {\chi _h } \right\}$ are $\chi _h$ [44] and $
b_p^\dagger b_h$$$0$\rangle$, $b_p^\dagger b_{p'}^\dagger
b_{h'} b_h$$$0$\rangle$, etc., are singly, doublyexcited, etc., determinants out of $$0$\rangle$. The $\left\{ {\psi _h , } {\psi _p } \right\}$ MSOs are constructed with conventional atomic basis set functions centered on the nuclear wave packets; therefore, the MSOs depend upon $\left\{ {R_i } \right\}$ and so do $\left {{\bf{z}};{\rm{ }}{\bf{R}}} \right\rangle$. The less conventional Thouless wavefunction provides a minimum and nonredundant representation of the most general singledeterminantal state $\left {{\bf{z}};{\rm{ }}{\bf{R}}} \right\rangle$ out of and nonorthogonal to the reference $$0$\rangle$ [44]. This reduces the number of electronic parameters $\left\{ {z_{{ph}} } \right\}$ and saves computational cost and also provides an invertible mapping between independent singledeterminantal 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 $\left {{\bf{z}};{\rm{ }}{\bf{R}}} \right\rangle$ is flexible to describe bond breaks as discussed in Section Ⅳ.
Having defined $\Psi _{\textrm{total}}$, the quantum Lagrangian $L$ [51] for such a trial function is:
$
L = \frac{{\left\langle {{\Psi _{{\rm{total}}}}\left[ {\frac{{i\hbar }}{2}\left( {\frac{{\vec \partial }}{{\partial t}}  \frac{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftarrow$}}
\over \partial } }}{{\partial t}}} \right)  \hat H} \right]{\Psi _{{\rm{total}}}}} \right\rangle }}{{\left\langle {{\Psi _{{\rm{total}}}}\left {{\Psi _{{\rm{total}}}}} \right.} \right\rangle }}
$

(3) 
where $\hat H$ is the total (nuclear and electronic) Hamiltonian of the system, and ${\vec \partial }/{\partial t}$ and ${\vec \partial }/{\partial t}$ take 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 zerowidth limit to the wave packets in $L$, the timedependent 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$ ($\delta S$=$\delta$$
\int_{t_1 }^{t_2 }$$L$($t$)d$t$=0). The resulted SLEND equations for the timedependent variational parameters $\{R_i, P_i, z_{ph}, z_{ph}^*\}$ are [4, 39]:
$
\left[ \begin{matrix}
i\mathbf{C} & 0 & i{{\mathbf{C}}_{\mathbf{R}}} & 0 \\
0 & i{{\mathbf{C}}^{*}} & i\mathbf{C}_{\mathbf{R}}^{*} & 0 \\
i\mathbf{C}_{\mathbf{R}}^{\dagger } & i\mathbf{C}_{\mathbf{R}}^{T} & {{\mathbf{C}}_{\mathbf{RR}}} & \mathbf{I} \\
0 & 0 & \mathbf{I} & 0 \\
\end{matrix} \right]\left[ \begin{matrix}
\frac{\text{d}\mathbf{z}}{\text{d}t} \\
\frac{\text{d}{{\mathbf{z}}^{*}}}{\text{d}t} \\
\frac{\text{d}\mathbf{R}}{\text{d}t} \\
\frac{\text{d}\mathbf{P}}{\text{d}t} \\
{} \\
\end{matrix} \right]=\left[ \begin{matrix}
\frac{\partial {{E}_{\text{total}}}}{\partial {{\mathbf{z}}^{*}}} \\
\frac{\partial {{E}_{\text{total}}}}{\partial \mathbf{z}} \\
\frac{\partial {{E}_{\text{total}}}}{\partial \mathbf{R}} \\
\frac{\partial {{E}_{\text{total}}}}{\partial \mathbf{P}} \\
{} \\
\end{matrix} \right]
$

(4) 
where the electronelectron ${\bf{C}}$ and nucleuselectron ${\bf{C}}_{\bf{R}}$ and ${\bf{C}}_{{\bf{RR}}}$ metric matrix coupling terms are [4, 39]:
$
\begin{eqnarray}
\begin{array}{*{20}c}
\displaystyle{{\bf{C}}_{ph:qg} = \left. {\frac{{\partial ^2 \ln S}}{{\partial z_{ph}^* \partial z_{qg} }}} \right_{{\bf{R' = R}}} }\\[4ex]
{\rm{ }}\left( {{\bf{C}}_{{\bf{X}}_{{\bf{ik}}} } } \right)_{ph} = \displaystyle{\left. {\frac{{\partial ^2 \ln S}}{{\partial z_{ph}^* \partial X_{ik} }}} \right_{{\bf{R' = R}}}} \\[4ex]
{\rm{ }}\left( {{\bf{C}}_{{\bf{XY}}} } \right)_{ij:kl} =\displaystyle{\left. {  2{\mathop{\rm Im}\nolimits} \frac{{\partial ^2 \ln S}}{{\partial X_{ik} \partial Y_{jl} }}} \right_{{\bf{R' = R}}}}
\end{array}
\end{eqnarray}
$

(5) 
where $E_{{\rm{total}}} $ is the total (nuclear and electronic) energy and $S$=$\left\langle {{\bf{z'}};{\rm{ }}{\bf{R'}}} \right{{\bf{z}};{\rm{ }}{\bf{R}}}\rangle$ 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 momenta $\left\{ {{\rm{ }}R_i^0 } \right\}$ and $\left\{ {{\rm{ }}P_i^0 } \right\}$ (cf. FIG. 1) and in the electronic state $\left\{ {z_{ph}^0 } \right. \left. {z_{ph}^{0^*} } \right\}$, to define the initial wavefunction: $
\Psi _{\textrm{total}}^i$. At the initial time, all $z_{ph}^0$ are zero ($z_{ph}^0$=0), so that $\left {{\bf{z}};{\rm{ }}{\bf{R}}} \right\rangle$=$$0$\rangle$, i.e. $\left {{\bf{z}};{\rm{ }}{\bf{R}}} \right\rangle$ is the ground UHF state corresponding to the incoming H$^+$+CO$_2$ reactants. A chemical reaction simulation is computed by integrating Eq.(4). As time elapses, the system will access the nonadiabatic regime where the Thouless parameters $\{ {z_{ph}( t )}\}$ cease to be zero and the excitedstate determinants in Eq.(2) arise. At final time, all the relevant properties of a chemical reaction (e.g. electrontransfer probabilities, differential cross sections, etc.) can be calculated from the final wavefunction $\Psi _{\textrm{total}}^f$. For instance, in the present system, an outgoing projectile can in practice capture up to two electrons: H$^+$+CO$_2$$\rightarrow$H$^{1n}$+CO$_2{^{n}}$, 0$\leq$$n$$\leq$2, because the probability of forming the unstable H$^{1n}$ ions with $n$$\geq$3 is negligible. Under these conditions, the probabilities of transferring $n$ electrons from the CO$_2$ target to the H$^+$ projectile $P_{n {\rm{\textrm{}ET}}}$, 0$\leq$$n$$\leq$2, are [52]
$
\begin{eqnarray}
P_{0 {\rm{\textrm{}ET}}} \hspace{0.15cm} &=&\hspace{0.15cm} \left( {1  N_\alpha } \right)\left( {1  N_\beta } \right) \\
P_{1 {\rm{\textrm{} ET}}} \hspace{0.15cm}&=&\hspace{0.15cm} N_\alpha \left( {1  N_\beta } \right) + N_\beta \left( {1  N_\alpha } \right) \\
P_{2 {\rm{\textrm{} ET}}} \hspace{0.15cm}&=&\hspace{0.15cm} N_\alpha N_\beta \end{eqnarray}
$

(6) 
where $N_\alpha$ and $N_\beta$ are the Mulliken population numbers of alpha and betaspin 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, 3942] that lacks dynamical correlation and has a limited amount of nondynamical correlation; one may wonder whether this situation may affect the accuracy of the presented results. Regarding dynamical correlation, we developed the SLEND/KohnShamdensityfunctionaltheory (SLEND/ KSDFT) [4, 41, 42] that incorporates such effect at the level of timedependent KSDFT [53]. SLEND/KSDFT performs better for molecular geometries and vibrational frequencies than SLEND but inherits from timedependent KSDFT a tendency to overestimate electrontransfer processes [4, 41, 42]. For that reason, we applied SLEND to H$^+$+CO$_2$, a system where electrontransfer processes are important. While lacking dynamical correlation, SLEND has provided good predictions of electrontransfer properties (e.g. integral and differential cross sections) in comparison with experimental data as documented in Refs.[4, 39]. Regarding nondynamical correlation (correct bond breaking), a multiconfiguration 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 nondynamical correlation effects on SLEND results. However, from previous results [4, 39], we can judge that the level of nondynamical correlation in the unrestricted singledeterminantal SLEND wavefunction provides reasonable bondbreaking 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 $$0$\rangle$ and all the singly, doublyexcited, etc., states out of $$0$\rangle$. This superposition creates a meanfield potential under which the nuclei evolve. In this regard, SLEND is similar to other meanfield 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 nonadiabatic 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 zerowidth 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 singularityfree. More details about SLEND in comparison with other dynamical methods are discussed in Refs.[4, 39]
Ⅲ. COMPUTATIONAL DETAILS
The nuclear parameters defining the reactants' initial conditions for the present calculations are depicted in FIG. 1. Nuclei in the CO$_2$ target are labeled: C, O(1), and O(2), respectively, whereas the nucleus making the incoming projectile H$^+$ is labeled H. The 631G basis set [56] is placed on the four moving nuclei to construct the electronic MSO $\left\{ {\psi _h , } \right.$
$\left. {\psi _p } \right\}$. Mediumsize basis sets like 631G have provided accurate results with SLEND without posing a high computational cost [33]. The CO$_2$ target is initially at rest ($\textbf{P}_\textrm{c}$=${\bf{P}}_{{\rm{O}}\left( {\rm{1}} \right)} $=${\bf{P}}_{{\rm{O}}\left( {\rm{2}} \right)} $=${{0}}$), with its center of mass at (0, 0, 0) (cf. FIG. 1) and in its electronic UHFSCF/631G [56] ground state. The angular orientation [$\alpha$, $\beta$] of the CO$_2$ target is determined by the polar $\alpha$ (0$^{\circ}$$\leq$$\alpha$$\leq$180$^{\circ}$) and azimuthal $\beta$ (0$^{\circ}$$\leq$$\beta$$\leq$360$^{\circ}$) angles of its internuclear CO(1) semiaxis (cf. FIG. 1). The H$^+ $ projectile is initially placed at ($b$, 0, 15 a.u.), where $b$ (0$\leq$$b$$\leq$+$\infty$) is the impact parameter. The incoming projectile initially travels along the positive $z$axis with momentum ${\bf{P}}_{\rm{H}}^{{\rm{in}}} $=(0, 0, $ + P_{{\rm{H }}z}^{\textrm{in}} $) corresponding to a kinetic energy $E_{\textrm{Lab}}$=30 eV. Benefiting from the D$_{\infty \textrm{h}} $ symmetry of CO$_2$, initial conditions not replicating symmetryequivalent trajectories can be generated by constraining $\alpha$ and $\beta$ in the ranges: 0$^{\circ}$$\leq$$\alpha$$\leq$180$^{\circ}$ and 0$^{\circ}$$\leq$$\beta$$\leq$90$^{\circ}$, respectively ($e.g.$ the [$\alpha$, $\beta$] and [$\alpha$, 180$^{\circ}$$\beta$] orientations generate equivalent trajectories, cf. FIG. 1). Two sets of trajectories starting from relevant target orientations [$\alpha$, $\beta$] are calculated. Set 1 contains target initial orientations [$\alpha$, 0$^{\circ}$] with 0$^{\circ}$$\leq$$\alpha$$\leq$180$^{\circ}$ in steps of $\Delta\alpha$=15$^{\circ}$ and with $b$ varying: 0.0$\leq$$b$$\leq$9.0 a.u., in steps of $\Delta b$=0.1 a.u. in the 0.0$\leq$$b$$\leq$7.0 a.u. subinterval, and in steps of $\Delta b$=0.2 a.u. in the 7.2$\leq$$b$$\leq$9.0 a.u. one. Set 2 contains collisions of H$^+$ projectile impacting orthogonally on the CO$_2$ axis and arising from targets at initial orientations [90$^{\circ}$, $\beta$] with $\beta$ in two ranges: 0$^{\circ}$$\leq$$\beta$$\leq$15$^{\circ}$ in steps of $\Delta \beta$=1$^{\circ}$, and 20$^{\circ}$$\leq$$\beta$$\leq$90$^{\circ}$ in steps of $\Delta \beta$=5$^{\circ}$, respectively, and with $b$ varying as in Set 1. These initial conditions produce a total of 42 symmetryindependent 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.
Ⅳ. RESULTS AND DISCUSSION
A. Reactive processes and trajectory analysis
Theoretically predicted reactions during the present simulations include:
(ⅰ) Nonchargetransfer scattering (NCTS), where the H$^+$ projectile scatters off after exciting the rotational, vibrational and electronic degrees of freedom of CO$_2$:
$
\begin{eqnarray}
&{\rm{H}}^{\rm{ + }}{\rm{ }}+{\rm{ O}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{C}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{O }}\left( {l_i, {\rm{ }}v_i , {\rm{ }}n_i } \right) \to \\
&{\rm{H}}^{\rm{ + }}+{\rm{ O}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{C}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{O }}\left( {l_f , {\rm{ }}v_f, {\rm{ }}n_f } \right){\rm{ }}\hspace{0.5cm}\left( {{\rm{NCTS}}} \right)
\end{eqnarray}
$

(7) 
(ⅱ) Chargetransfer scattering (CTS), where essentially one electron is transferred to the incoming projectile H$^+$ from the target CO$_2$:
$
\begin{eqnarray}
&{\rm{H}}^{\rm{ + }}+{\rm{ O}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{C}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{O }}\left( {l_i , {\rm{ }}v_i , {\rm{ }}n_i } \right) \to \\
&{\rm{ H}}^{\rm{0}} {\rm{ }}+{\rm{ }}\left[ {{\rm{O}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{C}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{O}}} \right]_{}^ \oplus {\rm{ }}\left( {l_f , {\rm{ }}v_f , {\rm{ }}n_f } \right)\hspace{0.5cm}{\rm{ (CTS)}}
\end{eqnarray}
$

(8) 
NCTS and CTS processes happen together on a given classical nuclear trajectory, with their respective probabilities given by $P_{0 \textrm{} {\rm{ET}}}$ and $P_{1 \textrm{} {\rm{ET}}}$ in Eq.(6).
(ⅲ) Collisioninduced C=O bond dissociation (C=O BD), where one C=O bond of the CO$_2$ molecule is cleaved by the incoming projectile H$^+$:
$
\begin{eqnarray}
&{\rm{H}}^{\rm{ + }}+{\rm{ O}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{C}}\hspace{0.1cm}=\hspace{0.1cm}{\rm{O }}\left( {l_i , {\rm{ }}v_i, {\rm{ }}n_i } \right) \to \\
&{\rm{H}}^q+{\rm{ O}}^{q'}+{\rm{C\hspace{0.1cm}=\hspace{0.1cm}O}}_{}^{q"} \left( {l'_f , {\rm{ }}v'_f , {\rm{ }}n'_f } \right)\hspace{0.3cm}{\rm{ (C}}\hspace{0.15cm}=\hspace{0.15cm}{\rm{O}}\hspace{0.15cm}{\rm{ BD)}} \\
& {\rm{ }}q + q' + q" = + 1
\end{eqnarray}
$

(9) 
Above, the quantum number sets: (${l_i}$, $v_i$, ${n_i }$) and (${l_f}$, $v_f$, ${n_f}$), 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$^+$+CO$_2$$\rightarrow$${\rm{H}}^q $+${{\text{C}}^{{{q}'}}} $+${\rm{O}}^{q"} $+${\rm{O}}^{q'''} $ ($q$+$q'$+$q"$+$q'''$=+1); CO/OH formation: H$^+$+CO$_2$$\rightarrow$${\rm{CO}}^q $+${\rm{OH}}^{q'} $ ($q$+$q'$= +1); and H$^+$ capture: H$^+$+CO$_2$$\rightarrow$HCO$_2$$^+$, 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$^{\circ}$, $\beta$] 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$^{\circ}$, $\beta$] C=O BD reactions are listed in Table Ⅰ, entered by orientations and impact parameter intervals within which the C=O BD reactions occur. The reactive channels listed in Table Ⅰ are also shown in the 2D polar plot in FIG. 2. Additional C=O BD simulations from Set 1 not listed in Table Ⅰ include: [75$^{\circ}$, 0$^{\circ}$]/$b$=1.01.4, and [105$^{\circ}$, 0$^{\circ}$]/$b$=0.3. The extension of the C=O BD reaction in CO$_2$ is smaller than those of analogous collisioninduced reactions in previously investigated targets ($e.g$. H$^+$+H$_2$ [28, 29], +CH$_4$ [46] and +C$_2$H$_2$ [30] at $E_{\textrm{Lab}}$=30 eV) likely due to the higher C=O bond strength.
Table Ⅰ C=O BD) reactions for representative [90$^{\circ}$, $\beta$$^{\circ}$] orientations $vs$. impact parameter $b$.
The two predicted reactions in the H$^+$+CO$_2$ system: C=O BD and NCTS/CTS are depicted in FIGs. 311 for simulations from the initial conditions: [90$^{\circ}$, 6$^{\circ}$]/$b$=0.8 a.u. (C=O BD, FIGs. 35), [90$^{\circ}$, 15$^{\circ}$]/$b$=2.0 a.u. (NCTS/CTS with negligible chargetransfer probability, FIGs. 68), and [90$^{\circ}$, 85$^{\circ}$]/$b$=1.6 a.u. (NCTS/CTS with high chargetransfer probability, FIGs. 911), respectively. FIGs. 311 share common features: FIGs. 3, 6, and 9 depict the four nuclei trajectories projected on appropriate planes ($xz$ plane in FIGs. 3 and 6 and $yz$ plane in FIG. 9) to obtain their clearest visualizations; FIGs. 4, 7, and 10 depict the relative distances between the target nuclear pairs: CO(1), CO(2) and O(1)O(2) $vs$. time; FIGs. 5, 8, and 11 depict the Mulliken population analysis of the alphaspin electrons on the nuclei: C, O(1), O(2) and H $vs$. time.
FIGs. 35 show a C=O reaction from the [90$^{\circ}$, 6$^{\circ}$]/$b$=0.8 a.u. initial conditions. FIG. 3 portrays the incoming projectile H$^+$ (H) initially approaching the CO(1) bond of CO$_2$ and 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 knotlike trajectory while O(1) apparently dissociates from CO$_2$. The C=O BD process suspected in FIG. 3 is verified in FIG. 4. Therein, the three relative distances between pairs of target nuclei: CO(1), CO(2), and O(1)O(2), initially remain constant at the values of the CO$_2$ equilibrium bond lengths; however, at the time of the projectileâ€™s closest approach (about 400 a.u.), both CO(1) and O(1)O(2) monotonically increase and remain so until the final time, whereas CO(2) becomes oscillatory. Those simultaneous motions indicate the dissociation of CO$_2$ into O(1) and CO(2) fragments, the latter arising vibrationally excited. FIG. 5 is also consistent with C=O BD and illustrates concomitant chargetransfer processes: the incoming H$^+$ (H) projectile starts with no electrons but when closest to CO$_2$, it transitorily acquires a Mulliken population number of alphaspin electrons from the CO$_2$ nuclei that is mostly lost after the H nucleus departs from the dissociating CO$_2$. On the other hand, the alphaMulliken populations on C, O(1), and O(2) initially remain constant at their CO$_2$ equilibrium 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 oscillationdriven electron transfers between C and O(2), and the latter dissociating from those two nuclei as indicated by its nonoscillatory alphaMulliken population that reflects no electron transfers through bonding.
FIGs. 68 show a NCTS/CTS reaction with negligible chargetransfer probability from the [90$^{\circ}$, 15$^{\circ}$]/$b$=2.0 a.u. initial conditions. FIG. 6 portrays the incoming projectile H$^+$ (H) initially approaching the O(1) terminus of the CO$_2$ molecule and eventually transferring momentum to O(1), C, and O(2); after that closest approach, the H nucleus repulsively scatters in the rightbackward direction, whereas the CO$_2$ molecule recoils, rotates, and vibrates without visible fragmentation. FIG. 7 confirms a NCTS/CTS reaction without target dissociation because the relative distances: CO(1), CO(2), and O(1)O(2) initially remain constant at the values of the CO$_2$ equilibrium 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 CO$_2$ remains 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 CO$_2$, it transitorily acquires a partial Mulliken population number of alphaspin electrons from the CO$_2$ nuclei that is mostly lost after the H nucleus departs from CO$_2$. FIGs. 911 show another NCTS/CTS reaction from the [90$^{\circ}$, 85$^{\circ}$]/$b$=1.6 a.u. initial conditions but with high chargetransfer probability. The interpretation of FIGs. 911 is similar to that of FIGs. 68. Worth noticing, FIG. 11 indicates that this NCTS/CTS reaction has a high chargetransfer probability. Therein, the incoming H$^+$ (H) projectile starts again with no electrons but when closest to CO$_2$, it acquires a Mulliken population number of alphaspin electrons through a process exhibiting a twojump pattern. Unlike the previous cases, that electron transfer becomes permanent after the H nucleus and CO$_2$ separate.
B. Scattering angles, rainbow and glory scattering angles
The scattering angle $\theta$ of the outgoing projectile ${\rm{H}}^q $ (0$\leq$$q$$\leq$+1) with respect to the incoming $z$direction is
$
\begin{eqnarray}
\sin \theta \hspace{0.15cm}&=&\hspace{0.15cm} \frac{{\left( {P_{{\rm{H }}x}^{{\rm{out }}2} + P_{{\rm{H }}y}^{\textrm{out}2} } \right)^{1/2} }}{{\left( {P_{{\rm{H }}x}^{\textrm{out}{\rm{ }}2} + P_{{\rm{H }}y}^{\textrm{out}{\rm{ }}2} + P_{{\rm{H }}z}^{\textrm{out}{\rm{ }}2} } \right)^{1/2} }}\\
&{\rm{ }}\left( {0 \le \theta \le 180^\circ } \right)
\end{eqnarray}
$

(10) 
where $\left\{ {P_{{\rm{H }}i}^{\textrm{out}{\rm{ }}} } \right\}$ ($i$=$x$, $y$, $z$) are the outgoing projectile momentum components at the final time. Scattering angle functions $\theta$($b$) $vs$. impact parameters $b$ from selected [90$^{\circ}$, $\beta$] target orientations in Set 2 are shown in FIG. 12 to illustrate the effect of the orientation angle $\beta$ on $\theta$($b$). In addition, values and locations of critical angles $\theta$($b$) (rainbow and glory scattering angles [57]) are listed in Table Ⅱ for selected [$\alpha$, $\beta$] orientations in Sets 1 and 2.
Table Ⅱ Position and values of primary (classical and uniform Airytreated) and secondary (classical) rainbow scattering angles, and positions of glory scattering angles from selected target orientations [$\alpha$, $\beta$].
In FIG. 12, $\theta$($b$) from the [90$^{\circ}$, 0$^{\circ}$] orientation exhibits a primary rainbow angle $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{ CL}}} $=$\theta\left( {b_R^\textrm{I} } \right)$=12.20$^{\circ}$ as a maximum ([d$\theta$($b$)/d$b$]$_{b = b_R^\textrm{I}} $) at the high impact parameter $b_R^\textrm{I} $=4.69 a.u. (cf. Table Ⅱ; 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 projectiletarget potential energy functions [57]. Therefore, accurate predictions of rainbow angles features will be a measure of the correctness of the projectiletarget interactions with SLEND. Away from the rainbow area, $\theta$($b$) becomes a decreasing function at impact parameters $b$ somewhat lower than $b_R^\textrm{I} $ and reaches a zero value $\theta _G $ (${b_G }$)=0$^{\circ}$ known as the glory scattering angle [57] at the intermediate impact parameter ${b_G }$ (cf. Table Ⅱ). Other zero and nearly zero scattering angles $\theta _S $$\rightarrow$0$^{\circ}$, 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$$ <$$b_G $, $\theta$($b$) becomes an increasing function. The observed features in $\theta$($b$) from the [90$^{\circ}$, 0$^{\circ}$] orientation result from the projectiletarget interactions that are overall repulsive and overall attractive in the intervals (1.0 a.u., $b_\textrm{G}$) and ($b_G$, $+\infty$), respectively, with the primary rainbow angle $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{ CL}}} $ corresponding to a maximum in the projectiletarget attraction; the glory scattering $\theta _G \left( {b_G } \right)$=0$^{\circ}$ occurs at $b$=$b_\textrm{G}$, i.e. at the transition from repulsive to attractive interactions where the net effect is no deflection. The $\theta$($b$) from the remaining [$\alpha$, 0$^{\circ}$] orientations exhibit exactly the same features as those discussed from the [90$^{\circ}$, 0$^{\circ}$] one; for clarity, those functions are not plotted in FIG. 12 but a selection of their primary rainbow and glory scattering angle data is listed in Table Ⅱ.
The scattering angle functions $\theta$($b$) from the [90$^{\circ}$, 0$^{\circ}$$ <$$\beta$$\leq$50$^{\circ}$] orientations in FIG. 12 also display primary rainbow angles $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{ CL}}} $ that do not differ much in value from the [90$^{\circ}$, 0$^{\circ}$] $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{ CL}}} $ and that tend to appear at lower impact parameters $b_R^\textrm{I} $ (cf. Table Ⅱ). More importantly, the glory scattering angle $\theta _G \left( {b_G } \right)$=0$^{\circ}$ from the orientation [90$^{\circ}$, 0$^{\circ}$] is replaced in the [90$^{\circ}$, 0$^{\circ}$$ <$$\beta$$\leq$50$^{\circ}$] ones for secondary rainbow angles $\theta _R^{\textrm{II}\hspace{0.05cm}{\rm{ CL}}} $=$\theta$$\left( {b_R^{\textrm{II}} } \right)$$\neq$0$^{\circ}$ that are minima in $\theta$(b) with [d$\theta$($b$)/d$b$]$_{b=b_R^{\textrm{II}} } $=0 at lower impact parameters $b_R^{\textrm{II}} $$ <$$b_R^\textrm{I} $ (cf. FIG. 12 and Table Ⅱ). Due to the initial reactants' symmetry, the four nuclei in [$\alpha$, 0$^{\circ}$] and [90$^{\circ}$, 90$^{\circ}$] simulations remain on the $xz$ plane of the initial conditions (cf. FIG. 1) during the entire scattering process. The lower symmetry in the [90$^{\circ}$, 0$^{\circ}$$ <$$\beta$$\leq$50$^{\circ}$] initial conditions causes the nuclei not to remain on the initial $xz$ plane but scatter off it; a secondary rainbow angle approximately corresponds to a maximum offplane scattering [57], a phenomenon also inferred experimentally and predicted with SLEND and other theories in other protonmolecule reactions ($e.g$. H$^+$+H$_2$ [5, 28, 29] and H$^+$+C$_2$H$_2$ [15, 30]).
An interesting effect displayed in FIG. 12 is that the values of the primary rainbow $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{ CL}}} $ and secondary rainbow/glory $\theta _R^{\textrm{II}\hspace{0.05cm}{\rm{ CL}}} $/$\theta _G^{} $ angles from the [90$^{\circ}$, $\beta$] orientations get closer as $\beta$ varies from 0$^{\circ}$ to 50$^{\circ}$ and so do their respective impact parameters $b_R^\textrm{I} $ and $b_R^{\textrm{II}} $/$b_G $. At some $\beta$ (50$^{\circ}$$ <$$\beta$$ <$55$^{\circ}$), the approaching primary (maximum) and secondary (minimum) rainbow angles finally coalesce into an inflexion point so that their rainbow angle features vanish. The $\theta$($b$) from the additional [90$^{\circ}$, 55$^{\circ}$$\leq$$\beta$$ <$90$^{\circ}$ orientations lie beyond that coalescence point and display no rainbow angle features (cf. FIG. 12). Like in H$^+$+C$_2$H$_2$ [30], the final value of $\beta$=90$^{\circ}$ brings back the rainbow angles features but in a complex pattern exhibiting: (ⅰ) a maximum primary rainbow angle $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{ CL}}} $=0.92$^{\circ}$ at $b_R^\textrm{I}$=5.08 a.u., (ⅱ) a maximum (not minimum) secondary rainbow angle $\theta _R^{\textrm{II}\hspace{0.05cm}{\rm{ CL}}} $=0.50$^{\circ}$ at $b_G^{\textrm{II}}$=2.95 a.u., and (ⅲ) two glory scattering angles at $b_G^\textrm{I}$=2.57 a.u. and $b_G^{\textrm{II}}$=3.56 a.u., respectively (cf. FIG. 12). Rainbow scattering angles manifest in conspicuous features of both the experimental and calculated differential cross section as discussed below.
C. Differential cross sections
The most important result from the Toennies group experiment on H$^+$+CO$_2$ at $E_{\textrm{Lab}}$=30 eV is the total NCTS and CTS differential cross sections (DCS) [13]. In the present context, quantum inelastic DCS [${\textrm{d}\sigma ^{i \to f} \left( \theta \right)}/\textrm{d}\Omega$] from initial ($i$) to final ($f$) states ($i$$\rightarrow$$f$) can be calculated as [57]:
$
\begin{eqnarray}
\left[ {\frac{{\textrm{d}\sigma ^{i \to f} \left( \theta \right)}}{{\textrm{d}\Omega }}} \right]\hspace{0.15cm} &=& \hspace{0.15cm}\frac{{k_f }}{{k_i }}\left {f^{i \to f} \left( \theta \right)} \right^2 \\
&=&\hspace{0.15cm} \frac{1}{{4k_i^2 }}\left {\sum\limits_{l = 0}^\infty {\left( {2l + 1} \right){\rm{ }}} T_{\textrm{quantum}}^{i \to f} \left( l \right)P_l \left( {\cos \theta } \right)} \right^2\\
\end{eqnarray}
$

(11) 
where $f^{i \to f} \left( \theta \right)$ is the scattering amplitude,
$k_i \left( {k_f } \right)$ is the projectile initial (final) wave vector, $l$ is orbital angular momentum quantum numbers, $T_{\textrm{quantum}}^{i \to f} \left( l \right)$ is the $T$matrix, and $P_l \left( {\cos \theta } \right)$ is Legrende polynomials. Since SLEND employs classical nuclear dynamics, the quantum DCS in Eq.(11) should be evaluated using semiclassical techniques [4, 57, 58] as follows. When $\theta$($b$) is not close to rainbow angles, the stationary phase approximation [57] gives the semiclassical DCS $\left[ {\textrm{d}\sigma _{\textrm{CS}}^{i \to f} \left( \theta \right)}/\textrm{d}\Omega\right]$ as [4, 57, 58]:
$
\begin{array}{l}
\left[ {\frac{{{\rm{d}}\sigma _{{\rm{SC}}}^{i \to f}\left( \theta \right)}}{{{\rm{d}}\Omega }}} \right] = {\left {\sum\limits_{m = 1}^N {f_m^{}\left( \theta \right)} } \right^2}\\
= {\left {\sum\limits_{m = 1}^N {\frac{{\sqrt {{b_m}} T_{{\rm{SLEND}}}^{i \to f}\left( {{b_m}} \right)}}{{\sqrt {\sin \left[ {\theta \left( {{b_m}} \right)} \right]{{\left[ {{\rm{d}}\theta \left( b \right)/{\rm{d}}b} \right]}_{b = {b_m}}}} }}} } \right^2}
\end{array}
$

(12) 
where $l$$\approx$$k_ib$, $T_{\textrm{quantum}}^{i \to f} $=$T_{\textrm{SLEND}}^{i \to f} $, and $m$(1$\leq$$m$$\leq$$N$) labels the $N$ classical projectile trajectories (branches) that from impact parameters $b_m$ lead to the same scattering direction $\theta$; in the present simulations, the maximum values of $N$ are: 3 from [$\alpha$, 0$^{\circ}$] and [90$^{\circ}$, 0$^{\circ}$$ <$$\beta$$\leq$50$^{\circ}$], 1 from [90$^{\circ}$, 55$^{\circ}$$ <$$\beta$$\leq$85$^{\circ}$], and 5 from [90$^{\circ}$, 90$^{\circ}$], respectively (cf. FIG. 12). In addition,
$
\begin{eqnarray}
T_{\textrm{SLEND}}^{i \to f} = \sqrt {P_{\textrm{SLEND}}^{i \to f} } \exp \left( {iS} \right)
\end{eqnarray}
$

(13) 
where $P_{\textrm{SLEND}}^{i \to f} $ is the probability of the process from a SLEND simulation and $S$ is the quantum action (cf. discussion in Section II). In the cases of NCTS and CTS ($i.e.$ in the cases of zeroelectron and oneelectron total transfers), $P_{\textrm{SLEND}}^{i \to f} $ is given by $
P_{0 {\rm{\textrm{}ET}}}$ and $P_{1{\rm{ \textrm{} ET}}}$ in Eq.(6), respectively. Notice that $S$ in $T_{\textrm{SLEND}}^{i \to f} $ produces a quantum oscillatory pattern in the semiclassical DCS of Eq.(12) via coherent sum of scattering amplitudes. If that sum is performed incoherently, a (predominantly) classical DCS $
\left[{\textrm{d}\sigma _{\textrm{CL}}^{i \to f} \left( \theta \right)}/\textrm{d}\Omega\right]$ is obtained:
$
\left[ {\frac{{{\rm{d}}\sigma _{{\rm{CL}}}^{i \to f}\left( \theta \right)}}{{{\rm{d}}\Omega }}} \right] = \sum\limits_{m = 1}^N {\frac{{{b_m}{{\left {T_{{\rm{SLEND}}}^{i \to f}\left( {{b_m}} \right)} \right}^2}}}{{\sin \theta \left {{{\left[ {{\rm{d}}\theta \left( b \right)/{\rm{d}}b} \right]}_{b = {b_m}}}} \right}}}
$

(14) 
that retains quantum probabilities through $T_{{\rm{SLEND}}}^{i \to f}\left( {{b_m}} \right){{\rm{}}^2}$ but lacks the mentioned oscillatory pattern; that classical DCS is the lowestlevel 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, ${\left[ {{\rm{d}}\theta \left( b \right)/{\rm{d}}b} \right]_{b = {b_m}}}$, 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$\theta$($b$)/d$b$=0$^{\circ}$; more sophisticated semiclassical techniques are required in those situations. The uniform Airy [57, 59] approximation is the most rigorous semiclassical treatment of an isolated rainbow angle that occurs in all [$\alpha$, 0$^{\circ}$] simulations. There, the scattering amplitudes $f_2$($\theta$) and $f_3$($\theta$) from two attractive trajectories near the primary rainbow angle direction (Eq.(12)) are replaced with a combined single amplitude $f_{23}^{{\rm{U}}{\rm{. Airy}}} \left( \theta \right)$:
$
\begin{array}{l}
f_{23}^{{\rm{U}}.{\rm{Airy}}}\left( \theta \right) = {\pi ^{1/2}}\exp \left[ {i\left( {{A_{23}}  \frac{1}{4}\pi } \right)} \right]\left\{ {\left[ {\left {\frac{{{\rm{d}}{\sigma _{{\rm{QC}}}}\left( \theta \right)}}{{{\rm{d}}\Omega }}} \right_2^{1/2} + \left {\frac{{{\rm{d}}{\sigma _{{\rm{QC}}}}\left( \theta \right)}}{{{\rm{d}}\Omega }}} \right_3^{1/2}} \right]\xi _{23}^{1/4}Ai\left( {  {\xi _{23}}} \right)  } \right.\\
\left. {i\left[ {\left {\frac{{{\rm{d}}{\sigma _{{\rm{QC}}}}\left( \theta \right)}}{{{\rm{d}}\Omega }}} \right_2^{1/2}  \left {\frac{{{\rm{d}}{\sigma _{{\rm{QC}}}}\left( \theta \right)}}{{{\rm{d}}\Omega }}} \right_3^{1/2}} \right]\xi _{23}^{  1/4}Ai'\left( {  \xi _{23}^{}} \right)} \right\}
\end{array}
$

(15) 
where $A_{23} \left( \theta \right)$ and $\xi _{23}^{} \left( \theta \right)$ are functions of the action $S$ along those trajectories [59], and $Ai$ and $Ai'$ are the Airy function and its first derivative [60]. The uniformAiry DCS from Eq.(15) evolves uniformly to the DCS in Eq.(12) when $\theta$ is away from the rainbow angle and to the transitional Airy approximation [59] DCS when $\theta$ is close to the rainbow angle. The uniform Pearcey [59] approximation is the most rigorous semiclassical treatment of the two adjacent rainbow angles that occur in all $[90^\circ {\rm{, }}$
$0^\circ $$ <$$\beta$$\leq$50$^\circ $] simulations but it is too cumbersome to apply. However, that approximation approaches the uniformAiry 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 [3032] is adopted herein. That approach involves approximating the Pearcey scattering amplitude $f_{123}^{{\rm{U}}{\rm{. Pearcy}}}$($\theta$)$\approx$$f_1$($\theta$)+$
f_{23}^{{\rm{U}}{\rm{. Airy}}}$($\theta$) and $f_{123}^{{\rm{U}}{\rm{. Pearcy}}}$($\theta$)$\approx$$f_3$($\theta$)+$f_{12}^{{\rm{U}}{\rm{. Airy}}}$($\theta$) near the rainbow angles $\theta _R^{\textrm{I}{\rm{ CL}}}$/$b_R^\textrm{I} $ and $\theta _R^{\textrm{II}{\rm{ CL}}} $/$b_R^{\textrm{II}} $, respectively, and switching between those amplitudes at the inflexion point $b_{}^{\textrm{IP}} $ of $\theta$($b$), $b_R^{\textrm{II}} $$ <$$b_R^{\textrm{IP}} $$ <$$b_R^{\textrm{I}} $. The above uniform approximations are valid for rainbow angles not close to 0$^{\circ}$; therefore, treatments of a few lowlying rainbow angle singularities are not possible with those techniques and are left uncorrected ($e.g$. the two rainbow angles: $\theta _R^{\textrm{I CL}} $=0.92$^{\circ}$ and $\theta _R^{\textrm{II CL}}$=0.5$^{\circ}$ in the [90$^{\circ}$, 90$^{\circ}$] simulation). Contrastingly, the [90$^\circ$, 55$^{\circ}$$ <$$\beta$$\leq$85$^{\circ}$] 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$^+$+CO$_2$ at $E_{\textrm{Lab}} $=30 eV from the individual target orientation [90$^{\circ}$, 0$^{\circ}$] are shown in FIG. 13, both at the classical (Eq.(14)) and the uniformAiry (Eqs. (12) and (15)) levels, and along with the experimental result [13]. The latter was reported unnormalized but is normalized in FIG. 13 to allow comparisons. The normalization criterion involves attaining a minimum deviation between experimental and uniformAiry results over all the investigated scattering angles $\theta$. The classical DCS shows an unphysical rainbow angle singularly at $\theta$=$\theta _R^{\textrm{I}{\rm{ CL}}} $=12.20$^{\circ}$, while the uniformAiry CS DCS shows a bounded rainbow angle maximum near $\theta$=10$^{\circ}$ in agreement with the experimental result. Except for that singularity, the classical DCS agrees fairly with the measured data. The uniformAiry DCS agrees better and even displays an oscillatory pattern produced by the coherent sums of the action terms in $T_{\textrm{SLEND}}^{i \to f}$, Eq.(13), and the Airy terms, Eq.(15), [57]. Those genuine quantum oscillations are artificially absent in the classical DCS (due to incoherent sum and lack 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 uniformAiry DCS as shown below. The comparison between classical (Eq.(14)) and uniformAiry (Eqs. (12) and (15)) DCS in FIG. 14 demonstrates the better performance of the latter; therefore, only uniformAiry 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 scattering angle $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{U.Airy}}} $$ <$$\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{ CL}}} $; $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{U.Airy}}} $ values provide a more accurate prediction of the experimental rainbow angle than $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{ CL}}} $; some $\theta _R^{\textrm{I}\hspace{0.05cm}{\rm{U.Airy}}} $ from selected orientations are listed in Table Ⅱ.
The final predicted total DCS for both NCTS and CTS processes in H$^+$+CO$_2$ at $E_{\textrm{Lab}}$=30 eV are averages of the DCSs from the simulated target orientations; they are shown in FIG. 14 at the uniformAiry level along with their corresponding experimental results [13]. The unnormalized 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 $\theta$. Conversely, the calculated CTS DCS in FIG. 14 agrees well with its measured counterpart at large scattering angles $\theta $$\geq$8$^{\circ}$ but it is lower in value than that counterpart at small scattering angles $\theta$$ <$8$^{\circ}$. Notably, both calculated NCTS and CTS DCS exhibit bounded rainbow angle features at $\theta$=10$^{\circ}$ in very good agreement with those in their experimental counterparts [13]. As stated previously, rainbow angle features are related to parameters defining model projectiletarget potential energy functions [57]. Therefore, the current accurate prediction of rainbow angles features is a measure of the correctness of the projectiletarget interactions with SLEND. The positions of the primary rainbow angle features in the calculated DCS are averages of the $\theta _R^{\textrm{I}{\rm{U}}.{\rm{ }}\textrm{Airy}} $ values per target orientation weighted with the individual DCS intensities; those positions distribute closely around $\theta$=10$^{\circ}$ (cf. Table Ⅱ) 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 $\theta$=7$^{\circ}$ [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 Ⅱ); 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 Ⅴ.
Ⅴ. CONCLUSION
The SLEND method [4, 39] has been applied to the H$^+$+CO$_2$ reactive system at $E_{\textrm{Lab}}$=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 singledeterminantal Thouless wavefunction [44]. The present investigation comprises a total of 3402 reaction simulations from 42 independent CO$_2$ orientations. The first part of this investigation (Section IV. A) focused on predicting the chemical reactions occurring in the H$^+$+CO$_2$ system; 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: nonchargetransfer scattering (NCTS), chargetransfer 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 chargetransfer processes (cf. FIGs. 311). H$^+$+CO$_2$ at $E_{\textrm{Lab}}$=30 eV turned out to be less reactive in terms of both extension and variety of reactions than similar protonmolecule systems simulated with SLEND (e.g. H$^+$+H$_2$ [28, 29], +CH$_4$ [46] and +C$_2$H$_2$ [30] at $E_{\textrm{Lab}}^{} $=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 CO$_2$ mass. The calculated scattering functions revealed complex scattering patterns in terms of rainbow and glory angles. Simulations from the [$\alpha$, 0$^{\circ}$] orientations ($i.e.$ those with the initial H$^+$ and CO$_2$ on the same plane) and the [90$^{\circ}$, 90$^{\circ}$] one ($i.e$. that with the initial projectile hitting orthogonally to a plane containing CO$_2$) keep the H$^+$ projectile on the same plane during scattering due to the symmetry of the projectiletarget 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$^{\circ}$, 0$^{\circ}$$ <$$\beta$$ <$90$^{\circ}$] (i.e. those where the above planar symmetries between H$^+$ and CO$_2$ at initial time are missing) do not keep the H$^+$ projectile on the same plane during scattering. As a result, simulations from the [90$^{\circ}$, 0$^{\circ}$$ <$$\beta$$\leq$50$^{\circ}$] 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 CO$_2$ and predicted with SLEND and other theories in other protonmolecule reactions (e.g. H$^+$+H$_2$ [5, 29] and H$^+$+C$_2$H$_2$ [15, 30]). Moreover, in simulations from the [90$^{\circ}$, 55$^{\circ}$$\leq$$\beta$$ <$90$^{\circ}$] 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$^{\circ}$$ <$$\theta $$\leq$18$^{\circ}$. On the other hand, SLEND CTS DCS agrees well with its experimental counterpart [13] at large scattering angles $\theta$$\geq$8$^{\circ}$ but is lower in value than its experimental counterpart at small scattering angles $\theta$$ <$8$^{\circ}$. 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 $\theta$=10$^{\circ}$ as do their experimental counterparts [13]. Rainbow angle features are related to the projectiletarget interactions; therefore, this good rainbow angle prediction indicates that SLEND reproduce proper projectiletarget 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 Ⅱ).
Ⅵ. 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).