The remarkable success of x-ray free-electron lasers and their ability to image biological macromolecules while outrunning secondary radiation damage due to photoelectrons, by using femtosecond pulses, raise the question of whether this can be done using pulsed high-energy electron beams. In this paper, we use excited state molecular dynamics simulations, with tabulated potentials, for rare gas solids to investigate the effect of radiation damage due to inelastic scattering (by plasmons, excitons, and heat) on the pair distribution function. We use electron energy loss spectra to characterize the electronic excitations responsible for radiation damage.

## I. INTRODUCTION

Radiation damage imposes fundamental limitations on the fidelity and resolution of imaging and diffraction methods. Recently, it has been suggested and demonstrated that damage may be overcome using femtosecond x-ray pulses to out-run the damage processes.^{1} Using these pulses, the serial crystallography approach has made it possible to solve the structure of protein micro-crystals, record the femtosecond dynamics of light-sensitive proteins in micro-crystals, and to observe protein dynamics at atomic resolution during enzyme catalysis^{2}. This ability to outrun damage is important because it circumvents the need to immobilize by freezing proteins to avoid damage in protein crystallography and so allows observations of dynamics under near-physiological conditions of temperature and chemical environment.^{3} The aim of this research is to determine if this can also be done using brief electron beam pulses for the collection of transmission electron diffraction data and, if not, to understand the severity of the damage processes that remain. As a preliminary, in this paper, we study electron-beam radiation damage in the simple Rare Gas Solids (RGSs) in order to avoid the complications of interatomic bonds and their excitations.

Differences between diffraction and damage processes using x rays and electrons have been compared by Henderson^{4} (continuous exposure) and Spence^{5} (pulsed mode). Unlike x rays, electrons are diffracted from both the nuclei (Rutherford scattering) and the electron density. While the photoelectron cross section greatly exceeds that for elastic scattering for x rays, elastic and inelastic scattering cross sections for electrons are comparable. Laser amplification is now available for hard x rays, but not for electron beams, which provide fewer electrons per pulse but have a much larger elastic cross section. The Coulomb repulsion between electrons makes focusing difficult simultaneously at the brief and high currents. Flux is the product of source brightness and emittance, with emittance proportional to the product of beam divergence and source size. For transmission Electron Diffraction (ED) from thin crystals, to avoid overlap of Bragg spots, this beam divergence must be less than the Bragg angle, which for typical photocathode emittance values results in a focused beam many micrometers wide. As a result, Ultrafast Electron Diffraction (UED) experiments are undertaken with very low fluence at the sample per shot, and samples, in the form of a thin slab, should be homogeneous over an area of several micrometers. A typical fluence^{6} (integrated flux) of ∼10^{−4} − 10^{−6} *e*/Å^{2} (or *dose* as defined in the Cryo-EM literature) is therefore much less than that used in Cryo-EM. These are the conditions for which our simulations are performed. Since the Cryo-EM *dose* of about <10 *e*/Å^{2} produces very little damage, our first conclusion is that Coulomb interactions in the beam and the resulting large beam diameters used with picosecond pulses in the UED field will produce negligible amounts of damage. Our simulations are therefore for much longer pulses than used in UED.

The simplest form of radiation damage and energy transfer from the beam to the sample occurs via direct Primary *Knock-on* Atom (PKA) elastic collisions between beam electrons and the nuclei (see Sec. II A). Other than PKA, radiation damage is caused by inelastic collisions between the sample and beam electrons, causing the sample to flow on a new potential energy surface. This change in nuclear–nuclear correlation may be captured by using a time-resolved Radial Distribution Function (trRDF), which may be computed using Molecular Dynamics (MD) simulations. In this paper, we describe a program that performs the required excited-state MD to simulate RGS irradiation for electrons (also x-ray) beam pulses. The specific aim of this paper is to use this program to compute the trRDF during irradiation by an electron beam pulse.

