# Why Did This Model Forecast This Future?
Closed-Form Temporal Saliency Towards Causal Explanations of Probabilistic
Forecasts

Chirag Raman Hayley Hung Marco Loog

Delft University of Technology, Delft, The Netherlands
{c.a.raman, h.hung<EMAIL_ADDRESS>
###### Abstract

Forecasting tasks surrounding the dynamics of low-level human behavior are of
significance to multiple research domains. In such settings, methods for
explaining specific forecasts can enable domain experts to gain insights into
the predictive relationships between behaviors. In this work, we introduce and
address the following question: given a probabilistic forecasting model how
can we identify observed windows that the model considers salient when making
its forecasts? We build upon a general definition of information-theoretic
saliency grounded in human perception and extend it to forecasting settings by
leveraging a crucial attribute of the domain: a single observation can result
in multiple valid futures. We propose to express the saliency of an observed
window in terms of the differential entropy of the resulting predicted future
distribution. In contrast to existing methods that either require explicit
training of the saliency mechanism or access to the internal states of the
forecasting model, we obtain a closed-form solution for the saliency map for
commonly used density functions in probabilistic forecasting. We empirically
demonstrate how our framework can recover salient observed windows from head
pose features for the sample task of speaking-turn forecasting using a
synthesized conversation dataset.

##  1 Introduction

The existence of multiple valid futures for a given observed sequence is a
crucial attribute of several forecasting tasks, especially surrounding the
dynamics of low-level human behavior. These tasks include the forecasting of
trajectories of pedestrians [1, 2, 3, 4, 5], vehicles [6, 7, 8, 9], and
autonomous robots [10, 11], or other more general nonverbal cues of humans
[12, 13, 14, 15, 16] and artificial virtual agents [17] in group conversation
settings. Consequently, rather than making point predictions, several machine
learning methods in these settings have attempted to forecast a distribution
over plausible futures [5, 12]. In this work, we introduce and address a
specific research question towards gaining domain-relevant insights into such
forecasts: how can we identify preceding observed sequences that are salient
for a given model towards making probabilistic forecasts over a particular
future window?

Similar questions have been recently posed for models making point forecasts
across various application domains [18, 19, 20]. The broader goal here has
been on designing architectures that learn to produce forecasts that are not
only accurate, but also interpretable. However, in what Lipton [21] terms The
Mythos of Interpretability, the notion of what renders these models
interpretable largely has non-overlapping motivations and attributes across
different works. Examples of these notions of interpretability111Barredo
Arrieta et al. [22] argue for the importance of distinguishing between
interpretability and explainability as different concepts, the latter denoting
any active action or procedure taken by a model with the intent of clarifying
or detailing its internal functions. From this perspective, what the cited
works term interpretability is closer to the notion of explainability. Also
see [22] for an overview of taxonomies of explainability approaches beyond
time-series forecasting. and the associated approaches for forecasting tasks
include: (i) injecting a suitable inductive bias into the model through a set
of basis functions and identifying how they combine to produce an output [18]
(an approach that has recently been applied to the probabilistic setting as
well [23]); (ii) employing a self-attention mechanism to learn temporal
patterns while attending to a common set of features [19]; and (iii) applying
the notion of saliency maps from computer vision [24] to time-series data as a
measure of how much each feature contributes to the final forecast [20].
Irrespective of the notion of interpretability, these methodologies are
underpinned by two common attributes: (a) the interpretability mechanism needs
explicit training as part of the model architecture, and (b) what constitutes
a good explanation is subject to the biases, intuition, or the visual
assessment of the human observer [25, 26]; a phenomenon we refer to as the
interpretation lying in the eye of the beholder. This is especially true for
when saliency maps have been used as tools for post-hoc explanations in
computer vision: the computed map may not measure the intended saliency, and
even be independent of both the model and data generating process Adebayo et
al. [26]. Indeed, things are further confounded by the lack of a common notion
of saliency. As Barredo Arrieta et al. point out, “there is absolutely no
consistency behind what is known as saliency maps, salient masks, heatmaps,
neuron activations, attribution, and other approaches alike” [22, Sec. 5.3].

###  1.1 Revisiting Explanations, and Drawbacks of Existing Methods

Against the backdrop of this research landscape, we begin by revisiting the
question: what constitutes an explanation? Miller [25] recently turned to the
vast body of research in philosophy, pscyhology, and cognitive science on the
topic and highlighted the importance of causality in explanation.
Specifically, in Table 1 we reproduce the simple categorization of explanatory
questions he proposed based on Pearl and Mackenzie’s Ladder of Causation [27].
Within this categorization, Miller argues that the why-questions are the most
challenging, because they use the most sophisticated reasoning and are
contrastive in nature.

Most of the existing interpretability approaches (involving injecting
inductive biases [18, 23] and attention-based mechanisms [19]) can be
considered to be associative within Miller’s categorization, since they reason
about the (unobserved) importance of input features conditioned on one
realization of the future prediction (the event). On the surface, the
perturbation-based saliency methods for interpretable forecasting can be
considered to be counterfactual in approach: Pan et al. [20] learn a saliency
mask that perturbs different parts of the input (simulating alternative
‘causes’) in a way that minimizes an error metric between the corresponding
prediction and ground-truth future (the ‘event’). So while saliency-based
approaches can be promising in answering such why-questions, we identify
several concerns in applying such an approach directly for causal reasoning:
(i) the alternative ‘causes’ generated by noisy perturbations to parts the
input may not correspond to valid real behaviors; (ii) the notion of saliency
in terms of minimizing an error metric is an arbitrary choice, and requires
the ground-truth future for optimizing the saliency mask for a specific
instance; (iii) the saliency map identifies instance-specific attributes
rather than aggregating temporal patterns across all instances; and (iv) the
map needs to be explicitly trained (for every instance in the case of [20]).

Table 1: Classes of Explanatory Question and the Reasoning Required to Answer. Reproduced from Miller [25, Table 3] Question | Reasoning |  Description
---|---|---
What? | Associative |  Reason about which unobserved events could have occurred given the observed events
How? | Interventionist |  Simulate a change in the situation to see if the event still happens
Why? | Counterfactual |  Simulating alternative causes to see whether the event still happens

###  1.2 Our Approach: Building upon a Unifying Framework of Bottom-Up
Saliency

In this work, our central idea for temporal saliency in forecasting settings
is as follows: the saliency of a preceding observed window is related to the
uncertainty over the future window of interest resulting from its observation.
Our framework addresses all the aforementioned concerns we identified with
saliency-based approaches in forecasting as follows: (i) the alternate
‘causes’ in our framework pertain to real preceding windows of actual observed
behavior; (ii) our information-theoretic perspective builds upon a fundamental
definition of saliency grounded in bottom-up preattentive human perception
[28, 29], and can be used to compute saliency for unseen test data where the
ground-truth future is unavailable; (iii) we utilize a distribution over the
futures that the model believes could have occurred over a future window
rather than the single future that did occur: what a model considers plausible
for a single forecast encapsulates the temporal patterns of behavior it has
learned from the entire dataset; and (iv) we derive a closed-form expression
for the saliency map that can be applied to any model that outputs a
probability distribution for its forecasts without additional training or
optimization.

It is reasonable to assume that a model that accurately estimates a
distribution over plausible futures for an input has captured structural
predictive dependencies between features in the data. Our long-term goal is to
enable causal insights using such forecasts. However, the current state of
explainable artificial intelligence is characterized by a lack of consensus
over whether explanations are truly meaningful for real-world applications.
Within this research landscape, rather than directly discovering the
behavioral causes of a forecast, we view the present work as part of a prudent
two-step methodology that incorporates a domain-expert in the loop. First, our
proposed framework here can be used to identify windows salient towards a
model’s forecasts as candidate causes. Second, these candidate causes can be
presented to the domain-expert for subsequent dedicated analysis to
investigate whether the salient relationships between the behaviors within the
observed and future windows are truly causal.

The rest of this paper is organized as follows: in Section 2 we establish the
general framework of bottom-up preattentive saliency which we build upon. In
Section 3 we formalize our proposed framework for temporal saliency in
forecasting settings, and demonstrate its empirical working in Section 4. We
review related literature in Section 5 and conclude with a discussion in
Section 6.

##  2 Background: Information Theoretic Preattentive Saliency in Images

Loog [28] developed a general closed-form expression for saliency using a
surprisal-based operational definition of bottom-up attention from the field
of computational visual perception. Let $L:\mathbb{R}^{n}\to\mathbb{R}^{d}$ be
a general $n$-dimensional $d$-valued image. A continuously differentiable
feature mapping $\phi:\mathbb{R}^{n}\to\mathbb{R}^{N}$ relates every image
location $x\in\mathbb{R}^{n}$ from $L$ to $N$ features. Further, let $p_{X}$
denote the probability density function over image locations $x$ and
$p_{\Phi}$ denote the probability density function over all feature vectors
$\phi(x)$. The distribution $p_{X}$ captures any prior knowledge that would
make one location more salient than another. In the absence of such prior
knowledge, $p_{X}$ is typically chosen to be uniform.

The saliency $S(x)$ of a location $x$ is then defined in terms of the amount
of information—or surprise—of its associated feature vector $\phi(x)$ relative
to the other feature vectors extracted from the same image. The intuition is
that the larger the information content of a certain combination of features,
$-\log p_{\Phi}(\phi(x))$, the more salient the location $x$:

$S(x)>S(x^{\prime})\iff-\log p_{\Phi}(\phi(x))>-\log
p_{\Phi}(\phi(x^{\prime})).$ (1)

This general definition unifies seemingly different definitions of saliency
encountered in the literature. It relates to salient observations being
considered unexpected, rare, or surprising [30, 31, 32, 33]. It also
incorporates the notion of saliency considered in decision-theoretic settings
[33, 34], where bottom-up saliency increases with an increase in the
information associated with feature vectors [35, 36].

Contrary to approaches that determine saliency maps through an explicit data-
driven density estimation [33, 37, 38, 36, 39, 40, 41, 42, 24, 43, 44, 45], a
closed form expression for saliency can be given once the feature mapping
$\phi$ is fixed. When $\phi$ is continuously differentiable, the information
content $-\log p_{\Phi}$ on $\phi(\mathbb{R}^{n})\subset\mathbb{R}^{N}$ over
all image features can be obtained from $\log p_{X}$ through a simple change
of variables [46] from $x$ to $\phi(x)$. The amount of information for the
feature vector $\phi(x)$ at every location $x$, and thereby the saliency
$S(x)$, is then given by the expression:

$\displaystyle-\log p_{\Phi}(\phi(x))$
$\displaystyle=-\log\frac{p_{X}(x)}{\sqrt{\det(J^{t}_{\phi}(x)J_{\phi}(x))}}$
(2) $\displaystyle=-\log
p_{X}(x)+\frac{1}{2}\log\det(J^{t}_{\phi}(x)J_{\phi}(x)),$ (3)

where $J_{\phi}:\mathbb{R}^{n}\to\mathbb{R}^{N\times n}$ denotes the Jacobian
matrix of $\phi$, and $\\_^{t}$ indicates matrix transposition.

A crucial implication of this result is that the computation of saliency can
be performed based on purely local measurements and without the need to refer
to previously observed data or any explicit density estimate of these.
Following Equation 3, Loog simplifies the definition of the saliency map to

$S(x)\coloneqq\det(J^{t}_{\phi}(x)J_{\phi}(x)),$ (4)

given that a monotonic transformation of the saliency map does not essentially
alter the map.

##  3 Closed-Form Temporal Saliency for Social Behavior Forecasting

###  3.1 Setup and Notation

Following the notation from the formulation of the task of Social Cue
Forecasting (SCF) [12], let
$\bm{t}_{\mathrm{obs}}~{}\coloneqq~{}[o1,o2,...,oT]$ denote a window of
consecutively increasing observed timesteps, and
$\bm{t}_{\mathrm{fut}}~{}\coloneqq~{}[f1,f2,...,fT]$ denote an unobserved
future time window, with $f1>oT$. Note that $\bm{t}_{\mathrm{fut}}$ and
$\bm{t}_{\mathrm{obs}}$ can be of different lengths, and
$\bm{t}_{\mathrm{fut}}$ need not immediately follow $\bm{t}_{\mathrm{obs}}$.
Given a set of $n$ interacting participants, let us denote their behavioral
cues over $\bm{t}_{\mathrm{obs}}$ and $\bm{t}_{\mathrm{fut}}$ respectively as
$\bm{X}\coloneqq[\bm{b}^{i}_{t};t\in\bm{t}_{\mathrm{obs}}]_{i=1}^{n}$ and
$\bm{Y}\coloneqq[\bm{b}^{i}_{t};t\in\bm{t}_{\mathrm{fut}}]_{i=1}^{n}$. The
vector $\bm{b}^{i}_{t}$ encapsulates the multimodal cues of interest from
participant $i$ at time $t$, such as head and body pose, facial expressions,
gestures, etc. The task of SCF is to forecast a distribution over possible
future behavioral cues for a given observed sequence $\bm{X}$ of the same
cues, denoted by the probability density function $p_{Y|X}$. Here $Y$ and $X$
denote the multivariate random variables associated with the future and
observed sequences respectively. While the original formulation of SCF dealt
with conversation settings, it is general enough to also apply to other social
settings such as pedestrian trajectory forecasting, where the social cue is
simply the location of an individual. In the rest of this work, for simplicity
we denote the feature array at an individual timestep $t$ in
$\bm{t}_{\mathrm{obs}}$ and $\bm{t}_{\mathrm{fut}}$ as $\bm{X}_{t}$ and
$\bm{Y}_{t}$ respectively, even though they correspond to the same behavioral
cues. The focus of this work is to compute the saliency
$S(\bm{t}_{\mathrm{obs}})$ of an observed $\bm{t}_{\mathrm{obs}}$ towards the
future occurring over a fixed choice of $\bm{t}_{\mathrm{fut}}$.

###  3.2 Defining $\phi$ in Terms of the Uncertainty over the Future Window
$\bm{t}_{\mathrm{fut}}$

The closed-form expression for information-based spatial saliency in Equation
3 makes it explicit that the choice of the feature mapping $\phi$ determines
the actual form of the saliency map. To extend the saliency framework to
forecasting settings, we begin with the following intuition: the saliency of
an observed window $\bm{t}_{\mathrm{obs}}$ is directly related to how certain
or uncertain it makes the resulting future over $\bm{t}_{\mathrm{fut}}$. This
idea is related to the notion of surprisal; an observation that changes the
certainty of the future is surprising, and thereby salient. We formalize this
intuition by mapping the window $\bm{t}_{\mathrm{obs}}$ to the differential
entropy of the predicted future distribution over $\bm{t}_{\mathrm{fut}}$,
conditioned on the observed features $\bm{X}$. That is, we define
$\phi:\bm{t}_{\mathrm{obs}}\mapsto h(Y|X=\bm{X})$, where the conditional
differential entropy of $Y$ given $\\{X=\bm{X}\\}$ is defined as

$h(Y|X=\bm{X})\coloneqq-\int p_{Y|X}(\bm{Y}|\bm{X})\log
p_{Y|X}(\bm{Y}|\bm{X})d\bm{Y}.$ (5)

Following Equation 5, our framework for computing the saliency map is as
follows. First we identify a specific $\bm{t}_{\mathrm{fut}}$ of interest.
This could correspond to a semantic event of forecasting interest into which a
domain-expert seeks to get some insight, such as a speaking-turn change [47,
48] an interaction termination [49, 50], or a synchronous behavior event [51].
Looking back in time before $\bm{t}_{\mathrm{fut}}$, we then compute
$h(Y|X=\bm{X})$ for different observed behaviors $\bm{X}$ corresponding to
different locations of a sliding $\bm{t}_{\mathrm{obs}}$ within a look-back
period determined to be reasonable for the specific application. The computed
differential entropy values are then inserted into Equation 4 to obtain the
saliency of the different behaviors over the $\bm{t}_{\mathrm{obs}}$ locations
towards forecasting the future over the chosen $\bm{t}_{\mathrm{fut}}$. The
general approach is summarized in Algorithm 1.

Algorithm 1 Temporal Saliency in Probabilistic Forecasting

1:The probability density function $p_{Y|X}$

2:A fixed $\bm{t}_{\mathrm{fut}}$ of interest, a sequence of $m$ preceding
observed windows
$O=[\bm{t}_{\mathrm{obs}}^{1},\ldots,\bm{t}_{\mathrm{obs}}^{m}]$, and the
behavioral features $\bm{X}^{j}$ for every $\bm{t}_{\mathrm{obs}}^{j}$

3:The saliency map $S(O)$ over the observed windows

4:

5:for each $\bm{t}_{\mathrm{obs}}^{j}\in O$ do

6: Compute the feature mapping $\phi(\bm{t}_{\mathrm{obs}}^{j})\leftarrow
h(Y|X=\bm{X}^{j})$

7:end for

8:

9:Compute the saliency map
$S(\bm{t}_{\mathrm{obs}})\leftarrow\det(J^{t}_{\phi}(\bm{t}_{\mathrm{obs}})J_{\phi}(\bm{t}_{\mathrm{obs}}))$

Differential entropy possesses certain favorable properties that make it a
suitable choice as $\phi$ for computing the saliency map in forecasting
settings. To begin, the scale of the forecast density does not affect the
resulting saliency map [see 52, Theorem 8.6.4]:

$\displaystyle h(aY)$ $\displaystyle=h(Y)+\log|a|,\text{for }a\neq
0,\text{and}$ (6) $\displaystyle h(\bm{A}Y)$
$\displaystyle=h(Y)+\log|\det(\bm{A})|,\text{when $\bm{A}$ is a square
matrix}.$ (7)

That is, scaling the distribution changes the differential entropy by only a
constant factor. So the saliency map resulting from inserting the entropy into
Equation 4 remains unaffected since the Jacobian term only depends on the
relative change in entropy across different choices of
$\bm{t}_{\mathrm{obs}}$. Similarly, translating the predicted density leaves
the saliency map unaffected [see 52, Theorem 8.6.3]:

$h(Y+c)=h(Y).$ (8)

###  3.3 Computing $h(Y|X=\bm{X})$

To compute the differential entropy of the future distribution, we need the
density function $p_{Y|X}$. Typically, this density is modeled as a
multivariate Gaussian distribution [12, 53, 54, 55]. When the decoding of the
future is non-autoregressive, the parameters of the distributions for all
$t\in\bm{t}_{\mathrm{fut}}$ are estimated in parallel. In these cases, the
differential entropy has a closed-form expression. Assuming a $d$-dimensional
predicted Gaussian distribution with mean $\bm{\mu}$ and covariance matrix
$\bm{K}$, the expression for the entropy of a multivariate Gaussian
distribution is given by [see 52, Theorem 8.4.1]

$h(Y|X=\bm{X})=h(\mathcal{N}_{d}(\bm{\mu},\bm{K}))=\frac{1}{2}\log[(2\pi
e)^{d}\det(\bm{K})].$ (9)

When $\bm{K}$ is diagonal, so that the predicted distribution is factorized
over participants and features, we can simply sum the $\log$ of the individual
variances to obtain the feature mapping $\phi$. Note that from Equation 9, for
a multivariate Gaussian distribution, the differential entropy only depends on
the covariance, or the spread of the distribution, aligning with the notion of
differential entropy as a measure of total uncertainty. (See 52, Tab. 17.1; 56
for closed-form expressions for a large number of commonly employed
probability density functions.)

Another common approach for inferring the future density function is the use
of probabilistic autoregressive decoders [12, 1, 57, 54]. Here, one possible
decoding approach [57, 1] involves taking a specific sample
$\widehat{\bm{Y}}_{t}$ from the density predicted at each
$t\in\bm{t}_{\mathrm{fut}}$, and passing it back as input to the decoder for
the estimation of the density at timestep $t+1$. Therefore, the density at
$t+1$ depends on the randomness introduced in sampling $\widehat{\bm{Y}}_{t}$.
This yields an arbitrarily complex form for the joint distribution $p_{Y|X}$
(over all $t\in\bm{t}_{\mathrm{fut}}$). Figure 1 illustrates the concept for
two timesteps. Here, a single forecast would only output the shaded red
distribution for $Y_{2}$. This precludes the direct computation of
$h(Y_{1},Y_{2})$ that requires the full joint distribution $p_{Y_{1},Y_{2}}$.
In such cases, we have two broad options: using a simplifying assumption to
retain computational simplicity, or approximating the differential entropy by
sampling.

The simpler option is to redefine our feature-mapping as
$\phi:\bm{t}_{\mathrm{obs}}\mapsto\sum_{t\in\bm{t}_{\mathrm{fut}}}h(Y_{t}|\widehat{\bm{Y}}_{<t},\bm{X})$,
i.e. we sum the differential entropy of the individual densities estimated at
each timestep as an approximation of the total uncertainty over the predicted
sequence. Note that following the chain rule for differential entropy ([see
52, Eq. 8.62]), the joint entropy can indeed be written as the sum of
individual conditionals. However, in general,

$h(Y|X=\bm{X})=\sum_{t\in\bm{t}_{\mathrm{fut}}}h(Y_{t}|Y_{<t},\bm{X})\neq\sum_{t\in\bm{t}_{\mathrm{fut}}}h(Y_{t}|\widehat{\bm{Y}}_{<t},\bm{X}).$
(10)

Figure 1: Illustrating predicted densities under the greedy autoregressive
decoding approach for two timesteps. For simplicity, we depict a joint
Gaussian distribution and omit the conditioning on $\bm{X}$ everywhere.

Employing this approximation (of summing the entropies across all timesteps)
relies on the observation that for autoregressive decoding, the parameters of
the predicted distribution for $Y_{t}$ is computed as a deterministic function
of the decoder hidden state. That is, $Y_{t}$ is conditionally independent of
$Y_{<t}$ given the hidden state of the decoder $\bm{s}_{t}$ at timestep $t$.
The underlying assumption is that for a well-trained decoder, $\bm{s}_{t}$
encodes all relevant information from other timesteps to infer the
distribution of $Y_{t}$. So at inference, despite being a function of the
single sample $\widehat{\bm{Y}}_{t-1}$, the predicted distribution conditioned
on $\bm{s}_{t}$ provides a reasonable estimate of the the uncertainty in
$Y_{t}$. The benefit of employing this assumption is that when each $Y_{t}$ is
modeled using a density function that has a closed-form expression for
differential entropy 52, Tab. 17.1; 56, every item in the sum can be computed
anlaytically, and we obtain a closed-form expression for the saliency map.
Apart from a multivariate Gaussian, the other common choice for modeling
$Y_{t}$ is using a Gaussian mixture [57, 1]. While a closed-form expression
for the differential entropy of a Gaussian mixture is not known,
approximations that approach the true differential entropy can be obtained
efficiently [58, 59, 60] to directly compute the feature mapping $\phi$.

What if the differential entropy for $Y_{t}$ does not have an analytical
expression or a computationally efficient approximation? In such cases,
$h(Y|X=\bm{X})$ can be estimated using sampling or other non-parametric
approaches [61, 62, 63, 64]. While computationally more expensive than
parametric methods, such sampling-based methods can provide approximations
that converges to the true entropy.

Another practical consideration arises from the fact that we consider
different candidate $\bm{t}_{\mathrm{obs}}$ to evaluate how their observation
changes the uncertainty over a fixed $\bm{t}_{\mathrm{fut}}$. So
$\bm{t}_{\mathrm{fut}}$ need not immediately follow $\bm{t}_{\mathrm{obs}}$.
As Raman et al. [12] point out, a model for forecasting human behavior ought
to support such a requirement given that domain experts are often interested
in behavioral phenomena that occur after an arbitrary offset [see 12, Sec.5].
In such cases, the typical approach which starts decoding immediately after
the end of $\bm{t}_{\mathrm{obs}}$ would entail discarding the predictions in
the gap between $\bm{t}_{\mathrm{obs}}$ and $\bm{t}_{\mathrm{fut}}$. This
could cascade prediction errors and lead to worse estimates of $p_{Y|X}$. One
approach to avoid such errors in practice is for the underlying model to
encode the offset $\Delta t=f1-oT$ as an input and commence decoding at
timestep $f1$ [12].

##  4 Empirical Demonstration Using Synthesized Variance over Futures

Figure 2: Illustrating the synthetic conversation dynamics dataset. Speakers
are denoted in orange and listeners in green. For a fixed
$\bm{t}_{\mathrm{fut}}$ we depict two preceding $\bm{t}_{\mathrm{obs}}$
windows. By construction, when observing a stable speaking turn over
$\bm{t}_{\mathrm{obs}}^{1}$, two valid futures are possible over
$\bm{t}_{\mathrm{fut}}$. These correspond to a turn handover to the immediate
left or right of the current speaker. Over $\bm{t}_{\mathrm{obs}}^{2}$, when a
listener nods to indicate the desire to take the floor, the future over
$\bm{t}_{\mathrm{fut}}$ becomes certain, corresponding to the listener
successfully taking over the speaking turn. Here $\bm{t}_{\mathrm{obs}}^{2}$
is consequently more salient than $\bm{t}_{\mathrm{obs}}^{1}$ towards
forecasting the turn change over $\bm{t}_{\mathrm{fut}}$.

The best a well-trained model can do is to predict the true variance over
plausible futures given an observation. In this section we use a synthetic
dataset of simplified conversation dynamics to demonstrate how our framework
empirically retrieves the windows that are salient by construction. Note that
the model does not need to have any conception of saliency, and as such our
framework doesn’t involve any training using a top-down notion of saliency.
The goal is to illustrate how our framework retrieves what would constitute as
bottom-up salient to a well-trained forecasting model.

###  4.1 The Dataset

For illustration, we choose a sample forecasting task that has received
extensive domain interest over the last decade: forecasting a turn change in a
multiparty conversation [65, 66, 67, 68, 47, 69, 70]. We synthesize a group
conversation following established patterns of social behavior. First, the
visual focus of attention of listeners is usually the speaker, while the
speaker might look at different listeners over the speaking turn [71]. Second,
head gestures and gaze patterns are predictive of the next speaker [65, 66,
70, 69]. We implement a simplified version of these dynamics as follows. The
speaker rotates towards the geometric center of the formation after acquiring
a speaking turn, while the listeners orient towards the speaker. A listener
nods their head to indicate a desire to acquire the floor, following which the
current speaker rotates to look at the listener and hands over the speaking
turn. We represent head pose with quaternions, given their common use in
representing human motion and pose [12, 72]. We simulate the turn changes to
occur once clockwise and once anticlockwise around the group so that each
participant yields the floor once each to the participant on their immediate
left and right. We provide the dataset and an animated visualization in the
Suppementary material.

Figure 3: Computing Saliency. We plot the quaternion dimensions qx and qy for
the listener that nods over $\bm{t}_{\mathrm{obs}}^{2}$ in Figure 2 (top). The
observation of the head-nod beginning in frame $138$ makes the future over
frames $183-228$ ($\bm{t}_{\mathrm{fut}}$) certain. This is reflected in the
reduction in the mean entropy over future values of all participants (middle).
The saliency map obtained using our framework (bottom) identifies the
beginning of the head-nod in frame $138$ as the salient timestep towards
forecasting the future over $\bm{t}_{\mathrm{fut}}$.

###  4.2 Computing Temporal Saliency

By construction, it is impossible to guess the future with certainty by simply
observing the current speaker speaking: there exist two valid examples of the
future head behavior in the dataset for such a sequence. Only the commencement
of a head-nod by a listener makes the future turn handover certain. Figure 2
illustrates this mechanism. Therefore, the observation of a head-nod is
salient towards forecasting a turn change in this dataset. Can our framework
identify the head-nod as salient? Let us assume a perfect model. Such a model
would predict the true variance over the two valid future quaternion
trajectories given a stable speaker turn. This variance vanishes when the
future becomes certain as soon as one of the listeners start nodding. We model
the future distribution using a Gaussian function for simplicity (setting std.
to $10^{-10}$ for the single future), but a more complex distribution that
predicts the appropriate change in variance would also work in practice.

We now implement Algorithm 1 as follows. We identify a window where a turn
change occurs in the data (frames $183$-$228$) and denote this $45$ frame
window as the $\bm{t}_{\mathrm{fut}}$ of interest. While we manually identify
an interesting event for illustration, such a window could also correspond to
an interesting prediction by a model. We generate a set of candidate
$\bm{t}_{\mathrm{obs}}$ by sliding a $30$ frame window over a horizon of $100$
frames prior to $\bm{t}_{\mathrm{fut}}$, with a stride of $1$ frame. For every
observed $\bm{t}_{\mathrm{obs}}$, we fit a Gaussian density to the quaternion
and speaking status features of all participants over the futures that can
occur during $\bm{t}_{\mathrm{fut}}$. We then set the entropy of this Gaussian
density as the feature $\phi$ for the $\bm{t}_{\mathrm{obs}}$. We obtain the
saliency map using Equation 4, and plot the listener features, entropy, and
saliency across timesteps in Figure 3. We observe that our framework
attributes higher saliency to the timesteps where the listener first begins
the head-nod in frame $138$. Once the nod is in motion, the saliency drops as
expected since the future is then already certain given the data.

##  5 Related Work: Explainable Methods for Time-Series Tasks Across Domains

The larger focus of explainability techniques involving time-series data has
been on the task of classifying time-series. The goal has been to attribute
relevance, to each input feature at a given timestep. Here saliency approaches
often overlap with techniques developed for image data, and can be categorized
into gradient-based, perturbation-based, and attention-based techniques. The
broad approach of gradient-based techniques involves evaluating the gradient
of the output class with respect to the input [73]. Several variants of this
idea have been proposed [74, 75, 76, 77, 78]. The idea behind perturbation-
based techniques is to examine how the output changes in response to some
perturbation of the input. Here, perturbations can be implemented by either
occluding contiguous regions of the input [79, 80]; performing an ablation of
the features [81]; or randomly permuting features [82]. Ismail et al. [83]
provide a benchmark of a subset of these techniques. Attention-based methods
incorporate an attention mechanism into the model that is trained to attribute
importance to previous parts of the input towards an output at each timestep.
Such techniques have been especially employed for healthcare data, with early
methods applying a reverse-time attention [84], and later methods applying the
attention to probabilistic state-space representations [85].

Some of these broad ideas have been applied to the regression setting, to make
interpretable forecasts of future values of the time-series features (see
Section 1). Lim et al. [19] proposed an attention-based architecture that
employs recurrent layers for local processing and interpretable self-attention
layers for capturing long-term dependencies. Pan et al. [20] recently proposed
a computing saliency as a mixup strategy with a learnable mask between series
images and their perturbed version. They view saliency in terms of minimizing
the mean squared error between the predictions and ground-truths for a
particular instance. Focusing on the univariate point-forecasting problem,
Oreshkin et al. [18] proposed injecting inductive biases by computing the
forecast as a combination of a trend and seasonality model. They argue that
this decomposition make the outputs more interpretable.

Developing explainable techniques for the probabilistic forecasting setting
remains largely unexplored. Rügamer et al. [23] transform the forecast using
an additive predictor using predefined basis functions such as Bernstein
polynomials. They relate interpretability to the coefficients of these basis
functions (a notion similar to that of Oreshkin et al. [18]). Panja et al.
[86] embed the classical linear ARIMA model into a non-linear autoregressive
neural network for univariate probabilistic forecasting. As before, the
explainability here also stems from the ‘white-box’ nature of the linear ARIMA
component. Li et al. [87] propose an automatic relevance determination network
to identify useful exogenous variables (i.e. variables that can affect the
forecast without being a part of the time-series data). To the best of our
knowledge, saliency based methods have not yet been considered within this
setting.

##  6 Discussion: Operationalizing Saliency for Forecasting and Implications
for Interdisciplinary Research

We have proposed a framework for computing the saliency of observed sequences
in probabilistic forecasting tasks, that results in a convenient closed-form
expression of the saliency map for commonly used probability density functions
to represent forecasts [57, 12, 1, 2]. Rather than defining saliency in a top-
down manner as a function of some task-specific error metric, we have started
from a more fundamental conception of bottom-up, or task-agnostic, saliency.
loogInformationTheoreticPreattentive2011’s original definition pertains to
preattentive saliency, which captures what is perceived to be subconsciously
informative before conscious (attentive) processing by the brain. Here, a
surprising or unexpected observation is salient. For instance, in a large
white image with a single black pixel, the black pixel is salient. The direct
application of this concept to time-series data would involve identifying
surprising task-agnostic temporal events. For instance, imagine viewing a
static landscape where a bird suddenly flies in. The entry of the bird into
the scene is unexpected, and therefore salient.

When applied to forecasting tasks, however, this idea of surprisal (or
unexpectedness or informativeness) that saliency represents needs to be tied
to the future outcome. The saliency computed by most methods working on point-
forecasting tasks deals with which past features are surprising given a
specific realization of the future. While not explicitly stated by these
works, we argue that this notion of saliency is related to the surprisal in
$p_{X|Y}$ for some specific $\bm{Y}$. We therefore interpret these methods as
being associative in nature within millerExplanationArtificialIntelligence2019
categorization in Section 1. In contrast, our approach is counterfactual
because we examine alternate future outcomes, while conceptualizing saliency
more naturally defined in terms of the changes in the uncertainty in $p_{Y|X}$
in response to different realizations of observed sequences. However, rather
than corresponding to random occlusions or perturbations of the input, the
different realizations of $X$ in our framework correspond to real features (or
behaviors) preceding a future, which is more suitable to present to domain
experts as candidate causes.

However, are we really capturing the salient predictive relationships in the
data? The distribution $p_{Y|X}$ is considered to capture the structural
predictive relationships between different features across the entire data. In
principle, when it is possible to have access to the true $p_{Y|X}$ as in
Section 4, the salient sequences identified by our framework reflect the true
predictive relationships between the behaviors contained in
$\bm{t}_{\mathrm{obs}}$ and $\bm{t}_{\mathrm{fut}}$. However, estimating this
density analytically would entail identifying the multiple futures in the data
corresponding to every occurrence of the same observed sequence of features.
In practice, subtle variations in behaviors and sensor measurement errors make
it infeasible to estimate $p_{Y|X}$ analytically, so a model is trained to
learn generalized patterns from the given data. In these cases, our framework
identifies the sequences that the model considers salient for its forecasts
given the data. Consequently, despite being grounded in a conceptually sound
definition of saliency, the measures obtained by our framework may not reflect
real-world truths about behavioral patterns. This is especially crucial to
consider in the healthcare and human behavior domains to avoid potential
prejudices against certain behaviors, or worse, misdiagnoses of affective
conditions. We therefore reiterate that the salient sequences retrieved by our
framework ought to be treated as candidate causes until subsequently examined
along with a domain-expert in the loop.

## Acknowledgments and Disclosure of Funding

This research was partially funded by the Netherlands Organization for
Scientific Research (NWO) under the MINGLE project number 639.022.606. Chirag
would like to thank Amelia Villegas-Morcillo and Yeshwanth Napolean for their
input on the manuscript, and the members of the TUDelft Pattern Recognition
Lab for the thoughtful discussions.

## References

  * Salzmann et al. [2021] Tim Salzmann, Boris Ivanovic, Punarjay Chakravarty, and Marco Pavone.  Trajectron++: Dynamically-Feasible Trajectory Forecasting With Heterogeneous Data.  _arXiv:2001.03093 [cs]_ , January 2021.
  * Rudenko et al. [2020] Andrey Rudenko, Luigi Palmieri, Michael Herman, Kris M Kitani, Dariu M Gavrila, and Kai O Arras.  Human motion trajectory prediction: A survey.  _The International Journal of Robotics Research_ , 39(8):895–935, 2020.
  * Huang et al. [2019] Yingfan Huang, Huikun Bi, Zhaoxin Li, Tianlu Mao, and Zhaoqi Wang.  STGAT: Modeling Spatial-Temporal Interactions for Human Trajectory Prediction.  In _2019 IEEE/CVF International Conference on Computer Vision (ICCV)_ , pages 6271–6280, Seoul, Korea (South), October 2019. IEEE.  ISBN 978-1-72814-803-8.  doi: 10.1109/ICCV.2019.00637.
  * Zhang et al. [2019] Pu Zhang, Wanli Ouyang, Pengfei Zhang, Jianru Xue, and Nanning Zheng.  SR-LSTM: State Refinement for LSTM towards Pedestrian Trajectory Prediction.  _arXiv:1903.02793 [cs]_ , March 2019.
  * Mohamed et al. [2020] Abduallah Mohamed, Kun Qian, Mohamed Elhoseiny, and Christian Claudel.  Social-STGCNN: A Social Spatio-Temporal Graph Convolutional Neural Network for Human Trajectory Prediction.  _arXiv:2002.11927 [cs]_ , February 2020.
  * Gilles et al. [2022] Thomas Gilles, Stefano Sabatini, Dzmitry Tsishkou, Bogdan Stanciulescu, and Fabien Moutarde.  THOMAS: Trajectory Heatmap Output with learned Multi-Agent Sampling.  _arXiv:2110.06607 [cs]_ , January 2022.
  * Zhao et al. [2020] Hang Zhao, Jiyang Gao, Tian Lan, Chen Sun, Benjamin Sapp, Balakrishnan Varadarajan, Yue Shen, Yi Shen, Yuning Chai, Cordelia Schmid, Congcong Li, and Dragomir Anguelov.  TNT: Target-driveN Trajectory Prediction.  _arXiv:2008.08294 [cs]_ , August 2020.
  * Carrasco et al. [2021] Sandra Carrasco, David Fernández Llorca, and Miguel Ángel Sotelo.  SCOUT: Socially-COnsistent and UndersTandable Graph Attention Network for Trajectory Prediction of Vehicles and VRUs.  _arXiv:2102.06361 [cs]_ , May 2021.
  * Zeng et al. [2020] Wenyuan Zeng, Shenlong Wang, Renjie Liao, Yun Chen, Bin Yang, and Raquel Urtasun.  DSDNet: Deep Structured self-Driving Network.  _arXiv:2008.06041 [cs]_ , August 2020.
  * Vemula et al. [2017] Anirudh Vemula, Katharina Muelling, and Jean Oh.  Modeling cooperative navigation in dense human crowds.  In _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 1685–1692. IEEE, 2017.
  * Ivanovic et al. [2021] Boris Ivanovic, Yifeng Lin, Shubham Shrivastava, Punarjay Chakravarty, and Marco Pavone.  Propagating state uncertainty through trajectory forecasting.  _arXiv preprint arXiv:2110.03267_ , 2021.
  * Raman et al. [2021] Chirag Raman, Hayley Hung, and Marco Loog.  Social Processes: Self-Supervised Meta-Learning over Conversational Groups for Forecasting Nonverbal Social Cues.  _arXiv:2107.13576 [cs]_ , July 2021.
  * Adeli et al. [2020] Vida Adeli, Ehsan Adeli, Ian Reid, Juan Carlos Niebles, and Hamid Rezatofighi.  Socially and Contextually Aware Human Motion and Pose Forecasting.  _arXiv:2007.06843 [cs]_ , July 2020.
  * Nguyen and Celiktutan [2022] Tan Viet Tuyen Nguyen and Oya Celiktutan.  Context-Aware Human Behaviour Forecasting in Dyadic Interactions.  _Understanding Social Behavior in Dyadic and Small Group Interactions_ , pages 1–19, 2022.
  * Barquero et al. [2022] German Barquero, Johnny Núñez, Sergio Escalera, Zhen Xu, Wei-Wei Tu, Isabelle Guyon, and Cristina Palmero.  Didn’t see that coming: A survey on non-verbal social human behavior forecasting.  _arXiv:2203.02480 [cs]_ , March 2022.
  * Yao et al. [2018] Taiping Yao, Minsi Wang, Bingbing Ni, Huawei Wei, and Xiaokang Yang.  Multiple Granularity Group Interaction Prediction.  In _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2246–2254, Salt Lake City, UT, June 2018. IEEE.  ISBN 978-1-5386-6420-9.  doi: 10.1109/CVPR.2018.00239.
  * Ahuja et al. [2019] Chaitanya Ahuja, Shugao Ma, Louis-Philippe Morency, and Yaser Sheikh.  To React or not to React: End-to-End Visual Pose Forecasting for Personalized Avatar during Dyadic Conversations.  _arXiv:1910.02181 [cs]_ , October 2019.
  * Oreshkin et al. [2020] Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  N-BEATS: Neural basis expansion analysis for interpretable time series forecasting.  _arXiv:1905.10437 [cs, stat]_ , February 2020.
  * Lim et al. [2020] Bryan Lim, Sercan O. Arik, Nicolas Loeff, and Tomas Pfister.  Temporal Fusion Transformers for Interpretable Multi-horizon Time Series Forecasting.  _arXiv:1912.09363 [cs, stat]_ , September 2020.
  * Pan et al. [2021] Qingyi Pan, Wenbo Hu, and Ning Chen.  Two Birds with One Stone: Series Saliency for Accurate and Interpretable Multivariate Time Series Forecasting.  In _Twenty-Ninth International Joint Conference on Artificial Intelligence_ , volume 3, pages 2884–2891, August 2021.  doi: 10.24963/ijcai.2021/397.
  * Lipton [2017] Zachary C. Lipton.  The Mythos of Model Interpretability, March 2017.
  * Barredo Arrieta et al. [2020] Alejandro Barredo Arrieta, Natalia Díaz-Rodríguez, Javier Del Ser, Adrien Bennetot, Siham Tabik, Alberto Barbado, Salvador Garcia, Sergio Gil-Lopez, Daniel Molina, Richard Benjamins, Raja Chatila, and Francisco Herrera.  Explainable artificial intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible AI.  _Information Fusion_ , 58:82–115, 2020.  ISSN 1566-2535.  doi: 10.1016/j.inffus.2019.12.012.
  * Rügamer et al. [2022] David Rügamer, Philipp F. M. Baumann, Thomas Kneib, and Torsten Hothorn.  Probabilistic Time Series Forecasts with Autoregressive Transformation Models.  _arXiv:2110.08248 [cs]_ , February 2022.
  * Dabkowski and Gal [2017] Piotr Dabkowski and Yarin Gal.  Real time image saliency for black box classifiers.  _Advances in neural information processing systems_ , 30, 2017.
  * Miller [2019] Tim Miller.  Explanation in artificial intelligence: Insights from the social sciences.  _Artificial Intelligence_ , 267:1–38, February 2019.  ISSN 0004-3702.  doi: 10.1016/j.artint.2018.07.007.
  * Adebayo et al. [2018] Julius Adebayo, Justin Gilmer, Michael Muelly, Ian Goodfellow, Moritz Hardt, and Been Kim.  Sanity checks for saliency maps.  _Advances in neural information processing systems_ , 31, 2018.
  * Pearl and Mackenzie [2018] Judea Pearl and Dana Mackenzie.  _The Book of Why: The New Science of Cause and Effect_.  Basic books, 2018.
  * Loog [2011] Marco Loog.  Information theoretic preattentive saliency: A closed-form solution.  In _2011 IEEE International Conference on Computer Vision Workshops (ICCV Workshops)_ , pages 1418–1424, Barcelona, Spain, November 2011. IEEE.  ISBN 978-1-4673-0063-6 978-1-4673-0062-9 978-1-4673-0061-2.  doi: 10.1109/ICCVW.2011.6130417.
  * vanderHeijden [1996] AHC vanderHeijden.  Perception for selection, selection for action, and action for perception, 1996.
  * Walker et al. [1998] Kevin N Walker, Timothy F Cootes, and Christopher J Taylor.  Locating salient object features.  In _BMVC_ , volume 98, pages 557–566. Citeseer, 1998.
  * Garcia et al. [2001] Jose A. Garcia, Joaquín Fdez-Valdivia, Xose R Fdez-Vidal, and Rosa Rodriguez-Sanchez.  Information theoretic measure for visual target distinctness.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 23(4):362–383, 2001.
  * Itti and Koch [2001] Laurent Itti and Christof Koch.  Computational modelling of visual attention.  _Nature reviews neuroscience_ , 2(3):194–203, 2001.
  * Torralba [2003] Antonio Torralba.  Modeling global scene factors in attention.  _JOSA A_ , 20(7):1407–1418, 2003.
  * Lingyun et al. [2007] Zhang Lingyun, Matthew H Tong, and Garrison W Cottrell.  Information attracts attention: A probabilistic account of the cross-race advantage in visual search.  In _Proceedings of the Annual Meeting of the Cognitive Science Society_ , volume 29, 2007.
  * Rosenholtz [1999] Ruth Rosenholtz.  A simple saliency model predicts a number of motion popout phenomena.  _Vision research_ , 39(19):3157–3163, 1999.
  * Gao and Vasconcelos [2007] Dashan Gao and Nuno Vasconcelos.  Bottom-up saliency is a discriminant process.  In _2007 IEEE 11th International Conference on Computer Vision_ , pages 1–6. IEEE, 2007.
  * Bruce [2005] Neil DB Bruce.  Features that draw visual attention: An information theoretic perspective.  _Neurocomputing_ , 65:125–133, 2005.
  * Itti and Baldi [2005] Laurent Itti and Pierre Baldi.  A principled approach to detecting surprising events in video.  In _2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05)_ , volume 1, pages 631–637. IEEE, 2005\.
  * Li et al. [2013] Xiaohui Li, Huchuan Lu, Lihe Zhang, Xiang Ruan, and Ming-Hsuan Yang.  Saliency detection via dense and sparse reconstruction.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pages 2976–2983, 2013.
  * Simonyan et al. [2013] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Deep inside convolutional networks: Visualising image classification models and saliency maps.  _arXiv preprint arXiv:1312.6034_ , 2013.
  * Li and Yu [2016] Guanbin Li and Yizhou Yu.  Visual saliency detection based on multiscale deep CNN features.  _IEEE transactions on image processing_ , 25(11):5012–5024, 2016.
  * Zhou et al. [2016] Bolei Zhou, Aditya Khosla, Agata Lapedriza, Aude Oliva, and Antonio Torralba.  Learning deep features for discriminative localization.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 2921–2929, 2016.
  * Fong and Vedaldi [2017] Ruth C Fong and Andrea Vedaldi.  Interpretable explanations of black boxes by meaningful perturbation.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pages 3429–3437, 2017.
  * Jiang et al. [2019] Lai Jiang, Zhe Wang, Mai Xu, and Zulin Wang.  Image saliency prediction in transformed domain: A deep complex neural network method.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 8521–8528, 2019.
  * Huang and Gao [2020] Kun Huang and Shenghua Gao.  Image saliency detection via multi-scale iterative CNN.  _The Visual Computer_ , 36(7):1355–1367, 2020\.
  * Boothby [1975] WM Boothby.  An introduction to differential geometry and riemannian manifolds, 1975\.
  * Keitel et al. [2015] Anne Keitel, Moritz M Daum, et al.  The use of intonation for turn anticipation in observed conversations without visual signals as source of information.  _Frontiers in psychology_ , 6(108):265–273, 2015\.
  * Levinson and Torreira [2015] Stephen C. Levinson and Francisco Torreira.  Timing in turn-taking and its implications for processing models of language.  _Frontiers in Psychology_ , 6, June 2015.  ISSN 1664-1078.  doi: 10.3389/fpsyg.2015.00731.
  * Bohus and Horvitz [2014] Dan Bohus and Eric Horvitz.  Managing Human-Robot Engagement with Forecasts and… um… Hesitations.  _Proceedings of the 16th International Conference on Multimodal Interaction_ , page 8, 2014.
  * van Doorn [2018] Felix van Doorn.  Rituals of Leaving: Predictive Modelling of Leaving Behaviour in Conversation.  _Master of Science Thesis, Delft University of Technology_ , 2018\.
  * Bilakhia et al. [2013] Sanjay Bilakhia, Stavros Petridis, and Maja Pantic.  Audiovisual Detection of Behavioural Mimicry.  In _2013 Humaine Association Conference on Affective Computing and Intelligent Interaction_ , pages 123–128, Geneva, Switzerland, September 2013. IEEE.  ISBN 978-0-7695-5048-0.  doi: 10.1109/ACII.2013.27.
  * [52] Thomas M Cover and Joy A Thomas.  _Elements of Information Theory_.
  * Salinas et al. [2019a] David Salinas, Michael Bohlke-Schneider, Laurent Callot, Roberto Medico, and Jan Gasthaus.  High-dimensional multivariate forecasting with low-rank Gaussian Copula Processes.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019a.
  * Salinas et al. [2019b] David Salinas, Valentin Flunkert, and Jan Gasthaus.  DeepAR: Probabilistic Forecasting with Autoregressive Recurrent Networks.  _arXiv:1704.04110 [cs, stat]_ , February 2019b.
  * Rangapuram et al. [2018] Syama Sundar Rangapuram, Matthias W Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, and Tim Januschowski.  Deep state space models for time series forecasting.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.
  * Lazo and Rathie [1978] A Verdugo Lazo and P Rathie.  On the entropy of continuous probability distributions (corresp.).  _IEEE Transactions on Information Theory_ , 24(1):120–122, 1978.
  * Ha and Eck [2017] David Ha and Douglas Eck.  A Neural Representation of Sketch Drawings.  _arXiv:1704.03477 [cs, stat]_ , May 2017.
  * Huber et al. [2008] Marco F. Huber, Tim Bailey, Hugh Durrant-Whyte, and Uwe D. Hanebeck.  On entropy approximation for Gaussian mixture random vectors.  In _2008 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems_ , pages 181–188, Seoul, August 2008. IEEE.  ISBN 978-1-4244-2143-5.  doi: 10.1109/MFI.2008.4648062.
  * Zhang and Luo [2017] Chenchen Zhang and Yuan Luo.  Approximating the differential entropy of gaussian mixtures.  In _GLOBECOM 2017-2017 IEEE Global Communications Conference_ , pages 1–6. IEEE, 2017.
  * Michalowicz et al. [2008] Joseph V. Michalowicz, Jonathan M. Nichols, and Frank Bucholtz.  Calculation of Differential Entropy for a Mixed Gaussian Distribution.  _Entropy_ , 10(3):200–206, September 2008.  ISSN 1099-4300.  doi: 10.3390/entropy-e10030200.
  * Ariel and Louzoun [2020] Gil Ariel and Yoram Louzoun.  Estimating differential entropy using recursive copula splitting.  _Entropy. An International and Interdisciplinary Journal of Entropy and Information Studies_ , 22(2):236, 2020.
  * Brewer [2017] Brendon J Brewer.  Computing entropies with nested sampling.  _Entropy. An International and Interdisciplinary Journal of Entropy and Information Studies_ , 19(8):422, 2017.
  * Ajgl and Šimandl [2011] Jiří Ajgl and Miroslav Šimandl.  Differential entropy estimation by particles.  _IFAC Proceedings Volumes_ , 44(1):11991–11996, January 2011.  ISSN 1474-6670.  doi: 10.3182/20110828-6-IT-1002.01404.
  * Beirlant et al. [1997] Jan Beirlant, Edward J Dudewicz, László Györfi, Edward C Van der Meulen, et al.  Nonparametric entropy estimation: An overview.  _International Journal of Mathematical and Statistical Sciences_ , 6(1):17–39, 1997.
  * Petukhova and Bunt [2009] Volha Petukhova and Harry Bunt.  Who’s next? Speaker-selection mechanisms in multiparty dialogue.  In _Workshop on the Semantics and Pragmatics of Dialogue_ , 2009.
  * De Kok and Heylen [2009] Iwan De Kok and Dirk Heylen.  Multimodal end-of-turn prediction in multi-party meetings.  In _Proceedings of the 2009 International Conference on Multimodal Interfaces_ , pages 91–98, 2009.
  * Ishii et al. [2013] Ryo Ishii, Kazuhiro Otsuka, Shiro Kumano, Masafumi Matsuda, and Junji Yamato.  Predicting next speaker and timing from gaze transition patterns in multi-party meetings.  In _Proceedings of the 15th ACM on International Conference on Multimodal Interaction_ , pages 79–86, 2013.
  * Rochet-Capellan and Fuchs [2014] Amélie Rochet-Capellan and Susanne Fuchs.  Take a breath and take the turn: How breathing meets turns in spontaneous dialogue.  _Philosophical Transactions of the Royal Society B: Biological Sciences_ , 369(1658), December 2014.  ISSN 0962-8436.  doi: 10.1098/rstb.2013.0399.
  * Ishii et al. [2017] Ryo Ishii, Shiro Kumano, and Kazuhiro Otsuka.  Prediction of Next-Utterance Timing using Head Movement in Multi-Party Meetings.  In _Proceedings of the 5th International Conference on Human Agent Interaction_ , HAI ’17, pages 181–187, New York, NY, USA, October 2017. Association for Computing Machinery.  ISBN 978-1-4503-5113-3.  doi: 10.1145/3125739.3125765.
  * Malik et al. [2020] Usman Malik, Julien Saunier, Kotaro Funakoshi, and Alexandre Pauchet.  Who speaks next? Turn change and next speaker prediction in multimodal multiparty interaction.  In _2020 IEEE 32nd International Conference on Tools with Artificial Intelligence (ICTAI)_ , pages 349–354. IEEE, 2020.
  * Hung et al. [2008] Hayley Hung, Dinesh Babu Jayagopi, Sileye Ba, Jean-Marc Odobez, and Daniel Gatica-Perez.  Investigating automatic dominance estimation in groups from visual attention and speaking activity.  page 233. ACM Press, 2008.  ISBN 978-1-60558-198-9.  doi: 10.1145/1452392.1452441.
  * Pavllo et al. [2018] Dario Pavllo, David Grangier, and Michael Auli.  QuaterNet: A Quaternion-based Recurrent Model for Human Motion.  _arXiv:1805.06485 [cs]_ , July 2018.
  * [73] David Baehrens, Timon Schroeter, Stefan Harmeling, Motoaki Kawanabe, and Katja Hansen.  How to Explain Individual Classification Decisions.  page 29.
  * Sundararajan et al. [2017] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In _International conference on machine learning_ , pages 3319–3328. PMLR, 2017.
  * Smilkov et al. [2017] Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise.  _arXiv preprint arXiv:1706.03825_ , 2017.
  * Shrikumar et al. [2017] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje.  Learning important features through propagating activation differences.  In _International conference on machine learning_ , pages 3145–3153. PMLR, 2017.
  * Lundberg and Lee [2017] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  _Advances in neural information processing systems_ , 30, 2017.
  * Bargal et al. [2018] Sarah Adel Bargal, Andrea Zunino, Donghyun Kim, Jianming Zhang, Vittorio Murino, and Stan Sclaroff.  Excitation Backprop for RNNs.  _arXiv:1711.06778 [cs]_ , March 2018.
  * Zeiler and Fergus [2014] Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  In _European Conference on Computer Vision_ , pages 818–833. Springer, 2014.
  * Mujkanovic et al. [2020] Felix Mujkanovic, Vanja Doskoč, Martin Schirneck, Patrick Schäfer, and Tobias Friedrich.  timeXplain – A Framework for Explaining the Predictions of Time Series Classifiers.  _arXiv:2007.07606 [cs, stat]_ , July 2020.
  * Suresh et al. [2017] Harini Suresh, Nathan Hunt, Alistair Johnson, Leo Anthony Celi, Peter Szolovits, and Marzyeh Ghassemi.  Clinical intervention prediction and understanding using deep networks.  _arXiv preprint arXiv:1705.08498_ , 2017.
  * Molnar [2020] Christoph Molnar.  _Interpretable machine learning_.  Lulu. com, 2020.
  * Ismail et al. [2020] Aya Abdelsalam Ismail, Mohamed Gunady, Héctor Corrada Bravo, and Soheil Feizi.  Benchmarking Deep Learning Interpretability in Time Series Predictions.  _arXiv:2010.13924 [cs, stat]_ , October 2020.
  * Choi et al. [2016] Edward Choi, Mohammad Taha Bahadori, Joshua A. Kulas, Andy Schuetz, Walter F. Stewart, and Jimeng Sun.  RETAIN: An Interpretable Predictive Model for Healthcare using Reverse Time Attention Mechanism.  _arXiv:1608.05745 [cs]_ , August 2016.
  * Alaa and van der Schaar [2019] Ahmed M Alaa and Mihaela van der Schaar.  Attentive state-space modeling of disease progression.  _Advances in neural information processing systems_ , 32, 2019.
  * Panja et al. [2022] Madhurima Panja, Uttam Kumar, and Tanujit Chakraborty.  An Interpretable Probabilistic Autoregressive Neural Network Model for Time Series Forecasting.  _arXiv:2204.09640 [cs, stat]_ , April 2022.
  * Li et al. [2021] Longyuan Li, Junchi Yan, Xiaokang Yang, and Yaohui Jin.  Learning Interpretable Deep State Space Model for Probabilistic Time Series Forecasting.  _arXiv:2102.00397 [cs, stat]_ , January 2021.


