# A real-time TDDFT scheme for strong-field interaction in Cartesian
coordinate grid

Abhisek Ghosal111Present Address: Department of Chemical Sciences, Tata
Institute of Fundamental Research (TIFR), Mumbai-400005, Maharastra, India
Amlan K. Roy  Email<EMAIL_ADDRESS><EMAIL_ADDRESS> Department of
Chemical Sciences
Indian Institute of Science Education and Research (IISER) Kolkata,
Nadia, Mohanpur-741246, WB, India.

###### Abstract

In this communication, we present a new approach towards RT-TDDFT through
time-dependent KS equations based on an _adiabatic eigenstate subspace_ (AES)
procedure. It introduces a second-order split operator technique in energy
representation to implement the approximate TD propagator in AES. Most of the
elements in TDKS matrix are directly computed in Cartesian coordinate grid
(CCG). To demonstrate the internal consistency of our proposed scheme, we
computed the TD dipole moment and high harmonic generation spectra using an
adiabatic local density approximation. The comparison with available
theoretical results ensures the feasibility of this proposed route.

Keywords: Strong-field interaction, adiabatic approximation, real time,
Cartesian grid, TDDFT.

##  I Introduction

Throughout the past several decades, time-dependent density functional theory
(TDDFT) has been found to be the most popular method for computing linear
response and excitation energies of many-electron systems, mostly on atoms,
molecules, and molecular clusters, due to its reasonable trade-off between
accuracy and efficiency (Casida and Huix-Rotllant, 2012; Marques et al., 2012;
Adamo and Jacquemin, 2013). It can be viewed as a reformulation of time-
dependent (TD) quantum mechanics where the basic variable is TD single-
particle electron density, $\rho(\mathrm{\mathbf{r}},t)$, instead of many-
electron wave function. This is based on a one-to-one correspondence between
TD external potential and $\rho(\mathrm{\mathbf{r}},t)$, for many-body systems
evolving from a given initial state (Runge and Gross, 1984), what is called
TDDFT. Now, the most widely recognized and implemented version of TDDFT is
that of linear-response theory (LR-TDDFT) based on Casida’s equation (Casida,
1995). It has been developed for computing excitation energies, mostly for
discrete spectrum. In comparison to the conventional wave function-based
formalisms, such as CI, MR-CI, complete active space self-consistent field,
equation-of-motion coupled cluster etc., LR-TDDFT has found wide applicability
in medium to large systems, for same obvious reason. This is evident from its
applicability to medium to large systems (Casida and Huix-Rotllant, 2012;
Marques et al., 2012; Adamo and Jacquemin, 2013). Though it is found to be the
most economically efficient method to compute low-lying excited states for
finite systems, it is limited to weak-field processes (in which the density is
perturbed slightly and remains close to ground state).

As an alternative to LR-TDDFT in frequency domain, one can directly solve the
time-dependent Kohn-Sham (TDKS) equation (Yabana and Bertsch, 1996) (or its
variant, the Liouville-von Neumann equation based on density matrix), the so-
called _real-time_ TDDFT (RT-TDDFT). This approach is non-perturbative _in
principle_ , and has the ability to capture all the non-linear effects in a
photo-chemical and physical process. These include non-linear optical
properties (Ding et al., 2013), photo-induced coupled electron-ion dynamics
(Falke et al., 2014) namely, photo-dissociation and photo-isomerisation,
electronic quantum control, and strong-field processes such as above-threshold
ionization (Schafer et al., 1993), high harmonic generation (HHG) (Chu and
Chu, 2001), multi-photon ionization (Sakai et al., 2003), charge-resonance
enhanced-ionization (Bocharova et al., 2011), etc. With the advancement in
atto-second science and relevant spectroscopic techniques, one can now probe
directly the electron dynamics with femto- or atto-second temporal resolution.
That is why further development of theoretical as well as sophisticated
computational methods are essential to complement these strong-field processes
(Ullrich and Bandrauk, 2012; Crawford-Uranga et al., 2014; Nisoli et al.,
2017; Kruchinin et al., 2018). Motivated by the rapid developments in ultra-
fast experimental techniques, RT-TDDFT is drawing more and more attention as a
reliable method of choice to simulate the strong-field phenomena.

In order to solve the TDKS equations, one, _in principle_ , needs a full time-
dependent propagator (TDP) which can be written as,

$\mathrm{\mathbf{U}}(t,0)=\mathcal{T}\mathrm{exp}\bigg{\\{}-i\int_{0}^{t}dt^{\prime}H_{\mathrm{s}}(\rho,t^{\prime})\bigg{\\}},$
(1)

where $\mathcal{T}\mathrm{exp}$ is the time-order exponential containing TDKS
Hamiltonian, $H_{\mathrm{s}}(\rho,t)$. But in practice, it has been used for a
short time of interval, and then the full-TDP can be rewritten as,

$\displaystyle\mathrm{\mathbf{U}}(t_{1},t_{2})=\mathrm{\mathbf{U}}(t_{1},t_{3})\mathrm{\mathbf{U}}(t_{3},t_{2})$
$\displaystyle\mathrm{\mathbf{U}}(t,0)=\prod_{i=0}^{N-1}\mathrm{\mathbf{U}}(t_{i}+\Delta
t_{i},t_{i}).$ (2)

In general, the TDP is unitary, $\mathrm{\mathbf{U}}(t+\Delta
t,t)=\mathrm{\mathbf{U}}^{-1}(t+\Delta t,t)$, and has the time-reversal
symmetry, $\mathrm{\mathbf{U}}(t+\Delta
t,t)=\mathrm{\mathbf{U}}^{-1}(t,t+\Delta t)$, for a Hermitian Hamiltonian,
which is obvious for $H_{\mathrm{s}}(\rho,t)$. As its exact form is non-
trivial to implement as yet, it has to be approximated, and accordingly
various _approximate time-dependent propagators_ (ATDP) have been introduced
in the literature, such as exponential midpoint and modified midpoint
algorithm, enforced time-reversal symmetry method, Magnus expansion, splitting
techniques, etc (Castro et al., 2004; Li et al., 2005; Watanabe and Tsukada,
2002). Recently, the so called classical based ATDPs have been utilized in RT-
TDDFT calculations which were largely ignored in the past (Gomez Pueyo et al.,
2018). This includes linear multi-step, Runge-Kutta, exponential Runge-Kutta
and their different versions. The accuracy and efficiency of a RT-TDDFT
calculation depends critically on the level of ATDP and spectral
characteristic of TDKS Hamiltonian. These have been implemented successfully
using plane waves (Schleife et al., 2012), numerical orbitals (Soler et al.,
2002; Hekele et al., 2021), Gaussian-type orbitals (Li et al., 2005; Lopata
and Govind, 2011; Nguyen and Parkhill, 2015), and on real-space grids (Andrade
et al., 2015) for a wide range of RT-TDDFT applications.

We, however, note that RT-TDDFT is not used widely as the method of choice in
literature due to its huge computational burden requiring a smaller time step
to propagate the density for a period of time (e.g., tens of femtoseconds at a
minimum). On the other hand, there are no such alternatives which can
accurately describe the non-linear electron dynamics of a large quantum
mechanical system. Besides, high-lying excitations are attractive targets for
RT-TDDFT which are very expensive for linear response algorithm due to the
iterative nature of implementation. In recent years, most developments have
been found to reduce the computational cost while maintaining the accuracy and
numerical stability of ATDP. In this direction, we offer a novel propagation
scheme based on _adiabatic eigenstate subspace_ (AES) technique that holds the
promise of a larger time step. This advocates the use of a split-operator
technique in energy representation (Son and Chu, 2009), to perform the ATDP in
AES. We have used a density predictor/corrector algorithm to maintain the
self-consistency in $\rho(\mathrm{\mathbf{r}},t)$ at each time step. Most of
the elements of TDKS matrix are directly computed in Cartesian coordinate grid
(CCG). The relevant TD properties are computed through TD dipole-moment
(TDDM). To demonstrate the internal consistency of our proposed scheme, we
consider the interaction of a molecule with a high-intensity ultra-short
strong laser field, taking H2 and N2 as representative. The TDDM and HHG
spectrum are calculated using _adiabatic local density approximation_ (ALDA),
along with a comparison of available theoretical references (Apra et al.,
2020). This would establish the main objective of present communication to
offer a general framework of RT-TDDFT of many-electron systems for strong-
field processes. Towards this goal, we present a RT-TDDFT implementation in
our in-house code InDFT (A. K. Roy, A. Ghosal and T. Mandal, 2019). The plan
of this article is as follows. In next section, we discuss the structure of
TDKS equation and its representation in AES, along with the propagation
scheme. Then, we highlight some of the TD quantities in CCG. In Sec. III, the
necessary computational and technical details are offered. In Sec. IV, we
demonstrate the internal consistency of this approach in strong-field regimes
considering TDDM and HHG spectrum. For H2, an _all-electron_ calculation is
performed, whereas a _pseudopotential_ approximation is invoked in case of N2.
Finally, a few conclusions as well as future and outlook are summarized in
last section.

##  II Methodology

###  II.1 Time-dependent KS equation

Similar to ground state KS-DFT (Kohn and Sham, 1965), the TD non-interacting
system is described by a single Slater determinant of fictitious orbitals
($\bm{\varphi}(t)=|\varphi_{1}(\mathrm{\mathbf{r}},t),\varphi_{2}(\mathrm{\mathbf{r}},t),\dots,\varphi_{j}(\mathrm{\mathbf{r}},t)|$)
and the interacting TD electron density $\rho(\mathrm{\mathbf{r}},t)$ is
represented as,

