# Interference and reflection from the event horizon of a quantum corrected
black hole

Sunandan Gangopadhyaya∗, Soham Sena† and Rituparna Mandala‡ a Department of
Astrophysics and High Energy Physics
S.N. Bose National Centre for Basic Sciences
JD Block, Sector III, Salt Lake, Kolkata 700 106, India

###### Abstract

In this work, we calculate the Hawking temperature for a quantum corrected
black hole geometry using the reflection from the horizon method. We observe
that quantum gravity corrections indeed show up in the Hawking temperature
formula of the quantum corrected black hole. It is important to notice that
the quantum gravity corrections arise in the Hawking temperature formula only
due to the underlying quantum gravity corrections to the lapse function of the
black hole metric rather than the semi classical methods used in the analysis.
We also substantiate our result by computing the Hawking temperature using the
tunneling approach.

## 1 Introduction

111
<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>
General theory of relativity developed by Albert Einstein is considered as the
most accurate theory describing the large scale structure of the universe[1,
2]. The Einstein’s field equations admit solutions with singularities which
are called black holes. In classical consideration any object that falls into
the event horizon of the black hole can never escape from it. Therefore
classically black holes act as perfect absorbers. This classical picture has
known limitations. From the viewpoint of thermodynamics, a black hole should
have entropy and a temperature. It was shown in [3] that if one tries to
extract energy from a Kerr black hole, there is a quantity which never
decreases and it was later found to be proportional to the area of the black
hole. Future investigations revealed that the area of the black hole is its
physical entropy[4, 5, 6]. Some ground breaking works also revealed that the
black holes have their own laws of thermodynamics[7]. It was first shown by
Stephen Hawking that if one considers quantum fluctuations in the curved
background, the black hole emits radiation[8, 9, 10]. Later it was implemented
by Stephen Hawking that the black hole radiation is identical to that of a
black body radiation with temperature $T=\frac{\kappa}{2\pi}$, where $\kappa$
is the surface gravity of the black hole. There have been several attempts to
give an alternative derivation of the Hawking radiation[11, 12, 13]. Some of
the very famous approaches to understand the origin of the Hawking temperature
revolves around the reflection from the horizon method[14] and the well known
tunneling approach[15, 16]. The radial null geodesic approach has been used
later to compute the Hawking temperature for different black hole
geometries[17, 18, 21, 20, 19]. The tunneling of Dirac particle through the
event horizon has also been considered in [22, 23, 24, 25, 26, 27].

It was shown in [14] that if a particle is falling into the event horizon of a
black hole, there is a finite probability that the particle will be reflected
back from the event horizon of the black hole. If one considers particle
trajectories around the black hole, classically there are ingoing trajectories
and the outgoing trajectories. The outgoing trajectories describe the
trajectories that leads a particle away from the centre of the black hole and
the ingoing trajectories on the other hand leads particles towards the
singularity of the black hole. If a particle enters the event horizon of a
black hole, there are no ways to escape classically. But if one consider
quantum mechanical analysis in the vicinity of the event horizon, it can be
observed that the solution of the incoming particle consists of both the
ingoing and outgoing trajectories to allow superposition between the two
trajectories. This outgoing part of the solution implies a finite probability
of a particle to be reflected back from the event horizon of the black hole
making the absorption cross-section finite which is proportional to the event
horizon area in the infrared region[28, 29]. Now the Hawking radiation process
includes radiations being emitted from the event horizon of the black hole. It
was analyzed in [14] that if one calculates the temperature of the black hole
for simple black hole geometries, it coincides with the Hawking temperature of
the black hole implying that this reflection from the horizon can work as an
alternative origin story for the Hawking radiation of a general black hole
geometry.

The tunneling approach serves as an alternate description of the source of the
Hawking radiation which describes the source of this radiation by tunneling of
scalar waves across the event horizon. The core idea of this approach lies in
the fact that there is an abrupt change in sign of the energy of the particle
when it crosses the event horizon of a black hole. If a pair is created
outside or inside of the event horizon of the black hole, after one member of
the pair has tunneled to the opposite site, it can materialize with zero total
energy. The mass of the black hole goes down as it radiates to maintain
conservation of energy. These black holes are therefore in highly excited
states from a quantum gravity point of view. The more popular approach to
treat Hawking radiation is by considering the black hole immersed in a thermal
bath in which equilibrium is possible. Hence, by the principle of detailed
balance there must be emission from the black hole itself.

In this work, our main motivation is to investigate the Hawking temperature of
quantum gravity corrected black hole geometries via both the reflection from
the horizon method and the tunneling approach. As we are working with quantum
gravitational black holes, one expects to find quantum gravitational
corrections to the classical form of the Hawking temperature. Here we have
considered the quantum corrected black hole geometry obtained from the flow of
the Newton’s gravitational constant and the Garfinkle-Horowitz-Strominger
(GHS) black hole.

In section 2 we have obtained the scalar field solution of the covariant
Klein-Gordon equation using the Hamilton-Jacobi approach. In sections 3 and 4,
we have obtained the Hawking temperatures for the quantum corrected black hole
and the GHS black hole using the the method of reflection from the horizon.

## 2 Hamilton-Jacobi method for obtaining the scalar field solution

In the $s$-wave approximation, the spacetime structure of a spherically
symmetric black hole is effectively $1+1$ \- dimensional. The generic metric
structure for such a static spherically symmetric black hole reads

$ds^{2}=-f(r)dt^{2}+g(r)^{-1}dr^{2}~{}.$ (1)

We shall now consider a scalar field in the above background. It satisfies the
Klein-Gordon quantum field equation for a scalar field with rest mass $m$ in
two spacetime dimensions

$-\frac{\hbar^{2}}{\sqrt{-g}}\partial_{\mu}\left(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\right)\Psi=m^{2}\Psi$
(2)

where $\sqrt{-g}=\sqrt{-\det{(g_{\mu\nu})}}=\sqrt{\frac{f(r)}{g(r)}}$ and
$\Psi$ is the scalar field. We will be considering massless scalar fields in
our analysis, and using the metric structure in eq.(1), we can express eq.(2)
as follows

$-\frac{\partial_{t}^{2}\Psi}{f(r)}+\left(\frac{f^{\prime}(r)}{2f(r)}+\frac{g^{\prime}(r)}{2g(r)}\right)g(r)\partial_{r}\Psi+g(r)\partial_{r}^{2}\Psi=0~{}.$
(3)

In order to obtain the solution of eq.(3), we take an ansatz of the form

$\Psi(t,r)=\exp\left(-\frac{i}{\hbar}I(t,r)\right)~{}.$ (4)

Substituting the ansatz in eq.(4) for $\Psi(t,r)$ in eq.(3), we obtain the
following equation involving $I(t,r)$

$\begin{split}&\frac{i}{f(r)}\left(\frac{\partial I}{\partial
t}\right)^{2}-\frac{\hbar}{f(r)}\frac{\partial^{2}I}{\partial
t^{2}}-ig(r)\left(\frac{\partial I}{\partial r}\right)^{2}\\\ &+\hbar
g(r)\frac{\partial^{2}I}{\partial r^{2}}+\hbar
g(r)\left(\frac{f^{\prime}(r)}{2f(r)}+\frac{g^{\prime}(r)}{2g(r)}\right)\frac{\partial
I}{\partial r}=0~{}.\end{split}$ (5)

It can be inferred from the forms of eq.(s)(3,5) that we can separate the
equations involving the time coordinate ($t$) and the radial coordinate ($r$).
Therefore, we can consider the form of $I(r,t)$ to be

$I(t,r)=\varepsilon t+\tilde{S}(r)$ (6)

where $\varepsilon$ denotes the energy of the particle.

Now substituting the form of $I(r,t)$ from eq.(6) in eq.(5), we get

$\begin{split}&\frac{i\varepsilon^{2}}{f(r)}-ig(r)\left(\frac{\partial\tilde{S}}{\partial
r}\right)^{2}+\hbar g(r)\frac{\partial^{2}\tilde{S}}{\partial r^{2}}\\\
&+\hbar
g(r)\left(\frac{f^{\prime}(r)}{2f(r)}+\frac{g^{\prime}(r)}{2g(r)}\right)\frac{\partial\tilde{S}}{\partial
r}=0~{}.\end{split}$ (7)

We now expand $\tilde{S}(r)$ in a power series in $\hbar$ as

$\tilde{S}(r)=S_{0}(r)+\hbar S_{1}(r)+\hbar^{2}S_{2}(r)+\cdots.$ (8)

Using this form of $\tilde{S}(r)$ from eq.(8), we can recast eq.(7) in the
following form

$\begin{split}&\frac{i\varepsilon^{2}}{f(r)}-ig(r)\biggr{(}\frac{\partial
S_{0}}{\partial r}\biggr{)}^{2}+\hbar g(r)\biggr{[}-2i\frac{\partial
S_{0}}{\partial r}\frac{\partial S_{1}}{\partial
r}+\frac{\partial^{2}S_{0}}{\partial r^{2}}\\\
&+\left[\frac{f^{\prime}(r)}{2f(r)}+\frac{g^{\prime}(r)}{2g(r)}\right]\frac{\partial
S_{0}}{\partial r}\biggr{]}+\hbar^{2}g(r)\biggr{[}-2i\frac{\partial
S_{0}}{\partial r}\frac{\partial S_{2}}{\partial r}\\\ &-i\left(\frac{\partial
S_{1}}{\partial r}\right)^{2}+\frac{\partial^{2}S_{1}}{\partial
r^{2}}+\left[\frac{f^{\prime}(r)}{2f(r)}+\frac{g^{\prime}(r)}{2g(r)}\right]\frac{\partial
S_{1}}{\partial r}\biggr{]}+\cdots=0~{}.\end{split}$ (9)

Equating the terms involving equal powers of the reduced Planck’s constant
($\hbar$) from eq.(9) to zero, we obtain the following set of equations

