# Nonlocality of the energy density for all single-photon states

M. Federico and H. R. Jauslin<EMAIL_ADDRESS>Laboratoire
Interdisciplinaire Carnot de Bourgogne UMR 6303, CNRS, Université de
Bourgogne, BP 47870, 21078 Dijon, France

###### Abstract

The nonlocality of single-photon states has been analyzed from several
different but interrelared perspectives. In this article, we propose a
demonstration based on the electromagnetic energy density observable and on
the anti-local property of the frequency operator $\Omega=c(-\Delta)^{1/2}$.
The present proof is based on the standard quantization of the electromagnetic
field, which can be formulated equivalently in the momentum representations or
in the position representations of Landau and Peierls [Z. Phys. 62, 188
(1930)] and of Białynicki-Birula [Progress in Optics, edited by E. Wolf
(Elsevier, Amsterdam, 1996)]. Our proof extends to all single-photon states
the results of Białynicki-Birula, which were formulated for two particular
classes of states, those involving a uniform localization [Phys. Rev. Lett.
80, 5247 (1998)] or alternatively states that are electrically or magnetically
localized [Phys.Rev. A 79, 032112 (2009)]. Our approach is formulated in terms
of Knight’s definition of strict localization [J. Math. Phys. 2, 459 (1961)],
based on the comparison of expectation values of single-photon states of local
observables with those of the vacuum.

††preprint: APS/123-QED

##  I Introduction

Localization of photons, together with the existence of a position operator
and consequently of a position wave function for photons, has maintained a
vivid debate in the physics community since the early days of quantum
mechanics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. The
question of the spatial localization of bosons in quantum field theory has
followed several different but interrelated lines.

One line is based on a definition of strictly localized states formulated by
Knight [3]. He showed that for states composed of a finite number of quanta, a
class of space- and time-dependent correlation functions of the Klein-Gordon
field cannot be zero anywhere. Only states involving the superposition of
infinitely many quanta, e.g., coherent states, can be spatially localized.
This result was reformulated and extended by De Bièvre [17, 18] in terms of
expectation values of Weyl unitary operators. In these type of approaches the
main property that leads to nonlocality is the antilocal property of the
frequency operator [19, 20, 21], defined as the unique positive self-adjoint
operator such that its square is minus the Laplacian:
$\Omega=c(-\Delta)^{1/2}$. Antilocality means that if for some square-
integrable field $\vec{v}(\vec{x})$ both $\vec{v}=0$ and $\Omega\vec{v}=0$ in
some finite volume, then $\vec{v}=0$ everywhere in ${\mathbb{R}}^{3}$. In [4,
5], following the results of Knight, Licht characterized the whole set of
strictly localized states by showing that they can be obtained by applying a
partial isometry $\hat{W}^{\dagger}\hat{W}=1$ on the vacuum
$\hat{W}^{\dagger}\ket{\varnothing}$.

Białynicki-Birula showed that for two particular classes of single-photon
states, the expectation value of the energy density cannot be zero in any
finite volume. The two classes considered are states with a spherically
uniform support [22], and states having either a magnetic or an electric
localization property, defined in [23]. The argument in [22] was constructed
using the Paley-Wiener theorem [24]: Single photons with an exponential radial
falloff of the form $\exp(-Ar)$, $A>0$, cannot exist but weaker falloffs are
allowed. An example corresponding to a falloff of the form $\exp(-A\sqrt{r})$
was given. In general, in order to fulfill the constraint from the Paley-
Wiener theorem, a quasiexponential localization is possible with a falloff
$\exp(-Ar^{\gamma})$, where $\gamma<1$. The advantage of the argument of [22]
is that it provides a localization limit and gives a concrete example of a
solution of Maxwell’s equations approaching that limit. In [23] the authors
introduced the notions of electrically and magnetically localized states and
using a proof of the nonlocality of the helicity operator $\Lambda$ (see
Appendix A), they showed that electrically or magnetically localized states
cannot be local if one uses the energy density observable. The interest of the
formulation of localization in terms of the expectation value of the energy
density is that it is an observable that can be measured in current quantum
optics experiments, using e.g. superconducting nanowire detectors [25]. We
show in Appendix B that Knight’s theorem does not imply the results on the
energy density since it does not apply to equal-time correlation functions.
The goal of the present article is to extend the results of Białynicki-Birula
et al. to all single-photon states.

We also mention [26, 27] where it was shown that there exist cylindrical
functions for which a Gaussian falloff is possible in the waist plane only,
making the localization stronger than the exponential limit shown in [22].
More recently in [28, 29, 30], some classes of strictly localized states were
constructed (that are not single-photon states) so that they approach single-
photon states as close as possible.

Following the experimental development of on-demand single-photon sources [31,
32, 33, 34, 35, 36] one can ask how it is possible to produce single photons
in a controlled way without contradicting causality, or how close
experimentally produced photons can be to perfect single-photon states.
Similar questions concerning spontaneous emission have been addressed in [37,
38]. A precise analysis of these experiments can provide a test of whether the
nonlocality predicted by the standard quantization is actually realized in
experiments. As stated, e.g., in [28] the answer can be expected to be that in
fact, those sources produce states that contain an infinity of multi-photon
components (similar but different from coherent sates) that yield local energy
densities, thus avoiding any conflict with causality. The presence of
multiphoton components could in principle be tested in experiments, although
one can expect that the multi-photon components will be quite small and will
depend on the type of source considered.

In a second line, Hegerfeldt [6, 7, 39, 40] established some general
properties of bosons: For states containing exclusively positive energies, if
the state function has a finite support at a given time, the time evolution
will spread it over all space at any later times (see also [41]). In a later
work, concerning causality in Fermi’s two-atom model [40], he showed that for
positive bounded observables $0\leq\hat{\mathcal{O}}\leq 1$, the expectation
value $\bra{\psi(t)}\hat{\mathcal{O}}\ket{\psi(t)}$ is either nonzero for
almost all $t$ or identically zero for all $t$. This opens the question of
compatibility with Einstein causality. Some recent works have explored the
possibility of constructing variations of quantum field theory involving
negative energy states that could avoid Hegerfeldt’s nonlocality [42, 43, 44,
45].

A link between the positivity of energies and the anti-locality of the
frequency operator $\Omega$ can be described as follows: The frequency
operator can be defined as the square root of minus the Laplacian because
$-\Delta$ is a positive operator, which is the reason why the energies are all
positive in the Coulomb gauge quantization of the electromagnetic field.

In the present article we consider the standard quantization of the
electromagnetic field in Fock space with a positive-energy spectrum, which is
most often formulated in a momentum representation [46, 47]. As it was shown
in [48], it can be formulated equivalently in two position representations:
the Landau-Peierls (LP) [1] and the Białynicki-Birula (BB) [11]
representations. The equivalence of these quantizations is formulated in terms
of isomorphisms of the Hilbert spaces with their corresponding scalar products
[48]. An advantage of the BB representation is that it is Lorentz covariant
and independent of a choice of gauge, since it is formulated directly in terms
of the electric and the magnetic fields. All our results are formulated in the
Schrödinger representation but they do not involve at all the time evolution.
The nonlocality of the expectation value of the energy density is proven at
any fixed time, for all single-photon states. It is not created by the time
evolution. This is an important difference from the works of Knight and
Hegerfeldt.

We recall the definition of a localized state that we will use as formulated
by Knight [3]: It is a state which cannot be detected by any means, outside
its volume of localization, i.e., for any local observable, a localized state
will give the same expectation value as the vacuum state, outside its volume
of localization. For example, let $\ket{\varphi}$ be a spatially localized
state in a volume $\mathcal{V}_{\text{s}}$ and
$\hat{O}(\mathcal{V}_{\text{d}})$ a local observable which can probe the state
$\ket{\varphi}$ over a volume $\mathcal{V}_{\text{d}}$. Since $\ket{\varphi}$
is localized, the expectation value of $\hat{O}(\mathcal{V}_{\text{d}})$ can
give the following results [illustrated in FIG. 1 (a) and (b1)]:

  * •

$\bra{\varphi}\hat{O}(\mathcal{V}_{\text{d}})\ket{\varphi}=\bra{\varnothing}\hat{O}(\mathcal{V}_{\text{d}})\ket{\varnothing}$
for $\mathcal{V}_{\text{d}}\cap\mathcal{V}_{\text{s}}=\o$,

  * •

$\bra{\varphi}\hat{O}(\mathcal{V}_{\text{d}})\ket{\varphi}\neq\bra{\varnothing}\hat{O}(\mathcal{V}_{\text{d}})\ket{\varnothing}$
for $\mathcal{V}_{\text{d}}\cap\mathcal{V}_{\text{s}}\neq\o$.

Figure 1: Sketch illustrating Knight’s definition of a localized state. We
consider a state $\ket{\varphi}$ with an associated localization volume
$\mathcal{V}_{\text{s}}$ and a detector $\hat{O}(\mathcal{V}_{\text{d}})$ with
finite volume $\mathcal{V}_{\text{d}}$ which probes the state. Two situations
can occur: (a) a general state where
$\mathcal{V}_{\text{d}}\cap\mathcal{V}_{\text{s}}\neq\o$ and the expectation
value of $\hat{O}(\mathcal{V}_{\text{d}})$ is not equal to that of the vacuum
or (b1) a localized state where
$\mathcal{V}_{\text{d}}\cap\mathcal{V}_{\text{s}}=\o$ and the expectation
value of $\hat{O}(\mathcal{V}_{\text{d}})$ is equal to that of the vacuum for
all observables $\hat{O}$. This means that the state cannot be “seen” outside
its volume of localization by any localized detector. (b2) A state is
considered to be nonlocalized if it can be probed outside its apparent volume
of localization by some observable $\hat{O}(\mathcal{V}_{d})$, i.e.,
$\bra{\varphi}\hat{O}(\mathcal{V}_{\text{d}})\ket{\varphi}\neq\bra{\varnothing}\hat{O}(\mathcal{V}_{\text{d}})\ket{\varnothing}$
for $\mathcal{V}_{\text{d}}\cap\mathcal{V}_{\text{s}}=\o$. In general, this
means that for such extended states, there is no $\mathcal{V}_{\text{s}}$
satisfying (a) and (b1) for any $\hat{O}(\mathcal{V}_{\text{d}})$.

The contraposition implies that a nonlocalized state is a state for which
there exists at least one local observable whose expectation value is not
equal to that of the vacuum at at least one point outside the volume of
localization [FIG 1 (b)].

In this article we will provide a demonstration of the nonlocalization of any
single-photon state, based on the measurement of the mean value of the local
energy density in any finite volume. It is inspired by the pioneering works of
Białynicki-Birula and Białynicki-Birula [22, 23] and the antilocal property of
the frequency operator $\Omega$ [19, 20, 21] (see Appendixes A and D).

Another important operator that will be used in the following is the helicity
operator $\Lambda$ [23, 48], defined by
$c\nabla\times=\Omega\Lambda=\Lambda\Omega$ (see Appendix A). Its spectrum for
transverse fields is $\\{\pm 1\\}$ and it allows us to decompose any
transverse field $\vec{v}$ into its positive and negative helicity parts
$\vec{v}=\vec{v}^{(h+)}+\vec{v}^{(h-)}$ where
$\vec{v}^{(h\pm)}=(1\pm\Lambda)\vec{v}/2$ and
$\Lambda\vec{v}^{(h\pm)}=\pm\vec{v}^{(h\pm)}$.

The article is structured as follows. We first recall how one can describe
single-photon states using position space representations. We then discuss the
notion of a local observable and explain why the energy density is relevant in
the context of local detection of single-photons. We finally formulate our
proof of the nonlocality of photons and illustrate it with some examples.

##  II Representation of a general single-photon state

In the quantum optics literature [46, 47, 49, 50, 48], photons are often
constructed using plane waves, which are not adapted to discuss their
localization in position space. To overcome this issue, one can construct
wavepackets from the general plane wave decomposition [46, 47, 49] or directly
quantize the field using pulses of arbitrary shapes [11, 22, 17, 18, 51, 50,
48, 52]. In this article we will follow the latter choice, and taking
advantage of the result shown in [48], we can switch back and forth between
two equivalent position representations of the quantum theory. The first one,
called the Landau-Peierls (LP) representation [1, 50, 48], is quite close to
the standard plane waves quantization since it also takes the electric field
and the vector potential as canonical variables, but in position space. The
second one is what we call the Białynicki-Birula (BB) representation [11, 22,
51, 48], which takes as canonical variables the electric and magnetic fields.

The LP representation is based on the complex field

$\vec{\psi}(\vec{x})=\sqrt{\frac{\varepsilon_{0}}{2\hbar}}\left[\Omega^{1/2}\vec{A}(\vec{x})-i\Omega^{-1/2}\vec{E}(\vec{x})\right],$
(1)

