Molecular dynamics simulations of inelastic x-ray scattering from shocked copper Open Access

The study of matter under ultra-high pressures and compressions has a history spanning more than a century.^1–9 An enduring and growing interest in this field is driven not only by a desire for a fundamental understanding of the behavior of matter under extreme conditions but also by the goal of linking such knowledge to a description of the interiors and dynamics of planets within our own solar system and beyond.^10–13 On a practical front, there also exists the prospect of finding and recovering new high-pressure phases that are metastable under ambient conditions, opening novel routes for materials synthesis.^14–18 

X-ray diffraction as a means of structure determination has been used since the earliest days of high-pressure research^19–23 and is now routinely used to map the phase diagram of matter subjected to high pressures within, for example, diamond anvil cells (DACs) used for static high-pressure experiments.^24–27 As well as such elastic scattering, the advent of high-brightness neutron and x-ray sources (such as synchrotrons) has led to the development of high-resolution inelastic scattering techniques,^28–33 whereby the exchange of momentum and energy of the scattered particle induced by the creation or destruction of a phonon has allowed phonon dispersion curves of matter under high pressure to be interrogated.^34–36 In the case of x rays, extremely high resolution is required for such measurements, given that phonons typically have energies up to a few tens of meV, while x-ray energies lie in the multi-keV domain.

As well as the above-mentioned method of using DACs to subject materials to high pressure, such extreme states can also be produced by rapid compression,^37–39 and, in particular, the past few years have seen a rise in interest in combining the technique of pressure applied by high-intensity nanosecond laser-ablation with a high-brightness transient x-ray source.^40–63 As local sound speeds generally increase with compression, in many cases, the rapid application of pressure causes the resulting compression wave to steepen into a shock, and the states in front of and behind the shock front are linked by the Hugoniot relations.^64,65 Shock compression is highly entropic, and as a result for most metals, the Hugoniot (the locus of points that can be reached by shock compression) crosses the melt curve at pressures of order 100–300 GPa.^66–68 However, with controlled laser pulses, shaped in time, or by compressing materials via multiple shocks (which can in principle be induced by appropriate target design⁶⁹), the application of high pressure can be achieved in a ramped manner,^57,70–74 without the compression wave steepening into a shock. This so-called “quasi-isentropic” (QI) compression has been shown to keep a material solid well into the TPa regime.^75–80 

Importantly, in conjunction with the development of these dynamic compression techniques, methods of producing x-ray diffraction patterns on the concomitant time-scales have also been developed, first by use of laser-plasma x-ray sources^{40,41,43,46,47,49,54,57,81} and subsequently with synchrotrons.^{59,61,82–85} More recently, significant advances have been made by using the output from so-called fourth-generation light sources—these are hard x-ray Free Electron Lasers (FELs)⁸⁶ with peak spectral brightnesses more than a billion times greater than those produced by any synchrotron. The use of such sources has led to a plethora of work that uses x-ray diffraction to study materials under shock and dynamic compression^{48,50–53,55,56,58,60,62,63,87} allowing structure and density to be deduced using x-ray pulses of a duration less than even the shortest phonon period in the system.

As noted above, in both shock and QI compression, density can be deduced via strong elastic x-ray scattering (i.e., Bragg or Laue diffraction). In combination with these structural measurements, pressure can be inferred from wave velocity profiles deduced from VISAR (Velocity Interferometer System for Any Reflector) measurements^88–90 and tabulated data about the equation of state (EOS) of the material. However, full information about the phase diagram remains elusive owing to the lack of any reliable temperature diagnostics. The final temperature induced during QI compression will be highly sensitive to the strength of the material—that is to say, its resistance to shear—as heat is generated by the plastic work performed during the relaxation of the induced shear-stresses. While it has been proposed that temperature could, in principle, be deduced via elastic scattering by monitoring the relative intensity of Bragg peaks and deducing the influence of the Debye–Waller effect on their intensity,^91,92 such an approach may be hampered by texturing effects and also relies on accurately knowing the effective Debye temperature of the system, which is itself a function of compression.