$\rho(\mathrm{\mathbf{r}},t)=\sum_{j=1}^{N}|\varphi_{j}(\mathrm{\mathbf{r}},t)|^{2},$
(3)

where $N$ is the total number of occupied orbitals. Now, these orbitals are
obtained by solving the TDKS equation as below (henceforth atomic units
employed unless otherwise stated),

$\displaystyle i\frac{\partial}{\partial
t}\varphi_{j}(\mathrm{\mathbf{r}},t)=\text{H}_{\mathrm{s}}(\mathrm{\mathbf{r}},t)\varphi_{j}(\mathrm{\mathbf{r}},t),$
(4)
$\displaystyle\text{H}_{\mathrm{s}}(\mathrm{\mathbf{r}},t)=\big{[}-\frac{1}{2}\nabla^{2}+v_{\mathrm{s}}(\mathrm{\mathbf{r}},t)\big{]},$
$\displaystyle
v_{\mathrm{s}}(\mathrm{\mathbf{r}},t)=v_{\mathrm{ext}}(\mathrm{\mathbf{r}},t)+v_{\mathrm{h}}(\mathrm{\mathbf{r}},t)+v_{\mathrm{xc}}(\mathrm{\mathbf{r}},t).$

Here $v_{\mathrm{s}}(\mathrm{\mathbf{r}},t)$ signifies the one-body TDKS
potential. It includes the external potential,
$v_{\mathrm{ext}}(\mathrm{\mathbf{r}},t)$, the classical Coulomb or Hartree
potential, $v_{\mathrm{h}}(\mathrm{\mathbf{r}},t)$ and the TD exchange-
correlation (XC) potential, $v_{\mathrm{xc}}(\mathrm{\mathbf{r}},t)$. We
restrict ourselves to spin-unpolarized systems for simplicity.

###  II.2 Adiabatic eigenstate subspace

The TDKS Hamiltonian, $\text{H}_{\mathrm{KS}}$ consists of both linear and
non-linear terms. The former including only the kinetic term, is considered
for possible _stiffness_ of TDKS equation. Once the system is excited from
ground state, the TD density will oscillate with a wide range of frequencies.
The Fourier components of those oscillations correspond to the Bohr
frequencies, including both valence-to-valence and core-to-valence
excitations. To maintain the numerical stability due to presence of such
highly oscillating Fourier components, Eq. $(4)$ has to be solved using a
sufficiently small time step, which is the main integrated part of ATDP.
Recently, AES approaches with varying flavor have been introduced in the
literature (Chu and Chu, 2001; Son and Chu, 2009; Russakoff et al., 2016) to
increase the time steps and hence reducing the computational cost.
Accordingly, the TDKS orbitals can be expanded as,

$\displaystyle\varphi_{j}(\mathrm{\mathbf{r}},t)=\sum_{l}a_{jl}(t)\phi_{l}(\mathrm{\mathbf{r}}),$
(5) $\displaystyle a_{jl}(t)=\delta_{jl}\quad\mathrm{at}\quad t=0.$

where $\phi_{l}(\mathrm{\mathbf{r}})$ denotes the solution of time-independent
KS (TIKS) equations at its ground state. This forms the basis of AES and the
expansion coefficients, $a_{jl}$(t) projects the TDKS orbitals into this
subspace.

In order to utilize the advantage of AES, we recast the Hamiltonian in
following way (Son and Chu, 2009),

$\text{H}_{\mathrm{s}}(\mathrm{\mathbf{r}},t)=\text{H}_{\mathrm{s}}^{\mathrm{st}}+\text{H}_{\mathrm{s}}^{\mathrm{dy}}.$
(6)

where
$\text{H}_{\mathrm{s}}^{\mathrm{st}}(\mathrm{\mathbf{r}})=-\frac{1}{2}\nabla^{2}+v_{\mathrm{s}}(\mathrm{\mathbf{r}},0)$
and
$\text{H}_{\mathrm{s}}^{\mathrm{dy}}(\mathrm{\mathbf{r}},t)=v_{\mathrm{ks}}(\mathrm{\mathbf{r}},t)-v_{\mathrm{ks}}(\mathrm{\mathbf{r}},0)$.
The first one is called stationary KS Hamiltonian corresponding to ground
state. This would remain constant throughout the propagation. The second one
stands for dynamic KS Hamiltonian which effectively contains the fluctuation
embedded in TD potential. It should be updated at each instant of time during
propagation. This kind of splitting may help us to utilize the outcome of AES
approach fully, and consequently may improve the computational efficiency.

We are now able to explore the effect of AES on TDKS equations in a well-
defined matrix form by exploiting Eqs. $(4)$ to $(6)$. After some mathematical
manipulation, we can recast Eq. $(4)$ in terms of TD expansion coefficients in
following fashion,

$i\frac{\partial a_{jk}}{\partial
t}=a_{jk}(t)\epsilon_{k}+\sum_{l}a_{jl}(t)\langle\phi_{k}|\text{H}_{\mathrm{s}}^{\mathrm{dy}}|\phi_{l}\rangle,$
(7)

where $a_{jk}$ is the $k$th coefficient of $j$th TDKS orbital and
$\epsilon_{k}$ denotes the eigenvalue of $k$th TIKS orbital from which the AES
is formed. Further, we can also set up all the equivalent equations of the
remaining coefficients for $j$th TDKS orbitals in a well-defined matrix form

$i\frac{\partial\mathrm{\mathbf{A}}_{j}}{\partial
t}=\big{(}\bm{\epsilon}+\mathrm{\mathbf{H}}_{\mathrm{s}}^{\mathrm{dy}}\big{)}\mathrm{\mathbf{A}}_{j}.$
(8)

Here, $\mathrm{\mathbf{A}}_{j}$ stands for a column matrix for the $j$th TDKS
orbital and $\bm{\epsilon}$ is a diagonal matrix containing the eigenvalues
corresponding to TIKS orbitals. In general,
$\mathrm{\mathbf{H}}_{\mathrm{s}}^{\mathrm{dy}}$ is a complex Hermitian matrix
and, therefore, it must be square symmetric (Lopata and Govind, 2011; Zhu and
Herbert, 2018). Again, it can be easily generalized to the spin-polarized
case, as given below,

$i\frac{\partial\mathrm{\mathbf{A}}_{j,\sigma}}{\partial
t}=\big{(}\bm{\epsilon}_{\sigma}+\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}\big{)}\mathrm{\mathbf{A}}_{j,\sigma},\quad\quad\sigma=\alpha~{}\mathrm{or}~{}\beta~{}\mathrm{spin}.$
(9)

All the quantities are defined exactly in the same way as Eq. $(8)$, but
contain the appropriate spin part. Similarly, the total TD electron density
would have contributions from both $\alpha$ and $\beta$ densities i.e.,
$\rho(\mathrm{\mathbf{r}},t)=\rho_{\alpha}(\mathrm{\mathbf{r}},t)+\rho_{\beta}(\mathrm{\mathbf{r}},t)$
and
$\rho_{\sigma}(\mathrm{\mathbf{r}},t)=\sum_{j}|\varphi_{j,\sigma}(\mathrm{\mathbf{r}},t)|^{2}$.
Now, Eq. $(9)$ is the real workhorse in our present spin-polarized RT-TDDFT
calculations.

###  II.3 Time-evolution through split-operator

In order to solve Eq. $(9)$, we use a second-order split-operator technique in
energy representation as an ATDP (Chu and Chu, 2001; Son and Chu, 2009).
Accordingly, it can be expressed as,

$\mathrm{\mathbf{U}}_{\sigma}(t+\Delta
t)=\bigg{[}\mathrm{exp}(-i\frac{\bm{\epsilon}_{\sigma}}{2}\Delta
t)\mathrm{exp}(-i\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}\Delta
t)\mathrm{exp}(-i\frac{\bm{\epsilon}_{\sigma}}{2}\Delta t)\bigg{]}+O(\Delta
t^{3}),$ (10)

where the dynamic KS potential term appears sandwiched between two TIKS
orbital energies. The reverse of this approach (KS orbital energy term
sandwiched between two dynamical KS potential terms) is also equally
legitimate, but it requires two times more computational cost than the former.
It is to be noted that this form of ATDP would be unitary and unconditionally
stable provided all the exponential terms are computed exactly. It provides a
very reliable source of second-order method. Now, the solution of Eq. $(9)$
can be expressed as,

$\displaystyle\mathrm{\mathbf{A}}_{j,\sigma}(t+\Delta
t)=\mathrm{\mathbf{U}}_{\sigma}(t+\Delta
t)\mathrm{\mathbf{A}}_{j,\sigma}(t),\quad$
$\displaystyle\mathrm{\mathbf{A}}_{j,\sigma}(t+\Delta
t)=\mathrm{\mathbf{E}}_{\sigma}\mathrm{\mathbf{S}}_{\sigma}\mathrm{\mathbf{E}}_{\sigma}\mathrm{\mathbf{A}}_{j,\sigma}(t),\quad$
(11)
$\displaystyle\mathrm{\mathbf{E}}_{\sigma}=\mathrm{exp}(-i\frac{\bm{\epsilon}_{\sigma}}{2}\Delta
t)\quad\mathrm{and}\quad\mathrm{\mathbf{S}}_{\sigma}=\mathrm{exp}(-i\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}\Delta
t).$