where $\Omega^{\pm 1/2}$ are self-adjoint operators constructed from the
frequency operator introduced before (see Appendix A), $\vec{A}$ is the vector
potential in the Coulomb gauge and $\vec{x}\in\mathbb{R}^{3}$. The LP field is
an element of the Hilbert space of square-integrable functions from which we
can construct a Fock space of states and annihilation and creation operators
directly on arbitrary pulse-shaped functions [50, 48]

$\displaystyle\hat{B}_{\vec{\psi}}\ket{\varnothing}$ $\displaystyle=0,$ (2a)
$\displaystyle\hat{B}^{\dagger}_{\vec{\psi}}\ket{\varnothing}$
$\displaystyle=\ket*{\vec{\psi}},$ (2b)

satisfying the general bosonic commutation relations [47, 49]

$\left[\hat{B}_{\vec{\psi}},\hat{B}^{\dagger}_{\vec{\psi}^{\prime}}\right]=\bra*{\vec{\psi}}\ket*{\vec{\psi}^{\prime}}_{LP}=\int_{\mathbb{R}^{3}}d^{3}x\
\vec{\psi}^{\star}\cdot\vec{\psi}^{\prime}.$ (3)

The state $\ket*{\vec{\psi}}$ constructed in (2) represents a single-photon
state carried by the classical pulse-shaped field $\vec{\psi}$. Its quantum
dynamics is determined by the classical dynamics of the classical pulse [50,
Appendix B] according to

$\ket*{\vec{\psi}(t)}=\hat{B}^{\dagger}_{\vec{\psi}(t)}\ket{\varnothing},$ (4)

where $\vec{\psi}(t)$ is a solution of the following complex representation of
Maxwell’s equations

$\displaystyle i\frac{\partial\vec{\psi}}{\partial t}$
$\displaystyle=\Omega\vec{\psi},$ (5a) $\displaystyle\nabla\cdot\vec{\psi}(t)$
$\displaystyle=0.$ (5b)

This formulation of the quantum theory is related to the standard quantization
using plane waves through the isomorphism which transforms $\vec{\psi}$ in the
momentum space representation $z(\vec{k},\sigma)$:

$\displaystyle\vec{\psi}(\vec{x})$
$\displaystyle=\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma=\pm}\vec{\phi}_{\vec{k},\sigma}(\vec{x})z(\vec{k},\sigma),$
(6a) $\displaystyle z(\vec{k},\sigma)$
$\displaystyle=\int_{\mathbb{R}^{3}}d^{3}x\
\vec{\phi}_{\vec{k},\sigma}^{\star}(\vec{x})\cdot\vec{\psi}(\vec{x}),$ (6b)

where $\vec{\phi}_{\vec{k},\sigma}$ are the plane waves of wave vector
$\vec{k}$ and polarization $\sigma$ (see Appendix A). Creation and
annihilation operators can thus be developed on the plane-wave basis

$\displaystyle\hat{B}^{\dagger}_{\vec{\psi}}$
$\displaystyle=\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma=\pm}z(\vec{k},\sigma)\hat{B}^{\dagger}_{\vec{\phi}_{\vec{k},\sigma}},$
(7a) $\displaystyle\hat{B}_{\vec{\psi}}$
$\displaystyle=\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma=\pm}z^{\star}(\vec{k},\sigma)\hat{B}_{\vec{\phi}_{\vec{k},\sigma}},$
(7b)

where $\hat{B}^{\dagger}_{\vec{\phi}_{\vec{k},\sigma}}$ and
$\hat{B}_{\vec{\phi}_{\vec{k},\sigma}}$ are the plane-wave creation and
annihilation operators, respectively, in the position representation, i.e.,
the analog of the standard $\hat{a}^{\dagger}_{\vec{k},\sigma}$ and
$\hat{a}_{\vec{k},\sigma}$ in momentum space [46, 47, 49].

The BB representation can be constructed from the LP representation using the
isomorphism $\mathcal{I}$ [48] as

$\vec{F}=\mathcal{I}\vec{\psi}=i\sqrt{\hbar}\Omega^{1/2}\vec{\psi}.$ (8)

Expressed in terms of the real electromagnetic variables, it takes the form

$\vec{F}=\sqrt{\frac{\varepsilon_{0}}{2}}\left(\vec{E}+ic\Lambda\vec{B}\right),$
(9)

where we have used the relation $\Omega=c\Lambda\nabla\times$ (see Appendix A)
and the definition of the vector potential $\nabla\times\vec{A}=\vec{B}$. The
BB vector $\vec{F}$ is related to the Riemann-Silberstein vector
$\vec{F}_{RS}=\sqrt{\varepsilon_{0}}(\vec{E}+ic\vec{B})/\sqrt{2}$ through the
helicity decomposition

$\vec{F}^{(h+)}=\vec{F}_{RS}^{(h+)},\quad\quad\vec{F}^{(h-)}=\left(\vec{F}_{RS}^{(h-)}\right)^{\star}.$
(10)

The construction of creation and annihilation operators for the BB
representation is directly given by the isomorphism $\mathcal{I}$

$\hat{C}_{\vec{F}}=\mathcal{I}\hat{B}_{\vec{\psi}}\mathcal{I}^{-1},\quad\quad\hat{C}_{\vec{F}}^{\dagger}=\mathcal{I}\hat{B}^{\dagger}_{\vec{\psi}}\mathcal{I}^{-1}.$
(11)

They satisfy the commutation relation

$\left[\hat{C}_{\vec{F}},\hat{C}^{\dagger}_{\vec{F}^{\prime}}\right]=\bra*{\vec{F}}\ket*{\vec{F}^{\prime}}_{BB}=\int_{\mathbb{R}^{3}}d^{3}x\
\vec{F}^{\star}\cdot\Omega^{-1}\vec{F}^{\prime}.$ (12)

We have used another letter to refer to the creation and annihilation
operators in the BB representation to emphasize that they do not act on the
same Hilbert space [48], as one can see with the two different scalar products
in the commutation relations (3) and (12), defining two Hilbert spaces

$\displaystyle\mathcal{H}_{LP}$ $\displaystyle=\left\\{\vec{\psi}\
\Big{|}\nabla\cdot\vec{\psi}=0,\bra*{\vec{\psi}}\ket*{\vec{\psi}^{\prime}}_{LP}<\infty\right\\},$
(13a) $\displaystyle\mathcal{H}_{BB}$ $\displaystyle=\left\\{\vec{F}\
\Big{|}\nabla\cdot\vec{F}=0,\bra*{\vec{F}}\ket*{\vec{F}^{\prime}}_{BB}<\infty\right\\}.$
(13b)

The BB scalar product (12) has the advantage that it is Lorentz invariant, and
thus the BB representation is Lorentz covariant. Furthermore, it does not
involve the choice of a gauge since the fields are expressed directly in terms
of the electric and magnetic fields.

##  III Detection model — Local energy observable

We remark that in the quantum field theory of the electromagnetic field the
localization properties cannot be established by looking only at the spatial
properties of the state functions. One has to consider the joint
representation of the states and the local observables. This can be done by
considering e.g. mean values or correlation functions. The fact that the LP
spatial properties do not correspond to the physically measurable properties
had already been stated by Pauli [53]; we will give an explicit example of
this particularity in Section V.

To show the nonlocality of single-photon states, it is enough to find one
particular local observable for which Knight’s localization criterion does not
hold. An operator is said to be a local observable if it represents a physical
measurement which can be made with an instrument well localized in space. Here
$\vec{\hat{E}}(\vec{x})$ and $\vec{\hat{B}}(\vec{x})$ are local since in
practice they can be measured by instruments involving e.g, localized charged
particles or magnetic moments and thus possibly designed as small as required.
Any operator that can be written as a point-wise function of
$\vec{\hat{E}}(\vec{x})$ and $\vec{\hat{B}}(\vec{x})$ is also considered to be
a local observable.

In this section, we introduce the local observable of the energy density that
we will use to show the nonlocality, and which can represent some actual
detectors. It was defined in [22] as

$\hat{\mathcal{E}}_{\text{em}}(\vec{x})=\frac{\varepsilon_{0}}{2}:\left(\vec{\hat{E}}^{2}(\vec{x})+c^{2}\vec{\hat{B}}^{2}(\vec{x})\right):,$
(14)

where $:\cdot:$ stands for the normal ordering.

If one considers an experimental setup close to what is done in [54], the
photons that are produced are “long”, i.e., carried by a pulse with a slowly
varying envelope. Consequently, the associated pulse described in space, is
much bigger than any actual detector. This means that the detector can probe
the photon field only partially, without “seeing” the full state at the same
time: $\mathcal{V}_{\text{s}}\gg\mathcal{V}_{\text{d}}$. Moreover, there exist
efficient single-photon detectors, e.g., superconducting nanowires [25], that
are sensitive to the electromagnetic energy. Those detectors are prepared at a
temperature that is slightly below the critical tempeature $T_{\text{c}}$ of
the superconducting nanowire, so it has no resistance. When a photon triggers
the detector, there is a local absoption of energy, heating up the detector
above $T_{\text{c}}$ and yielding a measurable resistance that signals the
detection of a photon.

What is measured in such experiments, can thus be modeled by the local energy
density $\hat{\mathcal{E}}_{\text{em}}(\vec{x})$ integrated over the volume of
the detector $\mathcal{V}_{\text{d}}$ i.e.

$\hat{\mathcal{E}}_{\mathcal{V}_{\text{d}}}=\int_{\mathcal{V}_{\text{d}}}d^{3}x\
\hat{\mathcal{E}}_{\text{em}}(\vec{x}).$ (15)

The expectation value of the energy density operator for a general single-
photon state in the LP representation,
$\ket*{\text{1ph}}=\hat{B}_{\vec{\psi}}^{\dagger}\ket{\varnothing}$ for any
$\vec{\psi}\in\mathcal{H}_{LP}$ or equivalently written [48] in the BB
representation
$\ket*{\text{1ph}}=\hat{C}_{\vec{F}}^{\dagger}\ket{\varnothing}$ for
$\vec{F}=i\sqrt{\hbar}\Omega^{1/2}\vec{\psi}\in\mathcal{H}_{BB}$, was computed
in [22]. It can be written as (see Appendix C for an indication of the
calculation)

$\displaystyle\langle\hat{\mathcal{E}}_{\text{em}}(\vec{x})\rangle_{\ket*{\text{1ph}}}$
$\displaystyle=\hbar\absolutevalue*{\Omega^{1/2}\vec{\psi}^{(h+)}(\vec{x})}^{2}+\hbar\absolutevalue*{\Omega^{1/2}\vec{\psi}^{(h-)}(\vec{x})}^{2}$
(16a)
$\displaystyle=\absolutevalue*{\vec{F}^{(h+)}(\vec{x})}^{2}+\absolutevalue*{\vec{F}^{(h-)}(\vec{x})}^{2}\geq
0.$ (16b)

We will show in the next section that this result implies the nonlocality of
single photons.

##  IV Proof of the nonlocality of the energy density for all single-photon
states

The result obtained above for the mean value of the energy density operator is
clearly greater than or equal to zero for any single-photon pulse, i.e., for
any function $\vec{\psi}$ or $\vec{F}$ representing the single-photon state.
In this section we will show that the average local energy is strictly
different from zero at any position in space

$\langle\hat{\mathcal{E}}_{\text{em}}(\vec{x})\rangle_{\ket*{\text{1ph}}}\neq
0\quad\forall\vec{x}\in\mathbb{R}^{3}.$ (17)

This result is a consequence of the following lemmas:

Lemma 1: For any field $\vec{v}(\vec{x})$ that is not identically zero,
$\Omega\vec{v}$ and $\vec{v}$ cannot be both zero in any open set of
$\mathbb{R}^{3}$ [19, 20, 21].

Lemma 2: Fields of $\pm$ helicity, i.e.
$\Lambda\vec{g}^{(h\lambda)}(\vec{x})=\lambda\vec{g}^{(h\lambda)}(\vec{x})$,
$\lambda=\pm$, have the property that either $\vec{g}^{(h\lambda)}$ is
identically zero or $\vec{g}^{(h\lambda)}\neq 0$ in any open set of
$\mathbb{R}^{3}$.

A proof of Lemma 1 is given in Appendix D and Lemma 2 can be shown directly as
follows: Let $\vec{g}$ be a transverse field,
$\vec{g}=\vec{g}^{(h+)}+\vec{g}^{(h-)}$. The helicity components are
eigenfunctions of the helicity operator
$\Lambda\vec{g}^{(h\pm)}=\pm\vec{g}^{(h\pm)}$. Using the definition of the
helicity operator $\Lambda=c\Omega^{-1}\nabla\times$, we can reformulate this
relation as

$c\nabla\times\vec{g}^{(h\pm)}=\pm\Omega\vec{g}^{(h\pm)}.$ (18)