,

#  Ab-initio semi-classical electronic transport in ZnSe: The role of
inelastic scattering mechanisms

Anup Kumar Mandia1, Renuka Patnaik2, Bhaskaran Muralidharan1 Seung Cheol Lee2
<EMAIL_ADDRESS>Satadeep Bhattacharjee2
<EMAIL_ADDRESS>1Indian Institute of Technology,
Mumbai-400076, India
2Indo Korea Science and Technology Center, Bangalore-560065, India

###### Abstract

We present a detailed ab-initio study of semi-classical transport in n-ZnSe
using Rode’s iterative method. Inclusion of ionized impurity, piezoelectric,
acoustic deformation and polar optical phonon scattering and their relative
importance at low and room temperature for various n-ZnSe samples are
discussed in depth. We have clearly noted that inelastic polar optical phonon
scattering is the most dominant scattering mechanism over most of the
temperature region. Our results are in good agreement with the experimental
data for the mobility and conductivity obtained at different doping
concentrations over a wider range of temperatures. Also we compare these
results with the ones obtained with relaxation time approximation (RTA) which
clearly demonstrate the superiority of the iterative method over RTA.

Rode iterative method, Relaxation Time Approximation, first principles
calculations

##  I Introduction

Semi-classical transport models ref5 ; ref6 ; ref7 ; refa ; refb ; refc ; refd
; refe ; reff ; refg ; refh  are still actively employed for transport
calculations across various materials in the bulk form. Typically, they are
employed to calculate the mobility, conductivity and eventually also the
thermoelectric coefficients ref5 ; ref6 ; ref7 ; refa ; refb ; refc ; refd ;
refe ; reff ; refg ; refh . Most of the prevalent semi-classical simulation
platforms employ the Boltzmann transport equation (BTE) within the relaxation
time approximation (RTA) and well known packages have been employed. The RTA
treats the relaxation time as constant or as a power law in energy and is
applicable only if the relaxation time is independent of the distribution
function, which is satisfied when a scattering mechanism is either elastic or
isotropic ref4 . For instance, polar optical phonon (POP) scattering is
neither elastic or isotropic and hence the RTA is inappropriate when one deals
with POP scattering. So one has to go beyond the RTA to capture POP scattering
mechanism properly. This brings a serious issue when it comes to materials
which possess strong optical phonon scattering mechanisms. Methods such as the
Rode iterative method ref5 ; ref6 ; ref7  or variational refc ; refd
techniques are often employed to properly account for POP scattering. In this
paper, we use the Rode iterative method for mobility and conductivity
calculation. The objective of this work is to use ZnSe as an example to
elucidate the limitations of the RTA while presenting a calculation that is
based on the Rode iterative scheme solution of the BTE.
ZnSe is a wide band gap semiconductor, which is employed for making blue light
emitting diode and lasers. The electron mobility in n-ZnSe crystals was first
analyzed in Ref. refi , where the electron mobility was measured at several
temperatures and it was shown that POP scattering is the most dominant
scattering mechanism near the room temperature region. There have been a few
more works ref6 ; ref7 ; reff ; refj ; refk ; refl ; refm  which have carried
out theoretical investigations on n-ZnSe electron mobility. However, most of
these calculations were based on the RTA method. Earlier works ref6 ; ref7 ;
reff  had also done theoretical investigations beyond RTA for n-ZnSe electron
mobility using an iterative and variational technique respectively and have
shown good agreement between theory and experiment. References ref6 ; ref7
had investigated n-ZnSe electron mobility various scattering mechanisms over a
wide temperature range with iterative technique, while Ref. reff  had
investigated the dependence on donor concentrations with variable compensation
ratios with variational technique. All of these works used parabolic or Kane
model for the mobility calculation and all the constants used were taken from
the experiments.

