Calculation of x-ray scattering patterns from nanocrystals at high x-ray intensity

With the advent of x-ray free electron laser (XFEL) sources,¹ studies of structural determination of biomolecules^2–5 have gained a new boost. XFELs provide intense radiation of a wavelength comparable to atomic scales. The characteristics of XFEL radiation and associated sample environments have triggered the development of new data collection methods such as serial femtosecond crystallography⁶ (SFX). The ultimate goal and dream is to perform atomic resolution single particle imaging.^7–11 Sample damage by x-rays and low signal to noise ratio at high photon momentum transfer limit the resolution of structural studies on non-repetitive structures such as individual biomolecules or cells.^7,12 Therefore, at high resolution, SFX is currently still a better option to use. XFELs deliver intense femtosecond pulses that promise to yield high-resolution diffraction data of nanocrystals (∼200 nm to 2 μm in size) before the destruction of the sample by radiation damage.^13,14 In SFX, a complete data-set can be obtained by exposing thousands of randomly oriented, individual crystals of proteins to the x-ray beam.

For imaging proteins and viruses at atomic resolution, one calls for high intensity and short x-ray pulses.^7,15–19 The shortcoming of high intensities is the rapid ionization of the atoms on the few femtosecond timescale, which affects the structure of the system. This radiation induced damage changes the atomic form factors^20,21 and may induce significant atomic displacement on longer times. Finally, radiation damage changes the scattering pattern. For a comprehensive theoretical study of signal formation in an SFX experiment, one needs to simulate (i) the radiation induced dynamics of the sample and (ii) pattern formation based on the dynamics. During the past decade, several models have been developed for studying the time evolution of small and large samples irradiated by XFEL pulses.^22–30 We use XMDYN,^31,32 a Monte-Carlo molecular-dynamics based code developed by the authors. In the theoretical study presented here, we consider a micron-size crystal in a 100 nm focus beam, a scenario where a nanocrystalline sample experiences fluences as high as to be used in single particle imaging experiments. As a consequence, the x-ray fluence is non-uniform throughout the sample. This may also have its imprint in the scattering pattern. The bottleneck one faces is that it is computationally not feasible to simulate a system with realistic size using tools which are capable to follow the dynamics of each atom, required for imaging studies. Therefore, we present an approach that involves the division of a crystal into smaller units (super-cells) and the calculation of their dynamics individually using periodic boundary conditions (PBC). In order to investigate the effect of inhomogeneous spatial fluence distribution, the super-cells are subjected to different fluences. Then we combine all the super-cells to form a nano-crystal and construct the scattering pattern under the influence of uniform (within the irradiated part of the sample) and non-uniform spatial beam profiles. We study and compare these two scenarios.

XMDYN^31–33 has been originally developed for modeling finite-size systems irradiated by an XFEL pulse. It unites a Monte-Carlo description of ionizations with a classical molecular-dynamical treatment of particle dynamics. XMDYN keeps track of the configuration of the bound electrons in neutral atoms and atomic ions. These configurations change dynamically because of different atomic processes like inner and outer-shell photoionization, Auger and fluorescence decay and collisional (secondary) ionization.

In order to treat x-ray-atom interactions, XMDYN uses the XATOM^21,34 toolkit, which is an ab-initio framework based on non-relativistic quantum electrodynamics and perturbation theory. XATOM provides rates and cross-sections of x-ray-induced processes such as photoionization, Auger decay, and x-ray fluorescence. XMDYN employs XATOM data, keeps track of all the ionization events along with the electron configuration of each atom, calculates impact ionization and recombination, and follows the trajectories of all the ionized electrons and atoms solving the classical equations of motion numerically. The framework based on these microscopic processes can describe complex many-body phenomena in ionized systems such as nanoplasma formation, charge screening, thermalization of electrons through collisions and thermal emission.³² In the current study, chemical bonds between carbon atoms are not considered. This is a good approximation when the fluence is high enough to cause severe ionization in the system early in the pulse. The immediate ionization of the atoms leads to fast bond breaking that allows their exclusion in simulations.³¹ 

The dimensions of the interaction volume are defined by the intersection of the x-ray beam and the crystal, therefore, its dimensions are determined by the focal area (∼100 × 100 nm²) and the thickness of the crystal along the beam propagation direction (μm). The number of atoms within this volume is of the order of 10⁹. This number is formidably large: it is not feasible to simulate the whole system by a single XMDYN run. In order to overcome this barrier, we propose the procedure of dividing the whole crystal into smaller units. These super-cells may contain several crystallographic unit cells. We follow the dynamics within each super-cell driven by the local fluence (assumed to be uniform throughout the super-cell) individually. For this purpose, we have developed an extension to XMDYN that applies PBC^35,36 to a super-cell, accounting also for the effect of the environment surrounding it.