Concurrently, dynamical (space–time) correlation may be captured by the Van Hove function,^{7} which is directly related, by the imaginary part of its Fourier transform, to the loss or the response function. This may be directly measured using the Electron Energy Loss Spectra (EELS). EELS spectra may be recorded from nanometer-sized regions in an electron microscope simultaneously with the elastic scattering. These spectra may be interpreted in terms of a frequency and momentum dependent dielectric function,^{8}

Here, *ε*(*q*, *ω*) is the dynamic dielectric function, with $\u2212\lambda =\lambda 2\pi $, where *λ* is the wavelength (in Bohr radius) of the incident particle, the scattering vector $q=2\u2061sin\u2061\theta \u2212\lambda $, and the energy loss *E* = *ℏω*. EELS spectra for the RGS which form at low temperature have been published and interpreted, identifying the specific crystal excitations to which the beam gives up energy as it traverses a thin slab of crystal. The results are shown in Fig. 1 of the study of Bernstorff and Saile.^{9}

Two radiative loss mechanisms also captured by EELS spectra are the Bremsstrahlung and Cherenkov effects. Both of these interactions emit light that can be assumed to escape but contribute negligibly to the dynamics of the sample. This is because RGSs are transparent to visible light due to the large band gap *E*_{G}, and to x rays, created by the Bremsstrahlung because of the small sample size and cross section. The nonradiative forms of energy loss by electrons beams in matter are valence electron excitation and ionization, inner-shell excitation and ionization, plasmons, excitons, and direct phonon excitations. Although inner shell excitations transfer more energy to the sample, they are much rarer events with smaller cross sections and are ignored in this model. Direct phonon excitations are similarly ignored due to their small cross sections relative to plasmons and excitons.

## II. METHODS

### A. The model

The dominant features seen in experimental EELS spectra, consist of two types: plasmon and exciton excitation. However, we first consider the ballistic PKA events, akin to Rutherford scattering. The maximum energy transfer from the beam electron to a nucleus occurs for a head-on collision and may result in irreversible nuclear displacement if this energy transfer exceeds the displacement energy, which depends on the atomic number (it is 20 eV for copper). The effect has been studied extensively in high-voltage (e.g., 1 MeV) TEM using a continuous electron beam. The total cross section for this process (about 10^{−6} Å^{2}) is much smaller than that for elastic scattering or inelastic excitation of valence electrons—for an atomic number of 20, these are both equal to about 10^{−2} Å^{2}. The volume fraction of nuclei displaced by a head-on collision is^{10} $Cd=\sigma jT|e|$, where *σ* is the cross section for displacement of a nucleus in a head-on collision (scattering angle 180°, largest energy transfer), *j* is the beam current density, and *T* the exposure time. For light elements at 1 MeV, Cd = 0.0125 with a continuous TEM current of 0.2 *μ*A and a beam diameter of 5 *μ*m after one minute. Since typical UED average currents are typically about 1 nA (for a bunch charge of a pA at a repetition rate of 1 KHz) and beam diameters, limited by space-charge, even larger, the volume fraction of atomic displacements in UED will be far less than this, and the effect can be neglected compared to other inelastic processes at these current densities.

Plasmon excitations are the dominate form of collective electronic excitation specially for higher *Z* RGS^{9} (also for proteins in ice^{11}), and they are highly localized and damped, as suggested by their broad spectral linewidth corresponding to a lifetime of ∼100 as. The products of plasmon decay are believed to be the main cause of radiation damage in Cryo-EM, on the basis of their EELS spectra.^{11} For simulations with a time step larger than 100 as, the plasmon excitation decay is simply modeled as a re-scaling of sample’s nuclear velocities by the plasmon energy. The plasmon energy is proportional to the square root of the density of valence electrons. Although a perhaps bad approximation, the RPA dielectric function is used in conjunction with Eq. (1) to determine the plasmon excitation cross section.

