Ultrafast electronic response of graphene to a strong and localized electric field

The way conduction electrons respond to ultrafast external perturbations in low dimensional materials is at the core of the design of future devices for (opto)electronics, photodetection and spintronics. Highly charged ions provide a tool for probing the electronic response of solids to extremely strong electric fields localized down to nanometre-sized areas. With ion transmission times in the order of femtoseconds, we can directly probe the local electronic dynamics of an ultrathin foil on this timescale. Here we report on the ability of freestanding single layer graphene to provide tens of electrons for charge neutralization of a slow highly charged ion within a few femtoseconds. With values higher than 1012 A cm−2, the resulting local current density in graphene exceeds previously measured breakdown currents by three orders of magnitude. Surprisingly, the passing ion does not tear nanometre-sized holes into the single layer graphene. We use time-dependent density functional theory to gain insight into the multielectron dynamics.

. Raman spectrum and STEM measurement a) Raman spectrum of a graphene sample; the FWHM of ∼30 cm −1 of the 2D mode peak ensures single layer graphene and the ratio of the integrated intensities of the D and G mode peaks points to a negligible amount of defects b) STEM measurements (left) show clean monolayer graphene regions with a lateral extension of up to hundreds of nm, surrounded by graphene regions covered with thin hydrocarbon adsorbates. The monolayer nature of the clean graphene regions was confirmed using atomic-resolution STEM imaging (middle). Using the roughly linear intensity scaling of the HAADF STEM signal with specimen thickness, and assuming as an upper limit the same density for the hydrocarbon adsorbates as for graphene, a thickness estimated for the hydrocarbon contamination was calculated as equivalent thickness in number of graphene layers (right). This analysis shows that the hydrocarbon deposits on the graphene have a thickness corresponding to ∼3 graphene monolayers with some regions on the suspended graphene basal planes reaching up to ∼6 layers.
TDDFT allows to consistently address the one-electron and the multi-electron processes, as well as to follow the real time evolution of the electron densities, currents, and fields in the system. In a number of recent publications, for the case of the negative, neutral, single and double charged projectiles, TDDFT has been shown to provide an ab initio description of the stopping power of ions in solids [9][10][11][12] , as well as the energy loss in transmission through graphene [13][14][15][16][17] .
The specific details on the present implementation of the TDDFT can be found elsewhere 9,18 . In brief, the time-dependent electron density of the many-body system n(r,t) is represented via the density of a fictitious system of non-interacting electrons, the Kohn-Sham system, given by: with the summation running over all the occupied KS orbitals ψ j (r,t). The time evolution of the KS orbitals is governed by the time-dependent Scrödinger equation: Here, T = − 1 2 ∆ is the kinetic energy operator, and V eff [n(r,t)] is the effective KS potential that depends on the electron density. The Hartree V H [n(r,t)] and exchange-correlation V xc [n(r,t)] potentials are generally time-dependent through the time-dependence of the electron density. The exchange-correlation potential is described within the adiabatic local density approximation (ALDA) 6,8,19 using the exchange-correlation kernel of O. Gunnarsson and B. I. Lundqvist 20 . At this point it is worth to stress that the ALDA with exchange-correlation kernels developed so far in quantum chemistry does not allow to describe the Auger processes involving strong binary electronelectron interactions with large energy exchange. Thus, strictly speaking, the present approach is well suited for the incoming part of the multicharged ion trajectory, where a strongly excited hollow atom is formed via multielectron capture from graphene towards loosely bound orbitals of the projectile. The relaxation of the hollow atom via Auger cascades towards inner electronic shells is beyond the reach of the present description. We further address this point below in this Supplementary note.
The last term in Supplementary Equation 3 describes the electron interaction with the projectile ion core. We consider a fully stripped ion of charge q in impinging at graphene along the straight line trajectory R(t) = R 0 − vt perpendicular to the graphene layer. R(t) defines position ion core. Since the calculated energy loss is much smaller than the kinetic energy, the projectile speed |v| is set constant during the collision. With these assumptions V ion (r,t) = V ion (|r − R(t)|). To avoid the divergence of V ion at r = R(t), the positive charge q in of the ion core is homogeneously distributed over the small sphere of radius R ρ (q in ) given by (R ρ (q in )/ρ) 3 = q in . This is equivalent to the description of the projectile as a small highly charged quantum dot or, equivalently, a fully ionised jellium metal (JM) cluster 21,22 characterised by the density parameter ρ. We have used ρ = 1 a 0 and ρ = 0.5 a 0 (a 0 = 0.53Å is a Bohr radius) so that (in units of Bohr radius): For example, when q in = 20 the corresponding radii are R 1 (20) = 1.43Å and R 0.5 (20) = 0.72Å, respectively so that the projectile considered here is very compact. Within the Jellium cluster model, the V ion potential is given by: where ζ = |r − R(t)|. It follows from Supplementary Equation 5 that R 1/0.5 (q in ) corresponds to the cutoff parameter for the Coulomb potential of the HCI.
This simple picture of a HCI suffices to model the electron capture onto the outer Rydberg states of the projectile with hollow atom formation, as well as the energy loss and the electronic response of graphene, which is the purpose of this paper. The inner shell structure of the Xe projectiles is not described with this model consistent with impossibility to describe electronic relaxation to these shells with available exchange-correlation potentials.
In the TDDFT simulation of the HCI interaction with graphene one might use the supercell geometry or finite size of a graphene patch as has been done in the ab initio studies of the energy loss of the singly and doubly charged projectiles [14][15][16] . However, the peculiarity of the projectile-surface charge transfer involving the HCI is that because of the large attractive potential, the projectile starts to capture the electrons far away from the graphene sheet. The strongly polarized hollow atom is formed with electrons occupying diffuse Rydberg orbitals. This requires that the computational mesh is large enough to cover an extended region of space in the direction perpendicular to the surface, and in the lateral direction. The former is needed to correctly capture the onset of the electron transfer and the latter guaranties that the periodic boundary conditions or the finite size effects do not alter the results of the calculations. Along with this difficulty, considering the large number of electrons transferred to the projectile and emitted into vacuum by ionization because of the HCI impact, the fraction of the graphene surface per projectile has to be sufficiently large to avoid charging effects which would block the electron transfer. All together these constraints render the calculations with a full atomistic representation of graphene computationally extremely demanding, if feasible at all. We, therefore, decided to adopt the JM model to describe the graphene valence electrons.
where Θ is the Heaviside step function. In order to impose the correct work function of graphene Φ = 4.6 eV, 29, 30 a constant attractive potential V stab = 13.6 eV is added to the V H [n(r,t)] potential inside the jellium disk. This is the so-called "stabilized jellium model" 31 .
With the KS orbitals represented on a mesh in cylindrical coordinates, Supplementary Equation 2 are solved using real-time propagation similar to that developed for the one-electron wave packet propagation in cylindrical coordinates and detailed elsewhere [32][33][34] . To avoid complications linked with the change of the number of electrons in the computation box we do not introduce the complex absorbing potentials (CAPS) 35 at the boundaries of the mesh. Since the mesh used in our calculations is large, the overall density change due to the emitted electrons distributed over the computational box remains very small.
Obviously, within the JM the motion of graphene lattice atoms is not considered. We thus calculate the electronic contribution to the projectile energy loss, which has been shown to be a dominant energy loss channel for the projectile speeds considered here 13 . The energy loss of the projectile ∆E P is calculated integrating the action of the Hellmann-Feynman force F z over the trajectory The Hellmann-Feynman force is obtained from the time dependent electron density as where Z is the projectile z-coordinate. Observe that, because of the symmetry of the problem, the only non zero component of the force is along z-axis. Our scheme preserves the total energy of the quantum-classical where ε XC is the exchange-correlation energy density 20,37 . As a consequence, ∆E P can be alternatively obtained as At the beginning of the collision the projectile is assumed to be fully ionized, and the initial conditions for the time dependent Supplementary Equation 2, ψ 0 j (r) ≡ ψ j (r,t = 0), correspond to the KS orbitals of the ground-state graphene disk. Therefore, prior to the TDDFT calculations, the standard density functional theory (DFT) study has been performed. The KS orbitals ψ 0 j (r) satisfy the stationary Kohn Sham equations: is time independent. We use superscript 0 to refer to the ground state of the system. Within the JM the electronic structure of the free-standing graphene layer is characterised by the two continua of states. The electron motion is quantized in z-direction perpendicular to the layer. The electron motion parallel to the layer is free, and it is characterized by the wave vector k || . The energies of the electronic states are given by: where E 1 = −31.95 eV and E 2 = −13.68 eV. These bands mimic the σ -and π-band of graphene, respectively. The electronic states ψ 0 j (r) of the model graphene disk originate from the quantization of these two free electron continua by the finite lateral size.
A remark is in order here concerning the JM used for graphene. Taking into account an extremely strong Coulomb potential of the ingoing projectile, the electronic structure of graphene is locally modified. Thus, even though the JM does not reproduce the details of the stationary band structure of pristine graphene including the Dirac cone at K-point, we believe that it allows to capture the essential physics of the dynamics of the HCI-surface charge transfer and energy loss.
To test the effect of the JM approximation on the actually measured quantities we performed calculation of the energy loss of He ++ projectiles and compared it with available TDDFT results obtained recently 6/14 with full atomistic description of the graphene layer 15,16,38 as we show in Supplementary Figure 2. In overall, the projectile energy dependence of the electronic stopping obtained in full atomistic TDDFT study is well reproduced within our model approach. In particular, we find the correct energy position of the energy loss maximum as well as the overall shape of the energy dependence. However, as follows from the results presented with open circles in Supplementary Figure 2, the JM systematically underestimates the energy loss value by ≈ 55 eV. We attribute the lower energy loss calculated with the present model to the higher values of the binding energies |E 1 | = 31.95 eV and |E 2 | = 13.68 eV states as compared to the σ (|E σ | = 20 eV) and π (|E π | = 7 eV) bands of graphene 16 . In particular, at low projectile velocities when weakly bound π electrons dominate the electronic stopping, this effect is more significant.
These arguments are further supported by the analysis of the final charge state of the He-projectile after the collision with graphene layer at normal incidence shown in Supplementary Figure 3. The mean number of electrons captured by the He ++ incident ion is determined on the exit part of the trajectory after the projectile has traversed the graphene layer. It is obtained integrating the electron density within the sphere of 2.65Å around the projectile. Even though this procedure is focused on the ground state population analysis and does not allow to account for the occupation of the diffuse excited orbitals, we apply it here for consistency with 16 . While at high projectile energies the JM and full atomistic TDDFT calculations agree quantitatively, qualitative differences emerge below 20 keV/amu. A decrease in the probability of electron capture with lowering energy is obtained in full atomistic calculations, while the present JM shows larger neutralisation rates at lower energies. This result is consistent with the higher binding energies of the electronic states of graphene in the present model. For example, the binding energy of the E 1 (k) band is resonant with the first ionization potential of He allowing the resonant electron capture by He + projectile, while the σ and π bands of graphene have lower binding energies and electron capture into inner projectile levels is then a nonadiabatic process requiring high enough projectile speed.
TDDFT Results for the highly charged ions.
The number of electrons captured by the HCIs N cap is shown in the left panel of Supplementary Figure 4 for the q in =10, 20 and 40 projectiles incident at the graphene layer represented with 1004 electrons jellium disk. N cap is determined on the outgoing trajectory path after the projectile traversed the target. It is obtained by integrating the electron density within the sphere centered at the projectile. The radius of the sphere is set as 12.5Å, i.e. it is sufficiently large to account for the occupation of excited states. Within the present approach, N cap appears to be quite independent of the projectile velocity. Consistent with experimental observations we obtain that the HCI with initial charge q in is almost neutralized in all cases with residual charge as small as q in /10. In the right panel of Supplementary Figure 4 we show the total number of electrons lost by the target N lost (ionized electrons). It is approximately twice larger than N cap so that N vac = N lost − N cap N cap electrons are emitted into the free-electron continuum upon impact. The slower is the projectile, the more adiabatic is the interaction. The electron emission N vac then decreases leading to the overall decrease of the number of ionized electrons for ion velocities lower than 0.44 nm/fs.
In the experiments, slower projectiles capture more electrons leading to the velocity dependent N cap . Indeed, the slow collision allows full relaxation of the projectile into the ground state via Auger processes. As discussed above, the present theoretical approach does not allow to include Auger processes. It is aimed at the description of the transfer of valence electrons between the HCI and graphene.
In Supplementary Figure 5 we analyse the dependence of the calculated energy losses on the model used to describe the HCI projectile, and on the initial charge of the projectile q in . As follows from the results shown in Supplementary Figure 5, reducing the cutoff parameters for Coulomb potential from R 1 (q in ) to R 0.5 (q in ) i.e. representing the projectile core as more compact, leads to the overall increase of the calculated energy loss. For the present case we calculated that a reduction of cutoff radius by a factor of two leads to the twice larger energy losses. This increase can be even more prominent at higher velocities and charge states. We argue that, together with higher binding energies of the graphene bands within the present jellium model, the sensitivity to the Coulomb potential cutoff is at the origin of the underestimation of the measured energy losses by present calculations as shown in Figure 3 of the paper. The dependence of the projectile energy loss ∆E P on the incident charge state q in is shown in right panel of Supplementary Figure 5. For the fixed model of the projectile, the quadratic dependence of ∆E P ∼ q 2 in is found in full accord with experimental data. The target size dependence of the calculated energy loss ∆E P and the number of captured electrons N cap is analysed in Supplementary Figure 6. The lateral size of the jellium disk D was varied by a factor of two (variation of the number of electrons by a factor of four from 500 to 2000). The results presented in Supplementary Figure 6 show that the calculations converged with respect to the finite size effects already for a 500 electrons jellium disk as far as the projectiles with q in = 10 and 20 incident charge states are concerned. For the highest incident charge state considered here, q in = 40, finite size effects are still present, although they do not affect the main conclusions of this work.
The convergence of the energy loss with respect to the finite size effect reflects convergence of the force exerted on the HCI, as can be seen in the left panel of Supplementary Figure 7. This is particularly the case in the incoming path of the trajectory that predominantly determines the energy loss (the ion is basically neutralized after crossing the graphene layer). Additionally, it is worth to mention that energy conservation