K. Balzer et al

*

# Neutralization dynamics of slow highly charged ions passing through graphene
nanoflakes–an embedding self-energy approach

K. Balzer  M. Bonitz  Rechenzentrum, Christian-Albrechts-Universität zu Kiel,
24098 Kiel, Germany  Institut für Theoretische Physik und Astrophysik,
Christian-Albrechts-Universität zu Kiel, Leibnizstraße 15, 24098 Kiel, Germany
<EMAIL_ADDRESS>
(xxx; xxx; xxx)

###### Abstract

We study the time-dependent neutralization of a slow highly charged ion that
penetrates a hexagonal hollow-centred graphene nanoflake. To compute the
ultrafast charge transfer dynamics, we apply an effective Hubbard nanocluster
model and use the method of nonequilibrium Green functions (NEGF) in
conjunction with an embedding self-energy scheme which allows one to follow
the temporal changes of the number of electrons in the nanoflake. We perform
extensive simulations of the charge transfer dynamics for a broad range of ion
charge states and impact velocities. The results are used to put forward a
simple semi-analytical model of the neutralization dynamics that is in very
good agreement with transmission experiments, in which highly charged xenon
ions pass through sheets of single-layer graphene.

###### keywords:

Highly charged ions, charge transfer, neutralization dynamics, single-layer
graphene, nonequilibrium Green functions, embedding self-energy

††articletype: ORIGINAL ARTICLE

##  1 Introduction

The interaction of heavy charged particles with matter is of fundamental
importance in various fields of physics. Particularly, the collision of ions
with surfaces of solids is of high current interest both, theoretically and
for technological applications in plasmas and plasma-surface interactions 1.
The ion-induced excitation and ionization of electrons in the surface material
as well as the capture of electrons by the ion (which changes the charge state
of the ion) is crucial for the electrostatic potential at the surface and in
the plasma sheath region. The theoretical description of charge transfer
processes in the plasma near the surface region has been the subject of
intense investigations, for recent studies and further references, see, e.g.,
Refs. 2, 3. In recent years, a number of studies was devoted to the energy
loss (stopping power) of charged particles in strongly correlated solids,
nanostructures and single-layered materials, such as graphene 4, 5, 6.

Another topic that caught attention has been the irradiation of solids by
highly charged ions (HCIs) 7, 8, in particular, the impact of HCIs on single-
layer graphene 9. Due to their high charge these ions carry a large amount of
potential energy 7, 8 and, depending on their kinetic energy, reveal manifold
and complex interactions with surfaces of solids. For slow HCIs with typical
velocities of $v\sim\\!10^{5}\,$m/s, i.e., with velocities less than the Bohr
velocity [the electron velocity in the first Bohr orbit], surface transmission
times are of the order of a few femtoseconds and are hence comparable to the
time scale of the electronic motion in the solid. In this regime, the exposure
of surfaces to individual or ensembles of highly charged ions allows for
interesting nanoscale applications including surface (re)structuring and
modification 10, 11, 8, 12 and surface diagnostics 13. Moreover, upon ion
impact, the extraction of electrons from the surface can lead to partial or
complete neutralization of the ion 9 as well as to secondary electron emission
14. The loss of charge carriers in the surface, together with the deposition
of potential energy, can ultimately culminate in the destruction of surface
domains in a Coulomb explosion-like sputtering scenario 15. How exactly a
surface and its electronic structure reacts to the impact of a highly charged
ion depends on multiple factors, in particular on the ion’s charge state, the
ion velocity and the material properties. Interesting questions, which have
recently been addressed in that regard, concern, for example, the formation of
correlated two-particle states (doublons) during ion impact 6, 16, 17 or the
influence of the band gap on the charge exchange 18. Of interest is also the
impact of multiple ions and possible correlations of their individual effects.
This requires to go beyond linear response theory and to perform a non-
adiabatic time-dependent correlated simulation, see Refs. 6, 17. In these
works a combined Nonequilibrium Green Functions(NEGF)-Ehrenfest approach19, 20
was developed 16, 5 that focused on ion stopping and energy loss spectra
(electronic stopping power) for closed surface models. However, the processes
related to charge transfer were not included in these simulations.

The goal of the present article, is to extend these time-resolved NEGF-
Ehrenfest simulations to include charge transfer processes. As a test case, we
simulate the electron emission from two-dimensional hexagonal hollow-centered
graphene nanoflakes (GNFs) during the impact of slow HCIs. In the course of
this, neutralizing charge transfer to the ion is described via an embedding
self-energy scheme, which mimics the formation of a highly excited hollow atom
in front of the GNF surface. The approach is capable to resolve the ultrafast
electronic response of the nanoflake at accurate exit charge states of the
ion, while the de-excitation and recombination processes inside the ion are
not explicitly considered. The article is organized as follows. Section 2
outlines the effective Hubbard-type model that is used to describe the
graphene nanoflake, including electronic correlations, and its interaction
with the ion. Further, we here motivate the spatio-temporal form of the charge
transfer amplitude. Section 3 presents the theory, derives the central
Keldysh-Kadanoff-Baym equations of the GNF and introduces the embedding self-
energy, which couples ion and surface states. Results of the NEGF calculations
are compiled in Sec. 4, and the conclusions are summarized in Sec. 5. For the
system of units used throughout the article, see below Eq. (2).

##  2 Model

Figure 1: (a) Hexagonal hollow-centered graphene nanoflakes (GNF) with a
finite number of honeycombs and $L=24$ and $54$ lattice sites. The red dot
($\odot$) in the center of the GNFs indicates the impact point of the
approaching ion. (b) Charge transfer couplings considered in the present
model, Eq. (2). Depending on the initial ion charge state $Z_{0}$, only a
selected set of central GNF sites $k\in M$ participate in the resonant
transfer of electrons at an effective single-particle energy
$\epsilon_{k}=\epsilon$.

The electrons in the covalent $\sigma$\- and $\pi$-bonds of the sp2-hybridized
carbon atoms of the GNF are described in terms of an effective Hubbard
nanocluster model with four equal decoupled bands at half-filling, cf. Fig.
1a for the lattice structure of GNFs with $L=24$ and $54$ sites. Using a
nearest-neighbor hopping $J$, a local Coulomb repulsion $U$ and an on-site
energy $E$, the overall Hamiltonian is given by
$\hat{H}(t)=4\hat{H}_{\textup{s}}(t)$, with the single-band contribution

$\displaystyle\hat{H}_{\textup{s}}(t)=$
$\displaystyle\sum_{i\sigma}\underbrace{\left\\{W_{i}(Z_{0};{\bf
S}(t))+E\right\\}}_{=\,E_{i}({\bf
S}(t))}\,\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}+\sum_{\langle
ij\rangle,\sigma}\left\\{W_{\langle ij\rangle}(Z_{0};{\bf
S}(t))-J\right\\}\,\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}+U\sum_{i}\left(\hat{n}_{i\uparrow}-\frac{1}{2}\right)\left(\hat{n}_{i\downarrow}-\frac{1}{2}\right)$
$\displaystyle+\sum_{k\in M,\sigma}\gamma_{k}[{\bf
S}(t)]\left(\hat{a}_{k\sigma}^{\dagger}\hat{c}_{k\sigma}+\hat{c}_{k\sigma}^{\dagger}\hat{a}_{k\sigma}\right)+\sum_{k\in
M,\sigma}\epsilon_{k}\,\hat{a}_{k\sigma}^{\dagger}\hat{a}_{k\sigma}\,.$ (1)

