Electronic stopping and proton dynamics in InP, GaP, and In$_{0.5}$Ga$_{0.5}$P from first principles
Cheng-Wei Lee André Schleife





Department of Materials Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA


schleife@illinois.edu
Department of Materials Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
Frederick Seitz Materials Research Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA




The phosphide-based III-V semiconductors InP, GaP, and In_0.5Ga_0.5P are promising materials for solar panels in outer space and radioisotope batteries, for which lifetime is a major issue.
In order to understand high radiation tolerance of these materials and improve it further, it is necessary to describe the early stages of radiation damage on fast time and short length scales.
In particular, the influence of atomic ordering, as observed e.g. in In_0.5Ga_0.5P, on electronic stopping is unknown.
We use real-time time-dependent density functional theory and the adiabatic local density approximation to simulate electronic stopping of protons in InP, GaP, and the CuAu-I ordered phase of In_0.5Ga_0.5P across a large kinetic energy range.
These results are compared to SRIM and we investigate the dependence on the channel of the projectile through the target.
We show that stopping can be enhanced or reduced in In_0.5Ga_0.5P and explain this using the electron-density distribution.
By comparing Ehrenfest and Born-Oppenheimer molecular dynamics, we illustrate the intricate dynamics of a proton on a channeling trajectory.


Electronic stopping and proton dynamics in InP, GaP, and In_0.5Ga_0.5P from first principles
    André Schleife
    December 30, 2023
============================================================================================





§ INTRODUCTION

Indium phosphide (InP) and In_0.5Ga_0.5P are well-suited materials for optoelectronic devices due to their direct (low-temperature) band gaps of 1.42 eV <cit.> and 1.99 eV <cit.>, respectively.
Gallium phosphide (GaP) has an indirect gap of 2.34 eV at low temperature <cit.>.
Using In_0.5Ga_0.5P, a tandem solar cell was demonstrated with an efficiency greater than 30 % in a double-junction <cit.> and over 40 % in a triple-junction <cit.> configuration.
In addition, In_0.5Ga_0.5P shows good resistance to energetic, charged-particle radiation, making it suitable for applications in extreme operational environments where lifetime is one of the major issues.
Examples include solar panels in outer space <cit.> and radioisotope batteries <cit.>.

Research devoted to analyzing degradation of solar panels caused by charged-particle radiation, typically relies on semi-classical models <cit.> derived from the Shockley-Read <cit.> and Hall<cit.> equation to describe recombination of electrons and holes in semiconductor devices.
This allowed attributing a gradual drop in efficiency of solar panels as fluence of radiation increases to decreased minority-carrier life times.<cit.>
In addition, radiation-induced defects in InP based solar devices were found to be annealed by injection of minority carriers<cit.> and the performance was partially recovered.
The enhanced annealing was attributed to the Bourgoin mechanism <cit.>, i.e., a change of the charge state of defects due to injection that leads to faster diffusion. 
These insights illustrate that the semi-classical approach is useful for optimizing the design of devices, however, it has no access to atomic-scale details of the interaction between the charged projectile ions and the target material.
Such details are essential for understanding the underlying atomistic mechanisms.
Achieving this goal requires modern first-principles simulations such as the ones described here.

Previous studies <cit.> showed that the defect dynamics in target materials exposed to charged-particle radiation differ between regions of bulk and interfaces, since interfaces can act as sink or source of defects.
Itoh reviewed the effect of interfaces on defect dynamics under the scenario of projectile kinetic energies that are too low to induce knock-on events <cit.>.
It was speculated that in this scenario, enhanced damage near interfaces can be attributed to stronger localization of excitons or slower recombination rates for Frenkel pairs <cit.>.
Furthermore, a recent study based on ab-initio molecular dynamics <cit.> for primary knock-on events under particle radiation shows that cations in the GaAs/AlAs superlattice are more likely to be displaced than cations in pure GaAs or AlAs.
Therefore, it is critical to model the effect of interfaces on radiation damage.