Given that the intensity of the Stokes and anti-Stokes peaks of inelastically scattered x rays (i.e., photons that during the scattering process have either created or annihilated a phonon) are determined by the principle of detailed balance for a system in thermodynamic equilibrium,^28,29,93,94 it has been suggested that inelastic scattering from transiently compressed matter could be used to determine temperature directly.^95,96 However, the transient nature of the compression process means that the crystal is only in the state of interest for a few nanoseconds at most. Owing in part to the very high resolutions required, inelastic x-ray scattering measurements are “photon-hungry,”^29,30 and it might seem that the prospect of retrieving information about the phonon modes on such timescales is unfeasible. Nevertheless, despite the considerable challenges that such measurements pose, progress is being made toward the development of inelastic scattering techniques on these short timescales. As noted above, FELs have an enormous peak spectral brightness. Indeed, recent work on both the Linac Coherent Light Source (LCLS)⁹⁵ and at the European XFEL (EuXFEL)^96,97 has demonstrated inelastic scattering from static, unshocked samples of almost perfect diamond single crystals. The temperature of samples at both ambient and elevated levels could be determined with an error of less than $8%$⁠. In these experiments, the photon numbers were low and just a few photons were collected for every FEL shot but this signal level can be increased significantly by seeding of the FEL.^98,99 Thus, when combined with a high-repetition-rate high-energy nanosecond optical laser for compressing samples, the prospect of obtaining meaningful inelastic scattering signals appears within reach. Indeed, this is one of the motivations for the building and placement of the DiPOLE laser¹⁰⁰ at the EuXFEL, DiPOLE being a 100 J diode-pumped laser capable of delivering pulses of a few nanoseconds duration at 10 Hz.

Thus, significant progress has been made experimentally in bringing about the prospect of inelastic scattering studies on the short timescales of relevance to dynamic compression. It is in the above context that we present here molecular dynamics (MD) simulations of x-ray inelastic scattering from the vibrational modes present within copper crystals shocked to pressures of order 70 GPa. While previous MD simulations have been post-processed to provide simulated high-pressure diffraction (i.e., elastic scattering) patterns,^{43,72,73,101–106} to our knowledge little work has been reported which attempts to investigate inelastic scattering in this regime. Our aim here is to provide a modest first step in this direction, assessing the feasibility of using x-ray inelastic scattering to gather information about the phonons present behind the shock front in a shocked single crystal. We note that the extremely high plastic strain rates present in such systems indicate that copious dislocations and other defects are homogeneously generated at the shock front, and these have certainly been seen in previous MD simulations, where, for example, dislocation densities of order $1013cm−2$ have been reported for shock pressures of order 100 GPa.^{44,101,107–111} Such a high density causes disruption to the perfect periodicity of the lattice on the length scale of tens of lattice spacings, raising the question as to the lifetime of coherent oscillations within the system and thus whether distinct phonon modes can be observed. Furthermore, such defects, while mobile during the time that they are relieving shear stress, will by interrupting the lattice give rise to a quasi-elastic scattering (QES) signal, even when scattering is occurring far from a Bragg peak. Experimentally, this quasi-elastic scattering must be differentiated from the inelastic scattering from the phonons, potentially placing limits on the energy resolution required.

While we will not attempt to address fully all of the above questions here, we will show for the conditions we have chosen that although copious defects are indeed generated when the crystal is shocked above the Hugoniot elastic limit (HEL), distinct phonon modes can still be discerned, though the phonon dispersion relation is, as expected, significantly modified by the compression process. Furthermore, the defects present (in this case mainly stacking faults) do indeed give rise to a quasi-elastic signal that competes with the inelastic features. As we shall show below, at least for these single-crystal simulations, the inelastic signal should still be discernible in the first Brillouin zone but care needs to be taken if measurements are made beyond this region of reciprocal space, as the stacking faults present will, at least along certain directions, completely overwhelm the inelastic signal. Finally, we note that given the level of energy resolutions achievable, in an experimental signal the particle velocity of the compressed portion of the sample will give rise to a potentially measurable Doppler shift of the scattered x rays that could be used as a further diagnostic of the system.

The paper is laid out in the following manner. First, in Sec. II, we describe the parameters of the simulations, both in terms of the MD simulations themselves and the Fourier transform methods employed. In Sec. III, we present our results, giving plots of the simulated scattered intensity as a function of wavevector and energy transfer and also lay out a pathway to obtaining the mean particle velocity behind the shock front via the Doppler effect. The implications of the results and proposals for further investigations are discussed in Sec. IV.

All simulations were performed using the open-source code Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)¹¹³ on a parallelized CPU cluster. The atomic response was modeled using the Cu EAM01 Mishin 2001 potential,¹¹⁴ which specifically was designed for high-pressure applications, accurately models dominant plasticity, e.g., stacking-fault energies, and hence is commonly used for shock simulations. The simulation box size was set to $50×50×455$ unit cells in the [100], [010], and $[001]$ directions, respectively. The total of $4.55×106$ copper atoms were arranged in a perfect fcc lattice that was aligned with the coordinate system. In the x- and y-directions, the boundary conditions were periodic, while in the z-direction the boundary condition was left non-periodic to allow for compression.

