# Gaussian Process for Trajectories

Kien Nguyen
Computer Science Department
University of Southern California
Los Angeles, CA 90007
<EMAIL_ADDRESS>
&John Krumm
Microsoft Research
Microsoft Corporation
Redmond, WA 98052
<EMAIL_ADDRESS>
&Cyrus Shahabi
Computer Science Department
University of Southern California
Los Angeles, CA 90007
<EMAIL_ADDRESS>

###### Abstract

The Gaussian process is a powerful and flexible technique for interpolating
spatiotemporal data, especially with its ability to capture complex trends and
uncertainty from the input signal. This chapter describes Gaussian processes
as an interpolation technique for geospatial trajectories. A Gaussian process
models measurements of a trajectory as coming from a multidimensional
Gaussian, and it produces for each timestamp a Gaussian distribution as a
prediction. We discuss elements that need to be considered when applying
Gaussian process to trajectories, common choices for those elements, and
provide a concrete example of implementing a Gaussian process.

_K_ eywords Gaussian process $\cdot$ trajectory interpolation $\cdot$
spatiotemporal data

##  1 Introduction

The availability of devices with location tracking capability has helped
generate a tremendous amount of trajectory data of humans, animals, vehicles,
and drones, which is valuable in various applications [1]. However,
trajectories may have missing locations between their measurements and
uncertainty in the measurements depending on how they were captured.
Interpolating missing locations between measurements, or making predictions
about future locations, is often necessary in order to derive high utility
from trajectories. This work presents Gaussian Process (GP) as a powerful and
flexible technique for trajectory interpolation and prediction.

A trajectory $S$ is defined as a sequence of time-ordered noisy location
measurements
$S=\\{\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots,\boldsymbol{x}_{|S|}\\}$,
where each measurement or data point $\boldsymbol{x}$ (bold symbol) includes
$\boldsymbol{x}.\textit{lon}$ and $\boldsymbol{x}.\textit{lat}$ as the
longitude and latitude, $\boldsymbol{x}.\textit{t}$ as the timestamp, and
$\boldsymbol{x}.\sigma$ as the accuracy or uncertainty. Previous work [2]
showed that it is reasonable to assume Gaussian noise for location
measurements such as GPS points. Therefore, $\boldsymbol{x}.\sigma$ can be
considered as the standard deviation of independent Gaussian noise of
longitude and latitude. Figure 1 shows an example of a trajectory with its
measurements shown as black crosses, uncertainty as the blue area, and some
predictions made by a GP, which is explained later.

Figure 1: Illustration of a Gaussian process trained from data (black cross)
making predictions (red cross).

GP has been studied and applied extensively [3]. There have been several works
that applied GP for different types of trajectories such as in robotics [4, 5,
6, 7, 8], video surveillance [9, 10], batch trajectories processing [11],
sensor trajectories in active sensing [12], travel-time prediction [13], and
human motion [14, 15]. This work describes how to apply GP for trajectory
interpolation and prediction, especially GPS-based trajectories.

##  2 Gaussian Process

A Gaussian Process (GP) is a generalization of the Gaussian probability
distribution. While a probability distribution describes scalar or vector
random variables, a stochastic process describes properties of functions,
i.e., instead of describing the probability of generating a point as a
Gaussian distribution, a GP describes the probability of generating a
function. If one can loosely consider a function as a very long vector, where
each entry in the vector indicates the function value $f(\textit{t})$ at a
specific input t (which is a timestamp in our case), then a GP implies that
any subset of entries in that long vector is distributed according to a
multidimensional Gaussian. We build a separate GP for latitude and another one
for longitude, assuming that their values change independently.

As an example, one wants to estimate locations of a user at a set of
timestamps
$\mathcal{T}=\\{\textit{t}_{1},\textit{t}_{2},\dots,\textit{t}_{|\mathcal{T}|}\\}$.
Thus, the input of the function $f(\textit{t})$ that a GP (for either latitude
or longitude) describes is a timestamp $\textit{t}\in\mathcal{T}$, and the
output is the mean and variance of the Gaussian distribution that describes
the distribution of $f(\textit{t})$. A GP implies that the set
$\\{f(\textit{t})|\textit{t}\in\mathcal{T}\\}$ is distributed according to a
$|\mathcal{T}|$-dimensional Gaussian. When $\mathcal{T}$ includes only one
timestamp t, the longitude (or latitude) of the estimated location is
distributed according to a univariate Gaussian distribution. Figure 1 shows an
example of a GP trained from data points (represented as black crosses) and
making predictions (represented as red crosses). The location axis can either
be longitude, latitude, or their transformed values such as converted to
another spatial coordinate system.

A GP is specified by its mean function $m(\textit{t})$ and covariance function
$k(\textit{t},\textit{t}^{\prime})$. A deterministic mean function
$m(\textit{t})$ is used to specify the mean of the multidimensional Gaussian,
defining the mean of the estimated locations. In the literature, GPs are often
studied with a zero mean function, which means $m(\textit{t})=0$ for all t,
because a deterministic mean function can be incorporated after the kernel
function is trained from training data.

The covariance matrix of the multidimensional Gaussian defines how points at
different values of the independent variable (time in our case) correlate with
each other. In general, timestamps that are closer would correlate more and
those that are further apart would correlate less, which means the correlation
decreases to zero as $|\textit{t}-\textit{t}^{\prime}|$ becomes larger. A GP
models such correlation using a scalar covariance kernel/function
$k(\textit{t},\textit{t}^{\prime})$. A GP can be very flexible by using and
combining different kernels, which will be discussed later. A kernel often has
some parameters that need to be trained from data, such as pertaining to the
uncertainty of the input data or how strong a correlation should be given the
time difference $|\textit{t}-\textit{t}^{\prime}|$.

When using a GP for trajectory interpolation and prediction, one would need to
specify its mean function $m(\textit{t})$ and covariance function
$k(\textit{t},\textit{t}^{\prime})$, train its parameters from data, then make
predictions. We discuss each of these elements in the next section.

##  3 Gaussian Process Elements

In this section, we describe the elements needed to implement a GP for
trajectories, including data preparation, mean function, covariance function
(or kernel), training, and inference. Our implementation uses the GPFlow
library [16], which is a package for building Gaussian process models in
Python, but the concepts also apply to other libraries.

###  3.1 Data Preparation

We first describe how to prepare data from a trajectory
$S=\\{\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots,\boldsymbol{x}_{|S|}\\}$ for
GP training and inference. As mentioned, two separate GPs are created for
longitude and latitude predictions. Input data for training each GP includes
the sequence $\\{(\textit{t}_{1},v_{1}),\dots,(\textit{t}_{|S|},v_{|S|})\\}$
where $\textit{t}_{i}=\boldsymbol{x}_{i}.\textit{t}$ and coordinate $v_{i}$
can be either $\boldsymbol{x}_{i}.\textit{lon}$ or
$\boldsymbol{x}_{i}.\textit{lat}$, and the measurement uncertainty
$\sigma_{m}=\boldsymbol{x}.\sigma$. From the training sequence, we create a
feature vector $\boldsymbol{X}=[\textit{t}_{1},\dots,\textit{t}_{|S|}]^{T}$
and label vector $\boldsymbol{y}=[v_{1},\dots,v_{|S|}]^{T}$. Similarly, for
interpolating or predicting locations of the user at timestamps in a timestamp
set $\mathcal{T}$, a feature vector
$\boldsymbol{X}_{*}=[{\textit{t}|\textit{t}\in\mathcal{T}}]^{T}$ is created.

In addition to $\boldsymbol{X},\boldsymbol{y},\boldsymbol{X}_{*}$ and
$\sigma_{m}$, one may also use a prior uncertainty $\sigma_{0}$, which
represents the prior information of the uncertainty of locations of the users
in terms of the standard deviation of a distribution. For example,
$\sigma_{0}$ can be the standard deviation of a large Gaussian distribution
covering an entire city area, which would indicate a high degree of prior
uncertainty.

###  3.2 Mean Function

As mentioned, the mean function $m(\textit{t})$ is a deterministic function
defining the mean of the estimated locations. A GP makes predictions based
mainly on parameters learned from training data. When a prediction timestamp
is far from training timestamps, the prediction tends to return to the mean
function.

A GP is often assumed to have zero mean, i.e.,
$m(\textit{t})=0,\forall\textit{t}$. However, when there is a clear pattern,
it can be beneficial to provide a more realistic mean to the model. For
example, a trajectory may have a linear trend or we may believe future
locations are likely to follow a linear extrapolation of two most recent
measurements.

Using the GPFlow library, common mean functions can be specified using
gpflow.mean_functions. For example, a constant mean function $m(\textit{t})=c$
with gpflow.mean_functions.Constant, a linear mean function
$m(\textit{t})=At+b$ with gpflow.mean_functions.Linear, or a custom mean
function can also be specified.

###  3.3 Kernel/Covariance Function

The kernel or covariance function $k(\textit{t},\textit{t}^{\prime})$ is the
core element that helps a GP capture trends in data and gives a GP flexibility
and efficiency. Different kernels capture different kinds of trends, and
kernels can also be combined together for more complex models. In this
section, we describe some common kernels and some common types of
combinations. Examples of these kernels are shown in Figure 2, and examples of
these combinations are shown in Figure 3. More detail about kernels can be
found in [3] and [17].

Note that parameters of a kernel can be set to be trained from data or fixed.
For example, if we know the measurement uncertainty $\sigma_{m}$, we can set
it to be fixed; otherwise, it can be learned from the data. If using the
GPFlow library, this can be set using the set_trainable method.

####  3.3.1 Common Kernel Types

##### Constant Kernel

The constant kernel outputs constant predictions, i.e., $f(\textit{t})=c$ with
$c\sim\mathcal{N}(0,\sigma^{2})$. The formula is

$k_{Constant}(\textit{t},\textit{t}^{\prime})=\sigma^{2}$ (1)

The constant kernel is often used to modify the mean of a GP or to scale the
magnitude of the other factor (kernel) of a kernel. The GPFlow library
provides this kernel by gpflow.kernels.Constant with parameter variance.

##### White Kernel

The white kernel explains/produces the noise-component of a measurement. The
formula is

$k_{White}(\textit{t},\textit{t}^{\prime})=\sigma^{2}\delta(\textit{t},\textit{t}^{\prime})$
(2)

where $\delta(\textit{t},\textit{t}^{\prime})=1$ when
$\textit{t}=\textit{t}^{\prime}$; otherwise,
$\delta(\textit{t},\textit{t}^{\prime})=0$. GPFlow provides this kernel with
gpflow.kernels.White with parameter variance. Since we know that the
measurement uncertainty is $\sigma_{m}$, we can set the variance value of this
kernel to $\sigma_{m}^{2}$ (measurement noise) and set this parameter as non-
trainable.

##### Linear Kernel

The linear kernel outputs linear predictions, i.e.,
$f(\textit{t})=c\textit{t}$ with $c\sim\mathcal{N}(0,\sigma^{2})$. The formula
is

$k_{Linear}(\textit{t},\textit{t}^{\prime})=\sigma^{2}\textit{t}\textit{t}^{\prime}$
(3)

The linear kernel is often used in combination with other kernels to capture
linear trends in data. This kernel is a non-stationary kernel, which means it
depends on the absolute locations of the two inputs. We will see later some
stationary kernels, which only depend on the relative positions of its two
inputs. The GPFlow library provides this kernel by gpflow.kernels.Linear with
parameter variance.

##### Squared Exponential (SE)/Radial Basis Function (RBF) Kernel

