# Gamma Radiation Source Localization for Micro Aerial Vehicles with a
Miniature Single-Detector Compton Event Camera

Tomas Baca1, Petr Stibinger1, Daniela Doubravova2, Daniel Turecek2,
Jaroslav Solc3, and Jan Rusnak3, Martin Saska1, Jan Jakubek2 1Authors are with
Czech Technical University in Prague, Technicka 2, 160 00, Prague.
<EMAIL_ADDRESS>are with Advacam s.r.o., U Pergamenky 12, 170
00, Prague.3Authors are with Czech Metrology Institute, Okruzni 31, 638 00,
Brno.

###### Abstract

A novel method for localization and estimation of compact sources of gamma
radiation for Micro Aerial Vehicles is presented in this paper. The method is
developed for a novel single-detector Compton camera, developed by the
authors. The detector is extremely small and weighs only 40 g, which opens the
possibility for use on sub-1 kg class of drones. The Compton camera uses the
MiniPIX TPX3 CdTe event camera to measure Compton scattering products of
incoming high-energy gamma photons. The 3D position and the sub-nanosecond
time delay of the measured scattering products are used to reconstruct sets of
possible directions to the source. The proposed method utilizes a filter for
fusing the measurements and estimating the radiation source position during
the flight. The computations are executed in real-time onboard and allow
integration of the detector info into a fully-autonomous system. Moreover, the
real-time nature of the estimator potentially allows estimating states of a
moving radiation source. The proposed method was validated in simulations and
demonstrated in a real-world experiment with a Cs137 radiation source. The
approach can localize a gamma source without estimating the gradient or
contours of radiation intensity, which opens possibilities for localizing
sources in a cluttered and urban environment.

##  I INTRODUCTION

Localization of ionizing radiation sources has recently become the domain of
autonomous mobile robots. Following the 2011 disaster at the Fukushima Daiichi
Nuclear Power Plant (FDNPP), considerable efforts went into minimizing the
exposure of human workers by sending Unmanned Ground Vehicles into potentially
hazardous areas. Mapping a radiation distribution within an area is
traditionally done with a dosimeter (measures intensity of radiation) while
scanning the whole area. On the other hand, the search for compact sources is
better done with sensors capable of estimating the direction to the source.
The direction estimation is traditionally achieved by pixel radiation
detectors (cameras) with sensor attachments, such as pinhole optics, coded
apertures, and collimators [1]. In recent years, Compton cameras made an
addition to the list of ways to estimate a set of directions to the source,
however, without an external attachment to the sensor. A _traditional_ Compton
camera provides the direction to a source by measuring products of a Compton
scattering reaction within two detectors (scatter and absorber). Information
on the products is used to calculate a set of possible directions to the
source in the form of a cone surface in 3D.

(a)

(b)

Figure 1: A miniature UAV, one of our experimental platforms, equipped with
the MiniPIX TPX3 Compton camera.

We aim to bring the Compton camera to the world of MAVs —
sub-$1\text{\,}\mathrm{kg}$ UAVs. With current state-of-the-art technologies,
MAVs can fly autonomously in cluttered environments, such as urban and indoor
environments and forests. An MAV can pursue an informed exploration of the
environment towards the estimated position of the source if such a real-time
estimate is available.

###  I-A Related work

Several UGVs have been deployed directly inside the damaged reactor buildings
of FDNPP under remote control. Various radiation detection methods have been
tested inside the power plant, including a coded aperture scintillator [2], a
semiconductor digital dosimeter [3], a Compton event camera composed of two
scintillators [4], and a time-of-flight gamma camera [5]. Ground-based robots
offer higher payload capacity and the ability to carry heavier sensory
equipment than most aerial vehicles. On the other hand, these robots tend to
be relatively bulky, struggle to navigate the cluttered corridors and
staircases inside the damaged buildings, and generally move slower than a
multirotor aircraft.

UAVs have been utilized to map the spread of the radioactive material outside
of the power plant. These range from large aircraft weighing more than
$90\text{\,}\mathrm{kg}$ equipped with heavy scintillation detectors [6, 7, 8]
to compact multi-rotors suitable for flying along a pre-defined trajectory
close to the ground [9, 10, 11]. Outside of Japan, several projects have
employed UAVs for radiation intensity mapping around uranium ore mines [12,
13, 14].

In [15], multiple fixed-wing UAVs equipped with miniature scintillators are
used for contour analysis of an irradiated area. Trajectory planning and data
processing are performed offline, contrary to our approach, which estimates
the position of the source in real-time during the flight. In [16], the
contour analysis is tackled using a single multirotor UAV. The contour
analysis uses a Gaussian mixture model to estimate multiple radiation sources’
positions with overlapping intensity fields. The projects mentioned above
utilize unmanned vehicles to deliver a radiation detector into a hazardous
environment. However, the approaches do not respond to measured data in real-
time and thus do not exploit the mobility of UAVs to improve the measurement.

Active path-planning driven by the onboard measurements has been shown in [7]
for an outdoor environment and in [17] for a GPS-denied indoor environment.
Both of these works rely on a scintillator detector to estimate the radiation
intensity in the UAV’s current position. As a result, the employed aerial
platforms have to be large with a payload capacity exceeding
$2\text{\,}\mathrm{kg}$. As in [7], the aerial platform is a
$90\text{\,}\mathrm{kg}$ unmanned aerial helicopter, which significantly
limits its deployment conditions due to personal safety and considerable
minimum distance to obstacles in the environment. Similarly, scintillating
detectors were used in [17] on a mobile ground robot to map a cluttered
environment and to search for radiation sources.

The lack of lightweight radiation detectors with immediate readout capability
severely limits the application potential of aerial dosimetry. However, a new
robotic methodology for motion planning and exploration needs to be developed
to accommodate the proposed measurement system’s specifics. The proposed
system is intended for outdoor and indoor environments, prompting a very
compact and self-sufficient vehicle. We build upon our previous experience
with indoor UAV deployment [18], where such radiation localization system
could be used during exploration of mine shafts and caves [19, 20], in
interiors of large buildings [21, 22], in forests [23] and in nuclear power
plants [24].

$\mathbf{x}$, $\bm{\alpha}$ vector, pseudo-vector, or tuple
$\mathbf{\hat{x}}$, $\bm{\hat{\omega}}$ unit vector or unit pseudo-vector
$\mathbf{\hat{e}}_{1},\mathbf{\hat{e}}_{2},\mathbf{\hat{e}}_{3}$ elements of
the _standard basis_ $\mathbf{X},\bm{\Omega}$ matrix
$x=\mathbf{a}^{\intercal}\mathbf{b}$ inner product of $\mathbf{a}$,
$\mathbf{b}$ $\in\mathbb{R}^{n}$ $\mathbf{x}=\mathbf{a}\times\mathbf{b}$ cross
product of $\mathbf{a}$, $\mathbf{b}$ $\in\mathbb{R}^{3}$
$\mathbf{x}=\mathbf{a}\circ\mathbf{b}$ element-wise product of $\mathbf{a}$,
$\mathbf{b}$ $\in\mathbb{R}^{n}$ $x_{[n]}$ $x$ at the sample $n$
$\mathbf{A},\mathbf{x},\mathbf{\Omega}$ LKF system matrix, state vector,
covariance $\Delta t$ time difference interval, [$\mathrm{s}$]
$\mathbf{x}_{\mathcal{B}},\mathbf{x}_{\mathcal{W}},\mathbf{x}_{\mathcal{C}}$
$\mathbf{x}$ in body-fixed, world, and camera frame
$\angle\left(\mathbf{a},\mathbf{b}\right)$ signed angle between vectors
(right-hand rule)

Table I: Mathematical notation, nomenclature, and notable symbols.

###  I-B Contributions

We present a novel method for gamma source radiation localization using a
single-detector Compton camera [25]. The method utilizes an onboard position
estimator of a UAV to estimate a real-time hypothesis of the radiation source
position. The proposed method is executed in real-time such that the results
can be used to drive an autonomous search algorithm. Unlike any existing
solution, the system is deployable onboard Micro Aerial Vehicles thanks to the
small form factor of the utilized Compton camera detector. The proposed
approach is intended as an addition to an intensity-based estimation approach
and the proposed method can, on its own, maintain a feasible hypothesis of a
moving radiation source. The used MiniPIX TPX3 detector with Cadmium Telluride
(CdTe) sensor can also detect $\beta$ and heavy-ion radiation and, therefore,
can be used without changes to estimate gradients and contours of radiation
intensity. Even more, particle tracks can be classified [1] to deduce the
nature of the radiation source.

##  II TIMEPIX3 COMPTON CAMERA

The process of recording and imaging ionizing gamma radiation is very
different from the near-visible light spectrum photography [1]. Refractive
optical lenses cannot be made for such energies, and reflective optics
(Kirkpatrick-Baez, Walter, Lobster-eye) are limited to energies not exceeding
$30\text{\,}\mathrm{keV}$ by the reflectivity of the used mirror materials
such as gold and silicon. Pinhole gamma cameras are a viable option. However,
their Field of View (FOV) is narrow, and the aperture shield absorbs a
valuable portion of the incoming particles.

Compton cameras use the Compton scattering effect (see Fig. 2) to reconstruct
a projection of a radioactive source [26]. Typically, a Compton camera
consists of two sensors. The first detector, a _scatterer_ , witnesses the
Compton scattering effect in its sensor by capturing the bi-product of the
reaction — the electron $e^{-}$. The electron is captured immediately upon its
creation. The second sensor is placed behind the first and is dedicated to
capturing the scattered photon $\lambda^{\prime}$ through the photo-electric
effect. When both detectors detect the scattering products, they
simultaneously measure the energy and the 3D position of the particles.
According to Compton, 1923 [27], the relation of energies of the particles to
the scattering angle $\Theta$ is modelled as

$E_{\lambda^{\prime}}=\frac{E_{\lambda}}{1+\left(E_{\lambda}/\left(m_{e}c^{2}\right)\right)\left(1-\cos\Theta\right)},$
(1)

where $E_{\lambda},E_{\lambda^{\prime}}$ are the energies of the incoming and
scattered photon ($\mathrm{J}$), $E_{e^{\scalebox{0.75}[1.0]{$-$}}}$ is the
energy of the electron ($\mathrm{J}$), $m_{e}=9.10938356\times
10^{-31}\,$$\mathrm{kg}$ is the electron rest mass, and $c=299792458$
$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ is the speed of light in vacuum. Since
Compton scattering is a radially symmetrical effect, the reconstruction of
$\Theta$ yields a set of possible incoming directions that forms a cone
surface.

$\lambda$$\lambda^{\prime}$$e^{\scalebox{0.75}[1.0]{$-$}}$$\Theta$

Figure 2:  Compton scattering: an incoming gamma photon $\lambda$ changes its
momentum as it scatters off a scattering center by a scattering angle
$\Theta$. The decrease in momentum creates an electron.