$\displaystyle\hbar^{0}:$
$\displaystyle~{}\frac{i\varepsilon^{2}}{f(r)}-ig(r)\biggr{(}\frac{\partial
S_{0}}{\partial r}\biggr{)}^{2}=0~{},$ (10) $\displaystyle\hbar^{1}:$
$\displaystyle~{}-2i\frac{\partial S_{0}}{\partial r}\frac{\partial
S_{1}}{\partial r}+\frac{\partial^{2}S_{0}}{\partial
r^{2}}+\left[\frac{f^{\prime}(r)}{2f(r)}+\frac{g^{\prime}(r)}{2g(r)}\right]\frac{\partial
S_{0}}{\partial r}=0~{},$ (11) $\displaystyle\hbar^{2}:$
$\displaystyle~{}-2i\frac{\partial S_{0}}{\partial r}\frac{\partial
S_{2}}{\partial r}-i\left(\frac{\partial S_{1}}{\partial
r}\right)^{2}+\left[\frac{f^{\prime}(r)}{2f(r)}+\frac{g^{\prime}(r)}{2g(r)}\right]\frac{\partial
S_{1}}{\partial r}\biggr{]}$ (12)
$\displaystyle+\frac{\partial^{2}S_{1}}{\partial r^{2}}=0~{}.$
$\displaystyle\vdots~{}~{}~{}$

In order to obtain the complete form of $\tilde{S}(r)$ in eq.(8), we need to
solve the complete set of equations given by eq.(s)(10-12) and the higher
order equations in $\hbar$ as well. From eq.(10), we obtain the following
solution

$\frac{\partial S_{0}}{\partial r}=\pm\frac{\varepsilon}{\sqrt{f(r)g(r)}}~{}.$
(13)

Integrating the above equation, with respect to the radial coordinate, we
obtain

$S_{0}=\pm\varepsilon\int^{r}\frac{dr}{\sqrt{f(r)g(r)}}~{}.$ (14)

Substituting the form of $\frac{\partial S_{0}}{\partial r}$ from eq.(13) in
eq.(11), we obtain the following relation

$\frac{\partial S_{1}}{\partial r}=0~{}.$ (15)

Eq.(15) implies that $S_{1}$=constant. Substituting eq.(15) back in eq.(12),
we obtain

$\frac{\partial S_{2}}{\partial r}=0~{}.$ (16)

Following a similar procedure, we find that

$\frac{\partial S_{n}}{\partial
r}=0~{},~{}\forall~{}n\in\\{1,2,3,\ldots\\}~{}.$ (17)

Eq.(17) indicates that $S_{n}~{}$= constant $\forall~{}n\in\\{1,2,3,\ldots\\}$
and therefore only $S_{0}$ survives while calculating the energy of the
system. Using eq.(14) in eq.(6), we obtain

$I(t,r)=\varepsilon t\pm\int\frac{\varepsilon}{\sqrt{f(r)g(r)}}~{}.$ (18)

Using the form of $I(t,r)$ from eq.(18), the solution for the massless scalar
field takes the form

$\Psi=\exp\left(-\frac{i\varepsilon
t}{\hbar}\mp\frac{i\varepsilon}{\hbar}\int\frac{dr}{\sqrt{f(r)g(r)}}\right)~{}.$
(19)

In the next section, we will be considering reflection of the massless scalar
field from the event horizon of a quantum corrected black hole.

## 3 Reflection from the horizon of a quantum corrected black hole

The metric of a quantum corrected black hole in $(3+1)$-dimensions (obtained
from the flow of the Newton’s gravitational constant) reads[30, 31, 32, 33,
34]

$ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta
d\phi^{2}$ (20)

where

$f(r)=1-\frac{2G(r)M}{r}$ (21)

and the form of $G(r)$ is given by (in natural units)

$G(r)=\frac{G}{1+\frac{\tilde{\omega}G}{r^{2}}}~{}.$ (22)

In eq.(22), $\tilde{\omega}$ is a constant denoting the quantum gravity
corrections to the black geometry . From the lapse function in eq.(21), we
obtain the inner and external radius of this quantum corrected black hole
geometry as follows

$r_{\pm}=GM\pm\sqrt{G^{2}M^{2}-\tilde{\omega}G}~{}.$ (23)

In case of the of the black holes with $f(r)=g(r)$, the scalar field solution
in eq.(19) reduce to the following form

$\Psi=\exp\left(-\frac{i\varepsilon
t}{\hbar}\mp\frac{i\varepsilon}{\hbar}\int\frac{dr}{f(r)}\right)~{}.$ (24)

In eq.(24), $\int^{r}\frac{dr}{f(r)}=r_{*}$ denotes the tortoise coordinate.
For the metric structure in eq.(21), we obtain the form of the tortoise
coordinate in the vicinity of the event horizon radius as

$\begin{split}\int_{r_{+}}^{r}\frac{dr}{f(r)}&\cong
GM\ln\left[r^{2}-2GMr+\tilde{\omega}G\right]+\\\
&\frac{G^{2}M^{2}}{\sqrt{G^{2}M^{2}-\tilde{\omega}G}}\ln\left[\frac{r-GM-\sqrt{G^{2}M^{2}-\tilde{\omega}G}}{r-GM+\sqrt{G^{2}M^{2}-\tilde{\omega}G}}\right]~{}.\end{split}$
(25)

Using the forms of $r_{\pm}$ from eq.(23), we can recast the above equation as

$\int_{r_{+}}^{r}\frac{dr}{f(r)}\cong\frac{2GM}{r_{+}-r_{-}}\left(r_{+}\ln\left[r-r_{+}\right]-r_{-}\ln\left[r-r_{-}\right]\right)~{}.$
(26)

The above form has two parts, one part indicates a solution inside the
external radius $r_{+}$ and the other outside it. Here we will be considering
the $\ln(r-r_{+})$ part of the solution describing the scalar field solution
outside the horizon radius $r_{+}$ as this is the only piece relevant for
obtaining the Hawking temperature. Hence, the complete solution of the
massless scalar field in the vicinity of the horizon radius $r_{+}$ is given
by (ignoring the second term in eq.(26))222A complete analysis keeping both
the terms in eq.(26) is given in an Appendix.

$\Psi(t,r)=e^{-\frac{i\varepsilon t}{\hbar}}\phi(r)$ (27)

where $\phi(r)$ has the form given by

$\phi(r)=\exp\left[{\mp\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\ln(r-r_{+})}\right]~{}.$ (28)

The ‘minus’ sign in the wavefunction in eq.(28) represents a radially
infalling massless scalar field that completely gets absorbed by the event
horizon of the black hole at $r=r_{+}$. In our case, we will be considering a
radially outgoing solution as well in order to allow interference between the
infalling and outgoing parts of the wavefunction. Therefore, the wavefunction
in the vicinity of the horizon radius $r_{+}$ can be written as

$\Phi(r)=e^{-\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}\ln(r-r_{+})}{r_{+}-r_{-}}\right)}+\mathcal{R}e^{\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}\ln(r-r_{+})}{r_{+}-r_{-}}\right)}.$ (29)

The first term in eq.(29) denotes an incoming wave whereas the outgoing wave
is described by the second term. In case of interference between the incoming
and the outgoing waves, the reflection coefficient $\mathcal{R}$ must have
non-zero value. If the black horizon acts as a perfect absorber then the
reflection coefficient $\mathcal{R}$ must be equal to zero. From the unitarity
condition, we know that $|\mathcal{R}|\leq 1$. The wave function in eq.(29) is
singular at the point $r=r_{+}$. In order to tackle this issue, we shall take
recourse to analytic continuation by considering $r-r_{+}=\mathcal{b}$ as a
complex number. Then $\mathcal{b}$ can be expressed as

$\mathcal{b}=r-r_{+}=\mathcal{a}e^{i\varphi}$ (30)

where $|\mathcal{b}|=\mathcal{a}$ is a real number and $\varphi$ denotes the
argument of the complex number $\mathcal{b}$. Using eq.(30), we can recast
eq.(29) as

$\Phi(r)=e^{-\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\ln\mathcal{b}}+\mathcal{R}e^{\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\ln\mathcal{b}}.$ (31)

In order to proceed further we will now use the analytical continuation method
[14]. In our analysis we are considering the massless scalar field very near
the outer horizon $r_{+}$ of the quantum corrected black hole. Therefore, we
can consider, $0<|\mathcal{b}|\ll 1$. We will now continue to rotate the
complex number $\mathcal{b}$ in the complex plane by an angle $2\pi$
(clockwise rotation). Under this clockwise rotation, the complex number
$\mathcal{b}$ is given by the following relation

$\mathcal{b}_{2\pi}=\mathcal{\alpha}e^{i(\varphi-2\pi)}~{}.$ (32)

Using modified complex number from the above equation, the modified wave
function $\Phi_{2\pi}(r)$ is given by

$\begin{split}\Phi_{2\pi}(r)=&e^{-\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\ln\mathcal{b}_{2\pi}}+\mathcal{R}e^{\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\ln\mathcal{b}_{2\pi}}\\\
=&\xi e^{-\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\ln\mathcal{b}}+\frac{\mathcal{R}}{\xi}e^{\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\ln\mathcal{b}}\end{split}$
(33)

where $\xi$ is given by

