Directional dependency of electronic stopping in nickel, projectile’s excited charge state and momentum transfer

We studied the directional dependency of electronic stopping power of swift light ions in nickel using real-time time-dependent density functional theory. We report a variation of electronic stopping for moving ions as the projectile probes different electronic densities of the host material. These results show that while the predicted magnitude stays in reasonable agreement with experiment, for v>2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v > 2$$\end{document}. a.u. simulating only low index crystallographic directions is not enough to sample the experimental average values. The ab initio simulations give us access to microscopic quantities, such as non-adiabatic forces, momentum transfer and transient excited state charges of the projectile and host ions, which are not available through other methods. We report these quantities for the first time.


Introduction
The study of energetic ions traversing condensed matter has important applications in areas of nuclear medicine [1,2], deep ion implantation [3][4][5], defect profiles [6], radiation protection [7] and atomic structure of surfaces [8,9]. The retarding force acting on the energetic ion (projectile) due to the interaction with the irradiated material, which results in a loss of the ion kinetic energy, is referred to as stopping power S. This force can be quantitatively described as the average energy deposited per path unit length by the ion to the host (target) material. The stopping power depends on the type of ion, its kinetic energy (or its velocity relative to the host material) and also on the properties of the material it passes through. In crystalline materials, not all trajectories are equivalent; there is yet an extra dependency on its orientations with respect to crystallographic directions. This dependency is the main topic of this manuscript. Atomistic simulations are in a unique position to explore this dependency in a systematic way, beyond what models [10] and the most accurate state-of-the-art experiments can do [11].
There are two main contributions to the stopping power phenomena: (i) electronic stopping power S e [12][13][14][15][16] which is mainly due to electronic excitations of the target electrons and (ii) nuclear stopping power S n [17][18][19][20] which is due to the elastic collision of the projectile with the target nuclei, where energy is dissipated as lattice defects or vibrations. Our discussion in this paper focuses on the first category of electronic stopping and pertains to the regime in which the phenomena occur a e-mail: correaa@llnl.gov (corresponding author) in the femtosecond time scale, with explicit simulation of electron dynamics.
Extensive experimental work on stopping of light ions on different targets at a wide range of energies has been performed over the years [21][22][23][24][25][26][27][28] (and references therein), obtaining a reasonably complete picture of the average behavior of electronic stopping power in a diversity of materials. Although there are numerous research works on stopping of ions at higher energies, there are little available data on intermediate and low energies regimes and the probing of the effects of crystalline structure of the target materials [29,30].
As in the general case, the electronic stopping power depends on the specific target material, the elemental projectile and its velocity. Apart from this, projectile direction with respect to the crystallographic directions of the target affects the electronic stopping power as well, since the projectile tends to lose energy more effectively as more electrons are available for excitation in the nearby region. For example, a higher electronic density tends to yield higher electronic stopping. This effect is typically averaged out over directions for polycrystalline samples and most experimental setups. In fact, in Ref. [31] we proposed that a random direction can capture this average. At low energies however, an explicit discrimination of low index directions is necessary to improve averaging methods [32].
Among these crystalline systems, nickel and nickelbased alloys are known for their optimal mechanical properties, thermal stability, withstanding to swelling and radiation tolerance, therefore making them good candidates for radiation resistant structural materials for energy applications [33][34][35][36][37] and aerospace applications [38][39][40]. Given the recent interest in this family of materials, our goal is to study S e of proton and of helium ions in a nickel target to describe the initial stages of particle radiation in these crystalline materials. We focus on crystalline effects essentially absent in historical theories until this point [41], i.e., beyond homogeneous electron gas and/or linear response approximations. Recent advances in simulation techniques [42] and computational power [43] enable this study.
White and Mueller [44] measured the electronic stopping for H and He at energies 30-140 keV and 50-120 keV, respectively, in Cr, Mn, Fe, Co and Cu targets, and in Ni target in particular. The measurements were made utilizing a high-resolution silicon surface-barrierdetector system. Their measured average stopping values (with a reported error of ± 3%) for 4 He (Fig. 2 in Ref. [44]) and 1 H (Fig. 3 in Ref. [44]) were in good agreement with previous measurements in Ref. [45].
Gertner et al. [46] investigated the electronic energy loss of ions in solids in the energy range 10 1 -10 4 keV/ nucleon both experimentally and theoretically. Experimentally, they measured the S e on C, Cr, Cu and Ni thin targets, where the energy loss for protons and deuterons was reported. Theoretically, they used the Lindhard-Winther (LW) [47] formalism for the homogeneous electron gas to compute the S e of ions in a freeelectron gas model (with no crystalline information). In the high-energy region (proton energies above 500 keV), they obtained relatively fair results between theory and experiment [46]. In this regime, the main contribution to S e is universal (linear response of the homogeneous electron gas); however, at lower energies, significant discrepancy is observed with other models where the major contribution to S e stems from the valence electrons only, and therefore, their lattice structure, electronic structure and directionality have to be modeled explicitly. Nonlinear response effects, which tend to be relatively stronger at low velocities, and explicit treatment of the target lattice are the main motivations to utilize ab initio simulations like the one presented here. These aforementioned historical results are now included indirectly in the SRIM database and model [10], which we include in the plots for comparison along with more recent experimental results by Roth et al. [48] (for H + ), Mikheev et al. [49] (for He) and Tran et al. [50] (for both H + and He projectiles in Ni).
Regarding directionality studies, Li et al. [51] calculated, from first principles, the electronic stopping power of slow H + and He 2+ ions in CdTe under channeling conditions. For protons, they obtained good agreement between calculated results and experimental results up to the stopping maximum along the 100 and 111 directions. For the case of helium ion, the 100 direction results agreed well with experimental results but in the 111 direction, they observed a transition between two velocities regimes at v = 0.4 a.u. which was attributed to extra energy loss due to charge resonance (in addition to electron-hole excitation). However, the results for 110 direction underestimated experimental data for v > 0.5 a.u. and v > 1.0 a.u. for protons and helium ions, respectively. Recently, Lee and Schleife [52] used time-dependent density functional theory to study the radiation tolerance of phosphide-based III-V semiconductor materials. To understand the high radiation tolerance of these materials, they calculated electronic stopping and proton dynamics in InP, GaP and In 0.5 Ga 0.5 P as an initial step in understanding collision cascades. They obtained a pronounced direction-dependence of electronic stopping along different ( 100 and 110 ) channels. This was attributed to the magnitude of electron density that interacts with the proton projectile in these crystalline channels.
In this paper, we present the electronic stopping power of nickel target under proton and helium irradiation, taking into account the directional dependence of the projectiles on the electronic stopping. The paper is structured as follows: We describe the simulation method, the theoretical framework and computational details in the next section; later, we present the results and their corresponding discussion, where we concentrate on the detailed analysis of the electronic stopping power and other microscopic quantities that can be extracted from postprocessing first principles simulations.