The simulations were run with a time step of 1 fs and a total duration of 30 ps. At first, the initialized fcc lattice is thermalized at 300 K for the first 2 ps under constant-NVE conditions, creating a realistic ambient crystal with temperature and pressure fluctuations on the order of $±1$ K and $±0.5$ GPa, respectively. The crystal is then shock-compressed by driving a rigid piston consisting of $50×50×15$ unit cells at the lower z-edge, at a fixed velocity, for 28 ps along $[001]$⁠. The final pressure reached ranged from 0 to 70 GPa depending on the velocity of the piston (i.e., the particle velocity behind the shock front). During the last approximately 8 ps of each of the eight simulations, the real-space positions of all of the atoms in a central, fully shocked region of $50×50×50$ unit cells extent and an atom count of approximately 500 000 are saved (“dumped”) every 32 fs; any further analysis will be only regarding these sets of atoms. The edge of the output box is approximately 50 unit cells distant from the piston, and its spatial extent is reduced for simulations with higher final pressure given the compression of the lattice, such that the number of dumped unit cells, and thus atoms, approximately stays the same.

The typical response of the crystal to the shock compression can be seen in Figs. 1 and 2, where the atoms in the crystal have been color-coded via a post-processing algorithm using the OVITO code,¹¹⁵ where the colors are assigned according to the common-neighbor-analysis (CNA)^116–118 and the Dislocation Extraction Algorithm (DXA),¹¹² with atoms in an fcc environment being colored green, and hcp atoms as red. The two plots shown correspond to pressures behind the shock front of 30 and 70 GPa, respectively. As has been found previously, for an initially perfect copper crystal, the Mishin potential predicts an onset of shock-induced plasticity¹¹⁹ (at the Hugoniot elastic limit or HEL) at around 30 GPa,¹²⁰ with entirely elastic behavior below this pressure. In Fig. 1, this onset of plasticity is evinced by the appearance of stacking faults, which appear as thin regions of hcp material in the ${111}$ planes. Furthermore, we note that for the Mishin potential employed here, a fcc-hcp phase boundary at a shock pressure of approximately 55 GPa has been found in a computational recrystallization study.¹²¹ However, while in shock simulations, the fcc-hcp fraction saturates at equal fcc-hcp phase fraction at the highest pressures explored in this work (60 and 70 GPa), and no full fcc to hcp phase transition is found.¹²² During a shock, the hcp material nucleated by plasticity can be thermodynamically similarly favorable as fcc and thus forms extended bands of hcp material several atomic layers thick that cannot be realized by plasticity alone, as seen in Fig. 2. Nonetheless, we shall refer to these hcp bands as “stacking faults” throughout since the fcc phase fraction always exceeds the hcp phase fraction, or is on the same order of, for the shock pressures up to 70 GPa.

Our approach to modeling the elastic and inelastic x-ray scattering from these shocked crystals is to take the space and time Fourier transform of their atomic coordinates. The intensity of x-rays scattered with a change of wavevector $k$ with an energy change of $ℏω$ is calculated via

with the atomic positions in real-space $ri$ and the simulation time steps $tj$⁠. Given the constraints of computational expense, the calculation of the inelastic scattering (i.e., that which is a function of $ω$⁠) is only undertaken along specific directions of interest in reciprocal space—i.e., in one dimension in reciprocal (⁠$k$⁠) space plus one dimension in frequency (⁠$ω$⁠) space. We can first identify such directions of interest using calculations of the elastic scattering alone, which can be modeled from the atomic coordinates dumped at a single moment in time but now calculated over a wide, three-dimensional volume of reciprocal space, as has previously been outlined by Warren¹²³ and is analogous to the calculations of Velterop et al.¹²⁴ and Kimminau et al.^102,125,126 Note that the atoms are modeled with neither extent nor form factor, i.e., as 3D delta-functions in real space. In reality, a form factor arising from the spatially extended electron distribution will affect the long-range behavior in reciprocal space, which in this work, however, is of lesser importance as for any given k-value, this intensity-reducing effect identically affects both the elastic and inelastic signals alike.

We note that strictly speaking, our approach cannot capture several aspects of the relevant physics, especially for the modeling of inelastic scattering. First, phonons are quanta of vibrations of the normal modes of the lattice, and their discrete nature is not captured in a classical simulation. Second, in merely taking the temporal Fourier transform, we completely ignore the creation and annihilation of phonons that occurs in the inelastic scattering process. However, we contend that this approach still provides useful information in that for modes with energies low compared with the thermal energy, the classical amplitude predicted by the MD should be a good indicator of the expectation value of the quantum amplitude. Indeed, the use of classical MD simulations to identify and study physical processes mediated by phonons is an established technique.^127,128 As phonon creation and annihilation are not taken into account, the intensities of the Stokes and anti-Stokes peaks will not be related by detailed balance, which should be the case in practice. If we wished to, we could mimic this by weighting the spectra as a function of $k$ and $ω$ by the factor²⁹ of

with the mode intensity $I$⁠, the momentum transfer $k$⁠, the mode energy $ℏω$⁠, the Boltzmann constant $kB$⁠, and the absolute temperature $T$ in the sample. As this is trivial and does not lead to any deeper understanding, this is not done in this work.