Existing first-principles studies that aim at unraveling the effect of interfaces are limited to the linear-response approximation and focus on optical properties<cit.> instead of electronic response to radiation.
Gumbs proposed an analytic expression for electronic stopping of a charged particle moving parallel to the surface of layered 2D free-electron gases, based on the random-phase approximation <cit.>.
However, this approach is limited by the linear-response approximation and the specific geometric setup used in the derivation.
In particular, the charged projectile moves outside of the heterostructure and parallel to the surface.
Recently, Cruz combined the Bethe stopping theory<cit.> with a model of quantum confinement that imposes boundary conditions on the system, to study the effect of interfaces on electronic stopping <cit.>.
Although this method is not limited to a specific geometric setup, it still suffers from the linear-response approximation and the assumption of a fully ionized projectile as well as quantum confinement.

For device applications, high-quality In_0.5Ga_0.5P is fabricated, using molecular-beam epitaxy or organo-metallic vapor-phase deposition.
This leads to well-defined, atomically ordered phases <cit.>, instead of random solid solutions, with the “CuAu-I” ordered phase<cit.> being one simple example.
These ordered phases have different electronic and phonon band structures compared to solid solutions and to bulk materials, giving rise to different optical, electronic, and thermal properties <cit.>.
As discussed above, there are a few studies exploring materials response to particle radiation for interfaces in heterostructures where the components are much thicker than monolayers that are observed in atomically ordered phases.
However, due to the different geometry, these existing approaches cannot be applied to ordered phases irradiated by fast ions.
To the best of our knowledge, there is no literature on how ordered phases with periodicities on the single-monolayer scale affect the ultrafast electronic response to particle radiation.
This is the focus of the present work.

Here we use real-time time-dependent density functional theory (RT-TDDFT) to study the electronic response of InP, GaP, and In_0.5Ga_0.5P to highly energetic protons.
We compute the electronic stopping power and dynamics of the proton projectile for the individual materials.
Our results indicate that interfaces in In_0.5Ga_0.5P give rise to both local enhancement as well as reduction of instantaneous stopping, compared to pure InP or GaP.
We attribute this behavior to the redistribution of electron density caused by the formation of the ordered phase.
In addition, we compare the dynamics of the proton projectile using Ehrenfest and Born-Oppenheimer molecular dynamics.
Their difference suggests the importance of including non-adiabatic and excited-electron effects.

In Sec. <ref> we summarize our computational approaches for ground-state calculations, real-time electron dynamics, and both average as well as instantaneous electronic stopping power.
In Sec. <ref> and <ref>, we report our results for average and instantaneous electronic stopping, respectively, for proton-irradiated InP, GaP, and In_0.5Ga_0.5P. 
In Sec. <ref>, we report the dynamics of a proton moving on a [100] channel using both Ehrenfest and Born-Oppenheimer molecular dynamics.
We compare the difference and discuss the importance to explicitly model electron dynamics. 
Lastly, we conclude and summarize our work in Sec. <ref>.



§ COMPUTATIONAL APPROACH




 §.§ Ground-state calculations




Using the Qb@ll code <cit.>, we performed ground-state density functional theory (DFT) <cit.> calculations for zinc-blende (zb) InP, zb-GaP, and the zb-based ordered CuAu-I phase <cit.> of In_0.5Ga_0.5P (see Fig. <ref>).
On GaAs(001) substrates, the “CuPt” type atomic ordering of InGaP is more commonly observed <cit.>, but the CuAu-I ordering was reported on GaAs(110) substrates before <cit.>.
Hence, even though the CuAu-I phase is not the most common atomic ordering of In_0.5Ga_0.5P, it is chosen here as a reasonable and computationally feasible test case.
Kohn-Sham (KS) wave functions are expanded into a plane-wave basis with cutoff energies of 50 hartree (E_H), 75 E_H, and 75 E_H for InP, GaP, and In_0.5Ga_0.5P, respectively, to obtain total energies converged to within 0.184 mE_H/atom.
The local-density approximation (LDA) is used to describe exchange and correlation <cit.> and the electron-ion interaction is described by norm-conserving Hamann, Schlüter, and Chiang pseudopotentials as modified by Vanderbilt <cit.>.
We use pseudopotentials with 4s^23d^104p^1, 5s^24d^105p^1, and
3s^23p^3 valence electrons for Ga, In, and P respectively.
The Brillouin zone is sampled using only the Γ point, which is justified for the 216-atom supercells  used here.