Method
We used time-dependent density functional theory (TDDFT) [53] for this study. TDDFT is one of the most practical frameworks used over the years to tackle the dynamics of many-body electronic systems [16,31,36,52,[54][55][56][57][58]. To obtain the electronic stopping, we solve the time-dependent Kohn-Sham (TDKS) equations. The TDKS equations describing the wavefunctions ϕ i of effective individual electrons are (in atomic units): (1) where the electronic density is given by The Kohn-Sham effective potential V KS is conceptually partitioned in three contributions: is the external potential due to nuclei core potential (including moving ions during the stopping process, hence the explicit tdependence). V H [n(t)](r) is the Hartree potential which describes the classical electrostatic interaction between electrons. V xc [n(t)](r) denotes the exchange-correlation (XC) potential and it is a placeholder for common approximations to the many-body effects of the electron problem. (Time dependence is implicit in n.) All the calculations are performed using the Qball code package [43] which contains time-dependent and other modifications over the Qbox code [59]. In this study, local density approximation (LDA) [60] XC functional was used within the adiabatic approximation. Note that the adiabatic label here refers to the fact that this XC term is taken to be an instantaneous function of the electronic density and not the density values at other (past) times; these simulations are still fully nonadiabatic with respect to ion motion.
For the current framework, a plane-wave pseudopotential scheme is used in solving the time-dependent Kohn-Sham (KS) equations. Electronic wavefunctions ϕ i are expanded in terms of plane waves. Coulombic ionic potentials are replaced by smoother versions with approximately the same scattering properties [61] that allows to use a finite set of plane-wave terms to simulate valence electrons and ignore core, deeply bound electrons.
We use a simulation cell consisting of 108 Ni atoms (1728 explicit valence electrons) in a cubic supercell with a simulation volume of (19.98 a B ) 3 forming an FCC lattice. Periodic boundary conditions along with Ewald-type summation for the Coulomb energy are used throughout this study. The ions are represented by the norm-conserving Troullier-Martins pseudopotentials [62], and the KS orbitals are expanded up to a 160 Ry energy cutoff.
As any time-explicit method, a particular instance simulation needs the specification of an initial condition, both for the ion positions and the electronic state. The initial condition corresponds to a perfect FCC lattice together with the projectile (either H or He) stationary at its initial position with electrons in the corresponding ground state. The electronic ground state (i.e., ϕ i (t = 0)) was obtained by performing a timeindependent density functional theory (DFT) calculation, by solving the time-independent KS equations for the target plus projectile system. Therefore, a timedependent DFT calculation was then carried out starting from the ground state, resulting in a time-evolving electronic wavefunction (ϕ i (t > 0)).
In each simulation, a specific velocity is given to the projectile and the electronic system is monitored for several femtoseconds. This was done to allow the moving projectiles to produce a dynamical steady state where the total energy of the electronic system increases at a well-defined rate as described in Ref. [42]. The remaining (host) ions are fixed at their perfect (unrelaxed) ion positions. This fixed-ions (host) strategy is not only a good approximation given the velocity of the projectile, but it is also a practical constrain to evaluate the electronic-only stopping (without spurious contributions from the nuclear-purely ionic-stopping).
The time propagation of Eq. 1 of the electronic system is achieved by the enforced time-reversal symmetry (ETRS) integration [63], implemented as described in Ref. [43]. A numerical integration time step of 0.024 attoseconds was used. Details of the convergence of the simulation parameters for nickel-based systems are discussed elsewhere [36,64].
The method described here, including the approximations regarding electronic interactions, proved in the past to achieve a good balance of accuracy while treating large complicated atomistic systems. The method may still lack the accuracy of many-body approaches, for example, those based on the dielectric response including the exact dynamic exchange-correlation treatment for low velocity projectiles [65].

Total energy increase of the target system
When the projectile is given a fixed velocity v, the total energy of the electronic system, initially at the ground state value, changes continuously as the projectile moves in the material. The total electronic energy E of the electronic system is given instantaneously by: where the terms are: the electronic kinetic energy, the electron-ion (external) potential V ext (r, t), the Hartree (Coulomb) interaction energy E H [n](r, t), the XC density-functional approximation E xc [n](r, t) and E ion−ion (t) is the interaction between the nuclei (ion). Figure 1 shows the increase in the total electronic energy as a function of the distance traveled by the projectile, resulting from a time-dependent simulation at different velocities in a specific direction.
As the projectile starts to travel in the crystal, a steady state is achieved at some point which is roughly marked by a vertical dotted line; x ≥ 5.0 a 0 . After this transient region, a value for the electronic stopping is then obtained as the slope fitting of the change in total energy versus projectile displacement (Eq. 5).
At finite velocities, the projectile produces electronic excitations, a non-adiabatic behavior that in turn increases the electronic energy. (In the real system, this happens at the expense of the ion kinetic energy; in the simulation, the ion kinetic energy is held constant for simplicity.) The smooth oscillations in the total energy, on top of an otherwise mostly linear increase, are correlated to the lattice periodicity of the target (host) crystal. The correlation is clear in the v = 0 case, in which the simulation should not show an overall increase in the energy (adiabatic behavior), just oscillations in the energy which are natural to the periodic crystal.
For lower velocities, a linear fit coupled with oscillatory function is used to extract the S e with a minimal fitting error. At higher velocities, extracting the slope is more straightforward since there is relatively less oscillation amplitude in the energy versus position curve, so a simple linear fit is enough to extract electronic stop-  Ref. [16]). The slopes are systematically extracted for different velocities and different trajectories for both the H and the He projectiles. This rate of change in E(t) as a function of penetration depth x(t) in the material is identified as the electronic stopping power S e and is given by the average over a trajectory;