Unlike plasmons, Frenkel exciton excitations in RGS are relatively long lived (∼25 ns),^{12} and are also prominent in experimental EELS spectra. In order to incorporate the collective effect on the structural dynamics, the complete static Coulomb interaction (explicit electrostatics between all nuclei and electrons) was computed for the sample. Excitons are then incorporated by adding a new set of Rydberg/Wannier states (*ϕ*_{n}(*r*)) in addition to the non-interacting Hartree–Fock (HF) states, given by the Wannier equation [Eq. (2)] with static dielectric constant *ε*(0) [*α*_{i}(0) being the static polarizability of atom *i* with Vol being the sample’s volume and *μ* being the exciton mass],

At every time step, *t*, of the stochastic simulation, all atoms are subject to the following operator:

Where {*w*^{†}, *w*} and {*s*^{†}, *s*} are creation and annihilation operators for Wannier orbitals and Slater (HF) orbitals, respectively. Both *X*(*t*) and *Y*(*t*) take random Boolean values, {0, 1}. The probability of 1 is dependent on the half-life^{12} of the Poisson process of decay (*X*(*t*)) and the excitation cross section and beam fluence at time *t* (*Y*(*t*)). Since the sample is homogeneous, excitons are supposed to be stationary, affecting only the total electrostatics and thereby dynamics.

In previous simulations of X-ray Free-Electron Lasers (XFEL) beam damage,^{13} secondary radiation damage (that is, atomic ionization due to ejected photo-electrons) was considered. This is very different in the case of electron beams, where the atomic electrons ejected by the beam mostly have energies of a few electron volts and where lower energy excitons and plasmons dominate the spectra. These low energy secondary electrons are detected and used for image-formation in the Scanning Electron Microscope (SEM). They have been detected and studied in detail in time co-incidence with the plasmons seen in EELS spectra in order to understand their origin.^{14} In this work, we have assumed that secondary electrons do not have sufficient energy to ionize neighboring RGS atoms. Due to the Ramsauer–Townsend effect, we can expect that because of the high band gap of rare gas solids, the inelastic mean free path of ejected secondary electrons will be greatly increased.^{15}

### B. Stochastic molecular dynamics

In order to study radiation damage, we wish to obtain the time-resolved structure factor and the radial distribution function (trRDF), both of which may be computed by MD simulations. MD was used to track the location of nuclei at every time step during irradiation. The MD was implemented by integrating Langevin’s equations with a heat bath parameter *γ* (in a.u.). In the limit *γ* → 0, we recover Newtonian dynamics. The primary difficulty with this implementation of MD is the need for constant modifications of the potential energy surface and the subsequent excited state dynamics as ionization events occur during the MD. These simulations should be performed on RGS clusters on time scales of typical electron pulses (∼10 ps) and length scales of typical proteins (∼100 Å).

Ideally, the *Ab Initio* MD (AIMD) simulation with dynamics based on the time-dependent Schrödinger equation would be desirable to model such a complex environment. However, due to size and time limitations, post-Hartree–Fock methods are not feasible. Similarly, dynamics based on Kohn–Sham Density Functional Theory (KS-DFT) is also not quite feasible. Although larger space–time scales may be accessed than those by post-HF methods, the desired space–time scale is still not achievable. In addition, KS-DFT as currently developed may not fully capture the dynamical electron correlations, which result in van der Waals interactions, the only form of cohesion in RGS.

Instead, the Hartree–Fock equations were solved using a STO-6G basis set for a range of inter-nuclear distances between two RGS atoms, with either or both excited (excited orbitals were solutions to Wannier’s equation). This formed a potential energy surface and was used to create a tabulated potential (implemented with resolution of 0.005 Å) on which nuclei may be propagated via classical MD. In addition to this tabulated potential, a London dispersion potential was included from the first Padé approximant of the Casimir–Polder interaction to incorporate long-range cohesion.

## III. RESULTS