The operator $\hat{c}^{\dagger}_{i\sigma}$ ($\hat{c}_{i\sigma}$) creates
(annihilates) an electron with spin $\sigma\in\\{\uparrow,\downarrow\\}$ on
site $i$ of the GNF oriented in the $x-y$-plane,
$\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}$ denotes the
local electron density, $\langle ij\rangle$ indicates nearest-neighbor
contributions, and the matrix elements $W_{i}(Z_{0};{\bf S}(t))$ and
$W_{\langle ij\rangle}(Z_{0};{\bf S}(t))$ describe the time-dependent Coulomb
interaction of the electrons with a highly charged ion, which has initial
charge state $Z_{0}$ and follows a normal incidence trajectory ${\bf
S}(t)=(0,0,z(t)=v_{z}t)$ with a constant velocity $v_{z}$, cf. Sec. 2.1.
Furthermore, the last two terms in Eq. (2) allow for charge transfer from a
set of selected lattice sites $k\in M$ onto initially unoccupied energy levels
$\epsilon_{k}$ of the external ion, where $a^{\dagger}_{k\sigma}$
($a_{k\sigma}$) are the corresponding creation (annihilation) operators, cf.
Sec. 2.2 for the discussion of the trajectory-dependent charge transfer
amplitude $\gamma_{k}[{\bf S}(t)]$ and the energies $\epsilon_{k}$. Depending
on the charge state $Z_{0}$, we couple per Hubbard band only a few of the
central GNF sites to the ion, cf. Fig. 1b. For a more detailed discussion
on this choice, see end of Sec. 3.

Throughout the article, we measure energies in units of
$J_{0}=\hbar^{2}/(ma_{0}^{2})\approx 3.78$ eV and times in units of
$t_{0}=\hbar/J_{0}\approx 0.17$ fs, where $m$ denotes the electron mass, and
$a_{0}=1.42$ Å is the molecular carbon-carbon bond length. The ion velocity
$v_{z}$ is given in units of $v_{0}=a_{0}/t_{0}\approx 0.82$ nm/fs$=8.2\cdot
10^{5}$ m/s.

###  2.1 Kinetic, potential, and interaction energy

To determine the single-particle matrix elements of the kinetic and potential
energies in Eq. (2), we use localized atomic orbitals of Gaussian form,
$\varphi_{i\sigma}({\bf r})=\frac{1}{\pi^{3/4}\mu^{3/2}}\textup{e}^{-({\bf
r}-{\bf R}_{i})^{2}/(2\mu^{2})}$, with width $\mu=\mu_{0}a_{0}$, and assume
the potential close to a lattice site with coordinate ${\bf R}_{i}$ to be
$v_{i}({\bf r})=-V\textup{e}^{-({\bf r}-{\bf R}_{i})^{2}/(2\nu^{2})}$, with
$\nu=\nu_{0}a_{0}$ and $V=V_{0}J_{0}$. Below we will use the parameters
$\mu_{0}$, $\nu_{0}$ and $V_{0}$ as free parameters of the model, their values
for the case of graphene are specified in Eqs. (9) and (10).

For sites located at positions ${\bf R}_{i}=(0,0,0)$ and ${\bf
R}_{j}=(a_{0},0,0)$, the kinetic part, $T_{\langle ij\rangle}$, of the
nearest-neighbor hopping matrix element, $-J=T_{\langle ij\rangle}+V_{\langle
ij\rangle}$, is given by

$\displaystyle T_{\langle ij\rangle}$
$\displaystyle=-\frac{\hbar^{2}}{2m}\int\textup{d}^{3}r\,\phi^{*}_{i\sigma}({\bf
r})\nabla^{2}\phi_{j\sigma}({\bf
r})=-\frac{1}{8\mu_{0}^{4}}(1-6\mu_{0}^{2})\,\textup{e}^{-1/(4\mu^{2})}J_{0}\,,$
(2)

whereas, for the potential part, $V_{\langle ij\rangle}$, we obtain

$\displaystyle V_{\langle ij\rangle}$
$\displaystyle=\int\textup{d}^{3}r\,\phi^{*}_{i\sigma}({\bf
r})\left\\{v_{i}({\bf r})+v_{j}({\bf r})\right\\}\phi_{j\sigma}({\bf
r})=-\frac{4\sqrt{2}\,\nu_{0}^{3}\textup{e}^{-(\mu_{0}^{2}+\nu_{0}^{2})/(2\mu_{0}^{2}(\mu_{0}^{2}+2\nu_{0}^{2}))}V_{0}}{(\mu_{0}^{2}+2\nu_{0}^{2})^{3/2}}J_{0}\,.$
(3)

For the on-site energy $E=T_{i}+V_{i}$, we have

$\displaystyle T_{i}$
$\displaystyle=-\frac{\hbar^{2}}{2m}\int\textup{d}^{3}r\,\phi^{*}_{i\sigma}({\bf
r})\nabla^{2}\phi_{i\sigma}({\bf r})=\frac{3}{4\mu_{0}^{2}}J_{0}\,,$ (4)
$\displaystyle V_{i}$
$\displaystyle=\int\textup{d}^{3}r\,\phi^{*}_{i\sigma}({\bf r})v_{i}({\bf
r})\phi_{i\sigma}({\bf
r})=-\frac{2\sqrt{2}\nu_{0}^{3}V_{0}}{(\mu_{0}^{2}+2\nu_{0}^{2})^{3/2}}J_{0}\,.$
(5)

Furthermore, the time-dependent Coulomb potential generated by the ion,

$\displaystyle w[Z_{0};{\bf r}-{\bf S}(t)]$
$\displaystyle=-\frac{Z_{0}e^{2}}{4\pi\epsilon_{0}}|{\bf r}-{\bf
S}(t)|^{-1}=-Z_{0}W_{0}a_{0}|{\bf r}-{\bf S}(t)|^{-1}\;,$ (6)

with electron charge $e$, vacuum permittivity $\epsilon_{0}$ and
$W_{0}=e^{2}/(4\pi\epsilon_{0}a_{0})\approx 10.14$ eV ($W_{0}/J_{0}\approx
2.68$), leads to a modification of the on-site energy and the nearest-neighbor
hopping. The corresponding matrix elements are treated as follows:

$\displaystyle W_{i}[Z_{0};{\mathbf{S}}(t)]$
$\displaystyle=\int\textup{d}^{3}r\,\phi^{*}_{i\sigma}({\bf r})w[Z_{0};{\bf
r}-{\bf S}(t)]\phi_{i\sigma}({\bf r})\approx w[Z_{0};{\bf R}_{i}-{\bf
S}(t)]\,,$ (7) $\displaystyle W_{\langle ij\rangle}[Z_{0};{\mathbf{S}}(t)]$
$\displaystyle=\int\textup{d}^{3}r\,\phi^{*}_{i\sigma}({\bf r})w[Z_{0};{\bf
r}-{\bf S}(t)]\phi_{j\sigma}({\bf
r})\approx\frac{\lambda}{2}\left(w[Z_{0};{\bf R}_{i}-{\bf S}(t)]+w[Z_{0};{\bf
R}_{j}-{\bf S}(t)]\right)\,,$ (8)

where $\lambda$ denotes the overlap integral $\lambda=\int_{\langle
ij\rangle}\textup{d}^{3}r\,\phi^{*}_{i\sigma}({\bf r})\phi_{j\sigma}({\bf
r})=\textup{e}^{-1/(4\mu_{0}^{2})}$, defined with respect to nearest neighbor
sites.

To determine the free parameters $\mu_{0}$, $\nu_{0}$ and $V_{0}$ in Eqs.
(2)-(5), we require

$\displaystyle-J=T_{\langle ij\rangle}+V_{\langle ij\rangle}$
$\displaystyle=-2.8\,\textup{eV},$ (9) $\displaystyle
E=T_{i}+V_{i}=-W_{\textup{f}}$ $\displaystyle=-4.6\,\textup{eV},$ (10)

where $J=2.8$ eV ($J/J_{0}\approx 0.74$) is the typical carbon-carbon hopping
energy in graphene 21, and $W_{\textup{f}}=4.6$ eV
($W_{\textup{f}}/J_{0}\approx 1.22$) is the work function of single-layer
graphene 22. This leads to a nonlinear system of equations, which has a
reasonable solution for $\mu_{0}=0.365$, $\nu_{0}=\lambda=0.153$ and
$V_{0}=51.274$. Finally, following Ref. 21, we set the on-site Coulomb
repulsion $U$ in Hamiltonian (2) to $U=1.6J=4.48$ eV ($U/J_{0}\approx 1.19$)
and, neglecting electron-electron correlations, use the Hartree approximation,