In undertaking these calculations, consideration must be given to the appropriate resolution in reciprocal space. Given that the largest spatial scale is that of the simulation box (which in reciprocal space will, as we see, give rise to diffraction features itself), we choose the spacing of the k-points in reciprocal space to be such that such size effects are resolved—i.e., in this work $Δk≤π/4Li$⁠, with $Li$ the simulation box extent in direction $x$⁠, $y$⁠, or $z$⁠. For the inelastic scattering calculations, the atomic positions were output every 32 fs for 8.192 ps and thus $Ntime=256$ data points were obtained at a sampling rate of 31.25 THz, yielding an energy resolution of 0.5 meV covering an energy range of $±64.62$ meV. To mitigate the influence of the finite time box, a periodic Hann window filter is applied. Note that the finite mean velocity of the atoms in the shocked crystal (i.e., the particle velocity) will introduce a phase shift into the Fourier transform. In the results given below, this phase shift is removed but we comment on the observation of such a shift (which corresponds to the Doppler effect) in Sec. III C.

An example of such a calculation for the intensity of the elastic (⁠$ω=0$⁠) scattering is shown in Fig. 3. Here, we plot the square modulus of the Fourier transform for a crystal shocked to 30 GPa (i.e., just above the HEL) as an intensity-threshold plot around the origin for a region in reciprocal space bounded by points with Miller indices $(±3,±3,±3.45)$⁠, where the indices are referenced to the ambient unshocked material. At this shock pressure, the crystal has been compressed to 87% of its original volume.

This elastic scattering calculation exhibits several features of interest. First, we note that from each of the intense points in reciprocal space, we observe streaks along the $⟨100⟩$ directions (i.e., along the coordinate axes). We highlight this feature at the origin of reciprocal space in the figure, by color-coding it red. This feature is the artifact of the finite size of the simulation box in reciprocal space, i.e., it is the well-known result that an orthorhombic shape in real space of dimensions $Lx$⁠, $Ly$⁠, and $Lz$ will have a “shape function” in reciprocal space that, in the vicinity of the reciprocal lattice vector $G$⁠, gives rise to a diffracted intensity

which has prominent “arms” along the $⟨100⟩$ directions.¹²³ Indeed, the intensity in these directions due to the shape function of the box of atoms is so great that it overwhelms any thermal diffuse scattering. For the shocked crystals, the same intensity pattern will dress every reciprocal lattice point, which will have shifted owing to the crystal compression.

However, as well as this artifact due to finite size effects, lines of high intensity in reciprocal space can be seen linking certain reciprocal lattice points along different directions. Specifically, there is considerable intensity between some of the reciprocal lattice points along $[1¯11]$ and $[111¯]$⁠. This quasi-elastic scattering is associated with the planar stacking faults seen in Figs. 1 and 2.^129,130 The fact that particular points in reciprocal space will have a large streak in intensity emanating from them along directions associated with the normals to stacking faults is well known. As shown by Warren¹²³ and has been seen in previous atomistic and MD simulations,^108,129,130 the phase shift, $Δϕ$⁠, introduced into an x ray as it traverses a stacking fault in the $(111)$ plane is given by

where $n111$ is the number of faulted planes. Analogous relations hold for the other stacking faults in the ${111}$ family of close-packed planes. We thus expect linear streaks in reciprocal space caused by the disruption to the real lattice by the stacking faults, along directions that are the normals to the stacking fault planes, for those reciprocal lattice points for which $Δϕ≠2mπ$⁠, where $m$ is an integer.

A study of Fig. 3 shows that the main features are consistent with the above picture: strong scattering due to the stacking faults is not seen at the origin (where only the shape function is visible, at least at the threshold intensity in the plot) or, for example, around the reciprocal lattice points corresponding to ${3,3,3}$ in the unshocked system. As a result, such effects are not prominent in the first Brillouin zone. In contrast, strong elastic scattering due to the stacking faults can be seen close to, for example, reciprocal lattice points corresponding to ${2,0,0}$ in the unshocked crystal. This reduction in scattering due to these planar defects around lattice points for which $Δϕ=2mπ$ with integer $m$ is referred to as the “invisibility condition.”¹²⁹ As we shall find below, the density of stacking faults found in these shocked samples is actually so high, and the lattice so distorted, that the effects of the faulting can still be discerned, although at a level considerably reduced compared to the case where the invisibility condition is not satisfied.