At first glance, Eq. $(11)$ seems to be a chain of simple matrix-vector
multiplication, but it involves exponential of $\bm{\epsilon}_{\sigma}$ and
$\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}$ in
$\mathrm{\mathbf{E}}_{\sigma}$ and $\mathrm{\mathbf{S}}_{\sigma}$
respectively. The exponential of $\bm{\epsilon}_{\sigma}$ is trivially
obtained by exponentiating every entry on the main diagonal. It becomes time-
independent, if uniform temporal spacing is used throughout the propagation.
On the other hand, it is not so trivial to obtain
$\mathrm{\mathbf{S}}_{\sigma}$, the exponential of
$\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}$. This can be performed
using a broad range of techniques such as diagonalization, power series
method, Lanczos algorithm, etc. Here, we have employed a spectral
decomposition of $\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}$ in
AES (i.e.,
$\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}=\mathrm{\mathbf{C}}_{\sigma}\mathrm{\mathbf{Q}}_{\sigma}\mathrm{\mathbf{C}}^{\dagger}_{\sigma}$)
and accordingly, it is found that,

$\mathrm{\mathbf{S}}_{\sigma}=\mathrm{exp}(-i\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}\Delta
t)=\mathrm{\mathbf{C}}_{\sigma}\mathrm{exp}(-i\mathrm{\mathbf{Q}}_{\sigma}\Delta
t)\mathrm{\mathbf{C}}^{\dagger}_{\sigma}.$ (12)

The diagonalization of $\mathrm{\mathbf{H}}_{\mathrm{s},\sigma}^{\mathrm{dy}}$
is straightforward to perform in AES through numerical routine, and that makes
this split-operator approach unitary and stable. We, however, note that it is
quite different from traditional split-operator techniques where
$\mathrm{\mathbf{S}}$ contains kinetic energy operator, usually solved in
momentum space.

###  II.4 Time-dependent quantities

Now we discretize various quantities like TD electron density, TIKS and TDKS
orbitals as well as dynamics KS Hamiltonian on a 3D simulation box having
$x,y,z$ axes, as below,

$\displaystyle r_{i}=r_{0}+(i-1)h_{r},\quad i=1,2,3,....,N_{r}~{},\quad
r_{0}=-\frac{N_{r}h_{r}}{2},\quad r\in\\{x,y,z\\},$ (13)

where $h_{r},N_{r}$ denote grid spacing and total number of points along each
direction. The TD electron density $\rho(\mathrm{\mathbf{r}},t)$ in this grid
may be simply written as (“g” symbolizes discretized grid),

$\rho(\mathrm{\mathbf{r}}_{g},t)=\sum_{j=1}^{N}f_{j}|\varphi_{j}(\mathrm{\mathbf{r}}_{g},t)|^{2},\quad\varphi_{j}(\mathrm{\mathbf{r}}_{g},t)=\sum_{l}a_{jl}(t)\phi_{l}(\mathrm{\mathbf{r}}_{g}).$

At this stage, the contributions of dynamic KS matrix are directly computed in
CCG,

$\displaystyle\mathrm{\mathbf{H}}_{\mathrm{s},ij}^{\mathrm{dy}}=\langle\phi_{i}(\mathrm{\mathbf{r}})|\text{H}_{\mathrm{s}}^{\mathrm{dy}}(\mathrm{\mathbf{r}},t)|\phi_{j}(\mathrm{\mathbf{r}})\rangle$
$\displaystyle\mathrm{\mathbf{H}}_{\mathrm{s},ij}^{\mathrm{dy}}=h_{x}h_{y}h_{z}\sum_{g}\phi_{i}(\mathrm{\mathbf{r}}_{g})\text{H}_{\mathrm{s}}^{\mathrm{dy}}(\mathrm{\mathbf{r}}_{g},t)\phi_{j}(\mathrm{\mathbf{r}}_{g}),$
(14)

and subsequently total TD norm, $N(t)$ is also computed in CCG as,

$N(t)=h_{x}h_{y}h_{z}\sum_{g}\rho(\mathrm{\mathbf{r}}_{g},t).$ (15)

The construction of various TD potentials in
$\text{H}_{\mathrm{s}}^{\mathrm{dy}}(\mathrm{\mathbf{r}}_{g},t)$ within the
adiabatic approximation at a given time can be calculated analogous to their
time-independent counterpart, and the details have been well documented in our
earlier works (Roy, 2008a, b, 2011; Ghosal and Roy, 2016; Ghosal et al., 2018;
Mandal et al., 2019; Ghosal et al., 2019; Roy et al., 2021a, b); hence not
repeated here. Here we need to consider only the TD potential arising from an
electric field.

Figure 1: Electric field of the laser pulse with carrier frequency,
$\omega_{0}=1.55$ eV and intensity, $I=10^{14}$ W cm-2.

The TD electric field potential, $v_{\mathrm{f}}(\mathrm{\mathbf{r}},t)$ is
introduced to represent the laser-matter interaction; in particular we are
interested in laser-molecule interaction. It is considered to be a TD
perturbation in $v_{\mathrm{ext}}(t)$ besides Coulombic attractive potential.
During such interaction, the system is excited by applied electric field,
$\mathrm{\mathbf{E}}(t)$ through transition dipole moment integral. Within the
“length gauge” (or dipole) approximation, the component of
$\mathrm{\mathbf{H}}_{\mathrm{s}}^{\mathrm{dy}}$ due to
$v_{\mathrm{f}}(\mathrm{\mathbf{r}},t)$ can be expressed as,

$\displaystyle
v_{\mathrm{f}}(\mathrm{\mathbf{r}},t)=-\mathrm{\mathbf{E}}(t)\bm{\cdot}\mathrm{\mathbf{r}},$
$\displaystyle\mathrm{\mathbf{v}}_{\mathrm{f},ij}=-\bigg{[}\mathrm{\mathbf{D}}_{ij}\bm{\cdot}\mathrm{\mathbf{E}}(t)\bigg{]},$
(16)

where $\mathrm{\mathbf{D}}$ stands for transition dipole moment integral in
AES. Here, we have mainly considered those properties which are derived from
TDDM, $\mu(t)$. This can be computed directly from the transition dipole
moment integral and TD density matrix, $\mathrm{\mathbf{P}}(t)$. Accordingly,
the $z$-polarized component can be expressed as,

$\mu_{z}(t)=\mathrm{Tr}[\mathrm{\mathbf{D}}_{z}\mathrm{\mathbf{P}}(t)].$ (17)

The elements of $\mathrm{\mathbf{D}}_{z}$ in AES can be written as follows,

$\displaystyle\mathrm{\mathbf{D}}_{z,ij}=\langle\phi_{i}(\mathrm{\mathbf{r}})|z|\phi_{j}(\mathrm{\mathbf{r}})\rangle,$
$\displaystyle\mathrm{\mathbf{D}}_{z,ij}=\sum_{\nu\lambda}C_{\nu i}C_{\lambda
j}\langle\chi_{\nu}(\mathrm{\mathbf{r}})|z|\chi_{\lambda}(\mathrm{\mathbf{r}})\rangle,$
(18)

where the first expression is further expanded in atomic orbital bases (or
contracted bases) due to the fact that each element involved in Eq. $(18)$ can
be computed analytically using well-known Obara-Saika recursion relation
(Obara and Saika, 1986). Moreover, the element of $\mathrm{\mathbf{P}}(t)$ in
AES can be expressed in the following way,

$\mathrm{\mathbf{P}}_{ij}(t)=\sum_{k}f_{k}a_{ik}(t)a_{jk}(t),$ (19)

where $a_{ik}$ is the $i$th coefficient of $k$th TDKS orbital.

Figure 2: The temporal dependence of TDDM along z-direction of H2 molecule
over the full duration of laser pulse, along with reference result (Set IX).

Before passing, it may be worth noting that the propagation scheme is not
restricted to TIKS orbitals. One can use different subspaces consisting of
atomic orbitals, numerical orbitals etc., but that may require additional cost
during diagonalization of $\text{H}_{\mathrm{s}}^{\mathrm{st}}$ on that
subspace. Here, our approach facilitates the removal of high-frequency
oscillations generated by static Hamiltonian through the split-operator
technique in energy representation. In this direction, it is worthwhile to
mention the work of (Jia et al., 2018) which is based on a parallel transport
gauge transformation of orbitals such that the transformed orbitals, more
specifically the residual vectors, are expected to oscillate less frequently
than the original ones. But they do need additional term in the Hamiltonian in
order to satisfy the desired transformation. On the other hand, the present
approach is rather simple and it works because the spurious high frequencies
of oscillation are damped by the
$\mathrm{exp}(-i\frac{\bm{\epsilon}_{\sigma}}{2}\Delta t)$ term. Therefore the
electron dynamics always follows the slowest possible (physical) oscillation.

##  III Computational details

For laser-molecule interaction, the electric field, $\mathrm{\mathbf{E}}(t)$
is conveniently written as,

$\mathrm{\mathbf{E}}(t)=\mathrm{\mathbf{n}}f(t)\ \mathrm{sin}(\omega_{0}t),$
(20)

where $\mathrm{\mathbf{n}}$ is the unit vector in the direction of electric
field polarization vector, $f(t)$ is the envelope function; $\omega_{0}$ is
the carrier frequency. Here, we employ a $\mathrm{cos}^{2}$ shape laser field,
as below,

$f(t)=\begin{cases}f_{0}~{}\mathrm{cos}^{2}\left(\frac{\pi}{2\kappa}(\kappa-t)\right)&\text{if
}|t-\kappa|\leq\kappa,\\\ 0&\text{otherwise}.\end{cases}$

where $\kappa$ signifies the width of envelope. We choose the particular laser
field with carrier frequency, $\omega_{0}=1.55$ eV (corresponding to a
Ti:sapphire laser), intensity, $I=10^{14}$ W cm-2, and $\kappa=10$ optical
cycle (o.c.), where $1$ o.c. $=\frac{2\pi}{\omega_{0}}$. Now, this particular
form of laser parameters are adopted from (Luppi and Head-Gordon, 2012)
because it permits us to neglect the ionization effects. Then, the propagation
is continued up to a maximum duration of $40$ o.c., corresponding to
t${}_{f}=4412.80$ a.u. ($106.74$ fs), where pulse duration is $20$ o.c. (see
Fig. $1$).