Directional effects
In this section, we concentrate on the behavior of S e for both, proton (hydrogen) and helium, in three different trajectories characterized by channeling directions.
Since the projectile explores different environments and electronic densities of the lattice with respect to its direction, we expect the total electronic energy E(t), and therefore the stopping power, to change with each direction.
In Fig. 2, our TDDFT simulation results for protons in bulk Ni are shown for three different projectile directions. For 100 and 110 directions, our calculated S e is below the empirical SRIM data [10] for all range of velocities considered. Ullah et al. [32] obtained a similar trend, an underestimation of experimental results when they calculated electronic stopping of H in Ge for directions with open channels and therefore low electronic density ( 001 and 011 projectile directions for v > 0.3 a.u. in the diamond structure). Our 111 direction result coincides with SRIM empirical data for velocities below 1.5 a.u. and for higher v, but mismatches the empirical results by more than 25% at higher velocities. No particular direction is expected to give the experimental result for experiments without directional information; in previous works, we developed an averaging method based on random directions (random high index) to compare with the bulk of these experiments.
The reason for the discrepancy between particular directions and random directions will be explained when we discuss the effective electronic density in an FCC Ni target along the three directions in the following section.
The calculated electronic stopping of He in nickel is shown in Fig. 2 for the same range of velocities and directions. As in the proton case, the 100 and 110 directions result in values below the experimental data. For the 111 direction, we find it to coincide with experimental data for most of the velocity regime. Although in the 1.5 < v < 3.0 a.u. range of our TDDFT results overshoots over SRIM data, the values still remain within 5%.
For v < 0.5 a.u., the S e results for the three different directions are converging into one curve. This means that as v → 0, the relative difference in magnitude of the electronic stopping decreases for the different directions. In previous simulations (Refs. [16,31,36,54]), a similar pattern is observed, where channeling (in this case 100 ) as well as off-channeling trajectories collapse into one curve at lower velocities, apparently showing independence of geometry on electronic stopping.
We observe that as the projectile velocity increases, there is significant difference in the electronic stopping for the different directions for both types of projectiles. For instance, in the case of helium, the difference between the S e for 100 and 110 directions is up to 10%, and between the 100 and 111 directions it can be a factor two. To understand in detail these direc-tional dependencies on the stopping power as depicted in Fig. 2, we plotted the ground state electronic density for the three directions ( 100 , 110 and 111 ) along the z-axis as shown in Fig. 3.
The 111 direction has the highest average density of 0.098 electrons/a 3 0 compared with the other two directions. This is consistent with the high electronic stopping analyzed in Fig. 2 for both projectiles, since more valence electrons are available to be excited in this path. The average densities yield 0.040 and 0.047 electrons/a 3 0 for 100 and 110 , respectively. This has translated into the relatively lower S e observed for both directions. The relatively small difference in the average electronic density for 100 and 110 directions also corresponds to slight differences in S e .