Contrary to recent works [4, 8, 28], we employ only a single sensor to measure
the scattering products [25]. Using cutting-edge innovations in the field of
particle imaging, we employ the Timepix3 detector [29] as part of the MiniPIX
TPX3 event camera with a thick $2\text{\,}\mathrm{mm}$ CdTe sensor
111https://advacam.com/camera/minipix-tpx3. The thick semiconductor sensor is
capable of causing the scattering effects, as well as capturing the products.
Thanks to the detector’s event camera operation mode, both bi-products are
timed with sub-nanosecond accuracy, and the depth of the event within the
sensor is calculated. In contrast with traditional frame-based cameras, event-
driven cameras output a continuous stream of data generated by the active
(hit) pixels. A similar trend emerged in the visible-light camera field [30].
Timepix3 promises a similar impact in the field of aerial radiation detection,
as its essential properties are perfectly suited for dynamically flying MAVs.

The coincidence detection window is set to $86\text{\,}\mathrm{ns}$, which is
the maximum time difference of two coinciding particles being measured at the
opposite side of the $2\text{\,}\mathrm{mm}$ CdTe at $450\text{\,}\mathrm{V}$.
The detection window is specific to this particular sensor and should be
adjusted for different thicknesses, materials, and bias voltages.

###  II-A Reconstructing the set of directions to the source

Figure 3:  An illustration of the Compton scattering effect within the CdTe
sensor on top of the Timepix3 detector. The ionized electron clouds are
produced by absorption $\lambda^{\prime}$ and $e^{-}$ and are gathered towards
the detector pixels by $450\text{\,}\mathrm{V}$ bias potential with a speed of
$\approx$ $23\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{ns}}^{-1}$.

Table I summarizes the mathematical notation used throughout this work and
Fig. 4 provides an illustration of the used coordinate systems. Ionizing
radiation separates electrons (ionization) from the CdTe sensor material in
the location of absorption. The newly created electron cloud is accelerated by
a $450\text{\,}\mathrm{V}$ electric potential towards one facet of the
$2\text{\,}\mathrm{mm}$ CdTe sensor where the charge is measured by individual
pixels of the Timepix3 detector. Figure 3 shows the charge-gathering process
within the sensor. The
$[\lambda^{\prime}_{x},\lambda^{\prime}_{y}]^{\intercal}_{\mathcal{C}}$ and
$[e^{\scalebox{0.75}[1.0]{$-$}}_{x},e^{\scalebox{0.75}[1.0]{$-$}}_{y}]^{\intercal}_{\mathcal{C}}$
coordinates of the events are obtained by calculating the centroids of the
respective pixel tracks (see Fig. 5). The z-coordinates $\lambda^{\prime}_{z}$
and $e^{\scalebox{0.75}[1.0]{$-$}}_{z}$ are unknown. However, it is possible
to estimate the relative difference along the $\mathbf{\hat{c}}_{z}$ using the
measured time-difference of the events. Although the absorption of
$\lambda^{\prime}$ and $e^{\scalebox{0.75}[1.0]{$-$}}$ happens virtually at
the same time222The speed of a gamma photon through CdTe material is near the
speed of light in vacuum and therefore can be neglected., the charge-gathering
causes a delay between the measured electron clouds. The relative z-distance
is obtained as

$\Delta z=2325.6\,(e^{\scalebox{0.75}[1.0]{$-$}}_{t}-\lambda^{\prime}_{t}),$
(2)

where $e^{\scalebox{0.75}[1.0]{$-$}}_{t}$ and $\lambda^{\prime}_{t}$ are the
times of arrival of the electron clouds from the electron and photon,
respectively, and $c_{s}$ = $23.256\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{ns}}^{-1}$ =
$2325.6\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ is the charge-gathering
speed through the CdTe sensor with $450\text{\,}\mathrm{V}$ bias333The charge
gathering speed is obtained empirically by analyzing tracks of muons from a
cosmic radiation.. The absence of a precise absolute z-axis coordinate of the
cone origin within the camera frame is inconsequential. It results only in a
shift of the resulting cone surface, which can be neglected compared to the
real-world distances involved. Only the difference between the scattering
products along $\mathbf{\hat{c}}_{z}$ is significant for estimating the cone
parameters. Therefore, coordinates of the electron and the secondary photon
events in the camera frame $\mathcal{C}$ are considered as

$\displaystyle\mathbf{e}^{\scalebox{0.75}[1.0]{$-$}}_{\mathcal{C}}$
$\displaystyle=[e^{\scalebox{0.75}[1.0]{$-$}}_{x},e^{\scalebox{0.75}[1.0]{$-$}}_{y},\Delta
z]^{\intercal},$ (3) $\displaystyle\bm{\lambda}^{\prime}_{\mathcal{C}}$
$\displaystyle=[\lambda^{\prime}_{x},\lambda^{\prime}_{x},0]^{\intercal}.$ (4)

Figure 4:  Compton camera coordinate frame $\mathcal{C}$ is attached to the
UAV body frame $\mathcal{B}$. The world coordinate frame $\mathcal{W}$ is
static. Transformations between the frames are provided by the UAV control
system [18] through ROS.

In order to reconstruct the scattering angle $\Theta$, we convert the measured
particle energies from ($\mathrm{keV}$) to ($\mathrm{J}$) as

$\displaystyle E_{e^{\scalebox{0.75}[1.0]{$-$}}}$
$\displaystyle=e_{e^{\scalebox{0.75}[1.0]{$-$}}}/6.242\times 10^{18},$ (5)
$\displaystyle E_{\lambda}$ $\displaystyle=e_{\lambda}/6.242\times 10^{18},$
(6)

and calculate the scattering angle using (1) as:

$\displaystyle\Theta=\cos^{-1}\underbrace{\left(1+m_{e}c^{2}\left(\frac{1}{E_{e^{-}}+E_{\lambda}}-\frac{1}{E_{\lambda}}\right)\right)}_{B},$
$\displaystyle\text{for}\ \scalebox{0.75}[1.0]{$-$}1<B<1$ (7)

Finally, the cone parameters are the cone origin
$\mathbf{o}_{\mathcal{C}}=\mathbf{e}^{\scalebox{0.75}[1.0]{$-$}}$, the
direction
$\mathbf{\hat{d}}_{\mathcal{C}}=\mathbf{e}^{\scalebox{0.75}[1.0]{$-$}}-\bm{\lambda}^{\prime}$,
and the inner angle $\Theta$.

|
---|---
photon trackelectron track $\lambda^{\prime}_{e}$ = | $394.22\text{\,}\mathrm{keV}$
---|---
$e^{\scalebox{0.75}[1.0]{$-$}}_{e}$ = | $315.70\text{\,}\mathrm{keV}$
$\Delta t$ = | $20.31\text{\,}\mathrm{ns}$
$\Delta z$ = | $0.47\text{\,}\mathrm{mm}$
$\Theta$ = | $1.13\text{\,}\mathrm{rad}$
Figure 5:  An example of a pair of Compton scattering products captured by the
Timepix3 detector. The events’ times are used together with the particles’
energies to reconstruct the scattering angle, $\Theta$.

##  III RADIATION SOURCE STATE ESTIMATION

Estimating the direction to the source is typically done under assumptions
that do not apply in scenarios involving mobile robots. The camera position is
often fixed during the whole measurement, and all the data are processed at
once after the measurement [26, 28]. In other cases, the camera needs to be
stationary to complete a single measurement set (tens of seconds) before it
can be moved further [4, 8]. However, the estimation should be conducted
continuously without stopping to utilize the advantages of multirotor
helicopters fully. Moreover, in all the previous works involving a mobile
robot, the used camera was a heavy device, which required the use of a large
UGV or UAV. We aim to bring the Compton camera detector to the domain of MAVs.
With appropriate onboard localization and mapping systems, such an MAV can
automatically perform the search for an ionizing radiation source. However,
this requires the MAV to automatically estimate the source’s position onboard
in real-time to incorporate newly acquired data on the fly. Here we present
our initial method developed for a single radiation source localization.

The proposed method relies on the Linear Kalman Filter (LKF) to estimate a
hypothesis and a covariance of the source’s position. However, the LKF is not
natively able to accept the measured cone surfaces as a vector of measurement.
To tackle this issue, we assume that the optimal correction would move the
hypothesis orthogonally in the cone surface direction. Such corrections behave
idempotently under ideal conditions. Applying such a correction requires a
novel approach of using LKF by defining an orthogonal projection to the cone
surface.

###  III-A Orthogonal projection to a cone surface

Projecting the LKF hypothesis $\mathbf{x}_{\mathcal{W}}$ orthogonally to the
cone surface to obtain ${\mathbf{x}^{\prime}}_{\mathcal{W}}$ is illustrated in
Fig. 6. Vector from the origin of the cone $\mathbf{o}_{\mathcal{W}}$ to the
given point $\mathbf{x}_{\mathcal{W}}$ is denoted as
$\mathbf{\hat{u}}=(\mathbf{x}-\mathbf{o})/\|\mathbf{x}-\mathbf{o}\|$. Firstly,
we find the angle between the cone direction $\mathbf{\hat{d}}$ and
$\mathbf{\hat{u}}$ as $\alpha=\angle(\mathbf{\hat{u}},\mathbf{\hat{d}})$.
Then, the angle between the cone surface and the desired orthogonal projection
of $\mathbf{x}$ to the cone surface is $\beta=\alpha-\Theta$. To find the
projection, we first find $\mathbf{x}_{r}$ by rotating $\mathbf{x}$ around the
vector $\mathbf{\hat{d}}\times\mathbf{\hat{u}}$ by the angle $\beta$. The
vector towards the projection is then constructed as
$\mathbf{\hat{v}}=(\mathbf{x}_{r}-\mathbf{o})/\|\mathbf{x}_{r}-\mathbf{o}\|$.
Finally, the orthogonal projection of $\mathbf{x}$ on the cone surface is:

$\mathbf{x}^{\prime}=\begin{cases}\mathbf{o}+\mathbf{\hat{v}}\cos\beta,&\text{if}\
\alpha<\pi/2,\\\ \mathbf{o},&\text{if}\ \pi/2\leq\alpha.\end{cases}$ (8)

Figure 6: Orthogonal projection of $\mathbf{x}$ to the cone surface.

###  III-B Linear Kalman Filter

LKFs filter uses a model with the hypothesis $\mathbf{x}_{[t]}$ and hypothesis
covariance $\mathbf{\Omega}_{[t]}$:

$\mathbf{x}_{[t]}=\mathbf{A}\mathbf{x}_{[t-1]},\leavevmode\nobreak\
\mathbf{A}=\left[\begin{smallmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{smallmatrix}\right],$ (9)

where $\mathbf{A}$ is designed for a stationary 3D target. Measurement noise
covariance with the correction being applied in the direction of the
$\mathbf{\hat{e}}_{3}$ axis is defined as

$\mathbf{R_{e}}=\left[\begin{smallmatrix}r&0&0\\\ 0&10^{9}&0\\\
0&0&10^{9}\end{smallmatrix}\right],$ (10)

where $r$ is used to weigh the correction in the direction to the cone
surface. The rest of the matrix is designed to mitigate corrections for other
directions. This canonical covariance $\mathbf{R}_{e}$ is then rotated every
time to align with the vector $\mathbf{x}^{\prime}-\mathbf{x}$ as
$\mathbf{R}=\mathbf{P}\mathbf{R_{e}}\mathbf{P}^{\intercal}$ where $\mathbf{P}$
is a rotation matrix corresponding to rotation around the axis
$\mathbf{\hat{e}}_{3}\times\left(\mathbf{x}^{\prime}-\mathbf{x}\right)$ by the
angle
$\angle\left(\mathbf{x}^{\prime}-\mathbf{x},\mathbf{\hat{e}}_{1}\right)$. An
obtained measurement is then the projection $\mathbf{x}^{\prime}$ and the
covariance $\mathbf{R}$.

Two variants of the filter were tested: the fully-3D approach presented above,
and a simplified 2D approach in which we assume that the radiation source is
located on the ground. In the latter case, we project the hypothesis onto the
ground plane after each correction.

Initialization of the LKF hypothesis cannot be done using the first arriving
measurement as typical. To tackle this issue, we first gather $N$ cone
surfaces with different points of origin. The gathering can be done using one
of many exploration schemes, e.g., a traditional gradient-based exploration
using particle intensity measurements as heuristics or flying a pre-determined
sweeping path. With the first $N$ cones measured, we initialize the filter by
optimizing for a point $\mathbf{p}$, which has the minimum square distance to
the cones’ surfaces. Firstly, we must define the distance to a cone surface.

###  III-C Distance to a cone surface

Given a point $\mathbf{p}$ such that $\text{if}\
\mathbf{\hat{d}}^{\intercal}_{i}\mathbf{p}<\mathbf{\hat{d}}_{i}^{\intercal}\mathbf{o}_{i}$
($\mathbf{p}$ is situated _behind_ the cone origin), the distance to the cone
surface is trivial:

$d(\mathbf{p},\mathbf{o},\mathbf{\hat{d}},\Theta)=\|\mathbf{o}-\mathbf{p}\|.$
(11)

Otherwise, the distance is expressed as

$d(\mathbf{p},\mathbf{o},\mathbf{\hat{d}},\Theta)=\|\mathbf{p}-\mathbf{o}\|\sin\left(\cos^{-1}\frac{\left(\mathbf{p}-\mathbf{o}\right)\circ\mathbf{\hat{d}}}{\|\mathbf{p}-\mathbf{o}\|}-\Theta\right).$
(12)

###  III-D Finding the closest point to multiple cone surfaces

The state estimator’s initialization is done by finding the closest point to
the first _N_ cone surfaces. The point $\mathbf{p}_{\mathcal{W}}$ is found by
minimizing the squared distance to each cone surface while satisfying
geometric constraints:

$\displaystyle\min_{\mathbf{p}\in\mathbb{R}^{3}}$
$\displaystyle\sum_{i=1}^{N}d(\mathbf{p},\mathbf{o}_{i},\mathbf{\hat{d}}_{i},\Theta_{i})^{2}$
(13) $\displaystyle\text{s.t.}\leavevmode\nobreak\
\mathbf{\hat{d}}^{\intercal}_{i}\mathbf{p}$
$\displaystyle\geq\mathbf{\hat{d}}_{i}^{\intercal}\mathbf{o}_{i},$
$\displaystyle\forall i$ $\displaystyle\in\\{1,\ldots,n\\},$ (14)
$\displaystyle\mathbf{p}^{\intercal}\mathbf{\hat{e}}_{3}$ $\displaystyle=0,$
$\displaystyle\forall i$ $\displaystyle\in\\{1,\ldots,n\\}.$ (15)

The constraint (15) is used optionally in the case of the 2D estimator to
enforce the search for a source near the ground plane.

This optimization problem is a nonlinear least-squares problem, which makes
finding a near-optimum solution impractical and time-consuming. However, we do
not need to solve this problem quickly or in fast succession, and the solution
is used only for initialization. Therefore, we utilize the Sequential Least-
Squares Programming optimization method. The Jacobian of the criterion is
calculated analytically. The problem takes approx. $1\text{\,}\mathrm{s}$ to
solve for $N=5$ on an average onboard computer.

##  IV EXPERIMENTS

The approach was verified in realistic simulations using the Gazebo/ROS
simulator. The Compton camera is modeled by combining real-time Monte-Carlo
ray tracing [1]. Randomized simulations suggest that the proposed method is
capable of initializing and further estimating the 3D position of the
radiation source with no prior knowledge. However, the UAV should conduct
motion orthogonally to the direction of the source. Otherwise, both the LKF
initialization and the LKF correction methods result in singular solutions
that fail to provide a full 3D position of the source. Similarly, the
estimator fails when the UAV is stationary. This can be explained using the
principles of the proposed method. When the UAV is stationary, all the cones
originate from the same point. Therefore, the information extracted from
multiple cones only provide data regarding the direction to the source. When
the UAV conducts a motion directly towards the radiation source, again, the
information extracted from the cones do not provide us with other than
direction data. However, if the UAV moves orthogonally to the direction to the
source, a non-trivial intersection of multiple cones can exist (if they
originate from different 3D positions) and therefore, the 3D position of the
source can be estimated. We can estimate the position of the source even in
situations where traditional gradient-based method might fail, e.g., in urban
environment. The proposed method could be used as an addition to the
intensity-based methods, that can utilize the same Timepix3 detector. The
orthogonal motion can be realized, e.g., by encircling the estimated
hypothesis.

The presented two practical tests were conducted in cooperation of _Advacam_
(developer of the MiniPIX TPX3 detector), the CTU (developer of the proposed
method and the aerial system) and the Czech Metrology Institute (provider of
the sample radiation sources). Only a few experiments were conducted due to
the autumn/winter conditions in our country and the current social distancing
rules, that have made such real-world tests increasingly difficult to
organize.

###  IV-A Simulations

(a)  Figure 7: The simulation model of the DJI F330 UAV.

(a)  (b)  (c)

Figure 8:  Top-down view of samples of simulations with a single radiation
source, moving at the speed of
$1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. First, the UAV searched for
the source by encircling the center of the designated $100\times 100$
$\mathrm{m}$ area. When the source for localized, the UAV then proceeded with
encircling the 2D LKF hypothesis by circles with $10\text{\,}\mathrm{m}$
radius. The method showed good tracking performance for all tested source
energies: $30\text{\,}\mathrm{GBq}$, $10\text{\,}\mathrm{GBq}$,
$3\text{\,}\mathrm{GBq}$, however, the initial search phase took increasingly
longer with a weaker source.

A real-time ray tracing simulation model proposed previously in [1] was used
to evaluate the proposed localization method in conjunction with the realistic
Gazebo/ROS robotics simulator. The DJI F330 quadrotor was modeled (see Fig. 7)
and provided as a part of the MRS UAV system [18]. Development of a research
strategy is not in the scope of this paper, hence we evaluate the localization
method using the following naive search strategy:

  1. 1.

the UAV encircles a designated area until a hypothesis is initialized by Eq.
(13),

  2. 2.

the UAV then proceeds to encircle the hypothesis, updating the LKF with the
incoming measurements,

  3. 3.

if more than 3 cones in a row exhibit an outlier condition, the hypothesis is
reset using Eq. (13).

Simulations show that the proposed method is able to localize a radiation
source even with the naive search strategy. Moreover, the proposed
localization method is even suitable for real time tracking of a moving
[137]Cs radiation source. As shown in Fig. 8, the UAV is able to track a
radiation source moving at the speed of
$1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. In the case of the
$3\text{\,}\mathrm{GBq}$ source, the average rate of detected Compton cones
was $1.7\text{\,}{\mathrm{s}}^{-1}$. In can be observed that a steady tracking
is achieved even with the $3\text{\,}\mathrm{GBq}$ source, despite a lack of a
specialized model in the LKF for a moving target. The simulation software
together with an example scenarios is provided as open source at
https://github.com/rospix.

###  IV-B Unmanned aerial vehicles

(a) (a)

(b) (b)

(c) (c)

Figure 9:  The first experimental flight: the T650 multi-rotor UAV (a) was
equipped with the MiniPIX TPX3 CdTe Compton camera, a 3D LiDAR, and an RGB
camera. SLAM localizes the UAV in 3D, which allows it to build a map of the
surrounding environment (b). The RGB camera (c) serves for visual re-
projection of the source’s position. The estimated position of the source is
shown in comparison with the ground truth position. Videos from the experiment
are available at http://mrs.felk.cvut.cz/radron-icuas.

The first verification of the MiniPIX detector onboard a UAV was conducted
using a general-purpose quadrotor helicopter built upon the Tarot 650 frame
and equipped with Intel NUC-i7 onboard computer. Figure 9a depicts the UAV in
mid-air. The Ouster OS1-16 3D LiDAR was used for ground truth localization and
mapping together with the A-LOAM SLAM [31]. Furthermore, the UAV was also
equipped with a wide-angle RGB camera to aid with data visualization. These
sensors are not required for the proposed approach. The only requirement is
the access to real-time 3D position and orientation state estimate of the UAV.
Such a state estimate is a requirement for any autonomous mission. It can be
accommodated by, e.g., a Global Navigation Satellite System (GNSS) receiver
coupled with an Inertial Measurement Unit (IMU) and a magnetometer, or by a
visual-inertial navigation system such as [32]. Therefore, the Compton camera
system can be employed on much smaller UAV thanks to its small weight and
size. The design and size of the UAV body is more limited by the localization
sensor suite rather the by the $40\text{\,}\mathrm{g}$ radiation detector. The
UAV utilized the open-source _MRS UAV system_ 444http://github.com/ctu-
mrs/mrs_uav_system[18] in ROS for automatic control and state estimation.

In the pursuit of miniaturization, a second experiment was conducted using the
DJI F330 UAV (see Fig. 11a), weighting approx. $1\text{\,}\mathrm{kg}$. The
UAV also conducted an autonomous flight, however, this time using a monocular
visual SLAM, the VINS mono [32]. A wide-angle greyscale camera was used for
localization in real-time.

###  IV-C Compton camera

The Compton camera is built upon the Timepix3 detector [29] with a
$2\text{\,}\mathrm{mm}$ CdTe sensor. CdTe sensors are exceptionally sensitive
to high-energy gamma rays. The sensitivity to high-energy $\gamma$ outweighs
the small size (therefore small effective volume) of the sensor since only
high-energy $\gamma$ can travel hundreds of meters through the air without
being absorbed. This allows sensing the sources from a further distance than
using traditional detectors. The Timepix3 detector with thick CdTe sensor can
detect photons ranging from 5 to $1000\text{\,}\mathrm{keV}$ with relatively
high quantum efficiency especially for Compton scattering events. The
miniaturized Timepix3-based detector MiniPIX TPX3 CdTe with incorporated data
acquisition, calibration and event reconstruction software layers was
developed by Advacam (see Fig. 1). The device does not require any cooling and
can operate at room temperature. The Compton camera payload weighs
$40\text{\,}\mathrm{g}$ and has dimensions of $80\times 20\times 15$
$\mathrm{mm}$. The processing software is implemented in ROS and is available
as open-source.

###  IV-D Experimental environment

These preliminary tests were conducted in $40\times 20$ $\mathrm{m}$ area with
a single [137]Cs source with activity of $187\text{\,}\mathrm{MBq}$
($187\times 10^{6}$ decay events per second). This was a relatively weak
source, sending only $\approx\,$
$30\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{cm}}^{-2}$ photons at
$10\text{\,}\mathrm{m}$ distance.