The spectrum of HHG is computed by taking the Fourier transform of the TDDM:

$H_{z}(\omega)=\Bigg{|}\frac{1}{t_{f}-t_{i}}\int_{t_{i}}^{t_{f}}\mu_{z}(t)e^{-i\omega
t}dt\Bigg{|}^{2}.$ (21)

We consider the applied TD field along the molecular axis ($z$-direction) and
the corresponding induced TDDM is measured along the same direction. In order
to perform discrete Fourier transform, standard FFTW3 package (Frigo and
Johnson, 2005) was invoked.

The dimension of $\mathrm{\mathbf{S}}$ is $N_{b,\mathrm{AES}}$ where
$N_{b,\mathrm{AES}}$ defines the total number of eigenstates in AES. The
computational cost of evaluation of exponential via diagonalization scales as
$N_{b,\mathrm{AES}}^{3}$. It incurs a negligible cost compared to
$\mathrm{\mathbf{H}}_{\mathrm{s}}^{\mathrm{dy}}$ formation, scaling as
$N_{b,\mathrm{AES}}^{2}N_{g}$, where $N_{g}$ defines the total number of grid
points in CCG. Moreover, Eq. $(11)$ involves a series of matrix-vector
multiplication scaling as $N_{b,\mathrm{AES}}^{2}$. Hence, construction of
dynamic KS matrix remains the most computationally expensive step than any
other operations. Further, it also requires twice the memory footprint of
ground-state DFT calculation as it becomes complex. Therefore, the numerous
build of dynamic KS matrix may enhance the computational cost. In order to
overcome this issue, we have used density predictor/corrector algorithm
additionally for improving the efficiency in terms of number of
$\mathrm{\mathbf{H}}_{\mathrm{s}}^{\mathrm{dy}}$ builds per unit time of
propagation. This is performed self-consistently by employing a convergence
criteria during time evolution i.e.,
$\mathrm{max}\\{|\rho^{\mathrm{new}}_{\sigma}-\rho_{\sigma}|\\}<\iota$, where
$\iota$ is about $1\times 10^{-7}$. The resulting generalized matrix-
eigenvalue problem and matrix diagonalization is solved through standard
LAPACK routine (Anderson et al., 1999) accurately and efficiently. These are
implemented in our in-house program InDFT (A. K. Roy, A. Ghosal and T. Mandal,
2019); all the results presented in this communication are computed from
there.

Table 1:  Numerical details of our proposed RT-TDDFT scheme in CCG (all in a.u.). See text for details. Set | $h_{r}$ | $\Delta t$ | Basis set | XC functional |  $\\{N_{x},N_{y},N_{z}\\}$(optimal)
---|---|---|---|---|---
I | 0.2 | 0.2 | 6-311G | L(S)DA | $\\{80,80,80\\}$
II | ,, | 0.1 | ,, | ,, | ,,
III | ,, | 0.05 | ,, | ,, | ,,
IV | ,, | 0.025 | ,, | ,, | ,,
V | 0.3 | 0.1 | SBKJC | ,, | $\\{52,52,60\\}$
VI | ,, | 0.08 | ,, | ,, | ,,
VII | ,, | 0.04 | ,, | ,, | ,,
VIII | ,, | 0.01 | ,, | ,, | ,,

At first, we solve the TIKS equations, for a given basis set and XC
functional. The corresponding TIKS orbitals form the AES which lies at the
heart of our propagation scheme. A systematic grid optimization technique from
(Ghosal and Roy, 2016; Ghosal et al., 2018; Mandal et al., 2019; Ghosal et
al., 2019; Roy et al., 2021a, b), provides an optimal grid which would be the
starting point of our RT-TDDFT calculations. It is well justified in terms of
current propagation scheme. We have performed the calculations at fixed
experimental geometry, taken from NIST database (Johnson, 2016), employing the
following finite atomic orbital basis sets: _all-electron_ basis $6$-$311$G
(includes only $s$-type of functions) for H2 and ECP basis SBKJC (Stevens et
al., 1984) for N2 respectively. For _all-electron_ calculations, we chose the
particular one to avid any node near the nucleus that would facilitate the
successful implementation of KS-DFT in CCG. These are adopted from EMSL Basis
Set Library (Feller, 1996). For TDXC functional, we employ spin-polarized LDA
in its ground state form within adiabatic approximation including Slater spin
exchange and Vosko-Wilk-Nusair (VWN)-V parametrized model for correlation
within a homogeneous electron gas approximation (Vosko et al., 1980).

Figure 3: The temporal dependence of TDDM along z-direction of H2 molecule
toward the peak duration of pulse, along with the reference result (Set IX).
See text for details.

##  IV Results and Discussion

In the present context, the primary objective is to introduce a simple CCG-
based RT-TDDFT scheme in dealing with laser-molecule interaction, where
perturbative models are not feasible. To the end, we would also discuss the
effectiveness of pseudopotential RT-TDDFT calculation in simulating strong-
field processes through a specimen case of N2. Let us begin with the case of
H2 molecule. At first, the computed TDDM along $z$ direction is displayed in
Figs. 2 and 3, for different temporal resolution. The first one gives it for
full duration of the simulation, while the second one provides it towards the
end of the pulse respectively. We project out the effect of temporal
resolution on the accuracy by comparing them with reference results. The
latter is identified as Set IX, referring to 6-311G basis set, spin-polarized
LDA XC functional, $\Delta t=0.1$ and other default parameters in NWChem (Apra
et al., 2020). The details regarding optimal grid, spatial resolution (or grid
spacing), temporal resolution, $\Delta t$ are given in Table I for
convenience. The polarization of the field was chosen to be the $z$-axis,
along molecular axis.

Figure 4: The temporal dependence of HHG spectra along z-direction of H2
molecule compared with reference result (Set IX). See text for details.

We begin our computation with a temporal resolution of $0.2$ a.u., in Set I
which is twenty times larger than the one used in (Luppi and Head-Gordon,
2012). This is then halved successively from Set II to Set IV. From Fig. $2$,
it is well manifested that TDDM coherently follows the $20$ o.c. of laser
pulse irrespective of the set of interest, which corroborates well with (Luppi
and Head-Gordon, 2012). Further, they also correctly follow the behavioral
pattern of TDDM from reference (Set IX), but the amplitude of oscillation at
different peak positions are somehow different. Figure 3 tells that, the
oscillation amplitude at the end of $20$ o.c. does not return back to zero. It
is quite possible due to non-adiabatic excitations during the propagation, as
discussed in (Schlegel et al., 2007), but in our current scenario, we do not
have enough flexibility in basis set to describe these excitations. We,
however, believe that this small fluctuation in amplitude is due to numerical
errors accumulated during the propagation. To analyze it further, we follow
Fig. $3$ very closely and it is observed that there is a sharp decrease in
amplitude of oscillations than the frequency of laser field, and remains
constant after Set II. Accordingly, we set the time step corresponding to Set
III as an optimal choice for comparing the results with reference (Set IX).

Figure 5: The temporal dependence of TD dipole moment along z-direction of N2
molecule over the full duration of the laser pulse.

Next, we present the HHG spectra from TDDM, in Fig. $4$. Clearly, the
distinctive shape of HHG spectra is obtained, which consists of a
characteristic shape: first a rapid decrease for the low-order harmonics,
followed by an extended plateau region, and finally a cut-off region. As
expected, only odd harmonics are observed, because a homo-nuclear diatomic
like H2 has a center of inversion. All the HHG spectra have roughly comparable
intensity profiles as expected from dipole. If we compare these spectra (Sets
I-IV), with reference (Set IX), the intensity profiles are different, but all
contain the same physics. The effects of smaller time step, as observed in
dipole moment plot, are also visible in the characteristic of cut-off region.
The results of Sets I and II do not correctly reproduce the cut-off region of
the spectra. To pursue further, we have calculated Keldysh parameter (Keldysh,
1965) theoretically for the given basis set and XC functional, which is $0.93$
and the corresponding maximum harmonic cut-off is $22$. At this point, it is
important to note that the calculated value of Keldysh parameter and maximum
HHG cut-off is based on the assumption of a single active electron in the
three-step model. It is well manifested by Sets I and II, having a close
resemblance with reference result (Set IX). But for Sets III and IV, the cut-
off position moves towards the higher harmonics, which is also found in recent
real-time wave function based calculation of H2 (Luppi and Head-Gordon, 2012).
We believe that the performance of our current calculations using spin-
polarized LDA XC functional along with small basis set such as $6$-$311$G in
describing HHG spectra may further be improved.

Figure 6: The temporal dependence of HHG spectra along z-direction of N2
molecule.

In order to analyze the stability of propagator and its dependence on time
step, we note that the norm of wave function was conserved in all Sets (I-IV)
during the propagation. Further, the exponential density predictor/corrector
algorithm took mostly two steps (sometimes three) to facilitate a larger time,
maintaining the accuracy and numerical stability. Moreover, it may be possible
to resolve the numerical artefacts by controlling spatial resolution to a
certain extent, but it has an adverse effect on computational cost, due to an
increase in number of grid points.