To elucidate the above observations of this quasi-elastic scattering signal and the effect that the stacking faults have upon it, we will consider the inelastic scattering along two high-symmetry paths in reciprocal space: the first path, which is labeled $P1$ in Fig. 4, connects reciprocal lattice points $(2,0,0)$ and $(1,1,1)$ along the $[1¯11]$ direction; the second path, $P2$⁠, connects the $(0,0,0)$ and $(1¯,1,1)$ peaks. For an fcc crystal with perfect translational symmetry, in which all Brillouin zones are exactly equivalent, the reciprocal-space intensity along paths $P1$ and $P2$ would, of course, be identical. However, in the presence of the copious shock-induced stacking faults, $P1$ may coincide with an intense streak of quasi-elastic scattering. This signal—when it exists—arises specifically from stacking faults of the type $(1¯11)$⁠, whose invisibility criterion is failed by both $(2,0,0)$ and $(1,1,1)$ (and indeed by every reciprocal lattice point along the line passing through them). Path $P2$⁠, by contrast, is permanently shielded from any such quasi-elastic signal because $(0,0,0)$ identically satisfies the invisibility criterion of every stacking-fault variant. The intensity of the quasi-elastic scattering one observes between the Bragg peaks thus depends strongly on which Brillouin zone one probes.

With this in mind, we now go on to examine the full spatiotemporal Fourier transforms of the crystal along paths $P1$ and $P2$ in reciprocal space. As shall be shown, the linear quasi-elastic features also have a signature that extends into the energy domain, which must be taken into account when considering inelastic scattering.

To illustrate the significant impact of the stacking-fault features found in the ${111}$-directions, we show in Fig. 5 the phonon modes along path $P1$ joining $(2,0,0)$ and $(1,1,1)$ for shock pressures between 0 and 70 GPa using the method described in Sec. II B. The color scale indicates the intensity of a mode of energy $ℏω$ with momentum $ℏk$⁠. At shock pressures below 30 GPa, distinct $(2,0,0)$ and $(1,1,1)$ Bragg peaks can clearly be observed along the zero-energy axis. Note that the effect of the star-like (triple-sinc-squared) shape function is only felt very close to the peaks because in the ${111}$ directions, it decays extremely quickly. As we ramp up the pressure, the Bragg positions move in k-space, which is attributable to the shock-induced lattice compression, and, from 30 GPa onward, to small plasticity-induced rotations^46,55,119 together with stacking-fault-induced phase shifts.^108,123 These effects give rise to increased distances in reciprocal space between specific reciprocal lattice points along $P1$ and $P2$ in addition to apparent rotations of these paths of up to $3°$—details of how paths $P1$ and $P2$ migrate due to compression and plasticity effects are shown for reference in Table I.

Above 30 GPa, there also exists between the Bragg peaks themselves a strong quasi-elastic signal along the zero-energy axis that arises from the presence of defects [namely, the stacking faults of the kind $(1¯11)$⁠, which do not undergo significant changes on the lattice vibration time scale] whose pressure-dependent intensity is shown in Fig. 6. Given that the distribution and abundance of ${111}$ stacking faults vary between simulations, we find no trivial pressure dependence of the plotted QES signal once we enter the plastic regime. While the distinct jump between 20 and 30 GPa is clearly related to the sudden increase in defect density, we can only resolve the binary relation that once stacking faults in $(1¯11)$ appear, the central QES signal increases by approximately five orders of magnitude relative to the ambient case; beyond this, the complex interaction of various defects and rotations means the signal strength along $P1$ in the plastic regime fluctuates by two orders of magnitude between different simulations.

As well as being connected by these quasi-elastic signals, the Bragg peaks are also linked up by arcs that span nonzero-energy regions of the spatiotemporal Fourier transform. These show the intensities of the classical Stokes- and anti-Stokes modes corresponding to the Cu $[qqq]$ longitudinal-acoustic (LA) phonon (upper arc) and the $[qqq]$ transversal-acoustic (TA) phonon (lower arc). Compression of the crystal along its $[001]$ direction causes its elastic moduli to increase via pressure hardening, which, in turn, causes the LA phonon energy to scale with shock pressure. This phenomenon is illustrated in Fig. 7, in which the maximum energy of the classical modes at $(32,12,12)$ is plotted as a function of shock pressure. We note in passing that the TA modes can be described by components oscillating in the $(001)$ plane and perpendicular to it in $[001]$⁠. Hence, upon elastic compression of the crystal along $[001]$⁠, the associated loss of isotropy causes the TA modes to split into one approximately ambient component related to the $(001)$ plane (lower TA) and one pressure-hardened component along $[001]$ (upper TA). At 20 GPa, one can clearly distinguish two distinct TA modes because of the large anisotropy of the unit cell at approximately $5%$ elastic shear strain. Above the elastic limit, the plastic flow mediated by the stacking faults permits some fraction of the deviatoric elastic shear strains to relax and the residual shear strain is only about $1%$⁠. This lowers the elastic anisotropy of the crystal to such a degree that the two TA modes can no longer easily be distinguished but instead form a continuum whose central frequency grows with pressure, much like the LA mode. The tendency of the phonon frequencies to increase with pressure should allow them to be more easily distinguished from the QES, a point to which we shall return in Sec. III B.