Channeling and off-channeling trajectories
In Fig. 4, we compare our results for the 111 and offchannel directions with SRIM database and experimental works. For v < 2.0 a.u., the 111 and off-channel simulation results have similar stopping power values with less than 1% difference, as seen in Fig. 4a, which is an indication of fewer geometric effects on the stopping power. The two direction results compare well with recent experimental results from Roth et al. [48] and Tran et al. [50] for projectile velocities less than 1.5 a.u. The maximum S e occurs at similar range in v for both directions with 111 direction having the lowest stopping power value. At higher projectile velocities, the offchannel results show better agreement with the SRIM database, with a gradual decrease in S e for v > 2.8 a.u. and a faster decrease in stopping values for 111 direction at the same velocities.
In Fig. 4b, there is an average of less than 10% difference between the 111 and off-channel simulation results for velocities less than 1.5 a.u. At this lower v, there is a slight overestimation and underestimation of the SRIM results for the 111 and off-channel directions, respectively; also, the 111 results are in good agreement with Tran et al. [50] and Mikheev et al. [49]. For velocities in the range 1.5-3.0 a.u., the offchannel results agree well with SRIM results, except near v = 3.5 a.u., where the 111 results agree better.
Geometric effects on the stopping power are significant for the helium ions even for v < 1.0 a.u., an observation which is absent in the proton case where the two direction ( 111 and off-channel) curves concur into a single curve at v < 1.0 a.u.