The SE/RBF kernel is often used as the default kernel for a GP. This is a
stationary and smooth kernel with the formula

$k_{SE}(\textit{t},\textit{t}^{\prime})=\sigma^{2}\exp{\Big{(}-\frac{(\textit{t}-\textit{t}^{\prime})^{2}}{2l^{2}}\Big{)}}$
(4)

The length scale $l$ indicates how much one data point affects another and is
often trained from data. The variance $\sigma^{2}$ is trained from data but
one can set this variance to $\sigma_{0}^{2}$ to provide prior information
about the variance of predictions. GPFlow provides this kernel by
gpflow.kernels.SquaredExponential with parameters variance and length_scale.
Also, instead of providing scalar values, prior distributions can also be used
as the prior for the length scale $l$ and variance $\sigma^{2}$.

##### Rational Quadratic (RQ) Kernel

Another common kernel that is often used as the default kernel for a GP is the
rational quadratic (RQ) kernel with the formula

$k_{RQ}(\textit{t},\textit{t}^{\prime})=\sigma^{2}\Big{(}1+\frac{(\textit{t}-\textit{t}^{\prime})^{2}}{2\alpha
l^{2}}\Big{)}^{-\alpha}$ (5)

The RQ kernel is equivalent to adding together many SE kernels with different
length scales $l$. It converges to the SE kernel when
$\alpha\rightarrow\infty$. Parameter $\alpha$ controls the relative weighting
of large-scale and small-scale variations. The GPFlow library provides this
kernel by gpflow.kernels.RationalQuadratic with parameters variance,
length_scale, and alpha.

##### Matérn Kernel

The Matérn kernel is another common kernel. It is a stationary kernel and a
generalization of the SQ kernel with an additional parameter $\nu$ controlling
the smoothness. It converges to the SE kernel when $\nu\rightarrow\infty$. The
formula for the Matérn kernel is rather complicated. However, there are three
common values of $\nu$, which are $\frac{1}{2},\frac{3}{2},\frac{5}{2}$, that
produce three common kernels called Matern12, Matern32, and Matern52. Their
formulas are

$\displaystyle k_{Matern12}(\textit{t},\textit{t}^{\prime})$
$\displaystyle=\sigma^{2}\exp{\Big{(}-\frac{1}{l}|\textit{t}-\textit{t}^{\prime}|\Big{)}}$
(6) $\displaystyle k_{Matern32}(\textit{t},\textit{t}^{\prime})$
$\displaystyle=\sigma^{2}\Big{(}1+\frac{\sqrt{3}}{l}|\textit{t}-\textit{t}^{\prime}|\Big{)}\exp{\Big{(}-\frac{\sqrt{3}}{l}|\textit{t}-\textit{t}^{\prime}|\Big{)}}$
(7) $\displaystyle k_{Matern52}(\textit{t},\textit{t}^{\prime})$
$\displaystyle=\sigma^{2}\Big{(}1+\frac{\sqrt{5}}{l}|\textit{t}-\textit{t}^{\prime}|+\frac{5}{3l}(\textit{t}-\textit{t}^{\prime})^{2}\Big{)}\exp{\Big{(}-\frac{\sqrt{5}}{l}|\textit{t}-\textit{t}^{\prime}|\Big{)}}$
(8)

Using the GPFlow library, these kernels are provided by
gpflow.kernels.Matern12, gpflow.kernels.Matern32 and gpflow.kernels.Matern52
with parameters variance and length_scale.

Figure 2: Examples of functions drawn from GPs with different kernels. Each
figure shows 3 functions drawing from the same GP with the blue line
indicating the mean of the GP and the blue area indicating the variance of the
GP. Predictions made from these kernels would have a similar shape.

####  3.3.2 Common Types of Kernel Combination

Kernels can be combined together to create more complex kernels, helping to
capture more complex trends in the data. For example, we may want to capture
both a wiggling movement with an SE kernel and measurement uncertainty with a
white kernel. Two common types of combining kernels are summing and
multiplying. Examples of these combinations are shown in Figure 3.

##### Summing Kernels

Summing two kernels $k_{1}$ and $k_{2}$ means summing their output:

$k_{Sum}(\textit{t},\textit{t}^{\prime})=k_{1}(\textit{t},\textit{t}^{\prime})+k_{2}(\textit{t},\textit{t}^{\prime})$
(9)

Roughly speaking, summing kernels is similar to an OR operation, which means
the output of $k_{Sum}$ would be higher if the output of either $k_{1}$ or
$k_{2}$ is higher. An example of summing kernels is to capture both a wiggling
movement with a SE kernel and measurement uncertainty with a white kernel:

$k(\textit{t},\textit{t}^{\prime})=k_{SE}(\textit{t},\textit{t}^{\prime})+k_{White}(\textit{t},\textit{t}^{\prime})=\sigma_{SE}^{2}\exp{\Big{(}-\frac{(\textit{t}-\textit{t}^{\prime})^{2}}{2l^{2}}\Big{)}}+\sigma_{White}^{2}\delta(\textit{t},\textit{t}^{\prime})$
(10)

An illustration of this combination is shown in Figure 3.

##### Multiplying Kernels

Multiplying two kernels $k_{1}$ and $k_{2}$ means multiplying their output:

$k_{Mul}(\textit{t},\textit{t}^{\prime})=k_{1}(\textit{t},\textit{t}^{\prime})\times
k_{2}(\textit{t},\textit{t}^{\prime})$ (11)

Roughly speaking, multiplying kernels is similar to an AND operation, which
means the output of $k_{Sum}$ would be higher if the output of both $k_{1}$
and $k_{2}$ is higher. An example of multiplying kernels is multiplying two
linear kernels to create a quadratic kernel or more to create higher-degree
kernels. An illustration of this combination is shown in Figure 3.

Figure 3: Examples of functions drawn from GPs with different types of kernel
combination. Each figure shows 3 functions drawing from the same GP with the
blue line indicating the mean of the GP and the blue area indicating the
variance of the GP. Predictions made from these kernels would have a similar
shape.

###  3.4 Training and Inference

After constructing a mean function and a kernel and setting trainable
parameters appropriately, training and inference processes are often
straightforward using provided methods in GP libraries. Using GPFlow, for
training, one can use existing optimizers in the package gpflow.optimizers
such as the Scipy optimizer with input $\boldsymbol{X},\boldsymbol{y}$. Once
trained, we can use the trained model to make predictions with predict_y,
which returns the mean and variance of predictions for new data points
$\boldsymbol{X}_{*}$.

Note that the training process minimizes the negative log marginal likelihood.
Since likelihood specifies the probability density of the
observations/measurements given the parameters, for trajectory data, we can
use Gaussian likelihood, which is the default in the GPFlow library. The
likelihood variance indicates the uncertainty learned from data. For
trajectory data, uncertainty in the movement of measurements can be captured
by the model with likelihood variance, and the uncertainty of each measurement
can be captured by a white kernel.

##  4 Gaussian Process Example

In this section, we present an example of using GP for trajectories. In this
example, we use given (training) timestamps
$\boldsymbol{X}=[11,12,\dots,20]^{T}$, given locations
$\boldsymbol{y}=[1,10,30,45,40,40,50,40,35,50]^{T}$, and seek computed values
at timestamps $\boldsymbol{X}_{*}=[21,22,\dots,24]^{T}$ (illustrated in Figure
1). For reproducibility, we use the GPFlow library and set the TensorFlow
random seed to 1. GPFlow also provides a utility method print_summary that
prints values of the parameters of a model.

(a) Not trained

(b) Trained

(c) Trained with White noise

(d) Same as (4(c)), but with fixed White noise

(e) Same as (4(d)), but with small likelihood variance

(f) Same as (4(e)), but with linear mean function

Figure 4: Example of GPs trained from the same training data (as black cross)
with different configuration of kernel and mean function.

The first model $m_{0}$ has a zero mean function and a Matern32 kernel with
default parameters, which is lengthscales $l=1$ and variance $\sigma^{2}=1$.
Figure 4(a) shows the predictions of this model before we trained its
parameters with $\boldsymbol{X}$. We can see that the predictions do not
appropriately capture measurements and give a quick return to 0.

After we train the model with the data to get a trained model $m_{1}$, we have
$l=7.35$, $\sigma^{2}=1346.36$, and likelihood variance (which indicates model
variance) $\sigma_{likelihood}^{2}=36.4$. Figure 4(b) shows the predictions of
this model after we trained with $\boldsymbol{X}$. It is clear that $m_{1}$
can capture the data and trend in the data significantly better than $m_{0}$.

Next, we consider model $m_{2}$ as a summed combination of a Matern32 kernel
and a white kernel, all with default parameters indicating default priors
(e.g., prior for measurement uncertainty captured by the white kernel is 1).
This also means all parameters can be trained from data. After we train the
model with the data, we have $l=7.35$, $\sigma_{Matern32}^{2}=1346.36$,
$\sigma_{White}^{2}=18.23$, and $\sigma_{likelihood}^{2}=18.23$. Figure 4(c)
shows that the model $m_{2}$ makes similar predictions compared to model
$m_{1}$. However, the white kernel captured some of the measurement noise.

Next, we consider model $m_{3}$, which similar to $m_{2}$, but we fix the
variance of the white kernel to variance $\sigma_{White}^{2}=4$, which means
we assume that each measurement has uncertainty $\sigma^{2}=4$. Fixing the
variance means the training process will not change this parameter. After we
train the model with the data, we have $l=7.35$,
$\sigma_{Matern32}^{2}=1346.36$, $\sigma_{White}^{2}=4$, and
$\sigma_{likelihood}^{2}=32.464$. Since we fixed the variance of the white
kernel, the likelihood variance was adjusted to capture more of the
uncertainty. Figure 4(d) shows that the model $m_{3}$ make similar predictions
compared to model $m_{2}$.

Next, we consider model $m_{4}$ similar to $m_{3}$ but in addition to fixing
the variance of the white kernel to $4$, we also set the likelihood variance
to a very small initial value $\sigma_{likelihood}^{2}=0.0001$, which
indicates that we are confident of our knowledge of measurement uncertainty.
After we train the model with the data, we have $l=4.11$,
$\sigma_{Matern32}^{2}=1328.03$, $\sigma_{White}^{2}=4$, and
$\sigma_{likelihood}^{2}=0.00015$. Figure 4(e) shows that, compared to model
$m_{3}$, the model $m_{4}$ make predictions closer to the training
measurements. Figure 1 shows predictions made by $m_{4}$ for
$\boldsymbol{X}_{*}$.

In the final example, we consider model $m_{5}$ similar to $m_{4}$ but using a
linear mean function starting with a guess $y=x+1$. After we train the model
with the data, we have $l=1.37$, $\sigma_{Matern32}^{2}=107.24$,
$\sigma_{White}^{2}=4$, $\sigma_{likelihood}^{2}=0.0001$, and a mean function
$y=4.84x-42.14$. The training process learned the linear line going through
the measurements and used that as the mean function. The model then made
predictions closer to that mean function and with smaller variance, as shown
in Figure 4(f).

##  5 Discussion

Gaussian process is a powerful and flexible technique that is capable of
capturing trends and uncertainty of location measurements. A GP can use a
simple or complex mean function, or even use output of other models, such as a
deep neural network model [18], as its mean function. A GP can use standard
kernels, combined kernels, or other complex kernels. These capabilities help
GP incorporate knowledge from other domains for better interpolation and
prediction.