In the present work, the electron mobility and conductivity is investigated
for n-type ZnSe at different doping concentrations and temperatures with ab-
initio inputs. For our theoretical analysis, all the parameters of band
structure, density of states, wave function admixture, dielectric constants,
piezoelectric constants, polar optical phonon (POP) frequencies, and acoustic
deformation potentials are calculated by using the density functional theory
(DFT). In ZnSe, since POP scattering is the dominant mechanism and is
inelastic, it modulates the electron energy, and hence a universal time
constant cannot be defined, making the RTA inappropriate. We thus use the Rode
iterative method with ab-initio inputs as previously discussed in Ref. ref1 .
In the Rode method, the perturbation in the distribution is obtained at low
electric field, keeping only the linear term. The perturbation in the
distribution function is used to calculate the mobility. Since all the
required inputs are calculated by using DFT, only the crystal structure is
given as input, it does not rely on experimental data. So this is the first
paper who calculated n-ZnSe electron mobility with abinito band structure and
inputs by using Rode method and a good agreement is obtained for different
concentrations over a wide temperature range. The same approach can be applied
to other new materials for which many constants are not determined
experimentally.
In this paper, we first compare the experimental data and RTA calculation, and
establish that the RTA results underestimate the mobility by more than a fifty
percent at room temperature. In the process of presenting the results here, we
also give a detailed insight about the various contributions of different
scattering mechanism in ZnSe at different temperatures and doping
concentrations. Since all the required inputs are calculated by using density
functional theory, our code can be further extended and used for mobility
calculations of new materials for which many constants are not known
experimentally.
This paper is organized as follows. In the following section, we first
describe in detail the methodology by first elaborating on the Rode scheme for
solving the Boltzmann transport equation. We specifically hint at how this
method scores over the RTA in capturing inelastic scattering processes. We
also detail the computational scheme of calculating various scattering rates
starting from ab-initio simulations. In Sec. III, we describe the results in
detail, where we first establish using the calculated scattering rates that
the inelastic POP scattering mechanism indeed dominates over a wide
temperature range in the case of n-type ZnSe. We then describe our theoretical
fits in comparison with relevant experimental works at different carrier
concentrations. In Sec. IV, we conclude with an outlook from this work.