The first real-world test555The first test video: https://youtu.be/oH4jMMHfGVA
was conducted under autonomous localization, mapping, and control, using the
Tarot 650 UAV. A human operator provided the UAV with control waypoints, and
the UAV followed them with $1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$
speed. The UAV first conducted a sweeping flight while gathering the first 5
cone surfaces for initialization. After the LKF was initialized, it started
encircling the current hypothesis. Such a strategy shows success in
simulations. After successfully localizing the source, the estimator can adapt
to the potentially changing position of the source. Figure 9 depicts the 3D
map of the environment with the hypothesis showed as a red cube surrounded by
a covariance ellipsoid. A yellow cube marks the ground truth position of the
source.

The real-world test666The second test video: https://youtu.be/t6ayhHJ9z5k was
conducted using the miniaturized UAV, the DJI F330. Again, the flight was
conducted autonomously. This time, the point of the experiment was to test an
updated version of the MiniPIX TPX3 detector API, developed specifically for
detecting the Compton coincidences within the embedded hardware of the
detector. Therefore, the UAV was tasked to autonomously encircle the already
known position of the radiation source. Figure 11 depicts the UAV in mid-air,
together with a snapshot of a 3D visualization of the recorded data.

Full-flight data from the experiments, including the RAW sensor data,
localization data, LiDAR point clouds and camera images are available on the
website http://mrs.felk.cvut.cz/radron-icuas in the form of _rosbag_ files.

###  IV-E Results

Figure 10: Plots of convergence of the 3D hypothesis and the 2D hypothesis in
comparison to the ground truth. Note that the vehicle maintained a nearly
constant height above ground during the experiment, hence the slow convergence
in the Z-axis.

Statistics of event occurrences observed during the experiments is shown in
Table II. Figure 10 shows plots of the estimated 3D and 2D position of the
gamma source in comparison to the ground truth. The hypothesis was correctly
estimated near the ground truth after gathering just 15 Compton events. The
experimental data show the estimator system working as expected. Future work
should emphasize on rejecting outliers. Outliers can occur due to the
omnipresent radiation background, even though occurrence of coincidences in
the radiation background is extraordinarily rare. The detector registers only
approx. 10 particles per second due to the radiation background, the chance of
two background particles coinciding within the time frame of approx.
$100\text{\,}\mathrm{ns}$ is extremely low. The largest source of outliers
lies in the uncertainty of origin of the already detected two coinciding
particles. In some cases their type (electron, scattered photon) can be
deduced from the measured energies using, e.g., the knowledge of the Compton
edge for a known isotope. However, if both measured energies lie in the
Compton continuum, the particle type is ambiguous. Moreover, even if the
incoming photon comes from a known direction, back-scattering can make the
measured by-products appear as incoming from an _opposite_ direction. However,
despite the occasional outliers, the estimator performed well despite the
tested radiation source’s low activity. It must be noted, that majority of the
measured signal, that can be seen in the provided videos and in the provided
datasets, is a valid measurement of the incoming gamma photons, not of a
background radiation. A majority of the signal originates from photons
measured by pure photo-absorption, therefore, the signal does not produce
coinciding particles used for Compton reconstruction. Nevertheless, the signal
could be used for traditional intensity-based estimation [17], since the CdTe
Timepix3 detector is fully capable to serve as a dosimeter. The method
proposed in this manuscript aims to provide an additional information that
could aid the traditional state-of-the-art methods, without the need for extra
hardware.

| Event count | Rate [s-1] | Relative share
---|---|---|---
Photoelectric | 6073 | 24.29 | 0.996
Compton | 18 | 0.072 | 0.003
Background | 8 | 0.032 | 0.001
Photoelectric | 5315 | 21.26 | 0.996
Compton | 18 | 0.072 | 0.003
Background | 4 | 0.016 | 0.001
Table II: Events captured by the Compton camera during the first 250 seconds
of a flight. Data from Experiment 1 are shown in the upper section, and
Experiment 2 in the lower section. Events exceeding energy of
$800\text{\,}\mathrm{keV}$ are labeled as background, as these are too
energetic to be released by the [137]Cs source.

(a) (a)

(b) (b)

(c) (c)

Figure 11:  The second experimental flight: the DJI F330 multi-rotor UAV (a),
(b) was equipped with the MiniPIX TPX3 CdTe Compton camera, and a monocular
grey-scale camera. VINS MONO SLAM was used to provide the UAV with 3D
odometry. A snapshot from the video http://mrs.felk.cvut.cz/radron-icuas is
shown, depicting a 3D plot of the UAV’s position, the latest Compton cone and
the estimated position of the radiation source.

##  V LESSONS LEARNED AND FUTURE WORK

The presented results are a pilot work for a project for autonomous
localization and tracking of compact ionizing radiation sources by Micro
Aerial Vehicles. Several significant challenges have been identified that are
yet to be tackled before the solution becomes fully autonomous and reliable.
The single-detector Compton camera should be compared to a multi-sensor stack.
Although the proposed single-detector variant offers a significantly simpler
hardware solution, the angular resolution is worse than with a sensor stack.
Moreover, reconstruction of cone surfaces using (7) should benefit from
outlier rejection using the values of
$E_{\lambda},E_{\lambda^{\prime}},E_{e^{-}}$. Priors of $E_{\lambda}$ can be
estimated in real-time by analyzing the spectrum, a type of incoming
particles. Also, the state estimation hypothesis should be generalized to a
mixture of distributions. In combination with the observed complex 3D map of
the environment, Monte-Carlo methods such as the particle filter should be
explored. If available, prior knowledge of the environment structure, e.g.,
wall material and thickness, should also be utilized by the estimator.

We pursue the use of MAVs equipped with monocular visual SLAM (see Fig. 1).
Autonomous planning methods should simultaneously explore an unknown
environment while maximizing the yield of information gathered by the Compton
camera. This task shares goals with the Defense Advanced Research Projects
Agency (DARPA) SubT challenge777https://www.subtchallenge.com in which the
authors participate. Simulations show that the multi-MAV system benefits
significantly from sharing the measured data between the MAV [33]. We aim to
utilize swarms of MAVs to increase the _baseline_ of the distributed sensor,
which increases the convergence speed of the estimator.

##  VI CONCLUSIONS

We have showcased a novel approach for localization and position estimation of
a gamma radiation source using a miniature single-detector Compton camera on
an Micro Aerial Vehicle. The proposed method of estimating the source’s
position utilizes the particle physics interactions within the CdTe sensor in
conjunction with a Linear Kalman Filter. A novel method of the filter
initialization is proposed, given the Compton camera detector’s unique form of
measurements. Moreover, a novel iterative-based correction method is proposed
to fuse the measurements in the form of a cone surface. The methods were
validated by real-world experiments showing promising results that were, until
now, the domain of heavy-lifting aircraft. Gazebo/ROS simulation plugin for
the Compton camera and the processing software is provided as open-
source888http://github.com/rospix to allow and motivate research on autonomous
systems using the Compton camera detector.

## ACKNOWLEDGMENT

This work was done on behalf of Medipix3 collaboration, was supported by TACR
project no. FW01010317, and by CTU grant no SGS20/174/OHK3/3T/13.