In a broader context, the focal point of this work is to deal with larger
systems and this is facilitated by the use of ECP in the realm of gaussian
basis set. The pseudopotential approximation would be a useful option for
successful implementation of RT-TDDFT in CCG. In such calculations (Foglia et
al., 2017), this may offer substantial advantages: (1) it reduces the
dimension of TDDM, which may significantly alleviate the cost of propagation
scheme (2) it also increases the value of required time-step, by eliminating
the need to integrate the motion of inner electrons, associated with high-
energy transition and fast oscillating frequencies. Thus, the relevance of
pseudopotential approximation may be significant when dealing with heavy
atoms, not only because of a reduction of dimension of TDDM, but also through
a considerable decrease in propagation time.

Towards this direction, we choose N2 molecule. At first, computation is
performed with a temporal resolution of $0.1$ a.u., with Set V which is ten
times larger than the one used in (Luppi and Head-Gordon, 2012), and then it
is gradually decreased from Set VI to Set VIII. From Fig. $5$, it is clear
that TDDM coherently follows the $20$ o.c. of laser pulse for all the sets of
interest. Also Fig. $5$ shows the same trend at the end of $20$ o.c., as found
in case of H2. Again, here also, the conclusion remains same and it is
possibly due to some numerical artefacts accumulated during propagation.
Furthermore, there is a sharp decrease in amplitude of oscillations than the
frequency of laser field. We have calculated the HHG spectra from TDDM, which
is reported in Fig. $6$. The distinctive shape of HHG spectra in low-order
harmonics and plateau region is well obtained, but the sharp cut-off region is
missing. In a closer look, odd harmonics are there in the lower and plateau
region. All the HHG spectra have roughly comparable intensity as expected from
TDDM, and the dependence of time step is not so prominent here. Moreover, the
numerical noise present in spectra is very hard to remove after filtering
through window function, and that causes considerable complexity in the
analysis of spectra. At this stage, it is difficult to distinguish the
structure and signal in the HHG spectra from the background. We have
calculated Keldysh parameter (Keldysh, 1965) theoretically for the given basis
set and XC functional, which is $1.04$ (tunneling ionization) and the
corresponding maximum harmonic cut-off is $21$. It is not sharply manifested
in all the sets. In a recent article (Foglia et al., 2017), it has been shown
that RT-TDDFT calculations for absorption spectra with pseudopotential
approximation can be quite sensitive to the choice of basis function. As most
of the basis sets designed for ECP are typically constructed to reproduce
ground-state properties, this may not be as flexible for TDDFT calculations,
which inherently requires both ground and excited states. Currently, we are
looking into it.

A few remarks may be made at this stage. The present work is a continuation of
our ongoing work to establish a viable alternative for DFT calculations, for
small to medium systems, through a CCG based pseudopotential KS-DFT framework
using LCAO-MO ansatz as implemented in InDFT (A. K. Roy, A. Ghosal and T.
Mandal, 2019). Here, we have made an attempt to establish a general framework
of RT-TDDFT for atoms/molecules allowing reasonable trade-off between accuracy
and efficiency. Therefore, it is our intention to offer the performance of our
scheme in CCG. It is to be noted that the present propagation scheme is not
restricted to CCG. One can use atom-centered grid such as Lebedev grid
(Lebedev and Vychisl, 1975, 1976) with Becke’s partitioning (Becke and
Dickson, 1988) to construct the TD Hamiltonian in accordance with the modern
computational code. However, as evident from the previous discussion, the
applicability of the present approach is crucially dependent on various
factors, especially basis set; however the present propagation scheme remains
quite reliable and accurate. Thus it is our firm belief that the current
propagation scheme is much more general than its current mode of
implementation, and has considerable promise for electron dynamics of strong-
field processes for atoms/molecules.

##  V Future and outlook

We have demonstrated the viability and suitability of an AES based RT-TDDFT
approach for electron dynamics in strong-field processes, through both all-
electron and pseudopotential approximation. This was applied for H2 and N2 as
prototypical cases. The properties derived from TDDM were offered. Besides
these, their performances on different time steps were also presented.

The most important conclusion is that in finite systems, it shows considerable
promise for electron dynamics of strong-field processes for larger systems.
The success of this approach relies on accurate implementation of our
numerical propagator in CCG as well as estimation of all dynamical KS matrix
elements. We primarily focused on the dependence of time step on TDDM and its
corresponding HHG spectra, for which very limited information is available in
CCG, so far.

The added flexibility of basis set is very important here. So far, we have
taken in to account only a particular kind of basis set designed for ECP
calculations. But, other basis sets derived from diffusive quantum Monte-
Carlo, may also throw some light on the quality of these spectra. This issue
leaves ample room for future code improvement.

##  VI Availability of data

Data available on request from the authors.

##  VII Acknowledgement

AG is grateful to UGC for a senior research fellowship. AKR sincerely
acknowledges funding from DST SERB, New Delhi, India (sanction order:
CRG/2019/000293). Partial financial assistance from BRNS, Mumbai (sanction
order: 58/14/03/2019-BRNS/10255) is acknowledged. A special thank goes to Dr.
N. Govind for numerous valuable and constructive comments. AG thanks Prof.
Vikas and Dr. E. Luppi for useful discussions.

## References

  * Casida and Huix-Rotllant (2012) M. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem. 63, 287 (2012).
  * Marques et al. (2012) M. A. L. Marques, F. M. S. Nogueira, E. K. U. Gross, and A. Rubio, _(Eds.) Fundamentals of Time-Dependent Density Functional Theory, Vol. 837_ (Springer-Verlag, Heidelberg, 2012).
  * Adamo and Jacquemin (2013) C. Adamo and D. Jacquemin, Chem. Soc. Rev. 42, 845 (2013).
  * Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
  * Casida (1995) M. E. Casida, in _Recent Advances in Computational Chemistry; D. P. Chong (Ed.)_ (World Scientific, NJ, 1995), pp. 155–192.
  * Yabana and Bertsch (1996) K. Yabana and G. F. Bertsch, Phys. Rev. B 54, 4484 (1996).
  * Ding et al. (2013) F. Ding, B. E. Van Kuiken, B. E. Eichinger, and X. Li, J. Chem. Phys. 138, 064104 (2013).
  * Falke et al. (2014) S. M. Falke, C. A. Rozzi, D. Brida, M. Maiuri, M. Amato, E. Sommer, A. De Sio, A. Rubio, G. Cerullo, E. Molinari, et al., Science 344, 1001 (2014).
  * Schafer et al. (1993) K. J. Schafer, B. Yang, L. F. DiMauro, and K. C. Kulander, Phys. Rev. Lett. 70, 1599 (1993).
  * Chu and Chu (2001) X. Chu and S.-I. Chu, Phys. Rev. A 64, 063404 (2001).
  * Sakai et al. (2003) H. Sakai, J. J. Larsen, I. Wendt-Larsen, J. Olesen, P. B. Corkum, and H. Stapelfeldt, Phys. Rev. A 67, 063404 (2003).
  * Bocharova et al. (2011) I. Bocharova, R. Karimi, E. F. Penka, J.-P. Brichta, P. Lassonde, X. Fu, J.-C. Kieffer, A. D. Bandrauk, I. Litvinyuk, J. Sanderson, et al., Phys. Rev. Lett. 107, 063201 (2011).
  * Ullrich and Bandrauk (2012) C. A. Ullrich and A. D. Bandrauk, in _Fundamentals of Time-Dependent Density Functional Theory; M. A. L. Marques, N. T. Maitra, F. M. S. Nogueira, E. K. U. Gross and A. Rubio (Eds.)_ (Springer, Berlin, 2012), pp. 351–371.
  * Crawford-Uranga et al. (2014) A. Crawford-Uranga, U. De Giovannini, E. Räsänen, M. J. Oliveira, D. J. Mowbray, G. M. Nikolopoulos, E. T. Karamatskos, D. Markellos, P. Lambropoulos, S. Kurth, et al., Phys. Rev. A 90, 033412 (2014).
  * Nisoli et al. (2017) M. Nisoli, P. Decleva, F. Calegari, A. Palacios, and F. Martín, Chem. Rev. 117, 10760 (2017).
  * Kruchinin et al. (2018) S. Y. Kruchinin, F. Krausz, and V. S. Yakovlev, Rev. Mod. Phys. 90, 021002 (2018).
  * Castro et al. (2004) A. Castro, M. A. Marques, and A. Rubio, J. Chem. Phys. 121, 3425 (2004).
  * Li et al. (2005) X. Li, S. M. Smith, A. N. Markevitch, D. A. Romanov, R. J. Levis, and H. B. Schlegel, Phys. Chem. Chem. Phys. 7, 233 (2005).
  * Watanabe and Tsukada (2002) N. Watanabe and M. Tsukada, Phys. Rev. E 65, 036705 (2002).
  * Gomez Pueyo et al. (2018) A. Gomez Pueyo, M. A. Marques, A. Rubio, and A. Castro, J. Chem. Theory and comput. 14, 3040 (2018).
  * Schleife et al. (2012) A. Schleife, E. W. Draeger, Y. Kanai, and A. A. Correa, J. Chem. Phys. 137, 22A546 (2012).
  * Soler et al. (2002) J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Condense. Matter 14, 2745 (2002).
  * Hekele et al. (2021) J. Hekele, Y. Yao, Y. Kanai, V. Blum, and P. Kratzer, J. Chem.Phys. 155, 154801 (2021).
  * Lopata and Govind (2011) K. Lopata and N. Govind, J. Chem. Theory Comput. 7, 1344 (2011).
  * Nguyen and Parkhill (2015) T. S. Nguyen and J. Parkhill, J. Chem. Theory Comput. 11, 2918 (2015).
  * Andrade et al. (2015) X. Andrade, D. Strubbe, U. De Giovannini, A. H. Larsen, M. J. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. J. Verstraete, et al., Phys. Chem. Chem. Phys. 17, 31371 (2015).
  * Son and Chu (2009) S.-K. Son and S.-I. Chu, Chem. Phys. 366, 91 (2009).
  * Apra et al. (2020) E. Apra, E. J. Bylaska, W. A. De Jong, N. Govind, K. Kowalski, T. P. Straatsma, M. Valiev, H. J. van Dam, Y. Alexeev, J. Anchell, et al., J. Chem. Phys. 152, 184102 (2020).
  * A. K. Roy, A. Ghosal and T. Mandal (2019) A. K. Roy, A. Ghosal and T. Mandal, _InDFT: A DFT Program for Atoms and Molecules in CCG_ (Theoretical Chemistry Laboratory, IISER Kolkata, India, 2019), This is based on the extension of an initial version of the code, established by A. K. Roy, in 2008, whose results were published in refs. [33-35].
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Russakoff et al. (2016) A. Russakoff, Y. Li, S. He, and K. Varga, J. Chem. Phys. 144, 204125 (2016).
  * Zhu and Herbert (2018) Y. Zhu and J. M. Herbert, J. Chem. Phys. 148, 044117 (2018).
  * Roy (2008a) A. K. Roy, Int. J. Quant. Chem. 108, 837 (2008a).
  * Roy (2008b) A. K. Roy, Chem. Phys. Lett. 461, 142 (2008b).
  * Roy (2011) A. K. Roy, J. Math. Chem. 49, 1687 (2011).
  * Ghosal and Roy (2016) A. Ghosal and A. K. Roy, _In Specialist Periodical Reports: Chemical Modelling, Applications and Theory; M. Springborg and J.-O. Joswig (Eds.) Vol. 13_ (Royal Society of Chemistry, London, 2016).
  * Ghosal et al. (2018) A. Ghosal, T. Mandal, and A. K. Roy, Int. J. Quant. Chem. 118, e25708 (2018).
  * Mandal et al. (2019) T. Mandal, A. Ghosal, and A. K. Roy, Theor. Chem. Acc. 138, 10 (2019).
  * Ghosal et al. (2019) A. Ghosal, T. Mandal, and A. K. Roy, J. Chem. Phys. 150, 064104 (2019).
  * Roy et al. (2021a) R. Roy, A. Ghosal, and A. K. Roy, Chem. Asian J. 16, 2729 (2021a).
  * Roy et al. (2021b) R. Roy, A. Ghosal, and A. K. Roy, J. Phys. Chem. A (in press) (2021b).
  * Obara and Saika (1986) S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986).
  * Jia et al. (2018) W. Jia, D. An, L.-W. Wang, and L. Lin, J. Chem. Theory and Comput. 14, 5645 (2018).
  * Luppi and Head-Gordon (2012) E. Luppi and M. Head-Gordon, Mol. Phys. 110, 909 (2012).
  * Frigo and Johnson (2005) M. Frigo and S. G. Johnson, Proceedings of the IEEE 93, 216 (2005).
  * Anderson et al. (1999) E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Dongarra, J. D. C. A. Greenbaum, S. Hammarling, A. A. McKenney, and D. Sorensen, _LAPACK Users’ Guide_ , vol. 9 (SIAM, 1999).
  * Johnson (2016) R. D. Johnson, _III (Ed.) NIST Computational Chemistry Comparisons and Benchmark Database, NIST Standard Reference Database, Number, Release 18_ (NIST, Gaithersburg, MD, 2016).
  * Stevens et al. (1984) W. J. Stevens, H. Basch, and M. Krauss, J. Chem. Phys. 81, 6026 (1984).
  * Feller (1996) D. Feller, J. Comp. Chem. 17, 1571 (1996).
  * Vosko et al. (1980) S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
  * Schlegel et al. (2007) H. B. Schlegel, S. M. Smith, and X. Li, J. Chem.Phys. 126, 244110 (2007).
  * Keldysh (1965) L. Keldysh, Sov. Phys. JETP 20, 1307 (1965).
  * Foglia et al. (2017) N. O. Foglia, U. N. Morzan, D. A. Estrin, D. A. Scherlis, and M. C. Gonzalez Lebrero, J. Chem. Theory Comput. 13, 77 (2017).
  * Lebedev and Vychisl (1975) V. I. Lebedev and Z. Vychisl, Mat. Mat. Fiz. 15, 48 (1975).
  * Lebedev and Vychisl (1976) V. I. Lebedev and Z. Vychisl, Mat. Mat. Fiz. 16, 293 (1976).
  * Becke and Dickson (1988) A. D. Becke and R. M. Dickson, J. Chem. Phys. 89, 2993 (1988).


