Pre-equilibrium stopping and charge capture in proton-irradiated aluminum sheets
PPre-equilibrium stopping and charge capture in proton-irradiated aluminum sheets
Alina Kononov
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
Andr´e Schleife ∗ Department of Materials Science and Engineering,University of Illinois at Urbana-Champaign, Urbana, IL 61801, USAMaterials Research Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA andNational Center for Supercomputing Applications,University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA (Dated: June 25, 2020)We present a first-principles study of pre-equilibrium stopping power and projectile charge capturein thin aluminum sheets irradiated by 6 – 60 keV protons. Our time-dependent density functionaltheory calculations reveal enhanced stopping power compared to bulk aluminum, particularly nearthe entrance layers. We propose the additional excitation channel of surface plasma oscillations asthe most plausible explanation for this behavior. We also introduce a novel technique to compute theorbital-resolved charge state of a proton projectile after transmission through the sheet. Our resultsprovide insight into the dynamics of orbital occupations after the projectile exits the aluminum sheetand have important implications for advancing radiation hardness and focused-ion beam techniques,especially for few-layer materials.
I. INTRODUCTION
Although ion irradiation represents a fundamentalmeans of manipulating materials and probing theirproperties, a detailed theoretical understanding ofthe interaction between an energetic charged particleand induced electronic excitations in a solid has provenchallenging. The Bethe formula models stoppingpower S = − dE/dz (i.e., the energy deposited per pen-etration depth) for fast projectiles, while the Fermi-Teller model applies to slow projectiles. The Lindhardstopping formula from linear-response theory and itsextensions accurately account for the small density per-turbations produced by fast, low- Z projectiles. However,these existing analytical models rely on ambiguousparameters such as projectile charge Z or participatingelectron density n which not only have multiple meaning-ful definitions, but may also evolve dynamically duringthe projectile’s transit. In addition, upon entering the material, quantitiessuch as stopping power and the projectile’s charge stateare not initially identical to the steady-state bulk valuespresumably achieved as the projectile moves through athick target. This leads to pre-equilibrium effects thatcannot be ignored when understanding electronic stop-ping near surfaces or in thin target materials of onlya few atomic layers. In these, the projectile may notreach an equilibrium charge state at all, since it mayneed to traverse many layers before doing so. While pre-equilibrium effects should occur even for light ions, theyshould be most prominent for highly charged projectileswith a large difference between initial and equilibratedcharge state.Indeed, several experiments on highly charged ions im-pacting carbon-based materials with thicknesses of only afew nm or less reported that the projectile’s initial charge influences energy and charge transfer distributions and even material damage characteristics.
Similarly,the response of aluminum foils to ion irradiation hasbeen shown to depend on both ion charge and foilthickness.
Such studies inferred projectile chargeequilibration time scales smaller than 10 fs and lengthscales shorter than 10 nm.
Sensitivity of elec-tron emission to incident ion charge was shown evenfor ∼
100 nm thick carbon foils and attributed to pre-equilibrium stopping and projectile charge. These experimental observations of pre-equilibrium ef-fects inspired exponential decay models for projectilecharge equilibration.
Since experiments cannot ac-cess the projectile’s charge state within the material,studies evaluating such models typically compare theirpredictions to measurements of the projectile’s charge af-ter transmission through the sample.
However, thisexit charge state may not be equivalent to the projectile’scharge state inside the material. Overall, the transitionand equilibration of the ion into the bulk regime is stillpoorly understood and requires further study.Conversely, for bulk materials under ion irradi-ation, stopping power has been fairly well-studiedexperimentally.
In addition, modern first-principlessimulations have provided a detailed description of theenergy and charge dynamics in bulk, as evident frommany recent studies of electronic stopping power inmetals, semiconductors, and insulators which were enabled by the rise of high performance com-puting. However, it remains difficult for experiments toaccess the poorly understood pre-equilibrium behaviorof stopping power near a surface or for materials thinnerthan the length scale of projectile charge equilibration.To improve our understanding of pre-equilibrium be-havior, computational modeling of thin materials underion irradiation is a promising alternative, as it offers a r X i v : . [ c ond - m a t . m t r l - s c i ] J un greater spatial and temporal resolution than currentlyachievable experimentally. However, modeling an ion’sinteraction with a surface requires the inclusion of a suf-ficiently large vacuum region and sufficiently large ma-terial surface, greatly increasing computational cost. Inaddition, extracting observables from the simulated time-dependent electron density to describe the charge dy-namics instigated by ion-induced electronic excitationspresents a further challenge. Since detailed understand-ing of pre-equilibrium behavior is currently absent, weaim to study these effects in the present work.To this end, we used a first-principles computationalapproach to calculate and analyze electronic stoppingand projectile charge state as protons traverse aluminumsheets. We focus on protons with kinetic energies of 6 –60 keV that move along a channeling trajectory through0.8 – 2.4 nm of aluminum, corresponding to 4 – 12 atomiclayers. The wealth of existing literature on the elec-tronic response of bulk aluminum to proton irradia-tion makes this an ideal system to study dynamicalbehavior near ion-irradiated surfaces. In particular,the stopping power of protons in bulk aluminum hasbeen studied extensively both experimentally andtheoretically. This existing wisdom provides op-portunities both to validate our results for the bulk limitand to clearly identify surface effects and pre-equilibriumbehavior. Moreover, light-ion irradiation is particularlywell-suited for first-principles studies because they ex-perience relatively little scattering in the host material,resulting in long, straight trajectories. Using a protonprojectile further allows us to accurately calculate thecharge captured by the projectile using analytic hydro-gen orbitals, as we discuss below.The remainder of this paper is structured as follows:In Section II we outline the first-principles computa-tional approach used here, and in Section III we describethe post-processing methods developed and used to ex-tract charge capture from time-dependent electron den-sity data. Section IV presents the results obtained forpre-equilibrium stopping power and charge capture, andSection V summarizes the conclusions drawn from thisstudy.
II. COMPUTATIONAL METHODS
We used real-time time-dependent density functionaltheory (TDDFT) as implemented in theQbox/Qb@ll code to simulate the real-time non-adiabatic electronic-ion dynamics as a proton traversesthin aluminum sheets. Different sheet thicknesses con-sisting of 4 – 12 atomic layers were studied, where onelayer corresponds to about 0.2 nm. The Kohn-Sham (KS)orbitals were expanded in a plane-wave basis with a cut-off energy of 50 Ry. The electron-ion interaction wasdescribed using a Troullier-Martins pseudopotential with 11 valence electrons for aluminum and a Hamann-Schl¨uter-Chiang-Vanderbilt pseudopotential for the proton projectile. The large simulation cells used inthis work, typically 38 × ×
150 a , allow reciprocal-space sampling using the Γ point only. The adiabaticlocal density approximation was used to describe ex-change and correlation. The electronic ground state ofthe aluminum sheet served as the initial condition of thetime-dependent calculations, and it was computed usingdensity functional theory with 100 empty states and aFermi electron temperature of 100 K in order to accel-erate self-consistent convergence of the metallic groundstate.Due to the few-fs time scale of our time-dependent sim-ulations, ions do not have enough time to respond to theforces they experience. Thus, we fix the positions of thealuminum ions and maintain a constant velocity of theproton projectile. The proton starts at a distance of atleast 25 a from the aluminum sheet and traverses it alonga [100] channeling trajectory (see Fig. 1). The enforcedtime reversal symmetry (ETRS) integrator with atime step of 0.042 at. u. (1.04 as) was used to propagatethe Kohn-Sham orbitals, a choice shown to produce highnumerical accuracy for similar systems. The cutoff en-ergy, treatment of semi-core electrons, target atom dy-namics, and channeling projectile trajectory used in thiswork were shown in Ref. 26 to produce accurate stop-ping power results for proton-irradiated bulk aluminumwithin the velocity regime presently considered. Accu-racy and convergence of the time propagation for theselarge simulation cells is addressed in the SupplementalMaterial of this paper (see Fig. S1). From our real-timeTDDFT simulations we obtain time-dependent electrondensities which we analyze further using the approachesdiscussed in the following section. Inputs and outputsfrom our simulations are available at the Materials DataFacility.
III. ELECTRON DENSITY ANALYSIS
It is a central goal of this work to quantify and analyzethe charge state of the projectile both inside and outsidethe aluminum target material. Information about theprojectile’s charge state could be estimated by analyz-ing coefficients of the KS orbitals in an atomic orbitalbasis, by projecting the KS orbitals onto atomicorbitals, or by considering probabilities that the KS or-bitals are localized within a given volume of interest. However, many of these methods have well-known short-comings, most notably that they rely directly on thesingle-particle, non-interacting KS orbitals which haveno rigorous physical meaning.Instead, a method that provides the projectile chargestate as a functional of only the total electron density—areal, physical quantity which, in principle, determines allobservables —would be preferable. Existing methodsthat extract atomic charges directly from the electrondensity typically rely on volume or charge parti-tioning, which either assume a definite boundary between z (a )0102030 x ( a ) − − − − − − − electron density (e/a ) FIG. 1. (Color online.) Illustration of our simulation cell at a time of 0.75 fs after a proton with a velocity of 1 at. u. impactsa 4-layer thick aluminum sheet. A slice of the electron density is overlaid with projected positions of aluminum atoms in grayand the projectile along with its trajectory in red. Electron density has been emitted into the vacuum and captured by theprojectile. Dashed yellow lines indicate the aluminum-vacuum boundaries used for calculating the sheet’s dipole moment, andthe 5 a projectile radius used to compare our charge capture method with volume partitioning is indicated by the red dashedcircle (see Section III). the captured electrons and nearby free electrons or ne-glect free electrons altogether. In this context, we ap-ply the DDEC6 charge partitioning method in thiswork to calculate the effective projectile charge state asthe proton traverses the aluminum sheet. However, thistechnique is not applicable in the presence of free elec-trons and, hence, we find that it overestimates the num-ber of electrons bound to the projectile once it emergesfrom the material into the vacuum containing emittedelectrons (see Fig. S6 of the Supplemental Material).Since the electron distribution emitted into vacuum,not including those captured by the projectile, shouldnot be assigned to any atoms, a different method isrequired to calculate charge captured by the projectileafter traversing the material. Once the projectile hasleft the target material, a common strategy simply inte-grates the electron density within a sphere centered atthe projectile. However, the radius chosen forthis sphere (for instance, 5 a in Ref. 55) defines an ar-tificial, discrete boundary between electrons captured bythe projectile and free electrons in the vacuum which, de-pending on the radius chosen and the occupied projectileorbitals, could either falsely include emitted electrons orfalsely exclude higher energy captured electrons.Instead, we present a novel, physically motivatedmethod of calculating the charge captured by a protonand also the orbital distribution of the captured chargeas a functional of the electron density. We first obtain theradial distribution n ( r, t ) of the electron density aroundthe projectile, in units of e/a , by integrating the electrondensity n ( r , t ) over spherical shells S , n ( r, t ) = 1∆ r (cid:90) S ( r, R ( t )) n ( r , t ) dr . (1) r (a )0 . . . . . . . . r a d i a l d i s tr i bu t i o n ( e /a ) simulation data1 s +uniform1 s +2 s +uniform1 s +2 s +2 p +uniform1 s contribution2 s contribution2 p contribution FIG. 2. (Color online.) Curve fitting used to determine theorbital occupations of the proton projectile and to calculatethe charge it captured. Red circles show simulation results forthe radial distribution of the electron density around the pro-jectile (1 at. u. of velocity) 1.2 fs after it exits from a 4-layeraluminum sheet. Least squares fits to this data using differentanalytic hydrogen orbitals and a uniform background chargeare shown in green, blue, and purple. Radial distributionsfor analytic 1 s , 2 s , and 2 p hydrogen orbitals are shown ingray; they are scaled by their respective occupations as ob-tained from the 1 s +2 s +2 p +uniform fit, which describes thesimulation results very well. Here, ∆ r is the smallest real-space grid spacing in eachsimulation and S ( r, R ( t )) is the spherical shell of thick-ness ∆ r and radius r , centered at the projectile’s position R ( t ). Integrating n ( r, t ) from Eq. (1) again over r wouldgive the number of electrons within a sphere around theprojectile. Using the electron densities from our simula-tions, represented on a real-space grid, we evaluate thefollowing discrete version of Eq. (1) to compute n ( r, t ): n ( r i , t ) = 1∆ r (cid:88) r i − < | r − R ( t ) |≤ r i n ( r , t ) ∆ V (2)where r i = i ∆ r ranges from 0 to 10 a and ∆ V is thevolume of each grid cell. We then fit the resulting ra-dial distribution to a linear combination of analytic ra-dial distributions of hydrogen orbitals (1 s , 2 s , and 2 p )and a uniform background density to account for nearbyfree electrons. The resulting fits capture the numericallycalculated radial distributions extremely well (see Fig. 2),with R values generally above 0.9. In the SupplementalMaterial, we show that the DFT orbitals of an isolatedH + ion closely match the analytical hydrogen orbitals(see Fig. S4) and that including even higher energy or-bitals has a negligible effect on the orbital occupationsand the total charge transfer (see Fig. S5).Finally, in order to analyze plasma oscillations inducedin the aluminum sheet, we compute the out-of-planedipole moment d ( t ) = (cid:90) V z n ( r , t ) dr (3)where ˆ z is normal to the aluminum surface and V is thevolume occupied by the sheet. In order to account for theelectron density extending from the aluminum surface,we include within V the region within 11 a of the alu-minum surface atoms (see yellow dashed lines in Fig. 1).Using this cutoff gives less than 6 × − electrons in thevacuum at the start of each calculation. IV. RESULTS AND DISCUSSIONA. Pre-equilibrium electronic stopping
A fast projectile impacting a bulk target or sufficientlythick foil reaches an equilibrium charge state within a fewfs and within the first few nm of material traversed.
The projectile then experiences equilibrium electronicstopping, and as a result its velocity decreases over atime scale much longer than the initial charge equilibra-tion. As the projectile slows down, its equilibrium chargestate also evolves as a function of velocity. However, thesituation is more complicated for thin target materials.Upon approaching and entering the surface, an ion dy-namically captures or loses electrons, leading to energytransfer dynamics which have been detected as a depen-dence of stopping power in thin foils on the initial chargestate of the projectile. This pre-equilibrium behavior inthe projectile charge and stopping power within the ma-terial surface would then also influence surface processessuch as electron emission. . . . . . . . s t o pp i n g p o w e r ( H a/a ) SRIMbulksheet . . . . in orange and SRIM data in blue). Inset: Elec-tronic stopping of a 25 keV proton (velocity of 1 at. u.) fordifferent aluminum sheet thicknesses. Dashed lines indicatebulk values from TDDFT (orange) and SRIM (blue). Foreach sheet, average stopping is computed across the two mid-dle layers as the most bulk-like region. Our results show clear pre-equilibrium effects in theenergy transferred from the projectile to the host mate-rial from the comparison of the stopping power in thinaluminum sheets to bulk aluminum: Figure 3 shows 13 –25 % greater stopping, depending on projectile velocityand sheet thickness, for H + in the aluminum sheets com-pared to theoretical and empirically fitted results forH in bulk aluminum. Since the stopping powers of Hand H + in bulk aluminum quickly converge toward thesame value, this comparison isolates surface effects inthe host material. Furthermore, the inset of Figure 3shows that the stopping power varies with sheet thick-ness and does not approach the bulk value even for a12-layer sheet, the thickest one we simulated.In order to explain why pre-equilibrium stopping islarger than bulk stopping, we first analyze the projec-tile charge state dynamics. This allows us to compare toanalytic models that predict that stopping power variesquadratically with projectile charge. Interestingly,we find from Figure 4 that the instantaneous projectilecharge state as calculated by the DDEC6 method isactually anticorrelated with the instantaneous stoppingpower: The proton experiences enhanced stopping de-spite lower effective charge within the first few atomiclayers, and the fluctuations arising from lattice periodic-ity are out of phase, with local maxima in stopping poweraligned with local minima in effective charge state andvice versa. Thus, the dynamical behavior of the projec-tile charge state appears incompatible with equilibriumstopping power models and does not explain the pre-equilibrium stopping behavior near the material’s sur- )0 . . . . . s t o pp i n g p o w e r ( H a/a ) . . . . . p r o j e c t il e c h a r g e s t a t e ( e ) FIG. 4. (Color online.) Instantaneous electronic stopping(red) and effective projectile charge state from DDEC6 analysis (blue) for a proton projectile with 1 at. u. of velocitytraversing aluminum sheets that are 4, 6, and 8 layers thick(solid, dashed, dotted, respectively). The entrance layer ofaluminum atoms is located 0 a .FIG. 5. (Color online.) Average electron density locatedwithin a 1 a radius of the proton’s trajectory for the ground-state initial condition (dashed black) and when the proton,approaching with 1 at. u. of velocity, is 3.5 a away from im-pacting a 4-layer aluminum sheet (solid red). Dotted blueshows the difference between the average densities just beforeimpact and in the ground state. The positions of the alu-minum atoms (proton projectile) are indicated in gray (gold). face. However, we also note that charge partitioningschemes may not be capable of accurately resolving therelatively small changes in the charge state occurringhere, as the captured electron density is superimposedwith the projectile’s wake in the host material.Another potential source of enhanced pre-equilibriumstopping, according to analytic models, would behigher electron density near the entrance surface. Weindeed observe such higher electron density, arising fromthe polarization induced on the aluminum sheet during hω (eV)0 . . . . . . . d i p o l e m o m e n t( e a ) bulksurface 4 layers6 layers8 layersFIG. 6. (Color online.) Fourier transform of the time-dependent out-of-plane dipole moment in the aluminumsheets after impact by a proton with 1 at. u. of velocity. Onlydata at least 0.7 fs after the projectile exits the material wasanalyzed in order to isolate plasma oscillations and excludecontributions from emitted electrons. The dashed (dotted)lines indicate the theoretical (experimental) energies of thebulk and surface plasmons. the proton’s approach (see Fig. 5). However, we foundthat before impact, this polarization is highly localizedto the very surface, without extending into the first twoatomic layers where the enhanced stopping is observed.Once the projectile enters the material, it again becomesimpossible to disentangle its wake from its captured elec-trons. Therefore, we cannot definitively attribute theenhanced surface stopping power to surface polarization.Finally, the higher stopping may be explained by sur-face or confinement effects of thin sheets: Aluminum sur-faces have an additional excitation channel in the formof surface plasmons. In addition, the plasmonic behaviorof atomically thin metal films has been shown to deviatefrom bulk and predicted to support multiple surface,subsurface, and bulk plasmon modes. Surface plasmonmodes have indeed been predicted to become increasinglyimportant and lead to higher stopping power for incidentelectrons as film thickness is reduced. To investigatethis mechanism, we performed a Fourier analysis of thetime-dependent out-of-plane dipole moment, Eq. (3), andwe indeed find indications for plasmon modes located be-tween the bulk and surface plasmon energies (see Fig. 6).While our frequency resolution is fairly low due to thefew-fs time-scale of our simulations, Figure 6 shows thatthe data for the 4-layer sheet hints at the possibility oftwo distinct plasmon peaks. Future studies with signif-icantly longer total propagation time would be neededto unambiguously distinguish specific bulk and surfacemodes. Nonetheless, surface plasmon excitations and/ornon-bulk plasmonic behavior represent plausible mecha-nisms for enhanced electronic stopping near an aluminumsurface.We also note that the de Broglie wavelength of elec- . . . . e l e c t r o n s (a) 0 . . . . . . . . total1 s s + 2 p FIG. 7. (Color online.) (a) Total charge captured by the pro-jectile, decomposed into occupations of the 1 s and 2 s + 2 p hy-drogen orbitals, after a proton with 1 at. u. of velocity impactsaluminum sheets of different thickness. Error bars indicatestandard deviations of these time-dependent values. (b) Thesame quantities are shown as a function of projectile velocityfor protons impacting 4-layer aluminum. trons at the Fermi surface in aluminum (6.8 a ) is com-parable to the thickness of the aluminum sheets (15 –92 a ), suggesting that quantum confinement may affectthe thinner sheets studied. However, since the instanta-neous stopping power remains almost identical for sheetswith different thicknesses until the projectile reaches theexit surface (see Fig. 4), we conclude that any quantumconfinement effects do not significantly influence elec-tronic stopping power in this system. B. Projectile charge capture
Experimental studies often infer information aboutpre-equilibrium behavior from charge state distribu-tions of the projectile after transmission through afoil.
Thus, we calculated the number of elec-trons captured by the projectile after emerging from thealuminum sheets by analyzing the electron density as de-scribed in Section III. Our method also provides infor-mation about the sub-fs real-time dynamics of orbitaloccupations of the captured electrons.First, we found that charge captured by the projectileremains nearly constant as a function of aluminum sheetthickness (see Fig. 7(a)), but decreases with higher pro-ton velocity across the entire velocity range considered(see Fig. 7(b)). The majority of the captured electrondensity occupies the 1 s orbital, with a smaller portionoccupying the 2 s and 2 p orbitals. Figure 7(b) also showsthat the 2 s and 2 p orbital occupation is largely inde-pendent of the velocity and the change in total capturedcharge can be attributed to the 1 s shell of the projectile.Hence, while slower projectiles capture more charge fromthe target and are close to their 1 s ground state, fasterprojectiles capture less charge but any captured electronsare more likely to occupy an excited state. .
50 0 .
75 1 .
00 1 .
25 1 . . . . . . . c a p t u r e d e l e c t r o n s insideafter exitstrippedFIG. 8. (Color online.) Average number of electrons cap-tured by the proton projectile inside (blue triangles) and afterexiting (purple squares) a 4-layer aluminum sheet. The dif-ference between these two quantities, the number of electronsstripped at the exit-side surface, is shown in orange exes. Er-ror bars indicate standard deviations of the time-dependentquantities. 0 . . . . . . . . . e l e c t r o n s c a p t u r e d s s +2 p total5a FIG. 9. (Color online.) Time-dependent hydrogen orbitaloccupations after a proton projectile with 1 at. u. of velocitytraverses a 4-layer aluminum sheet. The number of electronswithin 5 a of the projectile is included for reference. The lack of dependence of charge capture on sheetthickness (see Fig. 7(a)) is surprising given the pre-equilibrium stopping presented in Section IV A. In par-ticular, this behavior is a departure from experimentsemploying heavier/higher charge projectiles which ob-served that charge loss distributions depend on foilthickness.
Thus, our results indicate that light ionscan experience pre-equilibrium stopping power even afterprojectile charge equilibration. This conclusion is furthersupported by the finding that the projectile’s charge al-most reaches equilibrium even within the thinnest sheetstudied here (see Fig. 4) despite pre-equilibrium stoppingpower even for the thickest sheets (see inset of Fig. 3).Furthermore, this supports our interpretation that largerpre-equilibrium stopping for light ions is related to thetarget material, via its surface plasmons, rather than be-ing a property of the projectile charge state.We also note another interesting surface effect: Theprojectile emerges with a higher exit charge state thanits effective charge state within the sample. Our datashows this for a proton with a velocity of 1 at. u., whichequilibrates to a charge of about 0 . within the 8-layersheet (see Fig. 4), but retains only about 0 . exits with a charge of 0 . We find that the number ofelectrons that are stripped at the surface increases withproton velocity (see Fig. 8). These discrepancies betweenthe projectile’s charge state within the material and af-ter exiting into the vacuum lead us to advise caution indrawing conclusions about pre-equilibrium behavior frommeasurements of only the projectile’s exit charge state.Finally, our analysis of the femtosecond dynam-ics of projectile orbital occupations shows that for v (cid:38) .
75 at. u., the captured electrons fluctuate betweenthe 1 s and 2 s /2 p states. We show this explicitly for v = 1 . v = 1 .
25 at. u. and 1 . s occupation number features a strongpeak at a frequency ranging from (cid:126) ω = 10 . ± . . ± . n = 1 and n = 2hydrogen orbitals, the fluctuations with these frequenciessuggest a superposition of these two states. We also notethat for a projectile velocity of v = 0 .
75 at. u., the ampli-tude of the oscillations in the orbital occupation dynam-ics becomes very small, making them hard to interpret,and at lower velocities the oscillations disappear (see Fig.S8 of the supplemental material). In this low velocityregime, occupation of the 1 s orbital dominates. Thisagain indicates that slower projectiles are much closer totheir electronic ground state when exiting the aluminumsheet. V. CONCLUSIONS
Our first-principles simulations of proton-irradiatedaluminum sheets revealed detailed information aboutpre-equilibrium stopping power near a metal’s surface.We found higher average stopping power in the aluminumsheets compared to bulk aluminum and a pronouncedstopping power enhancement within the entrance layers.These deviations from bulk behavior are not adequatelyexplained by pre-equilibrium projectile charge, surfacepolarization, or quantum confinement; the most viable mechanism for surface stopping enhancement is surfaceplasmon excitations.We also presented a novel technique based on analyti-cal hydrogen orbitals to extract from the electron densitynot only the exit charge state of the projectile, but alsothe orbital occupations of the captured electrons. Theelectrons captured by the proton predominantly occupythe 1 s orbital, though for higher velocities the projec-tile emerges in a superposition of 1 s and 2 s /2 p orbitals.Moreover, the projectile’s exit charge state differs fromits equilibrium charge within the material, and the num-ber of electrons stripped from the projectile as it emergesfrom the exit-side surface increases with projectile veloc-ity.This work provides new details about pre-equilibriumbehavior in ion-irradiated surfaces and thin materials, of-fering unprecedented insight into the few-fs dynamics ofelectronic stopping power and projectile charge equilibra-tion. This study also has broad implications for advanc-ing radiation hardness and ion beam techniques for imag-ing, defect engineering, and nanostructuring, as theseapplications are chiefly concerned with energy deposi-tion near material surfaces. The electron density analysisframework developed in this work lays the foundation forfurther computational research on charge dynamics nearion-irradiated surfaces and 2D materials. ACKNOWLEDGMENTS
This material is based upon work supported by theNational Science Foundation under Grant No. OAC-1740219. Support from the IAEA F11020 CRP ”IonBeam Induced Spatio-temporal Structural Evolution ofMaterials: Accelerators for a New Technology Era” isgratefully acknowledged. This research is part of the BlueWaters sustained-petascale computing project, which issupported by the National Science Foundation (awardsOCI-0725070 and ACI-1238993) and the state of Illinois.Blue Waters is a joint effort of the University of Illi-nois at Urbana-Champaign and its National Center forSupercomputing Applications. An award of computertime was provided by the Innovative and Novel Compu-tational Impact on Theory and Experiment (INCITE)program. This research used resources of the ArgonneLeadership Computing Facility, which is a DOE Officeof Science User Facility supported under Contract DE-AC02-06CH11357. This work made use of the IllinoisCampus Cluster, a computing resource that is operatedby the Illinois Campus Cluster Program (ICCP) in con-junction with the National Center for SupercomputingApplications (NCSA) and which is supported by fundsfrom the University of Illinois at Urbana-Champaign.
VI. ACCURACY OF TIME EVOLUTION
Significant numerical error may accumulate when in-tegrating the time-dependent Kohn Sham (TDKS) equa-tions for thousands of time steps within the large su-percells involved in this study. To evaluate the accu-racy of the time propagation, we measured the error innet charge when evolving a ground state aluminum sheetwith no excitation (i.e., without a projectile). Ideally,such a test simulation should produce a static electrondensity since the initial wave function is an eigenstate ofthe time-independent Hamiltonian.We found that the popular fourth-order Runge-Kutta(RK4) propagator accumulates significant numericalerrors in net charge (see Fig. S1(a)), even when usinga small time step of dt RK4 = 0 .
35 as. This time stepwas shown in Ref. 26 to be sufficiently accurate in stop-ping power calculations for light ions traversing bulk alu-minum. However, accurately simulating an ion’s interac-tion with a surface requires much larger supercells con-taining sufficient vacuum to mitigate finite-size effects.Unfortunately, a non-negligible portion of the large er-ror in net charge manifests as a fictitious density in thevacuum region (see Fig. S1(b)), which would distort cal-culations of electron dynamics outside the material.In contrast, we found that the enforced time-reversalsymmetry (ETRS) propagator is far more accuratethan RK4 without increasing computational cost. Sinceeach step of the ETRS implementation involves twelve Hψ evaluations compared to RK4’s four Hψ evaluations,we retain a similar wall time per simulation time bychoosing a time step of dt ETRS = 3 dt RK4 = 1 .
04 as forETRS simulations. Even with this larger time step, wefind that ETRS still produces much smaller numericalerrors and essentially preserves the ground state densityduring the test simulations (see Fig. S1). ETRS has alsobeen shown to outperform more advanced Runge-Kuttamethods when applied to very similar TDKS problems. Thus, we used the ETRS integrator with a time step of1.04 as to produce almost all of the results reported inthe main text. The only exception is the stopping powerdata point for the 10-layer sheet, which was calculatedusing RK4 and a time step of 0.35 as; we find that thischoice does not influence stopping power, and we do notanalyze the subsequent charge dynamics for this case.During the simulations for which the electron density isanalyzed in the main text, the error in net charge neverexceeded 0.04 electrons.Finally, we confirmed the assumption that atomic mo-tion is negligible over the fs time-scale of the simulationsby comparing results from two simulations of a 25 keV(1 at. u. of velocity) H + projectile impacting a 2-bilayeraluminum sheet: one with fixed ionic velocities and onewith true Ehrenfest molecular dynamics. Allowing thenuclei to respond to Hellmann-Feynman forces changesthe average stopping power by less than 0.02% and thetotal charge captured by the projectile by less than 1%. − time (fs)10 − − − e rr o r i n e l e c t r o nnu m b e r (a) RK4ETRS0 25 50 75 100 125 150 z (a )10 − − e l e c t r o n s /a (b) RK4ETRSFIG. S1. (Color online.) (a) Error in net charge over timewhen propagating a ground state aluminum sheet (300 atoms)with no excitation using the ETRS and RK4 integrators withtime steps of 1.04 and 0.35 as, respectively. (b) Electron den-sity profiles from the same test simulations after propagatingfor 2.9 fs. The density profile of the initial ground state isshown in black. RK4 and ETRS accumulate about 0.5 and0.005 fictitious electrons in the vacuum region, respectively. VII. CONVERGENCE OF SUPERCELLDIMENSIONS
As shown in Fig. S2, we find that electronic stoppingpower does not depend strongly on the vacuum length,which we define as the difference between the length ofthe supercell and the thickness of the material (not ac-counting for finite surface widths). For channeling pro-tons with 1 at. u. of velocity, the average stopping powervaries by less than 0.2% within 4-layer aluminum sheetsamong the different vacuum lengths tested.To ensure convergence with respect to lateral materialdimensions, we simulated protons with 1 at. u. of veloc-ity irradiating 4-layer aluminum sheets of area 30 . a (128 atoms), 38 . a (200 atoms), and 45 . a (288atoms). Because of the high computational cost of thelargest sheet, the vacuum length was reduced to 70 a −
10 0 10 20 30 40projectile position (a )0 . . . . . s t o pp i n g p o w e r ( H a/a )
90 a vacuum140 a vacuum240 a vacuumFIG. S2. (Color online.) Instantaneous stopping power as achanneling proton with 1 at. u. of velocity impacts a 4-layeraluminum sheet is not sensitive to the vacuum length. Theentrance layer of aluminum atoms is located at 0 a .0 10 20 30projectile position (a )0 . . . . . s t o pp i n g p o w e r ( H a/a )
128 atoms200 atoms288 atomsFIG. S3. (Color online.) Instantaneous stopping power as aproton with 1 at. u. of velocity traverses a 4-layer aluminumsheet does not change as the lateral dimensions of the sheetare varied. The entrance layer of aluminum atoms is locatedat 0 a . for these tests. Also, instead of ETRS with a time stepof 0.043 at. u. as decided in Section VI, these aluminumcalculations used the RK4 time-stepper with a time stepof 0.014 at. u. We do not expect these parameter choicesto significantly influence the converged supercell dimen-sions. We found very good convergence of the instanta-neous stopping power with respect to the area of the ma-terial (see Fig. S3), with average stopping power withinthe material differing by only 0.5% among the differentsize aluminum samples tested. r ( a )0 . . . . . . P ( r )( e /a ) analyticalcalculatedFIG. S4. (Color online.) Radial distributions of the analyticalhydrogen orbitals and the ground-sate Kohn-Sham orbitals ofan isolated H + ion. VIII. CHARGE CAPTURE
First, we validated the applicability of our novel ap-proach for calculating the projectile’s exit charge state(see Section III) by examining the hydrogen orbitals cal-culated by DFT. We found that the radial distributionsof the densities corresponding to the ground-state Kohn-Sham orbitals φ j ( r ) of an isolated H + ion, i.e., P j ( r ) = (cid:90) π (cid:90) π | φ j ( r ) | r sin θ dθ dφ, match extremely well with the analytical hydrogen or-bitals (see Fig. S4). This suggests that the electrons cap-tured by the proton projectile in our TDDFT simulationscan be adequately described with analytical hydrogen or-bitals. Also, for heavier ions where analytical orbitalsare not available, our novel approach can still be appliedusing calculated radial distributions of the isolated pro-jectile.We also found that including higher energy orbitals(above 2 p ) in the orbital fitting does not significantlyalter the orbital occupations or the total charge capture(see Fig. S5). Therefore, we concluded that contributionsfrom orbitals above n = 2 are negligible and we consid-ered only the 1 s , 2 s , and 2 p orbitals when computing thecharge capture results presented in the main text.While the simple volume partitioning method of in-tegrating the electron density within a 5 a radius ofthe projectile produces qualitatively similar results forthe total charge captured as our novel approach, theDDEC6 charge partitioning method overestimates thetotal charge captured compared to the other two tech-niques (see Fig. S6). We attribute this overestimation tothis method’s inapplicability in the presence of a free elec-tron distribution, which causes erroneous assignment ofnearby emitted electrons to the projectile. We found that0 s s p s p d highest energy orbital included0 . . . . c a p t u r e d e l e c t r o n s total1 s s p s p d FIG. S5. (Color online.) Projectile orbital occupations andtotal captured electrons obtained from the novel orbital fittingtechnique (see Section III) as the number of included orbitalsvaries. The time-dependent occupations are averaged overtime after a proton with velocity 1 at. u traverses a 4-layeraluminum sheet, and their standard deviations are illustratedby error bars. The average number of electrons within 5 a of the projectile, shaded to indicate standard deviation, isincluded for reference.0 . . . . . . . . . e l e c t r o n s c a p t u r e d s s +2 p total 5a DDEC6
FIG. S6. (Color online.) Charge captured by a H + pro-jectile with 1.0 at. u. of velocity after traversing an 8-layeraluminum sheet. Results computed by DDEC6 analysis are compared to the novel fitting technique (see Section III)and the number of electrons within 5 a of the projectile. the DDEC6 results are particularly inaccurate when theprojectile is still near the exit-side surface, but becomequalitatively similar to the other methods over time.Finally, we present additional results for projectile or-bital occupation dynamics. Figure S7 illustrates similaroscillatory behavior at v = 1 . v = 1 . ∼ .
75 at. u., these oscillations disappear and oc-cupation of the 1 s orbital dominates (see Fig. S8). ∗ [email protected] J. Notte, B. Ward, N. Economou, R. Hill, R. Percival,L. Farkas, and S. McVey, AIP Conference Proceedings , 489 (2007). I. Utke, P. Hoffmann, and J. Melngailis, Journal ofVacuum Science & Technology B: Microelectronics andNanometer Structures , 1197 (2008). R. Ritter, R. A. Wilhelm, M. Stger-Pollach, R. Heller,A. Mcklich, U. Werner, H. Vieker, A. Beyer, S. Facsko,A. Glzhuser, and F. Aumayr, Applied Physics Letters ,063112 (2013). G. Hlawacek, V. Veligura, R. van Gastel, andB. Poelsema, J. Vac. Sci. Technol. B , 020801 (2014),10.1116/1.4863676. D. S. Fox, Y. Zhou, P. Maguire, A. ONeill, C. Coilein,R. Gatensby, A. M. Glushenkov, T. Tao, G. S. Duesberg,I. V. Shvets, M. Abid, M. Abid, H.-C. Wu, Y. Chen, J. N.Coleman, J. F. Donegan, and H. Zhang, Nano Letters ,5307 (2015). Z. Li and F. Chen, Applied Physics Reviews , 011103(2017). H. Bethe, Annalen der Physik , 325 (1930). . . . . . . . e l e c t r o n s c a p t u r e d s s +2 p total5a FIG. S7. (Color online.) Time-dependent hydrogen orbitaloccupations after a H + projectile with a velocity of 1.5 at. u.traverses a 4-layer aluminum sheet. The number of electronswithin 5 a of the projectile is included for reference.0 .
75 1 .
00 1 .
25 1 .
50 1 . . . . . . e l e c t r o n s c a p t u r e d s s +2 p total5a FIG. S8. (Color online.) Time-dependent hydrogen orbitaloccupations after a H + projectile with a velocity of 0.5 at. u.traverses a 4-layer aluminum sheet. The number of electronswithin 5 a of the projectile is included for reference. P. Sigmund,
Particle penetration and radiation effects:general aspects and stopping of swift point charges ,Springer series in solid-state science No. 151 (Springer,2008). E. Fermi and E. Teller, Physical Review , 399 (1947). J. Lindhard and A. Winther, Matematisk-fysiske Med-delelser , 1 (1964). P. Wang, T. M. Mehlhorn, and J. J. MacFarlane, Physicsof Plasmas , 2977 (1998). A. A. Correa, Computational Materials Science , 291(2018). C.-W. Lee, J. A. Stewart, S. M. Foiles, R. Dingreville, andA. Schleife, arXiv:2001.10162 [cond-mat] 2001.10162. M. Hattass, T. Schenkel, A. V. Hamza, A. V. Barnes,M. W. Newman, J. W. McDonald, T. R. Niedermayr, G. A.Machicoane, and D. H. Schneider, Physical Review Letters , 4795 (1999). R. A. Wilhelm, E. Gruber, R. Ritter, R. Heller, S. Fac-sko, and F. Aumayr, Physical Review Letters , 153201(2014). T. Schenkel, M. A. Briere, H. Schmidt-Bcking, K. Bethge,and D. H. Schneider, Physical Review Letters , 2481(1997). H. L. Sun, Y. C. Yu, E. K. Lin, C. W. Wang, J. L. Dug-gan, A. R. Azordegan, F. D. McDaniel, and G. Lapicki,Physical Review A , 4190 (1996). W. Brandt, R. Laubert, M. Mourino, andA. Schwarzschild, Physical Review Letters , 358(1973). R. A. Wilhelm, E. Gruber, J. Schwestka, R. Kozubek, T. I.Madeira, J. P. Marques, J. Kobus, A. V. Krasheninnikov,M. Schleberger, and F. Aumayr, Physical Review Letters , 103401 (2017). R. A. Wilhelm, E. Gruber, J. Schwestka, R. Heller,S. Fascko, and F. Aumayr, Applied Sciences , 1050(2018). P. Koschar, K. Kroneberger, A. Clouvas, M. Burkhard,W. Meckbach, O. Heil, J. Kemmler, H. Rothard, K. O.Groeneveld, R. Schramm, and H.-D. Betz, Physical Re-view A , 3632 (1989). P. Mertens and T. Krist, Nuclear Instruments and Methodsin Physics Research , 57 (1982). D. Semrad, P. Mertens, and P. Bauer, Nuclear Instru-ments and Methods in Physics Research Section B: BeamInteractions with Materials and Atoms , 86 (1986). S. P. Møller, A. Csete, T. Ichioka, H. Knudsen, U. I. Ug-gerhøj, and H. H. Andersen, Physical Review Letters ,042502 (2004). A. A. Correa, J. Kohanoff, E. Artacho, D. Snchez-Portal,and A. Caro, Physical Review Letters , 213201 (2012). A. Schleife, Y. Kanai, and A. A. Correa, Physical ReviewB , 014306 (2015). E. E. Quashie, B. C. Saha, and A. A. Correa, PhysicalReview B , 155403 (2016). M. Caro, A. A. Correa, E. Artacho, and A. Caro, ScientificReports , 2618 (2017). R. Ullah, E. Artacho, and A. A. Correa, Physical ReviewLetters , 116401 (2018). A. Lim, W. Foulkes, A. Horsfield, D. Mason, A. Schleife,E. Draeger, and A. Correa, Physical Review Letters (2016), 10.1103/PhysRevLett.116.043201. D. C. Yost and Y. Kanai, Physical Review B (2016),10.1103/PhysRevB.94.115107. D. C. Yost, Y. Yao, and Y. Kanai, Physical Review B ,115134 (2017). C.-W. Lee and A. Schleife, The European Physical JournalB , 222 (2018). K. Kang, A. Kononov, C.-W. Lee, J. A. Leveillee, E. P.Shapera, X. Zhang, and A. Schleife, Computational Ma-terials Science , 207 (2019). C.-W. Lee and A. Schleife, Nano Letters , 3939 (2019). J. M. Pruneda, D. Snchez-Portal, A. Arnau, J. I. Juar-isti, and E. Artacho, Physical Review Letters , 235501(2007). F. Mao, C. Zhang, J. Dai, and F.-S. Zhang, Physical Re-view A (2014), 10.1103/PhysRevA.89.022707. C.-K. Li, F. Wang, B. Liao, X.-P. OuYang, and F.-S.Zhang, Physical Review B , 094301 (2017). M. Pealba, A. Arnau, P. M. Echenique, F. Flores, andR. H. Ritchie, Europhysics Letters (EPL) , 45 (1992). M. Quijada, A. G. Borisov, I. Nagy, R. D. Muio, and P. M. Echenique, Physical Review A , 042902 (2007). M. A. Zeb, J. Kohanoff, D. Snchez-Portal, and E. Ar-tacho, Nuclear Instruments and Methods in Physics Re-search Section B: Beam Interactions with Materials andAtoms Proceedings of the 11th Computer Simulation ofRadiation Effects in Solids (COSIRES) Conference SantaFe, New Mexico, USA, July 24-29, 2012, , 59 (2013). E. Runge and E. K. U. Gross, Physical Review Letters ,997 (1984). M. Marques and E. Gross, Annual Review of PhysicalChemistry , 427 (2004). M. Marques,
Time-Dependent Density Functional Theory (Springer Science & Business Media). C. A. Ullrich,
Time-Dependent Density-Functional The-ory: Concepts and Applications (Oxford University Press,2011). C. A. Ullrich and Z.-h. Yang, Brazilian Journal of Physics , 154 (2014). A. Schleife, E. W. Draeger, Y. Kanai, and A. A. Correa,The Journal of Chemical Physics , 22A546 (2012). A. Schleife, E. W. Draeger, V. Anisimov, A. A. Correa,and Y. Kanai, Comput. Sci. Eng. , 54 (2014), specialIssue on ”Leadership Computing”. F. Gygi, IBM J. Res. Dev. , 137 (2008). E. W. Draeger and F. Gygi, “Qbox code, qb@ll version,”(2017), Lawrence Livermore National Laboratory. N. Troullier and J. L. Martins, Physical Review B , 1993(1991). D. Vanderbilt, Physical Review B , 8412 (1985). A. Zangwill and P. Soven, Physical Review Letters , 204(1980). A. Zangwill and P. Soven, Physical Review B , 4121(1981). S. Zhao, W. Kang, J. Xue, X. Zhang, and P. Zhang, Jour-nal of Physics: Condensed Matter , 025401 (2015). A. Castro, M. A. L. Marques, and A. Rubio, The Journalof Chemical Physics , 3425 (2004). E. W. Draeger, X. Andrade, J. A. Gunnels, A. Bhatele,A. Schleife, and A. A. Correa, Journal of Parallel andDistributed Computing , 205 (2017). B. Blaiszik, K. Chard, J. Pruyne, R. Ananthakrishnan,S. Tuecke, and I. Foster, JOM , 2045 (2016). A. Kononov and A. Schleife, “Dataset for ‘pre-equilibrium stopping and charge capture in proton-irradiated aluminum sheets’,” (2020). R. S. Mulliken, The Journal of Chemical Physics , 1833(1955). O. Knospe, J. Jellinek, U. Saalmann, and R. Schmidt, TheEuropean Physical Journal D - Atomic, Molecular, Opticaland Plasma Physics , 1 (1999). O. Knospe, J. Jellinek, U. Saalmann, and R. Schmidt,Physical Review A , 022715 (2000). C. M. Isborn, X. Li, and J. C. Tully, The Journal of Chem-ical Physics , 134307 (2007). Y. Miyamoto and H. Zhang, Physical Review B , 161402(2008). P.-G. Reinhard and E. Suraud, Journal of Cluster Science , 239 (1999). C. A. Ullrich, Journal of Molecular Structure:THEOCHEM , 315 (2000). F. Wang, X. C. Xu, X. H. Hong, J. Wang, and B. C. Gou,Physics Letters A , 3290 (2011). F. Wang, X. H. Hong, J. Wang, B. C. Gou, and J. G.Wang, Physics Letters A , 469 (2012). C.-Z. Gao, J. Wang, and F.-S. Zhang, Chemical Physics , 9 (2013). R. F. W. Bader,
Atoms in Molecules: A Quantum Theory (Clarendon Press, 1994). T. A. Manz and N. G. Limas, RSC Advances , 47771(2016). N. Gabaldon Limas and T. A. Manz, RSC Advances ,45727 (2016). P. Reinhard, E. Suraud, and C. Ullrich, The EuropeanPhysical Journal D - Atomic, Molecular, Optical andPlasma Physics , 303 (1998). Y. Miyamoto and H. Zhang, Physical Review B , 045433(2008). G. Avendano-Franco, B. Piraux, M. Grning, andX. Gonze, Theoretical Chemistry Accounts , 1289(2012). F. Mao, C. Zhang, C.-Z. Gao, J. Dai, and F.-S. Zhang,Journal of Physics: Condensed Matter , 085402 (2014). A. Ojanper, A. V. Krasheninnikov, and M. Puska, Physi-cal Review B , 035120 (2014). R. Herrmann, C. L. Cocke, J. Ullrich, S. Hagmann,M. Stoeckli, and H. Schmidt-Boecking, Physical ReviewA , 1435 (1994). N. Bohr,
The Penetration Atomic Particles Through Mat-ter (The Royal Danish Academy of Science, 1948). H. Ogawa, I. Katayama, I. Sugai, Y. Haruyama, M. Saito,K. Yoshida, M. Tosaki, and H. Ikegami, Nuclear Instru-ments and Methods in Physics Research Section B: BeamInteractions with Materials and Atoms , 80 (1993). J. F. Ziegler, M. D. Ziegler, and J. P. Biersack, NuclearInstruments and Methods in Physics Research Section B:Beam Interactions with Materials and Atoms 19th Interna-tional Conference on Ion Beam Analysis, , 1818 (2010). D. Pines, Reviews of Modern Physics , 184 (1956). C. J. Powell and J. B. Swan, Physical Review , 869(1959). G. Grosso and G. P. Parravicini, in
Solid State Physics(Second Edition) (Academic Press, 2014) pp. 287–331. R. A. Maniyara, D. Rodrigo, R. Yu, J. Canet-Ferrer, D. S.Ghosh, R. Yongsunthon, D. E. Baker, A. Rezikyan, F. J.Garca de Abajo, and V. Pruneri, Nature Photonics ,328 (2019). K. Andersen, K. W. Jacobsen, and K. S. Thygesen, Phys-ical Review B , 245129 (2012). R. H. Ritchie, Physical Review , 874 (1957). A. Brgi, M. Gonin, M. Oetliker, P. Bochsler, J. Geiss,T. Lamy, A. Brenac, H. J. Andr, P. Roncin, H. Laurent,and M. A. Coplan, Journal of Applied Physics , 4130(1993). L. Folkerts, S. Schippers, D. M. Zehner, and F. W. Meyer,Physical Review Letters74