A single simulation is demonstrated to show the capabilities of this program. In this simulation, a 1 MeV electron Gaussian pulse of 1 ns (0.5 ns FWHM) duration, 10 *e*/Å^{2} fluence, and 1 *μ*m spot size ($\pi 4\u2009\mu m2$ beam area) irradiated a solid Argon cluster, Ar_{4631}, centered and equilibrated at 50 K before the pulse. Visually, the results are shown on Fig. 1, and this figure compares a cluster at a pristine condition to after the pulse has left the sample. The energy of the molecular system is shown in Fig. 2 and demonstrates discrete jumps, both exciting and relaxing, in potential energy corresponding to the exciton vertical transitions [see Eq. (4)]. The difference in the total energy before and after the simulation (or pulse duration) constitutes the energy deposited in the sample (in this case 2.1 keV) or a dose (when divided by clusters mass) of 1.1 MGy.

Because of the stochastic nature of this simulation, averages over many simulations are needed to simulate experimental results. In this case, a set of 24 simulations of a 1 MeV electron pulse of 100 ps (50 ps FWHM) duration, 1 *e*/Å^{2} fluence, and 1 *μ*m spot size was used to irradiate the aforementioned cluster, Ar_{4631}. The results of the simulation-averaged atom-averaged trRDF are displayed in Fig. 3. This figure shows minimal radiation damage with this 1 *e*/Å^{2} fluence pulse. Energy calculations show an averaged deposited dose of 193 kGy.

Also investigated was the effect of the coupling, *γ*, of the sample to a Langevin thermostat from just plasmon excitations (ignoring exciton excitations). The results from varying this parameter to a 1 MeV electron pulse of 100 ps (50 ps FWHM) duration, 10 *e*/Å^{2} fluence, and 1 *μ*m spot size are shown in Fig. 4. This figure shows that weakly coupled systems, i.e., isolated clusters in vacuum, suffer from radiation damage to a higher degree than samples connected to a thermostat (which would effectively cool the sample).

## IV. CONCLUSION

This program allows the analysis of radiation damage mechanisms and excitations responsible by suppressing each of them in turn. It makes use of EELS spectra to identify important excitations for a given system. Damage mechanisms may be studied under a wide range of conditions. These include beam energy, exposure time, fluence (*dose*), and sample temperature. The program’s choice of RGS is well-suited as a first step in understanding radiation damage in soft-matter on these time and length scales. RGS are known for their small cohesive energies and thus radiation sensitivity. It is thus unsurprising that a dose of 1–5 MGy appears to cause dissociation of the cluster on these time-scales. By comparison, the maximum exposure/dose used in Cryo-EM (protein/soft-matter) samples are roughly <30 MGy.^{4} These samples in addition to being in the bulk (with neighboring unirradiated sample) have much larger cohesive energies (from covalent bonds) than these RGS clusters. Apart from covalent bonds, RGS and protein matter have similar EELS spectra^{11} in the form of excitons (albeit much less dominant) and plasmons. van der Waals interactions are important in both systems (RGS and protein crystals).

Although radiation damage from XFELs are largely irreversible, radiation damage from electron beams may be reversible. Excitons mostly undergo radiative decay, while plasmons decay into phonons producing heat. This heat may or may not be conducted away, depending on the sample geometry. Excitons are more damaging (potentially irreversible) when they are in close proximity in space and time (Fig. 1 shows an example of this effect). Distancing these electronic excitations in time (exciting them at different times and allowing them to decay) could explain the experimental results of ultrafast electron beams in paraffin showing increased resistance to damage.^{16}

The exposures required to resolve the structure of an individual protein (in single-particle mode) require at least one elastic event per non-hydrogen atom^{4} and therefore require a total exposure of ∼10^{3} *e*/Å^{2}. This is not possible in a single shot with pulsed photocathode sources due to space-charge limitations.

## SUPPLEMENTARY MATERIAL

See the supplementary material for additional specifics. The code performing all these simulations is freely available at https://github.com/jcandane.

## ACKNOWLEDGMENTS

This work was funded by NSF STC (Award No. 1231306) by the Helmholtz Association through the CFEL at DESY, the Swedish Research Council (VR), and the Swedish Foundation for International Cooperation in Research and Higher Education (STINT).