Relaxed atomic geometries are computed using fits to the Murnaghan equation of state <cit.>.
This yields lattice constants of 11.07 and 10.24 a_B for InP and GaP respectively.
For In_0.5Ga_0.5P, we first determine the a/c ration that gives similar pressure on all faces of the cell, and then scale the cell volume until the external pressure is below 0.5 GPa.
This yields cell dimensions a, b, and c of 10.71 a_B, 10.65 a_B, and 10.65 a_B, respectively.
All atomic positions are relaxed until forces are below 0.1 mE_H/a_B.

In order to isolate the effect of electronic excitations on ion dynamics, we also performed Born-Oppenheimer molecular dynamics (BOMD) simulations <cit.>.
Since the protons that represent particle radiation move very fast, smaller time steps compared to typical BOMD simulations were chosen.
This guarantees enough sampling points (210 points along the [100] trajectory) and conservation of energy. 
More specifically, a time step of 0.3 atomic units (at. u.) of time, 0.1 at. u., and 0.0375 at. u. is used for proton velocities of 0.5 at. u., 1.5 at. u., and 4.0 at. u., respectively.



 §.§ Real-time electron-ion dynamics


We study real-time electron-ion dynamics using the Ehrenfest molecular dynamics approach <cit.>.
Such simulations have become increasingly feasible even for solids <cit.>, both due to the commendable balance of accuracy and computational efficiency of TDDFT <cit.>, and due to the advent of modern supercomputers.
The electronic system is described by propagating time-dependent KS equations in real time using a fourth-order Runge-Kutta integrator <cit.>.
A time step of 0.0145 at. u. was used and we verified that the electronic stopping power extracted from these simulations changed by less than 0.02 % when the time step is halved.

Non-adiabatic electron-ion coupling is described by computing Hellman-Feynman forces from the time-dependent electron density <cit.>.
These simulations are carried out using the TDDFT implementation within the Qb@ll code <cit.>.



 §.§ Electronic stopping power


When charged particles travel through a target compound, they transfer kinetic energy to that material <cit.>.
The energy loss (dE) per penetration depth (dx) is known as stopping power S and has the unit of a force,

    S(x) = dE(x)/dx.

As indicated in Eq. (<ref>), stopping power is the instantaneous rate of energy transfer, e.g. from protons to the III-P compounds studied here.
In the low-kinetic energy regime, the projectile predominantly transfers energy to the ions of the target material (“nuclear stopping”).
However, for protons with kinetic energies higher than about 1 keV, more than 10 times as much energy is transferred from proton kinetic energy to the electronic system of the III-P target material than to the ions (“electronic stopping”).
This electronic-stopping regime is the focus of this work. 

In Fig. <ref> we compare electronic stopping for channeling, i.e. protons that travel on trajectories centered at [100] and [110] lattice channels, to off-channeling stopping geometries.
Our studies of off-channeling trajectories are motivated by experiment and enable us to compare to either amorphous or polycrystalline samples commonly used in practice.
Furthermore, even when the sample is a single crystal, experiment oftentimes studies off-channeling trajectories because standard Monte Carlo packages, such as SRIM <cit.>, fail to predict damage and distribution of defects in target materials under channeling conditions <cit.>.

In this work we follow the approach of simulating a random trajectory through the crystal, as devised in Ref. Schleife2015, to represent off-channeling protons.
For each velocity (projectile kinetic energy) we use a standard pseudo-random number generator to generate a random direction through the lattice.
In order to obtain results that are independent of the specific random direction, we ensure they are dissimilar from any lattice channel and each trajectory is simulated long enough to obtain convergence (see below).
We then fix the velocities of all atoms in the simulation, including the projectile, to exclude primary knock-on events <cit.>.
This also avoids numerical issues caused by very short distances between projectile and target atoms, for which large values of the Coulomb interaction would require much shorter time steps.
While this constitutes an approximation, it can be justified since the cross section for scattering between projectile and lattice atoms is very small for fast, light projectiles.
As discussed in detail in Ref. Schleife2015, this assumption of a frozen lattice is valid for high proton velocities such as the ones studied in this work, for which the time scale of interaction with the lattice is short.
This allows us to use the total-energy increase to compute electronic stopping for off-channeling protons <cit.>.
Full Ehrenfest dynamics simulations, where all ions are allowed to move according to Hellman-Feynman forces, are performed for channeling trajectories.