$\displaystyle
U\left(\hat{n}_{i\uparrow}-\frac{1}{2}\right)\left(\hat{n}_{i\downarrow}-\frac{1}{2}\right)\approx
U\sum_{\sigma}\left(\langle\hat{n}_{i\bar{\sigma}}\rangle_{\textup{s}}(t)-\frac{1}{2}\right)\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}\,,$
(11)

where
$\langle\hat{n}_{i\uparrow}\rangle_{\textup{s}}(t)=\langle\hat{n}_{i\downarrow}\rangle_{\textup{s}}(t)$
denotes the expectation value of the local electron density with respect to
$\hat{H}_{\textup{s}}(t)$. The computational approach outlined in Sec. 3 is,
however, not limited to a mean-field treatment and can be straightforwardly
extended to include electron-electron correlations via the many-body self-
energy, cf. Eq. (35) in Sec. 3.

###  2.2 Charge transfer amplitude

To get an idea of how the charge transfer amplitude $\gamma_{k}[\bf{S}(t)]$ of
Eq. (2) varies as function of time and ion-surface distance $z(t)$,
respectively, we start from the energy spectrum of a hydrogen-like ion with
charge (atomic number) $Z$, which is given by

$\displaystyle E_{n}(Z)=-E_{\textup{R}}\frac{Z^{2}}{n^{2}}\,,$ (12)

with Rydberg energy $E_{\textup{R}}=1$ Ry $\approx 13.61$ eV
($E_{\textup{R}}\approx 3.60J_{0}$) and $n=1,2,3,\ldots,\infty$. Equation (12)
represents a very simple approximation for the energy level structure of the
highly excited states of the impacting ion as function of the charge state
$Z$. The quantity $\textup{d}n/\textup{d}E$ is used in the following to
characterize the number of final states that can be reached by an emitted
electron during the neutralizing charge transfer. The number of final states
$dn$ per energy interval $dE$ is obtained by inversion of Eq. (12),

$\displaystyle\frac{\textup{d}n}{\textup{d}E}(Z,E)=\frac{Z}{2}\left(-\frac{E_{\textup{R}}}{E}\right)^{-1/2}\\!\left(\frac{E_{\textup{R}}}{E^{2}}\right)\,.$
(13)

For an electron with energy $E_{k}$ on site $k\in M$ of the graphene
nanoflake, the probability for resonant charge transfer onto the ion will then
be proportional to Eq. (13) and proportional to some distance-dependent matrix
element $\Omega_{k}$. Parametrizing all quantities via the ion trajectory,
$\textbf{S}(t)$, we define

$\displaystyle\gamma_{k}[\mathbf{S}(t)]=\gamma\,\cdot\,Z[\mathbf{S}(t)]\,\cdot\,\underbrace{\frac{1}{2}\left(-\frac{E_{\textup{R}}}{E_{k}}\right)^{-1/2}\\!\left(\frac{E_{\textup{R}}}{E_{k}^{2}}\right)E_{\textup{R}}}_{\displaystyle=\frac{\textup{d}n}{\textup{d}E}(Z,E_{k})\,\frac{E_{\textup{R}}}{Z}}\,\cdot\,\Omega_{k}[\mathbf{R}_{k}-\mathbf{S}(t)]\,,$
(14)

where
$E_{k}[\mathbf{S}(t)]=w[Z_{0};\mathbf{R}_{k}-\mathbf{S}(t)]-W_{\textup{f}}\,$,
(cf. Sec 2.1), and $\gamma=\gamma_{0}J_{0}$ is an overall proportionality
factor. Equation (14) describes the time- and distance-dependent probability
for an electron to undergo a resonant transition from a site $k$ of the
graphene nanoflake to the ion (or vice versa) and thereby contains information
on, both, the nature of the highly charged ion, as well as of the electronic
structure of the initial and final states via the overlap matrix element
$\Omega_{k}$.

From experiments it is known that the charge transfer amplitude (14) is peaked
in front of the surface, at a characteristic ion-surface distance. For this
reason, we choose an analytical form for $\gamma_{k}[\mathbf{S}(t)]$ which, on
the one hand, directly includes such a maximum and, on the other hand, is
still sufficiently general. A reasonable choice is a Gaussian which, aside
from the prefactor $\gamma_{0}$, involves only two further open parameters,
$z_{\textup{res}}$ and $d_{\textup{w}}$,

$\displaystyle\gamma_{k}[{\mathbf{S}}(t)]=\gamma_{0}\exp\\!\left(-\frac{(z(t)+z_{\textup{res}})^{2}}{2d_{\textup{w}}^{2}}\right)J_{0}\,,$
(15)

where $z_{\textup{res}}$ marks the point of resonance, i.e., the ion-surface
distance (point in time) where the charge transfer probability has its
maximum, and $d_{\textup{w}}$ sets a characteristic spatial (temporal) width
over which the charge transfer probability is nonzero. In the following, we
discuss how the parameters $z_{\textup{res}}$ and $d_{\textup{w}}$ can be
determined by comparing Eqs. (14) and (15). To this end, we have to make
reasonable assumptions for the dependencies of the charge state $Z$ as well as
the matrix element $\Omega_{k}$ in Eq. (14) on the ionic position (i.e. on
time). As a further simplification, we will take into account charge transfer
only from lattice sites which belong to the GNF’s innermost honeycomb and have
thus the lateral distance $a_{0}$ to the impact point of the ion, cf. Fig.
1. For the time dependence of the ionic charge state, $Z$, we use a
functional form which in a step-like but continuous behavior drops off from
the initial charge state $Z_{0}$ and approaches zero before the ion passes
though the graphene layer,

$\displaystyle Z[{\bf
S}(t)]=\frac{Z_{0}}{\textup{e}^{[z(t)+z^{*}]/\eta}+1}\,,$ (16)

which leads to an approximately constant neutralization rate around
$z(t)=z^{*}$ and ensures a smooth behavior around the onset of neutralization
as well as during the termination of the charge transfer. The parameters
$\eta=\eta_{0}a_{0}$ and $z^{*}=z_{0}^{*}a_{0}$ in Eq. (16) are chosen such
that neutralization starts visibly around a few lattice spacings
($\sim\\!5$-$7$ Å), in front of the GNF surface, and terminates at around
$z=-a_{0}$, which is motivated by the experimental results of Ref. 9. Further,
for the overlap matrix element, $\Omega_{k^{*}}$, where $k^{*}$ indicates a
site index on the central honeycomb, a reasonable choice is again a Gaussian,

$\displaystyle\Omega_{k^{*}}[{\bf R}_{k}-{\bf
S}(t)]=\textup{e}^{-z^{2}(t)/(4\zeta^{2})}\,,$ (17)

with $\zeta=\zeta_{0}a_{0}$, which exhibits a maximum when the ion penetrates
the $x$-$y$-plane of the graphene nanoflake.

In Fig. 2, we show the charge transfer amplitude $\gamma_{k^{*}}(z)$
[orange curves] as defined in Eq. (14) with the input of Eqs. (16) [red
curves] and (17) [green curves] under variation of the parameters $\zeta_{0}$
[panel (a)] and $Z_{0}$ [panel (b)] and compare the obtained form as function
of the ion-surface distance $z$ with the ansatz of Eq. (15). From both panels,
we observe that the charge transfer amplitude is indeed sharply peaked at some
characteristic ion-surface distance and only for large values of $\zeta_{0}$
extends to considerably larger ion-surface distances. Furthermore, we find
that the position of the peak in $\gamma_{k^{*}}$ depends only very weakly on
the ion’s initial charge state. The resonance point
$z_{\textup{res}}=-\sqrt{3}a_{0}$ thereby describes the peak position
reasonably well, cf. the black vertical lines in Fig. 2, and thus
determines the first open parameter in Eq. (15). To also fix the parameter
$d_{\textup{w}}$, we reproduce the width of the main peak in the orange curves
for $\gamma_{k^{*}}(z)$ and set $d_{\textup{w}}=0.6a_{0}$. The resulting form
of Eq. (15) is hence the black curve in Fig. 2a and 2b, which shows Eq.
(15) for $\gamma_{0}=1$ (scaled by a factor of $\times$ $0.65$).