Taking a given $\lambda=\pm$, if $\vec{g}^{(h\lambda)}$ is zero in an open set
$\mathcal{S}$, then $\nabla\times\vec{g}^{(h\lambda)}$ is zero in the same set
and finally (18) implies that $\Omega\vec{g}^{(h\lambda)}$ is also zero in
$\mathcal{S}$, which by Lemma 1 implies that $\vec{g}^{(h\lambda)}$ is zero
everywhere. We conclude that $\vec{g}^{(h\lambda)}$ is either identically zero
or non-zero for any open set as stated in Lemma 2.

The central result (17) can thus be shown by taking $\vec{g}$ to be the BB
representation $\vec{F}$ of any single-photon state, which is transverse, and
thus can be decomposed as $\vec{F}=\vec{F}^{(h+)}+\vec{F}^{(h-)}$. We conclude
then that $\vec{F}^{(h\pm)}$ is either identically zero or nonzero in any open
set. Moreover, since (16b) has two terms, even if one of them is identically
zero, the other cannot be zero too since it would mean that the single-photon
state itself is zero. This result is valid for any open set and therefore one
can extend it to any point $\vec{x}\in\mathbb{R}^{3}$, which completes the
proof.

Thus, since the zero point energy has been removed using the normal ordering
in (14),
$\langle\hat{\mathcal{E}}_{\text{em}}(\vec{x})\rangle_{\ket*{1\text{ph}}}$ is
never equal to the vacuum mean value, preventing Knight’s localization
criterion to be fulfilled for any single-photon state. In physical terms this
means that if the electromagnetic field is prepared in a single-photon state,
a detector, placed anywhere in space, which measures the energy in a finite
volume, has a non-zero probability of detecting the photon. The probability
can be small, but it is strictly non-zero.

##  V Illustration of the nonlocality

Figure 2: Illustration of the nonlocality of single-photon states for two
extreme cases. (a) The blue solid line shows the LP representation of a
single-photon state $\ket*{\psi^{LP}_{\text{comp}}}$ with compact support and
the green dotted line its BB representation. The red dashed line shows the
expectation value of the energy density computed for that state. The compact
support property of $\ket*{\psi^{LP}_{\text{comp}}}$ is lost for both the BB
representation and the energy. (b) The blue solid line shows the LP
representation of an extended single-photon state $\ket{\psi_{\text{ext}}}$
and the dotted green line its BB representation. The red dashed line shows the
expectation value of the energy density computed from that state. The extended
property of $\ket{\psi_{\text{ext}}}$ is visible from both representation and
for the energy. (c) The green dotted line shows the BB representation of a
single-photon state $\ket{F^{BB}_{\text{comp}}}$ with compact support and the
solid blue line its LP representation. The red dashed line shows the
expectation value of the energy density computed for that state. The compact
support property of $\ket{F^{BB}_{\text{comp}}}$ is lost for the LP
representation and the energy. (d), (e), and (f) are the same plots as (a),
(b), and (c), respectively, but with a logarithmic scale.

The nonlocality brought by the splitting into helicity components can be
illustrated through simple one-dimensional examples. We can compute the
expectation value
$\langle\hat{\mathcal{E}}_{\text{em}}(\vec{x})\rangle_{\ket{\text{1ph}}}$ for
single-photon states representing three extreme cases: First we consider a
state
$\ket*{\psi^{LP}_{\text{comp}}}=\hat{B}_{\psi^{LP}_{\text{comp}}}^{\dagger}\ket{\varnothing}$,
where $\psi^{LP}_{\text{comp}}\in\mathcal{H}_{LP}$ is a function of compact
support i.e. $\psi^{LP}_{\text{comp}}(x)=0$ outside an interval of size $L$;
then we consider a state
$\ket*{F^{BB}_{\text{comp}}}=\hat{B}_{F^{BB}_{\text{comp}}}^{\dagger}\ket{\varnothing}$,
where $F^{BB}_{\text{comp}}\in\mathcal{H}_{BB}$ is a function of compact
support, i.e., $F^{BB}_{\text{comp}}(x)=0$ outside an interval of size $L$;
and finally we consider a state
$\ket*{\psi_{\text{ext}}}=\hat{B}_{\psi_{\text{ext}}}^{\dagger}\ket{\varnothing}$,
where $\psi_{\text{ext}}\in\mathcal{H}_{LP}$ is extended over all space, i.e.,
$\psi_{\text{ext}}(x)\neq 0$ for any $x\in\mathbb{R}$. To construct
$\psi_{\text{ext}}$, we use real fields $E_{\text{comp}}(x)$ and
$A_{\text{comp}}(x)$ with support in the interval $[-L/2,L/2]$ of the form

$\displaystyle E_{\text{comp}}(x)$
$\displaystyle\propto\left\\{\begin{array}[]{ll}\sin^{2}(\frac{\pi}{L}x+\frac{\pi}{2})\quad\text{if
$x\in[-\frac{L}{2},\frac{L}{2}]$}\\\ 0\hskip
62.59596pt\text{otherwise},\end{array}\right.$ (19c) $\displaystyle
A_{\text{comp}}(x)$
$\displaystyle\propto\left\\{\begin{array}[]{ll}\sin^{2}(\frac{\pi}{L}x+\frac{\pi}{2})\quad\text{if
$x\in[-\frac{L}{2},\frac{L}{2}]$}\\\ 0\hskip
62.59596pt\text{otherwise},\end{array}\right.$ (19f)

and to build $\psi^{LP}_{\text{comp}}$ we take the extended fields

$\displaystyle E^{LP}_{\text{ext}}(x)$
$\displaystyle=\Omega^{1/2}E_{\text{comp}}(x),$ (20a) $\displaystyle
A^{LP}_{\text{ext}}(x)$ $\displaystyle=\Omega^{-1/2}A_{\text{comp}}(x).$ (20b)

The resulting $\psi^{LP}_{\text{comp}}$ and $\psi_{\text{ext}}$ are
represented as the solid blue lines in FIG. 2(a) and 2(d) and FIG. 2(b) and
FIG. 2(e), respectively. To construct a localized BB representation
$F^{BB}_{\text{comp}}$ we take the fields

$\displaystyle E^{BB}(x)$ $\displaystyle=E_{\text{comp}}(x),$ (21a)
$\displaystyle A^{BB}(x)$ $\displaystyle=\Omega^{-1}A_{\text{comp}}(x).$ (21b)

The resulting localized BB representation is shown by the green dotted line in
FIG. 2(c) and 2(f).

For these three examples, we compute the expectation value of the energy
density operator and obtain the results displayed as the red dashed lines in
FIG. 2. In general, one can see that the localization property of the LP or BB
representation does not give any information for the localization of the
energy density. Indeed, the compact support property of the states
$\psi^{LP}_{\text{comp}}$ and $F^{BB}_{\text{comp}}$ is not preserved for the
expectation value of the local energy [FIG 2(a), 2(d), and FIG. 2(c), and
2(f)], as expected. Moreover, a localized LP representation implies a
nonlocalized BB representation and vice versa due to the form of the
isomorphism $\mathcal{I}$. This is illustrated in FIG 2(a), and 2(d) and FIG.
2(c), and 2(f) where we have either a localized LP function and a nonlocalized
BB function or a nonlocalized LP function and a localized BB function. The
most general case is shown in FIG. 2(b) and 2(e), where neither the LP nor the
BB representation is localized and so neither is the energy density.

##  VI Concluding remarks

We point out that none of the statements in this article involve time
evolution. They are statements about the expectation value of the energy
density for any single photon state at any fixed time. This is an important
difference from the works of Knight and of Hegerfeldt. The nonlocality, as
probed by the energy density, is not created by the time evolution; it is
already present with the initial condition for any state as shown in FIG. 2.
Our point of view is that the nonlocality appears through the frequency
operator $\Omega$ and it manifests in several ways: in time-independent
expectation values like the energy density; in properties of the time
evolution, as described by Hegerfeldt; or in time-dependent correlation
functions, as treated by Knight. The time evolution of the states is
determined by the frequency operator $\Omega$, both in the Landau-Peierls and
in the Białynicki-Birula representations [48]:

$i\frac{\partial\vec{\psi}}{\partial t}=\Omega\vec{\psi},\quad\quad\quad
i\frac{\partial\vec{F}}{\partial t}=\Omega\vec{F}.$ (22)

The antilocality and the positivity of $\Omega$ are at the origin of the time-
independent and the time-dependent nonlocality properties of single-photon
states. This different aspects are complementary and they can contribute to a
better understanding of the nonlocality of single-photon states.

## Acknowledgments

We thank Jonas Lampart for many fruitful discussions. This work was supported
by the “Investissements d’Avenir” program, project ISITE-BFC/IQUINS
(ANR-15-IDEX-03), QUACO-PRC (ANR-17-CE40-0007-01) and the EUR-EIPHI Graduate
School (17-EURE-0002). We also acknowledge support from the European Union’s
Horizon 2020 research and innovation program under the Marie Sklodowska-Curie
Grant Agreement No. 765075 (LIMQUET).

##  Appendix A Plane wave basis, curl, frequency and helicity operators

We consider the generalized basis of transverse plane waves
$\\{\vec{\phi}_{\vec{k},\sigma}\\}$ defined by

$\vec{\phi}_{\vec{k},\sigma}(\vec{x})=(2\pi)^{-3/2}\vec{\epsilon}_{\sigma}(\vec{k})e^{i\vec{k}\cdot\vec{x}},$
(23)

where $\sigma=\pm$, and the circular polarization vectors can be chosen as

$\vec{\epsilon}_{+}(\vec{k})=\frac{1}{\sqrt{2}|\vec{k}|\sqrt{k_{x}^{2}+k_{y}^{2}}}\left[\begin{array}[]{c}-k_{x}k_{z}+i|\vec{k}|k_{y}\\\
-k_{y}k_{z}-i|\vec{k}|k_{x}\\\ k_{x}^{2}+k_{y}^{2}\end{array}\right],$ (24)

with $\vec{\epsilon}_{-}(\vec{k})=\vec{\epsilon}_{+}(\vec{k})^{\star}$. They
are eigenfunctions of the curl operator with eigenvalues
$\nabla\times\vec{\phi}_{\vec{k},\sigma}=\sigma|\vec{k}|\,\vec{\phi}_{\vec{k},\sigma}$.
To construct the frequency operator $\Omega$, we recall that it is defined as
the unique positive operator satisfying $\Omega^{2}=-c^{2}\Delta$ where
$\Delta$ is the Laplace operator which can be written for transverse fields as
$-\Delta=\nabla\times\nabla\times$. Therefore, the frequency operator
satisfies
$\Omega^{2}\vec{\phi}_{\vec{k},\sigma}=c^{2}\absolutevalue*{\vec{k}}^{2}\vec{\phi}_{\vec{k},\sigma}=\omega_{\vec{k}}^{2}\vec{\phi}_{\vec{k},\sigma}$.
The positive square root of $\Omega^{2}$ can thus be defined by its action on
the continuum basis of plane waves

$\Omega\vec{\phi}_{\vec{k},\sigma}=\omega_{\vec{k}}\vec{\phi}_{\vec{k},\sigma},\>\>\>\>\>\omega_{\vec{k}}>0,$
(25)

and the particular powers $\Omega^{\pm 1/2}$ used in the LP representation (1)
can similarly be defined as $\Omega^{\pm
1/2}\vec{\phi}_{\vec{k},\sigma}=\omega_{\vec{k}}^{\pm
1/2}\vec{\phi}_{\vec{k},\sigma}$. Here $\Omega^{2}$, $\Omega$ and $\Omega^{\pm
1/2}$ are all positive selfadjoint operators.

The helicity operator $\Lambda$ is then defined through a combination of the
curl and the frequency operators by

$c\nabla\times=\Omega\Lambda.$ (26)

It has the same eigenfunctions with eigenvalues
$\Lambda\vec{\phi}_{\vec{k},\sigma}=\sigma\vec{\phi}_{\vec{k},\sigma}$ and
therefore commutes with both $\nabla\times$ and $\Omega$. One can decompose
any transverse field $\vec{v}$ into a sum of a positive and a negative
helicity part $\vec{v}^{(h\pm)}$

$\vec{v}=\vec{v}^{(h+)}+\vec{v}^{(h-)},$ (27)

where $\Lambda\vec{v}^{(h\pm)}=\pm\vec{v}^{(h\pm)}$. The positive and negative
helicity parts can be constructed by applying the projectors
$\mathbb{P}^{(h\pm)}=(1\pm\Lambda)/2$, i.e.
$\vec{v}^{(h\pm)}=\mathbb{P}^{(h\pm)}\vec{v}$. We also remark that
$\Lambda^{2}=1$ and $\Lambda^{-1}=\Lambda$.

Helicity can be interpreted as the projection of the spin on the direction of
motion [48].

##  Appendix B A counterexample to the extension of Knigth’s construction for
equal-time correlation functions

In this appendix we show that Knight’s result on time-dependent correlation
functions cannot be extended to equal-time correlation functions, nor in
particular to the expectation value of the energy density. Since we are
interested in photons, we adapt Knight’s construction for the scalar Klein-
Gordon equation to the electromagnetic case but a similar argument can be made
for massive scalar fields. In fact the argument is based on the construction
that was made by Białynicki-Birula and Białynicki-Birula in [23] in order to
show that there exist single-photon states with either electric or magnetic
localization (but not both).

We start with a state
$\ket*{\vec{\psi}}=\hat{B}^{\dagger}_{\vec{\psi}}\ket{\varnothing}$ with the
properties

$\displaystyle\left[\hat{B}_{\vec{\psi}},\hat{B}^{\dagger}_{\vec{\psi}}\right]$
$\displaystyle=1,$ (28a)
$\displaystyle\left[\vec{\hat{A}}(\vec{x}_{j},t_{0}),\hat{B}^{\dagger}_{\vec{\psi}}\right]$
$\displaystyle=0\quad\forall\vec{x}_{j}\notin\mathcal{R},$ (28b)

where $\mathcal{R}$ is an open set of $\mathbb{R}^{3}$ and $t_{0}$ a given
time which we take as zero in the following without lost of generality. The
mean value of the equal-time correlation function of the potential vector
outside $\mathcal{R}$ for a state with such properties is thus

$\displaystyle\bra*{\vec{\psi}}\vec{\hat{A}}(\vec{x}_{1})\dots\vec{\hat{A}}(\vec{x}_{M})\ket*{\vec{\psi}}=\bra*{\varnothing}\hat{B}_{\vec{\psi}}\vec{\hat{A}}(\vec{x}_{1})\dots\vec{\hat{A}}(\vec{x}_{M})\hat{B}_{\vec{\psi}}^{\dagger}\ket*{\varnothing}$
$\displaystyle\quad\quad\quad=\bra*{\varnothing}\hat{B}_{\vec{\psi}}\hat{B}_{\vec{\psi}}^{\dagger}\vec{\hat{A}}(\vec{x}_{1})\dots\vec{\hat{A}}(\vec{x}_{M})\ket*{\varnothing}$
$\displaystyle\quad\quad\quad=\bra*{\varnothing}\left(1+\hat{B}^{\dagger}_{\vec{\psi}}\hat{B}_{\vec{\psi}}\right)\vec{\hat{A}}(\vec{x}_{1})\dots\vec{\hat{A}}(\vec{x}_{M})\ket*{\varnothing}$
$\displaystyle\quad\quad\quad=\bra*{\varnothing}\vec{\hat{A}}(\vec{x}_{1})\dots\vec{\hat{A}}(\vec{x}_{M})\ket*{\varnothing}.$
(29)

This calculation shows that outside $\mathcal{R}$, the mean value is that of
the vacuum. The key point of this result is to construct a state
$\ket*{\vec{\psi}}$ with the properties (28). To do so, the commutator (28b)
can be written as

$\displaystyle\left[\vec{\hat{A}}(\vec{x}_{j}),\hat{B}^{\dagger}_{\vec{\psi}}\right]$
$\displaystyle=\sqrt{\frac{\hbar}{2\varepsilon_{0}}}\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma=\pm}\omega_{\vec{k}}^{-1/2}\vec{\phi}_{\vec{k},\sigma}(\vec{x}_{j})\left[\hat{B}_{\vec{\phi}_{\vec{k},\sigma}},\hat{B}^{\dagger}_{\vec{\psi}}\right]$
$\displaystyle=\sqrt{\frac{\hbar}{2\varepsilon_{0}}}\Omega^{-1/2}\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma=\pm}\vec{\phi}_{\vec{k},\sigma}(\vec{x}_{j})\bra*{\vec{\phi}_{\vec{k},\sigma}}\ket*{\vec{\psi}}$
$\displaystyle=\sqrt{\frac{\hbar}{2\varepsilon_{0}}}\Omega^{-1/2}\vec{\psi}(\vec{x}_{j}),$
(30)