Charge fluctuations on projectiles
One of the advantages of using first principles atomistic methods is that it provides a lot of detailed information that can be used to contrast effective models. In particular, by using TDDFT for the stopping problem, the a priori unknown effective charge of the projectile (excited by collisions with the host) is not a simulation parameter like in model methods, but it is a natural outcome of the simulation. The main dynamical quantities in this framework are the Kohn-Sham orbitals (Eq. 1) and the total electronic density, which make the work of finding the (a posteriori) projectile charge a rather challenging task of postprocessing and physical interpretation. The challenge is that atom's 'charge' (or degree of ionization) is not well defined in the general framework and a sensible determination of it has to rely in separating particular bound and free electrons in the vicinity of the ion (projectile or host atom). In place of a complicated determination, we use an approximate density-based analysis method, called Bader analysis [67][68][69][70].
In Fig. 5, we show the electronic charge transfer values for the projectiles' trajectory as a function projectile displacement for v = 0.5 a.u., 1.0 a.u. and 2.0 a.u. along the 100 , 110 and 111 directions. The top panels and the bottom panels depict the charges on proton and helium projectiles, respectively. The charges gained or lost by the projectiles in real time are extracted from the time dependent electronic density, n(r, t), using the Bader Charge Analysis program [66]. The Bader analysis consists in assigning elec-  [10] (black solid lines) and recent experimental results ( [48] (H), [50] (H and He) and [49] (He)). The off-channel represents the average of random directions of the projectile in the nickel lattice and its data points are taken from Ref. [36] tronic charge to each atom by dividing space by surfaces defined by saddle points of the charge density. It is intuitively justified and practical in many molecular systems, although its application is extrapolated to timedependent problems and bulk systems. For example, the assigned charge on a particular ion can be discontinuous in time even if the density changes continuously (saddle points can disappear and appear suddenly); we see instances of this artifact in the analysis of our simulations of the H projectile in narrow channels. The method is practical enough to draw conclusions about relative changes in the nominal charges.
For the proton case, the initial charge value is nominally 1.09 and there are smooth oscillations of the projectile's charge, except in the 111 direction for which even at low velocities the projectile is completely depleted of its electron when the projectile passes through high density regions. High density regions can be identified in Fig. 3. In general, higher velocity implies more ionization regardless of other factors.
For the helium case, we obtained an initial electronic charge of 2.09 that approximately corresponds to the number of electrons on He. The smooth oscillation of  Fig. 5d and e, respectively. In contrast to H, the 111 direction did not show complete electronic charge depletion in Fig. 5f. For both projectiles, we observed alternate neutralization and re-ionization processes along their trajectories [51,71]. These differences in the electronic charge of the projectiles in Fig. 5 for different velocities and direction contribute to the change in electronic stopping observed in Fig. 2. The simulations give ionization cycles that are in the scale of atomic distances in the lattice, rather than longer distances related to purely atomic (to the projectile) processes such as Auger processes, although the lack of self-interaction corrections could be biasing this particular aspect of the results. (See Ref. [58] for a study of the relation between charge oscillation in a projectile and self-interaction effects in binary atomic collisions within the same simulation framework.)