There are some issues that need to be taken into account when one considers
using GPs for trajectories. One issue is that the length scale of a kernel,
which roughly indicates how far in time a measurement affects other
predictions, tends to be determined by the non-smooth region in the data. For
trajectories, these regions are often where we have dense measurements with
many direction changes. These regions can make the length scale become too
small, resulting in a very sharp change in variance. It can be beneficial to
create different GPs when the trajectory sampling rate changes significantly.
Other issues can be the assumption of independence between longitude and
latitude changes, choosing appropriate kernels, or computational complexity.
However, in general, GP is a useful technique for trajectory interpolation and
prediction.

## References

  * [1] Yu Zheng.  Trajectory data mining: an overview.  TIST, 6(3):1–41, 2015.
  * [2] Frank Van Diggelen.  Gnss accuracy: Lies, damn lies, and statistics.  GPS world, 18(1):26–33, 2007.
  * [3] Christopher K Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning.  MIT press Cambridge, MA, 2006.
  * [4] Gregory Cox, George Kachergis, and Richard Shiffrin.  Gaussian process regression for trajectory analysis.  In Proceedings of the annual meeting of the cognitive science society, volume 34, 2012.
  * [5] Tim D Barfoot, Chi Hay Tong, and Simo Särkkä.  Batch continuous-time trajectory estimation as exactly sparse gaussian process regression.  In Robotics: Science and Systems, volume 10. Citeseer, 2014.
  * [6] Lukas Hewing, Elena Arcari, Lukas P Fröhlich, and Melanie N Zeilinger.  On simulation and trajectory prediction with gaussian process dynamics.  In Learning for Dynamics and Control, pages 424–434. PMLR, 2020\.
  * [7] Gang Cao, Edmund M-K Lai, and Fakhrul Alam.  Gaussian process model predictive control of an unmanned quadrotor.  Journal of Intelligent & Robotic Systems, 88(1):147–162, 2017\.
  * [8] Elnaz Jahani Heravi and Sohrab Khanmohammadi.  Long term trajectory prediction of moving objects using gaussian process.  In 2011 First International Conference on Robot, Vision and Signal Processing, pages 228–232. IEEE, 2011.
  * [9] Kihwan Kim, Dongryeol Lee, and Irfan Essa.  Gaussian process regression flow for analysis of motion trajectories.  In 2011 International Conference on Computer Vision, pages 1164–1171. IEEE, 2011.
  * [10] David Ellis, Eric Sommerlade, and Ian Reid.  Modelling pedestrian trajectory patterns with gaussian processes.  In 2009 IEEE 12th International Conference on Computer Vision Workshops, ICCV Workshops, pages 1229–1234. IEEE, 2009.
  * [11] Mattias Tiger and Fredrik Heintz.  Online sparse gaussian process regression for trajectory modeling.  In 2015 18th International Conference on Information Fusion (Fusion), pages 782–791. IEEE, 2015.
  * [12] Jerome Le Ny and George J Pappas.  On trajectory optimization for active sensing in gaussian process models.  In Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference, pages 6286–6292. IEEE, 2009.
  * [13] Tsuyoshi Idé and Sei Kato.  Travel-time prediction using gaussian process regression: A trajectory-based approach.  In Proceedings of the 2009 SIAM International Conference on Data Mining, pages 1185–1196. SIAM, 2009.
  * [14] Jisoo Hong, Changmook Chun, Seung-Jong Kim, and Frank C Park.  Gaussian process trajectory learning and synthesis of individualized gait motions.  IEEE Transactions on Neural Systems and Rehabilitation Engineering, 27(6):1236–1245, 2019.
  * [15] Jack M Wang, David J Fleet, and Aaron Hertzmann.  Gaussian process dynamical models for human motion.  IEEE transactions on pattern analysis and machine intelligence, 30(2):283–298, 2007.
  * [16] Mark van der Wilk, Vincent Dutordoir, ST John, Artem Artemev, Vincent Adam, and James Hensman.  A framework for interdomain and multioutput Gaussian processes.  arXiv:2003.01115, 2020.
  * [17] David Duvenaud.  The kernel cookbook: Advice on covariance functions.  https://www.cs.toronto.edu/ duvenaud/cookbook/.  Last accessed: September 22, 2021.
  * [18] Vincent Fortuin and Gunnar Rätsch.  Deep mean functions for meta-learning in gaussian processes.  arXiv preprint arXiv:1901.08098, 2019.


#  Automated Attribute Extraction from Legal Proceedings

Subinay Adhikary<EMAIL_ADDRESS>Indian Institute of Science
Education and Research KolkataKolkataIndia ,  Sagnik Das
<EMAIL_ADDRESS>Indian Institute of Science Education and Research
KolkataKolkataIndia ,  Sagnik Saha<EMAIL_ADDRESS>Indian Institute
of Science Education and Research KolkataKolkataIndia ,  Procheta Sen
<EMAIL_ADDRESS>University of LiverpoolP.O. Box
1212LiverpoolUnited kingdom ,  Dwaipayan Roy<EMAIL_ADDRESS>Indian Institute of Science Education and Research KolkataKolkataIndia and
Kripabandhu Ghosh<EMAIL_ADDRESS>Indian Institute of Science
Education and Research KolkataKolkataIndia

(2018)

###### Abstract.

The escalating number of pending cases is a growing concern worldwide. Recent
advancements in digitization have opened up possibilities for leveraging
artificial intelligence (AI) tools in the processing of legal documents.
Adopting a structured representation for legal documents, as opposed to a mere
bag-of-words flat text representation, can significantly enhance processing
capabilities. With the aim of achieving this objective, we put forward a set
of diverse attributes for criminal case proceedings. We use a state-of-the-art
sequence labeling framework to automatically extract attributes from the legal
documents. Moreover, we demonstrate the efficacy of the extracted attributes
in a downstream task, namely legal judgment prediction.

Legal AI, Attribute Extraction

††copyright: acmcopyright††journalyear: 2018††doi:
XXXXXXX.XXXXXXX††conference: Make sure to enter the correct conference title
from your rights confirmation emai; June 03–05, 2018; Woodstock, NY††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06

##  1\. Introduction

The exponential increase in the number of pending cases (such as over $10$
million111https://njdg.ecourts.gov.in/njdgnew/index.php pending cases in
India) highlights the pressing need for efficient processing of legal
documents using AI tools. Legal experts require time to thoroughly review
lengthy documents. On average, a case document contains $3,829$ words.
Existing research (Anand and Wagh, 2022) proposed approaches to provide
different kinds of summaries for a legal document.

However, a structured representation of a document is better (Munir and Sheraz
Anjum, 2018; Horrocks, 2013) in terms of efficiency in processing any
downstream task. To address this, we propose to represent a legal document in
terms of different important concepts or attributes. By collaborating with
legal experts, we formulated a comprehensive set of concepts specifically for
criminal case documents. Law students were appointed to annotate $200$ legal
documents (significant annotations compared to other legal works in India
Bhattacharya et al.; Malik et al.) with $7$ different attributes. However, it
is not possible to manually annotate a large number of documents.
Consequently, we explore existing sequence labeling approaches using deep
neural networks to effectively extract the concepts or attributes from legal
documents.

Broadly speaking, the contributions of this paper are as follows.

  1. (1)

We propose a new framework to present each legal document in terms of
different important concepts or attributes.

  2. (2)

We show the effectiveness of existing sequence labeling approaches in
automatically extracting attributes from legal documents using a limited
amount of training data.

  3. (3)

We also demonstrate the effectiveness of structured representation in a
downstream task (i.e. judgment prediction).

##  2\. Related Work

Existing literature related to our research scope can be broadly categorized
into three areas: a) Legal NLP b) Entity Recognition. Each one of them is
described as follows.

Legal AI In legal cases, the documents often encompass lengthy and intricate
sentences, making it challenging and time-consuming to thoroughly read and
comprehend the entire content of a case document. To alleviate this extensive
effort, researchers have emphasized the extraction of noun phrases known as
catchphrases (Tran et al., 2018; Mandal et al., 2017) from the document. This
approach aims to capture the key elements and central themes of the text.
Additionally, summarization (Anand and Wagh, 2022; Shukla et al., 2022)
techniques have been proven to be effective in gaining a comprehensive
understanding of the document by condensing its content into a concise
summary. However, to get the thematic view of the document, we need to extract
fine-grained level information or concepts from the document. We plan to
capture the thematic view of a legal document in our proposed attribute
framework for legal documents.

Figure 1. Annotation of a Legal Document