which is zero outside $\mathcal{R}$ if one takes a function $\vec{\xi}$ with
compact support in $\mathcal{R}$ i.e. $\vec{\xi}(\vec{x}_{j})=0$ for all
$\vec{x}_{j}\notin\mathcal{R}$ and defines $\vec{\psi}=\Omega^{1/2}\vec{\xi}$.

##  Appendix C Derivation of the mean value of the energy density for a
general single-photon state

In order to have a self-contained proof of nonlocality discussed in this
article, we provide here an indication of the main steps for calculating the
mean value of the energy density given in (16). For simplicity, we will show
here only the derivation done with the BB representation of the states, but
one can equivalently compute it with the LP representation. The result can
anyway be expressed easily in both representations using the isomorphism
$\mathcal{I}$ [see eq. (8)].

We first remark that the normal-ordered energy density observable (14) can be
written as [22]

$\hat{\mathcal{E}}_{\text{em}}=\frac{\varepsilon_{0}}{2}:\left(\vec{\hat{E}}-ic\vec{\hat{B}}\right)\cdot\left(\vec{\hat{E}}+ic\vec{\hat{B}}\right):,$
(31)

and that electromagnetic field operators can be expressed as

$\displaystyle\vec{\hat{E}}$
$\displaystyle=\frac{1}{\sqrt{2\varepsilon_{0}}}\left(\vec{\hat{{\bf
F}}}+\vec{\hat{{\bf F}}}^{\dagger}\right),$ (32a) $\displaystyle\vec{\hat{B}}$
$\displaystyle=\frac{-i}{\sqrt{2\varepsilon_{0}c^{2}}}\Lambda\left(\vec{\hat{{\bf
F}}}-\vec{\hat{{\bf F}}}^{\dagger}\right),$ (32b)

in terms of the BB field operator $\vec{\hat{{\bf F}}}$ defined by [48]

$\vec{\hat{{\bf
F}}}(\vec{x})=\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma=\pm}\vec{g}_{\vec{k},\sigma}(\vec{x})\hat{C}_{\vec{g}_{\vec{k},\sigma}}.$
(33)

The set of functions $\\{\vec{g}_{\vec{k},\sigma}\\}$ introduced here is a
basis of the BB Hilbert space $\mathcal{H}_{BB}$ (see eq. (13b)). By inserting
equations (32) into (31) and using the decomposition $\vec{\hat{{\bf
F}}}=\vec{\hat{{\bf F}}}^{(h+)}+\vec{\hat{{\bf F}}}^{(h-)}$ which can be done
with the projectors $\mathbb{P}^{(h\pm)}=(1\pm\Lambda)/2$, we can rewrite the
energy density as

$\hat{\mathcal{E}}_{\text{em}}(\vec{x})=\ :\left(\vec{\hat{{\bf
F}}}^{(h+){\dagger}}+\vec{\hat{{\bf
F}}}^{(h-)}\right)\cdot\left(\vec{\hat{{\bf F}}}^{(h+)}+\vec{\hat{{\bf
F}}}^{(h-){\dagger}}\right):.$ (34)

Among the four terms of (34), only two will give non zero expectation values
and thus we have

$\langle\hat{\mathcal{E}}_{\text{em}}\rangle_{\ket*{1\text{ph}}}=\sum_{\lambda=\pm}\bra{\varnothing}\hat{C}_{\vec{F}}\vec{\hat{{\bf
F}}}^{(h\lambda){\dagger}}\cdot\vec{\hat{{\bf
F}}}^{(h\lambda)}\hat{C}_{\vec{F}}^{\dagger}\ket{\varnothing}.$ (35)

To work out the final result, we use the commutators

$\left[\vec{\hat{{\bf
F}}}^{(h\pm)}(\vec{x}),\hat{C}_{\vec{F}}^{\dagger}\right]=\vec{F}^{(h\pm)}(\vec{x}),$
(36)

and obtain (16) in the BB representation

$\langle\hat{\mathcal{E}}_{\text{em}}\rangle_{\ket*{1\text{ph}}}=\absolutevalue*{\vec{F}^{(h+)}}^{2}+\absolutevalue*{\vec{F}^{(h-)}}^{2}.$
(37)

##  Appendix D Anti-locality of the frequency operator

We provide in this appendix a proof of Lemma 2 used to show the nonlocality of
photons in Section IV. This result was shown in [19, 20, 21] and we will
sketch here the argument of [20]. We recall the statement.

Lemma 2: For any field $\vec{v}(\vec{x})$ that is not identically zero,
$\Omega\vec{v}$ and $\vec{v}$ cannot be both zero in any open set of
$\mathbb{R}^{3}$.

We are going to prove that if $\vec{v}$ and $\Omega\vec{v}$ are both equal to
zero in some open set $\mathcal{S}$, it imples that $\vec{v}(\vec{x})=0$
everywhere.

Proof. Since $\Omega$ is positive and selfadjoint, the operators
$U(t)=\exp(i\Omega t)$, $t\in(-\infty,+\infty)$ define a one-parameter family
of unitary operators. The field defined as
$\vec{u}(\vec{x},t)=U(t)\vec{v}(\vec{x})$ satisfies the wave equation

$\frac{\partial^{2}\vec{u}}{\partial t^{2}}=-\Omega^{2}\vec{u},$ (38)

with initial conditions

$\displaystyle\vec{u}(\vec{x},t=0)$ $\displaystyle=\vec{v}(\vec{x}),$ (39a)
$\displaystyle\frac{\partial\vec{u}}{\partial t}(\vec{x},t=0)$
$\displaystyle=i\Omega\vec{v}(\vec{x}).$ (39b)

Since the solutions of the wave equation propagate with a finite speed $c$,
the property of the initial conditions to be zero i.e. $\vec{v}(\vec{x})=0$ in
a set $\mathcal{S}$, implies that there is a $t_{0}>0$ and a non-empty open
subset $\mathcal{S}_{0}\subset\mathcal{S}$ such that $\vec{u}(\vec{x},t)=0$
for all $\vec{x}\in\mathcal{S}_{0}$ and $0\leq t<t_{0}$. Thus, for any
$\mathcal{C}^{\infty}$ field $\vec{\varphi}(\vec{x})$ with compact support in
$\mathcal{S}_{0}$

$\bra*{\vec{\varphi}}\ket*{\vec{u}(\cdot,t)}=\int_{\mathbb{R}^{3}}d^{3}x\
\vec{\varphi}(\vec{x})^{\star}\cdot\vec{u}(\vec{x},t)=0$ (40)

for $0\leq t<t_{0}$. We now consider the continuation of the variable $t$ into
the upper complex half plane and define the function

$f(z)=\bra*{\vec{\varphi}}\ket*{e^{i\Omega z}\vec{v}},$ (41)

for $\imaginary z\geq 0$ that has the following properties [20]:

  1. (i)

$f(z)$ is holomorphic for $\imaginary z>0$, and continuous for $\imaginary
z\geq 0$,

  2. (ii)

$f(t)=\bra*{\vec{\varphi}}\ket*{\vec{u}(\cdot,t)}$ when $t\in\mathbb{R}$, and
$f(t)\in\mathbb{R}$ for $t\in(0,t_{0})$,

  3. (iii)

$f(t)=0$ for $0\leq t<t_{0}$.

We remark that for $t>t_{0}$, $f(t)$ is not necessarily zero or real. The main
steps of a proof that $f(z)$ is holomorphic in the upper half-plane, by
showing the existence of the derivative $df/dz$, can be summarized as follows:
$\vec{v}$ can be developed in the basis $\vec{\phi}_{\vec{k},\sigma}$ of
continuum eigenfunctions of the Laplacian
$\vec{v}=\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma}\alpha_{\vec{k},\sigma}\vec{\phi}_{\vec{k},\sigma}$.
We can then write

$\displaystyle\frac{d}{dz}f(z)$
$\displaystyle=\frac{d}{dz}\int_{\mathcal{S}_{0}}d^{3}x\
\vec{\varphi}^{\star}\cdot e^{i\Omega z}\vec{v}$
$\displaystyle=\int_{\mathcal{S}_{0}}d^{3}x\ \vec{\varphi}^{\star}\cdot
e^{i\Omega z}i\Omega\vec{v}$
$\displaystyle=\int_{\mathcal{S}_{0}}d^{3}x\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma}\vec{\varphi}^{\star}\cdot\vec{\phi}_{\vec{k},\sigma}\alpha_{\vec{k},\sigma}i\omega_{\vec{k}}e^{i\omega_{\vec{k}}z}.$
(42)

The derivative can be brought inside the integral by the Weierstrass M-test
and the Lebesgue dominated convergence theorem, because the integral converges
absolutely since $e^{i\Omega z}$ is a contractive semi-group [55]. The
existence of the derivative is then obtained by exchanging the two integrals
by the Fubini-Tonelli theorem: Writing $z$ as $z=z_{r}+iz_{i}$, one obtains