##  II Methodology

###  II.1 The Boltzmann Transport Equation

The Boltzmann transport equation (BTE) describes the time-evolution of the
state of the charge carriers in both real and momentum co-ordinates via a
distribution function, f(r,k,t), and is given by ref2 ; ref3 ; ref4

$\frac{\partial\textit{f}}{\partial
t}+\textbf{v}\cdot\nabla_{r}f+\frac{\textbf{F}}{\hbar}\cdot\nabla_{\bf
k}f=\frac{\partial\textit{f}}{\partial
t}\Bigr{|}_{\begin{subarray}{c}coll\end{subarray}},$ (1)

where v is the carrier velocity, f describes the probability distribution
function of carriers in real and momentum space as a function of time, F is
the applied external force, $\frac{\partial\textit{f}}{\partial
t}\Bigr{|}_{\begin{subarray}{c}coll\end{subarray}}$ represents the change in
the distribution function with time due to collisions. The first term in (1)
represents the rate of change of the carrier distribution f with time. The
second term represents the diffusion due to a gradient in the carrier density
and the third term represents the change in f due to all external forces. In
the presence of an electric field E, the external force F is given by
$e\textbf{E}$, and the BTE then becomes

$\frac{\partial f}{\partial
t}+\textbf{v}\cdot\nabla_{r}f+\frac{e\textbf{E}}{\hbar}\cdot\nabla_{\bf
k}f=\frac{\partial f}{\partial
t}\Bigr{|}_{\begin{subarray}{c}coll\end{subarray}},$ (2)