# Latent-CF: A Simple Baseline for Reverse Counterfactual Explanations

Rachana Balasubramanian  Samuel Sharpe  Brian Barr  Jason Wittenbach  C. Bayan
Bruss

###### Abstract

In the environment of fair lending laws and the General Data Protection
Regulation (GDPR), the ability to explain a model’s prediction is of paramount
importance. High quality explanations are the first step in assessing
fairness. Counterfactuals are valuable tools for explainability. They provide
actionable, comprehensible explanations for the individual who is subject to
decisions made from the prediction. It is important to find a baseline for
producing them. We propose a simple method for generating counterfactuals by
using gradient descent to search in the latent space of an autoencoder and
benchmark our method against approaches that search for counterfactuals in
feature space. Additionally, we implement metrics to concretely evaluate the
quality of the counterfactuals. We show that latent space counterfactual
generation strikes a balance between the speed of basic feature gradient
descent methods and the sparseness and authenticity of counterfactuals
generated by more complex feature space oriented techniques.

##  1 Introduction

The GDPR provides a "right to be informed" about system functionality in
automated decision making processes. Though not explicit in the GDPR, the law
in spirit ostensibly encourages those creating algorithmic decision making
systems to build trust and increase transparency around these systems [21].
When considering the fair application of artificial intelligence in industrial
settings such as finance, it is important to consider the information and
agency afforded to those affected by the decisions of the AI system.

In particular, providing the affected party with a means to either change
their behavior to ensure a better outcome in the future, or to contest the
current outcome. One way to provide this is through counterfactual
explanations. Specifically, counterfactual explanations provide a plausible
change to a model’s input that would result in a change in the model’s
prediction.

Counterfactual explanations fall within the category of model-agnostic post-
hoc example based explanations. Counterfactual explanations trace their
lineage to philosophy and psychology and are deeply intertwined with analyses
of causality. Philosophers proposed the use of counterfactuals as a way to
understand causal arguments. If A caused B then if A had not occurred B would
not have occurred [11]. Much of the subsequent work in causal inference takes
this same approach to counterfactual analysis. Changes to antecedents should
result in changes to outcomes otherwise there is no causal relationship [17].

Researchers in psychology have noted the importance of counterfactual thinking
in human cognition, with development as early as two years old [4]. Unlike the
counterfactuals in philosophy, the functional theory of counterfactual
thinking suggests a reverse process. Namely, a person starts with a mismatch
between expectations and the present outcome. The counterfactual isn’t to
understand generic causal relationships, but rather it is goal oriented to
ensure a different outcome in the future. In this framing, one starts from the
undesirable outcome and works backward until a path opens up to the desirable
outcome. This fits more closely with the framework of reverse causal inference
[6]. Reverse causal inference starts from a deviation from expectations -
typically an anomaly - and performs model checking to understand the causes of
the deviation. Because this backward looking process can result in a number of
possible antecedent changes, philosophers and psychologists have begun to
explore the process by which individuals prune the possible candidate
counterfactuals referred to as counterfactual availability [5].

In the domain of model-based counterfactual explorations, the common approach
is aligned to the traditional forward looking causal structure. Namely,
counterfactual approaches explore alternative inputs that would lead to
changes in the model’s predictions. The objective of a counterfactual is that
given $Y=f(X)$, for a specific $x_{0}$ where $f(x_{0})=y_{1}$ generate
$x_{cf}$ where $f(x_{cf})=y_{2}$. What is clear from this formulation is that
there may be many possible paths to generating $x_{cf}$. As a result, forward
looking counterfactuals seek to place constraints on the possible changes both
for computational efficiency but also interpretability [21].

In this paper we propose a method for generating counterfactual explanations
in-line with the reverse causal inference framework and the functional theory
of counterfactual thinking. These frameworks are well-suited to providing
those affected by a decision with an explanation, especially when the decision
led to an undesired outcome. The explanations produced provide a method of
changing behavior to reach a desired goal in such cases.

We refer to our approach as Latent-CF. Given any black-box classifier and a
dataset that has been used to train the classifier, Latent-CF can generate
counterfactuals. Training a separate autoencoder on the same training data is
necessary for our approach. For any data point and desired class confidence of
the target counterfactual class, we traverse the latent space of the
autoencoder from the original data point’s encoded representation until the
desired class probability is reached and then use the decoder to generate the
corresponding counterfactual. Because the latent space is more compact,
starting from the desired outcome and searching for counterfactuals in that
space is more computationally efficient than doing so in feature space. In
practical terms, the approach is simple but provides a foundation for viewing
counterfactual explanations from a more functional perspective compared to
existing forward oriented counterfactuals.

In order to evaluate the quality of these counterfactuals we use a framework
similar to that of [12] by stating that a counterfactual explanation should
be:

  1. 1.

In distribution

  2. 2.

Sparse in the number of changes it makes in the features

  3. 3.

Computationally efficient

These criteria align also to notions of counterfactual availability making
them more actionable [5].