We compute averages of instantaneous electronic stopping for channeling projectiles by integrating over 2 lattice periods (unshaded area in Figs. <ref> and <ref>) after discarding the first half lattice period of a simulation, to avoid onset effects. 
Along the [100] and [110] directions, the 216-atom cell has three lattice periods but the length of the lattice period in [110] direction is by a factor of √(2) larger than that in the [100] direction. 
Onset effects are obvious, e.g. in Fig. <ref>(a), where stopping near the onset is much larger than at later stages of the simulation.
Discarding also the last half lattice period of the simulation, allows us to mitigate the impact of excited electrons that re-enter the simulation cell due to periodic boundary conditions<cit.>.

As discussed in Ref. Schleife2015, the average electronic stopping for off-channeling projectiles is calculated from the instantaneous value using the slope of a linear regression fit to the E(x) curve.
Initially, this result is sensitive to the trajectory length, however,
it eventually converges when the trajectory is long enough (approximately for a trajectory length of 200 a_B). 



§ RESULTS AND DISCUSSION




 §.§ Average electronic stopping




In Fig. <ref> we show the dependence of the electronic stopping power on proton kinetic energy as computed from our first-principles simulations.
This figure compares two channeling proton trajectories to the off-channeling configuration for GaP, InP, and In_0.5Ga_0.5P.
From this comparison, it becomes immediately clear that electronic stopping in all three III-P compounds depends strongly on the trajectory:
For all proton kinetic energies, the [110] channel leads to the smallest electronic stopping.
The [100] channel shows similar electronic stopping as the off-channeling trajectory before the stopping maximum, but also leads to smaller stopping close to and even more so after the peak of the curve.

The first observation of smaller stopping along the [110] channel can be explained by the effective electron density that the projectile interacts with along this trajectory.
When protons travel on a [110] channel, the average distance between the proton and first-nearest-lattice atoms is about 50 % longer than for protons on a [100] channel.
Since most of the electron density is located near the ions, protons on [110] channels interact with smaller electron density, leading to weaker electronic stopping.
This finding is consistent with a previous RT-TDDFT study of proton-irradiated germanium <cit.> and a study based on scattering theory for energy loss in a non-uniform electron gas <cit.>.

The second observation, that off-channeling projectiles lead to higher stopping than channeling projectiles, has been reported in the literature before and was attributed to stopping contributions from semi-core electrons <cit.>.
In order for semi-core electrons to contribute to electronic stopping, protons need to have high enough kinetic energy to excite the semi-core states.
In addition, these excitations require spatial proximity of the proton and semi-core wave functions, i.e. very small distances between proton and ions, which we only capture by random trajectories in our simulations.

Finally, comparison of our results to data that we computed using “The Stopping and Range of Ions in Matter” (SRIM) <cit.> shows good overall agreement and confirms our interpretation.
Since SRIM assumes an amorphous structure of the target material, the large range of electron density values that a projectile experiences as it traverses an amorphous target is most closely represented by our off-channeling trajectory.
Consequently, when comparing our results for off-channeling electronic stopping to SRIM, we find good agreement before the electronic-stopping peak, but deviations become significant especially for higher kinetic energies.
This behavior has been identified in the literature<cit.> before and one possible explanation invokes electronic-stopping contributions due to semi-core electrons that are missing in the pseudopotentials used here (see supplemental material).
Another limitation is the use of the adiabatic LDA in this work, and, while this a topic of ongoing research <cit.>, it is currently unknown how this quantitatively affects electronic stopping of protons.
Finally, the simulation cell size also affects the accuracy of plasma excitations since it limits the maximum wave length for a plasmon in the simulation <cit.>.



We also note that our results agree with SRIM regarding the relative magnitude of electronic stopping across the different materials.
Except for off-channeling projectiles with v=0.9 at. u.we consistently find stopping in GaP to be the largest,  in InP to be the smallest, and in In_0.5Ga_0.5P to be in between.
More specifically, we find that electronic stopping of In_0.5Ga_0.5P is very close to the average of stopping in InP and GaP.
The data in Table <ref> illustrates that the relative differences are below 1 % across most of the velocities for the [100] and [110] channels.
This also holds for the density of valence electrons (see Sec. <ref>) for these compounds:
That of In_0.5Ga_0.5P is 4.00 × 10^23 cm^-3, which is within 1.3 % of the average of 4.52 × 10^23 cm^-3 (GaP) and 3.58 × 10^23 cm^-3 (InP).
we assume to zeroth order that a proton moving on a channel through In_0.5Ga_0.5P interacts half of the time with InP-like electron density and half of the time with GaP-like electron density.
Within the Lindhard model, electronic stopping is proportional to the electron density <cit.>, hence, we conclude that this model and the linear approximation describe electronic stopping for channeling in the CuAu-I ordered phase of In_0.5Ga_0.5P very well.
We will refine this picture below, using the actual electron-density distribution in In_0.5Ga_0.5P.