$\xi=\exp\left[-\frac{4\pi\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\right]~{}.$ (34)

The real differential equation satisfied by $\Phi(r)$ is also satisfied by the
analytically continued function $\Phi_{2\pi}(r)$. This condition implies that
one of the coefficients must have an absolute value equal to unity. In case of
$\Phi(r)$ in eq.(31), the first coefficient has an absolute value one. Now for
the rotated wave function $\Phi_{2\pi}(r)$ the first coefficient
$\xi=\exp\left[-\frac{4\pi\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\right]<1$, therefore the
second coefficient must have an absolute value equal to one. Hence, we can
write

$\begin{split}\frac{|\mathcal{R}|}{\xi}&=1\\\
\implies|\mathcal{R}|&=\xi=\exp\left[-\frac{4\pi\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\right]~{}.\end{split}$ (35)

Therefore, we observe that the value of the reflection coefficient
$\mathcal{R}$ is non-zero. Hence, the probability of reflection from the even
horizon is given by

$\mathcal{P}=|\mathcal{R}|^{2}=\exp\left[-\frac{8\pi\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\right]~{}.$ (36)

Now using the principle of detailed balance we can write

$\mathcal{P}=\exp\left[-\frac{\varepsilon}{k_{B}T}\right]=\exp\left[-\frac{8\pi\varepsilon
GM}{\hbar c^{3}}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\right]$ (37)

where $k_{B}$ denotes the Boltzmann constant and $T$ denotes the Hawking
temperature of the black hole. From eq.(37), we can obtain the Hawking
temperature for this quantum corrected black hole to be

$\begin{split}T&=\frac{\hbar c^{3}}{8\pi
k_{B}GM}\left(1-\frac{r_{-}}{r_{+}}\right)\\\ &=\frac{\hbar c^{3}}{8\pi
k_{B}GM}-\frac{\hbar c^{3}}{8\pi
k_{B}GM}\left[\frac{1-\sqrt{1-\frac{\hbar\tilde{\omega}c}{GM^{2}}}}{1+\sqrt{1-\frac{\hbar\tilde{\omega}c}{GM^{2}}}}\right]~{}.\end{split}$
(38)

In general, the quantum gravity correction term $\tilde{\omega}$ is very
small, hence we can simplify the above result as

$T\cong\frac{\hbar c^{3}}{8\pi
k_{B}GM}\left(1-\frac{\hbar\tilde{\omega}c}{4GM^{2}}-\frac{\hbar^{2}\tilde{\omega}^{2}c^{2}}{8G^{2}M^{4}}+\mathcal{O}(\tilde{\omega}^{3})\right)~{}.$
(39)

From the form of eq.(39), we can observe that the Hawking temperature picks up
quantum gravity corrections only due to the underlying quantum gravitational
nature of the background spacetime geometry. Note that we have now considered
clockwise rotation only. In case of counter-clockwise rotation, we would have
arrived at the following relation

$|\mathcal{R}|=\exp\left[\frac{4\pi\varepsilon GM}{\hbar
c^{3}}\left(\frac{r_{+}}{r_{+}-r_{-}}\right)\right]~{}.$ (40)

Now from eq.(40), we observe that $|\mathcal{R}|>1$ which is in contradiction
with the unitarity condition ($|\mathcal{R}|\leq 1$). Hence, the counter-
clockwise rotation is forbidden in the current analysis.

In the next section we will calculate the Hawking temperature for a Garfinkle-
Horowitz-Strominger (GHS) black hole for which $f(r)\neq g(r)$.

## 4 Reflection from the horizon of a Garfinkle -Horowitz-Strominger black
hole

A generalized spacetime geometry arises in case of low energy string theories.
Among several low energy string theoretical actions, we consider the following
action

$\mathcal{S}_{\mathcal{A}}=\int
d^{4}x\sqrt{-g}e^{-2\phi}(-R-4(\nabla\phi)^{2}+F^{2})$ (41)

where $\phi$ is the dilaton field and $F_{\mu\nu}$ is a Maxwell field
associated with a $U(1)$ subgroup of $E_{8}\times E_{8}$ or
$\text{Spin}(32)/Z_{2}$. Among the family of solutions of the low energy
action given in eq.(41), its charged black hole solution in $(1+1)$-dimensions
is given by[35, 36, 37, 38]

$ds^{2}=-f(r)dt^{2}+g(r)^{-1}dr^{2}$ (42)

where $f(r)$ and $g(r)$ are given by

$\displaystyle f(r)$
$\displaystyle=\left(1-\frac{2Me^{\phi_{0}}}{r}\right)\left(1-\frac{Q^{2}e^{3\phi_{0}}}{Mr}\right)^{-1}~{},$
(43) $\displaystyle g(r)$
$\displaystyle=\left(1-\frac{2Me^{\phi_{0}}}{r}\right)\left(1-\frac{Q^{2}e^{3\phi_{0}}}{Mr}\right)$
(44)

where $\phi_{0}$ is the asymptotic constant value of the dilaton field. This
is the “Garfinkle-Horowitz-Strominger” (GHS) black hole in $(1+1)$-dimensions
[35]. The event horizon radius for a GHS black hole is given by

$r_{+}=2Me^{\phi_{0}}~{}.$ (45)

Using the metric structure from eq.(s)(43,44), we can compute $S_{0}$ from
eq.(14) in the vicinity of the horizon radius $r_{+}$ as

$\begin{split}S_{0}&=\pm\varepsilon\int^{r}_{r_{+}}\frac{dr}{\sqrt{f(r)g(r)}}=\pm\varepsilon\int^{r}_{r_{+}}\frac{dr}{1-\frac{2Me^{\phi_{0}}}{r}}\\\
&\cong\pm 2M\varepsilon e^{\phi_{0}}\ln(r-2Me^{\phi_{0}})=\pm\varepsilon
r_{+}\ln(r-r_{+})~{}.\end{split}$ (46)

As before, we will consider
$r-r_{+}=\mathcal{b}_{1}=\mathcal{a}_{1}e^{i\varphi_{1}}$ with
$\mathcal{a}_{1}=|\mathcal{b}_{1}|$ and
$\varphi_{1}=\text{arg}(\mathcal{b}_{1})$. Using this condition along with
eq.(46), we can write the radial part of the wave function given by
(considering reflection from the event horizon)

$\Phi^{1}(r)=e^{-\frac{i\varepsilon}{\hbar}r_{+}\ln\mathcal{b}_{1}}+\mathcal{R}e^{\frac{i\varepsilon}{\hbar}r_{+}\ln\mathcal{b}_{1}}~{}.$
(47)

Following the same analytical continuation method as before we do a clockwise
rotation by an angle $2\pi$ in the complex plane. Under this clockwise
rotation the radial part of the wave function $\Phi^{1}(r)$ transforms as
follows

$\begin{split}\Phi^{1}_{2\pi}(r)&=e^{-\frac{i\varepsilon}{\hbar}r_{+}\ln\left[\mathcal{b}_{1}e^{-2\pi
i}\right]}+\mathcal{R}e^{\frac{i\varepsilon}{\hbar}r_{+}\ln\left[\mathcal{b}_{1}e^{-2\pi
i}\right]}\\\
&={\xi_{1}}e^{-\frac{i\varepsilon}{\hbar}r_{+}\ln\mathcal{b}_{1}}+\frac{\mathcal{R}}{\xi_{1}}e^{\frac{i\varepsilon}{\hbar}r_{+}\ln\mathcal{b}_{1}}\end{split}$
(48)

where $\xi_{1}$ is given by

$\xi_{1}=\exp\left({-\frac{2\pi\varepsilon}{\hbar}r_{+}}\right)~{}.$ (49)

We observe that the coefficient of the first term in eq.(48) has a value less
than one, therefore by the arguments used in the earlier section we can
conclude the following

$|\mathcal{R}|=\xi_{1}=\exp\left[{-\frac{2\pi\varepsilon}{\hbar}r_{+}}\right]~{}.$
(50)

From eq.(50), we can calculate the probability of reflection from the horizon
of the GHS black hole as follows

$\mathcal{P}=|\mathcal{R}|^{2}=\exp\left[{-\frac{8\pi\varepsilon
GMe^{\phi_{0}}}{\hbar c^{3}}}\right]~{}.$ (51)

In order to obtain the Hawking temperature $T$, we need to use the principle
of detailed balance as before. The final form of the Hawking temperature for a
Garfinkle-Horowitz-Strominger black hole is given by

$T=\frac{\hbar c^{3}}{8\pi k_{B}GM}e^{-\phi_{0}}~{}.$ (52)

From eq.(52), we observe that for a GHS black hole the Hawking temperature has
the usual form of the temperature multiplied by an exponential term. Now for a
very small value of the asymptotic constant value of the dilaton field we can
use $e^{-\phi_{0}}\approx 1-\phi_{0}+\mathcal{O}(\phi_{0}^{2})$ . Hence, the
Hawking temperature takes the form as follows

$T\cong\frac{\hbar c^{3}}{8\pi k_{B}GM}-\frac{\hbar c^{3}\phi_{0}}{8\pi
k_{B}GM}+\mathcal{O}(\phi_{0}^{2})~{}.$ (53)

The form of the Hawking temperature in eq.(53) has the form similar to that of
a base temperature term along with higher order correction terms in
$\phi_{0}$.

Before concluding our discussion, we would like to point out that the same
result for the Hawking temperature can be obtained for a quantum corrected
black hole geometry using the tunneling method [15] using the form of the
action derived earlier in eq.(18). The solution of the covariant Klein-Gordon
equation from eq.(19) can be decomposed into one ingoing solution and one
outgoing solution as follows

$\displaystyle\Psi_{in}=$ $\displaystyle\exp\left[-\frac{i\varepsilon
t}{\hbar}-\frac{i\varepsilon}{\hbar}\int\frac{dr}{\sqrt{f(r)g(r)}}\right]$
(54) $\displaystyle\Psi_{out}=$ $\displaystyle\exp\left[-\frac{i\varepsilon
t}{\hbar}+\frac{i\varepsilon}{\hbar}\int\frac{dr}{\sqrt{f(r)g(r)}}\right]~{}.$
(55)

Here we are considering the case of tunneling across the event horizon of the
black hole. When the particle crosses the horizon, we observe an interchange
of signs between the space and time coordinates indicating an imaginary part
of the time coordinate for crossing the horizon. Two patches across the
horizon are generally connected by an imaginary time[39]. Hence, by using the
ingoing and outgoing wave functions in eq.(s)(54,55), we can express the
ingoing and outgoing solutions of the particle as [15]

$\displaystyle\mathcal{P}_{in}=$
$\displaystyle|\Psi_{in}|^{2}=\exp\left(\frac{2\varepsilon}{\hbar}\Im\left[t\right]+\frac{2\varepsilon}{\hbar}\Im\left[\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\right]\right)~{},$
(56) $\displaystyle\mathcal{P}_{out}=$
$\displaystyle|\Psi_{out}|^{2}=\exp\left(\frac{2\varepsilon}{\hbar}\Im\left[t\right]-\frac{2\varepsilon}{\hbar}\Im\left[\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\right]\right)~{}$
(57)

where $\Im[t]$ denotes the imaginary part of time. Using the principle of
detailed balance [15], we can write the following

$\begin{split}\frac{\mathcal{P}_{out}}{\mathcal{P}_{in}}=\exp\left[-\frac{4\varepsilon}{\hbar}\Im\left[\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\right]\right]=\exp\left[-\frac{\varepsilon}{T}\right]~{}.\end{split}$
(58)

From eq.(58), we can extract the form of the Hawking temperature to be

$T=\frac{\hbar}{4}\left(\Im\left[\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\right]\right)^{-1}~{}.$
(59)

The above equation gives the Hawking temperature for a large class of
spherically symmetric black holes from the viewpoint of tunneling across the
horizon. Substituting the form of the quantum corrected black hole geometry in
eq.(59), we recover eq.(38) thereby agreeing with the interference and
reflection method. It is very important to notice that if the lapse function
does not possess an explicit quantum gravity correction then the Hawking
temperature also has no such corrections.

## 5 Conclusion

In this paper, we have used the method of reflection from the horizon to
investigate the Hawking temperature of a quantum corrected black hole. At
first we have calculated the Hawking temperature of the quantum corrected
black hole using the consideration that there will be a finite probability of
a part of the scalar wave function being reflected from the event horizon of
the black hole. The computed temperature formula for the quantum corrected
black hole has structure similar to that of the Hawking temperature of a
classical black hole (Schwarzschild black hole) along with higher order
quantum gravity corrections. We have also computed the Hawking temperature for
a GHS black hole. For the next part of our analysis, we assume tunneling
through the event horizon of a black hole. We obtained the same result for the
Hawking temperature of the quantum corrected black hole as obtained earlier
using the method of reflection from the horizon. These findings showcase the
fact that quantum corrections arise in the Hawking temperature of the black
hole only due to underlying quantum gravitational effects in the spacetime
geometry. These results also indicate clearly that there should be no quantum
gravity corrections in the Hawking temperature of the black hole if there are
no such quantum gravitational effects in the geometry of the spacetime.

## Appendix: Hawking temperature with the complete solution in eq.(26)

In section (3), we have calculated the form of the Hawking temperature using
the part of the solution outside the event horizon of the black hole in
eq.(26). In this Appendix, we shall proceed with the same calculation keeping
both the terms. Following the discussion after eq.(27), we can write the form
of the radial part of the wavefunction in the vicinity of the event horizon
keeping both the terms in eq.(23) as follows

$\begin{split}\Phi(r)&=e^{-\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{(r_{+}-r_{-})}\ln(r-r_{+})-\frac{r_{-}}{(r_{+}-r_{-})}\ln(r-r_{-})\right)}\\\
&+\mathcal{R}e^{\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{(r_{+}-r_{-})}\ln(r-r_{+})-\frac{r_{-}}{(r_{+}-r_{-})}\ln(r-r_{-})\right)}~{}.\end{split}$
(60)

Again defining $r-r_{+}=\mathcal{b}=\mathcal{a}e^{i\varphi}$, we can recast
the above relation in the following form

$\begin{split}\Phi(r)&=e^{-\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{(r_{+}-r_{-})}\ln\mathcal{b}-\frac{r_{-}}{(r_{+}-r_{-})}\ln(\mathcal{b}+(r_{+}-r_{-}))\right)}\\\
&+\mathcal{R}e^{\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{(r_{+}-r_{-})}\ln\mathcal{b}-\frac{r_{-}}{(r_{+}-r_{-})}\ln(\mathcal{b}+(r_{+}-r_{-}))\right)}~{}.\end{split}$
(61)