$\displaystyle\absolutevalue{\frac{d}{dz}f(z)}$
$\displaystyle\leq\int_{\mathcal{S}_{0}}d^{3}x\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma}\absolutevalue*{\vec{\varphi}^{\star}\cdot\vec{\phi}_{\vec{k},\sigma}}\absolutevalue*{\alpha_{\vec{k},\sigma}}\omega_{\vec{k}}e^{-\omega_{\vec{k}}z_{i}}$
$\displaystyle=\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma}\absolutevalue*{\alpha_{\vec{k},\sigma}}\omega_{\vec{k}}e^{-\omega_{\vec{k}}z_{i}}\int_{\mathcal{S}_{0}}d^{3}x\absolutevalue*{\vec{\varphi}^{\star}\cdot\vec{\phi}_{\vec{k},\sigma}}$
$\displaystyle\leq
C\int_{\mathbb{R}^{3}}d^{3}k\sum_{\sigma}\absolutevalue*{\alpha_{\vec{k},\sigma}}\omega_{\vec{k}}e^{-\omega_{\vec{k}}z_{i}}<\infty,$
(43)

where $C$ is a constant. The last integral is finite for $z_{i}>0$ which
completes the argument for the existence of $df/dz$.

We will now use the Schwarz reflection principle [56, p.75], which in the
present context can be formulated as follows: If $f_{+}(z)$ satisfies the
properties that

  1. (S-i)

$f_{+}(z)$ is holomorphic in the open upper complex rectangle
$D_{+}=\left\\{\imaginary z>0,\real z\in(0,t_{0})\right\\}$,

  2. (S-ii)

$f_{+}(z)$ continuous in $D_{+}\cup(0,t_{0})$,

  3. (S-iii)

$f_{+}(t)$ is real in the interval $t\in(0,t_{0})$,

then $f_{+}(z)$ can be continued holomorphically trough the interval
$(0,t_{0})$ to the lower rectangle $D_{-}=\left\\{\imaginary z\leq 0,\real
z\in(0,t_{0})\right\\}$, by defining

$\displaystyle
f_{+}(z)=\left\\{\begin{array}[]{ll}f_{+}(z)\>\>\>\>\>\>\>\>\>\>\text{for
$z\in D_{+}$},\\\ (f_{+}(z^{\star}))^{\star}\>\>\text{for $z\in
D_{-}$}.\end{array}\right.$ (46)

This defines thus a holomorphic function $f_{+}(z)$ in the whole open set
$D_{+}\cup D_{-}=\left\\{\real z\in(0,t_{0})\right\\}$, which includes the
interval $(0,t_{0})$.

By applying the Schwarz reflection principle to the function $f(z)$ that we
combine with the properties (i) and (ii), it shows that $f(z)$ is analytic in
the union of the open upper half-plane and $D_{-}$. Property (iii) states that
$f(z)=0$ in the interval $z\in(0,t_{0})$, which implies that $f(z)=0$ in the
whole region where $f$ is holomorphic, in particular in the whole open upper
half-plane $\imaginary z>0$. Since according to (i), $f(z)$ is continuous for
$\imaginary z\geq 0$, this implies that $f(z)=0$ for $\imaginary z\geq 0$. In
particular, $f(t):\bra*{\vec{\varphi}}\ket*{\vec{u}(\cdot,t)}=0$ for
$-\infty<t<\infty$. Since $\vec{\varphi}$ is an arbitrary function, this
implies that $\vec{u}(\vec{x},t)=0$ for $-\infty<t<\infty$ and
$\vec{x}\in\mathcal{S}_{0}$. The unique continuation theorem for solutions
$\vec{u}(\vec{x},t)$ of the wave equation, proven e.g. in [57], states that if
$\vec{u}(\vec{x},t)=0$ for an open set for all $-\infty<t<\infty$, then
$\vec{u}(\vec{x},t)=0$ everywhere. In particular, for $t=0$ this implies that

$\vec{u}(\vec{x},t=0)=\vec{v}(\vec{x})=0,\forall\vec{x}\in\mathbb{R}^{3},$
(47)

which completes the proof.

## References

  * Landau and Peierls [1930] L. Landau and R. Peierls, Zeit. F. Phys 62, 188 (1930).
  * Newton and Wigner [1949] T. D. Newton and E. P. Wigner, Reviews of Modern Physics 21, 400 (1949).
  * Knight [1961] J. M. Knight, Journal of Mathematical Physics 2, 459 (1961).
  * Licht [1963] A. L. Licht, Journal of Mathematical Physics 4, 1443 (1963).
  * Licht [1966] A. L. Licht, Journal of Mathematical Physics 7, 1656 (1966).
  * Hegerfeldt [1974] G. C. Hegerfeldt, Physical Review D 10, 3320 (1974).
  * Hegerfeldt and Ruijsenaars [1980] G. C. Hegerfeldt and S. N. M. Ruijsenaars, Physical Review D 22, 377 (1980).
  * Cook [1982] R. J. Cook, Physical Review A 25, 2164 (1982).
  * Białynicki-Birula [1994] I. Białynicki-Birula, Acta Physica Polonica A 86, 97 (1994).
  * Sipe [1995] J. E. Sipe, Physical Review A 52, 1875 (1995).
  * Białynicki-Birula [1996] I. Białynicki-Birula, in  _Progress in Optics_, edited by E. Wolf (Elsevier, 1996).
  * Hawton [1999] M. Hawton, Physical Review A 59, 3223 (1999).
  * Keller [2005] O. Keller, Physics Reports 411, 1 (2005).
  * Smith and Raymer [2007] B. J. Smith and M. G. Raymer, New Journal of Physics 9, 414 (2007).
  * Keller [2016] O. Keller,  _Light - The Physics of the Photon_ (CRC Press, 2016).
  * Kiessling and Tahvildar-Zadeh [2018] M. K. H. Kiessling and A. S. Tahvildar-Zadeh, Journal of Mathematical Physics 59, 10.1063/1.5021066 (2018).
  * De Bièvre [2006] S. De Bièvre, in  _Large Coulomb Systems_ (Springer Berlin Heidelberg, 2006) pp. 15–61.
  * De Bièvre [2007] S. De Bièvre, in  _Contributions in Mathematical Physics_ (Hindustan Book Agency, 2007) pp. 123–146.
  * Segal and Goodman [1965] I. E. Segal and R. W. Goodman, Journal of Mathematics and Mechanics 14, 629 (1965).
  * Masuda [1972] K. Masuda, Publications of the Research Institute for Mathematical Sciences 8, 207 (1972).
  * Murata [1973] M. Murata, Journal of the Mathematical Society of Japan 25, 10.2969/jmsj/02540556 (1973).
  * Białynicki-Birula [1998] I. Białynicki-Birula, Physical Review Letters 80, 5247 (1998).
  * Białynicki-Birula and Białynicka-Birula [2009] I. Białynicki-Birula and Z. Białynicka-Birula, Physical Review A 79, 032112 (2009).
  * Paley and Wiener [1934] R. C. Paley and N. Wiener,  _Fourier Transforms in the Complex Domain_ (American Mathematical Society, New York, 1934).
  * Buller and Collins [2009] G. S. Buller and R. J. Collins, Measurement Science and Technology 21, 012002 (2009).
  * Saari  _et al._ [2005] P. Saari, M. Menert, and H. Valtna, Optics Communications 246, 445 (2005).
  * Saari [2012] P. Saari, in  _Quantum Optics and Laser Experiments_ (InTech, 2012).
  * Gulla and Skaar [2021] J. Gulla and J. Skaar, Physical Review Letters 126, 073601 (2021).
  * Ryen  _et al._ [2022] K. Ryen, J. Gulla, and J. Skaar, International Journal of Theoretical Physics 61, 10.1007/s10773-022-05133-7 (2022).
  * Gulla  _et al._ [2023] J. Gulla, K. Ryen, and J. Skaar, arXiv:2109.06472v2 (2023).
  * Kuhn  _et al._ [2002] A. Kuhn, M. Hennrich, and G. Rempe, Physical Review Letters 89, 067901 (2002).
  * McKeever  _et al._ [2004] J. McKeever, A. Boca, A. D. Boozer, R. Miller, J. R. Buck, A. Kuzmich, and H. J. Kimble, Science 303, 1992 (2004).
  * Scheel [2009] S. Scheel, Journal of Modern Optics 56, 141 (2009).
  * Eisaman  _et al._ [2011] M. D. Eisaman, J. Fan, A. Migdall, and S. V. Polyakov, Review of Scientific Instruments 82, 071101 (2011).
  * Wang  _et al._ [2019] H. Wang, Y.-M. He, T.-H. Chung, H. Hu, Y. Yu, S. Chen, X. Ding, M.-C. Chen, J. Qin, X. Yang, R.-Z. Liu, Z.-C. Duan, J.-P. Li, S. Gerhardt, K. Winkler, J. Jurkat, L.-J. Wang, N. Gregersen, Y.-H. Huo, Q. Dai, S. Yu, S. Höfling, C.-Y. Lu, and J.-W. Pan, Nature Photonics 13, 770 (2019).
  * Sinha  _et al._ [2019] U. Sinha, S. N. Sahoo, A. Singh, K. Joarder, R. Chatterjee, and S. Chakraborti, Optics and Photonics News 30, 32 (2019).
  * Debierre  _et al._ [2015] V. Debierre, I. Goessens, E. Brainis, and T. Durt, Physical Review A 92, 023825 (2015).
  * Debierre and Durt [2016] V. Debierre and T. Durt, Physical Review A 93, 023847 (2016).
  * Hegerfeldt [1994] G. C. Hegerfeldt, Physical Review Letters 72, 596 (1994).
  * Hegerfeldt [1997] G. C. Hegerfeldt, arXiv:quant-ph/9707016 (1997).
  * Karpov  _et al._ [2000] E. Karpov, G. Ordonez, T. Petrosky, I. Prigogine, and G. Pronko, Physical Review A 62, 012103 (2000).
  * Mostafazadeh and Zamani [2006] A. Mostafazadeh and F. Zamani, Annals of Physics 321, 2183 (2006).
  * Babaei and Mostafazadeh [2017] H. Babaei and A. Mostafazadeh, Journal of Mathematical Physics 58, 082302 (2017).
  * Hodgson  _et al._ [2022] D. Hodgson, J. Southall, R. Purdy, and A. Beige, Frontiers in Photonics 3, 10.3389/fphot.2022.978855 (2022).
  * Hawton [2023] M. Hawton, Physical Review A 107, 013711 (2023).
  * Mandel and Wolf [1995] L. Mandel and E. Wolf,  _Optical coherence and quantum optics_ (Cambridge university press, 1995).
  * Garrison and Chiao [2008] J. C. Garrison and R. Y. Chiao,  _Quantum Optics_ (Oxford University Press, 2008).
  * Federico and Jauslin [2023] M. Federico and H. R. Jauslin, Journal of Physics A: Mathematical and Theoretical 56, 235302 (2023).
  * Fabre and Treps [2020] C. Fabre and N. Treps, Reviews of Modern Physics 92, 035005 (2020).
  * Federico  _et al._ [2022] M. Federico, V. Dorier, S. Guérin, and H. R. Jauslin, Journal of Physics B: Atomic, Molecular and Optical Physics 55, 174002 (2022).
  * Białynicki-Birula and Białynicka-Birula [2013] I. Białynicki-Birula and Z. Białynicka-Birula, Journal of Physics A: Mathematical and Theoretical 46, 053001 (2013).
  * Fabre and Cortiñas [2023] C. Fabre and R. G. Cortiñas,  _Quantum Processes and Measurement_ (Cambridge University Press, 2023).
  * Pauli [1980] W. Pauli,  _General Principles of Quantum Mechanics_ (Springer Berlin Heidelberg, 1980).
  * Nisbet-Jones  _et al._ [2011] P. B. R. Nisbet-Jones, J. Dilley, D. Ljunggren, and A. Kuhn, New Journal of Physics 13, 103036 (2011).
  * Talvila [2001] E. Talvila, The American Mathematical Monthly 108, 544 (2001).
  * Cartan [1963] H. Cartan,  _Elementary Theory of Analytic Functions of One or Several Complex Variables_ (Addison-Wesley Educational Publishers Inc, 1963).
  * Masuda [1967] K. Masuda, Proceedings of the Japan Academy, Series A, Mathematical Sciences 43, 10.3792/pja/1195521603 (1967).


[1]Max Schwarzer These authors contributed equally to this work.

[2]Jesse Farebrother These authors contributed equally to this work.

[6]Sergei Kalinin [5]Igor Mordatch [1,5]Pablo Samuel Castro [4]Kevin M.
Roccapriore

1]Université de Montréal, Mila 2]McGill University, Mila

3]Reliant AI

4]Oak Ridge National Laboratory 5]Google DeepMind 6]Department of Materials
Science and Engineering, University of Tennessee, Knoxville, TN 37920

# Learning and Controlling Silicon Dopant Transitions in Graphene using
Scanning Transmission Electron Microscopy