Within the concept of PBC, a hypothetic crystal is constructed as a periodic extension of a selected super-cell. The total Coulomb interaction energy for a super-cell includes all the interactions within the given cell as well as pair interactions when one particle is in the cell while the other is in a periodic image within the super-cell based hypothetic lattice (PBC-crystal). Formally,

where N represents the total number of particles in the super-cell, q_i is the charge of the ith particle, ε₀ is the dielectric constant, L represents the dimension of the cell (here assumed to be a cube), $nL=n1c1+n2c2+n3c3$⁠, where $c1,c2,c3$ represent basis vectors of the PBC-crystal, and n₁, n₂, n₃ are integers indexing the periodic images. Hence, $|rij+nL|$ is the distance between the ith particle in the central super-cell (n = 0) and jth particle in the super-cell indexed by n. The symbol ′ represents the exclusion of the term j = i if and only if n = 0. The summation in Eq. (1) is not only computationally very expensive because of the formally infinite sum but is also conditionally convergent which states that the result depends upon the order of summation. To overcome this problem, we follow a route used often in the literature for spatially periodic systems, the method of minimum image convention.³⁷ According to the convention: (i) when evaluating Eq. (1), we do not use the same super-cell division of the PBC crystal for all particles, but we always shift the boundaries so that the selected particle appears in the center; (ii) we consider only n = 0 terms. The former choice ensures that no jump happens in the potential energy when a particle crosses a super-cell boundary and therefore “jumps” in the evaluation from one border of the cell to the opposite. The latter is a minimum choice considering interactions between a selected particle with the closest copy of the others only. Finally, one can assemble the entire real crystal from the individually simulated super-cells to model the whole dynamics. While in this way modeling becomes feasible even without the need of super-computers, we should also note a shortcoming of the approach: we do not allow particle transport, in particular, electron transport between the super-cells. For biologically relevant light elements, Auger and secondary electrons have energies E_kin ≲ 300 eV, which yields a short mean free path in a dense environment. Therefore, such electrons may travel only to neighboring super-cells experiencing similar fluences during the irradiation, so that the effect of net transport may be negligible. On the other hand, photoelectrons have an energy almost as high as the photon energy. Hence, they are fast and have a long mean free path: they can leave super-cells located at high fluences regions and can affect super-cells at larger distances experiencing lower fluences. We will overcome this shortcoming of the model in the future.

Although during a single shot experiment the sample may undergo significant changes, the scattering patterns are static: they accumulate diffracted signal over the whole pulse. Further, the signal may contain an imprint of a spatially non-uniform intensity profile. Formally, the scattering intensity at a specific reflection described by reciprocal vector Q, including the integration over time and the subdivision of the crystal volume into super-cells according to the approach introduced in Section II B, reads

In this equation, Q is the momentum transfer, $F={Fμ}$ is the x-ray fluence distribution throughout the crystal, the index μ runs over all super-cells, and ω is the photon energy. C(Ω) is a factor depending the polarization of the x-ray pulse, and g(t) represents the normalized temporal envelope. $fX,IX,jμ$ is the atomic form factor of the jth atom of species X in the μth super-cell, $IX,jμ$ is the associated electronic configuration, $I={IX,jμ}$ denotes a global electronic configuration, $rX,jμ$ represents the position vector of the jth atom of species X in the μth super-cell, and $r={rX,jμ}$ indicates the set of all atomic positions. N_X represents the total number of atoms for species X within a super-cell. P_I,r represents the probability distribution of electronic configuration I and atomic positions r, and R_μ represents the position of the μth super-cell. The atomic form factor

includes the dispersion corrections $fX,IX,jμ ′(ω)$ and $i fX,IX,jμ ″(ω)$⁠. This dispersion correction can be neglected when the applied photon energy is high above the ionization edges, which is fulfilled in our study. Note that the summation over $Fμ$ appears inside the modulus square in Eq. (2). The scattering amplitude from the μth super-cell is proportional to the x-ray field amplitude (⁠$∝Fμ$⁠) in that super-cell. A key assumption when performing the coherent sum in Eq. (2) is that the entire crystal is illuminated coherently, a condition that is fulfilled considering realistic XFEL beam parameters and crystal sizes.