Finally, as we focus on resonant charge transfer processes, the point
$z_{\textup{res}}=v_{z}t_{\textup{res}}$ also marks a specific GNF on-site
energy,
$E_{k^{*}}[\mathbf{S}(t_{\textup{res}})]=w[Z_{0};\mathbf{R}_{k^{*}}-\mathbf{S}(t_{\textup{res}})]-W_{\textup{f}}$,
which, in turn, provides a reasonable choice for the energy
$\epsilon_{k^{*}}$. The resonance condition is given by

$\displaystyle\epsilon_{k^{*}}\equiv
E_{k^{*}}[\mathbf{S}(t_{\textup{res}})]=-\frac{Z_{0}W_{0}}{\left(1+z_{\textup{res}}^{2}/a_{0}^{2}\right)^{1/2}}-W_{\textup{f}}\,,$
(18)

and, for the resonance point $z_{\textup{res}}=-\sqrt{3}a_{0}$, the resonance
condition (18) yields half the maximum induced field strength, i.e.,
$\epsilon_{k}=\epsilon=-Z_{0}W_{0}/2-W_{\textup{f}}$. With the above choice of
parameters, Eq. (15) leaves open only one parameter, the prefactor
$\gamma_{0}$, which will be used in Sec. 4 to adapt the obtained results to
xenon ion transmission experiments 9.

Figure 2: Analysis of the charge transfer amplitude (orange curves).
Evaluation of Eq. (14) for ions with different initial charge states $Z_{0}$
(as indicated in the figure), using the assumptions (16) and (17) with
$\eta_{0}=0.3$ and $z_{0}^{*}=-2.5$. In panel (a), the parameter $\zeta_{0}$
is varied between $\zeta_{0}=0.9$ and $\zeta_{0}=2.1$ in steps of $0.3$,
whereas in panel (b), we fix $\zeta_{0}=1.5$ but vary the initial charge
$Z_{0}$. In both panels, the red, blue and green lines are assigned to the
right ordinate and show the individual factors in Eq. (14). Furthermore, the
black curve marks the charge transfer amplitude as defined in Eq. (15), with
$d_{\textup{w}}=0.6a_{0}$ and $z_{\textup{res}}=-\sqrt{3}\,a_{0}$ as distance
of resonance, where the amplitude reaches its maximum.

##  3 Embedding self-energy approach

To compute the electron dynamics resulting from Hamiltonian (2), we use the
method of nonequilibrium Green functions (NEGF), see, e.g., Refs. 19, 23. The
central quantity, from which all relevant observables can be obtained, is the
one-particle nonequilibrium Green function of the graphene nanoflake,

$\displaystyle G^{c}_{ij\sigma}(t,t^{\prime})=-\frac{\textup{i}}{\hbar}\langle
T_{\cal
C}\,\hat{c}_{i\sigma}(t)\,\hat{c}^{\dagger}_{j\sigma}(t^{\prime})\rangle_{\textup{s}}\,,$
(19)

where $T_{\cal C}$ denotes time ordering on the Keldysh contour ${\cal C}$,
and $\langle T_{\cal C}\dots\rangle_{\textup{s}}=\textup{tr}[T_{\cal
C}\exp(A_{\textup{s}})\dots]/\textup{tr}[T_{\cal C}\exp(A_{\textup{s}})]$
indicates the ensemble average with
$A_{\textup{s}}=-\frac{\textup{i}}{\hbar}\int_{\cal
C}\textup{d}\bar{t}\,\hat{H}_{\textup{s}}(\bar{t})$. The density matrix and
the site occupations of the GNF are recovered from the NEGF (19) in the limit
of equal times,

$\displaystyle\langle\hat{\rho}_{ij\sigma}\rangle_{\textup{s}}(t)=\langle\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}\rangle_{\textup{s}}(t)$
$\displaystyle=-\textup{i}\hbar\,G^{c}_{ji\sigma}(t,t^{+})\,,$ (20)
$\displaystyle\langle\hat{n}_{i\sigma}\rangle_{\textup{s}}(t)$
$\displaystyle=-\textup{i}\hbar\,G^{c}_{ii\sigma}(t,t^{+})\,,$

where the notation $t^{+}$ means, that the time $t^{+}$ is infinitesimally
larger along the contour ${\cal C}$ than $t$. Further, the overall particle
number is given by

$\displaystyle\langle\hat{N}\rangle_{\textup{s}}(t)=-\textup{i}\hbar\sum_{i\sigma}G^{c}_{ii\sigma}(t,t^{+})\,,$
(21)

which needs to be multiplied by a factor of $4$ to obtain the total number of
electrons in all four Hubbard bands of the present GNF model, cf. Sec. 2.
Similar to the NEGF of the graphene nanoflake, we define the charge transfer
Green functions,

$\displaystyle G^{ca}_{ik\sigma}(t,t^{\prime})$
$\displaystyle=-\frac{\textup{i}}{\hbar}\langle T_{\cal
C}\,\hat{c}_{i\sigma}(t)\,\hat{a}^{\dagger}_{k\sigma}(t^{\prime})\rangle_{\textup{s}}\,,$
$\displaystyle G^{ac}_{kj\sigma}(t,t^{\prime})$
$\displaystyle=-\frac{\textup{i}}{\hbar}\langle T_{\cal
C}\,\hat{a}_{k\sigma}(t)\,\hat{c}^{\dagger}_{j\sigma}(t^{\prime})\rangle_{\textup{s}}\,,$
(22)

which couple the GNF sites to the energy levels $\epsilon_{k}$ of the external
ion (and vice versa). Finally, the NEGF of the ion is

$\displaystyle G^{a}_{kk^{\prime}\sigma}(t,t^{\prime})$
$\displaystyle=-\frac{\textup{i}}{\hbar}\langle T_{\cal
C}\,\hat{a}_{k\sigma}(t)\,\hat{a}^{\dagger}_{k^{\prime}\sigma}(t^{\prime})\rangle_{\textup{s}}\,,$
(23)

and the superscripts indicate the type and arrangement of the creation and
annihilation operators.

The equations of motion for the one-particle nonequilibrium Green functions
(19), (22) and (23) are the Keldysh-Kadanoff-Baym equations (KBE) 24, 23:

$\displaystyle\sum_{l}\left[\textup{i}\hbar\,\partial_{t}\delta_{il}-H^{c}_{il\sigma}(t)\right]G^{cc}_{lj\sigma}(t,t^{\prime})-\sum_{k\in
M}\delta_{ik}\gamma_{k}[{\mathbf{S}}(t)]\,G^{ac}_{kj\sigma}(t,t^{\prime})$
$\displaystyle=\delta_{\cal C}(t,t^{\prime})\,\delta_{ij}+\sum_{l}\int_{\cal
C}\textup{d}\bar{t}\,\Sigma^{c}_{il\sigma}(t,\bar{t})\,G^{c}_{lj\sigma}(\bar{t},t^{\prime})\,,$
(24)
$\displaystyle\sum_{l}\left[\textup{i}\hbar\,\partial_{t}\delta_{il}-H^{c}_{il\sigma}(t)\right]G^{ca}_{lk^{\prime}\sigma}(t,t^{\prime})-\sum_{k\in
M}\delta_{ik}\gamma_{k}[{\mathbf{S}}(t)]\,G^{a}_{kk^{\prime}\sigma}(t,t^{\prime})$
$\displaystyle=0\,,$ (25)
$\displaystyle\left[\textup{i}\hbar\,\partial_{t}-\epsilon_{k}(t)\right]G^{ac}_{kj\sigma}(t,t^{\prime})-\sum_{k^{\prime}\in
M}\delta_{kk^{\prime}}\gamma_{k^{\prime}}[{\mathbf{S}}(t)]\,G^{c}_{k^{\prime}j\sigma}(t,t^{\prime})$
$\displaystyle=0\,,$ (26)
$\displaystyle\left[\textup{i}\hbar\,\partial_{t}-\epsilon_{k}(t)\right]G^{a}_{kk^{\prime}\sigma}(t,t^{\prime})-\sum_{k^{\prime\prime}\in
M}\delta_{kk^{\prime\prime}}\gamma_{k^{\prime\prime}}[{\mathbf{S}}(t)]\,G^{ca}_{k^{\prime\prime}k^{\prime}\sigma}(t,t^{\prime})$
$\displaystyle=\delta_{\cal C}(t,t^{\prime})\,\delta_{kk^{\prime}}\,,$ (27)