In these shock simulations, the overall noise in the spatiotemporal Fourier transform increases at higher pressure due to a confluence of higher sample temperature, higher defect density, and also the emergence of the hcp phase (the volume fraction of which is shown in Table I). When taken into consideration along with the aforementioned stacking-fault-induced QES, these plasticity effects can cause the phonon modes to become far harder to resolve at greater pressures. Such signal reduction plays an important role in real experiments, which is why next, the limiting effects of realistic instrument functions (IFs) will be explored.

Dynamical compression experiments involving laser-driven shocked samples probed using XFELs are of far-reaching interest to the high energy density (HED) physics community. The inelastic x-ray scattering (IXS) technique represents another means of exploiting the power of the XFEL to diagnose material behavior under extreme conditions. The method does, however, place stringent conditions upon the energy resolution required of the instrument used to collect the inelastic signal. The resolution in the simulations presented in this paper is 0.5 meV, which is yet to be matched by dynamic experiments: at the IXS commissioning experiment of the high energy density (HED) endstation of the EuXFEL, the instrument function was on the order of 50 meV⁹⁷—this has been and will continue to be further improved upon in the future at both EuXFEL and LCLS. While the influence of the instrument function on the scattering signal is, in principle, deconvolvable, shot noise poses a limit due to finite integration time at experiments. In this section, we will use the spatiotemporal Fourier transforms calculated in Sec. III A to estimate the kind of energy resolution required to distinguish the phonon modes in shock-compressed copper from the background QES owed to the copious defects.

In IXS geometries, usually photons that correspond to a small range of momentum transfer are recorded, which would correspond to a vertical line-out in Fig. 5. For this reason, we show in Fig. 8 raw line-outs of the Fourier transforms along the energy axis at a scattering vector of $(32,12,12)$ (i.e., at the midpoint of $P1$⁠) for simulations at 0 and 30 GPa. Shown also are the same signals convolved with a normalized Gaussian of FWHM 10 (left) and 50 meV (right), whose purpose is to emulate the instrument function. The Stokes- and anti-Stokes peaks are, as expected, symmetrically spaced around the zero-energy modes. Given the classical nature of the simulations, aside from numerical fluctuations, their relative intensities are matched and are not subject to the detailed balance described by Eq. (2)—as mentioned previously, correcting for this detail is possible but unnecessary for the purposes of the present study. The key question we wish to address is whether these peaks can still be distinguished in the presence of the background QES.

At ambient pressure, peaks from the LA and TA phonons can be identified at approximately $±34$ and $±14$ meV respectively, while at 30 GPa (in the plastic regime) the QES arising from the stacking fault streaks overwhelmingly dominates the spectra, and only a faint TA signature can be made out at approximately $±15$ meV. For an instrument function with FWHM 50 meV, phonons can be resolved in neither the ambient nor shocked case, and the signal is simply proportional to the total scattered photons. For an instrument function of just 10 meV FWHM, the TA modes are resolvable, but only under ambient conditions; the ambient LA modes, meanwhile, greatly suffer from noise due to their being in the second Brillouin zone (where the effects of stacking faults are large) and so remain unresolvable. Once the plastic regime is entered at 30 GPa, not even the TA phonons can be resolved thanks to the extremely intense QES, even with such high energy resolution. These spectra indicate that the presence of stacking faults could render measurement of phonon modes in the second Brillouin zone unfeasible, at least along these specific directions in the second Brillouin zone in shock-loaded single-crystal copper.

However, as was illustrated in Fig. 4, the stacking faults should not adversely affect the first Brillouin zone, which is why it is instructive to consider the same spectra obtained midway along $P2$⁠, i.e., still along the $[1¯11]$ direction, yet starting at the origin rather than $(2,0,0)$⁠. The corresponding spatiotemporal Fourier transforms are illustrated in Fig. 9. Just as in Fig. 5, Bragg peaks are visible along the zero-energy axis [in this case, $(0,0,0)$ and $(1¯,1,1)$] that are again connected by LA and TA modes. Once more, a splitting of the TA modes into two branches is visible in the elastic regime, becoming less distinct from 30 GPa onward. Unlike those along path $P1$⁠, the TA modes here in the first Brillouin zone are noticeably fainter than are the LA modes. The most striking difference between Figs. 9 and 5, however, is in the contrast between the classical modes and the QES: along path $P2$⁠, the QES does not completely dominate the signal as it did along path $P1$⁠. This, as explained in Sec. III, is thanks to the origin satisfying the invisibility criterion for all four stacking-fault variants.