Using the same analytic continuation procedure [14] and rotating the complex
number $\mathcal{b}$ by an angle $2\pi$ in the clockwise direction on the
complex plane, we can recast the rotated form of the radial part of the
wavefunction as

$\begin{split}\Phi(r)&=e^{-\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{(r_{+}-r_{-})}\ln\mathcal{b}_{2\pi}-\frac{r_{-}}{(r_{+}-r_{-})}\ln(\mathcal{b}_{2\pi}+(r_{+}-r_{-}))\right)}\\\
&+\mathcal{R}e^{\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{(r_{+}-r_{-})}\ln\mathcal{b}_{2\pi}-\frac{r_{-}}{(r_{+}-r_{-})}\ln(\mathcal{b}_{2\pi}+(r_{+}-r_{-}))\right)}~{}.\end{split}$
(62)

where the form of $\mathcal{b}_{2\pi}$ is given by eq.(32). Since
$\ln(\mathcal{b}_{2\pi}+(r_{+}-r_{-}))=\ln(\mathcal{b}e^{-2\pi
i}+(r_{+}-r_{-}))=\ln(\mathcal{b}+(r_{+}-r_{-}))$, hence eq.(62) takes the
form

$\begin{split}\Phi(r)&=\xi e^{-\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{(r_{+}-r_{-})}\ln\mathcal{b}-\frac{r_{-}}{(r_{+}-r_{-})}\ln(\mathcal{b}+(r_{+}-r_{-}))\right)}\\\
&+\frac{\mathcal{R}}{\xi}e^{\frac{2i\varepsilon
GM}{\hbar}\left(\frac{r_{+}}{(r_{+}-r_{-})}\ln\mathcal{b}_{2\pi}-\frac{r_{-}}{(r_{+}-r_{-})}\ln(\mathcal{b}_{2\pi}+(r_{+}-r_{-}))\right)}\end{split}$
(63)

where the analytical form of $\xi$ is given in eq.(34). By continuing the
analysis in section (3) after eq.(34), the same form of the Hawking
temperature $T$ given in eq.(38) is obtained.

## References

  * [1] A. Einstein, Sitzungsber Preuss Akad Wiss (1915) 844.
  * [2] A. Einstein, Ann. der Physik 49 (1916) 769.
  * [3] D. Christodoulou, Phys. Rev. D 25 (1970) 1596.
  * [4] J.D. Bekenstein, Lett. Nuovo Cimento 4 (1972) 737.
  * [5] J. D. Bekenstein, Phys. Rev. D 7 (1973) 2333.
  * [6] J. D. Bekenstein, Phys. Rev. D 9 (1974) 3292.
  * [7] J. M. Bardeen, B. Carter and S. W. Hawking, Commun. Math. Phys. 31 (1973) 161.
  * [8] S.W. Hawking, Nature 248 (1974) 30.
  * [9] S.W. Hawking, Commun. Math. Phys 43 (1975) 199.
  * [10] S.W. Hawking, Phys. Rev. D 13 (1976) 191.
  * [11] J. B. Hartle and S. Hawking, Phys. Rev. D 13 (1976) 2188.
  * [12] G. W. Gibbons and S. Hawking, Phys. Rev. D 15 (1977) 2752.
  * [13] S. M. Christensen and S. A. Fulling, Phys. Rev. D 15 (1977) 2088.
  * [14] M. Yu. Kuchiev, Phys. Rev. D 69 (2004) 124031.
  * [15] M. K. Parikh and F. Wilczek, Phys. Rev. Lett. 85 (2000) 5042.
  * [16] K. Srinivasan and T. Padmanabhan, Phys. Rev. D 60 (1999) 024007.
  * [17] Q.-Q. Jiang, S.-Q Wu and X. Cai, Phys. Rev. D 73 (2006) 064003.
  * [18] Q.-Q. Jiang, S.-Q Wu and X. Cai, Erratum [Phys. Rev. D 73 (2006) 069902]: Phys. Rev. D 73 (2006) 069902.
  * [19] Y. Hu, J. Zhang and Z. Zhao, Mod. Phys. Lett. A 21 (2006) 2143.
  * [20] Z. Zu and B. Chen, Phys. Rev. D 75 (2007) 024041.
  * [21] C.-Z. Liu and J.-Y. Zhu, Gen. Relativ. Gravit. 40 (2008) 1899.
  * [22] R. Kerner and R. B. Mann, Class. Quant. Gravit. 25 (2008) 095014.
  * [23] D.-Y. Chen, Q.-Q. Jiang and X.-T. Zu, Class. Quant. Gravit. 25 (2008) 205022.
  * [24] R. Di. Criscienzo and L. Vanzo, Euro. Phys. Lett. 82 (2008) 60001.
  * [25] R. Li and J.-R. Ren, Phys. Lett. B 661 (2008) 370.
  * [26] R. Kerner and R. B. Mann, Phys. Lett. B 665 (2008) 277.
  * [27] D.-Y. Chen, Q.-Q. Jiang and X.-T. Zu, Phys. Lett. B 665 (2008) 106.
  * [28] S. A. Kaplan, Zh. Eksp. Teor. Fiz. 19 (1949) 951; arXiv: 2201.07971 [physics.hist-ph].
  * [29] Ya. B. Zel’dovich and I. D. Novikov, Dokl. Akad. Nauk. SSSR 155 (1964) 1033.
  * [30] M. Reuter, Phys. Rev. D 57 (1998) 971.
  * [31] C. Wetterich, Phys, Lett. B 301 (1993) 90.
  * [32] M. Reuter and C. Wetterich, Nucl. Phys. B 417 (1994) 181.
  * [33] F. Saueressig and M. Reuter, “Quantum Gravity and the Functional Renormalization Group: The Road towards Asymptotic Safety (Cambridge monographs on mathematical physics) ”, Cambridge University Press (2019).
  * [34] A. Bonanno and M. Reuter, Phys. Rev. D 62 (2000) 043008.
  * [35] D. Garfinkle, G. T. Horowitz, and A. Strominger, Phys. Rev. D 43 (1991) 3140; Erratum: Phys. Rev. D 45 (1992) 3888.
  * [36] G. W. Gibbons, Nucl. Phys. B 207 (1982) 337.
  * [37] G. W. Gibbons and K. Maeda, Nucl. Phys. B 298 (1988) 741.
  * [38] S. Gangopadhyay and S. Kulkarni, Phys. Rev. D 77 (2008) 024038.
  * [39] E. T. Akhmedov, T. Pilling and D. Singleton, Int. J. Mod. Phys. D 17 (2008) 2453.


