Atomic-scale diffractive imaging of sub-cycle electron dynamics in condensed matter

For interaction of light with condensed-matter systems, we show with simulations that ultrafast electron and X-ray diffraction can provide a time-dependent record of charge-density maps with sub-cycle and atomic-scale resolutions. Using graphene as an example material, we predict that diffraction can reveal localised atomic-scale origins of optical and electronic phenomena. In particular, we point out nontrivial relations between microscopic electric current and density in undoped graphene.

For interaction of light with condensed-matter systems, we show with simulations that ultrafast electron and X-ray diffraction can provide a time-dependent record of charge-density maps with sub-cycle and atomic-scale resolutions. Using graphene as an example material, we predict that diffraction can reveal localised atomic-scale origins of optical and electronic phenomena. In particular, we point out nontrivial relations between microscopic electric current and density in undoped graphene.
Any light-matter interaction is, at its fundamental level, defined by field-driven motion of charges in their atomic-scale environment. Depending on a material's structure, potential-energy landscape, correlation mechanisms and applied field strength, this motion may be highly complex, and it causes a rich variety of macroscopic optical phenomena, which are at the heart of current laser technologies 1 , quantum optics 2,3 , spectroscopy 4 , nonlinear microscopy 5,6 , frequency comb metrology 7,8 , optical information processing 9,10 , and light-driven electronics [11][12][13] .
In this paper, we take a real-space/sub-cycle perspective seeking the most direct way to comprehend the nature of optical and electronic phenomena at lightwave frequencies. We consider the following questions: Within a laser cycle, what fraction of the material's charge density is driven into which direction resulting in which macroscopic response? Are there certain electron densities that are particularly efficient for specific effects? And what experimental approach can possibly reveal such three-dimensional motion of charges in real space? Sub-cell and sub-cycle resolution in space and time, if feasible, would arguably provide the most straightforward access to the basic microscopic effects underlying complex light-matter interaction.