where $H^{\textup{c}}_{ij\sigma}(t)$ is the one-particle Hamiltonian of the
graphene nanoflake,

$\displaystyle
H^{\textup{c}}_{ij\sigma}(t)=\left\\{\begin{array}[]{cc}W_{i}[Z_{0};{\mathbf{S}}(t)]+E+U\left(\langle\hat{n}_{i\bar{\sigma}}\rangle-\frac{1}{2}\right)\,,&i=j\\\\[6.0pt]
W_{\langle ij\rangle}[Z_{0};{\mathbf{S}}(t)]-J\,,&ij=\langle
ij\rangle\\\\[6.0pt] 0\,,&\textup{otherwise}\end{array}\right.\,,$ (31)

$\delta_{\cal C}(t,t^{\prime})$ is the delta function on the Keldysh contour,
and $\Sigma^{c}_{ij\sigma}(t,t^{\prime})$ in Eq. (24) denotes the many-body
self-energy 25 of the GNF. As mean-field (Hartree) contributions of the system
(2) are already contained in the one-particle Hamiltonian (31), the self-
energy $\Sigma^{c}_{ij\sigma}(t,t^{\prime})$ accounts for electron-electron
correlations only. Moreover, we note that Eqs. (24)-(27) need to be
accompanied by their adjoint equations, containing the time derivative with
respect to $t^{\prime}$.

In the following, we perform a formal decoupling of the GNF part of the
Keldysh-Kadanoff-Baym equations (24–27) from the ion part 19, 17. To this end,
we assume that the NEGF of the ion can be written as
$G^{a}_{kk^{\prime}\sigma}(t,t^{\prime})=\delta_{kk^{\prime}}g_{k\sigma}(t,t^{\prime})$,
where the nonequilibrium Green function $g_{k\sigma}(t,t^{\prime})$ obeys the
ideal (noninteracting) KBE

$\displaystyle[\textup{i}\hbar\,\partial_{t}-\epsilon_{k}]g_{k\sigma}(t,t^{\prime})=\delta_{\cal
C}(t,t^{\prime})\,,$ (32)

which is solved by

$\displaystyle g_{k\sigma}(t,t^{\prime})$
$\displaystyle=\frac{\textup{i}}{\hbar}[n_{k\sigma}-\theta_{\cal
C}(t,t^{\prime})]\exp\left(-\frac{\textup{i}}{\hbar}\epsilon_{k}(t-t^{\prime})\right)\,,$
(33)

where $\theta_{\cal C}(t,t^{\prime})$ is the step function on the contour, and
$n_{k\sigma}=\langle\hat{a}_{k\sigma}^{\dagger}\hat{a}_{k\sigma}\rangle_{\textup{s}}$
denotes the occupation of the energy level $\epsilon_{k}$. Using this ansatz,
we find that the charge transfer Green function can be expressed as

$\displaystyle G^{ac}_{kj\sigma}(t,t^{\prime})=\int_{\cal
C}\textup{d}\bar{t}\,g_{k\sigma}(t,\bar{t})\,\gamma_{k}[{\mathbf{S}}(\bar{t})]\,G_{kj}^{c}(\bar{t},t^{\prime})\,,$
(34)

and insertion of (34) into Eq. (24) leads to the closed equation of motion

$\displaystyle\sum_{l}$
$\displaystyle\left[\textup{i}\hbar\,\partial_{t}\delta_{il}-H^{c}_{il\sigma}(t)\right]G^{c}_{lj\sigma}(t,t^{\prime})=\delta_{\cal
C}(t,t^{\prime})\delta_{ij}+\sum_{l}\int_{\cal
C}\textup{d}\bar{t}\,[\Sigma^{c}_{il\sigma}(t,\bar{t})+\Sigma^{\textup{emb}}_{il\sigma}(t,\bar{t})]\,G^{c}_{lj\sigma}(\bar{t},t^{\prime})\,,$
(35)

where we have defined the embedding self-energy

$\displaystyle\Sigma_{ij\sigma}^{\textup{emb}}(t,t^{\prime})=\left\\{\begin{array}[]{cc}\gamma_{i}[{\mathbf{S}}(t)]\,g_{i\sigma}(t,t^{\prime})\,\gamma_{i}[{\mathbf{S}}(t^{\prime})]\,,&i=j\in
M\\\\[6.0pt] 0\,,&\textup{otherwise}\end{array}\right.\,.$ (38)

Methodically, the additional embedding self-energy in the integration kernel
on the r.h.s. of Eq. (35) incorporates the effects of resonant charge transfer
into the GNF part (24) of the overall KBE and explicitly allows for temporal
changes of the electron number $\langle\hat{N}\rangle_{\textup{s}}(t)$ in the
graphene nanoflake. This enables us to characterize the loss of electrons from
the GNF and, hence, the degree of neutralization of the external ion.
Generally, the energy levels $\epsilon_{k}$ are unoccupied before the ion
interacts with the GNF, thus we set $n_{k\sigma}\equiv 0$ in Eq. (33).
Moreover, as each energy level $k$ can capture one electron per spin $\sigma$,
we introduce a number of ${\cal M}_{\sigma}=Z_{0}/2$ charge transfer channels
$\hat{c}_{k\sigma}\rightarrow a^{\dagger}_{k\sigma}$ in the Hamiltonian, i.e.,
${\cal M}_{\sigma}^{\textup{s}}={\cal M}_{\sigma}/4=Z_{0}/8$ transfer channels
in each of the four Hubbard bands of Eq. (2). For the case of $Z_{0}$ being
multiples of $8$, this implies the possibility of complete neutralization of
the ion during its traversal of the nanoflake. Furthermore, in the considered
scenario of normal incidence, charge transfer will occur predominantly from
the GNF center, i.e., from lattice sites that are closest to the ion and the
impact point. For this reason, we couple only sites of the nanoflake’s central
honeycomb. For the precise setup, which defines the set of indices $M$ in Eq.
(38), see Fig. 1b.

The Keldysh-Kadanoff-Baym equation (35) is the central equation which will be
solved in the next section for various parameter sets under the neglect
electron-electron correlations, i.e., for $\Sigma^{c}_{ij\sigma}\equiv 0$
(mean field approximation). However, we note that, even for
$\Sigma^{c}_{ij\sigma}\equiv 0$, the remaining embedding part of the self-
energy renders the equation of motion for $G^{c}_{ij\sigma}(t,t^{\prime})$
still an integro-differential equation on the full Keldysh contour ${\cal C}$.
Thus, the numerical effort is similar to solving a correlated problem, and the
inclusion of electron correlations, e.g., in a second-order Born or T-matrix
approximation, would be straightforward. To our advantage, the embedding self-
energy (38) is local in space
($\Sigma_{ij\sigma}^{\textup{emb}}\propto\delta_{ij}\Sigma_{i\sigma}^{\textup{emb}}$),
and, therefore, we use the auxiliary Hamiltonian approach developed in Ref. 26
to efficiently solve the KBE.

Figure 3: Charge transfer $\langle\hat{\rho}_{\textup{ct}}\rangle$, black
symbols, Eq. (39), as a function of the inverse ion velocity $v_{z}^{-1}$
obtained from the KBE (35) for a graphene nanoflake with $L=24$ sites. Panel
(a): NEGF results for an ion with initial charge state $Z_{0}=32$ and
different parameters $\gamma_{0}=1.0$, $1.5$, $2.0$ and $2.5$ in Eq. (15).
Panel (b): NEGF results for different charge states $Z_{0}=8$, $16$, $24$,
$32$, $40$, $48$ and $\gamma_{0}=2.12$. In both panels, the colored dots
indicate experimental data of Gruber et al., Ref. 9, for the transmission of
Xe${}^{Q_{\textup{in}}+}$ ions through single-layer graphene, where
$Q_{\textup{in}}=20$, $25$, $30$, $32$ and $35$, see legend in panel (a). The
colored solid lines are fits of the form (40) with the neutralization time
constant $\tau$ fitted to the experimental data, assuming a continuous
neutralization following an exponential function, compare with Ref. 9.

##  4 Results

In this section, we solve the KBE (35) with the charge transfer amplitude of
Eq. (15) and parameters $z_{\textup{res}}=-\sqrt{3}a_{0}$,
$d_{\textup{w}}=0.6a_{0}$ and energies
$\epsilon_{k}=\epsilon=-Z_{0}W_{0}/2-W_{\textup{f}}$ (setup of Fig. 1b as
described in Sec. 3) for ions with initial charge states $Z_{0}=8$, $16$,
$24$, $32$, $40$, $48$ and velocities $v_{z}$ in the range of $0.1$ and
$100v_{0}$, i.e., with velocities between $8.2\times 10^{4}$ and $8.2\times
10^{7}$ m/s. Due to the extremely large and time-dependent potential energy,
which is induced inside the graphene nanoflake by the presence of the highly
charged ion and which easily exceeds several hundred electron volts 7 (e.g.,
for $Z_{0}=32$, the maximum field strength on the central GNF lattice sites is
given by $-W_{0}Z_{0}\approx-325$ eV, all NEGF simulations are carried out on
a fine ($t,t^{\prime}$) time mesh. To resolve an extended time window of
$[-75t_{0},75t_{0}]$, where the ion passes through the GNF layer at time
$t=0$, we use $n_{\textup{t}}=12,000$ time steps, which corresponds to an
integration step size and temporal resolution of $\Delta t/t_{0}=0.0125$.
Initially, at time $t=-75t_{0}$, the graphene nanoflake is prepared in the
half-filled Hartree ground state of the paramagnetic phase (the filling is
$f=\frac{1}{L}\sum_{i\sigma}\langle\hat{n}_{i\sigma}\rangle_{\textup{s}}=0.5$
with
$\langle\hat{n}_{i\uparrow}\rangle_{\textup{s}}=\langle\hat{n}_{i\downarrow}\rangle_{\textup{s}}\,\forall\,i$).

We emphasize again, that, for the described model setup and, following the
self-energy approach of Sec. 3, the only remaining open parameter is the
amplitude $\gamma_{0}$ of the charge transfer function (15). Moreover, the
constraint of $Z_{0}$ being a multiple of $8$ is due to the restriction of the
model Hamiltonian to four equivalent Hubbard bands ($4$ bands $\times$ $2$
spin-projections) in combination with the choice of considered charge transfer
channels (cf. Sec. 3 and Fig. 1b).

###  4.1 Charge transfer for different ion velocities and charge states. Test
of the model

Our solutions of the KBE (35) reveal that the number of electrons that are
emitted from the graphene nanoflake sensitively depends on the interaction
time between the HCI and the GNF [which is set by the ion velocity $v_{z}$]
and on the ion’s charge state $Z_{0}$ [which determines the strength of the
local field inside the GNF lattice]. Figure 3 investigates, for a broad
range of parameters, the overall transfer of charge onto the approaching ion,
which is defined as

$\displaystyle\langle\hat{\rho}_{\textup{ct}}\rangle=4\left(L-\langle\hat{N}\rangle_{\textup{s}}(t\rightarrow+\infty)\right)\,,$
(39)

where $L$ denotes the number of lattice sites in the graphene nanoflake.

In Fig. 3a, we consider a GNF with $L=24$ sites, an ion with $Z_{0}=32$ and
vary the ion velocity $v_{z}$ as well as the amplitude $\gamma_{0}$. For small
ion velocities and a sufficiently large value of $\gamma_{0}$, we observe that
the charge transfer (black symbols) allows for a complete neutralization of
the ion, i.e., $\langle\rho_{\textup{ct}}\rangle\sim Z_{0}$. On the other
hand, for high ion velocities, the charge transfer rapidly decreases and tends
towards zero, which is a direct consequence of the reduced interaction time
(at constant electron mobility) and ultimately leaves the electronic system
undisturbed for $v_{z}^{-1}\rightarrow 0$. These two limiting cases are
recovered independently of the choice of $\gamma_{0}$. On the other hand, at
intermediate ion velocities, the amount of electrons transferred to the ion
and the dependence on $v_{z}$ sensitively depend on the choice of the charge
transfer amplitude $\gamma_{0}$. A key question of the present work is, how
accurately the simple model (15) describes the physical processes of the
interaction of the charged projectile and the graphene cluster. Fortunately,
there exists an extensive set of experimental data for the neutralization of
Xe${}^{Q_{\textup{in}}+}$ ions in single-layer graphene at similar impact
parameters 9. The measured results for ions with an initial charge
$Q_{\textup{in}}=32$ are shown by the orange dots in Fig. 3a and show that
the overall $v_{z}$-dependence is well captured by our model. Moreover,
adjusting a single open parameter, $\gamma_{0}$, allows us to bring the two
curves into complete agreement, for $\gamma\approx 2$. The best fit is
achieved for $\gamma_{0}=2.12$ which will be used in the following.

Figure 4: Time evolution of the overall particle number
$\langle\hat{N}\rangle_{\textup{s}}(t)$ and the site occupation
$\langle\hat{n}_{0\sigma}\rangle_{\textup{s}}(t)$ [central site of the GNF,
cf. site index A in Fig. 1a] during ion impact for a graphene nanoflake
with $L=24$ sites and (a) an initial ion charge state $Z_{0}=16$, (b)
$Z_{0}=32$ and (c) $Z_{0}=48$. The different colors refer to NEGF results for
different ion velocities that are given in the legend of the left figure. All
other parameters are as in Fig. 3b. In panel (b), the arrows indicate the
time $t_{\textup{res}}=z_{\textup{res}}/v_{z}$, when the ion passes through
the resonance point $z_{\textup{res}}=-\sqrt{3}a_{0}$, in front of the plane
of the graphene nanoflake.

The next question is whether this good agreement is just a coincidence for a
single charge state. To answer this, we display in Fig. 3b, all
experimental data available from Ref. 9, i.e., in addition to
$Q_{\textup{in}}=32$, also the charge transfer for $Q_{\textup{in}}=20$, $25$,
$30$, and $35$, cf. the colored dots and the fits (lines of the same color).
These data are compared to our simulations for different charge states that
are multiples of eight (see legend in the figure) where we now fix
$\gamma_{0}=2.12$ that was obtained from the case $Q_{\textup{in}}=32$ (left
figure). The striking result is that this single parameter value allows us to
reproduce the experimental curves for all available charge states, compare,
e.g., the blue line for $Q_{\textup{in}}=25$ with the simulations for
$Z_{0}=24$ (black circles). For other values agreement can be established by
interpolating between the available charge values. This gives strong support
for the present simple model (15). Moreover, this allows us also to predict
the charge transfer behavior for small impact velocities and additional charge
states for which no experimental data are available.

This good agreement with the experiments gives us confidence that our model,
combined with our KBE dynamics, captures the correct physics. Generally, we
find that ions with a lower initial charge state get fully neutralized already
at smaller inverse ion velocities and thus have smaller neutralization time
constants $\tau$, cf. also Ref. 9. From the figure and, using the fit formula

$\displaystyle\langle\hat{\rho}_{\textup{ct}}\rangle=Q_{\textup{in}}\left(1-\textup{e}^{-a_{0}/(v_{z}\tau)}\right),$
(40)

we furthermore extract neutralization time constants in the range of
$\sim\\!0.7\,t_{0}$ (for $Z_{0}=8$) to $\sim\\!2.6\,t_{0}$ (for $Z_{0}=48$),
with respect to the ion travel distance of one lattice spacing (carbon-carbon
bond length $a_{0}$).

###  4.2 Density response and finite size effects

After having verified the good quality of our model, in this section, we
analyze more details of the charge transfer. In particular, we investigate the
time-resolved electronic response of the graphene nanoflake during impact of
the HCI and perform simulations for GNFs of different size.

Figure 5: Electronic response to ion impact of different hexagonal hollow-
centred GNFs with $L=24$ (blue curves), $54$ (green), $96$ (orange) and $150$
(red) lattice sites. Ion parameters: Initial charge state $Z_{0}=32$ and
velocity $v_{z}$ as indicated. Upper three panels: time-dependent occupations
$\langle\hat{n}_{i\sigma}\rangle_{\textup{s}}(t)$ for sites
$i=\textup{A},\textup{B},\textup{C}$ as defined in Fig. 1a. Bottom panel:
time-resolved charge transfer $\langle\hat{\rho}_{\textup{ct}}\rangle(t)$, Eq.
(39), as well as the charge transfer function $\gamma_{k}[{\mathbf{S}}(t)]$
(violet peak). Furthermore, the black line in the third row indicates the
time-dependent potential energy $E_{k}[{\mathbf{S}}(t)]$ on site A and all
other sites of the central honeycomb (right ordinate).

In Fig. 4a-c, we monitor the total particle number in a single-band,
$\langle\hat{N}\rangle_{\textup{s}}(t)$, (lower panels) and the mean
occupation of the central site,
$\langle\hat{n}_{\textup{A}\sigma}\rangle_{\textup{s}}(t)$, (upper panels)
around the time of impact ($t=0$) for a GNF with $L=24$ sites and three
different values of $Z_{0}$ at various ion velocities. Generally, the
electrons are attracted by the approaching ion and thus move inside the
graphene nanoflake, i.e., along the edges of the honeycomb lattice, towards
the impact point and accumulate in the GNF center (cf. also Ref. 9). This
explains the increase of
$\langle\hat{n}_{\textup{A}\sigma}\rangle_{\textup{s}}$, for $t<0$, before the
onset of electron emission given by the decline of
$\langle\hat{N}\rangle_{\textup{s}}$. The accumulation of electrons in the GNF
center is thereby more pronounced for large ion velocities, as in this case
less electrons are simultaneously emitted from the central honeycomb.
Particularly for $v_{z}^{-1}<1.0$, this even leads to a clearly enhanced
double occupation
$\langle\hat{n}_{\textup{A}\uparrow}\hat{n}_{\textup{A}\downarrow}\rangle\approx\langle\hat{n}_{\textup{A}\uparrow}\rangle\langle\hat{n}_{\textup{A}\downarrow}\rangle$
(Hartree approximation) in the graphene nanoflake and, thus, potentially to
the formation of doublons 6. Moreover, the emission of electrons leads to the
emergence of fast Bloch oscillations of the central density, the amplitude
(frequency) of which decreases (increases) with the field strength generated
by the ion. This is seen in the upper panels of Fig. 4 a to c where the ion
charge is varied form 16 to 48.

For $t>0$, the ion quickly departs from the nanoflake. Since in the present
model no further significant charge transfer occurs, the total particle number
remains constant. At the same time, the density on the central site finally
starts to relax in an oscillatory manner where details depend on the remaining
filling of the GNF. Also, the period of the oscillation increases with
increasing distance of the ion due to the decay of its Coulomb field. The
equilibration process is considered in some more detail in Fig. 5, which
shows the time evolution of the electron density along a path
$\overline{\textup{A}\textup{B}\textup{C}}$ of adjacent sites (see Fig. 1a)
away from the GNF center and discusses the influence of the cluster size $L$.
On the innermost site A, from which electron emission occurs, we observe that,
independently of the ion velocity, the time-dependent occupation around the
time of ion impact converges for $L\gtrsim 96$, whereas for times essentially
before and after the impact as well as on sites B and C the dynamics still
depends on the cluster size. On sites B and C, we clearly see the accumulation
of electrons during the impact of the ion, particularly for cluster sizes
$L\geq 54$, where the maximum site occupation increases with the cluster size
and approaches values
$\langle\hat{n}_{\textup{C}\sigma}\rangle_{\textup{s}}\gtrsim 0.9$, for the
largest nanoflake, cf. in particular the top panel in Fig. 5b, where the
site occupation remains large for several femtoseconds. For $L=24$, we note
that site C is located at the edge of the nanoflake, and electron transfer
towards the center is not compensated from more distant sites. Therefore, site
C exhibits a density decrease, in striking contrast to the other sites.
Moreover, due to the electron capture and the finite nature of the graphene
nanoflakes, the ion always leaves behind nanoclusters with a filling $f<0.5$.
This deviation from the half-filled case is of course most pronounced for the
smallest GNF containing the smallest amount of electrons, see the blue curves
which, for $t\rightarrow\infty$, indicate final occupations considerably less
than $0.5$.

Finally, the bottom panels of Fig. 5 show that the overall charge transfer
does in the present model not depend significantly on the size of the graphene
nanocluster (all curves are virtually on top of each other). We note, that
this is the case even though the occupations on site A clearly differ for the
different nanoflakes directly before the onset of the charge transfer.

##  5 Discussion and conclusions

In summary, we have studied the time-dependent neutralization dynamics of
highly charged ions which resonantly capture electrons during transmission
through a two-dimensional atomic monolayer. This was done by extending our
recently developed Nonequilibrium Green functions-Ehrenfest dynamics approach
5, 6 to include charge transfer processes. As test systems we have chosen
finite, hexagonal graphene nanoflakes. Our primary focus has been on the
electronic response of the GNF layer, while the treatment of the ion was
confined to the occupation of excited states with an effective energy
$\epsilon_{k}=\epsilon$ that correspond to a transient configuration of a
“hollow atom (ion)”. From the information about the occupied level
$\epsilon_{k}$ we reconstructed the one-particle NEGF of the ion [Eq. (23)]
which allowed us to compute the two-time embedding self-energy
$\Sigma^{\textup{emb}}_{ij\sigma}(t,t^{\prime})$ [Eq. (38)]. Due to the
appearance of the embedding self-energy in the integration kernel of the
Keldysh-Kandaoff-Baym equation (35), the number of electrons in the graphene
nanoflake is no longer conserved, even for a conserving approximation of the
regular many-body self-energy $\Sigma^{c}_{ij\sigma}(t,t^{\prime})$, and thus,
by construction, the change of particle number corresponds to the number of
electrons emitted from the GNF. This procedure resembles the embedding schemes
that were used before in NEGF calculations of transient photoabsorption 27, 28
or of the transient dynamics in quantum transport setups 29. The NEGF
simulations were combined with an analytical parametrization of the charge
transfer amplitude. Comparison to available experimental data for highly
charged ions revealed excellent agreement giving strong support for our
approach.

Let us discuss limitations and possible future improvements. The first
limitation of our approach is that the initially formed hollow ion state does
not undergo further de-excitation as a function of time. These processes might
influence the neutralization as well as the charge transfer. For fast impact
velocities, and hence small GNF-ion interaction times, such de-excitation
effects should, however, not essentially affect the emission of electrons from
the graphene nanoflake. Secondly, we have assumed that the potential energy
induced inside the GNF [Eqs. (6) and (7)] only depends on the HCI’s initial
charge $Z_{0}$ during the whole interaction process, although the capture of
electrons modifies the charge state as function of time. This effect could be
incorporated into the model via an effective, time-dependent charge state
$Z(t)$, which is constructed from the number of electrons emitted from the
nanoflake. For the equation of motion of the nanoflake’s NEGF, such a feedback
loop would, however, necessitate an iterative solution of the problem. It
remains to be checked whether the KBE will converge in such a procedure. On
the other hand, we are confident that the current ansatz is again justified
for sufficiently fast ions, where electrons move into energetically high-lying
and loosely bound hollow state wave packets, with which the definition of a
clear ion charge state is delicate anyway. Furthermore, the excellent
agreement of our NEGF results with the experimental data for a reasonably
motivated choice of the charge transfer amplitude [Sec. 2.2] gives us
confidence that the approach includes most of the relevant physics. A thorough
comparison of the results to predictions obtained from a classical over-the
barrier model 30 is left for future work. The same concerns the extension of
the simulations to include charge transfer from lattices sites outside the
innermost honeycomb ring. Moreover, an interesting topic is to investigate the
influence of the charge transfer on the stopping power and the modification of
the charge transfer by the initial acceleration of the ion 5.

The main advantage of the embedding self-energy approach is that the numerical
effort is not substantially larger compared to a NEGF simulation for the
closed, isolated graphene nanoflake system and that it requires only minor
adjustments of existing NEGF codes. In addition, the calculations can
straightforwardly be extended beyond the mean-field level to incorporate
electron-electron correlations, e.g., using the (local) second-order Born
approximation, as in Ref. 5. More advanced self-energies (including GW or
T-matrix), larger systems and more complex and realistic Hamiltonians should
become accessible in conjunction with the current development of more powerful
and efficient KBE solvers 31, 32 based on the generalized Kadanoff-Baym
ansatz. Finally, despite the above mentioned limitations, the approach offers
further interesting applications such as the investigation of charge transfer
for electron- or hole-doped graphene nanoflakes or the study of correlation
effects and doublon formation for systems with substantially larger Coulomb
interaction strengths and, hence, pronounced band gaps.

## Acknowledgments

We acknowledge Lasse Wulff’s contributions to the code development at the
initial stage. Further, we acknowledge computing time at the Computing Centre
of Kiel University.

## References

  * 1 Michael Bonitz, Alexey Filinov, Jan-Willem Abraham, Karsten Balzer, Hanno Kählert, Eckhard Pehlke, Franz X. Bronold, Matthias Pamperin, Markus Becker, Dettlef Loffhagen, Holger Fehske, Frontiers of Chemical Science and Engineering 2019, 13 (2), 201–237.
  * 2 J. Marbach, F. X. Bronold, H. Fehske, The European Physical Journal D 2012, 66 (4), 106.
  * 3 M. Pamperin, F. X. Bronold, H. Fehske, Phys. Rev. B 2015, 91, 035440.
  * 4 Shijun Zhao, Wei Kang, Jianming Xue, Xitong Zhang, Ping Zhang, Journal of Physics: Condensed Matter 2014, 27 (2), 025401.
  * 5 Karsten Balzer, Niclas Schlünzen, Michael Bonitz, Phys. Rev. B 2016, 94, 245118.
  * 6 Karsten Balzer, Maximilian Rodriguez Rasmussen, Niclas Schlünzen, Jan-Philip Joost, Michael Bonitz, Phys. Rev. Lett. 2018, 121, 267602\.
  * 7 Friedrich Aumayr, Hannspeter Winter, e-Journal of Surface Science and Nanotechnology 2003, 1, 171–174.
  * 8 S Facsko, R Heller, A S El-Said, W Meissl, F Aumayr, Journal of Physics: Condensed Matter 2009, 21 (22), 224012.
  * 9 Elisabeth Gruber, Richard A. Wilhelm, Rémi Pétuya, Valerie Smejkal, Roland Kozubek, Anke Hierzenberger, Bernhard C. Bayer, Iñigo Aldazabal, Andrey K. Kazansky, Florian Libisch, Arkady V. Krasheninnikov, Marika Schleberger, Stefan Facsko, Andrei G. Borisov, Andrés Arnau, Friedrich Aumayr, Nature Communications 2016, 7 (1), 13948.
  * 10 U Kentsch, S Landgraf, M Schmidt, H Tyrroff, G Zschornack, F Grossmann, V.P Ovsyannikov, F Ullmann, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 2004, 216, 196 – 201, Proceedings of the E-MRS 2003 Symposium E on Ion Beams for Nanoscale Surface Modifications.
  * 11 F. Aumayr, A.S. El-Said, W. Meissl, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 2008, 266 (12), 2729 – 2735, Radiation Effects in Insulators.
  * 12 Robert Ritter, Richard A. Wilhelm, Michael Stöger-Pollach, Rene Heller, Arndt Mücklich, Udo Werner, Henning Vieker, Andre Beyer, Stefan Facsko, Armin Gölzhäuser, Friedrich Aumayr, Applied Physics Letters 2013, 102 (6), 063112.
  * 13 T. Schenkel, A. V. Hamza, A. V. Barnes, M. W. Newman, G. Machicoane, T. Niedermayer, M. Hattass, J. W. McDonald, D. H. Schneider, K. J. Wu, R. W. Odom, Physica Scripta 1999, T80 (A), 73.
  * 14 Walter Meissl, Daniel Winklehner, Friedrich Aumayr, Martin C. Simon, Rainer Ginzel, José; R. Crespo López-Urrutia, Joachim Ullrich, Beate Solleder, Christoph Lemell, Joachim Burgdörfer, e-Journal of Surface Science and Nanotechnology 2008, 6, 54–59.
  * 15 Masahide Tona, Satoshi Takahashi, Kazuo Nagata, Nobuo Yoshiyasu, Chikashi Yamada, Nobuyuki Nakamura, Shunsuke Ohtani, Makoto Sakurai, Applied Physics Letters 2005, 87 (22), 224102.
  * 16 N. Schlünzen, K. Balzer, M. Bonitz, L. Deuchler, E. Pehlke, Contributions to Plasma Physics 2019, 59 (6), e201800184.
  * 17 Michael Bonitz, Karsten Balzer, Niclas Schlünzen, Maximilian Rodriguez Rasmussen, Jan-Philip Joost, physica status solidi (b) 2019, 256 (7), 1800490.
  * 18 S. Creutzburg, J. Schwestka, A. Niggas, H. Inani, M. Tripathi, A. George, R. Heller, R. Kozubek, L. Madauß, N. McEvoy, S. Facsko, J. Kotakoski, M. Schleberger, A. Turchanin, P. L. Grande, F. Aumayr, R. A. Wilhelm, Phys. Rev. B 2020, 102, 045408.
  * 19 G. Stefanucci, R. van Leeuwen, Nonequilibrium Many-Body Theory of Quantum Systems, Cambridge: Cambridge University Press, 2013.
  * 20 K. Balzer, M. Bonitz, Nonequilibrium Green’s Functions Approach to Inhomogeneous Systems, Springer, Berlin Heidelberg, 2013.
  * 21 M. Schüler, M. Rösner, T. O. Wehling, A. I. Lichtenstein, M. I. Katsnelson, Phys. Rev. Lett. 2013, 111, 036601.
  * 22 Young-Jun Yu, Yue Zhao, Sunmin Ryu, Louis E. Brus, Kwang S. Kim, Philip Kim, Nano Lett. 2009, 9, 3430.
  * 23 L. Kadanoff, G. Baym, Quantum Statistical Mechanics, New York: Benjamin, 1962.
  * 24 L.V. Keldysh, Soviet Phys. JETP 1965, 20, 1018, (Zh. Eksp. Teor. Fiz. 47, 1515 (1964)).
  * 25 N Schlünzen, S Hermanns, M Scharnke, M Bonitz, Journal of Physics: Condensed Matter 2019, 32 (10), 103001.
  * 26 Karsten Balzer, Martin Eckstein, Phys. Rev. B 2014, 89, 035148\.
  * 27 E. Perfetto, A.-M. Uimonen, R. van Leeuwen, G. Stefanucci, Phys. Rev. A 2015, 92, 033419.
  * 28 E Perfetto, G Stefanucci, Journal of Physics: Condensed Matter 2018, 30 (46), 465901.
  * 29 Petri Myöhänen, Adrian Stan, Gianluca Stefanucci, Robert van Leeuwen, Phys. Rev. B 2009, 80, 115107.
  * 30 Joachim Burgdörfer, Peter Lerner, Fred W. Meyer, Phys. Rev. A 1991, 44, 5674–5685.
  * 31 Niclas Schlünzen, Jan-Philip Joost, Michael Bonitz, Phys. Rev. Lett. 2020, 124, 076601.
  * 32 Jan-Philip Joost, Niclas Schlünzen, Michael Bonitz, Phys. Rev. B 2020, 101, 245101.