Ab initio electronic stopping power and threshold effect of channeled slow light ions in HfO_{2}
Chang-Kai Li, Feng Wang, Bin Liao, Xiao-Ping OuYang, Feng-Shou Zhang
aa r X i v : . [ phy s i c s . a t m - c l u s ] J un Ab initio electronic stopping power and threshold effect of channeled slow light ions inHfO Chang-Kai Li , , Feng Wang , Bin Liao , , Xiao-Ping OuYang and Feng-Shou Zhang , , , ∗ The Key Laboratory of Beam Technology and Material Modification of Ministry of Education,College of Nuclear Science and Technology,Beijing Normal University Beijing 100875, China Beijing Radiation Center, Beijing 100875, China School of Physics, Beijing Institute of Technology,Beijing 100081, China Center of Theoretical Nuclear Physics,National Laboratory of Heavy Ion Accelerator of Lanzhou,Lanzhou 730000, China Northwest Institute of Nuclear Technology,Xi’an 710024, China (Dated: June 29, 2017)We present ab initio study of the electronic stopping power of protons and helium ions in aninsulating material, HfO . The calculations are carried out in channeling conditions with differ-ent impact parameters by employing Ehrenfest dynamics and real–time, time–dependent densityfunctional theory. The satisfactory comparison with available experiments demonstrates that thisapproach provides an accurate description of electronic stopping power. The velocity-proportionalstopping power is predicted for protons and helium ions in the low energy region, which conformsthe linear response theory. Due to the existence of wide band gap, a threshold effect in extremelylow velocity regime below excitation is expected. For protons, the threshold velocity is observable,while it does not appear in helium ions case. This indicates the existence of extra energy loss chan-nels beyond the electron–hole pair excitation when helium ions are moving through the crystal. Toanalyze it, we checked the charge state of the moving projectiles and an explicit charge exchangebehavior between the ions and host atoms is found. The missing threshold effect for helium ions isattributed to the charge transfer, which also contributes to energy loss of the ion. PACS numbers: 61.85.+p,31.15.A − ,61.80.Az,61.82.Ms I. INTRODUCTION
Interaction of intruding ion with matter is of contin-uing interest, the deceleration force involved is mainlycharacterized by stopping power (SP) [1] or energy de-position per unit distance traveled through the mate-rial. This quantity can generally be decomposed intotwo parts. The first one is the electronic stopping power S e , which arises from the excitation of the electrons ofthe target. The second one is the nuclear stopping power S n due to elastic Coulomb collisions with the nuclei ofthe target. For particles with velocities below the Fermivelocity of the target, nuclear and electronic stoppingare both relevant, and the result of the interaction is acollision cascade [2]. Shortly before the particle stopseventually, a global stopping maximum denoted Braggpeak occurs. Thus, studying the energy transferred fromslow ions to the target material is at the heart of moderntechnologies, such as nuclear fission/fusion reactors [3, 4],semiconductor devices for space missions [5] and cancertherapy based on ion beam radiation [6, 7].The interest in modeling the SP of ions with velocitiesbetween 0.1 and 1 atomic units (a.u. hereafter) has fu- ∗ Corresponding author. [email protected] eled a huge amount of research [8]. In this regime theelectronic component is generally dominant and the elec-tronic energy loss is predominantly due to interactionwith valence electrons [9]. The electronic stopping poweris a crucial quantity for ion irradiation: it governs thedeposited heat, the damage profile, and the implantationdepth. Ever since the discovery of S e , various modelsand theories have been proposed to calculate S e depend-ing on the energy regime of the ion. For swift ions, basedon the assumption that the atoms of the target are clas-sical oscillator, Bohr [10, 11] suggested a formula for the S e ; employing the first Born approximation, Bethe [12]has introduced the first calculations of inelastic and ion-ization cross section; the Bloch correction [13] provides aconvenient link between the Bohr and the Bethe scheme.On the other hand, for slow particles, Fermi and Teller[14] using electron gas models had reported S e for varioustargets; a more general treatment of the S e applicable tothe whole velocity regime was later developed by Lind-hard [15] through linear dielectric theory.Energetic ions penetrate great depth along channelsbetween low–index crystallographic planes, moderatingthrough collisions with electrons, until finally they hita host atom initiating a cascade of atomic displace-ments. This channeling phenomenon has been exploitedin many important applications such as ion implantationand depth profiling. Glancing collisions with host atomsconfine the trajectory of a channelling ion, so most of itsenergy is lost through electronic excitation. Based on thefree electron gas (FEG) model, S e is predicted to be S e ∝ v for v < sp –bonded metals [16–19].A different behavior is expected in materials that havea finite minimum excitation energy T min , such as noblegas and wide band–gap insulators, given the finite energyrequired to excite outmost electrons. This finite mini-mum excitation energy is expected to suppress the energydispassion due to electron–hole pair excitation at ion en-ergies of several keV [14], which results in a thresholdeffect of SP with respect to the ion velocity [20]. Instead,however, no threshold effect was originally observed inmost systems [21–23], with the exception that a thresh-old velocity of v ≃ . v ≤ S e behaviorof slow ions shooting through a large band–gap insulatorHfO ( E g ≈ S e canbe extracted from the change of kinetic projectile energyusing the thickness of the target.This article is outlined as follows. In Sec. II, we brieflyintroduce the theoretical framework and the computa-tional details. Results are presented and discussed inSec. III, where we concentrate on the analysis of S e andits threshold effect. In the end, conclusions are drawnin Sec. IV. Henceforth, if no special reservation is made,Hartree atomic units with m = | e | = ~ = 1 are used inthis paper. II. MODEL AND METHODS
The collision behaviors of intruding ions with the hostnuclei and electrons are characterized by the Ehren-fest coupled electron–ion dynamics combined with time–dependent density–functional theory (ED–TDDFT) [29–34]. In this model, electronic degrees of freedomare treated quantum–mechanically within the time–dependent Kohn–Sham (KS) formalism, while the ionsare handled classically. This method allows for an excitedelectronic state ab initio molecular dynamics (AIMD)simulation. The ED–TDDFT can, in general, be definedby the following coupled differential equations: M I d ~R I ( t ) dt = − Z Ψ ∗ ( x, t )[ ∇ I ˆ H e ( ~r, ~R ( t ))]Ψ( x, t ) dx −∇ I X I = J Z I Z J | ~R I ( t ) − ~R J ( t ) | , (1) i ∂ Ψ( x, t ) ∂t = ˆ H e ( ~r, ~R ( t ))Ψ( x, t ) , (2)here M I and Z I denote the mass and charge of the I thnuclei, respectively, and ~R I ( t ) describes the correspond-ing ionic position vector. Ψ( x, t ) is the many–body elec-tron wave function in the time domain, for which wedefine x ≡ { x j } Nj =1 , with x j ≡ ( ~r j , σ j ), where the coordi-nates ~r j and the spin σ j of the j th electron are implicitlytaken into account. Here N is the number of electrons ofthe system.The electronic Hamiltonian is expressed as ˆ H e ( ~r, ~R ( t )),which depends on the instantaneous distribution of thepositions of all the nuclei, ~R ( t ) ≡ { ~R ( t ) , · · · , ~R M ( t ) } ( M is the number of nuclei of the system), and of allthe electrons ~r ; thus, it basically consists of the kineticenergy of electrons, the electron–electron potential andthe electron–nuclei potential, which can be formulatedas ˆ H e ( ~r, ~R ( t )) = − N X j ∇ j + X i
FIG. 1. Electronic stopping power for protons (a) and he-lium ions (b) versus velocity along the middle axes of < > , < > and < > channels of HfO . The lines are guides tothe eye. The inset in (a) is enlargement of the main figure. A. Threshold effect in S e of HfO To investigate the different threshold effects betweenthe protons and helium ions, we have explored the pos-sible energy dissipation channels. In the present work,we find S e is related to charge transfer, which is an ad-ditional energy dissipation channel besides electron–holepair excitation. Many mechanisms may contribute tothe charge transfer. Besides the direct transitions suchas excitation, ionization, and capture [53], the Auger pro-cess between the host atoms and ions also plays a pro-nounced role, in which an electron jumps from the va-lence band of the host atom to an ion bound state andvice versa. The energy released in such transitions isbalanced by an electronic excitation in the medium oron the projectile [54]. Other possible mechanisms areresonant charge transfer and radiative decay processesof the projectiles. Since pseudopotentials are adopted in the present TDDFT simulation, Auger processes follow-ing the inner–shell vacancy can not be considered. Nev-ertheless, because the kinetic ion energy in the presentwork is restricted to 25 keV/u and lower, according tothe interpretation in Ref. [55], direct transition mecha-nisms are dominant in such a low velocity regime, Augerprocesses following the inner–shell vacancy make a minorcontribution to charge transfer.Generally, both the neutralization and re-ionization ofthe ion included in the charge transfer process contributeto the decreasing of ion’s kinetic energy, due to the pro-motion of electronic states of either the host atoms or theprojectile itself. Pe´nalba et al. [56] reported that chargetransfer is an important energy loss channel, especiallyfor projectile around stopping maximum. For protonswith v = 1 a.u. in Aluminum, charge transfer accountsfor 15% of total SP. However, it has not been consideredin the original SP theory that accounts for linear velocitydependence. The different threshold effect between pro-ton and helium ion in Fig. 1 is studied by checking theircharge transfer behaviors in low velocity regime.As a first step, the time evolution of a helium ion mov-ing through the < > channel for a given velocity of1.0 a.u. is visualized in Fig. 2. Four snapshots coveringthe entire collision process are presented. Before enter-ing the crystal ( t = 0.059 fs), the helium ion projectileis a bare ion [Fig. 2(a)]; when the helium ion is gettingclose to the crystal ( t = 0.095 fs) and penetrating alongthe channel ( t = 0.486 fs), it exchanges charge with thehost atoms [Figs. 2(b) and (c)]; after traversing the HfO film ( t = 0.628 fs); the exiting ion still retains some cap-tured electrons [Fig. 2(d)]. The example in Fig. 2 couldbe qualitative evidence for the charge transfer during thecollisions. The gray region is the change in electron dis-tribution caused by the intruding ion.The electron density change induced by the heliumion in real time is quantified by integrating the valencecharge density within the volume around the projectileion with a radius of 1.26 ˚A in the time–dependent calcula-tion, from which the ground state electron density of thetarget in the corresponding volume has been subtractedand we thus obtained the number of induced electronsin real time. In the present work, we deem this quan-tity as the charge transferred to the projectiles. A pointshould be noted is that free electrons caused by electronscattering process may also be included in such approx-imation, where electrons pile up close to projectile dueto the attractive interaction between electrons and theions [57]. The choice of 1.26 ˚A as the integration radiusis a compromise between various factors. In this study,we are interested in learning the real–time electron occu-pying the intruding ion orbitals, and we get it throughthe discrepancy of density, i.e., the change in electrondistribution around the ion between the time–dependentand the ground state calculations. In theory, a largerintegration radius can be more effective to fully take avariety of mechanisms and also the highly occupied or-bits into account. However, at the same time, it may (a) (b)(c) (d) FIG. 2. (Color online) Snapshots of time evolution of theelectron density change caused by a He ion moving throughHfO crystal with v = 1.0 a.u. along the middle axis of < > channel (side view). (a) t = 0 .
059 fs, the ion is above thecrystal. (b) t = 0 .
095 fs, the ion is entering the channel.(c) t = 0 .
486 fs, the ion is penetrating along the channel.(d) t = 0 .
628 fs, the leaving ion retains some electrons afterthe collision. The gray regions are the change of electrondistribution caused by the intruding ions. H + HfO (a) -10 -8 -6 -4 -2 0 2 4 6 8 100.00.51.01.5 He HfO (b) I ndu ce d c h a r g e ( un it s o f | e - | ) z (¯) FIG. 3. (Color online) Electrons captured by protons (a) andhelium ions (b) versus z coordinates in low velocity regimealong the middle axis of < > channel. The vertical dashedlines are the positions of O layers, the vertical solid lines arethe positions of Hf layers. include more free electrons and excited state electrons ofthe host atoms that do not belong to the ion and alsomore excited state electrons caused by the former steps,as shown in Fig. 2.Electrons induced by the H + and He ions with verylow velocities moving along the middle axis of < > channel are presented in Fig. 3; the periodic variation ininduced electron reflects the periodicity of the crystal.As can be seen, the charge exchange between the pro-jectile and host atoms takes place alternately along theirtrajectories. Generally, protons show less active chargeexchange behavior during the collision than helium ions.In the present work, we found charge transfer of protonsis sensitive to ions velocity, it becomes less and less obvi-ous as the velocity decreases in the extremely low veloc-ity regime (not shown). While, charge transfer of heliumions shows a weak correlation with velocity. Consider-ing charge transfer accounts for a noticeable share of S e ,in the extremely low velocity regime the contribution to S e from electron–hole pair excitation becomes feeble andcharge transfer may be the dominant energy loss chan-nel. The different threshold effect shown in Fig. 1 can beattributed to the different charge transfer behaviors be-tween two kinds of ions in extremely low velocity regime.To depict the electron–hole pair excitation, we calcu-lated the electron occupation number distribution in theconduction band for protons at a given velocity of 0.01a.u. along the middle axis of < > channel and 0.5 a.u.along the middle axis of < > channel. Such choices aremade with the consideration to compare the cases withand without explicit electron–hole excitation, because0.01 a.u. in < > channel is lower than the thresholdvelocity and 0.5 a.u. in < > channel is higher thanthreshold velocity. Results are shown in Fig. 4. The ex-citation energy of the crystal has been extended to 7 eVin the conduction band by adding empty KS obtials inthe simulation. We can see clearly that the excited elec-trons are distributed broadly in the considered energyrange for the proton at velocity of 0.5 a.u. in < > channel. However, the electrons occupation in conduc-tion band is implicit for the proton at velocity of 0.01a.u. in < > channel, and the number of excited elec-trons in this energy range is much smaller than that ofthe 0.5 a.u. proton.As already alluded to in the previous paragraph, elec-trons excitation from valence band to the conductionband are suppressed by the energy gap. It should benoted that charge transfer discussed above also con-tributes to electrons excitation. The implicit electronsexcitation for 0.01 a.u. proton also demonstrates themissing of charge transfer in very low velocity regime.This means the projectile ion does not lose any energy tothe target electron subsystem when the velocities are be-low a certain threshold, resulting in a threshold effect inthe S e dependence on the velocity. Although the emptyKS level does not provide the exact description of thereal excited states of the nanostructure, this can be con-sidered as the first approximation to the true excitationsof the system [58, 59]. FIG. 4. Electron occupation number distribution after thecollision ends. The left side of the vertical lines shows theoccupation of the ground state and the right side shows theexcited states. The black line and red one are cases of protonsat 0.01 a.u. in < > channel and 0.5 a.u. in < > channel,respectively. See more details in the text. B. impact parameter dependence of chargetransfer and S e Since the occupied He-1 s level is strongly affected bythe interaction distance [60–63]. So charge transfer be-havior may differ for different trajectories. To investigatethe effect of impact parameter on charge exchange andthreshold velocity, we show in Figs. 5 and 6 the S e behav-ior and the corresponding charge transfer profiles versus z coordinates for a given velocity of 0.3 a.u. along twotrajectories in < > channel and three trajectories in < > channel. The trajectories are chosen to sampledifferent impact parameters within the channel. In thepresent work, the value of impact parameter is definedas the closest distance to any of Hf atoms along the iontrajectory. The trajectories along the middle axes of thechannel have the highest impact parameters and trajecto-ries close to Hf atoms have the lowest impact parameters.For the five trajectories in Fig. 5, the impact parametersare in order 1.279 ˚A, 0.904 ˚A, 1.809 ˚A, 1.357 ˚A and 0.904˚A.As can be seen, no threshold effect can be found in < > channel in both the center and off–center chan-neling cases shown as trajectories 1 and 2, respectively.While, threshold effect is observed along the trajecto-ries with low impact parameters in the < > channel,i.e. trajectories 4 and 5, the threshold velocity are 0.05 a.u. and 0.15 a.u., respectively. We notice that thecharge transfer behavior shown in Fig. 6 is active for thetrajectory in < > channel with low impact parame-ter. While, the amplitudes of charge transfer are muchsmaller along the trajectories in < > channel with lowimpact parameters, and trajectory 5 shows a even lessevident charge transfer behavior than trajectory 4. Thisresult consist with the threshold effect shown in Fig. 5,which indicates that an active charge transfer may causevanish of the threshold. He HfO <100> (a) He HfO <110> (b) d E / dx ( e V / ¯ ) Velocity (a.u.)
FIG. 5. (Color online) Electronic stopping power of heliumions in HfO as a function of the velocity along two trajec-tories in the < > channel (a) and three trajectories in the < > channel (b). The lines are guides to the eye. The insetshows a sectional view of the < > channel and the trajecto-ries. The gray circles and the blue ones represent host atomsin different transverse planes (defining the channel), while theblack circles show the projectile positions for different impactparameters. For the five trajectories impact parameters arein order 1.279 ˚A, 0.904 ˚A, 1.809 ˚A, 1.357 ˚A and 0.904 ˚A. The instantaneous energy loss rate of the projectileion in a condensed matter system often depends stronglyon the specific path taken by the ion and its proximityto atoms and bonds over the course of the trajectory.In the present work, we find S e is remarkably impactparameter–dependent, results along different trajectorieshave different magnitudes. To understand the depen-dence of the S e on the impact parameter, we show inFig. 7 the ground state electronic density and the axialforce along along three different ion trajectories (shownas black circles in Fig. 5(b)) in < > channel for heliumions at a given velocity of 0.3 a.u.. The density values areobtained by averaging electron density of cylindrical iontrack with a radius of 0.36 ˚A step by step. The trend ofthe axial force is similar to that of electron density, andthere is a proportional relation between the two to some HfO <100>(a) I ndu ce d c h a r g e ( un it s o f | e - | ) -10 -8 -6 -4 -2 0 2 4 6 8 100.00.51.01.5 HfO <110>(b) z (¯) FIG. 6. Electron capture by He along trajectories withdifferent impact parameter in the < > channel (a) and inthe < > channel (b) for a given velocity of 0.3 a.u.. Thecurves are corresponding to trajectories in Fig. 5. The verticaldashed lines and solid lines in (a) are the positions of O layersand Hf layers, respectively. The vertical dash-dotted lines in(b) are positions of crystalline layers comprising both O andHf atoms. extent, which is in accordance with the density functionaltheory (DFT) results [64] for FEG. The average densityvalues obtained along the whole trajectories in < > , < > , and < > channels are 0.35, 0.33, and 0.28electrons/˚A , respectively, which consists with the am-plitude of S e shown in Fig. 1(b). This suggests that the S e in channeling conditions is related to the average den-sity along the projectile’s trajectory, corroborating theinterpretation in the literatures [27, 65–68].The trend of force also have a direct correlation withthe threshold effect. As can be seen, the applied force onhelium ion varies like a sine function with respect to thedisplacement, and trajectories with the low impact pa-rameters have a larger variation. For trajectories 3 and4, the force is generally above zero, i.e. positive z direc-tion, which means the ions are predominantly exposedto drag force along the channel. For trajectory 5, theforce exerting on the ion turns up and down periodically(in the z axis direction), leading to little net energy lossin the following path, which is consistent with the rela-tively small threshold velocity of trajectory 5 shown inFig. 5(b).It should be noted that, since projectile have distinctcharge transfer behaviors along different trajectories, andcharge transfer contribute significantly to S e . So the ef-fect of impact parameter on S e is also embodied throughthe charge transfer. <110> channel D e n s it y ( ¯ - ) (a) -10 -8 -6 -4 -2 0 2 4 6 8 10-80-60-40-200204060 HfO <110> (b) A x i a l f o r ce ( e V / ¯ ) z (¯) FIG. 7. (Color online) Electron densities (a) and axial force(b) exerted on helium ion for a given velocity of 0.3 a.u. versus z coordinates in HfO along three trajectories in < > . Thecurves are corresponding to trajectories in Fig. 5(b). The ver-tical dashed lines are positions of crystalline layers comprisingboth O and Hf atoms. IV. CONCLUSIONS
Theoretical study from first principles the electronicstopping power of slow light projectiles in HfO has beenpresented. The velocity–proportional S e of HfO is pre-dicted. A quantitative agreement between the experi-mental data and our results is achieved. Threshold effectis found when proton is channeled in the HfO thin film,while, the expected threshold effect was not found for he-lium ion channelling the crystal, which was interpreted asa consequence of charge transfer. We have learned thatthe S e is sensitive to the impact parameter due to thedifferent electron density experienced by the ions, whichis consistent with the assumptions form FEG model. Ourresults shed light on describing the interaction betweenthe ions and the target electrons without restricting theelectrons to the adiabatic surface. To obtain a deeperunderstanding of the effect of charge transfer on inelas-tic energy loss, a thorough theoretical analysis of possiblecharge exchange mechanism in combination with suitableexperimental studies is highly desirable. ACKNOWLEDGEMENTS
One author wants to thank Dr. F. Mao for thefruitful discussions. This work was supported bythe National Natural Science Foundation of China un-der Grant Nos. 11635003, 11025524 and 11161130520,National Basic Research Program of China underGrant No. 2010CB832903, and the European Commis- sions 7th Framework Programme (FP7-PEOPLE-2010-IRSES) under Grant Agreement Project No. 269131. [1] T. L. Ferrell and R. H. Ritchie, Phys. Rev. B , 115(1977).[2] R. S. Averback and T. D. D. L. Rubia, Solid State Phys , 281 (1997).[3] G. R. Odette and B. D. Wirth, in Handbook of MaterialsModeling: Methods (Springer Netherlands, Dordrecht,2005) p. 999.[4] A. P. Horsfield, A. Lim, W. M. C. Foulkes, and A. A.Correa, Phys. Rev. B , 245106 (2016).[5] H. Zeller, Solid-State Electronics , 2041 (1995).[6] O. I. Obolensky, E. Surdutovich, I. Pshenichnov,I. Mishustin, A. V. SolovYov, and W. Greiner, Nucl.Instrum. Methods B , 1623 (2008).[7] K. G. Reeves, Y. Yao, and Y. Kanai, Phys. Rev. B ,041108 (2016).[8] R. Cabrera-Trujillo, P. Apell, J. Oddershede, andJ. Sabin, A I P Conference Proceedings CP680 , 86(2003).[9] A. Mertens and H. Winter, Phys. Rev. Lett. , 2825(2000).[10] N. Boh, Philos. Mag. , 581 (1915).[11] N. Boh, Philos. Mag. , 10 (1913).[12] H. Bethe, Annalen der Physik , 325 (1930).[13] F. Bloch, Annalen der Physik , 285 (1933).[14] E. Fermi, E. Teller, and V. Weisskopf, Phys. Rev. ,399 (1947).[15] J. Lindhard, Mat. Fys. Medd. Dan. Vid. Selsk, , 41(1954).[16] C. P. Race, D. R. Mason, M. W. Finnis, W. M. C.Foulkes, A. P. Horsfield, and A. P. Sutton, Rep. Prog.Phys. , 116501 (2010).[17] J. E. Vald´es, G. A. Tamayo, G. H. Lantschner, J. C.Eckardt, and N. R. Arista, Nucl. Instrum. Methods B , 313 (1993).[18] G. Mart´ınez-Tamayo, J. C. Eckardt, G. H. Lantschner,and N. R. Arista, Phy. Rev. A , 3131 (1996).[19] J. M. Pitarke and I. Campillo, Nucl. Instrum. MethodsB , 147 (1999).[20] S. N. Markin, D. Primetzhofer, and P. Bauer, Phys. Rev.Lett. , 113201 (2009).[21] K. Eder, D. Semrad, P. Bauer, R. Golser, P. Maier-Komor, F. Aumayr, M. Pe˜nalba, A. Arnau, J. M. Ugalde,and P. M. Echenique, Phys. Rev. Lett. , 4112 (1997).[22] M. Pe˜nalba, J. I. Juaristi, E. Zarate, A. Arnau, andP. Bauer, Phys. Rev. A , 012902 (2001).[23] S. P. Møller, A. Csete, T. Ichioka, H. Knudsen, U. I.Uggerhøj, and H. H. Andersen, Phys. Rev. Lett. ,042502 (2004).[24] C. Auth, A. Mertens, H. Winter, and A. Borisov, Phys.Rev. Lett. , 4831 (1998).[25] L. N. Serkovic, E. A. S´anchez, O. Grizzi, J. C. Eckardt,G. H. Lantschner, and N. R. Arista, Phys. Rev. A ,040901 (2007).[26] L. N. Serkovic Loli, E. A. S´anchez, O. Grizzi, and N. R.Arista, Phys. Rev. A , 022902 (2010).[27] R. Ullah, F. Corsetti, D. S´anchez-Portal, and E. Arta-cho, Phys. Rev. B , 125203 (2015). [28] F. Mao, C. Zhang, J. Dai, and F.-S. Zhang, Phys. Rev.A , 022707 (2014).[29] M. A. Zeb, J. Kohanoff, D. S´anchez-Portal, A. Arnau,J. I. Juaristi, and E. Artacho, Phys. Rev. Lett. ,225504 (2012).[30] A. P. Horsfield, D. R. Bowler, A. J. Fisher, T. N.Todorov, and M. J. Montgomery, J. Phys. Condens. Mat-ter , 3609 (2004).[31] E. K. U. Gross, J. F. Dobson, and M. Petersilka, “Den-sity functional theory of time-dependent phenomena,” in Density Functional Theory II: Relativistic and Time De-pendent Extensions , edited by R. F. Nalewajski (SpringerBerlin Heidelberg, Berlin, Heidelberg, 1996) pp. 81–172.[32] F. Calvayrac, P. G. Reinhard, E. Suraud, and C. A.Ullrich, Phys. Rep. , 493 (2000).[33] J. L. Alonso, X. Andrade, P. Echenique, F. Falceto,D. Prada-Gracia, and A. Rubio, Phys. Rev. Lett. ,096403 (2008).[34] J. L. Page, D. R. Mason, and W. M. C. Foulkes, J. Phys.Condens. Matter , 611 (2008).[35] A. Castro, M. Isla, J. I. Mart´ınez, and J. Alonso, Chem.Phys. , 130 (2012).[36] J. P. Perdew and Y. Wang, Phys. Rev. B , 13244(1992).[37] X. Andrade, A. Castro, D. Zueco, J. L. Alonso,P. Echenique, F. Falceto, and n. Rubio, J. Chem. TheoryComput. , 728 (2009).[38] G. Avenda˜no-Franco, B. Piraux, M. Gr¨uning, andX. Gonze, Theore. Chem. Acc. , 1289 (2012).[39] N. Troullier and J. L. Martins, Phys. Rev. B , 1993(1991).[40] T. Otobe, M. Yamagiwa, J.-I. Iwata, K. Yabana,T. Nakatsukasa, and G. F. Bertsch, Phys. Rev. B ,165104 (2008).[41] Y. Jiao, F. Wang, X. Hong, W. Su, Q. Chen, andF. Zhang, Physics Letters A , 823 (2013).[42] K. Yabana and G. F. Bertsch, Phys. Rev. B , 4484(1996).[43] A. A. Shukri, F. Bruneval, and L. Reining, Phys. Rev.B , 035128 (2016).[44] A. Castro and M. A. L. Marques, Propagators for theTime-Dependent Kohn-Sham Equations (Springer BerlinHeidelberg, 2004) pp. 197–210.[45] M. A. L. Marques, A. Castro, G. F. Bertsch, and A. Ru-bio, Comput. Phys. Commun. , 60 (2003).[46] A. Castro, H. Appel, M. Oliveira, C. A. Rozzi, X. An-drade, F. Lorenzen, M. A. L. Marques, E. K. U. Gross,and A. Rubio, Phys. Status Solidi B , 2465 (2006).[47] F. Wang, X. C. Xu, X. H. Hong, J. Wang, and B. C.Gou, Phys. Lett. A , 3290 (2011).[48] A. A. Correa, J. Kohanoff, E. Artacho, D. S´anchez-Portal, and A. Caro, Phys. Rev. Lett. , 213201(2012).[49] R. W. G. Wyckoff,
Crystal Structures Second edition (In-terscience Publishers, 1963) pp. 239–444.[50] D. Primetzhofer, Nucl. Instrum. Methods B , 100(2014). [51] A. Schleife, Y. Kanai, and A. A. Correa, Phys. Rev. B , 014306 (2015).[52] E. E. Quashie, B. C. Saha, and A. A. Correa, Phys. Rev.B , 155403 (2016).[53] G. Schiwietz and P. L. Grande, Phys. Rev. A , 052703(2011).[54] R. D. Mui˜no, Nucl. Instrum. Methods B , 8 (2003),14th International Workshop on Inelastic Ion-SurfaceCollisions.[55] C. L. Zhang, X. H. Hong, F. Wang, Y. Wu, and J. G.Wang, Phys. Rev. A , 032711 (2013).[56] M. Pe´nalba, A. Arnau, P. M. Echenique, F. Flores, andR. H. Ritchie, Europhys. Lett. , 45 (1992).[57] P. Echenique, R. Nieminen, and R. Ritchie, Solid StateCommun. , 779 (1981).[58] H. Appel, E. K. U. Gross, and K. Burke, Phys. Rev.Lett. , 043005 (2003).[59] A. Wasserman and K. Burke, Phys. Rev. Lett. , 163006(2005).[60] N. P. Wang, E. A. Garc´ıa, R. Monreal, F. Flores, E. C. Goldberg, H. H. Brongersma, and P. Bauer, Phys. Rev.A , 012901 (2001).[61] S. Wethekam, D. Vald´es, R. C. Monreal, and H. Winter,Phys. Rev. B , 075423 (2008).[62] R. Monreal, E. Goldberg, F. Flores, A. Narmann,H. Derks, and W. Heiland, Sur. Sci. Lett. , 271(1989).[63] D. Primetzhofer, D. Goebl, and P. Bauer, Nucl. Instrum.Methods B , 8 (2013).[64] P. M. Echenique, R. M. Nieminen, J. C. Ashley, andR. H. Ritchie, Phys. Rev. A , 897 (1986).[65] N. P. Wang and I. Nagy, Phys. Rev. A , 4795 (1997).[66] N. P. Wang, I. Nagy, and P. M. Echenique, Phys. Rev.B , 2357 (1998).[67] H. Winter, J. I. Juaristi, I. Nagy, A. Arnau, and P. M.Echenique, Phys. Rev. B , 245401 (2003).[68] L. Martingondre, G. A. Bocan, M. Blancorey, M. Al-ducin, J. I. Juaristi, and R. D. Mui˜no, J. Phys. Chem.C117