As described above, for off-channeling projectiles we use different random trajectories for the different velocities and III-P compounds.
Due to the statistical nature of this approach, convergence is computationally challenging:
While each trajectory converges to good accuracy after about 200 a_B, a given random trajectory may represent a good cell average only after much longer lengths.
We observe this for v=0.9 at. u., where Fig. <ref> shows a different relative ordering for the different materials, compared to the other velocities.
The InP trajectory in this case more often samples close proximity to semi-core electrons and, thus, higher stopping (see supplemental material).
Much longer runs would be required to eliminate this influence from the final stopping result.

From previous electronic-structure calculations<cit.> it is expected that formation of an ordered phase results in breaking of translational symmetry and, therefore, splitting and energy shifting of bands and states inside the band gap.
However, our results indicate, that after averaging over instantaneous stopping, the local electronic structure of In_0.5Ga_0.5P has a very minor influence on electronic stopping.
We attribute this to the large projectile velocities studied in this work, compared to the changes in the electronic structure.
The situation is different for instantaneous stopping, which we discuss next.



 §.§ Instantaneous electronic stopping






Our RT-TDDFT results unambiguously show that instantaneous electronic stopping reveals a dependency on the local environment.
Since all the III-P compounds have slightly different cell parameters, we use the normalized cell length for InP, GaP, and In_0.5Ga_0.5P in order to help visualization and comparison (see Figs. <ref>, <ref>, and <ref>).
This ensures that the same local environment is compared for all the III-P compounds.
Figure <ref> illustrates that instantaneous electronic stopping of protons moving with three different velocities on a [110] channel in In_0.5Ga_0.5P oscillates between InP-like and GaP-like behavior.
When the proton is near the InP layer of In_0.5Ga_0.5P, it locally follows the curve of InP and, similarly, when it is near the GaP layer it follows GaP stopping.
After averaging instantaneous stopping along the trajectory as discussed above, we then find that average stopping in In_0.5Ga_0.5P is very close to the average of GaP and InP electronic stopping.

For protons on a [100] channel, however, we find a totally different behavior and even a velocity dependence, as illustrated in Fig. <ref>.
As can be seen from this figure, for velocities of 0.5 at. u. and 4.0 at. u., the instantaneous stopping of In_0.5Ga_0.5P is locally larger or smaller than that of InP and GaP.
For these two velocities, the ordered phase of In_0.5Ga_0.5P gives rise to local enhancement and reduction of electronic stopping.
However, in the case of a proton with a velocity of 1.5 at. u. the stopping is again within the boundaries defined by InP and GaP, similar to what we discussed above for the [110] channel.
We attribute this velocity dependence to electronic states that appear in the ordered In_0.5Ga_0.5P phase and that lead to the observed behavior.



The ground-state electron density allows to analyze this in more detail and we find that its spatial distribution in In_0.5Ga_0.5P contributes to the local enhancement and reduction.
To illustrate this, Fig. <ref> shows the difference of the electron-density average along the [001] direction between In_0.5Ga_0.5P and GaP as well as InP as a 2D plot.
The top panel shows that in In_0.5Ga_0.5P there is less charge around Ga ions than in GaP, and the bottom panel shows that there is more charge around In atoms in In_0.5Ga_0.5P, compared to InP.
The difference for P atoms is negligible.
Comparing this to the data in Fig. <ref> illustrates that enhanced stopping occurs near In atoms and reduced stopping is observed near Ga atoms for a proton on a [100] channel, which matches the behavior of the electron density near these atoms.
Contrary, the proton on a [110] channel is further away from these atoms and does not sample these electron-density differences.
Hence, no local enhancement or reduction of electronic stopping is observed in Fig. <ref>.
Our observation that not all proton velocities lead to enhancement or reduction of electronic stopping cannot be understood in this model.
Instead, we conjecture that this is related to the specific electronics states in In_0.5Ga_0.5P that are responsible also for the electron-density differences discussed above.



 §.§ Dynamics of a proton on a [100] channel in In_0.5Ga_0.5P