Under steady state and spatially homogeneous conditions, the above equation
can be rewritten as:

$\frac{e\textbf{E}}{\hbar}\cdot\nabla_{k}f=-\int[s(k,k^{\prime})f(1-f^{\prime})-s(k^{\prime},k)\;f^{\prime}(1-f)]dk^{\prime},$
(3)

where $e$ is the electronic charge, $s(k,k^{\prime})$ represents the
transition rate from a state $k$ to a state $k^{\prime}$. At low electric
fields, the solution to the BTE is given by the distribution function ref5 ;
ref6 ; ref7

$f(k)=f_{0}[\epsilon(k)]+g(k)cos\theta,$ (4)

where $f_{0}[\epsilon(k)]$ is the equilibrium distribution function, and
$\cos\theta$ is the angle between applied electric field and $k$. Here we have
neglected the higher order terms, since we are calculating mobility under low
electric field conditions. Now we have to calculate the perturbation in the
distribution function $g(k)$ for calculating the low-field transport
properties. The perturbation $g(k)$ is given by ref5 ; ref6 ; ref7

$g_{k,i+1}=\frac{S_{i}(g_{k},i)-v(k)(\frac{\partial f}{\partial
z})-\frac{eE}{\hbar}(\frac{\partial f}{\partial
k})}{S_{o}(k)+\frac{1}{\tau_{el}(k)}},$ (5)