Charge deficiency on host atoms in the 100 direction
In the same way as for the projectile, the charge analysis can be applied to the (stationary) host atoms (Fig. 6). Persistent charges in the host atoms lead to impulsive Coulomb-like forces that can impart radial momentum [57]. The charge deficiency (relative to its initial value) on a particular neighboring Ni atom due to the incoming projectiles (H and He) is depicted in Fig. 6. We choose an atom that is first neighbor to the trajectory channel. Initially, the projectiles are ∼ 3.33 a 0 from the Ni atom shown in Fig. 6a. Using the Bader analysis, the charges on the Ni atom are again extracted from the time-dependent electronic charge density, except that this time the ion center has a fixed position. As the proton and helium get closer to the neighboring Ni atom, the charges decrease gradually, indicating that the particular atom is losing charge as the projectile approaches. As the projectile passes the particular host Ni atom (after x > 8 a 0 ), its charge gradually stabilizes due to less interaction with the projectile. Interestingly, the cases above v > 1a 0 do not recover to the initial value for the duration of the simulation. Although outside the practical time-scales of this simulation framework, a persistent charge could lead to large momentum transfer, as discussed later.

Ratio of helium/hydrogen electronic stopping
To understand our S e results, we compared with linear response theory (LRT) [47]. We evaluate linear response stopping S L based on the RPA Lindhard dielectric function for a homogeneous electron gas for effective electronic density n [72] as: where Z 1 is interpreted as the effective charge of the projectile, v is the velocity of the projectile and ε RPA is a model of the dielectric function for the linear response approximation (RPA) at frequency ω and momentum transfer k. Regardless of the dielectric model, a characteristic of linear response is that S L has a quadratic dependence on the projectile ion charge S L ∝ Z 2 1 , that is, an exact ratio of 4 for fully ionized hydrogen and helium. Therefore, deviations from this factor could be used as proxy of nonlinear effects in nickel (including bound charge effects).
The ratio of S e for helium ion to the S e for proton are depicted in Fig. 7 as a function of the projectile velocity, and the results are compared with linear response theory (LRT) [47] for full ionization and SRIM [10]. The TDDFT ratio surpasses 4 for v ≥ 2.5 a.u, which is beyond the maximum stopping for H. Their ratio should eventually approach 4 again for even higher velocities. Yost et al. [56] observed a similar trend when they utilized the same theoretical approach to study protons and α particles stopping in silicon carbide.
This is an indication that an ideal factor for 4 is achieved as both the projectile becomes fully stripped and the linear response approximation becomes more accurate at large velocities. At lower velocities, for LRT to correctly predict the maximum stopping, effective ion charge models and additional higher-order Z need to be incorporated [56,73,74]. Figure 8 shows the non-adiabatic radial forces exerted on a neighboring target atom by projectiles (H + and He) for a channeling trajectory, respectively. Initially, the distance of the closest Ni target atom to the projectile is about ∼ 3.33 a 0 . For zero velocity (adiabatic limit), a maximum force on a host atom is obtained exactly when the projectile is closest. The magnitude of this adiabatic force is higher for a He projectile. This maximum value of the force is due to electronic excitation [16,31,54] at the closest proximity of the projectile to the target atom. As the projectile travels further away, the force decays smoothly to zero (its original value) for low velocities. (Note that the non-adiabatic curves have been shifted for clarity.)

Radial force
At higher velocities v > 0.7 a.u., there is persistent oscillation in the forces' curves, as the projectile passes (a) (b) Fig. 8 Radial force exerted on host atom (Ni atom closest to the projectile's channel 100 trajectory) for a H + and for b He, versus projectile position for different velocities. The non-adiabatic curves have been shifted upwards for clearer view. Notice how for v > 0.7 a.u. the maximum shifts relative to the position of maximum approach and, instead of going back to zero, the forces remain oscillating even for large distances Fig. 9 Radial momentum transferred to host atom (Niatom closest to projectile channel 100 trajectory) vs. projectile velocity and the position of the maximum force with respect to their low velocity (near-adiabatic) counterparts shifts. These oscillations can be attributed to persistent plasmonic charge oscillations in the electronic system [54].