We show in Fig. 10 the corresponding raw line-outs and synthetic experimental spectra at scattering vector $(1¯2,12,12)$⁠, where the spectra are shown again for 0 and 30 GPa. The LA (TA) modes in this case are found at $±34$ (⁠$±14$⁠) meV for 0 GPa and at $±43$ (⁠$±16$⁠) meV for 30 GPa. A relatively weak yet still significant background QES is present, which may be due to the small but finite strains associated with the network of copious stacking faults and the general increased disorder in the plastic region behind the shock front. For the ambient crystal, the LA modes can be resolved with an instrument function as wide as 50 meV FWHM due to their higher energies compared with the dominant TA modes at $(32,12,12)$⁠. However, the QES still dominates the synthetic spectra once we enter the plastic regime, and at 30 GPa, no phonons can be resolved once again. For an instrument function of only 10 meV FWHM, however, the LA modes in both ambient and shocked copper can be measured, and at 30 GPa, even the effect of the TA modes on the spectrum can be resolved. Ignoring the detrimental effect of the increased background QES that inevitably accompanies it, we are assisted to a great extent here by pressure hardening: the inelastic peaks at higher pressure can be more easily identified due to the increased separation of the phonons in energy space—a crucial benefit conferred by the compression itself. In any case, these spectra indicate that, even within the first Brillouin zone, energy resolutions better than 50 meV will be required to perform IXS measurements of shock-loaded single-crystal copper, at least along the directions studied.

As noted above, the current resolutions that have been demonstrated at fourth-generation light sources are in the region of 45 meV at a photon energy of 7.5 keV—i.e., a resolution better than one part in $105$⁠. In future developments, it is envisaged that the use of even higher photon energies—with even greater resolution—will be possible. Indeed, resolutions of order just a few meV (e.g., 3 meV at 17.79 keV) have been demonstrated in synchrotron settings.³² It is important to note that the shock or ramp compression process will impart a bulk velocity shift to the compressed portion of the crystal—in the parlance of the shock community this is the “particle velocity” (as opposed to the shock velocity at which the compression front moves). With the above spectral resolution, this bulk velocity will exhibit itself as a Doppler shift of the scattered signal whenever the collection geometry differs from pure forward scattering. In this final section, we discuss the prospects of exploiting this Doppler shift to experimentally measure the particle velocity behind the shock front.

Typically, for shocks in the 100 GPa range, the particle velocity $Up$ is several $1000ms−1$⁠, resulting in $β=Up/c$ of the order of $10−5$⁠, with $c$ being the speed of light in vacuum. This factor is sufficiently large that the relativistic Doppler effect will need to be taken into account when analyzing the scattered signal: indeed, it will contain a direct measure of the particle velocity within the material. In our simulations, the bulk sample velocity $Up$ behind the shock front equals that of the driving piston in $[001]$⁠, which we were able to “correct for” by transforming into the center-of-mass frame of the shocked material. Without this correction, any spatiotemporal line-out with a component in $[001]$ would experience a Doppler shift.

In general, for a given x-ray photon with energy $ℏω0$⁠, the scattered frequency after absorption and re-emission is

with the angle of absorption $Ωa$⁠, and the angle of emission $Ωe$ in the rest frame of the source and of the detector.¹³¹ For a sample shocked to $Up$ in direction $n^s$ and probed by x rays with incident and scattered wavevectors $k0$ and $k$⁠, the angle of absorption follows as $Ωa=arccosn^s⋅k0k0$⁠, and the angle of emission as $Ωe=arccosn^s⋅kk$⁠. Given that for x-ray Thomson scattering the magnitude of the scattering vector $q$ follows as $q=2k0sinΘ2$⁠, this Doppler shift will hence be strongly dependent on the scattering angle.¹²³ Without loss of generality, setting

but subsequently confining scattering to the $yz$-plane, $k$ reduces to

This is illustrated in Fig. 11 and in Table II, a selection of such Doppler shift cases is shown. Note that this shift tends toward unity for direct forward scattering, for which the effects of absorption and re-emission cancel out. Typically, $α$ is small, and thus the energy shift is most significant for backscatter. Considering the example of copper at 100 GPa with $Up=1714ms−1$ at $α=10°$⁠, $Θ=30°$⁠, and $ℏω0=7.5$ keV, as is typical for the current EuXFEL HED configuration, this would cause a shift of just 2 meV. For backscatter, however, at $Θ=170°$⁠, this increases to 82 meV, significantly shifting the signal and providing the opportunity to directly extract $Up$ even with the current modest spectral resolutions of 45 meV. If higher resolutions in the sub-5-meV range could be obtained at FELs at the greater photon energies cited above, then particle velocities should be directly measurable with accuracies well below $100ms−1$⁠. While this Doppler-shift approach to measuring particle velocity cannot compete with a standard VISAR diagnostic,^88–90 it would give a velocity measurement within the material, without need for observation of the motion of a boundary, and without making the assumption that the velocity at a free surface is twice the particle velocity. Furthermore, as it does not rely on optical reflection from a surface, such an approach would not be subject to “blanking” due to poor reflectivity. These considerations illustrate that the Doppler effect in IXS could yet play an important role in the diagnosis of materials rapidly compressed by laser-ablation.