We conduct the following experiments to evaluate the quality of the
counterfactuals generated by Latent-CF. First, using the criteria above, we
compare Latent-CF to approaches that generate counterfactual explanations
utilizing gradient descent in feature space varying loss and clipping
strategies to encourage the explanations to be in-sample and sparse. Second,
we conduct a visual analysis of the interpretability of the counterfactual
explanations. Through these experiments, we are able to show that Latent-CF
simple method provides a strong baseline to compare current (and future)
methods in reverse counterfactual generation.

##  2 Previous Work

There is a large body of existing work on post-hoc explainability, see [14]
for an overview. Many methods treat the creation of explanations as a problem
of credit assignment - for each input sample, provide the relative importance
of the input features to that samples prediction. Early methods had roots in
computer vision, where a common approach is to calculate sensitivities and
show a heat map depicting which pixels are responsible for an image’s
classification. These methods, e.g. [19], rely on using a zero information
baseline. The need for a baseline can become problematic if applied to tabular
or other data where the definition of a baseline may not exist.

The use of local surrogate models, such as LIME [18] and its extensions, can
overcome the need for global baselines to be established. However, the
perturbation method used does not ensure that the local surrogate is built on
in-sample data.

Game theoretic approaches [20], [1], [13] view the credit assignment task as a
coalitional game amongst features. These methods baseline a feature’s
contribution against the average model prediction. Game theoretic approaches
are challenged by the exponential number of coalitions that must be evaluated.

Other researchers have developed methods that express the commonality of
features the must be present or missing. Dhurandhar et al. [2] try to find
perturbations in feature space utilizing an autoencoder reconstruction loss to
keep the explanations in sample. Their method seeks to find contrastive
explanations, which describe what minimal set of features or characteristics
must be missing to explain why it does not belong to an alternate class.

Counterfactual explanations avoid the pitfalls of previous local attribution
methods (no need for baselines, no approximation to game theory constructs, no
need for universality in features). All they require is a sample in need of
explanation and a search method to find the decision boundary. Lash et al. [9]
maintain sparsity in their inverse classification methods by partitioning the
features into mutable and immutable categories, and imposing budgetary
constraints on allowable changes to the mutable features. Laugel et al. [10]
advocate a sampling approach with their growing spheres method.

Wachter et al. [21] proposed the use of counterfactual explanations to help
the individuals impacted by a model based decision why a particular decision
was reached, to provide grounds to contest adverse decisions, and to
understand what would result in an alternative decision. van Looveren [12]
explore loss functions to search for perturbations in feature space and
introduce prototype loss which encourages the counterfactual to be close to
the average representation in latent space of $K$ training samples from the
target class. In contrast to our method which searches for a single
counterfactual instance, Mothilal et al [15] suggest generating a set of
diverse counterfactuals. McGrath et al. [7] used counterfactuals to explain
credit application predictions. They introduced weights, based on either
global feature importance or nearest neighbors in an attempt to make them more
interpretable.

Some recent counterfactual generation techniques search directly in a latent
space. Pawelczyk et al. [16] construct a model agnostic technique that sample
observations uniformly in $l_{p}$ spheres around the original point’s
representation in latent space. Joshi et al. [8] utilize a very similar
approach to ours using a VAE or other generative models to learn $p(x)$ and
traverse the latent space with gradient descent to find counterfactuals. They
include an additional loss term in the feature space to encourage sparse
changes. We choose to use the most simple manifestation of latent space
generated counterfactuals—a basic autoencoder and loss including only score of
a decoded latent representation—to convincingly illustrate the benefits of
searching in a smaller representative space.

##  3 Methods

###  3.1 Pretrained Models

The classifier is built using Keras in Tensorflow 2. It uses four stacked
layers of convolution with a 3x3 kernel size, 2x2 pooling kernels and filter
depths = [32, 16, 8, 1] followed by a dense layer with sigmoid activation as
the output. It is trained on the standard train/test split of MNIST filtered
to only 4’s and 9’s.

###  3.2 Latent-CF

In feature space perturbation methods, a sample $x_{0}$ is fed through a
classifier $f:X\rightarrow[0,1]$ and iteratively updated until $f(x)$ is close
to $0.5$ (decision boundary) or some other target probability $p$ where the
closeness to the boundary is measured by a loss function $\mathcal{L}$.
Approaches in [12] and [2] incorporate an autoencoder in $\mathcal{L}$ to
guide the search for counterfactual examples. However, Latent-CF, similar to
[8], directly perturbs the latent representation of the encoder, $z=E(x)$,
until the probability of the decoded sample, $f(D(z))$, is close to $p$. Our
algorithm is detailed in Algorithm 1 and the architecture of Latent-CF is
illustrated in Figure 1.

parameter : $p$ probability of target counterfactual class (0.5 for decision
boundary), $tol$ tolerance

Input: Instance to explain $x_{0}$, classifier $f$, encoder $E$, decoder $D$,
$\mathcal{L}$ loss function

Output: $x_{cf}$ the counterfactual explanation

Encode sample to latent space $z=E(x_{0})$

Calculate $loss=\mathcal{L}(p,f(D(z))$

while  _$loss >tol$_ do

$z\leftarrow z-\alpha\nabla_{z}loss$

Calculate $loss=\mathcal{L}(p,f(D(z))$

$x_{cf}=D(z)$

Algorithm 1 Latent-CF Figure 1: Comparing approaches for counterfactual
generation between gradient descent in feature space (left) and Latent-CF
gradient descent in the autoencoder latent space (right)

###  3.3 Baseline Methods

Our first baseline method, Feature GD, uses gradient descent to directly
perturb the feature space minimizing the $\ell_{2}$ distance to the decision
boundary. Two other methods introduce some small changes to make iterative
improvements over Feature GD. First, we add feature clipping after every
gradient step in Feature GD + clip to ensure pixel values stay within the
correct domain.

We implement Feature GD + MAD loss, which encourages in-sample counterfactuals
and sparse changes with Median Absolute Deviation (MAD) scaling loss, as
developed by Wachter et al [21]. Instead of the $\ell_{2}$ loss, this method
uses the $\ell_{1}$ norm weighted by the inverse median absolute deviation,
such that $MAD_{k}$ of a feature $k$ over the set of points $P$ is as shown in
equation 1.

$MAD_{k}=median_{j\in{P}}(|X_{j,k}-median_{l\in{P}}(X_{l,k})|)$ (1)

This results in a distance metric $d$ for between a synthetic data point
$x^{\prime}$ and original data point $x$ as described by equation 2.

$d(x,x^{\prime})=\sum_{k\in{F}}\frac{|x_{k}-x^{\prime}_{k}|}{MAD_{k}}$ (2)

This distance metric encourages changing only the pixels that vary in the
training set, and punishes changing the border pixels that are rarely
nonzero.The architecture for these three methods is illustrated on the left of
Figure 1.

We compare a final feature gradient descent method from [12] which we label
Prototype-CF. The authors include five different loss terms in their objective
to achieve desirable properties of their counterfactuals.

$L=cL_{pred}+\beta L_{1}+L_{2}+L_{AE}+L_{proto}$

$L_{pred}$ is designed to encourage counterfactuals of a different class.
$L_{1}$ and $L_{2}$ are combined to form an elastic-net regularizer on the
feature perturbations for sparse changes. They include $L_{AE}$ as used by
Dhurandhar et al. in [3] to ensure in-sample reconstructions of
counterfactuals. Finally, they guide the search for counterfactuals by
introducing prototypes, which are the average euclidean representation of a
class in the latent space defined by the $K$ closest encoded samples to
$E(x_{0})$. Specifically, for a class $i$ the prototype is defined as

$proto_{i}=\frac{1}{K}\sum_{k=1}^{K}E(x_{k}^{i})$

where
$\left\|E(x_{k}^{i})-E(x_{0})\right\|_{2}\leq\left\|E(x_{k+1}^{i})-E(x_{0})\right\|_{2}$.
$L_{proto}$, defined as the distance to the closest prototype, effectively
tries to speed up the search for counterfactuals by pushing $x_{0}+\delta$ to
the closest prototype.

###  3.4 Evaluation Metrics

In order to evaluate the quality of the counterfactual, we use three
evaluation metrics:

1) In-distribution
The counterfactual should not change the image such that proposed features
will have a low probability of occurring. We used per-pixel intensity kernel
density estimation (KDE) to measure the extent to which our counterfactuals
are in-sample. Specifically, we estimate the density over intensity values for
each pixel across the target class population. We compute the probability for
each pixel intensity given the pixel specific KDE and take the average over
all pixels. An example of this can be seen in Figure 2 comparing probabilities
of pixels in counterfactuals as well as probabilities of the original pixels.

Figure 2: Histogram of Kernel Density Estimate probabilities per pixel
computed on original images and counterfactuals. In this case the
counterfactual images generated by baseline methods have more low probability
pixels than the original images, signifying out of distribution changes.

2) Sparsity
The counterfactual should be sparse in the number of changes it makes in the
features. We compute the percent of features that are changed in the
generation of the counterfactual.

3) Computational Efficiency
Generating the counterfactual must be computationally efficient. We measure
the time it takes to generate a counterfactual.

###  3.5 Experiments

In order to generate a baseline comparison of the counterfactual methods, we
analyze the MNIST digit classification dataset. We performed a series of
experiments. We first compare each method outlined above as well as the method
proposed in [12] on the four criteria proposed. As mentioned in Section 3.1,
we use a convolutional neural network for the binary classification task of
distinguishing fours and nines. Given this classifier, we generate
counterfactuals for fours given nines and vice versa. We also conduct a visual
comparison of the counterfactual images generated by each model type, and use
these along with our evaluation metrics to draw conclusions about the quality
of each method.

##  4 Results and Discussion