# Mechanism of electron-beam manipulation of single dopant atoms in silicon

Alexander Markevich  [<EMAIL_ADDRESS>Bethany M Hudak  [
Jacob Madsen  [  Jiaming Song  [  Paul C Snijders  [  Andrew R Lupini  [
<EMAIL_ADDRESS>Toma Susi  [<EMAIL_ADDRESS>
###### Abstract

The precise positioning of dopant atoms within bulk crystal lattices could
enable novel applications in areas including solid-state sensing and quantum
computation. Established scanning probe techniques are capable tools for the
manipulation of surface atoms, but at a disadvantage due to their need to
bring a physical tip into contact with the sample. This has prompted interest
in electron-beam techniques, followed by the first proof-of-principle
experiment of bismuth dopant manipulation in crystalline silicon. Here, we use
first principles modeling to discover a novel indirect exchange mechanism that
allows electron impacts to non-destructively move dopants with atomic
precision within the silicon lattice. However, this mechanism only works for
the two heaviest group V donors with split-vacancy configurations, Bi and Sb.
We verify our model by directly imaging these configurations for Bi, and by
demonstrating that the promising nuclear spin qubit Sb can be manipulated
using a focused electron beam.

univie] University of Vienna, Faculty of Physics, Boltzmanngasse 5, 1090
Vienna, Austria NRL]Naval Research Laboratory, Material Sciences and
Technology, 4555 Overlook Ave SW, Washington, DC 20375, USA univie]University
of Vienna, Faculty of Physics, Boltzmanngasse 5, 1090 Vienna, Austria
NWU]School of Physics, Northwest University, 1 Xuefu Avenue, Xi’an, Shaanxi
710127, China ORNL1]Materials Science and Technology Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37830, USA ORNL2]Center for
Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge,
Tennessee 37830, USA univie]University of Vienna, Faculty of Physics,
Boltzmanngasse 5, 1090 Vienna, Austria DISTRIBUTION A: Approved for public
release, distribution is unlimited.

Notice: This manuscript has been authored by UT-Battelle, LLC, under Contract
No. DE-AC0500OR22725 with the U.S. Department of Energy. The United States
Government retains and the publisher, by accepting the article for
publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for the
United States Government purposes. The Department of Energy will provide
public access to these results of federally sponsored research in accordance
with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-
access-plan).

##  1 Introduction

Quantum information processing has surged in recent decades from a speculative
idea 1 into a technological race for quantum supremacy 2. Multiple viable
architectures with various strengths and weaknesses have been proposed 3,
including solid-state qubits. Due to our reliance and technological
capabilities in silicon manufacturing, the nuclear spins of positively charged
donor atoms within crystalline silicon have attracted particular interest 4,
5. However, two factors make such qubits challenging: the precise placement of
dopants into an ordered array, and their individual addressing for quantum
logic operations. The mere detection of single atoms inside a semiconductor is
extremely difficult, and there are only a few techniques capable of this
feat—particularly, atom probe tomography 6 and aberration-corrected scanning
transmission electron microscopy 7 (STEM).

Phosphorus (31P) as a dopant has attracted the greatest interest due to its
long nuclear spin coherence lifetime, enabling it to function as a quantum
memory 8. Typically, traditional non-deterministic ion implantation and post-
hoc lithography have been used for fabricating devices 9. A degree of
deterministic placement control has been enabled by scanning tunneling
microscopy based hydrogen depassivation lithography of Si surfaces 10,
followed by molecular beam epitaxy growth to complete the crystal over a
chemically introduced donor 11. However, the spin dephasing time for the
hyperfine-coupled electron on P sites is rather short 12, hindering its
electric-field control. Very recently, antimony (Sb) has been suggested as a
promising alternative, with two orders of magnitude higher dephasing time and
robust electrical control via its nuclear quadrupole moment 13.

Despite some advances, dopant placement remains the major hurdle for further
progress. Electron-beam manipulation of covalently bound lattice impurities
was recently proposed as a new kind of atomically precise manipulation tool
14, 15. Initially, this was limited to impurities in carbon nanomaterials
including graphene 16, 17, 18, 19 and single-walled carbon nanotubes 20, but a
similar capability was recently also demonstrated for heavy Bi dopants in
silicon 21, 22. However, the atomistic mechanism for their movement within the
crystal was not clear, nor was it known whether other dopants could be
affected, which beam energies are optimal, nor whether the surrounding lattice
is irreparably damaged.

The electron-beam energy in a (S)TEM is typically in the range of 30–-300 keV,
while displacement threshold energies are on the order of 10 eV. However,
conservation of momentum implies that only a very small fraction of the
kinetic energy can be transferred from the fast electron to a nucleus in a
single elastic collision. Electronic interactions or multiple collisions might
also contribute depending on the material and imaging conditions. In an
insulating sample, ionization or “radiolysis” would complicate the picture 23,
but the electrical conductivity of highly doped Si should mitigate such
damage. Moreover, in a recent study of atomic displacements in MoS2 under
electron irradiation, 24 it was demonstrated that electronic excitations
strongly affected displacement threshold energies at lower accelerating
voltages but made only a minor contribution at a primary beam energy of 80 kV.
The relatively low beam current typical for aberration-corrected STEM suggests
that sample heating is negligible because of the thermal conductivity of the
sample 25. The dominant beam interaction mechanism in these samples is
therefore expected to be elastic collisions transferring between 10–20 eV from
single fast incident electrons to individual nuclei. The knock-on threshold
energy for bulk Si has been calculated 26 at about 12.5 eV, corresponding to
an electron beam energy of about 140 keV.

In this study, we combine density functional theory molecular dynamics
(DFT/MD) and aberration-corrected STEM to systematically study the behaviour
of group V dopants in silicon under electron irradiation. Similar to
heteroatom-doped carbon systems 27, 14, 19, momentum-conserving electron
scattering at modest primary beam energies can only affect the lighter Si
neighbors of the impurity atom, allowing its movement to be controlled by
directing the electron beam at a chosen lattice site. Since any beam-induced
dynamics occur on timescales that are far below the experimental time
resolution 14, modeling is required to understand the atomistic details of the
process.

We reveal that dopant manipulation in crystalline silicon proceeds via a novel
type of nondestructive mechanism we call _indirect exchange_ : the impacted Si
is knocked into a metastable interstitial configuration, with the impurity
taking its original lattice position in the subsequent recombination. Contrary
to the process in graphene, the primary knock-on atom does not end up as a
neighbor of the impurity, but rather as a next-nearest neighbor, displacing
another Si atom during the dynamics.

Our simulations show that this mechanism only works with the heavier two
pnictogen dopants, Sb and Bi, which induce greater lattice distortion, and
whose split-vacancy configuration in the case of Bi we are occasionally able
to image directly, whereas for the lighter two, As and P, the ejected Si
interstitial recombines with a vacancy without exchange of atom positions. To
confirm these predictions, we demonstrate for the first time the manipulation
of Sb impurities in a thin crystalline slab of silicon, finding that similarly
to Bi these can be successfully manipulated using focused electron
irradiation.

##  2 Methods

###  2.1 Computational Methods

####  2.1.1 DFT simulations

All simulations were performed using density functional theory (DFT) as
implemented in the GPAW package 28. We used the PBE exchange-correlation
functional 29, a localized (dzp) basis set 30, and a real-space grid spacing
of 0.2 Å. A cubic supercell of 216 atoms was used in all calculations. The
Brillouin zone was sampled using a $3\times 3\times 3$ k-point mesh according
to the Monkhorst-Pack scheme 31.

It should be noted that while manipulation of dopants has been performed in Si
slabs of $\sim$10-–20 nm thickness, for simulations we used a periodic
supercell representing an infinite crystal. In a sub-surface layer of a real
sample, lattice strains induced by the surface layer reconstruction can affect
the dynamic behavior of crystal defects. Considerable lattice strain has been
shown to extend for five atomic layers into the crystal lattice, 32, 33 which
in the case of Si is $\sim$0.5 nm. In our experiments on dopant manipulation,
we have estimated typical depths of the dopant atoms that have been probed at
about 3-–5 nm from the surface. 21 At such a distance from the surface, the
lattice strain induced by the surface reconstruction should be negligible and
the phenomena observed experimentally are characteristic for a bulk crystal
volume.

To simulate the dynamics of the atoms induced by electron irradiation, we used
Velocity-Verlet (NVE) molecular dynamics 34 with a time step of 0.5 fs. The
effect of momentum-conserving electron scattering on Si neighbours of a dopant
was modelled by assigning a momentum that corresponds to a specific
transferred energy to one Si atom. Threshold energies were estimated by
varying the amount of energy transferred to the Si atom at 0.1 eV intervals. A
detailed description of the methodology can be found in Ref. 35.

Heating of systems with As and P dopants and interstitial Si was performed
using Langevin (NVT) dynamics 36 up to 1200 K or until any structural changes
occurred. Random snapshots of the system were then taken at temperatures in
the range of 800–1200 K and six Velocity-Verlet simulations were performed for
each dopant until the recombination of the defect occurred. The energy barrier
for the recombination of the Bi-V complex with a Si interstitial was
calculated with the climbing-image nudged elastic band (cNEB) method 37 using
nine system images with a spring constant of 0.1 eV/Å, relaxed until all
atomic forces were below 0.02 eV/Å.

####  2.1.2 Image simulations

All image simulations were performed using the multislice code abTEM 38. Based
on the experimental HAADF contrast, we assumed a slab thickness of 15 nm and
that the dopants were near the center of the slab. The size of the structure
prohibits fully optimizing all the positions at the _ab initio_ level of
theory; hence, smaller supercells around the defects were optimized and
inserted into an idealized slab. The acceleration voltage, beam convergence
semi-angle, and detector semi-angles were set to the experimental values. The
remaining unknown imaging parameters were estimated by matching the image
contrast through a grid search. The optimal parameters were a defocus of 180
Å, a spherical aberration of 40 µm, beam tilt in $x$ and $y$ of 13 and 9 mrad,
and a source spread of 0.8 Å (FWHM). We used the PRISM algorithm 39 to quickly
simulate a large number of required images for the grid search. The final
image was simulated using the standard multislice algorithm with thermal
diffuse scattering included in the frozen phonon approximation.