In summary, this work has made two modest steps toward understanding the physics of shock experiments that use high-resolution IXS as a diagnostic. First, we have shown how and why quasi-elastic scattering can be a challenge when measuring phonons via inelastic scattering, particularly when intending to reliably determine the absolute temperature of shock-compressed material by trying to resolve the Stokes and anti-Stokes peaks. Second, we demonstrated that through the shock-induced Doppler shift of the inelastic scattering signal, it should be possible to directly extract the yet-inaccessible net particle velocity behind the shock front.

The simulations presented here provide insight into a rich selection of physics. Upon shock, the phonon energies rise due to pressure hardening. The resulting enhanced separation of the Stokes- and anti-Stokes peaks proves experimentally to be advantageous, as it relaxes the energy resolution required to distinguish them. However, we have also shown that collecting photons between the Bragg peaks (in order both to maximize this phonon energy and to avoid elastic scattering signal) can still be a challenging task when looking at defective matter traversed by stacking faults, which give rise to QES. In particular, the 3D spatial Fourier transforms show that certain regions in k-space are exceptionally affected, via linear features connecting certain Bragg peaks in higher Brillouin zones. Indeed, due to the overall increased disorder of the shock-compressed system, even within the first Brillouin zone, there exists a background QES that could impair phonon measurements if the diagnostic’s energy resolution were insufficiently great, at least along the particular directions we have studied. The simulations herein have provided an indication of the kind of resolution required to distinguish the inelastic peaks from the quasi-elastic signal threatening to swamp them (i.e., a few tens of meV).

To generate the spatiotemporal Fourier transforms shown in Figs. 5 and 9, we moved the reciprocal-space paths along which we were probing “in sync” with the compression, in order that our paths always coincided with the Bragg peaks. For real experiments, however, it is more common to commit to one small, predetermined k-range. This could mean that the streaks drift through the detector range, potentially saturating any phonon signal in temporally integrating setups. That being said, while the $⟨111⟩$-streaks caused by stacking faults can impair inelastic measurements, deliberately probing higher Brillouin zones so as to maximize such linear features could be utilized to make statements about the plasticity-induced microstructure of the shocked material.

This work on copper dealt with a single fcc crystal shocked along the $[001]$ direction using the extensively tested Mishin potential. While the ${111}$ family of stacking faults is most strongly generated when shocking along $[001]$⁠, the behavior upon shocking in other directions of high symmetry with different defect structures are prime candidates of future simulation work, such as $[111]$ and $[110]$ that produce more isotropic plasticity.¹³² Similarly relevant will be to investigate more elaborate crystal structures and different materials; of particular importance are Debye–Scherrer-type experiments, in which the scattering signals (inelastic and quasi-elastic) result from an average over signals from every grain in a full polycrystalline aggregate. Simulating the signals obtained from polycrystalline targets represents a natural avenue for further work but would be computationally far more intensive.

In future experiments, exploiting the Doppler shift of the inelastic signal could be an interesting way of directly monitoring the particle velocity, complementing VISAR measurements; the relative novelty of XFELs means using IXS in such a way (in environments characterized by very large particle velocities) has never been done before. We expect mean velocities of the order of several $1000ms−1$ that cause energy shifts of tens of meV, which would make direct measurements of particle velocities within the interior of shocked systems possible. The simulations performed here all entailed compression via a largely over-driven plastic shock wave, in which there exists just one near-discontinuity in particle velocity. However, other loading scenarios exist in which the wave structure is more complex, giving rise to distinct signatures in the IXS. For example, upon ramp compression, a gradient in the particle velocities would broaden the energy-resolved diffraction signal and hence should provide valuable in situ information on the ramp itself, such as a three-wave structure that can arise in certain crystal system and shock configurations, which is composed of an elastic wave, a plastic wave, and a phase-transition wave at the same time, each of which causes a jump in the particle velocity.

O.K., P.G.H., and J.S.W. gratefully acknowledge support from AWE via the Oxford Centre for High Energy Density Science (OxCHEDS), and O.K. acknowledges further support from the Engineering and Physical Sciences Research Council (EPSRC). J.S.W. is also grateful for support from the EPSRC under Grant No. EP/S025065/1. R.E.R.’s work was performed under the auspices of the U.S. Department of Energy (DOE) at the Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 and the computing resources were provided in part by the LLNL Institutional Computing Grand Challenge program.

The authors declare no competing financial interests. UK Ministry of Defence^© Crown Owned Copyright 2021/AWE.

The data that support the findings of this study are available from the corresponding author upon reasonable request.