## References

  * [1] T. Baca, M. Jilek, P. Manek, _et al._ , “Timepix Radiation Detector for Autonomous Radiation Localization and Mapping by Micro Unmanned Vehicles,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2019, pp. 1–8.
  * [2] K. Ohno, S. Kawatsuma, T. Okada, _et al._ , “Robotic control vehicle for measuring radiation in Fukushima Daiichi Nuclear Power Plant,” in _2011 IEEE International Symposium on Safety, Security, and Rescue Robotics_. IEEE, 2011, pp. 38–43.
  * [3] K. Nagatani _et al._ , “Emergency response to the nuclear accident at the Fukushima Daiichi Nuclear Power Plants using mobile rescue robots,” _Journal of Field Robotics_ , vol. 30, no. 1, pp. 44–63, 2013.
  * [4] Y. Sato, Y. Terasaka, W. Utsugi, _et al._ , “Radiation imaging using a compact Compton camera mounted on a crawler robot inside reactor buildings of Fukushima Daiichi Nuclear Power Station,” _Journal of Nuclear Science and Technology_ , vol. 56, no. 9-10, pp. 801–808, 2019.
  * [5] H. Kinoshita, R. Tayama, E. Y. Kometani, _et al._ , “Development of new technology for Fukushima Daiichi nuclear power station reconstruction,” _Hitachi Review_ , vol. 63, no. 4, pp. 183–190, 2014.
  * [6] Y. Sanada _et al._ , “Aerial Radiation Monitoring Around the Fukushima Dai-ichi Nuclear Power Plant Using an Unmanned Helicopter,” _Journal of environmental radioactivity_ , vol. 139, pp. 294–299, 2015.
  * [7] J. Towler, B. Krawiec, and K. Kochersberger, “Radiation mapping in post-disaster environments using an autonomous helicopter,” _Remote Sensing_ , vol. 4, no. 7, pp. 1995–2015, 2012.
  * [8] J. Jiang, K. Shimazoe, Y. Nakamura, _et al._ , “A prototype of aerial radiation monitoring system using an unmanned helicopter mounting a GAGG scintillator Compton camera,” _Journal of Nuclear Science and Technology_ , vol. 53, no. 7, pp. 1067–1075, 2016.
  * [9] J. MacFarlane, O. Payton, A. Keatley, _et al._ , “Lightweight aerial vehicles for monitoring, assessment and mapping of radiation anomalies,” _Journal of environmental radioactivity_ , vol. 136, pp. 127–130, 2014.
  * [10] G. Christie, A. Shoemaker, _et al._ , “Radiation search operations using scene understanding with autonomous UAV and UGV,” _Journal of Field Robotics_ , vol. 34, no. 8, pp. 1450–1468, 2017.
  * [11] P. G. Martin, S. Kwong, N. Smith, Y. Yamashiki, O. D. Payton, _et al._ , “3D unmanned aerial vehicle radiation mapping for assessing contaminant distribution and mobility,” _International Journal of Applied Earth Observation and Geoinformation_ , vol. 52, pp. 12–19, 2016.
  * [12] O. Salek, M. Matolin, and L. Gryc, “Mapping of radiation anomalies using UAV mini-airborne gamma-ray spectrometry,” _Journal of environmental radioactivity_ , vol. 182, pp. 101–107, 2018.
  * [13] A. Keatley, P. Martin, K. Hallam, O. Payton, R. Awbery, F. Carvalho, J. Oliveira, L. Silva, M. Malta, and T. Scott, “Source identification of uranium-containing materials at mine legacy sites in Portugal,” _Journal of Environmental Radioactivity_ , vol. 183, pp. 102–111, 2018.
  * [14] P. Martin, O. Payton, J. Fardoulis, _et al._ , “The use of unmanned aerial systems for the mapping of legacy uranium mines,” _Journal of environmental radioactivity_ , vol. 143, pp. 135–140, 2015.
  * [15] J. Han, Y. Xu, L. Di, and Y. Chen, “Low-cost multi-UAV technologies for contour mapping of nuclear radiation field,” _Journal of Intelligent & Robotic Systems_, vol. 70, no. 1-4, pp. 401–410, 2013.
  * [16] A. A. R. Newaz, S. Jeong, H. Lee, H. Ryu, and N. Y. Chong, “UAV-based multiple source localization and contour mapping of radiation fields,” _Robotics and Autonomous Systems_ , vol. 85, pp. 12–25, 2016.
  * [17] F. Mascarich, T. Wilson, C. Papachristos, and K. Alexis, “Radiation source localization in GPS-denied environments using aerial robots,” in _2018 IEEE ICRA_. IEEE, 2018, pp. 6537–6544.
  * [18] T. Baca, M. Petrlik, M. Vrba, V. Spurny, R. Penicka, _et al._ , “The MRS UAV System: Pushing the Frontiers of Reproducible Research, Real-world Deployment, and Education with Autonomous Unmanned Aerial Vehicles,” _Journal of Intelligent & Robotic Systems_, vol. 102, no. 26, pp. 1–28, May 2021.
  * [19] T. Roucek, M. Pecka, P. Cizek, T. Petricek, _et al._ , “DARPA Subterranean Challenge: Multi-robotic exploration of underground environments,” in _IEEE MESAS_ , vol. 11995, 2019, pp. 274–290.
  * [20] M. Petrlik, T. Baca, D. Hert, M. Vrba, T. Krajnik, and M. Saska, “A Robust UAV System for Operations in a Constrained Environment,” _IEEE Robotics and Automation Letters_ , vol. 5, 4 2020.
  * [21] V. Kratky, P. Petracek, _et al._ , “Autonomous reflectance transformation imaging by a team of unmanned aerial vehicles,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 2302–2309, 4 2020.
  * [22] M. Saska, D. Hert, T. Baca, V. Kratky, and T. Nascimento, “Formation Control of Unmanned Micro Aerial Vehicles for Straitened Environments,” _Autonomous Robots_ , pp. 1573–7527, 2020.
  * [23] P. Petracek, V. Walter, T. Baca, and M. Saska, “Bio-Inspired Compact Swarms of Micro Aerial Vehicles without Communication and External Localization,” _Bioinspiration & Biomimetics_, 2020.
  * [24] V. Kratky, P. Petracek, T. Baca, and M. Saska, “An Autonomous Unmanned Aerial Vehicle System for Fast Exploration of Large Complex Indoor Environments,” _Submitted to Journal of Field Robotics_ , 2020, available online: http://mrs.felk.cvut.cz/data/papers/darpa_urban.pdf.
  * [25] D. Turecek, J. Jakubek, E. Trojanova, and L. Sefc, “Single layer Compton camera based on Timepix3 technology,” _Journal of Instrumentation_ , vol. 15, no. 01, p. C01014, 2020.
  * [26] ——, “Compton camera based on Timepix3 technology,” _Journal of Instrumentation_ , vol. 13, no. 11, p. C11022, 2018.
  * [27] A. H. Compton, “A quantum theory of the scattering of X-rays by light elements,” _Physical review_ , vol. 21, no. 5, p. 483, 1923.
  * [28] F. Terzioglu, P. Kuchment, and L. Kunyansky, “Compton camera imaging and the cone transform: a brief overview,” _Inverse Problems_ , vol. 34, no. 5, p. 054002, 2018.
  * [29] T. Poikela, J. Plosila, _et al._ , “Timepix3: a 65K channel hybrid pixel readout chip with simultaneous ToA/ToT and sparse readout,” _Journal of instrumentation_ , vol. 9, no. 05, p. C05013, 2014.
  * [30] B. Kueng _et al._ , “Low-latency visual odometry using event-based feature tracks,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2016\.
  * [31] J. Zhang and S. Singh, “LOAM: Lidar Odometry and Mapping in Real-time,” in _Robotics: Science and Systems_ , vol. 2, no. 9, 2014.
  * [32] T. Qin, P. Li, and S. Shen, “Vins-mono: A robust and versatile monocular visual-inertial state estimator,” _IEEE Transactions on Robotics_ , vol. 34, no. 4, pp. 1004–1020, 2018.
  * [33] P. Stibinger, T. Baca, and M. Saska, “Localization of Ionizing Radiation Sources by Cooperating Micro Aerial Vehicles With Pixel Detectors in Real-Time,” _IEEE Robotics and Automation Letters_ , vol. 5, pp. 3634–3641, 4 2020.

  *[FDNPP]: Fukushima Daiichi Nuclear Power Plant
  *[UAV]: Unmanned Aerial Vehicle
  *[MAVs]: Micro Aerial Vehicle
  *[UAVs]: Unmanned Aerial Vehicle
  *[MAV]: Micro Aerial Vehicle
  *[UGVs]: Unmanned Ground Vehicle
  *[CdTe]: Cadmium Telluride
  *[FOV]: Field of View
  *[ROS]: Robot Operating System
  *[UGV]: Unmanned Ground Vehicle
  *[LKF]: Linear Kalman Filter
  *[LKFs]: Linear Kalman Filter
  *[LiDAR]: Light Detection and Ranging
  *[SLAM]: Simultaneous Localization And Mapping
  *[GNSS]: Global Navigation Satellite System
  *[IMU]: Inertial Measurement Unit
  *[DARPA]: Defense Advanced Research Projects Agency


supplement

# A first-principles study of bilayer 1T’-WTe2/CrI3: A candidate topological
spin filter

Daniel Staros  Department of Chemistry, Brown University, Providence, RI 02912
Brenda Rubenstein  Authors to whom correspondence should be addressed: Brenda
Rubenstein<EMAIL_ADDRESS>Daniel Staros,
<EMAIL_ADDRESS>and Panchapakesan Ganesh<EMAIL_ADDRESS>Department
of Chemistry, Brown University, Providence, RI 02912  Panchapakesan Ganesh
Center for Nanophase Materials Sciences Division, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA

###### Abstract

The ability to manipulate electronic spin channels in 2D materials is crucial
for realizing next-generation spintronics. Spin filters are spintronic
components that polarize spins using external electromagnetic fields or
intrinsic material properties like magnetism. Recently, topological protection
from backscattering has emerged as an enticing feature through which the
robustness of 2D spin filters might be enhanced. In this work, we propose and
then characterize one of the first 2D topological spin filters: bilayer
CrI3/1$T$’-WTe2. To do so, we use a combination of Density Functional Theory
and maximally localized Wannier functions to demonstrate that the bilayer
satisfies the principal criteria for being a topological spin filter; namely
that it is gapless, exhibits charge transfer from WTe2 to CrI3 that renders
the bilayer metallic despite the CrI3 retaining its monolayer ferromagnetism,
and does not retain the topological character of monolayer 1$T$’-WTe2. In
particular, we observe that the atomic magnetic moments on Cr from DFT are
approximately 3.2 $\mu_{B}$/Cr in the bilayer compared to 2.9 $\mu_{B}$/Cr
with small negative ferromagnetic moments induced on the W atoms in
freestanding monolayer CrI3. Subtracting the charge and spin densities of the
constituent monolayers from those of the bilayer further reveals spin-
polarized charge transfer from WTe2 to CrI3. We find that the bilayer is
topologically trivial by showing that its Chern number is zero. Altogether,
this evidence indicates that BL 1$T$’-WTe2/CrI3 is gapless, magnetic, and
topologically trivial, meaning that a terraced WTe2/CrI3 bilayer
heterostructure in which only a portion of a WTe2 monolayer is topped with
CrI3 is a promising candidate for a 2D topological spin filter. Our results
further suggest that 1D chiral edge states may be realized by stacking
strongly hybridized ferromagnetic monolayers, like CrI3, atop 2D $non$magnetic
Weyl semimetals like 1$T$’-WTe2.

## Introduction

As signs continue to suggest that Moore’s Law has plateaued after many
decades, researchers have begun to seek new routes to designing faster,
smaller, more energy-efficient, and more versatile electronic devices. The key
to realizing such devices will be discovering, characterizing, and designing
novel nanoscale quantum electronic components whose many electronic degrees of
freedom, including their electron spinWolf  _et al._ (2001); Ahn  _et al._
(2020) and momenta,Zhang  _et al._ (2021, 2022a) can be manipulated to enable
faster, more energy-efficient operations on denser data.

Along these lines, nanoscale spintronics have been hailed as extremely
promising routes towards denser data storage and potentially faster and more
efficient reading and writing. Unlike conventional electronics which harness
the charge of an electron, spintronic materials store information in
electrons’ two possible spin states,Wolf  _et al._ (2001) which can be
manipulated more rapidly and with less energy than electrons’ charges.El-
Ghazalya  _et al._ (2020) Spintronic devices are also less volatile than
conventional electronic devices because they can preserve their spin even in
the absence of electric power.Wolf  _et al._ (2001); Barla  _et al._ (2021)
Moreover, one of the primary advantages of spintronic devices is that they can
be readily integrated into modern CMOS-based circuits.Barla  _et al._ (2021)

Since the birth of spintronics with the discovery of the giant
magnetoresistive effect,Wolf  _et al._ (2001) the world of spintronic device
components has expanded to include various spintronic analogues to traditional
resistors and transistors, as well as new components unique to controlling
spin currents like spin filters and spin injectors.Hirohataa  _et al._ (2020)
Many of these components take advantage of the properties of magnetic
materials, in which spins are already selectively ordered. For example, two
spintronic analogues to traditional resistors, spin valves and magnetic tunnel
junctions (MTJs), typically consist of two ferromagnetic layers separated by
an insulating layer.Barla  _et al._ (2021); Hirohataa  _et al._ (2020) Varying
the magnetic orientation of one of the magnetic layers and keeping the other
fixed allows the resistance to spin currents to be changed by taking advantage
of spin-selective quantum tunneling as in MTJs, or the giant magnetoresistive
or spin-transfer torque effects as in spin valves. Perhaps the most
fundamental spintronic device component, however, is that which enables the
generation of spin current in the first place: the spin filter.

Figure 1: (a) Edge conductance around ML WTe2. (b) Proposed spin-polarized
edge conductance/1D chiral edge state of the CrI3/WTe2 bilayer studied here.

Spin filters are devices that generate spin-polarized currents from
unpolarized electric currents by selectively transmitting electrons with a
particular spin and blocking those with the opposite spin. In general, such
devices have most often taken advantage of the inherent spin polarization in
ferromagnetic or multiferroic materials,Inomata  _et al._ (2008); Bea  _et
al._ (2008); Zhang  _et al._ (2022b) spin-selective quantum tunneling using
barrier materials with different spin-dependent transmission
probabilities,Massarotti  _et al._ (2015) or spin orbit coupling in Rashba-
type spin filters to achieve this.Koga, Nitta, and Takayanagi (2002);
Cummings, Akis, and Ferry (2006) Spin filters can additionally be designed by
placing ferromagnetic insulators in close proximity to superconducting
junctions since the presence of an out-of-plane magnetic field can break time-
reversal symmetry, causing the electrons to form spin-polarized currents
through the spin-transfer torque effect.Hasan and Kane (2010); Mak, Shan, and
Ralph (2019); Yoshida (2019)