Concept for sub-cycle imaging in real space
The example system considered in this work is graphene, which is well known for its unusual electronic and optical properties 14 and high potential for future high-speed electronics. Here, graphene is chosen due to its relatively simple structure and the large number of valence electrons per nucleus. The imaging concept is ultrafast pump-probe diffraction with electron or X-ray pulses of sub-laser-cycle duration, short enough to freeze out the relevant electronic and resultant atomic motion driven by light. Such pulses are achievable with free-electron lasers and, potentially, also with electron pulses, as we discuss below. A series of diffraction snapshots recorded under such conditions for different delays with respect to the driving field can provide real-time insight into rearrangement of electronic charge in condensed matter with spatial and temporal resolution. This idea was introduced earlier 15 and modelled in detail for electronic superposition states in atoms and diatomic molecules [16][17][18][19][20] , as well as for polyatomic molecules on much slower time scales 21 . Seminal experiments with femtosecond X-rays 22-24 revealed sensitivity of ultrafast diffraction to light-field-driven electron dynamics, albeit in a cycle-averaged experiment without sufficient temporal resolution. Our main goal here is to demonstrate the feasibility and importance of fully resolving electronic motion driven by the oscillating light fields in space and time. We also show how graphene's complex relations between atomic-and macroscopic-scale phenomena may be revealed by ultrafast pump-probe diffraction. Figure 1(a) depicts a schematic of the proposed experiment. A material (here graphene) interacts with a 10-fs laser pulse carried at a frequency of 140 THz, i.e., a wavelength of 2.1 μm. For such a pulse, the direct excitation of phonons is negligible, and graphene's response is mostly electronic. The 7-fs optical period of a 2.1-μm pulse is short enough to neglect, in a first approximation, lattice motion; at the same time, it allows for realistic sub-cycle diffraction measurements with few-fs electron and X-ray pulses. As a peak electric field, we use 0.5 V/Å, which is high enough for the onset of nonlinear and strong-field effects 25 while still below graphene's damage threshold 26 . Current laser technology readily provides such pulses 1,27 . The optical electric field in our calculations is polarised along graphene's Γ -K direction, i.e. in real space from the hexagon centre towards the middle of a C-C bond. Snapshots of the moving charges are taken by pump-probe diffraction with ultrashort electron or X-ray pulses of sub-cycle duration for later slow-motion replay.
Real-space electron dynamics and its complexity. First we report a simulation of graphene's electron dynamics in real space and discuss its relation to macroscopically observable optical phenomena. The method employed is a numerical solution of the time-dependent Schrödinger equation for a phenomenological lattice potential; see Methods. Figure 1(b) shows the static density of graphene's valence electrons without the optical field. From a chemical perspective, graphene is a sp 2 -hybridised structure; each carbon atom, which has electron configuration 1s 2 2s 2 2p 2 , contributes three of its four valence electrons to three sigma bonds and one to a delocalised electron density above and below the plane of the nuclei.
Hence, there are two electrons per sigma bond and two delocalised electrons per unit cell. In the computed ground-state charge-density map of Fig. 1(b), these are evident: the sigma bonds form density maxima between carbon atoms, and the p z electrons contribute to an additional delocalised density at larger z . For sufficiently large z , for example, at > . z 0 6 Å, there are local density maxima above/below each C-C bond, not above/below the nuclei. This is a consequence of the overlapping p z orbitals in conjunction with remaining parts of the bond-forming sp 2 densities. In other words, the sigma bonds and the delocalised electron density overlap significantly in real space. We note that there is also some remaining density of the valence electrons concentrated around the nuclei.
Electric current in undoped graphene can only exist if some electrons are excited into conduction bands. Let us consider a gedankenexperiment where we transfer valence electrons within a 2-eV energy range from the Fermi level and below into the lowest conduction band conserving the quasi-momentum. Note that this includes the electrons occupying the vicinities of the K points in reciprocal space. Figure 1(c) depicts the change of the electron density upon such excitation. There is no change in the graphene plane, but a maximum change at a distance of ≈ . z 0 6 Å above and below it. Remarkably, the electron density is transferred from two localised stripes (dark blue) into a central area (yellow) with a peak (red) at the hexagon centre. This real-space representation of the static band-population change provides the first sign of electronic motion's complexity in real space.
Next we discuss electron dynamics driven by the laser field. Generally, a material's response can be separated into localised charge displacements and delocalised charge motion. At low, perturbative fields, the response of dielectrics for frequencies below the bandgap is due to localised charge displacement causing their polarisation. In contrast, electric current in conducting media (metals, doped semiconductors, conjugated polymers, etc.) is caused by delocalised charges. In general, especially in complex materials, both localised and delocalised charges shape the optical response.
In undoped graphene, microscopic (atomic-scale) electron and current densities are, to a significant degree, complementary to each other in real space. This is evident from Fig. 1(d) where we plot magnitude of current density averaged over time. The result looks dramatically different from both the ground-state electron density of Fig. 1(b) and also the excited-state density change of Fig. 1(c). Although the electron density transferred to the conduction band is concentrated in the unit cell centre [ Fig. 1(c)], the current is concentrated at its periphery and ≈ . 0 6 Å above and below the atomic plane where the delocalised and overlapping p z orbitals are situated. Hence, the electric current predominantly flows through regions where these orbitals form continuous chains along the laser polarisation. Note that little or no current flows at the hexagon centre, where carbon orbitals have a negligible overlap. Also, there is no appreciable contribution to graphene's local or non-local currents originating from any C-C bond electrons in the atomic plane. These combined observations further suggest value of a real-space visualisation approach. Figure 1(e) identifies the spatial regions associated with optical absorption by depicting the residual change of electron density at t = 20 fs, after the excitation pulse has passed. The deposited energy per pulse is, in our model, 1.9 eV per unit cell. This causes no density change in the plane of the nuclei, but above and below it, charge is transferred towards the carbon hexagon centre where the initial electron density is low. This effect of the intense laser pulse is qualitatively different from that shown in Fig. 1(c) where electrons from the top of the valence band are artificially transferred to the conduction band. This difference is yet another indication of significance of real-space electron distributions.
To further elucidate the importance of charge movement in real space, we introduce some regions of interest, labelled A-G in Fig. 2(a). These regions are chosen to encompass distinguishable areas of electron density. The volumes associated with the C-C sigma bonds form two Bravais sub-lattices, one of electrons per unit cell, close to the overall number of 8 valence electrons of graphene. This indicates that all relevant charges are covered by our assignment.
In each region, the time-dependent effects of the laser field on the charge density are decomposed into two components: local displacement and redistribution between different volumes. We quantified the displacement by evaluating the first momentum (centre of mass) in each volume. For the C-C bond perpendicular to the laser polarisation (region A), Fig. 2(b) compares the ground-state electron density (solid) to that at the instance of maximum field strength (dashed). At this moment, the main effect of the field is displacement of the charge within the bond, though some decrease of the charge is also visible.
Temporal dynamics of displacements and charges are shown in Fig. 2(c-e) for the two types of C-C bonds (volumes A and B) and for the centre of the graphene hexagon (volume D). The field-induced local displacement of the electrons in volumes A and B is about one picometre; see Fig. 2(c). At all times, the shift is approximately proportional to the laser field. The C-C bond represented by volume B also experiences a shift that is proportional to the electric field [ Fig. 2(c), green]. The direction of the shift is shown by black arrows in Fig. 2(a) and it is not parallel to the laser field, as one could naively expect. Instead, the centre of mass of volume B moves predominantly along the bond. This is due to the polarisability of a bond being larger along its direction than normal to it. The valence charges around the nuclei (volume C) shift by only ~0.07 pm in the laser field and give/receive a maximum of only ~0.002 electron. These charges are, therefore, negligible for graphene's optical and electronic dynamics. . This reflects the symmetry of the system with respect to reflection → − x x in conjunction with pseudo-spin conservation implying charge conservation within each sub-lattice. As each of the two type-A bonds receives its peak charge increase of ~0.03 electrons, each of the four type-B bonds loses ~0.015 electrons. Our results show that this charge transfer mainly takes place at ≈ . z 0 6 Å, i.e., significantly out of the plane of graphene.
Unlike the charge displacement, the amount of charge in a certain volume is modulated in time with twice the laser frequency, and it is symmetric with respect to the sign of the laser field. Also, the phase of this field-induced charge migration between different types of bonds indicates that the vector potential, rather than the electric field, determines these dynamics: There is little charge transferred at the field extrema, while a maximum amount of charge is transferred between volumes A and B at the zero electric field. The fact that the magnitude of the electric current is also maximal at the zero field suggests that the atomic-scale charge redistribution shown in Fig. 2(d) is a prerequisite for inducing a significant electric current in undoped graphene. The effective charge transfer between regions A and ( , ′) B B is also largely symmetric with respect to → − y y and, consequently, negligibly contributes to net polarisation and current. This charge transfer effect is therefore hidden from far-field spectroscopic characterisation. Figure 2(e) shows time-dependent charge Q D in the centre of the unit cell, volume D. This charge has local maxima at the crests of the laser field and minima at its zero crossings. On a few-femtosecond time scale, there is gradual accumulation of charge in volume D that persists after the pulse end, which is also evident in Fig. 1(e). Physically, this implies that the charge accumulation in the unit cell centre, D, has both an adiabatic component (the local maxima following the absolute strength of the laser field) and a non-adiabatic component caused by resonant absorption. This absorption in graphene is broadband 28 , causing rapid dephasing (Landau damping), leading to the stationary residual electronic excitation.