In order to construct the scattering pattern, Eq. (2) cannot be used directly as the P_I and r configuration space is too large. However, by calculating realizations of super-cell dynamics with XMDYN, a Monte-Carlo sampling of the distribution $PI,r(F,ω,t)$ represented in Eq. (2) becomes feasible. To construct the time evolution of the crystal through global configurations and to calculate patterns, we used the following strategy, implemented in the code XSINC (x-ray scattering in nano-crystals).

We discretize the fluence space and calculate many super-cell trajectories for each fluence value with XMDYN. XSINC selects randomly a trajectory for each super-cell within the crystal (a local realization), so that the corresponding fluence values are matching the best. These trajectories describe the local time evolution of the super-cells and together they form a global realization of the crystal. Then, taking into account the spatial and temporal pulse profiles, XSINC calculates the scattering amplitudes and intensities for the global configuration at different times based on the corresponding snapshots. Finally, the incoherent sum of these patterns corresponds to a time integrated pattern measured at in a single-shot experiment. In our calculation, we perform a dense sampling of the fluence space. As a consequence, two neighboring super-cells experience very similar fluence. Therefore, it is a good approximation to take into account the direct effect of the neighboring cells by applying periodic boundary conditions and this construction leads to a realistic global trajectory. In the scheme above, several parameters are convergence parameters of the method (Table I). Results are considered converged when characteristic properties of the Bragg peaks, such as the width and height of the intensity distribution in reciprocal space, converge during monotonic increase (or decrease) of the parameter. As an example, Figures 1(a) and 1(b) illustrate the convergence of the time integrated peak intensity as a function of the number of local (super-cell) realizations per fluence point for the reflection (1 1 1) for the Gaussian and flattop spatial profile cases. We note that convergence implicitly depends on the total number of different realizations used to build a global realization. Therefore, in the Gaussian case, where 350 different fluence points are used, convergence starts at a much smaller value.

In our investigations, we consider a diamond cube of a size of 1 μm. We investigate the cases of flattop and Gaussian beam profiles (Fig. 2). Other parameters of the pulses are the same in both cases: photon energy is 10 keV, total number of photons per x-ray pulse is 1 × 10¹², the temporal pulse envelope is Gaussian with a duration of 10 fs FWHM, and focus size is 100 × 100 nm² FWHM. The size of the diamond unit cell is a = b = c = 3.57 Å containing 8 carbon atoms. The parameter choices listed in Table I yield converged results.

The coherent scattering patterns depend on the presence of the atomic bound electrons as well as on the atomic positions. The XMDYN and XATOM simulations allow to analyze their change due to radiation damage for both diamond and for the isolated carbon atom cases. Radiation damage is initiated by atomic photoionization events. In case of isolated carbon atoms, Auger decays contribute approximately to the same extent to the overall ionization. At the maximum fluence in our study, ∼35% of the atoms are photoionized (Fig. 3(a)). Although the absorbed energy is 10 keV per photon, almost all of this energy is taken away from the atom by the high energy photoelectron. The picture is different when the atom is embedded in a crystal environment (Fig. 3(c)). The high energy photoelectrons stay within the medium and distribute their energy by causing further ionization via secondary ionization events. As a consequence, neutral atoms disappear early in the pulse and by the end even fully stripped carbon ions (C⁶⁺) appear. Many electrons are promoted to (quasi-)free states within the sample. This also illustrates the importance of secondary ionizations in the progress of radiation damage in a dense environment.^38–40 In the center of focus, the sample absorbs 3.5 keV energy per atom that heats up the plasma electrons beside the ionizations. Despite the high charge states, recombination remains negligible during the pulse (number of events less than 1% per atom in the simulation) due to the extreme conditions.

Figure 4 represents the mean displacement of the carbon atoms during the pulse. The average atomic displacement is much below the maximum achievable resolution, ∼1.2 Å at 10 keV, even at the highest fluence. This suggests that the patterns are affected predominantly due to the bound-electron loss through the modification of atomic scattering form factors. Despite the heavy ionization, atomic displacements remain negligible during the ultrashort pulse duration due to the highly symmetrical sample environment. We note here again that in our calculations we neglected the chemical bonds. In low fluence regions bonds may survive and stabilize the structure against the emerging Coulomb forces. As the observed displacements are far below the resolution even without any stabilization due to bonds, bondless modeling of the current scenario is applicable.