###  2.2 Experimental Materials and Methods

####  2.2.1 Sample growth

Sb-doped Si films were grown on Si(100) substrates. The Si(100) surface was
prepared with a (2$\times$1) reconstruction using standard ultra-high-vacuum
degassing and flashing procedures, and surface quality examined by scanning
tunneling microscopy and low-energy electron diffraction. On the
Si(100)-2$\times$1 surface, approximately 0.5 monolayers (ML) of Sb atoms was
subsequently deposited with the substrate held at 550 °C, using a standard
effusion cell held at 400 °C. This Sb dose, higher than desirable for qubit
applications, was chosen in order to have a sufficient doping concentration
for good visibility and easy dopant manipulation. Next, a nominally $\sim$5 ML
Si capping layer was grown with the substrate held at room temperature and
annealed at 450 °C for 5 s. A final capping layer of 20 ML Si was grown with
the substrate held at room temperature and annealed at 720 °C for 1 s to
crystallize.

####  2.2.2 Specimen preparation

Cross-sectional STEM samples were prepared by taking the as-grown sample,
cutting it in half to produce a sandwich, and then slicing it into sections
that were glued to a support grid. The samples were mechanically polished to
approximately 10 µm, then ion-milled to electron transparency. Immediately
before imaging, the samples were baked at a nominal 140 °C in vacuum overnight
to reduce surface contamination. Thickness for typical images was measured by
electron energy-loss spectroscopy as about 0.1–0.2 mean free paths (about
10–20 nm), although both thicker and thinner regions were examined. Images
taken far away from the film and glue indicate the existence of an amorphous
oxide layer about 1 nm in thickness. After ion-milling, the Si surface was
chemically and structurally disordered. We found that there was often apparent
motion in the images, due to the surface oxides and overlying amorphous
carbon. Repeated intense illumination of the same area usually resulted in
milling of the surface Si.

####  2.2.3 Scanning transmission electron microscopy

Aberration-corrected STEM was performed in a Nion UltraSTEM 200 operated at
160 keV at a low beam current of 20$\pm$2 pA and a nominal beam convergence
semiangle of 30 mrad. We acquired images with a high-angle annular dark-field
(HAADF) detector with a nominal semiangular range of 60–400 mrad under a
variety of dose conditions, including field of view and scan speed. An ImageJ
smoothing filter that replaces each pixel with the average of its 3$\times$3
neighbors was applied to Fig. 3 to enhance contrast. Sb impurities are
identified by their characteristic $Z$-contrast, and their presence in the
sample was also verified by electron energy loss spectroscopy collected at 200
keV with a dispersion of 0.5 eV/ch and an acquisition time of 1 s/pixel. The
additional minor EEL intensity in Fig. 3 starting at $\sim$766 eV likely
corresponds to the Sb $M_{3}$ edge, while the more intense features starting
at $\sim$931 eV can likely be attributed to the $L_{2,3}$ edge of Cu
contamination near the disordered surface layer.

##  3 Results and discussion

### Uncovering the manipulation mechanism

To uncover the mechanism of directed dopant movement in silicon, we performed
DFT/MD simulations to model the effect of momentum-conserving elastic electron
scattering (ie. knock-on impacts) on the lattice atoms. We distinguish two
types of lattice jumps for the impurities: either diagonally or laterally
between the atomic columns as viewed along the [110] zone axis. While such
jumps are equivalent from the perspective of lattice symmetry, we observe
slight differences in terms of their electron-beam induced dynamics. Similar
to the case of graphene, directed manipulation is made possible by the fact
that beam-induced dynamics are triggered by an electron impact on one of the
lattice atoms neighboring the impurity.

Let us first consider the case of the heaviest stable group V dopant, bismuth.
In the original experimental work on Bi manipulation in crystalline silicon,
the authors used electron energies of 60, 100, 160 and 200 keV 21, 22. For 60
and 100 keV, almost no changes in the structure were observed, while for 160
keV, the positioning of dopants was possible by rapidly dragging a stationary
electron probe over multiple atomic columns. This clearly indicated that a
knock-on process was involved; indeed, originally the movement was thought to
occur via the sputtering of Si atoms by the beam, followed by impurity
migration via the created vacancies. As we show below, the actual mechanism is
somewhat more complicated—but also more favorable for electron-beam
manipulation.

For 160 keV primary beam energy, the maximum kinetic energy that an
elastically back-scattering electron can transfer to a stationary Si atom is
$\sim$14.5 eV, with the transferred momentum directed along the beam
direction. Figure 1a-f shows snapshots from a DFT/MD trajectory where a
diagonal Si neighbour of Bi is assigned a velocity that corresponds to 14.5 eV
of kinetic energy in the [110] direction. The impacted Si atom moves to an
interstitial configuration (Fig. 1d), leaving behind a bismuth-vacancy (Bi-V)
defect. This configuration is metastable, with the impurity atom further
taking over the lattice position originally occupied by the impacted Si atom,
while the latter becomes its second-nearest neighbour, displacing another Si
atom to fill the vacant lattice site. In many of our MD trajectories, this
occurred within the simulation time of 1–2 ps. We further estimate a nudged
elastic band (NEB) energy barrier of only $\sim$70 meV for the recombination.

Figure 1: Indirect exchange of Bi and Si lattice positions induced by an
electron impact. Snapshots from DFT/MD trajectories showing the jump of a Bi
impurity, shown in pink, by one lattice site through the indirect exchange
mechanism. The Si atom shown in black has been assigned a momentum
corresponding to a kinetic energy of 14.5 eV in the [110] lattice direction as
indicated by the arrow. The Bi impurity is shown in pink color and Si atoms
within the bulk in light grey. We distinguish two types of impurity jumps:
either diagonally (a-f) or laterally (g-l) between the atomic columns as
viewed along the [110] zone axis of the experimental geometry (here, the
dynamics are visualized from the side for clarity). During the dynamics, the
Bi moves to the lattice position originally occupied by the impacted Si, while
the latter becomes its second-nearest neighbour upon recombination, pushing
the Si neighbor of the impurity shown in dark grey into the lattice site
vacated by the dopant. (See also Supporting Information.)

Notably, no atoms are—nor need to be—sputtered from the structure. Thus, this
process is similar to the direct exchange mechanism discovered earlier for
substitutional dopants in graphene 14, 19, where the impurity atom exchanges
places with a neighboring lattice C atom. However, in the case of bulk
silicon, the impacted Si neighbour of the impurity instead ends up as its
second-nearest neighbour, and thus we call this process indirect exchange.

The dynamics discussed thus far result in the Bi dopant moving diagonally
between the columns as viewed along the [110] zone axis, and such lattice
dynamics occur for transferred energies above 13.6 eV. Below this energy, the
impacted Si atom does not reach the interstitial configuration and quickly
returns to its initial position. Above $\sim$22 eV, the dynamics becomes more
complex, with multiple atoms displaced from their lattice sites.

For lateral Bi jumps across the columns (Fig. 1g-l), the dopant has two Si
neighbours that lie in the same atomic column along the probe direction. Our
simulations show that only impacts on the Si positioned closer to the probe
and thus before the plane of the impurity can result in the movement of the
dopant between the columns through the indirect exchange mechanism. The
calculated threshold energy for this process is 13.9 eV, slightly higher than
that for jumps between the columns.

It should be noted that lateral jumps also result in the change of the dopant
position perpendicular to the (110) plane by one lattice site. This provides
the possibility to control the depth of a dopant within the slab by repeated
lateral jumps, but only in the direction towards the electron source (although
measuring this is challenging). Conversely, manipulation strictly within the
same plane is only possible for diagonal jumps.

### Differences amongst group V dopants

To explore the potential of the electron-manipulation technique further, we
systematically modelled the remaining group V dopants in silicon. Our
simulations show that Sb impurities can be manipulated in the same way as Bi,
with lattice jumps occurring through the indirect exchange mechanism.
Threshold energies for diagonal and lateral Sb movement were calculated to be
13.0 eV and 13.6 eV, respectively. However, in the case of the lighter
impurities, As and P, no indirect exchange process was observed in our
simulations. Instead, in some of the trajectories where the impacted Si atom
moved to an interstitial configuration, recombination occurred with no change
in the dopant’s lattice position. We further performed simulations at elevated
temperatures in the range of 800–1200 K with a starting configuration where
the Si atom is at the interstitial site. In all simulations, the recombination
occurred without the movement of the dopant atom to the target lattice site.

Such a clear difference in the behaviour of larger and smaller dopant atoms
can be explained by the atomic structure of their vacancy complexes. Oversized
impurities, such as Bi and Sb, form a split-vacancy structure, whereas As and
P adopt a substitutional dopant-vacancy geometry instead 40. Similar
structures are formed during the indirect exchange process. When a Si
neighbour is knocked into an interstitial configuration, Bi and Sb dopants
move towards the newly vacant lattice site (Fig. 2a), opening an easy route
for recombination, where only one additional lattice Si atom needs to be
displaced. In contrast, As and P remain close to their original lattice
positions even when their neighbor is displaced (Fig. 2b). MD simulations at
elevated temperatures suggest that in this case, the recombination of the
interstitial Si, which involves the displacement of two additional lattice
atoms, turns out to be energetically more favourable than the jump of the
dopant.

Remarkably, in the Bi-doped sample, we were occasionally able to record
intense atomic contrast between the diagonal Si columns, far too high to
correspond to a Si atom, which, even if stationary, would not be visible
compared to the lattice contrast. However, this atomic contrast is consistent
with a single Bi atom as verified by quantitative image simulations (Fig.
2c-d). We thus believe that the impacted Si interstitial can occasionally
diffuse away from the site instead of recombining, allowing the split-vacancy
configuration to be directly recorded. This idea is supported by the fact that
such contrast was observed more frequently at 200 keV, where the electron beam
can transfer more energy to the impacted Si atom, facilitating its escape.