Experimental feasibility of sub-cycle electron diffraction.
Here we propose an experiment for direct visualisation of purely electronic dynamics in graphene, and discuss its feasibility. The following three considerations are critical. First, the spatial resolution must be atomic. Second, temporal resolution must be sub-cycle. Third, there should be sufficient signal-to-noise ratio for the expected charge density changes, which manifest themselves in diffraction as time-dependent changes of Bragg peak intensities. In our case, sub-cycle means a quarter cycle of the laser, i.e. ~1.75 fs for the 2.1-μm driving laser considered here, in order to temporally resolve the largely nonlinear electron motion driven by a strong field, including, e.g., dynamics responsible for third harmonic generation.
Two currently emerging technologies can provide these prerequisites. Free-electron lasers (FELs) can deliver pulses with attosecond substructures 29,30 , and significant efforts are being made to scale up the photon energy toward wavelengths that would be sufficiently short for atomic-scale imaging 31,32 . The low scattering cross section of X-rays is compensated by the superior average photon flux of FELs; note that a particular peak flux is not required since the experiment is reversible. Some ionisation of the graphene sample eventually caused by a pump or a probe pulse will be rapidly replenished 33 . The second approach is ultrafast electron diffraction 34,35 . The de Broglie wavelength of electrons is as small as several picometres for kinetic energy of tens of kiloelectronvolts. The scattering efficiency from graphene is high, ~10 −3 from the central beam into a typical Bragg spot; see Methods. Ultrashort pulses can be produced by microwave compression 36,37 and, in the case of only one or a few electrons per pulse, simulations predict attosecond pulse durations 38,39 . Indeed, there is indirect experimental evidence 40 of an attosecond electron pulse train 15 .
Here we consider electron diffraction at 30 keV, an energy low enough for efficient scattering and high enough to provide a 7-pm de Broglie wavelength. Diffraction is modelled neglecting inelastic scattering and the electron transit time through the sample, which is a few attoseconds; see Methods. At a realistic 300 kHz pump-probe repetition rate 41 , about 10 9 single-electron pulses are incident on the sample within one-hour integration time. With a scattering probability of ~10 −3 , this allows realistically low shot noise in the Bragg peak intensities.
The inset in Fig. 3 shows the expected Bragg pattern, where we have indicated Miller indices for the two primary reciprocal-space vectors. For the first two diffraction orders, Fig. 3(b,c) show the scattering probabilities as functions of the delay between the laser field [ Fig. 3(a)] and a 1-fs electron pulse. Notably, the atomic-scale electron dynamics modulate the scattering probabilities by up to 2% within the first diffraction order and up to 0.2% within the second one. Such values are sufficiently high to be experimentally accessible with state-of-the-art pump-probe electron diffraction 42,43 .
The time-dependent scattering probabilities shown in Fig. 3 have symmetry that is specific for our experimental geometry: → − k k Retrieval of time-dependent charge-density maps. Which details can be retrieved? Bragg-peak intensities represent the magnitudes of complex scattering amplitudes, whose phases are not measured. Many advanced techniques exist to circumvent this problem using prior information, but these techniques may require some adaptation to reveal spatially delocalised changes in electron density, which are in contrast to localised atomic-type scatterers in crystallography.
Therefore, we refrain here from algorithmic procedures and report on a reconstruction specific to the symmetry in our example. First, we assume that the initial microscopic scattering potential is known, for example from ground-state calculations or state-of-the-art crystallography 44,45 . Second, we assume that the scattering amplitudes maintain their initial phases during the interaction with the laser pulse; see Methods. In this way, we reconstruct the part of the time-dependent electronic motion that is symmetric with respect to the reflection → − y y, which is hidden from spectroscopy in dipole approximation. We used only the lowest two orders of Bragg scattering as depicted in the inset in Fig. 3. Figure 4 shows the results. In Fig. 4(a), we depict four representative times in the laser cycle and, additionally, a time ≈ t 20 fs max , after the pulse end. The electron densities are integrated over z because z-dependent information is absent in diffraction geometry at normal incidence. Figure 4(b) shows changes in the simulated electron density, Δ σ, and Fig. 4(c) shows the Δ σ (s) that is symmetric with respect to reflection → − y y, i.e. the part hidden from far-field spectroscopy. The outcomes of the reconstruction are shown in Fig. 4(d). The reconstructed charge-density maps are nearly identical to the y-symmetric parts of the actual electronic motion. Figure 4(e) shows the time-dependent charges in regions A and D, which are reliably reconstructed. A similar quality of reconstruction is also achieved for regions B and C (not shown). The spatial resolution can be further improved by incorporating more Bragg peaks in the reconstruction, provided that the signal-to-noise ratio is sufficiently large.
General value of the real-space approach. We conclude by returning to our initial conjecture on the importance of time-dependent real-space recording of electronic motion. A periodic solid can be fully treated in either reciprocal or real space. However, the simplicity of the description or the feasibility of an experiment can define a preferred perspective. In graphene, our real-space analysis has revealed that electrons in different spatial regions respond to an external field in distinctly different ways. In particular, current flows in stripes above and below the plane of the chemical bonds, absorption populates the hexagon centre, the electrons localised at the bonds are periodically displaced along their structural extensions, and charges are periodically transferred in and out of the conducting stripes. These spatial peculiarities of the light-induced electron dynamics suggest the real-space perspective to be a worthwhile concept exposing features that are not easily discernible otherwise.
Real-space measurements will probably contribute most when applied to problems where spectroscopy does not provide sufficient information. They may give new insights into electronic origins of light-induced structural deformations, for example in graphite 46 , vanadium dioxide 15 , charge-density-wave materials 47 or nanoparticles 42 . Knowing the mechanisms of charge transfer on an atomic scale may be indispensable for understanding strongly-correlated systems, i.e. materials where electron-electron interactions play a major role, as well as the build-up of electron-boson interaction, which has recently been shown to be an ultrafast process in high-temperature superconductors 48 . Identifying structure-function relations for optical phenomena can facilitate design of novel materials and ultrafast electronic devices. Ultrafast electron and X-ray diffraction bear a promise to become important techniques to study light-matter interaction in condensed matter.