Effect of the PBC approach on the dynamics. While ionic motion is negligible during the pulse, fast photoelectrons can travel long distances. However, PBC confines all plasma electrons artificially within the supercell they have been created in. Neglecting particle transport may lead to error in (i) local plasma electron density and (ii) local energy density. Whenever a photoelectron is ejected it leaves behind a positive charge located on an ion. If we consider Coulomb interaction only, a positive space charge would build up in a central cylinder because of photoelectron escape. Photoelectron trapping within the interaction volume would start early in the pulse, at an average ion charge as low as +0.005. An analogous phenomenon was discussed for finite samples in the literature.²⁹ However, photoelectrons cause secondary ionization as well, so an atomic bound electron is promoted to a low energy continuum state. If this slow electron was created in an outer region, it can efficiently contribute to the screening of the space charge the photoelectron left behind. Based on these arguments we can conclude that (i) considering the interaction region to be neutral is a good approximation and (ii) in all regions we overestimate the energy density by confining fast photoelectrons within a supercell. Similarly, as the Coulomb forces are the driving forces of the ionic motions, we may also overestimate the atomic/ionic displacements. In our study eventually the effect on the scattering signal is relevant, as will be discussed in Sec. III C.

In this section, we analyze the changes of the Bragg peak intensity profiles in reciprocal space due to the severe radiation damage. In Figs. 5(a) and 5(b), snapshots of the 1D Bragg peak profiles in reciprocal space are depicted for the reflection Q = (1 1 1) for Gaussian and flattop spatial beam profiles, respectively. Two apparent features can be seen, valid for other reflections as well.

• The width of the Bragg peak does not change during the pulse. This is consistent with the expectation based on the negligible ion displacements: no visible Debye-Waller-like broadening occurs. However, the widths are different for the Gaussian and flattop cases. The reason is the difference between the size of the illuminated parts of the crystal. In the flattop profile case, the focus size defines strictly the region exposed. On the other hand, a Gaussian profile has no sharp edge and therefore illuminate a larger region, yielding a narrower Bragg peak and a larger effective crystal size.

• Snapshots of the Bragg peak intensities behave differently for flattop and Gaussian beams. The snapshots of the Bragg intensities depend not only on the scattering power of the sample but also on the instantaneous x-ray intensity. However, as the instantaneous x-ray intensities are equal at the same time before and after the maximum of the pulse, a direct comparison of the corresponding snapshots of the Bragg profiles reflects exclusively the effect of different damage extents. In the Gaussian profile case, these corresponding curves show small difference only, indicating that a significant contribution is coming from regions in the crystal suffering little damage (Fig. 5(a)). In contrast, applying a flattop pulse profile, the scattering pattern is formed only from extensively ionized parts of the crystal. A consequence of the loss of atomic bound-electrons is the decrease of the atomic form factors yielding significant signal drop for longer times (Fig. 5(b)).

The above findings are reflected by the time integrated signals that correspond to the situation one would encounter in an experiment (1D cut: Figs. 5(c) and 5(d); 2D cut: Fig. 6). Note that for the Gaussian spatial profile, there is only a small decrease of the signal compared to the ideal (no damage) case.

Effect of the PBC approach on the x-ray scattering patterns. In Section III B, we discussed that the PBC approximation overestimates ionization and atomic displacements, and therefore radiation damage throughout the sample. It means that the method gives an upper bound to the effect of radiation damage on the scattering patterns. A trivial lower bound is the case without any radiation damage.

We presented a methodology for the simulation of x-ray scattering patterns from serial femtosecond crystallography experiments with a high intensity x-ray beam. Our approach includes the simulation of radiation damage within the sample with the codes XMDYN and XATOM as well as the calculation of the patterns using the code XSINC. In the approach the crystal is divided into smaller units. The time evolution (the radiation damage process) of these units is calculated using periodic boundary conditions. Finally, a nanocrystal is assembled from the small units for the calculation of the patterns integrated over the pulse.

As a demonstration, we investigated spatial pulse profile effects on the Bragg peaks for a diamond nanocrystal. We found that if a Gaussian profile is used (assuming realistic XFEL parameters such as tight focus and ultrashort pulse duration), the time integrated signal intensity is reduced only by a small amount compared to the damage-free case. For a flattop profile, the decrease is much more significant. The intensity reduction is due to the change of the form factors caused by ionizations. In both cases, the width of the Bragg peak was connected to the size of the illuminated region in the crystal, but was not affected by damage. We analyzed the shortcoming of the periodic boundary condition approach. The method overestimates radiation damage in the interaction region, so it gives an upper bound to the effect of radiation damage on the patterns. In the future, the simulation method developed here is to be applied to more complex scenarios.

This work has been supported by the excellence cluster “The Hamburg Centre for Ultrafast Imaging (CUI): Structure, Dynamics, and Control of Matter at the Atomic Scale” of the Deutsche Forschungsgemeinschaft.