Figure 2: Group V dopant-vacancy complexes in Si. (a,b) Atomic structures of
dopant-vacancy complexes with a Si interstitial optimised with DFT for (a)
heavier Bi with a split-vacancy configuration (Sb appears similar), and for
(b) lighter As with substitutional dopant-vacancy configuration (P appears
similar). (c) HAADF image (Gaussian filtered) and (d) image simulation of a Bi
split-vacancy near a Bi substitution in Si(110) at 160 keV. The overlaid line
profiles indicate the intensity plotted over (e) the Si column, Bi split
vacancy, and Si column direction (red) and (f) the Bi split vacancy and Si
column with Bi substitution direction (blue), with an excellent quantitative
match obtained. (See also Supporting Information.)

### Electron-beam manipulation of antimony

To confirm our theoretical predictions, we attempted to experimentally
manipulate Sb impurities using the focused electron probe of an aberration-
corrected Nion UltraSTEM200 operating at 160 keV. By ’dragging’ the beam along
atomic columns adjacent to the Sb-containing column, we can controllably
direct the position of the Sb dopant in the same manner that Bi atoms were
directed 21. Since the beam cannot be simultaneously scanned and manually
directed, an image (Fig. 3c) is acquired and used as a survey image. The beam
is manually positioned over a series of neighboring atomic columns in the
survey image, and then a second scanned image is acquired (Fig. 3d), revealing
the new location of the Sb atom. Fig. 3e-j illustrate the path traced by the
beam. As in the case of Bi 21, the Sb has moved several columns away from its
initial position, following the path that the operator directed with the
electron beam. (The additional bright contrast near the endpoint is likely due
to dopants at different depths being simultaneously dragged into position,
which is a challenge for 3D manipulation.) White circles in Fig. 3 indicate Sb
atoms that remain stationary during the process, highlighting that the
movement is not simply an overall drift. Thus, we are able to guide the
position of the Sb dopant at 160 keV, in agreement with our DFT/MD prediction.

Figure 3: Antimony manipulation. Identification and electron-beam manipulation
of Sb dopants in Si(110). (a) HAADF image of the Sb-doped Si sample acquired
at 200 keV, with a green overlaid rectangle corresponding to the electron
energy-loss spectrum image shown in the middle. The spectrum in (b) was summed
over the 15$\times$10 pixels indicated by the red rectangle to enhance the
signal to noise ratio, showing a prominent Sb $M_{4,5}$ edge starting at
$\sim$528 eV (see Methods). (c) HAADF close-up of the region at the lower end
of the green rectangle in (a) captured before manipulation at 160 keV. Yellow
circles indicate the targeted Sb atom, and white circles mark Sb atoms that
remain stationary during the process. (b) Image acquired after direct Sb
positioning, with the manipulation path traced as a white line. (e-j) Screen-
captured frames showing the position of the electron beam as it was used to
guide the Sb dopant. Scale bar is 1 nm. (See also Supporting Information.)

### Influence of the atom emission angle

Our simulations indicate that due to the Pauli repulsion between their
electron clouds, the impacted Si atom moving along the [110] lattice direction
glances off either a neighbouring Si atom (Fig. 1b,c) or the impurity (Fig.
1h,i), which results in the change of its direction of motion and significant
loss of kinetic energy. This loss can be reduced if the initial velocity of
the impacted Si is directed slightly away from the neighbouring lattice site.
Indeed, our simulations show that for the case of diagonal Bi jumps, a change
of the emission angle $\Theta$ from 0 to 15° (see inset of Fig. 4) results in
a significant decrease of the threshold energy for indirect exchange—from 13.6
to 9.4 eV. It further drops to 8.3 eV for an emission angle of 30°, but
remains the same for 45°. Similar behavior is also observed for lateral jumps,
with differences between 0–0.6 eV between the two directions and for the Sb
dopant (see Table 1).

Table 1:  Calculated threshold energies (in eV) for the diagonal and lateral jumps of the Bi and Sb dopants as a function of the emission angle of the impacted Si atom. $\Theta$ (°) | 0 | 5 | 10 | 15 | 30 | 45
---|---|---|---|---|---|---
Bi (diag.) | 13.6 | 11.6 | 10.3 | 9.4 | 8.3 | 8.3
Bi (lat.) | 13.9 | 11.8 | 10.3 | 9.4 | 8.2 | 8.7
Sb (diag.) | 13.0 | 11.2 | 10.1 | 9.3 | 8.3 | 8.3
Sb (lat.) | 13.6 | 11.6 | 10.3 | 9.4 | 8.2 | 8.6

However, increase of the emission angle $\Theta$ also implies that the
electron scattering angle is reduced 41, $\theta=\pi-2\Theta$, and therefore
the amount of kinetic energy transferred to the Si by the probe electron
becomes smaller. In Fig. 4 we compare calculated threshold energies of Bi and
Sb dopants averaged between diagonal and lateral directions for different
emission angles (points) with the maximum kinetic energy that can be
transferred to a stationary Si atom in the corresponding direction of motion
at different electron energies (solid lines). For small emission angles, the
decrease of the threshold energy is significantly faster than the decrease in
the transferred energy. This suggest that an electron energies around 120 keV
should be sufficient for manipulation, and might present a useful future
trade-off between speed and control. The direction of the energy transfer and
thus manipulation probability could potentially also be controlled by tilting
the beam 19.

We should also address the question of surface sputtering. Our simulations for
the $(1\times 1)$ Si surface reconstruction show that the threshold energy for
sputtering Si atoms in the direction normal to the surface is only $\sim$8.5
eV. We note, however, that surface sputtering depends on the surface
reconstruction, oxidation and contamination, and therefore this value is
merely an idealized estimate. Surface sputtering can be reduced by decreasing
the electron energy, but it seems difficult to completely avoid this process
for energies at which efficient manipulation can be performed.