Where we have considered the electric field along the z-direction. The
transition rate in Eq.3 can be split into two contributions: One from from
inelastic (in) and other being the elastic (el) component respectively, i.e,
$s(k,k^{\prime})=s(k,k^{\prime})_{in}+s(k,k^{\prime})_{el}$. The expression
for $S_{i}$ and $S_{o}$ which appear in Eq.5 are given by ref7

$\
S_{o}(k)=\int[s_{in}(k,k^{\prime})(1-f^{\prime})+s_{in}(k^{\prime},k)f^{\prime}]dk^{\prime}$
(6) $\ S_{i}(g_{k},i)=\int
Xg_{k^{\prime},i}[s_{in}(k^{\prime},k)(1-f)+s_{in}(k,k^{\prime})f]dk^{\prime}$
(7)

where X is the cosine of the angle between the initial and the final wave
vectors. The elastic part of the scattering rate is given by,
$\frac{1}{\tau_{el}}=\int(1-X)S_{el}(k,k^{\prime})dk^{\prime}$.
As a initial guess take $g_{k,0}=0$, it will give $S_{i}(g_{k},0)=0$. Then
$g_{k,1}$ is given by

$g_{k,1}=\frac{-v(k)(\frac{\partial f}{\partial
z})-\frac{eE}{\hbar}(\frac{\partial f}{\partial
k})}{S_{o}(k)+\frac{1}{\tau_{el}(k)}},$ (8)

The value of $g_{k,1}$ is used to calculate to $S_{i}(g_{k},1)$, then
$S_{i}(g_{k},1)$ is used to calculate $g_{k,2}$.

$g_{k,2}=\frac{S_{i}(g_{k},1)-v(k)(\frac{\partial f}{\partial
z})-\frac{eE}{\hbar}(\frac{\partial f}{\partial
k})}{S_{o}(k)+\frac{1}{\tau_{el}(k)}}.$ (9)

These steps are repeated until $g_{k,i+1}$ converges. This convergence is
exponential ref7 , so normally it requires a few iteration for convergence.
Typically, five iterations are required for the perturbation to converge for
polar optical phonon scattering in ZnSe. Now, the calculated perturbation is
used to calculate mobility and Seebeck coefficient. The effect of inelastic
POP scattering is included through the terms $S_{i}(g)$ and $S_{o}$ and the
effect of elastic scattering is included through the relaxation time
$\frac{1}{\tau_{el}(k)}$ term. The second term in the numerator of (5)
represents the thermal driving force and the third term in the numerator
represents the electrical driving force. For mobility calculations, the
thermal driving force is set to zero and only the electrical driving force is
considered. $\frac{1}{\tau_{el}(k)}$ is the sum of the momentum relaxation
rates of all elastic scattering process, which is given by

$\
\frac{1}{\tau_{el}(k)}=\frac{1}{\tau_{ii}(k)}+\frac{1}{\tau_{pz}(k)}+\frac{1}{\tau_{ac}(k)},$
(10)

where the subscripts el,ii , pz and ac are used for elastic, ionized impurity,
piezoelectric and acoustic deformation potential scattering processes
respectively. The carrier mobility $\mu$ is then given by ref1 ; ref5 ; ref6 ;
ref7

$\mu=\frac{1}{3E}\frac{\int
v(\epsilon)D_{s}(\epsilon)g(\epsilon)d\epsilon}{\int
D_{s}(\epsilon)f(\epsilon)d\epsilon},$ (11)

where $D_{S}(\epsilon)$ represents density of states. The carrier velocity is
calculated directly from the ab-initio band structure by using

$\ v(k)=\frac{1}{\hbar}\frac{\partial{\epsilon}}{\partial{k}}.$ (12)

From these, we can evaluate the electrical conductivity given as

$\ \sigma=ne\mu_{e},$ (13)

where $n$ is the electron carrier concentration, and $\mu_{e}$ is the electron
mobility.

###  II.2 Computational method

The electronic structure calculations ref26  are performed using first-
principles methods within the frame-work of DFT with Perdew-Burke Ernzerhof
exchange correlation energy functionalpbe  based on a generalized gradient
approximation. We used a projector augmented wave method as implemented in
Vienna ab-initio simulation package (VASP)ref8 ; ref9 ; ref10 . The Kohn-Sham
wave functions of the valence electrons were expanded in plane wave basis with
energy cut-off of 500 eV. Ionic relaxation was performed using conjugate-
gradient method, until forces were reduced to within 0.01 eV/Angstrom. The
Brillouin zone sampling was carried out using Monkhorst Pack grid of 11x11x11
k-points. The band structure is computed along the high-symmetry k-points in
the irreducible Brillouin zone, with 100 k-points between each pair of high-
symmetry points. Computed band structure with the self-consistent density of
states (DOS) is shown in the Fig.1. Since as an input for the transport
calculation within Rode’s method, only band structure for one valley is
needed, we have performed non-self consistent calculations of the band
energies in a special k-point mesh around the $\Gamma$ point with 8531
k-points. Using such a dense mesh we have obtained very accurate group
velocity and effective mass.

###  II.3 Scattering Mechanisms

In this work, for calculating the mobility of ZnSe, four different types of
scattering mechanisms are included (i) Ionized impurity scattering, (ii)
Acoustic deformation potential scattering, (iii) Piezoelectric scattering and
(iv) Polar optical phonon scattering. For ionized impurity scattering, the
formulation provided by the Brooks-Herring approach is used. The ionized
impurity is a significant scattering mechanism at low temperatures and high
doping concentrations. The momentum relaxation rate for this mechanism with
its abinito counterpart is given by following expression ref1 ; ref7