Entity Recognition Named Entity recognition (Tian et al., 2020) is a well-
established task in NLP. Entities are defined as objects which have an
independent existence. Examples of named entities are persons and
organizations. Domain-specific entities (e.g. biomedical entity) are in
general defined by domain experts through an ontology structure. Initially
rule-based systems (Sari et al., 2010), and syntactical structure (Zhang et
al., 2014) were used to automatically extract entities from a text. With the
advancement of machine learning and deep learning techniques existing research
has used conditional random fields, sequence-to-sequence models, and large
language models for Named entity recognition (Curran and Clark, 2003;
Hammerton, 2003; Sang and Meulder, 2003; Sun et al., 2018). More specifically
legal name entities extraction (Skylaki et al., 2020; Kalamkar et al., 2022)
has been implemented using Conditional Random Fields (CRF by used a sequence
tagger, specifically the BiLSTM-CRF model. To the best of our knowledge, this
is the initial endeavor to automate attribute extraction from legal
proceedings.

##  3\. Automatic Attribute Extraction

In this work, we particularly focused on criminal case documents from the
Supreme Court of India222Available from https://indiankanoon.org/. We propose
seven important legal concepts or attributes corresponding to a criminal case
document with the help of law experts. The different concepts and the
corresponding descriptions are as follows.

  1. (1)

ExpertWittest: Witness testimony from Expert people includes forensic and
ballistic experts.

  2. (2)

Wittest: Testimony from non-Expert witness testimony has been mentioned during
the judgment.

  3. (3)

Assault: Hurt by a sharp weapon.

  4. (4)

Riot: Unlawful enterprise in a violent manner.

  5. (5)

Homicide: Homicide amounting to murder.

  6. (6)

Imprisonment: Sentenced to life imprisonment.

  7. (7)

Evidence: Evidence of the crime was found.

In each document, legal experts highlighted the text span corresponding to a
concept present in that legal document. Figure 1 shows an example annotation.
Legal experts annotated $200$ legal documents with gold standard annotations.

###  3.1. Problem Definition

It is not possible to manually annotate all the legal documents. However, to
exploit the attributes corresponding to each case in a downstream task we
would need a lot of annotated documents. As a result of this, we focus on the
automatic annotation of legal documents with the attributes or concepts. We
cast the automatic annotation problem as a sequence labeling problem (Rei and
Søgaard, 2018). In this context, each token in the input text is assigned a
corresponding label. Sequence labeling has been extensively employed in NLP
research to tackle various challenges. The advantage of casting the automatic
annotation of a document into a sequence labeling problem is that along with
sequence labeling along with all the attributes mentioned in a case, we can
also obtain the highlighted sequence of tokens responsible for generating that
attribute.

As described in Section 3, there are $7$ important concepts or labels in a
legal document. All the tokens in a legal document may not be useful for
annotation. Consequently, we introduced an additional attribute named ‘NoTag’
to complement the aforementioned seven concepts (i.e., $8$ in total).

###  3.2. Methodology

Here we describe the working principle of a sequence labeling model. If $x$ is
the input sequence (i.e. $x=x_{1},x_{2},\ldots x_{n}$) and $y$ is the desired
output sequence (i.e. $y=y_{1},y_{2},\ldots y_{n^{\prime}}$) then
$n=n^{\prime}$. The sequence labeling problem learns a mapping
$H:X\xrightarrow{}Y$ using a state-of-the-art BiLSTM-CRF (Huang et al., 2015).
The BiLSTM-CRF network essentially computes the maximum likelihood of each
sequence label as follows.

(1) $P(y_{0,\ldots n}|r_{0\ldots n})=\prod_{i=1}^{n}f(y_{i-1},y_{i},x)$

In Equation 1, $f(y_{i-1},y_{i},x)$ estimates the likelihood of an individual
label (i.e. $y_{i}$) based on its neighboring sequence. The hidden state
output vectors ($r_{i}$) obtained from a BiLSTM model are used to estimate the
likelihood in $f(y_{i-1},y_{i},x)$.

##  4\. Experiment Framework

Dataset: For sequence labeling input only the highlighted spans from the
annotators are used for training. Two annotators manually annotated $200$
documents. While highlighting a span, annotators were only asked to choose the
sequence of words which are relevant. The highlighted text may not span over a
whole sentence or it may span across a sentence. The details of the dataset
are shown in Table 1.

For judgment prediction downstream task, we used a
dataset333https://github.com/Exploration-Lab/CJPE consisting of $5,082$
documents for training and $94$ documents for testing.

Dataset | Tags
---|---
Type | Statistics | ExpWittet | Wittest | Homicide | Assault | Imprisonment | Riot | Evidence
Train | #sentences | 31 | 35 | 50 | 39 | 50 | 32 | 24
Train | #tokens | 502 | 930 | 501 | 734 | 509 | 797 | 580
Test | #sentences | 372 | 162 | 78 | 48 | 34 | 20 | 10
Test | #tokens | 7789 | 3841 | 1036 | 755 | 344 | 410 | 257

Table 1. Dataset Statistics

| Accuracy
---|---
Method | ExpWitTest | WitTest | Homicide | Assault | Imprisonment | Riot | Evidence | Overall
SeqLabel (Roberta) | 0.85 | 0.50 | 0.18 | 0.39 | 0.40 | 0.78 | 0.24 | 0.67
SeqLabel (BERT) | 0.61 | 0.28 | 0.14 | 0.35 | 0.36 | 0.41 | 0.18 | 0.60
SeqLabel (LegalBERT) | 0.84 | 0.51 | 0.35 | 0.36 | 0.50 | 0.69 | 0.15 | 0.68

Table 2. Automated Tag Extraction Results using Flair

Preprocessing As a preprocessing step, each legal case document was split into
sentences using a sentence tokenizer from NLTK
444https://www.nltk.org/api/nltk.tokenize.html. Then each highlighted sentence
was split into tokens based on whitespace. While preparing the dataset for the
training of the sequence labeling model, the words not highlighted within a
sentence are tagged with NoTag. The sequence labeling approach is implemented
in Flair model555https://pypi.org/project/flair/ and is used for all our
experiments.

As suggested in (Malik et al., 2021) we utilized BERT-uncased to generate
embeddings for the last $510$ tokens of each document for judgment prediction
task. The embeddings were eventually used as input to a logistic regression
classifier. For the ‘with tag’ version in Table 3, we appended all the
extracted attributes of each test document to its corresponding content. Table
3 shows the method with tag overperforms the without tag version.

|  | Metric
---|---|---
Method | Acc (Class 0) | Acc (Class 1) |
Embedding | Input Format | Precision | Recall | Precision | Recall | Acc
BERT | Text | 0.50 | 0.34 | 0.54 | 0.70 | 0.53
BERT | Text+Tag | 0.54 | 0.52 | 0.59 | 0.62 | 0.56
BERT | Text+Span | 0.61 | 0.25 | 0.56 | 0.86 | 0.57
InLegalBERT | Text | 0.57 | 0.36 | 0.57 | 0.76 | 0.57
InLegalBERT | Text + Tag | 0.63 | 0.31 | 0.58 | 0.84 | 0.60
InLegalBERT | Text +Span | 0.75 | 0.40 | 0.62 | 0.88 | 0.66
Table 3. Judgement Prediction Results

##  5\. Results

We can observe from Table 2 that from all the different variations of the
Flair model, the Roberta version (i.e. SeqLabel (LegalBERT)) performed the
best in terms of overall accuracy. The accuracy for NoTag was more than $80\%$
for all the approaches. Consequently, we didn’t report it in Table 2.

To investigate the effectiveness of the legal document attributes in a
downstream task, we explored both ‘Text+Tag’ and ‘Text+span’ approaches in
Table 3. For the ‘Text+Tag’ version, we appended each test document’s
extracted attributes to its corresponding content. Similarly for the
‘Text+Span’ version, we appended the highlighted texts corresponding to each
tag to the content of the document. Table 3 shows that for both BERT and
InLegalBERT variations the ‘Text+Tag’ and ‘Text+Span’ versions performed
better compared to using only the text version. It is also interesting to
observe that the ‘Text+Span’ version is performing better compared to the
‘Text+Tag’ version. Consequently, we can conclude that the highlighted
portions corresponding to each tag are more meaningful in legal document
representation compared to the tag names.

##  6\. Conclusion

In this paper, we propose a novel attribute framework to capture the thematic
view of a legal document in a structured way. We also show the effectiveness
of a state-of-the-art sequence labeling approach in automatically predicting
the attributes of a legal document. We also demonstrate the usefulness of the
attributes in a downstream task. It is always difficult to manually annotate a
reasonable number of training data for a new task. In the future, we would
like to explore advanced deep-learning approaches to address the issue of not
having enough training data for automatic attribute extraction model training.

## References

  * (1)
  * Anand and Wagh (2022) Deepa Anand and Rupali Wagh. 2022.  Effective deep learning approaches for summarization of legal texts.  _Journal of King Saud University - Computer and Information Sciences_ 34, 5 (2022), 2141–2150.  https://doi.org/10.1016/j.jksuci.2019.11.015
  * Bhattacharya et al. (2019) Paheli Bhattacharya, Shounak Paul, Kripabandhu Ghosh, Saptarshi Ghosh, and Adam Wyner. 2019\.  Identification of Rhetorical Roles of Sentences in Indian Legal Judgments. In _Legal Knowledge and Information Systems - JURIX 2019: The Thirty-second Annual Conference, Madrid, Spain, December 11-13, 2019_ _(Frontiers in Artificial Intelligence and Applications, Vol. 322)_ , Michal Araszkiewicz and Víctor Rodríguez-Doncel (Eds.). IOS Press, 3–12.  https://doi.org/10.3233/FAIA190301
  * Curran and Clark (2003) James R. Curran and Stephen Clark. 2003.  Language Independent NER using a Maximum Entropy Tagger. In _Proceedings of CoNLL-2003_ , Walter Daelemans and Miles Osborne (Eds.). Edmonton, Canada, 164–167.
  * Hammerton (2003) James Hammerton. 2003\.  Named Entity Recognition with Long Short-Term Memory. In _Proceedings of CoNLL-2003_ , Walter Daelemans and Miles Osborne (Eds.). Edmonton, Canada, 172–175.
  * Horrocks (2013) Ian Horrocks. 2013\.  _What Are Ontologies Good For?_ 175–188.  https://doi.org/10.1007/978-3-642-34997-3_9
  * Huang et al. (2015) Zhiheng Huang, Wei Xu, and Kai Yu. 2015.  Bidirectional LSTM-CRF Models for Sequence Tagging.  arXiv:1508.01991 [cs.CL]
  * Kalamkar et al. (2022) Prathamesh Kalamkar, Astha Agarwal, Aman Tiwari, Smita Gupta, Saurabh Karn, and Vivek Raghavan. 2022\.  Named entity recognition in indian court judgments.  _arXiv preprint arXiv:2211.03442_ (2022).
  * Malik et al. (2021) Vijit Malik, Rishabh Sanjay, Shubham Kumar Nigam, Kripa Ghosh, Shouvik Kumar Guha, Arnab Bhattacharya, and Ashutosh Modi. 2021\.  ILDC for CJPE: Indian legal documents corpus for court judgment prediction and explanation.  _arXiv preprint arXiv:2105.13562_ (2021).
  * Mandal et al. (2017) Arpan Mandal, Kripabandhu Ghosh, Arindam Pal, and Saptarshi Ghosh. 2017.  Automatic Catchphrase Identification from Legal Court Case Documents. In _Proceedings of the 2017 ACM on Conference on Information and Knowledge Management_ (Singapore, Singapore) _(CIKM ’17)_. Association for Computing Machinery, New York, NY, USA, 2187–2190.
  * Munir and Sheraz Anjum (2018) Kamran Munir and M. Sheraz Anjum. 2018.  The use of ontologies for effective knowledge modelling and information retrieval.  _Applied Computing and Informatics_ 14, 2 (2018), 116–126.
  * Rei and Søgaard (2018) Marek Rei and Anders Søgaard. 2018.  Zero-Shot Sequence Labeling: Transferring Knowledge from Sentences to Tokens. In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_. Association for Computational Linguistics, New Orleans, Louisiana, 293–302.  https://doi.org/10.18653/v1/N18-1027
  * Sang and Meulder (2003) Erik F. Tjong Kim Sang and Fien De Meulder. 2003.  Introduction to the CoNLL-2003 Shared Task: Language-Independent Named Entity Recognition.  _CoRR_ cs.CL/0306050 (2003).  http://arxiv.org/abs/cs/0306050
  * Sari et al. (2010) Yunita Sari, Mohd Fadzil Hassan, and Norshuhani Zamin. 2010\.  Rule-based pattern extractor and named entity recognition: A hybrid approach. In _2010 International Symposium on Information Technology_ , Vol. 2. 563–568.
  * Shukla et al. (2022) Abhay Shukla, Paheli Bhattacharya, Soham Poddar, Rajdeep Mukherjee, Kripabandhu Ghosh, Pawan Goyal, and Saptarshi Ghosh. 2022\.  Legal Case Document Summarization: Extractive and Abstractive Methods and their Evaluation. In _Proceedings of the 2nd Conference of the Asia-Pacific Chapter of the Association for Computational Linguistics and the 12th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_. Online only.
  * Skylaki et al. (2020) Stavroula Skylaki, Ali Oskooei, Omar Bari, Nadja Herger, and Zac Kriegman. 2020.  Named entity recognition in the legal domain using a pointer generator network.  _arXiv preprint arXiv:2012.09936_ (2020).
  * Sun et al. (2018) Peng Sun, Xuezhen Yang, Xiaobing Zhao, and Zhijuan Wang. 2018\.  An Overview of Named Entity Recognition. In _2018 International Conference on Asian Language Processing (IALP)_. 273–278.  https://doi.org/10.1109/IALP.2018.8629225
  * Tian et al. (2020) Yuanhe Tian, Wang Shen, Yan Song, Fei Xia, Min He, and Kenli Li. 2020.  Improving Biomedical Named Entity Recognition with Syntactic Information.  https://doi.org/10.21203/rs.3.rs-21994/v1
  * Tran et al. (2018) Vu D. Tran, Minh Le Nguyen, and Ken Satoh. 2018.  Automatic Catchphrase Extraction from Legal Case Documents via Scoring using Deep Neural Networks.  _CoRR_ abs/1809.05219 (2018).
  * Zhang et al. (2014) Xiantao Zhang, Dongchen Li, and Xihong Wu. 2014.  Parsing named entity as syntactic structure. In _Interspeech_.  https://api.semanticscholar.org/CorpusID:33212550


|  Low-energy pathways lead to self-healing defects in CsPbBr3†
---|---
|  Kumar Miskin,∗a Yi Cao,b Madaline Marland,b Jay I. Rwaka,b Farhan Shaikh,b
David T. Moore,c John A. Marohnd and Paulette Clancy b
|  Self-regulation of free charge carriers in perovskites via Schottky defect
formation has been posited as the origin of the well-known defect tolerance of
metal halide perovskite materials that are promising candidates for
photovoltaic applications, like solar cells. Understanding the mechanisms of
self-regulation, here for a representative of more commercially viable all-
inorganic perovskites, promises to lead to the fabrication of better-
performing solar cell materials with higher efficiencies. We investigated
different mechanisms and pathways of the diffusion and recombination of
interstitials and vacancies (Schottky pairs) in CsPbBr3. We use Nudged Elastic
Band calculations and ab initio-derived pseudopotentials within Quantum
ESPRESSO to determine energies of formation, migration, and activation for
these defects. Our calculations uncover defect pathways capable of producing
an activation energy at or below the value of 0.53 eV observed for the slow,
temperature-dependent recovery of light-induced conductivity in CsPbBr3. Our
work reveals the existence of a low-energy diffusion pathway involving a
concerted “domino effect” interstitial mechanism, with the net result that
interstitials can diffuse more readily over long distances than expected.
Importantly, this observation suggests that defect self-healing can be
promoted if the “domino effect” strategy can be engaged.

##

††footnotetext: a Department of Materials Science and Engineering, Johns
Hopkins University, USA. E-mail<EMAIL_ADDRESS>b Department
of Chemical and Biomolecular Engineering, Johns Hopkins University, USA
††footnotetext: c National Renewable Energy Lab, USA††footnotetext: d
Department of Chemistry and Chemical Biology, Cornell University,
USA††footnotetext: † Electronic Supplementary Information (ESI) available: See
DOI: 00.0000/00000000.

##  1 Introduction

Metal halide perovskites (MHPs) are one of the most exciting materials classes
at the forefront of photovoltaic materials development. They have demonstrated
performance comparable to market-leading solar cell materials, as evidenced by
their dramatic rise in power conversion efficiency (PCE) since emerging about
a decade ago.19 MHPs, like all perovskites, adopt the general formula ABX3,
where A is a monovalent cation (MA^+, FA^+, Cs^+), B is a divalent metal
cation (Sn^2+, Pb^2+), and X is a halide (Br^-, I^-, Cl^-). Early work on MHPs
centered on methylammonium (MA) lead iodide (or MAPbI3). More recently, the
field is migrating towards all-inorganic perovskites, like CsPbX3, with band
gaps in the range 1.7-2.3 depending on the halide. These perovskites have
proven to be somewhat more durable in conditions (UV, heat, and humidity) that
can degrade MAPbI3 and other hybrid organic-inorganic perovskites over time.
Due to its importance, we focus our defect studies on CsPbBr3.

Defect migration and recombination strongly affect device performance and
long-term durability and are generally unwanted in photovoltaic devices. 40 In
that regard, metal halide perovskites (MHP) are known to be far more defect-
tolerant than their silicon solar cell counterparts.39, 34, 14 Surprisingly,
the underlying atomic-level mechanisms and prevalent pathways for defect
creation and remediation remain largely unknown. 4 The most widely accepted
theory posits that the most probable perovskite defects form “shallow” traps
near the conduction or valence band and therefore still allow long lifetimes
for charge carriers. 29, 14 Conversely, experimental work argues that deeper
traps must be present 12 and that ionic defects, which may be created by light
16, 28 and voltage,28 are a source of device hysteresis and instability.15, 24

Experimental measurements of defect activation energies are rare, especially
for the all-inorganic perovskite studied in this paper. Here, we obtain
activation energies, Ea, computationally: The activation energy is determined
from the sum of the energy to form the defect, Ef, plus the energy to move the
defect in the lattice, Em. Both Ef and Em can be determined computationally,
here using density functional theory calculations, and hence activation
energies can be estimated.

In terms of extant defect studies, especially computational studies, a larger
volume of prior work exists for hybrid organic-inorganic systems, with
methylammonium (MA+) and formamidinium (FA+) A-site cations in the lattice,
compared to fewer studies for inorganic MHP systems. 1, 23, 36 Meggiolaro and
De Angelis provided a fairly comprehensive review of ab initio-generated
formation energies for MAPbI3 and offered “best practices” for such
calculations. 21 Oranskaia et al. found lower hydrogen-bonding strength in
MAPbBr3 than FAPbBr3 and, consequentially, lower halide migration and defect
formation energies (DFE). 26 This finding suggests that careful A-site cation
selection may reduce perovskite ion migration.

While the available literature reports defect energies, mechanistic
understanding is currently limited, especially for all-inorganic MHP
materials. Evarestov et al. found that perfect CsPbX3 structures follow the
stability trend CsPbBr3 > CsPbI3 > CsPbCl3, corroborating CsPbBr3 defect
formation energies from Kang and Wang. 5, 14 In CsPbI3 systems, Xue et al.
found that an anti-site defect of Pb on a Cs site forms with comparable ease
to other defects. 37 Balestra reported significantly lower activation energy
for vacancy migration for mixed halide CsPb(Br_xI_1-x)3 ($x=1/3$, $2/3$) than
for homogeneous halide systems in which all three X atoms are the same halide.
2 Ten Brinck et al. calculated defect formation energy values for
interstitial, vacancy, and anti-site defects in non-periodic 3 $nm$ CsPbBr3
nanoparticles and concluded that the most stable defect position is on the
nanoparticle surface. 30

Our calculations are motivated by Tirmzi _et al._ ’s measurements of the slow
recovery of light-induced conductivity in CsPbBr3.32 In these measurements, a
charged atomic-force microscope cantilever was brought near the surface of a
CsPbBr3 film and experienced a non-contact friction depending on the sample’s
total conductivity, electronic plus ionic.31 After irradiating the sample with
above-bandgap light, Tirmzi _et al._ observed a prompt increase in friction,
as expected from the presence of photo-generated electrons and holes.
Electrons and holes recombine on the microsecond timescale in perovskites.
After the light was turned off, however, sample friction took $10$ s to
recover at room temperature. Tirmzi _et al._ concluded that the sample’s
conductivity was dominated by ionic conductivity that that light was
_creating_ mobile ionic defects:

$\mathrm{nil}\ce{->[$h\nu$]}\mathrm{Br}_{\mathrm{Br}}^{\mathrm{x}}+\mathrm{h}^{\bullet}+e^{\prime}\ce{->}\mathrm{V}_{\mathrm{Br}}^{\bullet}+\mathrm{Br}_{\mathrm{i}}^{\mathrm{x}}+e^{\prime},$
(1)

in Kröger-Vink notation. The friction recovery time in the dark was
temperature dependent, with an activation energy of $E_{\mathrm{a}}=0.53\pm
0.03\>\mathrm{eV}$ (95% confidence interval). The associated defect-recovery
reaction in the dark is

$\mathrm{V}_{\mathrm{Br}}^{\bullet}+\mathrm{Br}_{\mathrm{i}}^{\mathrm{x}}+e^{\prime}\ce{->}\mathrm{nil}.$
(2)

For this reaction we would expect the measured $E_{\mathrm{a}}$ to be
dominated by the _migration_ energy governing the diffusion of the neutral
interstitial $\mathrm{Br}_{\mathrm{i}}^{\mathrm{x}}$.

In this paper, we have modeled a number of potential recombination pathways
for interstitial, vacancy, and anti-site defects using Nudged Elastic Band –
Density Functional Theory (NEB-DFT). Unlike previous studies in the
literature, we have calculated both formation and migration energies for a
range of defect pathways, allowing us the ability to calculate activation
energies, an experimentally accessible property. Coupling Nudged Elastic Band
calculations with ab initio pseudopotentials provides accurate estimates for
the migration energies. The conductivity-recovery rate observed by Tirmzi and
coworkers 32 was not noticeably different near grain boundaries; they were
observing a bulk effect. This observation justifies our ignoring ground-
boundary effects in the following computations.

We are looking to find a match between computational and experimental
activation and migration energies and, as a result, uncover the dominant
defect(s) and defect pathways that are present. With this information in hand,
we hope to better inform device synthesis and offer a different bulk-defect
passivation strategy to improve overall MHP performance.

##  2 Results and Discussion

Fig. 1: Flowchart of defect studies in CsPbBr3. This begins with an initial
lattice set-up (top left; blue box). Then a defect is created, whether a
vacancy or an interstitial (top right; pink box). Converged calculations yield
a formation energy. Additional calculations using NEB produce the migration
energy (bottom right; green box), for a set of considered pathways (middle;
orange boxes). Converged NEB calculations produce Em. Then the activation
energy can be calculated (dark green cube).

To validate our NEB calculations, we first predicted an X-ray diffraction
(XRD) spectrum generated from our simulated pristine supercell and VESTA
software25 and compared it against that found experimentally for a thin film
MHP by our collaborators. As seen in the comparison, found in the
Supplementary Information, the simulated peaks line up essentially perfectly,
with the density functional theory-generated observations matching the
experimental results. The relaxed cell parameters for the pristine CsPbBr_3
supercell ($3\times 3\times 2$) (with Pnma spacegroup) were found to be 8.267,
8.177 and 11.768 Å. These correspond very favorably to experimentally reported
values:18 8.23, 8.18 and 11.73 Å.

###  2.1 Defect Formation Energies

First, we examined defect formation energy values in CsPbBr3, focusing on the
consideration of various kinds of interstitial sites, as mentioned in the
Methods section. The Br_i atom was placed in a chosen interstitial position
and allowed to relax in the crystal. We explored different Br interstitials
based on their local electronic environment as explained in more detail later
in Section 2.2.

For the Br_i-Br_i split-interstitial system, a Br interstitial shares a
lattice site with another Br atom. For the Br_i-Pb_i split-interstitial
system, a Br interstitial shares a lattice site with a Pb atom. The defect
formation energy (DFE) for this configuration varies depending on the axis
along which the Pb-Br pair is aligned. As depicted in Figure S2 of the
Supplementary Information, the DFE is found to be $0.83$ eV, on average, when
the pair is aligned along the $z$-axis. In contrast, the DFE value increases
to $1.14$ eV when the split interstitial is aligned along the $y$-axis. This
difference of approximately 0.4 eV between the two orientations is significant
and consistent, indicating a pronounced anisotropic phenomenon in defect
formation.

In the Pb_Br and Br_Pb anti-site systems, the defect formation energy varied
with the location of the introduced anti-site. It is observed that, as the
anti-site is placed closer to the center of the supercell system, the DFE
comes closer to ten Brinck’s experimental values. 30 This can be explained by
the effect of periodic boundary conditions on the simulation, wherein it is
observed that defects near the edges of the supercell tend to result in
distorted energy values. All the values are reported in the Supplementary
Information. Our calculations resulted in large formation energies of 6.08 eV
and 1.30 eV for the Pb_Br and Br_Pb anti-site defects, respectively. These
values are similar to those reported by Kang et al. 14 Pb_Br shows a high DFE
due to the introduction of a larger Pb atom in place of a smaller Br, causing
a considerable shift in neighboring atomic positions.

Table 1 reports DFE values for all the defects that were studied in this work.

###  2.2 Defect Migration Energies

Our first task here was to identify possible interstitial sites for an extra
bromine atom relaxed at different sites in the crystal. We then looked at the
neighborhood of these different relaxed bromine interstitials to map out
possible sites. Figure 2 lists the neighboring interatomic distances from the
Bri atom for each CsPbBr3 relaxed structure without a paired vacancy. These
correspond to different sites to which the interstitial Br relaxed in the
supercell. The letters A,B,C.. etc. represent different neighborhoods for the
relaxed Br interstitial, further visualized in Figure S.2.1.

Fig. 2: Locations of first, second, third, and fourth neighbor interatomic
distances (Å) from a reference neutral interstitial Bri in (a) relaxed
$2\times 1\times 1$ CsPbBr3 structures and (b) neighbor distances found at the
NEB transition state, sorted from lowest Em to highest on the $x$-axis. The
color key for both graphs is shown inset on the right-hand side. The letters
A,B,C.. etc. represent different neighborhoods for the relaxed Br
interstitial. This figure demonstrates that Em is dependent not only on the
initial and final structure local electronic environments but also on the
transition structure local electronic environment.

A change in the coordinates of the pre-relaxation Bri of as little as 10% in
one direction resulted in different identities of its neighbors and
corresponding distances between the reference Bri atom and the neighbors, as
shown between systems C and D. We observed that, while Bri systems C and D
exhibited different sets of neighbors after structural relaxation, the
resulting migration energies from a reference system to systems C and D,
respectively, were remarkably similar (Em = 0.78 eV for system A to system C,
versus Em = 0.80 eV for system A to D) in Figure 3. Those differences are
within the uncertainty of the predictions themselves. These small differences
indicate that considering only the local electronic environment of the initial
or final system does not fully explain the origin of the $E_{m}$ value. In
contrast, systems C and D show similar transition state structures, i.e., at
the saddle point, along the NEB pathway from A to C and A to D, respectively;
Figure 3b. Thus, the local Bri environment in both the initial structure and
the transition structure influences the value of Em.

Bringing all these works together, we calculate the more experimentally
accessible activation energies for all the defects we studied. Table 3 reports
the activation energy ($E_{a}$) associated with Bri returning to a vacant
lattice site through direct recombination (creating a perfect crystal).

###  2.3 Potential Mechanisms for Defect Migration

We studied NEB calculations of the movement of interstitial sites occurring
via two distinct mechanisms: The first involved the assistance of vacancies in
the lattice, while the second relied on a coordinated movement of the
interstitials themselves.

####  2.3.1 Vacancy-Assisted Mechanism

To delve deeper into the intricacies of migration dynamics, we strategically
introduced a vacancy into the lattice structure along with an interstitial.
This approach allowed us to quantify the activation energies necessary for an
interstitial atom to migrate towards, and subsequently recombine with, the
vacancy (resulting in a healed perfect crystal). Table 3 reports the
activation energy ($E_{a}$) associated with Bri returning to a vacant lattice
site through direct recombination.

In the context of vacancy-assisted migration, our investigations unveiled a
“domino effect” strategy that can assist the translocation of an interstitial
across the lattice, a phenomenon particularly pronounced in interactions
involving second and third nearest neighbors. Rather than migrating directly
into the vacancy site, the interstitial atom induces a successive movement of
one or two adjacent atoms along the pathway, effectively achieving the final
state through a series of intermediate steps. This phenomenon is clearly
demonstrated in Figure 4, which illustrates how the NEB-generated activation
energy can be ’flattened’ by this domino effect, providing an easier migration
pathway for the atoms. In Figure 4, we can clearly see that the barrier for
long-distance jumps of the interstitial was significantly reduced by the
domino effect. We are unaware of any other references to such a novel and low-
energy pathway in the literature, marking this out as an important finding.

We calculated the average mean square displacement (MSD) of each atom from its
lattice site. This was computed for each pathway at the saddle point of
diffusion moves that involve recombination of defects. Looking at the MSD for
(a) and (c) in Figure 4, we can see that it follows a linear regime in the
direct jump scenario: The activation energy increases with increasing MSD as
might be expected. In contrast, in the domino effect strategy, (b) and (d) in
Figure 4, the activation energy tapered off for sufficiently large MSD values
of the interstitials. This shows that, for any migration more than a second or
a third neighbor distance away, a domino movement is far more probable than a
direct jump.

A summary of the NEB findings can be found in Figure 8, offering insights into
the energetics of the vacancy-assisted migration process.

We examined, in some detail, one of the recombination pathways that involved a
third neighbor jump. Here, we start from the direct jump and sequentially add
one atom at a time to see how the number of intermediate atoms affects the
final migration energy. As depicted in Figure 4, the direct jump has two local
energy peaks in the migration process, indicating two relatively high energy
barriers to overcome. When an intermediate atom is introduced, only one peak
remains in the migration pathway and the migration energy is reduced. We posit
that the involvement of an intermediate atom assisted in helping the original
interstitial atom bypass one of the barriers, although one barrier still
remains. Introduction of yet another intermediate atom flattens the left of
the two barriers and the resulting migration energy is effectively reduced to
zero. This suggests that, even if the defect is relatively distant from the
vacancy site, it can still migrate and recombine simultaneously. This suggests
that some defects can self-heal without additional energy input if the “domino
effect” strategy can be engaged.

Fig. 3: Neutral interstitial Bri Em (eV) in CsPbBr3 structures.

Table 3 reports the activation energy ($E_{a}$) associated with Bri returning
to a vacant lattice site through direct recombination to a perfect crystal.

Fig. 4: Migration energy of a Br_i-Pb_i split interstitial following two
distinctly different pathways and its effect on the first, second and third
nearest neighbors of the moving Br atom (first three columns, respectively).
(a, top) The energy barrier for the migration of a Br_i-Pb_i split
interstitial via a direct jump to a neighboring vacancy site. Energy profiles
for the first, second, and third nearest neighbors are grouped based on
distance to the Br interstitial site. Each line stands for a different
migration pathway which is color-coded by their activation rank order. (b,
bottom) The energy barrier associated with deploying a “domino effect”
strategy for inter-lattice diffusion. (c) Maps the correlation between
activation energy and mean squared displacement (MSD) for the direct jump
mechanism. Color-coded dots, based on activation energy, are rank-ordered with
shaded areas indicating confidence intervals for the regression curve. This
shows that Ea increases with MSD, as expected. (d) Displays a similar
correlation to (c), except that this relates to the “domino effect” diffusion
mechanism. Here, unlike in (c), the activation energy is pretty flat
regardless of the MSD value. Fig. 5: Summary of the activation energy values
in this study for defect recombination events plotted as a function of the
Mean Squared Displacement for a perfect crystal involved in the pathway and
type of defect studied. Color key for the data shown is given on the top RHS
of the figure. The plot shows that forming a variety of Pb-Br split
interstitials in the presence of a Br vacancy (blue dots) offers many
opportunities for low activation energy events (< 0.4 eV). The single
experimental activation energy for CsPbBr3 from Tirmzi et al. is shown as a
red dashed line with the error bar for measurement shown as shaded blue. Fig.
6: Migration energy differences for Br_i migration with one (red line) or two
(orange line) intermediate Br atoms participating in the motion compared to
zero atoms in the defect transition (purple line). The differences in
migration energies are lower as one, and then two, intervening atoms are
involved, representative of the “domino effect” strategy.

####  2.3.2 Interstitial Migration Mechanism

In addition to the previously described vacancy-assisted migration, our
investigations have revealed that Br interstitials are quite capable of
migrating to alternate interstitial sites. This migration energy is influenced
by the local environment of the interstitials in all three structures
(initial, transition and final) as we saw in Section 2.2. Notably, when the
vacancy is not the closest neighbor to the interstitial, the “domino effect”
strategy appears to facilitate this migration as well.

The NEB findings, in the form of migration energy barriers and a visualization
of the pathway mechanism, as captured in the initial, transition, and final
states of the system, are summarized in Figure 7. This Figure offers a compact
overview and associated insight into the energetics and mechanisms of the
defect migration process. These results clearly indicate that the most
favorable migration pathway amongst the ones we studied features a collective
and concerted movement of several atoms, consistent with the interstitial
“domino effect” in vacancy-assisted migration.

Figure6 accurately captures this lowering of the activation energy through the
“domino effect.” We see that, for an interstitial to migrate, having one more
intermediate atom involved reduced the activation barrier for migration.
Adding a second intermediate atom reduced this barrier almost to zero. This
shows that for any defect travel larger than the nearest neighbor, the
migration barrier can be reduced by a domino effect of atoms moving in a
concerted exchange for metal halide perovskites.

Table 1: Defect Formation Energies in CsPbBr3 Defect444Br-rich Environment | Formation Energy (eV) | Other Work (eV)
---|---|---
V_Br | 2.71 | 2.6714
Br_i | 0.75 | 0.7014
Pb_Br | 6.08 | 6.2814
Br_Pb | 1.30 | 1.4014
Br_i-Br_i split | 1.55 | N/A
Br_i-Pb_i split | 1.19 | N/A
Cs_Pb | 0.67 | 0.8114
Pb_Cs | 2.87 | 2.44 14
Table 2: Defect Migration Energies in CsPbBr3 Defect Pathway | Migration Energy (eV)
---|---
V_Br | 0.89
V^._Br | 0.75
Bri | 0.53-0.81
Br_i-Br_i split | 0.62
Br_i-Pb_i split | 0.48
Table 3: Activation Energies for Defect Recombination in CsPbBr3 Recombination Pathway | Activation Energy (eV)
---|---
Bri | 1.14
Pb_Br= Br_Pb | 0.17
Br_i-Br_i split | 0.70
Br_i-Pb_i split | 0.10$-$0.14
Fig. 7: Summary of Defect Migration Energies for four defects: From top to
bottom, a neutral Br vacancy, a positively charged Br vacancy, a Br
interstitial (no vacancy) and a split interstitial involving lead and bromine
atoms. Left column: The NEB-generated energy barrier. The right three columns
show a visualization of the initial state, the transition state, and the final
state, as indicated. Atoms not engaged in the pathway are shown with reduced
intensity. Those actively engaged in the pathway are shown in more intense
colors. Key: Red = Br; Black = Pb. In all images, the red tail, shown as a
solid line, indicates the path taken by the Br atom during the pathway Fig. 8:
Summary of Recombination Energies for five defects: From top to bottom, a Pb-
Br anti-site, a Br-Br split interstitial, and three Pb-Br split-interstitials
that are aligned in the $x-$, $y-$ and $z-$ directions, respectively in the
lattice. Left column: The NEB-generated energy barrier. The right three
columns show a visualization of the initial state, the transition state, and
the final state, as indicated. Atoms not engaged in the pathway are shown with
reduced intensity. Those actively engaged in the pathway are shown in more
intense coloring. Key: Red = Br; Black = Pb. In all images, the red tail,
shown as a solid line, indicates the path taken by the Br atom during the
pathway

A summary of the NEB-derived formation energies, migration energies and
activation energies investigated in this work can be found in Tables 1,2 & 3,
respectively.444Br-rich environment

##  3 Methodology

Nudged Elastic Band (NEB) calculations are a well-established methodology for
finding the minimum energy pathway that follows the transition of one
arrangement of atoms (the initial state) to a different arrangement (the final
state). 22, 3, 11, 10 These calculations give rise to an estimation of the
maximum potential energy needed to effect this transition, the “saddle point
energy,” which is assumed to be synonymous with the activation energy barrier
for the transition.

In the calculations here, we estimate both the defect’s formation energy and
migration energy, the energy required to move the defect from point A to point
B in the supercell. Both of these energy barriers contribute to the overall
activation energy for CsPbBr3, as explained in the Introduction. We use a
$3\times 3\times 2$ supercell of CsPbBr3 containing a total of 360 atoms for
the calculation of defect formation energies. This is close to the practical
limit of system size for a reasonably accessible computational resource
effort. Perdew-Burke-Ernzerhof (PBE) calculations were performed using
ultrasoft pseudopotentials33 with a plane wave cut-off of 40 Ryd (400 Ryd on
the charge density). These pseudopotentials were chosen because they
accurately represent defects in these kinds of systems.21 The well-known
electronic structure code, Quantum ESPRESSO, 7, 6, 8 was used for all density
functional theory (DFT) calculations reported below. The migration energies
and pathways for recombination of these defects were estimated using NEB. Due
to the high computational cost of running Nudged Elastic Band simulations, we
chose a smaller $2\times 1\times 1$ supercell of 40 atoms for the migration
energy and the activation barrier for recombination. Details of all the
parameters used to set up the Quantum ESPRESSO calculations, as well as sample
scripts, are provided in the Supplementary Information.

Based on Kang et al.’s work and that of others, we selected point defects that
are thought to be particularly prevalent in lead halide perovskites. Some were
chosen for their low defect formation energies.14 While much of our focus
involved looking at simple point defects (V, I), we also explored the
migration and recombination of additional ones not considered previously,
including “split interstitial” defects in which two atoms share a lattice
site, as well as more concerted pathways.

In the ab initio calculations of defect energetics, we introduced a specific
defect into a perfect lattice, fixed the cell parameters, and then allowed the
atoms to relax. The workflow for such calculations is shown in Figure 1. The
defect formation energy (DFE) was then computed based on the total energies of
the relaxed structures, using the following equation: 35

$E^{f}(D^{q})=E(D^{q})-E_{\text{bulk}}-\sum_{i}\Delta
n_{i}\mu_{i}+q(E_{V}+E_{F})+E_{\text{corr}},$ (3)

where $E^{f}(D^{q})$ represents the defect formation energy of a defect $D$ in
charge state $q$ and $E_{\text{bulk}}$ the total energy of the pristine bulk
supercell. This equation 3 incorporates the changes in the number of atoms
($\Delta n_{i}$) & the chemical potentials ($\mu_{i}$). The Fermi level
($E_{F}$) of a semiconductor is treated as an independent variable that can
assume any value within the band gap ($E_{V}$), the valence band maximum (VBM)
of the bulk material. The correction term ($E_{\text{corr}}$) for finite size
effects and periodic boundary conditions.

For the NEB calculations of the migration energy, we began by choosing initial
and final configurations of a defective lattice and allowed them to relax. We
created ten images between the initial and final states for each NEB-derived
energy barrier to investigate a manually selected recombination or migration
defect pathway.

As well as studies of simple vacancies and interstitials and anti-site
defects, we also considered some more exotic, concerted defect pathways. A
particularly interesting example involves a defect that “moves” in a concerted
manner involving multiple atoms. It moves such that the interstitial atom you
start with is not the same one that eventually moves across the lattice — a
“domino effect.” As we shall show, having other bromine interstitial atoms
participate in the migration mechanism drastically reduces the associated
activation barrier. A more comprehensive explanation of the “domino effect” in
CsPbBr3 can be found in the Supplementary Information. We identified the
intermediate atoms along the “domino effect” pathway in the NEB calculation
and compared the migration energy with that of a direct jump, as well as a
second neighbor jump (involving two Br atoms) and a third neighbor jump
(involving three Br atoms). These comparisons allowed us to examine the
contribution of each additional atom to the change in migration energy.

We now consider the creation of a set of defective systems studied in this
paper.

###  3.1 Br vacancies

Halide point defects, especially vacancies, are expected to be one of the
most, if not the most, readily moveable defect in MHPs.1, 13 As there is a
single unique site for bromine atoms, a vacancy can be formed by removing a
bromine atom from the simulation cell and letting the system relax. We also
looked at the effect of migration of a vacancy in both a charged and a neutral
supercell. The migration of a vacancy to its nearest neighbor site and the
resulting activation energy associated with the creation and migration of
these vacancies are given in Table 1.

###  3.2 Br interstitials

To create an interstitial Br in the system (Bri), we placed a Br atom at a
location offset from a bulk crystal reference site by a user-specified amount,
typically 10-25% in the $x$-, $y$-, $z$\- directions, to test the tolerance of
the resulting relaxed configuration to its starting point. Then we performed a
relax calculation on the system. This is a different approach to other ways to
predict the Bri locations, e.g., from a Voronoi analysis or other methods. 9
We characterized six unique interstitial locations, determined by their
differing nearest neighbor distances and atom types, Figure 2. Their reference
coordinates can be found in the Supplementary Information, Table S2.1. To
determine Em, different Bri systems were selected to be the first and final
NEB images, with the resultant “forward” direction selected as the lower Em
value. For recombination, the activation energy ($E_{a}$) is the energy
associated with Bri returning to a vacant lattice site through direct
recombination to create a perfect crystal.

###  3.3 Split interstitials

A split interstitial consists of two self-interstitials sharing a single
substitutional lattice site. These are some of the most common interstitials
found in semiconductor materials, 20 especially amongst the halide defects27.
We explored the defect migration energies of split interstitials involving two
bromine atoms sharing one site (Br-Br) and a bromine and lead atom sharing a
site (Br-Pb). Both of these are located at the bromine lattice sites. Table 2
summarizes the formation energy for creating these split interstitials.

###  3.4 Anti-sites

An anti-site refers to the presence of a “foreign” (unexpected) atom on a
lattice site. Here, we explore two types of anti-sites: The first involves
lead on a bromine site, and bromine on a lead site (Pb_Br and Br_Pb), followed
by simulation of the recombination of a pair of anti-sites Br_Pb= Pb_Br into
the null (perfect) crystal system. The second involves Pb on a Cs site, and
vice versa.

The atom sites were carefully selected to ensure that they are located at the
center of the supercell. To simulate the recombination of a pair of anti-sites
into the bulk crystal, the two displaced atoms had to be introduced within a
first nearest neighbour distance. Twenty such candidate sites were identified.
Amongst these, the lowest energy system was selected to further investigate
the recombination of a pair of anti-sites into the perfect crystal. Defect
formation energies were calculated for single anti-site defect and tabulated
in Table 1.

##  4 Conclusions

We have studied a broad range of defects in the CsPbBr_3 metal halide
perovskite, including the defect formation energies (Ef) and migration
energies (Em) studied using accurate ab initio Nudged Elastic Band
calculations. The latter will be a helpful addition to the literature. The
formation energies reproduce those of Kang and Wang 14. Unlike other studies,
the combination of both Ef and Em allows us to calculate the activation energy
of the various defects, which could be compared to experiment. Above, we
highlighted Tirmzi _et al._ ’s experiment measuring the slow, activated
recovery of light-induced conductivity in CsPbBr3. We find that the calculated
$E_{m}$ for $\mathrm{V}_{\mathrm{Br}}^{\bullet}$ agrees well with the 0.53 eV
activation energy for relaxation of light-induced conductivity observed by
Tirmzi and co-workers. Beyond the more well-studied intrinsic defects, namely,
a Br vacancy and a Br interstitial defect, we have included additional defects
like the family of split interstitials (Br-Br and Pb-Br) as well as anti-site
defects. Finally, we have identified a formerly unexplored diffusion pathway
which we have labeled the “domino effect” strategy involving multiple proximal
Br interstitials. The merits of this latter pathway are described below.

The low activation energies that we observed for some of the recombination
mechanisms can help explain the self-regulation process seen in perovskites.
38 Defects like vacancies, interstitials and anti-sites help the perovskite
“heal” by adjusting the concentration of charge carriers. The pathways that
the 10-window NEB calculations allowed us to explore give us accurate energy
barriers as well as insight into the different mechanisms of self-regulation
that metal halide perovskites can undergo. The domino effect of recombination
aided by other nearby bromine atoms lowered the activation threshold, a
phenomenon we saw across most of the pathways we studied.

Looking at a holistic picture of all the recombination pathways we studied and
their associated energy barriers, the lowest value of the activation energy
that we observed (see Figure 5) was 0.10 eV and a number of Pb-Br split
interstitials exhibited energy barriers below 0.4 eV. Br-Br split
interstitials in the presence of a vacancy can also heal with a barrier under
0.7 eV. An experimentally measured activation energy will be an ensemble
average of a number of pathways, and must surely include many of the ones we
have uncovered here.

Investigating the pathway that results in the lowest barrier (0.1 eV), this
mechanism is governed by the process V_Br + Br_i-Pb_i split -> Br_Br + Pb_Pb,
with the interstitial being located at the third nearest neighbor to the
vacancy and the reaction being aided by two other proximal bromine atoms in
the lattice. Because of its low formation energy and the ease of movement for
split interstitial defects in the system, the activation energy is lowered by
quite a large margin. Adding to this lower energy, the participation of other
bromine atoms to create a “domino effect” drastically decreases the energy
barrier by reducing the total displacement that each bromine has to travel to
return to the defect-free bulk state. As the defect moves, it affects other
intermediary atoms, transferring energy and reducing the overall activation
energy needed. Essentially, the intermediary atom acts as a stepping stone,
making defect movement more efficient. This “domino effect” strategy,
involving a collective movement of defects in inorganic metal halide
perovskites, like CsPbBr3, enables a drastic reduction of the activation
energy needed to traverse the lattice.

The calculations of defect pathways as well as Ef, Em, and Ea values provided
in this paper can have an impact on the design of better materials and,
ultimately, devices. Here are two pertinent suggestions: First, the links we
provide between structure and their corresponding formation energies can
assist the the interpretation of X-ray diffraction patterns. While
challenging, it is possible to use this paper’s data to uncover the type and
concentration of defects that correspond to a particular choice of processing
conditions. The conjunction of this paper’s data and the XRD patterns could be
used to guide the design of processing protocols to improve defect management
which is a key design feature for the future for perovskites17.

Second, this paper has suggested that there are low-energy pathways that can
self-heal defects. How can we exploit this observation? We hypothesize that
proximal (NN to 3rd NN) defects generated by illumination will readily heal.
It is uncertain what will transpire in more defective systems, putting a
premium on starting with more perfect crystals. This is a testable hypothesis,
experimentally and computationally, which we will pursue in future work.

## Author Contributions

Kumar Miskin : conceptualization, methodology, data curation, visualization
and writing. Yi Cao, Madaline Marland, Farhan Shaikh: data curation,
visualization and writing. Jay Rwaka : data curation. David Moore : funding
acquisition, writing - reviewing and editing. John Marohn : conceptualization,
funding acquisition and writing. Paulette Clancy : conceptualization, funding
acquisition, project administration, resources, supervision and writing.

## Data Availability

All scripts as detailed in SI can be found on the github link :
https://github.com/pclancy-lab/perovskite_defect_study/

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

The authors acknowledge support from the U.S. Department of Energy (DOE) Basic
Energy Sciences award DE-SC0022305. Miskin thanks Johns Hopkins University for
his support in AY 22-23. Computational resource support was provided by the
petascale Hopkins facility, Advanced Research Computing at Hopkins (ARCH)
(rockfish.jhu.edu), supported by National Science Foundation award OAC
1920103. Partial funding for ARCH’s infrastructure was originally provided by
the State of Maryland.

## Notes and references

  * 1 Azpiroz, J. M., Mosconi, E., Bisquert, J., and De Angelis, F. Defect migration in methylammonium lead iodide and its role in perovskite solar cell operation.  Energy and Environmental Science 8, 7 (7 2015), 2118–2127.
  * 2 Balestra, S. R., Vicent-Luna, J. M., Calero, S., Tao, S., and Anta, J. A. Efficient modelling of ion structure and dynamics in inorganic metal halide perovskites.  Journal of Materials Chemistry A 8, 23 (6 2020), 11824–11836.
  * 3 Berne, B. J., Ciccotti, G., and Coker, D. F. Classical and Quantum Dynamics in Condensed Phase Simulations Euroconference on "Technical advances in Particie-based Computational Material Sciences". Coordinator: M. Mareschal LERICI, Villa Marigola Yfe World Scientific.  Tech. rep.
  * 4 Brandt, R. E., Poindexter, J. R., Gorai, P., Kurchin, R. C., Hoye, R. L., Nienhaus, L., Wilson, M. W., Polizzotti, J. A., Sereika, R., Žaltauskas, R., Lee, L. C., Macmanus-Driscoll, J. L., Bawendi, M., Stevanović, V., and Buonassisi, T. Searching for "defect-Tolerant" Photovoltaic Materials: Combined Theoretical and Experimental Screening.  Chemistry of Materials 29, 11 (6 2017), 4667–4674.
  * 5 Evarestov, R. A., Kotomin, E. A., Senocrate, A., Kremer, R. K., and Maier, J. First-principles comparative study of perfect and defective CsPbX3 (X = Br, I) crystals.  Physical Chemistry Chemical Physics 22, 7 (2 2020), 3914–3920.
  * 6 Giannozzi, P., Andreussi, O., Brumme, T., Bunau, O., Buongiorno Nardelli, M., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Cococcioni, M., Colonna, N., Carnimeo, I., Dal Corso, A., de Gironcoli, S., Delugas, P., DiStasio, R. A., Ferretti, A., Floris, A., Fratesi, G., Fugallo, G., Gebauer, R., Gerstmann, U., Giustino, F., Gorni, T., Jia, J., Kawamura, M., Ko, H.-Y., Kokalj, A., Küçükbenli, E., Lazzeri, M., Marsili, M., Marzari, N., Mauri, F., Nguyen, N. L., Nguyen, H.-V., Otero-de-la Roza, A., Paulatto, L., Poncé, S., Rocca, D., Sabatini, R., Santra, B., Schlipf, M., Seitsonen, A. P., Smogunov, A., Timrov, I., Thonhauser, T., Umari, P., Vast, N., Wu, X., and Baroni, S. Advanced capabilities for materials modelling with Quantum ESPRESSO.  Journal of Physics: Condensed Matter 29, 46 (11 2017), 465901.
  * 7 Giannozzi, P., Baroni, S., Bonini, N., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Chiarotti, G. L., Cococcioni, M., Dabo, I., Dal Corso, A., de Gironcoli, S., Fabris, S., Fratesi, G., Gebauer, R., Gerstmann, U., Gougoussis, C., Kokalj, A., Lazzeri, M., Martin-Samos, L., Marzari, N., Mauri, F., Mazzarello, R., Paolini, S., Pasquarello, A., Paulatto, L., Sbraccia, C., Scandolo, S., Sclauzero, G., Seitsonen, A. P., Smogunov, A., Umari, P., and Wentzcovitch, R. M. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials.  Journal of Physics: Condensed Matter 21, 39 (9 2009), 395502.
  * 8 Giannozzi, P., Baseggio, O., Bonfà, P., Brunato, D., Car, R., Carnimeo, I., Cavazzoni, C., de Gironcoli, S., Delugas, P., Ruffino, F. F., Ferretti, A., Marzari, N., Timrov, I., Urru, A., and Baroni, S. Quantum ESPRESSO toward the exascale.
  * 9 Goyal, A., Gorai, P., Peng, H., Lany, S., and Stevanović, V. A computational framework for automation of point defect calculations.  Computational Materials Science 130 (4 2017), 1–9.
  * 10 Henkelman, G., and Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points.  Journal of Chemical Physics 113, 22 (12 2000), 9978–9985.
  * 11 Henkelman, G., Uberuaga, B. P., and Jónsson, H. Climbing image nudged elastic band method for finding saddle points and minimum energy paths.  Journal of Chemical Physics 113, 22 (12 2000), 9901–9904.
  * 12 Heo, S., Seo, G., Lee, Y., Lee, D., Seol, M., Lee, J., Park, J. B., Kim, K., Yun, D. J., Kim, Y. S., Shin, J. K., Ahn, T. K., and Nazeeruddin, M. K. Deep level trapped defect analysis in CH3NH3PbI3 perovskite solar cells by deep level transient spectroscopy.  Energy and Environmental Science 10, 5 (2017), 1128–1133.
  * 13 Kaiser, W., Hussain, K., Singh, A., Alothman, A. A., Meggiolaro, D., Gagliardi, A., Mosconi, E., and De Angelis, F. Defect formation and healing at grain boundaries in lead-halide perovskites.  Journal of Materials Chemistry A 10, 46 (11 2022), 24854–24865.
  * 14 Kang, J., and Wang, L. W. High Defect Tolerance in Lead Halide Perovskite CsPbBr3.  Journal of Physical Chemistry Letters 8, 2 (1 2017), 489–493.
  * 15 Kim, G. Y., Senocrate, A., Wang, Y., Moia, D., and Maier, J. Photo-Effect on Ion Transport in Mixed Cation and Halide Perovskites and Implications for Photo-Demixing**.  Angewandte Chemie International Edition 60, 2 (1 2021), 820–826.
  * 16 Kim, G. Y., Senocrate, A., Yang, T.-Y., Gregori, G., Grätzel, M., and Maier, J. Large tunable photoeffect on ion conduction in halide perovskites and implications for photodecomposition.  Nature Materials 17, 5 (5 2018), 445–449.
  * 17 Lee, Y. M., Maeng, I., Park, J., Song, M., Yun, J.-H., Jung, M.-C., and Nakamura, M. Comprehensive understanding and controlling the defect structures: An effective approach for organic-inorganic hybrid perovskite-based solar-cell application.  Frontiers in Energy Research 6 (2018).
  * 18 López, C. A., Abia, C., Alvarez-Galván, M. C., Hong, B. K., Martínez-Huerta, M. V., Serrano-Sánchez, F., Carrascoso, F., Castellanos-Gómez, A., Fernández-Dĺaz, M. T., and Alonso, J. A. Crystal Structure Features of CsPbBr3 Perovskite Prepared by Mechanochemical Synthesis.  ACS Omega 5, 11 (3 2020), 5931–5938.
  * 19 Manser, J. S., Christians, J. A., and Kamat, P. V. Intriguing Optoelectronic Properties of Metal Halide Perovskites, 11 2016.
  * 20 Masri, P., Harker, A. H., and Stoneham, A. M. Nature and diffusion of the self-interstitial in silicon.  Journal of Physics C: Solid State Physics 16, 18 (6 1983), L613–L616.
  * 21 Meggiolaro, D., and De Angelis, F. First-Principles Modeling of Defects in Lead Halide Perovskites: Best Practices and Open Issues.  ACS Energy Letters 3, 9 (9 2018), 2206–2222.
  * 22 Mills, G., and Jónsson, H. Quantum and thermal effects in ${\mathrm{h}}_{2}$ dissociative adsorption: Evaluation of free energy barriers in multidimensional quantum systems.  Tech. rep., Feb 1994.
  * 23 Minns, J. L., Zajdel, P., Chernyshov, D., Van Beek, W., and Green, M. A. Structure and interstitial iodide migration in hybrid perovskite methylammonium lead iodide.  Nature Communications 8 (2017), 5 p.
  * 24 Moia, D., and Maier, J. Ion Transport, Defect Chemistry, and the Device Physics of Hybrid Perovskite Solar Cells.  ACS Energy Letters 6, 4 (4 2021), 1566–1576.
  * 25 Momma, K., and Izumi, F. VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data.  Journal of Applied Crystallography 44, 6 (12 2011), 1272–1276.
  * 26 Oranskaia, A., Yin, J., Bakr, O. M., Brédas, J. L., and Mohammed, O. F. Halogen Migration in Hybrid Perovskites: The Organic Cation Matters.  Journal of Physical Chemistry Letters 9, 18 (9 2018), 5474–5480.
  * 27 Park, J.-S., Calbo, J., Jung, Y.-K., Whalley, L. D., and Walsh, A. Accumulation of Deep Traps at Grain Boundaries in Halide Perovskites.  ACS Energy Letters 4, 6 (6 2019), 1321–1327.
  * 28 Senocrate, A., and Maier, J. Solid-State Ionics of Hybrid Halide Perovskites.  Journal of the American Chemical Society 141, 21 (5 2019), 8382–8396.
  * 29 Shi, H., and Du, M. H. Shallow halogen vacancies in halide optoelectronic materials.  Physical Review B - Condensed Matter and Materials Physics 90, 17 (11 2014), 174103.
  * 30 Ten Brinck, S., Zaccaria, F., and Infante, I. Defects in lead halide perovskite nanocrystals: Analogies and (many) differences with the bulk.  ACS Energy Letters 4, 11 (11 2019), 2739–2747.
  * 31 Tirmzi, A. M., Christians, J. A., Dwyer, R. P., Moore, D. T., and Marohn, J. A. Substrate-Dependent Photoconductivity Dynamics in a High-Efficiency Hybrid Perovskite Alloy.  The Journal of Physical Chemistry C 123, 6 (2 2019), 3402–3415.
  * 32 Tirmzi, A. M., Dwyer, R. P., Hanrath, T., and Marohn, J. A. Coupled Slow and Fast Charge Dynamics in Cesium Lead Bromide Perovskite.  ACS Energy Letters 2, 2 (2 2017), 488–496.
  * 33 Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism.  Physical Review B 41, 11 (4 1990), 7892–7895.
  * 34 Walsh, A., Scanlon, D. O., Chen, S., Gong, X. G., and Wei, S. Self-Regulation Mechanism for Charged Point Defects in Hybrid Halide Perovskites.  Angewandte Chemie International Edition 54, 6 (2 2015), 1791–1794.
  * 35 Wang, J., Lukose, B., Thompson, M. O., and Clancy, P. Ab initio modeling of vacancies, antisites, and Si dopants in ordered InGaAs.  Journal of Applied Physics 121, 4 (1 2017), 045106.
  * 36 Xue, H., Brocks, G., and Tao, S. First-principles calculations of defects in metal halide perovskites: A performance comparison of density functionals.  Physical Review Materials 5, 12 (12 2021), 125408.
  * 37 Xue, H., Vicent-Luna, J. M., Tao, S., and Brocks, G. Compound Defects in Halide Perovskites: A First-Principles Study of CsPbI3.  Journal of Physical Chemistry C 127, 2 (1 2023), 1189–1197.
  * 38 Yao, Z., Xu, Z., Zhao, W., Zhang, J., Bian, H., Fang, Y., Yang, Y., and Liu, S. Enhanced Efficiency of Inorganic CsPbI3-xBrx Perovskite Solar Cell via Self-Regulation of Antisite Defects.  Advanced Energy Materials 11, 23 (6 2021), 2100403.
  * 39 Yin, W.-J., Shi, T., and Yan, Y. Unusual defect physics in CH3NH3PbI3 perovskite solar cell absorber.  Applied Physics Letters 104, 6 (2 2014).
  * 40 Yuan, Y., and Huang, J. Ion Migration in Organometal Trihalide Perovskite and Its Impact on Photovoltaic Efficiency and Stability.  Accounts of Chemical Research 49, 2 (1 2016), 286–293.