Recently, the concept of topological spin filters has been put forth as one
promising option for improving the robustness of spin filters at higher
temperatures by taking advantage of quantum anomalous Hall conductance, which
is topologically protected from backscattering and could minimize dissipation
as a result. By extension, edge-conductance dominated by one spin channel
amounts to nearly dissipationless conductance of one direction of quantum
spin, also known as topological spin filtering, which can manifest as a chiral
edge state along the edge of a partially exposed topologically trivial bilayer
and nontrivial monolayer (see Figure 1).Howard  _et al._ (2021) The magnetic
Weyl semimetal Co3Sn2S2 has recently been discussed as a potential avenue
towards realizing such higher temperature chiral conducting edge states,Hsieh
_et al._ (2016); Howard  _et al._ (2021) which could in principle also become
spin-polarized. Along the same vein, a topological spin filter may be
constructed by placing a ferromagnet near topologically nontrivial 2D
materials such as 1$T^{\prime}$ transition metal dichalcogenides,Qian  _et
al._ (2014) which could also give rise not only to spin-polarized currents,
but spin-polarized helical edge modes that are topologically protected from
backscattering.

While many such three-dimensional spin filters have been proposed, two-
dimensional materials and their heterostructures possess a larger design space
advantageous for engineering new spintronic devices.Sierra  _et al._ (2021)
Two-dimensional heterostructures can be designed to exhibit a wide array of
emergent properties by mixing and matching the properties of their constituent
monolayers,Novoselov  _et al._ (2016); Pham  _et al._ (2022) twisting them,He
_et al._ (2021); Xu  _et al._ (2022); Morissette  _et al._ (2023); Fumega and
Lado (2023) straining them,Dai, Liu, and Zhang (2019) or placing them in
proximity to electric and/or magnetic fields.Maximenko  _et al._ (2022);
Marmodoro  _et al._ (2022); Tao  _et al._ (2022); Iordanidou  _et al._ (2023)
Despite this, researchers have only recently made significant strides towards
truly 2D spin filters which promise to be smaller, more tunable, and ideally
more efficient than their 3D counterparts. Graphene is one 2D material that
originally garnered spintronic interest when it was predicted to exhibit
nearly perfect spin filtering when interfaced with only a couple of layers of
a ferromagnetic metal.Karpan  _et al._ (2007) However, experimental attempts
to realize such a graphene-based spin filter fell short, initially showing
tunnel magnetoresistance ratios of 0.4% for graphene/NiFe, with additional
attempts increasing this ratio to no more than 5%.Ahn  _et al._ (2020);
Piquemal-Banci  _et al._ (2017); Cobas  _et al._ (2012) More recently, the
pivotal discovery of giant magnetoresistance in bilayer CrI3 in 2018 has
reinvigorated the search for better 2D spin filters based on atomically thin
magnets,Song  _et al._ (2018); Wang  _et al._ (2018); Klein  _et al._ (2018);
Zhu  _et al._ (2021) and a slew of inspired studies have since been published
that take advantage of their properties.Zhu  _et al._ (2021); Luo  _et al._
(2023); Li  _et al._ (2023a, b) These examples suggest that, with the right
combination of 2D monolayers, 2D spin filters, and even 2D topological spin
filters, should also be within reach.

Notably, there are few studies which consider the proximity effects of 2D
magnets stacked atop a monolayer of $1T^{\prime}$-WTe2, which is the only MX2
monolayer that exists in the $1T^{\prime}$ phase in its ground state and the
only such member which is topological as a freestanding monolayer.Tang  _et
al._ (2017) Until fairly recently, the only such example consisted of one
layer of $1T^{\prime}$-WTe2 interfaced with one layer of permalloy (Ni80Fe20)
to form a film with several-nm thickness that exhibited out-of-plane magnetic
anistropy .Feng  _et al._ (2018); Li  _et al._ (2018) More recently,
proximity-induced magnetic order was observed in monolayer $1T^{\prime}$-WTe2
placed onto antiferromagnetic trilayer CrI3, where edge conductance jumps were
observed upon switching of CrI3’s magnetic state.Chang (2020); Zhao  _et al._
(2020) Most recently, proximity-induced half-metallicity and complete spin-
polarization was predicted in bilayer $1T^{\prime}$-WTe2/CrBr3 and attributed
to strong orbital hybridization and charge transfer at the interface of the
heterostructure.Bora  _et al._ (2023) Nonetheless, to the best of our
knowledge, no investigation of the topological properties of a
1$T^{\prime}$-WTe2/CrX3 bilayer has yet been performed, let alone with the
goal of realizing a new type of topological spin filter. Taken together, these
discoveries point towards bilayer 1$T^{\prime}$-WTe2/CrI3 as a strong
potential candidate for topological spin filtering which could leverage the
perfect spin filtering of a 1$T^{\prime}$-WTe2/CrX3 hetrostructure in
proximity to the dissipationless edge states of $1T^{\prime}$-WTe2.

Thus, in this manuscript, we use ab initio simulations to identify bilayer
1$T^{\prime}$-WTe2/CrI3 as a promising candidate for a 2D topological spin
filter. 1$T^{\prime}$-WTe2 is a nonmagnetic Weyl semimetal which exhibits
topological edge conductance in its monolayer form,Qian  _et al._ (2014) while
monolayer CrI3 is a ferromagnetic Mott insulator.Huang  _et al._ (2017) One
can thus imagine that, by placing these two materials in proximity, the CrI3’s
magnetism could potentially polarize WTe2’s edge conductance, forming a
topological spin filter. To determine whether this is in fact the case, we
predict the band structures, topological invariants, and interlayer charge and
magnetization density transfer for BL CrI3/1$T^{\prime}$-WTe2 with and without
spin-orbit coupling. In so doing, we unequivocally demonstrate that the
proximity of CrI3 to WTe2 foremost results in strong interlayer coupling
between the two layers, spin-polarization on the WTe2, and an overall trivial
BL topology. These considerations, taken together with the metallic nature of
the bilayer and previous evidence for spin-polarized helical edge modes in
monolayer $1T^{\prime}$-WTe2, provide convincing evidence for the possibility
of realizing chiral edge states at the interface of a terraced
$1T^{\prime}$-WTe2/CrI3 bilayer. Specifically, our results imply that electric
current injected into the metallic bilayer portion of the terraced
heterostructure would become spin-polarized before transferring to WTe2 and
exiting via spin-polarized edge conductance in chiral edge states around the
monolayer WTe2 portion; this terraced, strained $1T^{\prime}$-WTe2/CrI3
bilayer is then likely a strong candidate for a highly robust, ultra-thin spin
filter with 1D chiral edge states.

## Methods

In order to determine the ground electronic states of the constituent
monolayers and 1$T^{\prime}$-WTe2/CrI3 bilayer, as well as the charge transfer
present in the bilayer, we calculated ground electronic states and charge and
spin densities using self-consistent Density Functional Theory. We then
Wannierized the DFT orbitals into a form allowing for the calculation of Chern
numbers, which determine whether a material is topologically nontrivial. As
the different rhombohedral angles of monolayer $R3$ CrI3 and
1$T^{\prime}$-WTe2 do not lend themselves to the simple construction of a
commensurate supercell without the introduction of different strain to both
layers, the layers were strained slightly by hand such that they could share a
common cell. We therefore begin with a description of our methodology by first
discussing the determination of the appropriate strain and interlayer distance
of the bilayer before going into greater detail about the DFT calculations and
computation of the Chern numbers.

### Bilayer Supercell Construction

To construct our bilayer, we used a highly-accurate DMC-optimized monolayer
CrI3 structure containing 8 atoms in its unit cell and exhibiting triclinic
($R3$) symmetry.Staros  _et al._ (2022) A monolayer $1T$’-WTe2 structure was
obtained from the Materials Project website.Jain  _et al._ (2013) During the
simulation of these monolayers, more than 20 Å of vacuum was added to both
structures to prevent spurious self interactions.

Figure 2: (Left) Collinear PBE+$U$ spin density difference (yellow = positive)
$s_{BL}-s_{CrI3}-s_{WTe2}$ with an isosurface level value of 0.0025 shows
positive magnetization on the Cr atoms. The WTe2 layer is magnetized, but the
isosurface level does not resolve it. (Right) Noncollinear PBE+$U$
magnetization density difference (yellow = positive)
$m_{BL}-m_{CrI3}-m_{WTe2}$ with an isosurface level value of 0.0025 shows
enhanced magnetization relative to the collinear non-SOC bilayer.

In stacking the layers, special consideration was given to how to align them
since their monolayer structures are incommensurate. In particular, the WTe2
cell was rotated such that the original lattice constants for the monolayer
cells were strained as little as possible. The result of this process was a
bilayer with a lattice constant of 7.01 Å, which means that our DMC-optimized
ML CrI3 structure is stretched by 2.5% relative to the monolayer, and WTe2 (a
= 3.505 Å) is stretched by 0.8% relative to the experimental bulk
1$T^{\prime}$-WTe2 value of 3.477 Å.Mar, Jobic, and Ibers (1992) Additionally,
our monolayer 1$T^{\prime}$-WTe2 lattice constant is close to the value of
3.502 Å previously obtained using DFT structural relaxation.Xiang  _et al._
(2016)

### Density Functional Simulations

All simulations of structural and electronic properties were performed using
DFT as implemented within the Quantum ESPRESSO package.Giannozzi  _et al._
(2009, 2017) The PBE and PBE+$U$ functionalsPerdew, Burke, and Ernzerhof
(1996) with $U=2$ eV on the chromium atoms were selected to model these
materials because previous studies demonstrated that a trial wavefunction
utilizing a Hubbard $U$ value of 2-3 eV minimizes the fixed-node error in DMC
calculations of CrI3.Staros  _et al._ (2022); Ichibha  _et al._ (2021) Our
calculations used norm-conserving, scalar-relativistic Cr and relativistic I
pseudopotentials and recently developed spin–orbit relativistic effective W
and Te pseudopotentials.Guangming  _et al._ (2022) We employed a Monkhorst-
Pack k-point mesh with dimensions $10\times 10\times 1$ and a plane wave
energy cutoff of 300 Ry.

### Calculating Topological Invariants

To calculate the Chern numbers for monolayer WTe2 and the bilayer
heterostructure, subsets of DFT single-particle Bloch functions were
bijectively rotated onto sets of maximally-localized Wannier functions
(MLWF’s)Marzari and Vanderbilt (1997) starting from selected columns of the
density matrix from DFT via the SCDM-k method.Damle, Lin, and Ying (2017) This
mapping was performed for the isolated set of 31 monolayer 1$T^{\prime}$-WTe2
bands ranging from -10 eV below to 0.6 eV above the monolayer Fermi level, and
for entangled sets of bilayer 1$T^{\prime}$-WTe2/CrI3 bands as detailed in the
Supplementary Information. All of the obtained MLWF’s were well-localized and
replicated the DFT band structure well over the span of bands involved in
calculating topological invariants. Next, the hopping terms and correction
terms for the lattice vectors of the hopping terms output by Wannier90 were
used as input to the tight-binding model for the Z2Pack software for
calculating topological invariants. Z2Pack is capable of calculating the
evolution of hybrid Wannier charge centers across the surface defined by an
explicit Hamiltonian $H(\textbf{k})$, a tight-binding model, or an explicit
first-principles calculation.Soluyanov and Vanderbilt (2011); G.  _et al._
(2017) Thus, with our tight-binding model as input, we used Z2Pack to
calculate the hybrid Wannier center evolution of the MLWF’s corresponding to
the bands up to and including the two orbitals involved in WTe2’s spin-orbit-
induced gap openingQian  _et al._ (2014) on a small k-space sphere with a
radius of 0.001 centered at the $\Gamma$-point of the first Brillouin zone.
All of the Z2Pack calculations passed the line and surface convergence checks
to within the default tolerances of the Z2Pack software.Soluyanov and
Vanderbilt (2011); G.  _et al._ (2017)