In the following, we provide deeper insight into the intricate dynamics of a proton on a [100] channel in In_0.5Ga_0.5P.
In particular, we disentangle the influence of electronic excitations on the dynamics by comparing Ehrenfest to BOMD.
To this end, Fig. <ref> shows both the forces and the resulting displacement of the proton as it travels through the material for three different velocities.
The force acting on the projectile in different locations in the material in BOMD simulations does not depend on the projectile velocity, as confirmed by the solid lines in Fig. <ref>.
This changes in Ehrenfest dynamics, and in the following we discuss the three different components of that force (see  Fig. <ref> for definition of a, b, and c).

Most notably, the force component parallel to the a direction differs strongly between Ehrenfest and BOMD simulations.
This difference is completely expected and corresponds to electronic stopping, as discussed above.
As such, it is entirely attributed to non-adiabatic excitations that are captured by Ehrenfest dynamics, but not by BOMD simulations, for which the oscillations around zero force integrate to zero.

As shown in Fig. <ref>, we also find non-zero forces for the b and c direction, but only non-zero displacement for b direction within BOMD simulations.
The initial position of the proton at the center of the channel is not the equilibrium position in b direction since In_0.5Ga_0.5P breaks the symmetry along the b axis of InP and GaP.
As the proton moves through the material, it interacts with first-nearest-lattice atoms that repeat in the order In, P, Ga, and P in the directions b, c, -b, and -c, respectively.
It experiences repulsion from all of these atoms, but only the repulsion from P is oscillatory around zero.
In b direction, the repulsion from In is larger than that from Ga, resulting in the displacement shown in Fig. <ref>.

Furthermore, Fig. <ref> illustrates that these force components acting on the proton parallel to the b and c directions become significantly larger and depend on the proton velocity within Ehrenfest dynamics.
While the overall shape of the force parallel to b still strongly resembles the BOMD force, it becomes slightly more asymmetric, leading to larger displacements of the proton along this direction (see Fig. <ref>).
More importantly, the force along c significantly deviates from the BOMD force, both qualitatively and quantitatively, and even shows a different frequency of the oscillatory behavior.
Since these oscillations are not entirely symmetric around zero force, this leads to velocity-dependent displacements of the proton along c that are absent in BOMD simulations (see Fig. <ref>, in particular for v=1.5 at. u.).

Limited by the computational cost of Ehrenfest dynamics, we only report a trajectory of about 30 a_B.
However, even for this short trajectory we clearly identify an effect of electronic excitations on the trajectory of the proton projectile.
While BOMD predicts deviations from an ideal trajectory along the center of the [100] lattice channel in In_0.5Ga_0.5P, this is amplified and becomes velocity-dependent in Ehrenfest dynamics, due to the excitation of electrons.
Velocity-dependent non-adiabatic forces caused by electronic excitations were identified before using RT-TDDFT <cit.>. 

Oscillations of projectiles have been reported before for channeling, however, the effect of electronic excitations is generally neglected <cit.>.
In this work, we accurately quantify this effect and while we find that the magnitude is small, our first-principles results provide the first direct quantitative evidence of an electronic contribution to such oscillations.
In particular, we show that these excited-electron contributions cause non-zero forces even for cases where BOMD finds zero force and, thus, significantly affect the dynamics of fast protons as it moves through the material.
More computational work and, ideally, longer Ehrenfest trajectories are necessary to further study this behavior.



§ CONCLUSIONS


We reported on RT-TDDFT first-principles simulations to investigate  electronic stopping of protons in InP, GaP, and the CuAu-I ordered phase of In_0.5Ga_0.5P.
We compare our results from this parameter-free approach to data that we obtained using SRIM and find very good agreement for proton kinetic energies below about 25 keV.
The agreement is worse for higher kinetic energies, potentially due to core electronic states that were not treated as valence electrons in our pseudopotential implementation.
Nevertheless, we find a pronounced direction-dependence of electronic stopping along different channels and explain this using the magnitude of the electron density the proton projectile interacts with.
We also find a clear indication of local enhancement and reduction of stopping for the [100] channel, and explain this by local enhancement and reduction of the ground-state electron density.
The dependence of this effect on the proton velocity underscores its non-adiabatic character.