$\frac{1}{\tau_{ii}(k)}=\frac{e^{4}N}{8\pi\epsilon_{0}^{2}\hbar^{2}k^{2}v(k)}[D(k)ln(1+\frac{4k^{2}}{\beta^{2}})-B(k)],$
(14)

where $\epsilon_{0}$ is the dielectric constant, $\hbar$ is the reduced Planck
constant, and $\beta$ is the inverse screening length given by

$\beta^{2}=\frac{e^{2}}{\epsilon_{0}k_{B}T}\int
D_{s}(\epsilon)f(1-f)d\epsilon,$ (15)

where N is the concentration of ionized impurity and it is given by

$\ N=N_{A}+N_{D}$ (16)

where $N_{A}$ and $N_{D}$ are the acceptor and donor concentrations
respectively. The expressions for D(k) and B(k) are taken from equations (91)
and (92) of Ref. ref7 .

$\ D(k)=1+\frac{2\beta^{2}c^{2}}{k^{2}}+\frac{3\beta^{4}c^{4}}{4k^{4}}$ (17)

$\
B(k)=\frac{4k^{2}/\beta^{2}}{1+4k^{2}/\beta^{2}}+8\frac{\beta^{2}+2k^{2}}{\beta^{2}+4k^{2}}c^{2}+\frac{3\beta^{4}+6\beta^{2}k^{2}-8k^{4}}{(\beta^{2}+4k^{2})k^{2}}c^{4}$
(18)

Acoustic deformation occurs due to the coupling of electrons with non-polar
acoustic phonons. The momentum relaxation rate for acoustic deformation
potential scattering is given by ref1 ; ref7

$\frac{1}{\tau_{ac}(k)}=\frac{e^{2}k_{B}TE_{D}^{2}k^{2}}{3\pi\hbar^{2}c_{el}v(k)}[3-8c^{2}(k)+6c^{4}(k)],$
(19)

where c(k) is the contribution of the p-type orbital to the wave function of
the conduction band, $c_{el}$ is spherically averaged elastic constant and
$E_{D}$ is acoustic deformation potential and is given by conduction band
shift (in eV) per unit strain due to acoustic waves. For ab-initio
calculations, the wave function admixture $c(k)$ is obtained through
projecting the Kohn-Sham wavefunctions onto the spherical harmonics which are
non-zero only within spheres centering the ions and is implemented in VASP
package.
Chemical bonds in compound semiconductors such as ZnSe are partly ionic in
nature. Zn atom has a slight positive and Se atom has a slight negative
charge. The magnitude of this charge is determined by the degree of the ionic
nature of the bond, and it is a fraction of electronic charge ref24 . The
vibrations of atoms cause changes in the lattice constant. This perturbs the
dipole moment between the atoms that eventually scatter the electrons. The of
polar scattering due to the long-wavelength acoustic phonons is called
piezoelectric scattering and the polar scattering due to optical phonons is
called polar optical phonon (POP) scattering. Piezoelectric scattering is
important at low temperatures and at low doping densities in polar materials.
Since ZnSe is polar, it necessary to include piezoelectric scattering to get a
better fit for the mobility output at low temperature. The momentum relaxation
rate for piezoelectric scattering with ab-inito parameters as input is given
by ref1 ; ref7

$\frac{1}{\tau_{pz}(k)}=\frac{e^{2}k_{B}TP^{2}}{6\pi\epsilon_{0}\hbar^{2}v(k)}[3-6c^{2}(k)+4c^{4}(k)]$
(20)

where $P$ is dimensionless piezoelectric coefficient and it is given by ref7

$\
P^{2}=h_{14}^{2}\epsilon_{0}\frac{[(\frac{12}{c_{l}})+(\frac{16}{c_{t}})]}{35}$
(21)

where $h_{14}$ is one element of piezoelectric stress tensor and $c_{l}$
$c_{t}$ are the spherically averaged elastic constant for longitudinal and
transverse modes respectively and are given by ref7

$\ c_{l}=(3c_{11}+2c_{12}+4c_{44})/5$ (22)

$\ c_{t}=(c_{11}-c_{12}+3c_{44})/5$ (23)

where $c_{11}$, $c_{12}$ and $c_{44}$ are three independent elastic constants.
The POP scattering is the most dominant scattering mechanism near room
temperature and in the higher temperature regime. Since POP scattering is
inelastic and anisotropic, the Rode iterative scheme is used to directly
evaluate the scattering rates so that the momentum relaxation rate is given by
ref5 ; ref7

$S_{o}=(N_{po}+1-\textit{f}^{-})\lambda^{-}_{0}+(N_{po}+\textit{f}^{+})\lambda^{+}_{0}$
(24)
$\lambda^{+}_{o}=\beta^{+}[(A^{+})^{2}ln\mid\frac{k^{+}+k}{k^{+}-k}\mid-A^{+}cc^{+}-aa^{+}cc^{+}]$
(25) $\beta^{+}=\frac{e^{2}\omega_{po}k^{+}}{4\pi\hbar
kv(k^{+})}(\frac{1}{\epsilon_{\infty}}-\frac{1}{\epsilon_{0}})$ (26)
$A^{+}=aa^{+}+\frac{(k^{+})^{2}+k^{2}}{2k^{+}k}cc^{+},$ (27)

where the subscript plus denotes the scattering out by absorption so it is to
be evaluated at an energy $\epsilon+\hbar\omega_{po}$ and the subscript minus
denotes scattering out by emission so that it is to be evaluated at energy
$\epsilon-\hbar\omega_{po}$. If energy is less than $\hbar\omega_{po}$, then
the emission of phonons is not possible and hence $\lambda^{-}_{0}$ is to be
considered to be zero, $N_{po}$ is the number of phonons and is given by

$N_{po}=\frac{1}{exp(\hbar\omega_{po}/k_{B}T)-1}.$ (28)

The in scattering operator for POP is given by ref5 ; ref7

$S_{i}=(N_{po}+\textit{f})\lambda^{-}_{i}g^{-}+(N_{po}+1-\textit{f})\lambda^{+}_{i}g^{+}$
(29)
$\lambda^{+}_{i}(k)=\beta^{+}[\frac{(k^{+})^{2}+k^{2}}{2k^{+}k}(A^{+})^{2}ln\mid\frac{k^{+}+k}{k^{+}-k}\mid-(A^{+})^{2}-\frac{c^{2}(c^{+})^{2}}{3}]$
(30)

###  II.4 Ab-initio Inputs

The computed band structure and density of states for ZnSe is shown in Fig. 1.
As already mentioned, we have calculated the band structure using the density
functional theory using a three dimensional $k$ mesh around the conduction
band minimum (CBM). The conduction band is expressed as the average energy of
the electrons as a function of k from the CBM. In order to evaluate the group
velocities, we have first calculated the derivatives of the conduction band
energy with respect to $k$, performed an analytical fitting of the conduction
band with a six degree polynomial and divided the conduction band into four
segments ref1 ; ref11  to obtain a smooth curve for both mobility and
conductivity. The Fermi level is obtained by calculating the carrier
concentration using equation (31) and matching it to the given concentration:

$n=\frac{1}{V_{0}}\int_{\epsilon_{c}}^{\infty}{D_{S}(\epsilon)f(\epsilon)d\epsilon},$
(31)

where $D_{S}(\epsilon)$ represents density of states at energy $\epsilon$,
where $\epsilon_{c}$ represents the bottom of conduction band and $V_{0}$
represents the volume of the cell. We have calculated the low and high
frequency dielectric constants ref12 ; ref13 , polar optical phonon
frequencies $\omega_{po}$ ref14 , the elastic constant ref24 ; ref25  and the
piezoelectric constant by using the density functional perturbation theory as
implemented in VASP. Calculated abintio value of low and high frequency
dielectric constant is 7.45 and 3.44 respectively. We have obtained abinitio
polar optical phonon frequency of 5.88 THz. The computed phonon band structure
with phonon density of states for ZnSe is shown in figure 2. The acoustic
deformation potential is obtained by mimicking a uniform lattice deformation
due to acoustic phonons which results a shift of the conduction band minimum.
The deformation potential is given by

$E_{D}=-V(\frac{\partial E_{CBM}}{\partial V})|_{V=V_{0}}$ (32)

Where ECBM is the conduction band minimum and V0 is the equilibrium volume. We
have obtained abinitio value of 0.574 nm of lattice constant, 1.17 eV of band
gap, 12 eV of acoustic deformation potential, 0.0392 of Piezoelectric
coefficient and -0.08162 $C/m^{2}$ of Piezoelectric constant($e_{14}$).
Calculated abinitio value of elastic constant $c_{11}$, $c_{12}$ and $c_{44}$
are $7.99\times 10^{10}N/m^{2}$, $4.54\times 10^{10}N/m^{2}$ and $3.71\times
10^{10}N/m^{2}$ respectively.

Electronic correlation effects often play a major role in the carrier
transport for $\mathrm{3d-4p}$ systems. An efficient approach for the
inclusion of such correlation effect is multiband mean field Hubbard model
(DFT+U). Karazhanov et al. ref27  used this approach to study the effect of
correlation in ZnX (X = O, S, Se, Te) and concluded that the effects of
correlation in ZnO are more important than in ZnS, ZnSe and ZnTe. The spin-
orbit interaction (SOI) is another effect that we have not included in our
calculation. Usually for the semiconductors with zinc blend structures such as
ZnSe, the effect of SOI is to lift the six-fold degeneracy of the valence band
edge at $\Gamma$-point. However, since we are only interested in electronic
transport (not the hole) and our formalism does not depend on the exact
position or the nature of the edge of the valence band, we have neglected the
effect of the SOI.

###  II.5 Simulation Flowchart

The flowchart for our mobility and conductivity calculations using the Rode
method with the ab-initio inputs is shown in Fig. 3. First, we have to
calculate all the required inputs using the first principles methods detailed
in the previous section. Then we have to perform the analytical fitting of the
band structure, to obtain smooth curves in order to calculate the group
velocities of the carriers with $k$. Then the smoothed band structure is used
to calculate the Fermi level using Eq. (31). The calculations are sensitive
with respect to the calculated Fermi level and hence must be done carefully
with good precision. Next we have to calculate various scattering rates for
different scattering mechanisms by using (14), (19), (20), and (24). These
scattering rates are then used to calculate the perturbation in the
distribution function using (5). Now the obtained $g(k)$ is used to calculate
the desired transport coefficients. If $g(k)$ is obtained right after the
first order iterate, it will replicate the RTA results of the transport
coefficients.

##  III Results

###  III.1 Scattering Rate vs Energy

For a comparison between the different scattering mechanisms in n-ZnSe, the
scattering rate as a function of energy is plotted in Fig. 4 and Fig. Ab-
initio semi-classical electronic transport in ZnSe: The role of inelastic
scattering mechanisms, for different doping concentrations $N_{D}=1\times
10^{10}cm^{-3}$ and $N_{D}=1\times 10^{15}cm^{-3}$ at temperatures 30 K and
300 K. At low doping and lower temperatures, piezoelectric scattering is
considered the most dominant scattering mechanism for low energy carriers. At
a temperature of 30 K, the average energy of carriers is
$\frac{3}{2}k_{B}T=0.0038eV$, and hence most of the carriers are in the low
energy region. At low temperatures, around the average energy of carriers,
both piezoelectric and acoustic deformation potential scattering mechanisms
are approximately equal. So at lower doping and at lower temperatures,
piezoelectric scattering and acoustic deformation potential scattering are the
most dominant scattering mechanisms for ZnSe.
At higher doping ($1\times 10^{15}$) for low temperatures, ionized impurity
scattering is the most dominant for low energy carriers. In Fig. 4 and Fig.
Ab-initio semi-classical electronic transport in ZnSe: The role of inelastic
scattering mechanisms there is a sudden change in POP scattering rate after a
particular energy, which is due to the fact that if an electron energy is
smaller than the POP energy $\hbar\omega_{po}=0.024eV$, then it can scatter
only by the absorption of an optical phonon, but if the energy of an electron
is larger than the optical phonon energy, then it can scatter by both emission
and absorption of optical phonons. At 300 K, POP scattering is the most
dominant scattering mechanism for all electron energies for both higher and
lower doping making it the most important scattering mechanism at temperatures
ranging from room temperature to higher temperature. Since all working devices
operate in this temperature region, it is necessary to properly include POP
scattering for a good theoretical mobility calculation. Since RTA is not valid
for inelastic scattering mechanisms, other abinito RTA based codes such as the
BoltzTrap ref20  code may not produce the desired results in the case of
n-ZnSe.