## Results and Discussion

### Interlayer Charge Transfer and Magnetic Induction

As a first step toward understanding the physics of our bilayer, we began by
examining how the layers influence each other’s electronic structure. To do
so, we analyzed the difference between the bilayer and individual monolayers’
charge and spin densities. If proximity effects are truly at play, we would
expect to see significant differences in their bilayer charge and spin
densities relative to the separate monolayer densities. That said, when a
monolayer of CrI3 is stacked on a monolayer of 1$T^{\prime}$-WTe2, the DFT-
predicted charge density difference between the bilayer and monolayers,
$\rho_{BL}-\rho_{CrI3}-\rho_{WTe2}$, clearly shows charge accumulating near
the CrI3/WTe2 interface as ithe charge is drawn downwards (see Figure 3). This
suggests that CrI3 is a charge acceptor and WTe2 is a charge donor in the
bilayer. We see this charge transfer effect both with (non-collinear
calculations) and without (collinear calculations) including SOC in our DFT
calculations (see Figure S3), suggesting that it is a robust feature of the
bilayer.

Interestingly, charge also accumulates in between and around the bilayer with
charge accumulating near the CrI3 within the vdW interface, and withdrawn from
the portion of CrI3 which is facing away from the interface. The existence of
significant charge density within the bilayer gap confirms the strong
hybridization between the iodine and tellurium atoms, which is also reflected
in the significant atomic overlaps of all four atomic species in the partial
densities of states or PDOS (see Figure 5). Additionally, the metallicity of
the bilayer is reflected in the PDOS occupations of all four atomic species at
and near the Fermi level, indicating that this hybridization causes CrI3 to
lose its Mott insulating nature when it is interfaced with WTe2. The Lowdin
charges of the constituent monolayers and bilayer are tabulated in
Supplementary Tables I and III to quantify the extent of charge transfer in
this bilayer. Summing the individual atomic charges of the monolayers and the
bilayer yields an electron transfer of approximately 0.06 $e$ per primitive
bilayer cell.

Figure 3: (Left) Collinear PBE charge density of monolayer WTe2 in the absence
of CrI3 with an isosurface level value of 0.8. (Right) Collinear PBE+$U$
charge density difference between the bilayer and individual monolayers,
$\rho_{BL}-\rho_{CrI3}-\rho_{WTe2}$, with an isosurface level value of 0.0025.
Charge transfer from WTe2 to CrI3 is evident. Yellow indicates a positive
charge density and light blue indicates a negative charge density. Chromium
atoms are colored dark blue and iodine atoms are purple, while tungsten atoms
are grey and tellurine atoms are beige.

Lastly, we consider the spin polarization that accompanies these charge
transfer effects by evaluating the collinear spin density difference
$s_{BL}-s_{CrI3}-s_{WTe2}$ and noncollinear magnetization density difference
$m_{BL}-m_{CrI3}-m_{WTe2}$ between the bilayer and individual monolayers. Both
quantities point to spin-polarized charge transfer, resulting in induced
magnetization on the WTe2 as exhibited by the positive magnetization density
(yellow) on WTe2 in Figure 2. As can be observed from Figure 2, spin-orbit
coupling enhances this polarization over the collinear case because spin-orbit
coupling increases the time-reversal symmetry-breaking in CrI3, which in turn
influences the WTe2. As before, we tabulate the DFT-predicted atomic magnetic
moments of ML CrI3 and bilayer 1$T^{\prime}$-WTe2/CrI3 in Supplementary Tables
SII and SIV to quantify the spin-polarized charge transfer. The atomic
magnetic moments on the Cr atoms are enhanced by WTe2, with net negative
atomic moments of -0.1 $\mu_{B}$ induced on the W atoms, which are of the same
order as the iodine moments in the freestanding monolayer. The induced spin
polarization is thus stronger than the electron transfer.

### Topological Properties of ML 1$T^{\prime}$-WTe2 and BL WTe2/CrI3

To assess the potential for topology that can give rise to chiral edge states
in our bilayer, we first examined the band structures of the individual
monolayers and combined bilayer structure. Previous modeling has shown that
1$T^{\prime}$-WTe2 possesses a band crossing below the Fermi level, which
gives rise to its nontrivial topology.Tang  _et al._ (2017) As a first step,
we thus determined the band structures of our strained WTe2 monolayer with and
without spin-orbit coupling (SOC) to verify that our slight distortion does
not change the bands significantly.

Figure 4: (Top) Collinear PBE band structure (black) and partial density of
states (PDOS) of 0.8%-strained 1$T^{\prime}$-WTe2 (relative to the bulk
1$T^{\prime}$-WTe2 monolayer) with the interpolated MLWF band structure
overlain (light blue). The corresponding PDOS of the W $d$-orbitals and Te
$p$-orbitals are shown in grey and red on the right. (Bottom) Noncollinear PBE
band structure with spin-orbit coupling.

Indeed, as shown in Figure 4, a band crossing occurs mid-way between the
$\Gamma$ and X high-symmetry points as is also observed in simulations of
pristine WTe2 without spin-orbit coupling. Additionally, the density of states
in this region has slightly more W $d$-orbital character than Te $p$-orbital
character and is consistent with the previous assignment of these bands to the
hybridized W 5$d_{xz}$ and 5$d_{z^{2}}$ orbitals.Tang  _et al._ (2017) These
band structures confirm that the strain applied to the monolayer did not alter
its topology. Integration of the Berry curvature over the fiber bundle of
MLWFs up to the band circled in Figure 4 yields a Chern number of 0, meaning
the topology is trivial as expected, since spin-orbit coupling was turned off.
The same analysis of the noncollinear WTe2 monolayer when spin-orbit coupling
is included (bottom of Figure 4) exhibits a band gap within the 5$d_{xz}$ and
5$d_{z^{2}}$ bands about halfway between the $\Gamma$ and X high-symmetry
points. Integration of the Berry curvature over the fiber bundle of these
bands yields a Chern number of 1, verifying that this monolayer is
topologically nontrivial.Tang  _et al._ (2017)

Figure 5: (Top) Collinear PBE+$U$ band structure (black) and partial density
of states (PDOS) of bilayer 1$T^{\prime}$-WTe2/CrI3 with the interpolated MLWF
band structures overlain (light blue). The PDOS of the W $d$-orbitals, Te
$p$-orbitals and $d$-orbitals, and I $p$-orbitals are shown respectively in
grey, red, dark blue, and purple. (Bottom) Noncollinear PBE+$U$ band structure
with spin-orbit coupling.

Next, we calculated the band structure, partial density of states, and Chern
number for the bilayer composed of ML 1$T^{\prime}$-WTe2 and ML CrI3 (Figure
5) with and without SOC. This bilayer loses the Mott insulating behavior of
CrI3, with the bilayer exhibiting a metallic band structure and finite density
of states at and around the Fermi level. The bilayer also maintains monolayer
1$T^{\prime}$-WTe2’s band crossing below the Fermi level, which opens when
spin-orbit coupling is introduced as is visible in the upper panel of Figure
5. The integration of the Berry curvature over the fiber bundle of MLWF’s up
to and including the W 5$d_{xz}$ and 5$d_{z^{2}}$ orbitals yields a Chern
number of 0, when spin-orbit coupling is not included in the DFT calculation,
while the same process for the spin-orbit-coupled MLWF’s also yields a Chern
number of 0, indicating that this bilayer loses the topological character of
monolayer 1$T^{\prime}$-WTe2 when CrI3 is stacked on top of it. Since the
constituent 1$T^{\prime}$-WTe2 layer is itself topological, this observation
lends itself to the possibility of realizing a chiral edge state in a terraced
1$T^{\prime}$-WTe2/CrI3 bilayer due to the change in the Chern number from 0
in the bilayer to 1 as the CrI3 layer ends but the topological
1$T^{\prime}$-WTe2 layer continues. The step acts as the boundary where the
chiral edge state should exist due to the bulk-boundary correspondence. In
addition, the strong hybridization of the WTe2 with the ferromagnetic CrI3
should break the degeneracy of the two chiral edge states of the ML WTe2,
causing the step edge of CrI3 to host only a spin-polarized chiral edge state.

## Conclusions

In conclusion, we have used a combination of Density Functional Theory and
maximally-localized Wannier function-based tight-binding models to demonstrate
that bilayer 1$T^{\prime}$-WTe2/CrI3 is a topologically nontrivial metallic
material which exhibits enhanced Cr magnetic moments and spin-polarized charge
transfer from WTe2 to CrI3. Most notably, the topologically nontrivial
monolayer WTe2 becomes trivial when a monolayer of CrI3 is placed on top of
it, and this is reflected in a Chern number of 1 for monolayer
1$T^{\prime}$-WTe2 and 0 for the bilayer. When taken in conjunction with the
spin-polarized and metallic nature of this bilayer, our findings suggest that
it may be a possible candidate for realizing a chiral conducting edge state
via the construction of a terraced bilayer composed of a sheet of monolayer
1$T^{\prime}$-WTe2 that is partly covered by a monolayer of CrI3.
Additionally, this is the first evidence for this type of behavior in a system
composed of a $non$magnetic Weyl semimetal placed next to an atomically thin
magnet, thus expanding the concept of such terraced chiral edge states beyond
magnetic Weyl semimetal materials.

## Acknowledgements

This study was funded by The U.S. Department of Energy through the Office of
Science Graduate Student Research (SCGSR) Program. The funder played no role
in study design, data collection, analysis and interpretation of data, or the
writing of this manuscript. This research was partly conducted as part of a
user project at the Center for Nanophase Materials Sciences (CNMS), which is a
US Department of Energy, Office of Science User Facility at Oak Ridge National
Laboratory. This project used resources of the Oak Ridge Leadership Computing
Facility and the National Energy Research Scientific Computing Center, both
U.S. Department of Energy Office of Science User Facilities operated under
Contract No. DE-AC0205CH11231. P.G. initiated and co-supervised the research.
Simulations and analysis were performed by D.S.; the manuscript was prepared
by D.S. and B.R. All authors read, edited, and approved the final manuscript.
All authors declare no financial or non-financial competing interests.

## Data availability

The datasets used and/or analyzed during the current study are available from
the corresponding author upon request.

## References