Figure 4: Manipulation energetics. Calculated threshold energies for the
indirect exchange for Bi (black squares) and Sb (open triangles) dopants
(average of values for diagonal and lateral jumps (see Table 1), as a function
of the impacted Si atom (colored black in the inset) emission angle $\Theta$.
Solid lines show the maximum kinetic energy that can be transferred to a
stationary Si atom in the corresponding emission direction by probe electrons
at different accelerating voltages.

## Conclusions

Solid-state quantum computers using nuclear spin qubits in silicon require the
fabrication of a large ordered arrays of precisely positioned single donor
atoms with a spacing of tens of nm, buried tens of nm below the crystal
surface. This is a daunting challenge for traditional ion implantation and
lithography techniques and remains the greatest hurdle for further progress.
Deterministic ion implantation has been proposed as a possible solution 42,
but the achievable positional precision seems limited 43. Hydrogen-resist
lithography in combination with scanning tunneling microscopy can achieve the
required precision 11, but qubit fabrication is complex and the choice of
donor species is limited. Direct electron-beam manipulation of donor elements
with atomic precision within the crystal lattice is thus a highly promising
alternative. With process optimization both in terms of sample synthesis and
tuning the primary beam energy, the manipulation of Sb dopants that we
demostrate here may open a path to the deterministic and potentially scalable
fabrication of solid-state nuclear spin qubits in silicon.

The authors thank J. Meyers for the cross-section specimen preparation and J.
Kotakoski for useful discussions.

A.M., J.M., and T.S. were supported by the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(Grant agreement No. 756277-ATMEN), and acknowledge computational resources
provided by the Vienna Scientific Cluster. Electron microscopy work (B.M.H.
and A.R.L.) was supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences, Division of Materials Sciences and Engineering, and
sample growth (P.C.S. and J.S.) by the Laboratory Directed Research and
Development Program of Oak Ridge National Laboratory, managed by UT-Battelle,
LLC for the U.S. Department of Energy, and performed at the Oak Ridge National
Laboratory’s Center for Nanophase Materials Sciences (CNMS), a U.S. Department
of Energy Office of Science User Facility.

Supporting Information accompanying the article include a screen capture movie
file of the Sb manipulation sequence shown in Fig. 3, and a rendered video
illustration of the indirect exchange mechanism following Fig. 1. Open data of
the simulated molecular dynamics trajectories, thermal MD trajectories, and
relaxed split vacancy atomic configurations, as well as raw experimental
scanning transmission electron microscopy images shown in Fig. 2 can be found
on the University of Vienna institutional repository Phaidra
(doi:10.25365/phaidra.275). Open code to reproduce the image simulations in
Fig. 2 can be found on GitHub.

## References

  * Benioff 1980 Benioff, P. The computer as a physical system: A microscopic quantum mechanical Hamiltonian model of computers as represented by Turing machines. _Journal of Statistical Physics_ 1980, _22_ , 563–591
  * Arute et al. 2019 Arute, F.; Arya, K.; Babbush, R.; Bacon, D.; Bardin, J. C.; Barends, R.; Biswas, R.; Boixo, S.; Brandao, F. G. S. L.; Buell, D. A. et al. Quantum supremacy using a programmable superconducting processor. _Nature_ 2019, _574_ , 505–510
  * DiVincenzo 2000 DiVincenzo, D. P. The Physical Implementation of Quantum Computation. _Fortschritte der Physik_ 2000, _48_ , 771–783
  * Kane 1998 Kane, B. E. A silicon-based nuclear spin quantum computer. _Nature_ 1998, _393_ , 133–137
  * Morton et al. 2011 Morton, J. J. L.; McCamey, D. R.; Eriksson, M. A.; Lyon, S. A. Embracing the quantum limit in silicon computing. _Nature_ 2011, _479_ , 345–353
  * Kelly and Miller 2007 Kelly, T. F.; Miller, M. K. Atom probe tomography. _Review of Scientific Instruments_ 2007, _78_ , 031101
  * Nellist et al. 2004 Nellist, P. D.; Chisholm, M. F.; Dellby, N.; Krivanek, O. L.; Murfitt, M. F.; Szilagyi, Z. S.; Lupini, A. R.; Borisevich, A.; Sides, W. H.; Pennycook, S. J. Direct Sub-Angstrom Imaging of a Crystal Lattice. _Science_ 2004, _305_ , 1741
  * Morton et al. 2008 Morton, J. J. L.; Tyryshkin, A. M.; Brown, R. M.; Shankar, S.; Lovett, B. W.; Ardavan, A.; Schenkel, T.; Haller, E. E.; Ager, J. W.; Lyon, S. A. Solid-state quantum memory using the 31 P nuclear spin. _Nature_ 2008, _455_ , 1085–1088
  * Pla et al. 2012 Pla, J. J.; Tan, K. Y.; Dehollain, J. P.; Lim, W. H.; Morton, J. J. L.; Jamieson, D. N.; Dzurak, A. S.; Morello, A. A single-atom electron spin qubit in silicon. _Nature_ 2012, _489_ , 541–545
  * Randall et al. 2009 Randall, J. N.; Lyding, J. W.; Schmucker, S.; Von Ehr, J. R.; Ballard, J.; Saini, R.; Xu, H.; Ding, Y. Atomic precision lithography on Si. _Journal of Vacuum Science & Technology B: Microelectronics and Nanometer Structures Processing, Measurement, and Phenomena_ 2009, _27_ , 2764–2768
  * Fuechsle et al. 2012 Fuechsle, M.; Miwa, J. A.; Mahapatra, S.; Ryu, H.; Lee, S.; Warschkow, O.; Hollenberg, L. C. L.; Klimeck, G.; Simmons, M. Y. A single-atom transistor. _Nature Nanotechnology_ 2012, _7_ , 242–246
  * Muhonen et al. 2014 Muhonen, J. T.; Dehollain, J. P.; Laucht, A.; Hudson, F. E.; Kalra, R.; Sekiguchi, T.; Itoh, K. M.; Jamieson, D. N.; McCallum, J. C.; Dzurak, A. S. et al. Storing quantum information for 30 seconds in a nanoelectronic device. _Nature Nanotechnology_ 2014, _9_ , 986–991
  * Asaad et al. 2020 Asaad, S.; Mourik, V.; Joecker, B.; Johnson, M. A. I.; Baczewski, A. D.; Firgau, H. R.; Madzik, M. T.; Schmitt, V.; Pla, J. J.; Hudson, F. E. et al. Coherent electrical control of a single high-spin nucleus in silicon. _Nature_ 2020, _579_ , 205–209
  * Susi et al. 2014 Susi, T.; Kotakoski, J.; Kepaptsoglou, D.; Mangler, C.; Lovejoy, T. C.; Krivanek, O. L.; Zan, R.; Bangert, U.; Ayala, P.; Meyer, J. C. et al. Silicon–Carbon Bond Inversions Driven by 60-keV Electrons in Graphene. _Phys. Rev. Lett._ 2014, _113_ , 115501
  * Kalinin et al. 2016 Kalinin, S.; Borisevich, A.; Jesse, S. Fire up the atom forge. _Nature_ 2016, _539_ , 485–487
  * Susi et al. 2017 Susi, T.; Meyer, J.; Kotakoski, J. Manipulating low-dimensional materials down to the level of single atoms with electron irradiation. _Ultramicroscopy_ 2017, _180_ , 163–172
  * Dyck et al. 2017 Dyck, O.; Kim, S.; Kalinin, S. V.; Jesse, S. Placing single atoms in graphene with a scanning transmission electron microscope. _Applied Physics Letters_ 2017, _111_ , 113104
  * Tripathi et al. 2018 Tripathi, M.; Mittelberger, A.; Pike, N. A.; Mangler, C.; Meyer, J. C.; Verstraete, M. J.; Kotakoski, J.; Susi, T. Electron-Beam Manipulation of Silicon Dopants in Graphene. _Nano Letters_ 2018, _18_ , 5319–5323
  * Su et al. 2019 Su, C.; Tripathi, M.; Yan, Q.-B.; Wang, Z.; Zhang, Z.; Hofer, C.; Wang, H.; Basile, L.; Su, G.; Dong, M. et al. Engineering single-atom dynamics with electron irradiation. _Science Advances_ 2019, _5_ , eaav2252
  * Mustonen et al. 2019 Mustonen, K.; Markevich, A.; Tripathi, M.; Inani, H.; Ding, E.-X.; Hussain, A.; Mangler, C.; Kauppinen, E. I.; Kotakoski, J.; Susi, T. Electron-Beam Manipulation of Silicon Impurities in Single-Walled Carbon Nanotubes. _Advanced Functional Materials_ 2019, _29_ , 1901327
  * Hudak et al. 2018 Hudak, B. M.; Song, J.; Sims, H.; Troparevsky, M. C.; Humble, T. S.; Pantelides, S. T.; Snijders, P. C.; Lupini, A. R. Directed Atom-by-Atom Assembly of Dopants in Silicon. _ACS Nano_ 2018, _12_ , 5873–5879
  * Jesse et al. 2018 Jesse, S.; Hudak, B. M.; Zarkadoula, E.; Song, J.; Maksov, A.; Fuentes-Cabrera, M.; Ganesh, P.; Kravchenko, I.; Snijders, P. C.; Lupini, A. R. et al. Direct atomic fabrication and dopant positioning in Si using electron beams with active real-time image-based feedback. _Nanotechnology_ 2018, _29_ , 255303
  * Susi et al. 2019 Susi, T.; Meyer, J. C.; Kotakoski, J. Quantifying transmission electron microscopy irradiation effects using two-dimensional materials. _Nature Reviews Physics_ 2019, _1_ , 397–405
  * Kretschmer et al. 2020 Kretschmer, S.; Lehnert, T.; Kaiser, U.; Krasheninnikov, A. V. Formation of Defects in Two-Dimensional MoS2 in the Transmission Electron Microscope at Electron Energies below the Knock-on Threshold: The Role of Electronic Excitations. _Nano Letters_ 2020, _20_ , 2865–2870
  * Börner et al. 2016 Börner, P.; Kaiser, U.; Lehtinen, O. Evidence against a universal electron-beam-induced virtual temperature in graphene. _Phys. Rev. B_ 2016, _93_ , 134104
  * Holmström et al. 2008 Holmström, E.; Kuronen, A.; Nordlund, K. Threshold defect production in silicon determined by density functional theory molecular dynamics simulations. _Phys. Rev. B_ 2008, _78_ , 045202
  * Susi et al. 2012 Susi, T.; Kotakoski, J.; Arenal, R.; Kurasch, S.; Jiang, H.; Skakalova, V.; Stephan, O.; Krasheninnikov, A. V.; Kauppinen, E. I.; Kaiser, U. et al. Atomistic Description of Electron Beam Damage in Nitrogen-Doped Graphene and Single-Walled Carbon Nanotubes. _ACS Nano_ 2012, _6_ , 8837–8846
  * Enkovaara et al. 2010 Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Dulak, M.; Ferrighi, L.; Gavnholt, J.; Glinsvad, C.; Haikola, V.; Hansen, H. A. et al. Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method. _J. Phys. Condens. Matter_ 2010, _22_ , 253202
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Larsen et al. 2009 Larsen, A. H.; Vanin, M.; Mortensen, J. J.; Thygesen, K. S.; Jacobsen, K. W. Localized atomic basis set in the projector augmented wave method. _Phys. Rev. B_ 2009, _80_ , 195112
  * Monkhorst and Pack 1976 Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. _Phys. Rev. B_ 1976, _13_ , 5188–5192
  * Appelbaum and Hamann 1978 Appelbaum, J. A.; Hamann, D. R. Theory of reconstruction induced subsurface strain — application to Si(100). _Surface Science_ 1978, _74_ , 21–33
  * Kamiyama et al. 2012 Kamiyama, E.; Sueoka, K.; Vanhellemont, J. Ab initio study of vacancy and self-interstitial properties near single crystal silicon surfaces. _Journal of Applied Physics_ 2012, _111_ , 083507, _eprint: https://doi.org/10.1063/1.4703911
  * Swope et al. 1982 Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. _The Journal of Chemical Physics_ 1982, _76_ , 637–649
  * Susi et al. 2016 Susi, T.; Hofer, C.; Argentero, G.; Leuthner, G. T.; Pennycook, T. J.; Mangler, C.; Meyer, J. C.; Kotakoski, J. Isotope analysis in the transmission electron microscope. _Nature Communications_ 2016, _7_ , 13040
  * Schneider and Stoll 1978 Schneider, T.; Stoll, E. Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions. _Physical Review B_ 1978, _17_ , 1302–1322
  * Henkelman et al. 2000 Henkelman, G.; Uberuaga, B.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. _J. Chem. Phys._ 2000, _113_ , 9901
  * Madsen and Susi 2021 Madsen, J.; Susi, T. abTEM: transmission electron microscopy from first principles. _Open Research Europe_ 2021, _1_ , 13015
  * Ophus 2017 Ophus, C. A fast image simulation algorithm for scanning transmission electron microscopy. _Advanced Structural and Chemical Imaging_ 2017, _3_ , 13
  * Höhler et al. 2005 Höhler, H.; Atodiresei, N.; Schroeder, K.; Zeller, R.; Dederichs, P. H. Vacancy complexes with oversized impurities in Si and Ge. _Physical Review B_ 2005, _71_ , 035212
  * Egerton 2013 Egerton, R. Beam-Induced Motion of Adatoms in the Transmission Electron Microscope. _Microscopy and Microanalysis_ 2013, _19_ , 479–486
  * Jamieson et al. 2005 Jamieson, D. N.; Yang, C.; Hopf, T.; Hearne, S. M.; Pakes, C. I.; Prawer, S.; Mitic, M.; Gauja, E.; Andresen, S. E.; Hudson, F. E. et al. Controlled shallow single-ion implantation in silicon using an active substrate for sub-20‐keV ions. _Applied Physics Letters_ 2005, _86_ , 202101
  * Pacheco et al. 2017 Pacheco, J. L.; Singh, M.; Perry, D. L.; Wendt, J. R.; Ten Eyck, G.; Manginell, R. P.; Pluym, T.; Luhman, D. R.; Lilly, M. P.; Carroll, M. S. et al. Ion implantation for deterministic single atom devices. _Review of Scientific Instruments_ 2017, _88_ , 123301