###  III.2 Electron Mobility

By using the Rode iterative method based on the flowchart detailed above, the
electron mobility is calculated for three different experimentally
characterized n-ZnSe samples of different doping concentrations over a wide
temperature range. The donor and acceptor concentrations for different
experimentally characterized n-ZnSe samples taken from different references is
shown in Tab. 1. Figure 5 shows a comparison between the experimentally
measured and theoretically calculated mobility using our Rode scheme as well
as the RTA method. Overall there is good agreement between the theoretical and
experimental curves when considering the Rode scheme results and score much
better than RTA results. The RTA underestimates the mobility for all samples,
which is due to POP scattering, since POP scattering is inelastic and
anisotropic and hence changes the distribution function such that a constant
or simple power law form for the relaxation time cannot be defined. The
calculated mobility for sample (a), (b) and (c) using the Rode method have an
average relative error of $18.44\%$ and $21.63\%$ and $4.08\%$ respectively.
While the mobility calculated using the RTA have average errors of $113.46\%$,
$55.69\%$ and $42.77\%$ for sample (a), (b) and (c) respectively.
For sample (a), the Rode calculation shows an error of $-15.03\%$ and
$-7.27\%$ at 200 K and 300 K respectively, while the RTA shows an error of
$-133.45\%$ and $-117.88\%$ at 200 K and 300 K respectively. At lower
temperatures, even RTA shows a smaller error since at lower temperatures, POP
scattering becomes insignificant, since the POP energy of
$\hbar\omega_{po}=0.024eV$ is required. Around 40 K or lower, the electrons
have smaller probability of scattering by the absorption or emission of 0.024
eV of optical phonon energy. The Rode method shows an error of $-6.45\%$ at
300 K for sample (b), while the RTA shows an error of $-54.26\%$ at 300 K with
respect to the experimental curve. For sample (c), the Rode calculation shows
an error of $-3.86\%$ and $-1.05\%$ at 200 K and 300 K respectively, while the
RTA shows an error of $-51.49\%$ and $-57.85\%$ at 200 K and 300 K
respectively. Therefore, all in all, the RTA results are quite inappropriate
for mobility calculations of polar materials.
Figure 6 shows the calculated conductivity with respect to the experimental
one for sample (b). Sample (b) contains shallow as well as deep donors. With
the increasing temperature conductivity deceases in experimental curve, but
around 230 K there is sudden change in conductivity curve slope, this is due
to ionization of deep donors. Deep donors have a doping concentration of
$N_{deep}=4.5\times 10^{15}cm^{-3}$ ref22  and have ionization energy of
$E_{deep}=130meV$ ref22 . The temperature dependent free electron
concentration for the semiconductor having deep donors is given by equation

$\frac{n(n-n_{0})}{N_{deep}+n_{0}-n}=\frac{N_{C}}{g}exp(\frac{-E_{deep}}{k_{B}T})$
(33)

where $g$ is the degeneracy of deeper level, here $g=2$ is taken for
calculation and $n_{0}$ is the concentration of the electrons activated from
shallow donors in the region of their exhaustion. The conductivity calculated
by the Rode iterative method is again much better than that calculated by
using the RTA method. The calculated conductivity shows a very good
qualitative and quantitative agreement with respect to the experimental data.
Figure 7 shows the overall mobility and the mobility by considering only one
type of scattering mechanism at a time for sample (c). The lower mobility is
most dominant in deciding the overall mobility. From this figure it is clear
that POP is the most dominant scattering mechanism, from 70 K to higher
temperatures. Figures 8 and 9 show the mobility for different doping
concentrations at 77 K and 300 K by assuming a compensation ratio of unity.
With increasing doping concentration, the curves show a decrease in mobility
due to an increase in the number of ionized centers.

##  IV Conclusion

We presented an ab-inito semi-classical transport calculation for the mobility
and conductivity of n-ZnSe by using Rode iterative method in order to
conclusively illustrate the role of inelastic scattering processes. We
recognized that inelastic polar optical phonon scattering is the most dominant
scattering mechanism over most of the temperature region. A good agreement
with various experiments was observed for different doping concentrations over
a wider range of temperatures. In comparing our results against that obtained
using the relaxation time approximation method, we clearly noted the
discrepancy in explaining experimental results, thereby pointing out the need
to advance semi-classical transport calculations beyond the relaxation time
approximation. Further work will be extended p-type semiconductor and
multivalley transport.

## References

  * (1) D. L. Rode, Phys. Rev. B 2, 1012(1970).
  * (2) D. L. Rode, Phys. Rev. B 2, 4036(1970).
  * (3) D. L. Rode, Semiconductors and Semimetals (Academic Press, New York, 1975), Chapter 1.
  * (4) A. T. Ramu, L. E. Cassels, N. H. Hackman, H. Lu, J. M. O. Zide, and J. E. Bowers, Journal of Applied physics 107, 083707 (2010).
  * (5) A. T. Ramu, L. E. Cassels, N. H. Hackman, H. Lu, J. M. O. Zide, and J. E. Bowers 109, 033704 (2011).
  * (6) D. J. Howarth and E. H. Sondheimer, Proc. R. Soc. Lond. A 219, 53(1953).
  * (7) H. Ehrenreich, Phys. Rev. 120, 1951 (1960).
  * (8) H. E. Ruda and B. Lai, Journal of Applied physics 68, 1714 (1990).
  * (9) H. E. Ruda, Journal of Applied Physics 59, 1220(1986).
  * (10) A. Sztein, J. Haberstroh, J. E. Bowers, S. P. DenBaars, and S. Nakamura, Journal of Applied Physics 113, 183707 (2013).
  * (11) A. Popescu, A. Datta, G. S. Nolas, and L. M. Woods, Journal OF Applied Physics 109, 103709 (2011).
  * (12) M. Lundstrom, Fundamentals of Carrier Transport, 2nd ed. (Cambridge University Press, Cambridge, UK, 2009).
  * (13) M. Aven and B. Segall, Phys. Rev. 130, 81 (1963).
  * (14) D. D. Nedeog, Phys. stat. sol.(b) 80, 369(1977).
  * (15) 0\. V. Emelyanenko, G.N.Ivanova, T.S. Lagunova, D. D. Nedeoglo, G. M. Shmelev and A. V. Simashkevich, phys. stat. sol. (b) 96, 823(1979).
  * (16) Y. Fukuda and M. Fukai, J.Phys. Soc. Japan 23, 902(1967).
  * (17) B. R. Sethi, Phys. stat. sol. (a) 42, 791 (1977).
  * (18) A. Faghaninia, J. W. Ager III and C. S. Lo, Phys. Rev. B 91, 235123 (2015).
  * (19) Jasprit Singh, Electronic and Optoelectronic Properties of Semiconductor Structure, (Cambridge University Press, Cambridge, UK, 2003).
  * (20) D. K. Ferry, Semiconductor Transport, (Taylor & Francis, London, 2000).
  * (21) Y.Petroff, M.Balkanski, John P. Walter and Marvin L.Cohen, Solid State Communications 7, 459(1969).
  * (22) Perdew, J.P, Burke, K., Ernzerhof, M., Physical Review Letters 77, 3865 (1996)
  * (23) G. Kresse and J. Hafner, Physical Review B 49, 14251 (1994)
  * (24) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
  * (25) G.Kresse and J. Furthmullerb, Computational Materials Science 6, 15 (1996).
  * (26) B. K. Ridley, Quantum processes in semiconductors, (Clarendon Press, Oxford,1999).
  * (27) A. Faghaninia, J. W. Ager and C. S. Lo, 2015 IEEE 42nd Photovoltaic Specialist Conference (PVSC), New Orleans, LA, 2015, pp. 1-4.
  * (28) X. Gonze and Changyol Lee, Phys. Rev. B 55 10355 (1997).
  * (29) B. Liu, M. Gu, Z. Qi, X. Liu, S. Huang and C. Ni, Phys. Rev. B 76, 064307 (2007).
  * (30) P. Giannozzi and S. D. Gironcoli, Phys. Rev. B 43, 7231 (1990).
  * (31) B. H. Lee, Journal of Appl. Phys. 41, 2984(1970).
  * (32) C. G. Hodgins and J. C. Irwin, Physica Status Solidi (a) 28, 647(1975).
  * (33) A. Kartsev, D. Feya Oleg, N. Bondarenko and A. G. Kvashnin, Phys. Chem. Chem. Phys. 21, 5262(2019).
  * (34) Karazhanov et al., JOURNAL OF APPLIED PHYSICS 100, 043709 (2006)
  * (35) G. K. H Madsena and D. J. Singh, Computer Phy. Comm. 175, 67(2006).
  * (36) R.J. Nelmes, M.I. McMahon, Semiconductors and Semimetals (Academic Press, New York, 1998), Chapter 3.
  * (37) K. C. Agarwal, B. Daniel, T. Hofmann, M. Schubert, C. Klingshirn, M. Hetterich, Phys. Status Solidi B. 243, 914 (2006).
  * (38) D. T. F. Marple, Journal Of Applied Physics 35, 1879 (1964).
  * (39) C. G. Hodgins and J. C. Irwin, Phys. Stat. Sol. 28, 647 (1975).
  * (40) J. D. Zook, Phys. Rev. 136, A869 (1964).
  * (41) M. Aven, Journal of Applied Physics 42, 1204(1971).
  * (42) N. D. Nedeoglo, V. P. Sirkeli, D. D. Nedeoglo, R Laiho and E Lahderanta, Journal of Phy.: Condensed Matter 18, 8113(2006).
  * (43) A. N. Avdonin, D. D. Nedeoglo, N. D. Nedeoglo, and V. P. Sirkeli, Phys. Stat. Sol. (b) 238, 45(2003).

Figure 1: Band structure and density of states for ZnSe, Fermi level is set to
zero at valence band maximum

Figure 2: Phonon band structure with phonon density of states for ZnSe

Figure 3:  Flowchart for transport coefficient calculation from ab-initio
inputs

(a) 30 K (b) 300 K

Figure 4: Scattering Rates vs Energy for $N_{D}=1\times 10^{10}cm^{-3}$

(a) $30K$

(b) $300K$ (c) Scattering Rates vs Energy for $N_{D}=1\times 10^{15}cm^{-3}$

Table 1:  Doping of different experimentally fabricated n-ZnSe Samples Sample | $N_{D}-N_{A}(cm^{-3})$ | Donor, $N_{D}(cm^{-3})$ | Acceptor, $N_{A}(cm^{-3})$
---|---|---|---
a ref21  | $1\times 10^{15}$ | $2.9\times 10^{15}$ | $1.9\times 10^{15}$
b ref22  | $1.1\times 10^{16}$ | $6\times 10^{16}$ | $4.9\times 10^{16}$
c ref23  | $6.3\times 10^{15}$ | $7.5\times 10^{15}$ | $1.2\times 10^{15}$

Figure 5: Calculated and experimental measured mobility with temperature
variation for ZnSe at Different Doping. More detail about donor and acceptor
concentration is shown in table 1

Figure 6: Calculated and experimental measured conductivity with temperature
variation for sample b

Figure 7: Contribution of mobility from different scattering mechanisms for
Sample (c)

Figure 8: Calculated mobility for different doping concentration at 77 K

Figure 9: Calculated mobility for different doping concentration at 300 K