While these effects will be difficult to observe directly in experiment, we conjecture that they significantly contribute to the dynamics of charged ions in semiconductor materials.
To investigate this further, we directly study the dynamics of a proton moving through In_0.5Ga_0.5P, using Ehrenfest and BOMD.
This comparison reveals an influence of electronic excitations both on force and displacement of the proton.
Even though the trajectories reported here are very short, they nevertheless illustrate that excited electronic states can trigger dynamics that is absent in a solid in its ground state.
We believe that these effects contribute to oscillations of charged projectiles as they move through a material.
Excited electronic states need to be taken into account in order to understand radiation damage on an atomistic level, and the use of TDDFT in an Ehrenfest MD scheme is a particularly appealing approach to do so, striking a desirable balance between accuracy and computational cost.


C.-W. L. acknowledges support from the Government Scholarship to Study Abroad from the Taiwan Ministry of Education.
Financial support from Sandia National Laboratories through the Sandia-UIUC collaboration is gratefully acknowledged (SNL grant no. 1736375).
An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program.
This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. Data used in this work are available at the Materials Data Facility at http://dx.doi.org/10.18126/M2R337.




§ AUTHOR CONTRIBUTION STATEMENT


A. S. conceived the original idea and supervised the project. A. S. and C.-W. L. planned the simulations and C.-W. L. carried out all simulations. A. S. and C.-W. L. analyzed and discussed the results and both contributed to the manuscript.







§ SUPPLEMENTAL MATERIAL




 §.§ Estimation of threshold velocity for excitation due to fast charged particle


As discussed in Ref. Lim:2016, the threshold velocity, below which no electronic stopping is allowed, can be estimated by Planck's constant (h), distance between equivalent lattice position (λ), and band gap (Δ),

    v_th=λΔ/h.

We extend Eq. (<ref>) to estimate the threshold velocity to excite electrons from each shell to the conduction band minimum by replacing the band gap with the corresponding energy difference, calculated by subtracting electron affinity from ionization energy. 
The distance between equivalent lattice positions is 1/2 of a lattice period, i.e., 1.55, 1.43, and 1.50 a_B for InP, GaP, and In_0.5Ga_0.5P, respectively.
Since this is only an estimate, 1.50 a_B is used for all the calculations.
The electron affinity for InP, GaP, and In_0.5Ga_0.5P is 0.16, 0.14, and 0.15 E_H, respectively, and 0.15 E_H is used for all the calculations.
Ionization energy and threshold velocity (kinetic energy) for each shell are shown in Table <ref>.
Note that since the estimation is based on atomic spectral data and intra-band excitations within valence electrons are not considered, it can only serve as rough estimation and all-electron calculation is ultimately needed to study the contribution of semi-core electrons.





 §.§ Calculation of electronic stopping for an off-channeling trajectory




The energy transfer from proton to the target material depends on the local environment and is trajectory dependent.
When a proton is closer to nuclei of the target material, it has higher probability to excite electrons since the electron density is higher. 
The proton also has higher chance to excite core electrons for the same reason.
Therefore, the energy transferred from proton to target materials is larger when proton travels near positions of nuclei.
While the shortest distance between proton and target ion leads to the sharpest peaks in Fig. <ref>, we also note that these results are affected by the cutoff radius of the pseudopotentials used here.
For this reason, we only use the average to extract stopping, as explained in detail in Ref. Schleife2015.

Nevertheless, counting the peaks in each trajectory, we clearly find that there are much more and higher peaks for the trajectory, on which proton travels in InP, than the other two trajectories.
This indicates that for the short trajectory we used to calculate the electronic stopping of InP, the proton happens to experience region of higher electron density.
Therefore, we predict higher electronic stopping than fully converged value for InP.
Decreasing height of the peaks for InP also suggest that the trajectory start to explore region of lower electron density.
Therefore, we expect a much longer trajectory can have better sampling of the target materials and thus predict electronic stopping closer to converged value.
 


 §.§ Error estimate for channeling projectiles




