The role of multiple ionization of H2O in heavy ion collisions

Collisional ionization processes involving H2O molecules and C6+, O8+, Si13+ ions are studied by means of the classical trajectory Monte Carlo method using molecular orbital calculations to define the ionization stages of the water molecule. Net total and single-differential cross sections in energy and angle are obtained by using a newly developed model that goes beyond the commonly applied one-active electron approximation. This model allows us to access the fraction of electron emission arising from single and multiple electron ionization. Calculated cross sections are contrasted and benchmarked against available experimental data at impact energies in the MeV/u range. The present results highlight the important role of multiple ionization in the emission of electrons where we find the majority of electrons emitted with energies greater than ~50 eV arise from multiple ionization collisions.


Introduction
Ion beams utilizing protons, alpha particles, and carbon ions are being used for radiotherapy (Njkoo et al 2006, Suit et al 2010, Ebner et al 2017, Sokol et al 2017, Solov'yov 2017, with the number of such facilities that is rapidly expanding. The number of facilities under construction or in planning stages worldwide will double by 2020 from those operational in 2016 (Dosanjh 2018). This growth is the result that ion-beam therapy delivers a localized dose of energetic electrons that allows the treatment of deep seated tumours while minimizing the irradiation of the surrounding healthy tissue which is particularly important in areas near the brain or spinal cord. As such, this technique complements electron-beam impact therapies (particularly suited for superficial tumours, like skin cancer) and x-ray irradiation. In contrast to electrons and photons, heavy ions transfer most of their energy at the point where they almost stop in the spatial region known as the Bragg peak.
The ion-beam irradiation process can be considered a multi-scale process. Collision processes leading to electron production via direct ionization, secondary electron transport and Auger processes occur in timescales (<10 −15 s) and represent the start of the process. Fragmentation and recombination characterize the subsequent chemical stage (10 −15 -10 −12 s). At this point, the biological damage over the target zone has been achieved and biochemical and biological stages follow. In these stages the affected cells will try to repair their damage with variable success depending on the extent and origin of the damage. At present, it is widely accepted that double strand breaks in DNA leads to enhanced cell death. However, the physical mechanisms that lead to double and multiple strand breaks remain unknown. Therefore, a clear understanding of the electron emission mechanisms acting during the collision process seems, from a physical point of view, a necessary step to improve a basic understanding on this topic.
On the experimental side, considerable data exists for the ionization of water by protons from the pioneering work of Rudd, Toburen and others during the 1970s and 1980s (Toburen et al 1980, Rudd et al 1985, Bolorizadeh and Rudd 1986, Rudd et al 1992. Similar investigations have now been extended to a solid water target by experimentalists at the East Carolina University (Toburen et al 2010, Shinpaugh et al 2012. Studying water is of considerable importance since it comprises approximately 60% of the human body. Only very recently have a limited set of data become available for the heavier ions important to radio-therapy and deep space applications (Ohsawa et al 2005, Dal Capello et al 2009, Bhattacharjee et al 2017, Bhattacharjee et al 2018. These data play an The role of multiple ionization of H 2 O in heavy ion collisions N Bachi 1,3 , S Otranto 1 , G S Otero 1 and R E Olson 2 important role as benchmarks to any calculation. However, none of these data display ionized electron spectra that consider the importance of multiple ionization. Moreover, from the theoretical point of view most of the calculations have been performed within the one-active electron formulation (Illescas et al 2011, Liamsuwan and Nijkoo 2013, Bhattacharjee et al 2017, Bhattacharjee et al 2018 and have mainly directed at the net cross section level with multiple ionization estimates limited to double electron processes. Similar studies have been carried out for molecular targets of interest to the ion-therapy program such as the four nitrogen bases of DNA. In these cases, atomic cross sections for the conforming atoms have been added up with population coefficients provided by quantum mechanical structure codes (Galassi et al 2012, Champion et al 2015, Sarkadi 2016. Up to our knowledge, the only theoretical model that considered a multielectronic treatment of the collision process was the classical trajectory Monte Carlo-classical overbarrier model introduced by Abbas et al (2008). In this model, virtual electrons are created as the projectile approaches the target and the overbarrier condition to remove multiple electrons is fulfilled. Molecular targets are simplified to single-center atomic-like targets, and the virtual electrons interact with the ionic center and the other electrons by means of pure Coulomb potentials.
In this paper we present cross sections for collisions of ions of charge states Z P = 6+, 8+ and 13+ with H 2 O target for an energy relevant to the Bragg peak where maximum damage is done to biological media and for which specific experimental data are available. The ions chosen were those commonly used for hadron-cancer radiotherapy. Our calculations were performed by using molecular structure data for H 2 O and its ions. The scattering calculations were done using the classical trajectory Monte Carlo method that includes all valence electrons for the H 2 O molecule, and as such, this calculation goes beyond the one-electron level so that the neglected important multiple ionization processes can be addressed. Moreover, since the energies of interest near the Bragg peak energy have a large Bohr parameter, Z p /v p > 0.5, easily applied one-electron perturbative theoretical methods such as the Born approximation and variants thereof are usually considered to be outside their region of validity.
Present calculations predict multiple ionization will be an important mechanism for the ejected electron emission, particularly for carbon-beam irradiation. The multiple electron removal in a single collision can represent an important fraction of the ionized electrons and dominates the ejected electron spectra at specific energies and angles such as around the binary peak. These energetic electrons will also have the longest range in biological materials.

Theory
To predict the role of multiple ionization, we have developed an adaptive twelve-body CTMC model in its microcanonical formulation which we denote Ad-CTMC. Hamilton's equations, are numerically solved for the projectile, the three molecular centers and eight non-interacting electrons by means of a step adaptive Runge-Kutta-Gill method. The three nuclear charges of the molecular target are arranged according to its equilibrium configuration, with a distance between hydrogen and oxygen atoms of r O-H = 1.814 a.u. and an angle between the two O-H bonds equal to 106.5° quantities that are kept fixed throughout the simulation. The eight electrons interact with the three centers via a multi-center potential given by, which represents the sum of terms corresponding to the interaction of each electron with the O and the H centers: This multi-center potential describes the interaction among one electron and the mean field created by the nuclei and the rest of the electrons of the molecule. In the above equation, r represents the electron distance to the respective molecular center. The corresponding parameters N O = 8.18, N H = 1.91, α O = 3.06 and α H = 5.07 are determined by a least squares fitting procedure over the molecular anisotropic potential derived from the molecular orbitals ϕ j (r) provided by Moccia (1964) and that is given by, In the above expression, Z l represents the nuclear charges of the atom conforming the molecule. The sum index j runs over the valence molecular orbitals (1b1, 3a1, 1b2 and 2a1) and N ij = 2 for i ≠ j and N ij = 1 for i = j. The present least squares procedure provides the best overall fit that minimizes the overall dispersionχ 2 = χ 2 1b1 + χ 2 3a1 + χ 2 1b2 + χ 2 2a1 in a cubic box of 40 a.u. per side centered in the O-atom.
Throughout the collision, the projectile interacts with the eight electrons, the two protons and the O 6+ core. For the latter, we employ the potential model developed by Green et al from Hartree-Fock calculations (Green et al 1969) and later on generalized by Garvey et al (1975): This potential model consists in the sum of the asymptotic Coulomb potential in which the nuclear charge Z j is screened by N − 1 electrons, and a short range term which improves the description of the inner region where the electron screening is weaker. Using this functional form, these authors tabulated the ξ and η parameters that minimize the energy for different atoms and positive ions with nuclear charges values up to 54. The corresponding parameters in the present case for the O 6+ ion are Z O = 8, N = 3, ξ = 7.046 and η = 5.161 (Garvey et al 1975).
As a result, the projectile asymptotically sees the O-core as a +6 charge. As the relative distance decreases, this charge-value increases up to +8. Electrons are initially sorted by randomly sorting the polar angles (θ r , ϕ r ) which characterize the position vector r of the electron with respect to the O-atom center. Then, the maximum radial distance r 0 at that particular orientation is obtained from the classical return point condition E i − V e−H2O (r 0 , θ r , ϕ r ) = 0. This distance r 0 varies according to the binding energy E i of each electron. In our case, these correspond to the ionization potentials of the 1b1, 3a1, 1b2 and 2a1 valence molecular orbitals which are indicated in table 1.
The values of r for this particular orientation are then confined to the interval 0 < r < r 0 and are obtained following the standard microcanonical procedure (Reinhold and Falcón 1986). Energy conservation is used to determine the electrons momenta p i . At this point the target has been initialized and the simulation stage begins.
As the projectile approaches the target interacts separately with the particles conforming the H 2 O molecule (8 electrons plus the two H + and the O 6+ ion). During every step of the simulation, the energy of each of the electrons relative to the target ion center of mass is checked. If one of the electrons acquires a positive energy, the remaining bound electrons are automatically re-sorted with the ionization potentials corresponding to H 2 O + ion by assuming a vertical transition for the molecule. These electrons retain their momenta and are spatially re-located to match the new ionization potential. By so doing, we can keep track of the momentum transferred by the projectile to the recoiling target ion. The N O parameter of the re-sorted electrons is also diminished in one unit to let them see an asymptotic +2 charge. This procedure is repeated as new electrons are emitted from the target. In table 1 we show the electronic structure of the H 2 O q+ ions and the ionization potentials of its electrons as provided by the GAMESS code using a linear combination of Gaussian functions centered at the different atomic centers conforming the molecule (the denominated 3-21G base) in an unrestricted Hartree-Fock (UHF) treatment (Schmidt et al 1993).
It is important to notice that the present methodology pretends to model the multiple ionization dynamics in a collisional time scale (~10 −17 s) that is much shorter than the time scale needed for molecular fragmentation or dissociation (~5 × 10 −15 s). In this sense, the (H 2 O) q+ ionic states in table 1 do not intend to mean final stable ionic states, but refer to the target charge state during the sequential vertical transitions the target goes through. At the large impact energies considered in this work, the direct removal of q electrons take place in a time scale in which the interatomic distances in the molecule have no time to re-adapt and throughout the simulation they are kept fixed. Recent studies have married Coulomb explosion models together with CTMC dynamics to gain insight on the energetic spectra of the resulting fragments following collisions of H 2 O molecules with heavy ions (Wolff et al 2016). These type of studies are beyond the scope of the present work.
Once the projectile recedes and causes no further significant perturbation, the integration of the Hamilton's equations is stopped and the emitted electrons' energies and angles are calculated from the asymptotic value of their momenta. This set of data is collected and used to calculate the single differential cross sections (SDCS) in energy and angle (Reinhold and Falcón 1986).
CTMC models which contain n-electrons have been previously used in the literature (Olson et al 1987, 1989, Cariatore and Otranto 2013. Collision models based on the independent electron (IE) scheme, consider electrons only bound with fixed ionization potentials which represent the ground state valence molecular orbital. In this case, the electron density seen by the projectile is correct, but the energy deposition required for multiple electron removal is underestimated. A scheme which tries to correct this deficiency is the one that considers sequential electrons (SEQ). In this case, multiple-electron removal is achieved with the correct energy deposition by the projectile. However, the electron density of the molecular ground state is underestimated in this case. The present methodology is expected to reproduce the correct electron density for the valence molecular orbitals, and at the same time require the proper energy deposition by the projectile for multiple electron removal.
To analyze the capabilities of the methodology hereby introduced, in the next section we compare the results obtained in pure ionization studies to the available experimental data.

Results
In figure 1 we show the single-differential cross sections (SDCS) in energy and angle for 4 MeV/u C 6+ and Si 13+ and 3.75 MeV/u O 8+ collisions on H 2 O. Present calculations are benchmarked against the experimental and the CDW-EIS data of Bhattacharjee et al (2017Bhattacharjee et al ( , 2018. As expected, charge exchange and transfer ionization have been found to be negligible at these large impact energies. At low electron emission energies, and for the three projectiles under consideration, the spectrum is clearly dominated by the single ionization process. This energy region is mainly populated with electrons arising in soft, large impact parameter collisions. Present theoretical results are expected to underestimate those contributions at these large impact energies based on their classical nature as will be detailed below. As we move to the intermediate to large emission energy regions we find that the multiple ionization dominates beginning at energies of 70 eV, 50 eV and 40 eV respectively. Regarding the angular distributions of the emitted electrons, single ionization provides the leading contribution although multiple ionization is quite significant. The binary peak structure, clearly visible in the 60°-90° range, is present in both mechanisms for C 6+ . However, as the projectile charge increases we observe that the multiple ionization process becomes responsible of the peak structure. In this case, single ionization provides a smooth step down behavior for angles greater than about 80°. Experimental results from Bhattacharjee et al (2017) suggested that the total net cross sections σ net measured, (σ net = σ 1ion + 2σ 2ion + 3σ 3ion + . . .) followed a scaling law ~Z n p with n < 2. In fact, it was suggested for n a value of 1.7 ± 0.1. Present Ad-CTMC results for σ net fitted over a wide range of Z P values predict an n-value of 1.7 in good agreement with the experimental observations. This is shown in figure 2, where we notice that light projectiles evidence a slight deviation from this functional form.
In order to examine more closely the implications of this scaling obtained at the total cross section level, in figure 3(a) we show our scaled SDCS in energy for the three projectiles under consideration. These cross sections are compared to the reported SDCS. Since the O 8+ data have been measured at an impact energy of 3.75 MeV/u, the cross sections are expected to be larger than those corresponding to 4 MeV/u. Hence, we estimate a multiplicative correction factor of 0.58 for the O 8+ data by interpolating the reported C 6+ and Si 13+ σ net cross sections at 4 MeV/u. The separate contributions of the single and multiple ionization mechanisms are analyzed in figure 3(b). It can be seen that both mechanisms exhibit a clear dependence on Z P despite the fact that these differences are washed out once added up. These results indicate that multiple ionization gains relevance against single ioniz ation as Z P increases, especially in collision events which lead to energetic electron emission.
The impact parameter dependence for the single and multiple ionization processes predicted by the Ad-CTMC model is shown in figure 4 for the three projectiles under study. Single and multiple ionization are seen to dominate different impact parameter regions. The former can be associated to larger impact parameters for which the projectile momentum transfer is low. In contrast, the latter is associated to inner impact parameters for Table 1. Molecular orbitals and ionization potentials for H 2 O q+ ions following vertical transitions calculated with GAMESS using the 3-21G UHF (Schmidt et al 1993 which the projectile transfers a much larger momentum either to the emitted electrons or to the molecular recoil ion. Note that unitarity is preserved in our calculations at small impact parameters. Such would not be the case for one-electron calculations that simply add up the contributions from the separate molecular orbitals.    In figure 5 we show the molecular recoil ion charge state fraction of the water ion or its fragments predicted by the Ad-CTMC model. As expected, the single ionization process provides the dominant contribution. Multiple ionization contributions decrease with the degree of ionization in all cases. However, it seems clear that as the projectile charge increases a shoulder-like structure shows up. This trend has been previously observed in experiments regarding the multiple ionization of Ar atoms by highly stripped ions (Muller et al 1986). Argon has similar ionization energies and number of valence electron in its m-shell compared to H 2 O. Results for N 7+ and Fe 15+ ion impact at 1.4 MeV/u are also displayed in figure 5. To trace the physical origin of this structure we inspect again the impact parameter distributions for the three projectiles. As figure 4 suggests, as Z P increases the impact parameter range within which multiple ionization takes place increases. Besides, double ionization, which is clearly dominant over higher multiple ionization mechanisms in the whole b-range for low Z P values, is seen to peak and decrease in intensity for inner impact parameters. As a result, higher ionization degrees dominate the multiple ionization in the small b-region. In terms of recoil charge state fractions, the shoulder-like structure evidences the increasing role of multiple ionization mechanisms with increasing Z P . Although there are no data to directly compare our predictions for the recoil charge distributions of the water target, the Ar data of (Muller et al 1986) displays similar tends to the calculated values to lend credence to our conclusions.
Finally, we address the energy spectra of secondary electrons generated for the collisions systems under study. This is probably the most relevant contribution that the atomic and molecular collisions field can provide towards an accurate physical simulation of a medical irradiation procedure. The needed information in this case is the number of electrons produced per unit length of the ion's trajectory and their energy distributions. Accurate estimations for these quantities are expected to lead to improved numerical models regarding the transport stage and, consequently, a more realistic description of the extension of the dose.
The average energy of the secondary electrons is given by, where W is the secondary electron energy and E 0 is the impact energy. Semi-empirical expressions for proton impact have been previously used to provide an estimate for dσ/dW (Rudd et al 1992, Surdotovich et al 2009, Solov'yov 2017. In these estimates for single ionization the projectile charge enters as a multiplicative factor Z 2 P which for fully stripped projectiles can either represent the nucleus charge or some effective value to indirectly account for charge transfer processes (Barkas 1963). As a result, W ave is independent of the projectile charge and the larger cross section is related to a larger number of electrons emitted. Present results contradict this picture. As it has been stated above, single and multiple electron emission can be associated to different impact parameters and, as a result, to different momentum transfers by the projectile. Low energy electrons are characteristic of large impact parameter collisions involving low momentum transfers, while multiple ionization is provided by ion collisions corresponding to small impact parameters, for which the projectile momentum transfer maximizes, and with such, the chance of energetic electron emission. These statements are corroborated by the single and multiple ionization contributions to the SDCS shown in figures 1 and 3. While the present treatment is expected to be more appropriate at intermediate impact energies, the available experimental data corresponds to the MeV/u regime. In that context, the present classical model predicts a 1/E 0 behavior for the cross sections instead of the quantum mechanical log(E 0 )/E 0 . The difference is due to the classical underestimation of the CTMC ionization probability at low projectile momentum transfers compared to the first born approximation (FBA) (Kaganovich et al 2006). As a result, the low energy electron emission is underestimated as can be inferred from figures 1 and 3(a). To correct this underestimation of the single ioniz ation channel in our SDCS, we propose the following correction term, where N e = 2 is the number of electrons in each molecular orbital. This correction term incorporates the difference between the FBA and the CTMC SDCS considering hydrogenic targets with effective charges for the different molecular orbitals.  The obtained results are shown in figure 6 and are seen to provide an improved description of the experimental data compared to the pure classical predictions, especially for low projectile momentum transfers. As the emitted electron energy (and consequently the projectile momentum transfer) increases, this correction term tends to zero, clearly indicating that the CTMC and FBA tend to be in concordance. We have verified that the correction term of equation (7) fulfils the Z 1.7 P scaling for the SDCS in energy giving support to our previous statements. In table 2 we compare the effect of including this correction at the total ionization cross section σ net level against the reported experimental values and the predictions of the CDW-EIS model in its prior and post forms (Bhattacharjee et al 2017).
The resulting mean secondary electron emission energies W ave for 4 MeV/u C 6+ , and Si 13+ and 3.75 MeV/u O 8+ collisions on H 2 O are shown in table 3. The contributions arising from the corrected single and multiple ionization processes are presented together with those obtained for the resulting total ionization. The experimental W expt ave obtained from the reported data (Bhattacharjee et al 2017) is also included for comparison at the total ionization level. As the projectile charge increases, the W ave for single ionization decreases while that corresponding to multiple ionization increases. In terms of the total ionization, and based on the limited number of projectile charges explored, our results do not suggest a clear trend as a function of Z P . The reported experimental data lead to average energies for the secondary electrons that are lower than those theoretically predicted by either the present CTMC models as well as Rudd's semi-empirical expression. For the latter, a W ave value of about 46 eV for C 6+ projectiles at 4 MeV/u (Rudd et al 1992). It is worth noting that W ave for total ionization can increase by a factor of 2 if events corresponding to E e < 10 eV are not considered. This highly sensitive energetic region requires special care from the experimental point of view, since the use of pre-acceleration voltages to collect slow electrons can result in a lensing effect that increases the acceptance solid angle, and as a result the cross section for low energy electrons. Provided that the determination of accurate values for W ave is probably one of the major contributions that the atomic and molecular collisions community can provide to the ion-therapy program, a precise determination of these experimental uncertainties would allow implementing convolution processes over the available theoretical results. Further work is needed in order to provide a detailed evaluation of this statement.

Conclusions
In this work we have presented a new CTMC adaptive model which explicitly considers the eight electrons of water's valence molecular orbitals. This model has been conceived to compensate for the well-known limitations of the IE and sequential models regarding the energy deposition required in multiple ionization events and the electron distribution in the different molecular orbitals.
Single-differential cross sections in energy and angle have been contrasted against recently reported experimental data. Present results suggest that single and multiple ionization processes lead to distinctive electron energy emission spectra which can be related to different of impact parameters. Consequently, the average energy of secondary electrons corresponding to multiple ionization events is expected to be larger to that corresponding to single ionization. These trends are expected to accentuate for increasing Z P values.
Additional tests of the present model at lower impact energies, in which charge transfer and transfer ionization become relevant are currently under way. These studies are expected to shed light on the role of the different collisional mechanisms in the Bragg peak region.
We hope that the present results will stimulate further work in this area, given their particular relevance for numerical codes that simulate the irradiation procedure in ion therapy. Table 3. Average energies of the secondary electrons (in eV) following 4 MeV/u C 6+ , Si 13+ and 3.75 MeV/u O 8+ collisions with H 2 O. The Expt. values in the table have been obtained by numerically integrating the reported data of Bhattacharjee et al (2017Bhattacharjee et al ( , 2018