<EMAIL_ADDRESS><EMAIL_ADDRESS>Joshua Greaves
<EMAIL_ADDRESS>Ekin Dogus Cubuk<EMAIL_ADDRESS>Rishabh Agarwal
<EMAIL_ADDRESS>Aaron Courville<EMAIL_ADDRESS>Marc G.
Bellemare<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [  [  [

###### Abstract

We introduce a machine learning approach to determine the transition dynamics
of silicon atoms on a single layer of carbon atoms, when stimulated by the
electron beam of a scanning transmission electron microscope (STEM). Our
method is data-centric, leveraging data collected on a STEM. The data samples
are processed and filtered to produce symbolic representations, which we use
to train a neural network to predict transition probabilities. These learned
transition dynamics are then leveraged to guide a single silicon atom
throughout the lattice to pre-determined target destinations. We present
empirical analyses that demonstrate the efficacy and generality of our
approach.

###### keywords:

Microscopy, Machine Learning

##  1 Introduction

Sub-atomically focused electron beams in scanning transmission electron
microscopes (STEMs) can induce a broad spectrum of chemical changes, including
defect formation, reconfiguration of chemical bonds, and dopant insertion.
Several groups have shown the feasibility of direct atomic manipulation via
electron beam stimulation, which holds great promise for a number of
downstream applications such as material design, solid-state quantum
computers, and others [16, 29, 3, 30, 4]. One of the challenges for advances
in this space is that these types of atomic manipulation rely on manual
control by highly-trained experts, which is expensive and slow.

The ability to accurately automate this type of beam control could thereby
result in tremendous impact on the feasibility of atomic manipulation for real
use cases. A critical requirement for this automation is accurate estimation
of the transition dynamics of atoms when stimulated by focused electron beams.
To date, the microscopy community has relied on heuristic estimates for these
transition dynamics, with anecdotal evidence of their accurateness. Indeed,
the common practice has been to use the physically intuitive, but heuristic,
assumption that the optimal beam position is directly on a neighboring atom.

In this paper we present a technique for estimating these atomic transition
dynamics using machine learning techniques on collected observations. Our
approach consists in a sequence of steps ultimately resulting in a probability
distribution over possible beam positions relative to a particular atom,
conditioned on the atom’s prior position as well as the electron beam’s
location and dwell time. We demonstrate the effectiveness of our approach by
using it to estimate the transition dynamics of a silicon atom on a carbon
lattice (graphene).

As evidence that this process can improve beam control, we use the estimated
transition probabilities to automate the sequential control of a silicon atom
on graphene towards a pre-specified target position. Modern STEMs are capable
of this type of automation, and our work paves the way for future advances in
automated atomic control.

##  2 Related work

Atomic manipulation was first demonstrated by Stroscio and Eigler [26] by
using the tip of a scanning tunneling microscope (STM) to position individual
Xenon atoms on the surface of a single crystal surface to form the IBM company
logo. Further demonstrations, such as quantum corrals and molecular cascades,
have demonstrated the potential of the method. Perhaps the application that
has attracted the most interest is in using tip-induced atomic motion as an
enabling tool for the fabrication of P and other atoms in Si qubits, the
building blocks for quantum computers.

Despite the feasibility of manual control with STM tips, this type of atomic
manipulation is limited to metallic/conducting surfaces. On the other hand,
while STEMs can manipulate atoms embedded within a several layer thick
specimen, it is still a rather haphazard and unpredictable process relative to
using STM tips. To date, electron beam induced effects (with a STEM) have been
studied purely by human operation, most typically by scanning a raster pattern
(where the electron dose tends to be concentrated non-uniformly on one side of
the image) in a selected field of view. The more sophisticated experiments
involve manual positioning of the electron beam by a human, but this kind of
motion is unpredictable and unreliable, and useful statistics are challenging
(if not impossible) to glean from experiments conducted in this manner. We
note that other experiments have been performed which control the electron
beam in non-standard trajectories, effectively performing direct-write beam
patterning processes – but with the critical point that the atomic landscape
(i.e., position of atoms) is not considered [6, 7, 8].

The potential of the electron beams of scanning transmission electron
microscopes to affect matter on the atomic level has been recognized since the
early days of the technique. Most of these effects have been generally
classified as a beam damage, denoting unwanted changes in materials structure
induced by the beam. Indeed, minimization of beam damage, along with the need
to increase spatial and energy resolutions, remains one of the three primary
drivers behind STEM development, having spurred the high-voltage machines of
the 1980s and 1990s and the aberration corrected low-voltage machines of the
last two decades. It was also discovered that electron beam effects can be far
more subtle, including crystallization and amorphization of oxides and
semiconductors [19, 33, 12, 22, 9].

The emergence of aberration-corrected STEMs have made the atomic-resolution
imaging relatively routine, and spurred a new wave of electron beam matter
manipulation on the atomic level. The electron beam was shown to be able to
deposit single atoms from chemisorbed species [2] and form ordered vacancy
arrays [14]. Similarly, electron beams have been shown to induce direct atomic
motion and creation of functional defects [1, 34, 27].

The combination of simple beam control and feedback systems has enabled the
direct assembly of crystalline materials with a single unit plane precision
via directed crystallization and amorphization [15]. These systems have also
demonstrated potential for direct single atom dopant movement [16], and
finally, the controlled manipulation of Bismuth dopants in bulk silicon [11].
In 2016, it was proposed that the combination of machine learning with
electron beam manipulation can become a third paradigm for direct atomic
construction [17]. In 2017, Dyck et al [3], Susi et al [28] and Susi et al
[29] demonstrated single atom manipulation and insertion experiments for
silicon in graphene [29, 3, 30, 4], an approach soon extended to direct atomic
assembly of homo- [4] and hetero-atomic artificial molecules [5].

A number of theories for beam manipulation have been proposed, including those
based on phonon-assisted knock-on and electronic excitations. However the
causal relationship of the electron beam position relative to the silicon atom
has only been suggested and demonstrated anecdotally. Physical intuition
dictates that the damage mechanism is primarily through momentum transfer or
so-called “knock-on” processes; therefore, the ideal placement of the electron
beam would seem to be positioned exactly centered on a carbon (first)
neighbor. While this is intuitive and appears to have been a successful route
by multiple groups, damage mechanisms tend to be complex and are dictated by
more than one process. For example, ionization or sputtering processes may be
occurring as well, meaning it is unclear if the suggested beam position is
actually the ideal one for inducing the most efficient transition of a silicon
hop. Moreover, the anecdotal but physically intuitive rule of placing the
electron beam on the center of a carbon neighbor is mostly valid only for a
direct Si substitution (i.e., 3-fold coordinated silicon). For any other
configuration, the rules are already not the same, and the optimal beam
position for causing a transition event is not clear.

##  3 Problem description

Our system consists of graphene: a single layer of carbon atoms arranged in
3-fold configuration (i.e. every carbon atom is connected to three other
carbon atoms). On this lattice, a single silicon atom (hereafter referred to
as the dopant) has taken the place of one of the carbon atoms. We focus an
electron beam on a position in the area spanned by the dopant and its three
carbon neighbours for a specified amount of time (referred to as the dwell
time). This electron beam stimuli can result in the dopant moving to one of
its neighbours (by trading places with the respective carbon atom), or in the
configuration remaining unchanged111Note that for long dwell-times this
”unchanged” outcome can be a result of the dopant moving twice: once to its
neighbour and then back to its original position.. This configuration has been
extensively explored, and hence is an ideal system for exploring this type of
automation [32, 3, 20, 21]. We provide an illustration of this configuration
in the left panel of Figure 1.

Our objective is to learn a probability distribution over the position of the
dopant, conditioned on its current position, beam location, and beam dwell
time. If accurate, we can use this distribution to determine the optimal beam
location and dwell time so as to induce the dopant to move to one of its
neighbouring positions. In the center panel of Figure 1 we can see a heat map
depicting the transition probabilities for each of the dopant’s neighbours
(differentiated with colors) for varying beam locations.

Equipped with these probability maps, we can repeatedly induce transitions of
the dopant to neighbouring atoms, resulting in a full trajectory. In other
words, we can use a fully greedy strategy to move the dopant to any pre-
specified target position on the lattice via a simple shortest-distance path.
In the right panel of Figure 1 we depict one such possible trajectory.

Figure 1: Left: Illustration of pristine graphene with a single dopant. Carbon
atoms depicted with black circles, while the silicon atom in blue; Center:
Contour map learned by our method, depicting the probability of transitioning
to each of the neighbours for different beam positions. The distributions for
each of the neighbours are differentiated using three colours, and the numbers
indicate neighbour ordering as discussed in Section 4.3; Right: Example
trajectory of dopant towards a goal position.

##  4 Description of Methodology

###  4.1 Data collection

Our methodology relies on real data collected with a STEM device, so it is
important that the data gathered is informative for the task at hand. Since we
are concerned with the transition dynamics of the dopant, our data collection
approach is as follows:

  1. 1.

Acquire an initial image of the graphene. For this, we used a field of view of
3 nm with the dopant in the center. A sample image is shown on the bottom-left
panel of Figure 3.

  2. 2.

Sample a position uniformly within a 2.84 Å radius222Chosen to be the
cumulative length of two carbon-carbon bonds. of the atom. In the top-center
panel of Figure 3 we display the normalized beam positions.

  3. 3.

Focus the electron beam at that position for a dwell time drawn from a
distribution typically sampled from 1 to 10 seconds. The top-left corner of
Figure 3 displays the distribution of dwell times.

  4. 4.

Acquire a final image of the graphene.

Figure 2: An overview of the full pipeline for learning the transition
probabilities.

In this section we detail the methodology we used to learn the transition
dynamics of electron-beam induced atomic manipulation, illustrated in Figure
2. We used a Nion UltraSTEM 100, which grants us access to nearly every
microscope control via a Python API. The silicon dopant atoms have been
inserted into the lattice in a previous step at a higher (100 kV) accelerating
voltage, which is described elsewhere [25].

There are a few considerations that are worth mentioning. First, it is
important to gather multiple samples using the same beam position and dwell
time, as transitions are probabilistic. Second, image acquisition is done
using the same electron beam, which implies that imaging itself can cause a
transition; to mitigate this, the imaging electron dose should be minimized.
Third, placing the electron beam in a known position relative to the silicon
should be conducted in a controlled and automated fashion – for this, atomic
coordinates must be known with high confidence in as close to real time as
possible, and flexible control of the beam position is needed. This last
consideration is addressed by using deep ensemble neural networks to identify
the atomic species and positions both quickly and reliably [10].

###  4.2 Atomic alignment

In the bottom-left panel of Figure 3 we display one of the raw images captured
with the STEM device, and to the right of it the processed output. The
processed image is an “idealized” configuration as illustrated in Figure 1,
which are easier to operate on. While it is relatively simple to identify the
atoms in a single image, there are some challenges that arise when using more
than one raw image, as discussed in Section 4.1: the graphene sheet may have
physically moved between image acquisition steps (specimen drift), and there
may be aberrations (such as warping) caused by the electron beam. To be able
to use the images acquired for learning transition dynamics, we need not only
to detect the atoms in each of the separate images, but also map each of the
atoms from one image to the next.

We considered three solutions to this problem. The classical approach is to
take the cross-correlation between the two scans and estimate the drift to be
the arg max of this cross-correlation. Unfortunately, as images of this system
are generally dominated by the very bright silicon dopants, this has the net
effect of always aligning the dopant positions between time steps – leading to
a conclusion that the dopants never move, which is known to be false. A
second, more sophisticated alternative is to use the iterative closest points
(ICP) algorithm on the extracted atom positions (a similar technique as used
by Roccapriore et al [23]). This allows the algorithm to equally weigh the
dopant and non-dopant atoms, simply aligning the lattices together. This
method leads to acceptable alignments in most cases, but we found it to be
quite sensitive to failures in atom detection.

Figure 3: Top row: Histogram of dwell times (left) and standardized beam
positions (center) used for data acquisition; beam position data augmented for
three neighbours (right). In the beam position plots, grey and colored circles
represent negative and positive transitions, respectively. The colored circles
with a black border represent the neighbouring atoms. Bottom left: Raw image
acquired from the STEM (left), and processed image after atom detection
(right). Bottom right: The result of conducting multiple alignment iterations
with our trained convolutional aligner. Orange dots reflect estimated atom
positions in the previous scan, blue dots reflect atom positions in the
current scan, and red dots are dopant positions. As can be observed, the
discrepancies diminish with additional iterations.

To address these issues, we instead use a denoising convolutional neural
network [13] that solves the alignment problem directly from scans of the
system. Given a stack of image observations $(o_{1},...,o_{n})$, the network
is trained to predict the drift $d_{n}$ between $o_{n-1}$ and $o_{n}$, as a
two-dimensional vector. Historical observations $o_{1},\ldots,o_{n-2}$ serve
to provide context and noise reduction, but are not used in the loss; as they
have already been approximately aligned with $o_{n-1}$, they provide
additional information about the needed shift. After the network is applied,
we take its prediction $\hat{d}_{n}$ and shift the current observation by it
to align it to the prior observations, and add it to the array of historical
observations.

To train this network, we generated a dataset of synthetic trajectories, each
consisting of sequential image observations of a doped graphene system under
random, correlated drift. To add robustness, we applied both trajectory-wide
augmentation by randomly dropping atoms, introducing regions of bright
contamination, and adding synthetic large holes to the system. We then
simulated drift, treating the direction and magnitude of the drift as a
temporally correlated random variables. Finally, we took a series of synthetic
scans with the simulated cumulative drifts applied; for each scan, we also
randomly perturbed the system, occasionally dropping or moving atoms. We
parameterized the network as six convolutional layers followed by
downsampling, followed by a single fully-connected layer. We trained drift
correction prediction with mean squared error.

Note that the estimation of ${d}_{n}$ is a denoising task. Given this, we can
iteratively apply our network for more precise drift correction. We find that
this has a very large impact on performance. Qualitatively it is clear that
the alignment quality increases with additional iterations (see bottom-right
panel of Figure 3). Quantitatively, a transition model trained on only single-
step-aligned data produces markedly different estimates of the optimal beam
position, which are not able to successfully cause transitions when applied to
the greedy controller (see “Past Neighbor” in Figure 4).

###  4.3 Data filtering, augmentation, and structure

As there are many possible sources of noise and error in our data, we apply
several types of filtering prior to training. Out of a starting data set size
of $6,754$ examples, we discard transitions where: (1) There is no recorded
beam position ($793$ examples); (2) There is not exactly one (1) detected
dopant atom before and after the transition ($4$ examples); (3) The dopant
does not have the expected number of neighbors (3) before and after the
transition, either because they are absent, indicating non-pristine graphene,
or were not detected ($3,593$ examples); and (4) the neighbors after the
transition did not approximately align with the neighbors before the
transition (chosen as an average distance $>$ 0.71Å, $411$ examples).

These led us to discard approximately 80% of the data we received from the
microscope, resulting in a final dataset of $1,953$ examples. Once filtering
has been done, we further post-process our observations to create a
uniformized training set. To express beam positions in a consistent format, we
translated them to a frame-of-reference relative to the current position of
the dopant to be moved, with this atom at the origin. We then label the
neighbor closest to the beam position as neighbor 1; we rotate the system,
including the beam, such that this atom lies on the x-axis. The neighboring
atoms are then numbered accordingly in counter-clockwise order; these are the
indices we predict in our classification. We denote “no movement” as index 0
in the classification task. This labeling is illustrated in the top-center
panel of Figure 1.

In the structure noted above, there is no systematic difference between the
three neighboring atoms; only their distance from the beam separates them. As
we expect our system to be invariant to both rotation and reflection, we
enforce this by adding data augmentation. To do this, we first reflect across
the x-axis with 50% probability. This exchanges the second and third
neighbors. We then apply 0°, 120°or 240°rotations with equal probability,
rotating the neighbor indices accordingly. This leads to an effective sixfold
increase in data coverage. The top-right panel of Figure 3 visualizes the beam
positions post-augmentation.

###  4.4 Learning the transition dynamics

Experimentally, we modeled our problem as a classification task, estimating
$P(S^{\prime}|s_{0},a)$, where $S^{\prime}$ is the position of the silicon at
the next time step, $s_{0}$ is its current position, and $a$ is the beam dwell
action chosen by the user, specified as a two-dimensional coordinate $x$ and a
duration in seconds $\Delta t$. To constrain our predictions to respect the
observed physical reality (i.e., that the probability of moving to another
state should be monotonically increasing in $\Delta t$), we further formulate
the problem as predicting transition rates [31]. We decompose this as a total
rate $\lambda$ and a categorical distribution over possible next states $y$;
this corresponds to decomposing $P(S^{\prime}|a)$ into a probability of some
transition occurring, $P(S^{\prime}\neq s_{0}|s_{0},a)$, and a distribution
over which next state is chosen if a transition occurs,
$P(S^{\prime}|s_{0},a,S^{\prime}\neq s_{0})$. We parameterize
$P(S^{\prime}\neq s_{0}|s_{0},a)$ as $1-e^{-\lambda\Delta t}$ – i.e., as an
exponential cumulative distribution function. If desired, per-neighbor rates
are simply given by $\lambda y$. We then formulate maximum likelihood loss
functions for the total rate and distribution over neighbors:

$\displaystyle\mathbf{J}_{rate}$ $\displaystyle=-\mathbb{I}_{(S^{\prime}\neq
s_{0})}\cdot(-\lambda\Delta
t)-\mathbb{I}_{(S^{\prime}=s_{0})}\log(1-\exp(-\lambda\Delta t)$ (1)
$\displaystyle\mathbf{J}_{next}$ $\displaystyle=-\mathbb{I}_{(S^{\prime}\neq
s_{0})}(S^{\prime}\cdot\log y)$ (2) $\displaystyle\mathbf{J}_{total}$
$\displaystyle=\mathbf{J}_{rate}+\mathbf{J}_{next}\,.$ (3)

Given the above loss function, we trained a three-layer neural network using
Adam with weight decay [18] and ReLU hidden layer nonlinearities, trained with
the cross-entropy loss for 500 epochs using batch size 256. We predicted
$\lambda$ with a softplus activation and $y$ with a softmax. We normalize all
inputs to the network prior to training, for stability.

To improve the robustness of our transition model to initializations, we
trained an ensemble of transition predictors of bootstrap-resampled datasets.
In addition to giving us the ability to esimate uncertainties, this
significantly improved our overall accuracy and robustness To support more
rapid inference, we also distilled this ensemble to a single transition
predictor. By using widely-sampled random beam positions and a large number of
training steps, this distillation could be made to fairly precisely match the
ensemble predictions.

We display an example of the probability contours found by our system in the
center panel of Figure 1. Note that we are overlaying three different
probability distributions (one for each neighbour), where the colours are used
to distinguish them. Our main finding confirms what was anecdotally held to be
true by the microscopy community:

To induce the dopant to transition to one of its neighbours, the optimal beam
placement is directly on the neighbour, with a 50.0% (95% CI 0.246 to 0.754)
probability of causing a transition with a five second dwell time, using a
nominal beam current of 90 pA.

Perhaps more relevant than dwell time is the number of electrons emitted,
which is a function of both dwell time and the beam current. For our
experiments we used a beam current of 90 pA, resulting in approximately 3
billion electrons in a five second period.

##  5 Empirical evaluation

|
---|---

Figure 4: Left: The proportion of transitions that induced the intended
transition in under a $5$ second dwell. Right: The beam placement strategies
considered in our experiments.

While our main finding is consistent with previously held beliefs in the
community, we do not have ground-truth data for the learned transition
probabilities to quantitatively assess the accuracy of our predictions.
However, as discussed in Section 3, the purpose of learning these transition
probabilities is to be able to automate atomic manipulation, so in this
section we evaluate the efficacy of our learned transition functions for this
purpose. Specifically, our experiments are conducted as follows:

  1. 1.

Start from a configuration with one dopant and 3-fold connections to
neighbours, with the field-of-view (FOV) centered at the dopant.

  2. 2.

Pick an arbitrary carbon atom as the goal position.

  3. 3.

Focus the beam on a position dictated by one of the strategies (defined
below).

  4. 4.

Acquire an image to determine if we caused a transition.

  5. 5.

Repeat above steps until the dopant has arrived at the goal position, or until
we have reached the maximum allowable attempts.

In order to determine that on-neighbour is in fact the optimal beam placement,
we use a few different strategies for beam placement, detailed below and
illustrated in Figure 4 (Right).

  * •

On neighbour: our proposed optimal strategy.

  * •

Short of neighbour: place the beam in between the dopant and its neighbour.

  * •

Past neighbour: place the beam beyond the neighbour atom

  * •

Perpendicular offset: offset the beam perpendicularly from neighbour.

  * •

Learned dynamics (less data): we ran the learning method detailed in Section
4.4, but with approximately half of the collected data. When doing so, the
resulting ”optimal” beam placement was just past the neighbour.

The first controller was chosen as the optimal policy according to our most
recent learned transition models; the second, third, and fourth were chosen to
test how deviations from this policy affected results, in a systematic way.
The fifth and final controller represents the optimal policy according to an
early iteration of our learned models, and serves to validate the progress
made as additional data was introduced.

For consistency, we used a constant $5$s dwell time at a nominal beam current
of 90 pA for all agents. We measured the number of times the electron beam was
used to try and induce a dopant transition and report the findings in Figure
4. We observe that the on neighbor strategy induces the intended transition on
average approximately $50\%$ of the time, whereas the other approaches are
below $25\%$.

Figure 5: Experimental realization of automation of silicon motion throughout
graphene using learned transition rates. HAADF-STEM image snapshots of five
steps are shown in a, while corresponding predicted images are shown below in
b. Images are aligned to the initially selected goal position, denoted by red
arrows, and beam positions shown in green circles. The silicon is the bright
atom in a due to its higher atomic number, and is shown by an orange circle in
the predicted images. The path summary of the silicon atom between each step
is shown in c. Scale bars all 2 Angstroms.

Automated silicon manipulation was experimentally realized on a Nion
UltraSTEM100 microscope, and an example using the on neighbor protocols is
shown in Figure 5. Here, a goal location - shown by red arrows in Figure 5a \-
is randomly chosen relative to the silicon, and the controller directs the
electron beam on carbon neighbors to drive the silicon to the goal position.
To accomplish this task automatically, the atomic coordinates must be known in
as close to real time as possible; here, deep ensemble neural networks were
used which have been shown to be both fast and robust even to out of
distribution effects [10, 24]. These predictions are shown in Figure 5b where
the silicon is depicted by a larger orange circle. Fast access to the atom
positions allows the chosen controller to then place the beam at appropriate
positions relative to the silicon to fully automate this process. Note that
the first step in Figure 5 has caused the silicon to jump multiple sites,
whereas the remaining steps cause single atomic jumps. Real-time detector
feedback was not used to indicate if the silicon transitioned earlier than the
designated dwell time, and therefore this was one limitation in the current
work, which can explain the (albeit lower probability) longer initial jump. We
also note that this behavior can be exploited for more efficient beam
manipulation; i.e., one beam position can intentionally induce more than one
atomic jump.

##  6 In-depth empirical analyses

###  6.1 Synthetic data

To test the learning behavior of our algorithm beyond our relatively small
real-world dataset, we generated many datasets of simulated microscope
interactions. We use these datasets to evaluate the scaling behavior of our
algorithm, demonstrate that models ultimately converge to accurate estimates
with sufficient data, and to evaluate our algorithmic and hyperparameter
decisions.

|
---|---

Synthetic Rate Functions

Figure 6: Left:  Number of positive transitions in dataset versus the sum of
squared prediction error. Right:  A visualization of some of the synthetic
rate functions overlayed on a sheet of graphene. The contours represent the
rate of transitioning to the associated neighbor.

To generate these datasets, we sampled synthetic transition probability
distributions as mixtures-of-Gaussians,333Chosen because they generalize the
family of rate functions we expect to see in the real world by permitting
multimodality, with the Gaussian nature of the distributions being sensible as
the electron beam is approximately a Gaussian beam. each giving the non-
normalized transition rates to neighboring states given a certain beam
position (equivalent to the predicted per-state rates $\lambda y$ in our
neural network model). Figure 6 (Right) is representative of the type of
synthetic rate functions used throughout this section. Doing so grants us
ground-truth data for transition probabilities, allowing us to quantitatively
assess the accuracy of our learned transition model.

When creating a dataset, we first generate a simulated graphene sheet with a
single dopant. We then uniformly sample random actions within 2Å of the
silicon and simulate the transitions of the silicon according to the synthetic
rate function, continuing until a certain number of positive transitions have
been observed; as positive transitions are generally far rarer than negative
transitions, they are the critical determiner of effective dataset size. To
simplify comparisons between datasets, we enforce that the synthetic rate
functions for each dataset have the same maximum value, preventing us from
sampling datasets that are entirely positive transitions or almost entirely
negative transitions.

We can then use a synthetic dataset to evaluate a learning algorithm by
training on it and directly comparing its predictions (here, the predicted
rates $\lambda y$) to the known true values across a large grid of beam
positions (a uniformly-spaced 2d grid surrounding the dopant). This provides a
measure of the total accuracy of the learned rate function across a full range
of possible beam locations, even those that would not normally be chosen when
attempting to induce a transition.

###  6.2 Data scaling

To test the data scaling of our model – and show that it converges to near-
perfect predictions in the limit of large datasets – we use this evaluation
procedure at a range of data scales. We start by collecting 30 datasets with
at least 1,000 positive transitions each. We then train our models on a range
of subsampled scales from these datasets, and report a bootstrapped confidence
intervals over our our 30 synthetic datasets at each scale in Figure 6 (Left).
As expected, we find that our model has very poor performance with limited
amounts of data, but rapidly improves once hundreds of positive transitions
are available, matching our own experience of the model’s improvement as
additional data was collected.

To evaluate the utility of these learned models for control, we measure the
regret, or total suboptimality caused by acting according to our learned model
compared to the true rate function. To do this, we take the argmax of a
learned model $\hat{R}$, and evaluate the true synthetic rate function $R$ at
this point. We then define the regret as the difference between the true
maximum of the rate function and this value:

$\displaystyle\text{Regret}(\hat{R})\triangleq\max_{a}R(a)-R(\arg\max_{a}\hat{R}(a))$
(4)

Intuitively, this measure captures how much slower control according to a
learned rate function would be than control according to the true rate
function.

|
---|---

Figure 7: Left:  Positive transitions in dataset versus the incurred regret
for our rate model versus a $k$-nearest neighbors and a random forest
classifier. Right:  Comparing the rate of convergence when measuring maximum
regret versus mean squared error.

Figure 7 (Right) depicts the relative rate of convergence between the mean-
square error and the regret of the true rate function. We find that the model
tends to converge more quickly to the global shape of the rate function (as
measured by the mean squared error) as compared to true argmax (as measured by
the regret). This matches our empirical experience where the learned rate
functions can look similar although their argmax can differ (e.g., Learned
Rates (Less Data) vs Learned Rates in Figure 4 (Right)).

###  6.3 Hyperparameter and algorithmic evaluation

Beyond scaling experiments, we leverage this synthetic evaluation framework to
validate several hyperparameter decisions that we found to be important in the
real data setting. Firstly, in Figure 8 we examine the importance of training
an ensemble of transition prediction models and then distilling them. We find
that distillation significantly improves regret at all data scales on
synthetic data; we visualize examples of distilled and non-distilled models
trained on real data, to illustrate how this might take place. We find that
the result is largest in the intermediate-data regime when meaningful
transition models can be learned but inaccuracies are common, matching our
experience in the progression of real data experiments.

We also evaluate other possible machine learning paradigms, including families
of classical algorithms such as K-nearest-neighbors and random forests in
Figure 7 (Left). We find that these algorithms are also capable of learning
reasonable transition models, but generally require significantly more data to
do so and ultimately converge to worse solutions, suggesting that the
inductive biases provided by modern neural networks are important to this
task.

|
---|---

Distilled Rate Functions
---
Non-Distilled Rate Functions

Figure 8: Comparisons of transition prediction models trained a single models
(top) and as a distilled ensemble of multiple models (bottom) on a suite of
synthetic datasets (left) and on real data (right). Note that the distilled
ensemble has both better overall regret on synthetic data and shows more
regular predicted rates on real data.

##  7 Conclusion

The last note left by Richard Feynman stated “What I cannot create, I do not
understand.” Building solid state quantum computers, creating nano-robots, and
designing new classes of biological molecules and catalysts alike requires the
capability to manipulate and assemble matter atom by atom, probe the resulting
structures, and connect them to the macroscopic world; all this necessitates
accurate estimates of the transition dynamics induced by sub-atomically
focused electron beams. Until now, the elements of relevant knowledge have
been limited to a few research groups, and atomic manipulation has been
performed via direct control by human operator one beam positioning at a time.
The characteristic timescale of human-operated experiments vastly exceeds the
intrinsic latency of the electron microscope, for which hundreds of
fabrications steps per second should be possible. Similarly, human control
necessarily lacks precision, reproducibility, and systematic error correction
capabilities. While sufficient for a proof of concept, atomic scale
fabrication with the precision and throughput necessary for applications such
as nanopore fabrication for protein sequencing, molecule screening platforms
for physics and biology, and particularly quantum communication, sensing, and
computing devices requires moving beyond the current human control paradigm.

Our work is a robust first step for determining transition probabilities via
machine learning, and paves the way for further advances in this space. The
scenario we considered in this work is somewhat idealized: we limit ourselves
to single dopant and 3-way lattices for our learned dynamics. Nevertheless,
these settings allowed us to confirm, via a data-driven approach, the commonly
held belief that placing the electron beam directly on the neighbour has the
highest probability of inducing a transition of the dopant. Going forward, we
will be exploring broader settings: multiple dopants, graphene with 4-way
connections and aberrations (such as holes).

Supplementary information

The code we used for this work is publicly available at
https://github.com/google/putting-dune.

Acknowledgments

Autonomous STEM research was supported by the Center for Nanophase Materials
Sciences (CNMS) proposal [CNMS2022-B-01647], which is a US Department of
Energy, Office of Science User Facility at Oak Ridge National Laboratory. The
authors would also like to thank their colleagues at Google DeepMind for their
useful feedback on this work.

## References

  *   * Cretu et al [2012] Cretu O, Rodríguez-Manzo J, Demortière A, et al (2012) Electron beam-induced formation and displacement of metal clusters on graphene, carbon nanotubes and amorphous carbon. Carbon 50(1):259–264
  * van Dorp et al [2012] van Dorp WF, Zhang X, Feringa BL, et al (2012) Molecule-by-molecule writing using a focused electron beam. ACS Nano 6(11):10076–10081
  * Dyck et al [2017] Dyck O, Kim S, Kalinin SV, et al (2017) Placing single atoms in graphene with a scanning transmission electron microscope. Applied Physics Letters 111(11):113104. 10.1063/1.4998599, URL https://doi.org/10.1063/1.4998599, https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.4998599/14503703/113104_1_online.pdf
  * Dyck et al [2018] Dyck O, Kim S, Jimenez-Izal E, et al (2018) Building structures atom by atom via electron beam manipulation. Small 14(38):1801771
  * Dyck et al [2019] Dyck O, Ziatdinov M, Lingerfelt DB, et al (2019) Atom-by-atom fabrication with electron beams. Nature Reviews Materials 4(7):497–507
  * Dyck et al [2023a] Dyck O, Lupini AR, Jesse S (2023a) Atom-by-Atom Direct Writing. Nano Letters (23):2339–2346
  * Dyck et al [2023b] Dyck O, Lupini AR, Jesse S (2023b) A platform for atomic fabrication and in situ synthesis in a scanning transmission electron microscope. Small Methods n/a(n/a):2300401
  * Dyck et al [2023c] Dyck O, Yeom S, Lupini AR, et al (2023c) Top-down fabrication of atomic patterns in twisted bilayer graphene (adv. mater. 32/2023). Advanced Materials 35(32):2370228
  * Frantz et al [2001] Frantz J, Tarus J, Nordlund K, et al (2001) Mechanism of electron-irradiation-induced recrystallization in si. Phys Rev B 64:125313. 10.1103/PhysRevB.64.125313, URL https://link.aps.org/doi/10.1103/PhysRevB.64.125313
  * Ghosh et al [2021] Ghosh A, Sumpter BG, Dyck O, et al (2021) Ensemble learning-iterative training machine learning for uncertainty quantification and automated experiment in atom-resolved microscopy. npj Computational Materials 7(1). 10.1038/s41524-021-00569-7, URL https://doi.org/10.1038/s41524-021-00569-7
  * Hudak et al [2018] Hudak BM, Song J, Sims H, et al (2018) Directed atom-by-atom assembly of dopants in silicon. ACS Nano 12(6):5873–5879
  * I. Jencic [1995] I. Jencic IMRMAKM. W. Bench (1995) Electron‐beam‐induced crystallization of isolated amorphous regions in si, ge, gap, and gaas. Journal of Applied Physics 78(2):974–982
  * Jain and Seung [2008] Jain V, Seung S (2008) Natural image denoising with convolutional networks. In: Koller D, Schuurmans D, Bengio Y, et al (eds) Advances in Neural Information Processing Systems, vol 21. Curran Associates, Inc.
  * Jang et al [2017] Jang JH, Kim YM, He Q, et al (2017) In situ observation of oxygen vacancy dynamics and ordering in the epitaxial lacoo3 system. ACS Nano 11(7):6942–6949
  * Jesse et al [2015] Jesse S, He Q, Lupini AR, et al (2015) Atomic-level sculpting of crystalline oxides: Toward bulk nanofabrication with single atomic plane precision. Small 11(44):5895–5900
  * Jesse et al [2018] Jesse S, Hudak BM, Zarkadoula E, et al (2018) Direct atomic fabrication and dopant positioning in si using electron beams with active real-time image-based feedback. Nanotechnology 29(25):255303
  * Kalinin et al [2016] Kalinin SV, Borisevich A, Jesse S (2016) Fire up the atom forge. Nature 539(7630):485–487
  * Kingma and Ba [2015] Kingma D, Ba J (2015) Adam: A method for stochastic optimization. In: International Conference on Learning Representations (ICLR), San Diega, CA, USA
  * Lulli and Merli [1993] Lulli G, Merli PG (1993) Comparison of results and models of solid-phase epitaxial growth of implanted si layers induced by electron- and ion-beam irradiation. Phys Rev B 47:14023–14031. 10.1103/PhysRevB.47.14023, URL https://link.aps.org/doi/10.1103/PhysRevB.47.14023
  * Markevich et al [2020] Markevich A, Hudak B, Lupini A, et al (2020) Uncovering the mechanism for electron-beam manipulation of dopants in silicon. Microscopy and Microanalysis 26(S2):2560–2561. 10.1017/S1431927620022035
  * Markevich et al [2021] Markevich A, Hudak BM, Madsen J, et al (2021) Mechanism of electron-beam manipulation of single-dopant atoms in silicon. The Journal of Physical Chemistry C 125(29):16041–16048
  * Robertson and Jenčič [1996] Robertson I, Jenčič I (1996) Regrowth of amorphous regions in semiconductors by sub-threshold electron beams. Journal of Nuclear Materials 239:273–278. https://doi.org/10.1016/S0022-3115(96)00485-0, URL https://www.sciencedirect.com/science/article/pii/S0022311596004850, proceedings of the JIM ’95 Fall Annual Meeting (117th) on Lattice Defects and Radiation Induced Phenomena
  * Roccapriore et al [2021] Roccapriore KM, Creange N, Ziatdinov M, et al (2021) Identification and correction of temporal and spatial distortions in scanning transmission electron microscopy. Ultramicroscopy
  * Roccapriore et al [2022] Roccapriore KM, Boebinger MG, Dyck O, et al (2022) Probing electron beam induced transformations on a single-defect level via automated scanning transmission electron microscopy. ACS Nano 16(10):17116–17127. 10.1021/acsnano.2c07451, URL https://doi.org/10.1021/acsnano.2c07451
  * Roccapriore et al [2023] Roccapriore KM, Schwarzer M, Greaves J, et al (2023) Discovering the electron beam induced transition rates for silicon dopants in graphene with deep neural networks in the stem
  * Stroscio and Eigler [1991] Stroscio JA, Eigler DM (1991) Atomic and molecular manipulation with the scanning tunneling microscope. Science 254(5036):1319–1326. 10.1126/science.254.5036.1319, URL https://www.science.org/doi/abs/10.1126/science.254.5036.1319, https://www.science.org/doi/pdf/10.1126/science.254.5036.1319
  * Susi et al [2014] Susi T, Kotakoski J, Kepaptsoglou D, et al (2014) Silicon–carbon bond inversions driven by 60-kev electrons in graphene. Phys Rev Lett 113:115501. 10.1103/PhysRevLett.113.115501, URL https://link.aps.org/doi/10.1103/PhysRevLett.113.115501
  * Susi et al [2017a] Susi T, Kepaptsoglou D, Lin YC, et al (2017a) Towards atomically precise manipulation of 2d nanostructures in the electron microscope. 2D Materials 4(4):042004
  * Susi et al [2017b] Susi T, Meyer JC, Kotakoski J (2017b) Manipulating low-dimensional materials down to the level of single atoms with electron irradiation. Ultramicroscopy 180:163–172
  * Tripathi et al [2018] Tripathi M, Mittelberger A, Pike NA, et al (2018) Electron-beam manipulation of silicon dopants in graphene. Nano Letters 18(8):5319–5323
  * Voter [2007] Voter AF (2007) Introduction to the kinetic monte carlo method. In: Radiation effects in solids. Springer, p 1–23
  * Wang et al [2014] Wang WL, Santos EJ, Jiang B, et al (2014) Direct observation of a long-lived single-atom catalyst chiseling atomic structures in graphene. Nano letters 14(2):450–455
  * Yang et al [1997] Yang X, Wang R, Yan H, et al (1997) Low energy electron-beam-induced recrystallization of continuous gaas amorphous foils. Materials Science and Engineering: B 49(1):5–13
  * Yang et al [2014] Yang Z, Yin L, Lee J, et al (2014) Direct observation of atomic dynamics and silicon doping at a topological defect in graphene. Angewandte Chemie International Edition 53(34):8908–8912