On the prediction of electronic stopping power from first principles

Abstract

The prediction of electronic stopping power in materials from first principles using real-time time-dependent density functional theory, as implemented in the SIESTA code, combines classical motion of the nuclei with evolution of the electronic wavefunctions according to the time-dependent Kohn-Sham equation. The wavefunctions are numerically integrated according to the Crank-Nicolson integrator algorithm. In this work a correction to the existing Crank-Nicolson algorithm was implemented in parallel, based on using the correct time-dependent Kohn-Sham equation involving the covariant derivative of the wavefunction for a moving basis, and was shown to reproduce expected experimental results including oscillations of the Kohn-Sham energy for a projectile starting to move from a stationary position, which were lacking with the original integrator. The new integrator also predicted the electronic stopping power for a proton in graphite in line with experimental results, at velocities between 0.1 a.u. and 10 a.u., when the original integrator had been demonstrated to be inaccurate at velocities above 1 a.u.. The accuracy, convergence with timestep, and parallelisation of the new integrator was examined using example simulations of a proton projectile in graphite, and two helium atoms. The new integrator is not strictly unitary, unlike the previous implementation, and therefore requires a smaller timestep for accuracy. The removal of the basis changing step, which dealt with the moving basis in the original integrator via Lowdin orthonormalisation, resulted in a significant increase in speed using the new integrator, as a result of the removal of the computationally expensive matrix inversion in the basis change. The implementation of this new integrator extends the velocity range over which electronic stopping power can be calculated using the SIESTA code in any material above 1 a.u.; this allows exploration of higher projectile energy ranges. The prediction of electronic stopping power for a proton moving in graphite was investigated using both integrators, finding a threshold velocity below which electronic stopping power falls to zero when the projectile is travelling perpendicular to the graphitic planes, in the direction in which the material has a band gap and is therefore an insulator. No such threshold was found for a projectile travelling parallel to the graphitic layers; these results were in line with previous work finding a velocity threshold for electronic stopping power in insulators and semiconductors but not in metals. The electronic stopping power at low projectile velocities in a large simulation of 1000 diamond atoms was calculated and found to be in line with experimental data, with a threshold velocity found in all three directions. The valence electron density and volume of the excited electron cloud for projectiles moving at different velocities were also compared, investigating the duration of the transient response at the start of the simulation, and the extent to which the projectile interacts with the already excited electron density in the neighbouring simulation box with periodic boundary conditions, which at high velocities begins to occur slightly before the projectile re-enters the box. The velocities and simulation sizes for which the calculation of electronic stopping power can be accurately carried out was explored. Known non-adiabatic force terms for a moving basis were derived, and partially implemented into SIESTA; the relationship between these forces and the intrinsic curvature of the manifold was demonstrated, and in addition, the interpretation of the Pulay forces as the result of the extrinsic curvature of the manifold was also shown.