## References

  * Wolf  _et al._ (2001) S. Wolf _et al._ , “Spintronics: A spin-based electronics vision for the future,” Science 294 (2001).
  * Ahn  _et al._ (2020) E. Ahn _et al._ , “2d materials for spintronic devices,” npj 2D Materials and Applications 4 (2020).
  * Zhang  _et al._ (2021) S. Zhang _et al._ , “Two-dimensional heterostructures and their device applications: progress, challenges and opportunities—review,” Journal of Physics D: Applied Physics 54 (2021).
  * Zhang  _et al._ (2022a) T. Zhang _et al._ , “2D spontaneous valley polarization from inversion symmetric single-layer lattices,” npj Computational Materials 8 (2022a).
  * El-Ghazalya  _et al._ (2020) A. El-Ghazalya _et al._ , “Progress towards ultrafast spintronics applications,” Journal of Magnetism and Magnetic Materials 502 (2020).
  * Barla  _et al._ (2021) P. Barla _et al._ , “Spintronic devices: a promising alternative to cmos devicess,” Journal of Computational Electronics 20 (2021).
  * Hirohataa  _et al._ (2020) A. Hirohataa _et al._ , “Review on spintronics: Principles and device applications,” Journal of Magnetism and Magnetic Materials 509 (2020).
  * Inomata  _et al._ (2008) K. Inomata _et al._ , “Highly spin-polarized materials and devices for spintronics,” Science and Technology of Advanced Materials 9 (2008).
  * Bea  _et al._ (2008) H. Bea _et al._ , “Spintronics with multiferroics,” Journal of Physics: Condensed Matter 20 (2008).
  * Zhang  _et al._ (2022b) L. Zhang _et al._ , “Electrical switching of spin-polarized current in multiferroic tunneling junctions,” npj Computational Materials 8 (2022b).
  * Massarotti  _et al._ (2015) D. Massarotti _et al._ , “Macroscopic quantum tunnelling in spin filter ferromagnetic josephson junctions,” Nature Communications 6 (2015).
  * Koga, Nitta, and Takayanagi (2002) T. Koga, J. Nitta, and H. Takayanagi, “Spin-filter device based on the rashba effect using a nonmagnetic resonant tunneling diode,” Physical Review Letters 88 (2002).
  * Cummings, Akis, and Ferry (2006) A. Cummings, R. Akis, and D. Ferry, “Electron spin filter based on rashba spin-orbit coupling,” Applied Physics Letters 89 (2006).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, “Colloquium : Topological insulators,” Reviews of Modern Physics 82 (2010).
  * Mak, Shan, and Ralph (2019) K. F. Mak, J. Shan, and D. Ralph, “Probing and controlling magnetic states in 2D layered magnetic materials,” Nature Reviews Physics 1 (2019).
  * Yoshida (2019) N. Yoshida, “Theory of the spin-filtering effect in ferromagnet/ferromagnetic insulator/superconductor junctions,” Journal of Physics Communications 3 (2019).
  * Howard  _et al._ (2021) S. Howard _et al._ , “Evidence for one-dimensional chiral edge states in a magnetic weyl semimetal co3sn2s2,” Nature Communications 12 (2021).
  * Hsieh  _et al._ (2016) T. Hsieh _et al._ , “Bulk topological proximity effect,” Physical Review Letters 116 (2016).
  * Qian  _et al._ (2014) X. Qian _et al._ , “Quantum spin hall effect in two-dimensional transition metal dichalcogenides,” Science 346 (2014).
  * Sierra  _et al._ (2021) J. F. Sierra _et al._ , “Van der Waals heterostructures for spintronics and opto-spintronics,” Nature Nanotechnology 16 (2021).
  * Novoselov  _et al._ (2016) K. S. Novoselov _et al._ , “2d materials and van der waals heterostructures,” Science 2 (2016).
  * Pham  _et al._ (2022) P. Pham _et al._ , “2d heterostructures for ubiquitous electronics and optoelectronics: Principles, opportunities, and challenges,” Chemical Reviews 122 (2022).
  * He  _et al._ (2021) F. He _et al._ , “Moiré Patterns in 2D Materials: A Review,” ACS Nano 15 (2021).
  * Xu  _et al._ (2022) Y. Xu _et al._ , “Coexisting ferromagnetic–antiferromagnetic state in twisted bilayer CrI3,” Nature Nanotechnology 17 (2022).
  * Morissette  _et al._ (2023) E. Morissette _et al._ , “Dirac revivals drive a resonance response in twisted bilayer graphene,” Nature Physics (2023).
  * Fumega and Lado (2023) A. Fumega and J. Lado, “Moiré-driven multiferroic order in twisted CrCl ${}_{\textrm{3}}$ , CrBr ${}_{\textrm{3}}$ and CrI ${}_{\textrm{3}}$ bilayers,” 2D Materials 10 (2023).
  * Dai, Liu, and Zhang (2019) Z. Dai, L. Liu, and Z. Zhang, “Strain Engineering of 2D Materials: Issues and Opportunities at the Interface,” Advanced Materials 31 (2019).
  * Maximenko  _et al._ (2022) Y. Maximenko _et al._ , “Nanoscale studies of electric field effects on monolayer 1t’-wte2,” npj Quantum Materials 7 (2022).
  * Marmodoro  _et al._ (2022) A. Marmodoro _et al._ , “Electric field control of magnons in magnetic thin films: Ab initio predictions for two-dimensional metallic heterostructures,” Physical Review B 105 (2022).
  * Tao  _et al._ (2022) W. Tao _et al._ , “Multiband superconductivity in strongly hybridized 1T -WTe2/NbSe2 heterostructures,” Physical Review B. 105 (2022).
  * Iordanidou  _et al._ (2023) K. Iordanidou _et al._ , “Electric Field and Strain Tuning of 2D Semiconductor van der Waals Heterostructures for Tunnel Field-Effect Transistors,” ACS Applied Materials & Interfaces 15 (2023).
  * Karpan  _et al._ (2007) V. Karpan _et al._ , “Graphite and graphene as perfect spin filters,” Physical Review Letters 99 (2007).
  * Piquemal-Banci  _et al._ (2017) M. Piquemal-Banci _et al._ , “2d-mtjs: introducing 2d materials in magnetic tunnel junctions.” Journal of Physics D: Applied Physics 50 (2017).
  * Cobas  _et al._ (2012) E. Cobas _et al._ , “Graphene as a tunnel barrier: graphene-based magnetic tunnel junctions.” Nano Letters 12 (2012).
  * Song  _et al._ (2018) T. Song _et al._ , “Giant tunneling magnetoresistance in spin-filter van der waals heterostructures,” Science 360 (2018).
  * Wang  _et al._ (2018) Z. Wang _et al._ , “Very large tunneling magnetoresistance in layered magnetic semiconductor cri3,” Nature Communications 9 (2018).
  * Klein  _et al._ (2018) D. Klein _et al._ , “Probing magnetism in 2d van der waals crystalline insulators via electron tunneling,” Science 360 (2018).
  * Zhu  _et al._ (2021) Y. Zhu _et al._ , “Giant tunneling magnetoresistance in van der waals magnetic tunnel junctions formed by interlayer antiferromagnetic bilayer cobr2,” Physical Review B. 103 (2021).
  * Luo  _et al._ (2023) B. Luo _et al._ , “Stacking order, perfect spin polarization, and giant magnetoresistance in zigzag graphene/h-bn heterobilayer nanoribbons,” Physical Review Applied 19 (2023).
  * Li  _et al._ (2023a) Q. Li _et al._ , “Stacking- and strain-dependent magnetism in janus crste bilayer,” Applied Physics Letters 122 (2023a).
  * Li  _et al._ (2023b) Y. Li _et al._ , “Nonvolatile electrical control of spin polarization in the 2d bipolar magnetic semiconductor vsef,” npj Computational Materials 9 (2023b).
  * Tang  _et al._ (2017) S. Tang _et al._ , “Quantum spin hall state in monolayer 1t’-wte2,” Nature Physics 13 (2017).
  * Feng  _et al._ (2018) J. Feng _et al._ , “Excellent spin-filtering and giant tunneling magnetoresistance in a dual-electrode van der waals magnetic tunnel junction based on ferromagnetic crse2,” Nature Physics 611 (2018).
  * Li  _et al._ (2018) P. Li _et al._ , “Spin-momentum locking and spin-orbit torques in magnetic nano-heterojunctions composed of Weyl semimetal WTe2,” Nature Communications 9 (2018).
  * Chang (2020) C. Chang, “Marriage of topology and magnetism,” Nature Materials 19 (2020).
  * Zhao  _et al._ (2020) W. Zhao _et al._ , “Magnetic proximity and nonreciprocal current switching in a monolayer WTe2 helical edge,” Nature Materials 19 (2020).
  * Bora  _et al._ (2023) M. Bora _et al._ , “Adaptive half-metallicity via magnetic proximity in an electrically sensitive 1T’-WTe2/CrBr3 vdW heterostructure,” Applied Surface Science 623 (2023).
  * Huang  _et al._ (2017) B. Huang _et al._ , “Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit,” Nature 546 (2017).
  * Staros  _et al._ (2022) D. Staros _et al._ , “A combined first principles study of the structural, magnetic, and phonon properties of monolayer cri3,” Journal of Chemical Physics 156 (2022).
  * Jain  _et al._ (2013) A. Jain _et al._ , “Commentary: The materials project: A materials genome approach to accelerating materials innovation,” APL Materials 1 (2013).
  * Mar, Jobic, and Ibers (1992) A. Mar, S. Jobic, and J. Ibers, “Metal-metal vs tellurium-tellurium bonding in wte2 and its ternary variants tairte4 and nbirte4,” Journal of the American Chemical Society 114 (1992).
  * Xiang  _et al._ (2016) H. Xiang _et al._ , “Quantum spin hall insulator phase in monolayer wte2 by uniaxial strain,” AIP Advances 6 (2016).
  * Giannozzi  _et al._ (2009) P. Giannozzi _et al._ , “Quantum espresso: a modular and open-source software project for quantum simulations of materials,” Journal of Physics: Condensed Matter 21 (2009).
  * Giannozzi  _et al._ (2017) P. Giannozzi _et al._ , “Advanced capabilities for materials modelling with quantum espresso,” Journal of Physics: Condensed Matter 29 (2017).
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Physical Review Letters 77 (1996).
  * Ichibha  _et al._ (2021) T. Ichibha _et al._ , “Cri3 revisited with a many-body ab initio theoretical approach,” Physical Review Materials 5 (2021).
  * Guangming  _et al._ (2022) W. Guangming _et al._ , “A new generation of effective core potentials from correlated and spin–orbit calculations: Selected heavy elements,” The Journal of Chemical Physics 157 (2022).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, “Maximally localized generalized Wannier functions for composite energy bands,” Physical Review B 56 (1997).
  * Damle, Lin, and Ying (2017) A. Damle, L. Lin, and L. Ying, “Scdm-k: Localized orbitals for solids via selected columns of the density matrix,” Journal of Computational Physics 334 (2017).
  * Soluyanov and Vanderbilt (2011) A. Soluyanov and D. Vanderbilt, “Computing topological invariants without inversion symmetry,” Physical Review B 83 (2011).
  * G.  _et al._ (2017) D. G. _et al._ , “Z2pack: Numerical implementation of hybrid wannier centers for identifying topological materials,” Physical Review B 95 (2017).