Method | In-Distribution | Sparsity (%) | Time (s)
---|---|---|---
Feature GD | 0.098 | 77.5 | 0.346
Feature GD + clip | 0.111 | 48.7 | 0.401
Feature GD + MAD loss | 0.117 | 38.4 | 0.364
Prototype-CF | 0.133 | 16.8 | 15.400
Latent-CF | 0.135 | 24.2 | 1.134
Table 1: Mean and standard deviation for quality metrics given a target
counterfactual class probability of 0.5: in-distribution, sparsity, and time
for all methods in this study.

In Table 1, we present a comparison of each method across all of the proposed
metrics. Each of these metrics were calculated using a fixed probablity,
$p=0.5$, as the decision boundary. We will investigate each metric in further
detail later on, but this serves to point out some immediate trends. First, we
can notice that by far Prototype-CF method takes much longer than the Feature
GD methods. Secondly, we notice that Latent-CF and the Prototype-CF method
change a much smaller percentage of the pixels, making their generated
counterfactuals sparser. Finally, Prototype-CF and Latent-CF result in less
out-of-distribution changes to the image resulting in in-distribution scores
close the original images themselves (0.137).

(a) Counterfactual run time for each method

(b) Fraction of pixels changed by each method

(c) In-distribution probability over all test image pixels

Figure 3: Average and 95% CI of time, sparsity, and in-distribution
probability across all test images calculated for various counterfactual
confidence levels.

In Figure 3 we show these performance metrics when pushing beyond the decision
boundary of 0.5 to higher levels confidence. We only include results for
Prototype-CF at 0.5 since their formulation uses a penalty to encourage
bounding the confidence from below. Figure 3(c) shows how the probability of
being in sample changes for different confidence values. Models that perform
gradient descent in feature space tend to drift more out of sample as the
classifier seeks to obtain higher confidence in the counterfactual image,
while the Latent-CF method stays relatively stable. We notice that for most
Feature GD counterfactuals, the changes are somewhat adversarial in that they
make a significant number of perturbations around the border of the image.
These changes get larger in magnitude as we target higher confidences, thus
decreasing average pixel in-distribution probability.

To get a quantitative evaluation of how each model and method performed in
terms of sparsity, we calculated the fraction of the 784 pixels that were
changed to get from the original image to the counterfactual image. Figure
3(b) demonstrates how the proportion of pixels changed varies with increasing
confidence.

(a) Original Image

(b) Latent-CF

(c) feature gradient descent

(d) feature gradient descent with clipping

(e) MAD loss

(f) Prototype-CF

Figure 4: Visual inspection of counterfactual generation by each of the four
methods proposed here. The heat map corresponds to the changes in intensities
of each pixel. Note smaller scales for feature gradient descent methods.

Finally, efficiency is important when generating counterfactuals for
production systems. Table 1 shows the time in seconds that it took for each
model type to generate a counterfactual with a given method. Unsurprisingly,
the feature gradient descent methods tend to run the fastest and methods that
add additional transformations, architectures, or losses increases the
computation necessary and slows down counterfactual generation. We also see
these effects across levels of counterfactual confidence in Figure 3(a).
Latent-CF strikes a balance between all the elements of a quality
counterfactual generation method. Latent-CF demonstrates the advantages of
Prototype-CF while also resulting in an order of magnitude speed up.

We also provide a visual examination of counterfactuals in Figure 4. We start
with an image classified as a 4 (Figure 4(a)) and generate counterfactual
images with each method. We include the counterfactual image and heat map
which illustrates the addition or subtraction of luminosity for a pixel in the
counterfactual image compared to the original. These are presented for each of
the four described methods in Figures 4(b) to 4(f). From the heat maps, we can
observe that the initial attempt using feature gradient descent is noisy and
does not seem to be confined to any features of the 4. Once we add clipping,
we do observe a less noisy heat map, with aspects of the 4 becoming visible,
but there are still many noisy changes in surrounding pixels. Using MAD loss
results in similar changes as clipping. Finally, we can observe that Latent-CF
and Prototype-CF actually begin to modify the shape of the four, adding
curvature and closing the top of the of the four as could be expected in a
nine. We can see a significant change in the heat map as Latent-CF directly
manipulates the shape of the four as changes are concentrated on the border of
the four while Prototype-CF makes some extraneous changes inside and outside
the main shape.

##  5 Conclusion

We demonstrate the benefits of performing perturbations in a representative
latent space compared to various methods in feature space for counterfactual
generation. We show that these methods benefit from sparsity and in-sample
perturbations that simple feature based methods lack. Latent based generation
also exhibits significant speed up compared to more complex feature based
methods designed to remedy simple methods’ difficulty generating sparse and
in-sample counterfactuals.

## Broader Impact

A positive impact of this work is an increased focus in future research on
latent perturbation methods to produce higher quality counterfactuals for the
benefit of transparency and parties affected by model decisions. Possible
negative outcomes of this work is the presence of bias in data and/or
classifiers that are not addressed in this approach, and thus, may persist in
the counterfactual explanations themselves. The failure of counterfactual
systems to provide accurate and fair explanations could lead to
counterproductive decisions by affected parties further worsening their
standing in the eyes the decision system.

## References

  * [1] A. Datta, S. Sen, and Y. Zick.  Algorithmic transparency via quantitative input influence: Theory and experiments with learning systems.  In 2016 IEEE Symposium on Security and Privacy (SP), pages 598–617, May 2016.
  * [2] Amit Dhurandhar, Pin-Yu Chen, Ronny Luss, Chun-Chen Tu, Paishun Ting, Karthikeyan Shanmugam, and Payel Das.  Explanations based on the missing: Towards contrastive explanations with pertinent negatives.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS’18, page 590–601, Red Hook, NY, USA, 2018. Curran Associates Inc.
  * [3] Amit Dhurandhar, Vijay Iyengar, Ronny Luss, and Karthikeyan Shanmugam.  Tip: Typifying the interpretability of procedures, 2018.
  * [4] Kai Epstude and Neal J Roese.  The functional theory of counterfactual thinking.  Personality and social psychology review, 12(2):168–192, 2008.
  * [5] Sina Fazelpour.  Norms in counterfactual selection.  Philosophy and Phenomenological Research, 06 2020.
  * [6] Andrew Gelman and Guido Imbens.  Why ask why? forward causal inference and reverse causal questions.  Working Paper 19614, National Bureau of Economic Research, November 2013\.
  * [7] Rory Mc Grath, Luca Costabello, Chan Le Van, Paul Sweeney, Farbod Kamiab, Zhao Shen, and Freddy Lécué.  Interpretable credit application predictions with counterfactual explanations.  CoRR, abs/1811.05245, 2018.
  * [8] Shalmali Joshi, Oluwasanmi Koyejo, Warut Vijitbenjaronk, Been Kim, and Joydeep Ghosh.  Towards realistic individual recourse and actionable explanations in black-box decision making systems, 2019.
  * [9] Michael T. Lash, Qihang Lin, W. Nick Street, and Jennifer G. Robinson.  A budget-constrained inverse classification framework for smooth classifiers.  In Proceedings of The IEEE International Conference on Data Mining (ICDM), ICDM’17, pages 1184–1193, New Orleans, Louisiana, 2017. IEEE.
  * [10] Thibault Laugel, Marie-Jeanne Lesot, Christophe Marsala, Xavier Renard, and Marcin Detyniecki.  Comparison-based inverse classification for interpretability in machine learning.  In Jesús Medina, Manuel Ojeda-Aciego, José Luis Verdegay, David A. Pelta, Inma P. Cabrera, Bernadette Bouchon-Meunier, and Ronald R. Yager, editors, Information Processing and Management of Uncertainty in Knowledge-Based Systems. Theory and Foundations, pages 100–111, Cham, 2018. Springer International Publishing.
  * [11] David Lewis.  Causation.  The Journal of Philosophy, 70(17):556–567, 1973.
  * [12] Arnaud Van Looveren and Janis Klaise.  Interpretable counterfactual explanations guided by prototypes, 2019.
  * [13] Scott M. Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In Isabelle Guyon, Ulrike von Luxburg, Samy Bengio, Hanna M. Wallach, Rob Fergus, S. V. N. Vishwanathan, and Roman Garnett, editors, Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, 4-9 December 2017, Long Beach, CA, USA, pages 4765–4774, 57 Morehouse LaneRed HookNYUnited States, 2017. Curran Associates Inc.
  * [14] Christoph Molnar.  Interpretable Machine Learning.  2019\.  https://christophm.github.io/interpretable-ml-book/.
  * [15] Ramaravind Kommiya Mothilal, Amit Sharma, and Chenhao Tan.  Explaining Machine Learning Classifiers through Diverse Counterfactual Explanations.  Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency, pages 607–617, January 2020.
  * [16] Martin Pawelczyk, Klaus Broelemann, and Gjergji Kasneci.  Learning model-agnostic counterfactual explanations for tabular data.  Proceedings of The Web Conference 2020, Apr 2020.
  * [17] Judea Pearl.  The algorithmization of counterfactuals.  Annals of Mathematics and Artificial Intelligence, 61(1):29, 2011\.
  * [18] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  Why should i trust you?: Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 1135–1144, New York, NY, USA, 2016\. ACM, Association for Computing Machinery.
  * [19] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje.  Learning important features through propagating activation differences.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 3145–3153. JMLR.org, 2017.
  * [20] Erik Ã…Â trumbelj and Igor Kononenko.  An efficient explanation of individual classifications using game theory.  Journal of Machine Learning Research, 11(1):1–18, 2010.
  * [21] Sandra Wachter, Brent Mittelstadt, and Chris Russell.  Counterfactual explanations without opening the black box: Automated decisions and the gdpr.  Harv. JL & Tech., 31:841, 2017.