Methods
The model. We designed a numerical model adapted for studying ultrafast electron diffraction on an atomically thin sample exposed to a strong field. We solve the time-dependent Schrödinger equation (TDSE) in three spatial dimensions for an ensemble of electrons moving in a local "effective" lattice potential V(r) that is periodic in any plane parallel to the atomic plane: V(r + a 1,2 ) = V(r). Here, a 1 and a 2 are the lattice unit vectors, both of which are orthogonal to e z . Due to this periodicity, the Bloch theorem applies. We introduce the in-plane crystal momentum k , which is a conserved quantity in the velocity gauge. Henceforth we use atomic units and write the TDSE as Here, A(t) is related to the external electric field by ′ ( ) = − ( ) F A t t. Equation (1) implies the dipole approximation; thus, we neglect the effect of the magnetic field.
Let Σ G denote a sum over vectors of the reciprocal lattice. For each z, we decompose both the wave function ψ k and the lattice potential in Fourier series: For a function f(r) that has the periodicity of the lattice, the decomposition coefficients are This equation is solved for a set of initial conditions, each of which is defined by k and an initial valence band n. One of the most important physical observables for electron diffraction is the electron density ρ ( , ) r t , which we decompose into plane waves as well: Equations (2) and (4) yield The integral here is taken over the first Brillouin zone (BZ), and the contributions from all initial bands are added together; ρ ( , ) z t G is guaranteed to be non-negative. Our phenomenological lattice potential V(r) is evaluated as a sum of atomic potentials , and ξ = . 0 54 3 . The band structure for this potential is shown in Fig. 5; it is similar to that obtained in ab initio simulations 49 . The Fermi velocity in our model is equal to = . × / v 0 98 10 m s F 6 in the direction orthogonal to Γ -K, which is close to the literature value: . × / 1 0 10 m s 6 50 .
We note that, unlike most pseudopotential methods, our lattice potential supports core states, which is important for electron diffraction. However, to reduce the required number of plane waves in expansion (2), we allow the 1s states to be less tightly bound as they are in a carbon atom: the 1s state in the potential V at (r) has an energy of − 121 eV, while the K-edge of carbon is at 282 eV. This allows us to consider only the reciprocal-space vectors with ≤ . G 8 256. While the presence of the 1s state in the atomic potential is essential for obtaining correct nodal structures in the 2s and 2p states, dynamics of inner-shell electrons in our calculations have a negligible effect on Bragg peaks. Therefore, solving the TDSE, we consider only those initial states that correspond to valence bands.
The laser pulse in our simulations has a central wavelength of λ π ω = / = . μ c 2 21 m L L , a full width at half maximum (FWHM) duration of 10 fs, and a peak field of = . F 0 5 0 V/Å, which corresponds to a peak intensity of × / 3 10 W cm 12 2 . We define the pulse by its vector potential: , a = 4.65 atomic units. Discretising the crystal momentum, we use a Monkhorst grid 51 with 32 × 32 = 1024 nodes. An 8-th order scheme allows us to evaluate ∂ z and ∂ z 2 on a relatively coarse z-grid: Δ = . z 0 2. The grid contains 256 nodes, and we apply periodic boundary conditions. For solving the TDSE, we use the 4-th order Magnus propagator 52,53 , evaluate matrix exponentials using Expokit 54 , and perform the required sparse-matrix operations with the aid of CXSparse 55 .
Electron diffraction. We are interested in the intensities of Bragg peaks; their positions do not change as long as the lattice motion can be neglected because the interaction with a laser field preserves lattice periodicity. The first step in modelling elastic scattering of energetic electrons is the evaluation of the scalar potential induced by atomic nuclei, as well as core-and valence-band (VB) electrons. The spatial Fourier transform of the total potential is given by  Fig. 1(d) is evaluated as