Radial momentum
Although eventually the radial force returns to its original value, the integral of the radial force in time does not vanish and depends on the velocity. A net force acting during a period of time produces a momentum change. For forces that depend only on the position (adiabatic forces), the momentum transfer simply decays inversely proportional to the velocity of the perturbation (as the interaction time shortens). This is the expected behavior for any classical potential for projectile-ion interaction, in particular any binary collision model [75] but deviations from it correspond here to increasing electronic non-adiabatic effects. Figure 9 shows the lateral momentum transfer that the host atom incurs due to the non-adiabatic forces, calculated Δp = F[x(t)] dt = F(x)dx/v, by the projectile to the target atoms. This is the inverse of the projectile velocity (1/v) multiplied by the time integral of the force on the host atom over some specified projectile distance (Fig. 8). At lower velocities, the adiabatic approximation gives a good estimation of the momentum transfer up to v = 0.2 a.u. and v = 0.15 a.u. for proton and alpha projectiles, respectively. At this adiabatic limit, the momentum transfer is proportional to the function k/v, where k is the integral over the adiabatic force curve as a function of projectile displacement (see Fig. 8a, b, for v = 0). For higher velocities, (v > 0.2 a.u.), there is a complete deviation from their adiabatic counterparts. This is the point from which considering electronic effects would affect the results of models based purely on the molecular dynamics and interatomic potentials (in atomic time-scales).
The momentum transfers are surprisingly proportional (parallel curves in the log scale) for the proton and the He projectiles at different velocities; this proportionality extends to the non-adiabatic regime of momentum transfer. This points to a possible universal behavior to the momentum imparted by a projectile of arbitrary charge and across host materials (compare to Fig. 4 in Ref. [54]).

Conclusion
In summary, we presented first-principles calculations of electronic stopping power in nickel crystals for protons and helium particles in three different channel directions using a real-time electron dynamics formalism based on TDDFT. Our TDDFT results are comparable with SRIM database for most of the velocity range considered for the 111 and random directions.
We have learned that the electronic stopping is sensitive to projectile direction with higher density regions yielding the highest S e value due to stronger electronion interaction. The lower density regions (most open channels) showed the lowest stopping power.
We have observed a rapid charge oscillation on the projectiles at narrower channels with a complete pro-jectile ionization (depletion of electrons) for hydrogen in the 111 direction at regions during its trajectory in the crystal. Complete charge depletion for He is not achieved in any channel at any velocity studied, which is attributed to the enhanced binding of the remaining electron in the single ionized state of the projectile.
The non-adiabatic lateral momentum transfer turns to follow an approximately proportionality pattern between protons and alpha particles, both in the adiabatic limit and in the non-adiabatic limit. This scaling could help develop phenomenological models of nonadiabatic forces of swift ions in solids [76]. Very similar momentum transfer curve shapes were found for other materials, hinting to a possible universal behavior [54].
While absolute values of stopping power can be obtained by many methods, such as ab initio (this work), semianalytic methods [77] or experiments [78], it is clear that ab initio methods can access microscopic quantities such as non-adiabatic forces, specific ion charges and momentum transfer that are otherwise out of the scope of other methods.
Acknowledgements This work was performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. The work was supported by the US Department of Energy, Office of Science, Materials Sciences and Engineering Division. Computing support for this work came from the Lawrence Livermore National Laboratory Institutional Computing Grand Challenge program.

Author contributions
E.Q. carried out all the simulations and wrote most of the paper. X.A. helped running the simulations and revised the manuscript. A.A.C. conceived the directional and charge analysis of the results, guided the simulations and contributed to writing the paper.

Data Availibility Statement
This manuscript has no associated data or the data will not be deposited [Authors' comment: All the results presented in the paper can be reproduced with the information provided in the manuscript. The authors can be contacted for further information.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.