# Simulating two-dimensional correlation spectroscopies with third-order
infrared and fifth-order infrared–Raman processes of liquid water

Hideaki Takahashi  0000-0001-6465-2049 Yoshitaka Tanimura  0000-0002-7913-054X
<EMAIL_ADDRESS>Department of Chemistry, Graduate School of
Science, Kyoto University, Kyoto 606-8502, Japan

(Last updated: )

###### Abstract

To investigate the possibility of measuring the intermolecular and
intramolecular anharmonic coupling of balk water, we calculate third-order
two-dimensional (2D) infrared (IR) spectra and fifth-order 2D
IR–IR–Raman–Raman spectra expressed in terms of four-body correlation
functions of optical observables. For this purpose, a multimode Brownian
oscillator model of four interacting anharmonic oscillators strongly coupled
to their respective heat baths is employed. The nonlinearity of the
system–bath interactions is considered to describe thermal relaxation and
vibrational dephasing. The linear and nonlinear spectra are then computed in a
non-Markovian and nonperturbative regime in a rigorous manner using the
discretized hierarchical equations of motion in mixed Liouville–Wigner space
(DHEOM–MLWS). The calculated 2D spectra for stretching–bending,
bending–librational, stretching–librational, and stretching–translational
modes consist of various positive and negative peaks exhibiting essential
details of the intermolecular and intramolecular mode–mode interactions under
thermal relaxation and dephasing at finite temperature.

##  I INTRODUCTION

The delineation of the spectroscopic line shapes of liquid water into
vibrational modes and their interactions has been the subject of numerous
experimental and theoretical studies.[1, 2, 3, 4] The uniqueness of water
arises from its high-frequency intramolecular modes that promote bond forming
and bond breaking through complex hydrogen bonding (HB)[5, 6, 7, 8] and its
low-frequency intermolecular modes that realize irreversible nuclear
motion.[9, 10, 11, 12] The interplay between these modes plays key roles in
chemical reaction processes.[13, 14]

Experimentally, infrared (IR)[15, 16] and third-order off-resonant Raman
spectroscopies[17, 18, 19] have been utilized for the analysis of the
vibrational motions of liquid water. The observables of these methods are
defined in terms of the dipole $\hat{\bm{\mu}}$ and polarizability
$\hat{\bm{\Pi}}$ as
$R_{\mathrm{IR}}^{(1)}(t_{1})=i\langle[\hat{\bm{\mu}}(t_{1}),\hat{\bm{\mu}}]\rangle/\hbar$
and
$R_{\mathrm{Raman}}^{(3)}(t_{1})=i\langle[\hat{\bm{\Pi}}(t_{1}),\hat{\bm{\Pi}}]\rangle/\hbar$,
respectively, as functions of a single time duration $t_{1}$, which is
referred to as one-dimensional (1D) spectroscopy.[20] The IR process has been
utilized not only for investigations of the intramolecular OH stretching
motion ($\sim$3700 cm-1) and the intramolecular HOH bending motion ($\sim$1600
cm-1), but also those of the HB intermolecular translation (vibrational)
motion (33$\sim$100 cm-1) and inter-HB molecular vibrational motion
(33$\sim$400 cm-1), which are also regarded as THz processes. The Raman
process, which exploits the large polarizability of water molecules, has been
used to probe both intramolecular and intermolecular modes. This approach is
particularly useful for investigations of intermolecular modes because of its
superior time resolution and intensity in such a low-frequency region.
Although 1D spectroscopic approaches are useful for the classification of
vibrational line shapes of liquids into intermolecular and intramolecular
modes, it is unclear whether the peaks that we investigate are mutually
coupled, and whether the width of the peaks that we measure are of an
inhomogeneous origin or a homogeneous origin, as these peaks are usually broad
and overlap.

Two-dimensional (2D) experiments, in which an additional time correlation is
imprinted on the system response, are capable of determining anharmonic
coupling among vibrational modes and homogeneous linewidths that are needed to
understand the degree of coupling to the bath (i.e., the surrounding molecules
with respect to the excited transition).[21, 20, 22, 23, 24] Recent
developments of experimental techniques have significantly augmented the study
of both intermolecular and intramolecular modes in the 0$\sim$4000 cm-1
frequency range through 2D IR,[25, 26, 27, 28, 29] 2D THz–Raman
spectroscopies,[30, 31, 32, 33] and 2D THz–IR–visible (TIV) spectroscopy,[34,
35] which is equivalent to 2D IR–IR–Raman spectroscopy.[36, 37] Although 2D IR
and 2D THz–Raman experiments provide information about intramolecular and
intermolecular degrees of freedom, respectively, neither approach can directly
measure the anharmonic coupling between intermolecular and intramolecular
modes. In contrast, 2D TIV, which utilizes light sources to excite both
intermolecular and intramolecular modes, allows us to investigate their
mode–mode coupling mechanism. Note that the analysis of both 2D IR-Raman and
2D THz–Raman spectroscopies defined in terms of three-body response functions
such as
$R_{\mathrm{IIR}}^{(3)}(t_{2},t_{1})=-\langle[[\hat{\bm{\Pi}}(t_{2}+t_{1}),\hat{\bm{\mu}}(t_{1})],\hat{\bm{\mu}}(0)]\rangle/\hbar^{2}$
is far more complex than that of third-order 2D IR spectroscopy defined in
terms of four-body correlation functions such as
$R_{\mathrm{2DIR}}^{(3)}(t_{3},t_{2},t_{1})=-i\langle[[\hat{\bm{\mu}}(t_{3}+t_{2}+t_{1}),[[\hat{\bm{\mu}}(t_{2}+t_{1}),\hat{\bm{\mu}}(t_{1})],\hat{\bm{\mu}}(0)]]\rangle/\hbar^{3}$
because the contributions from the nonlinear polarizability and dipole moments
involved in the laser interactions give rise to additional peaks, and
nondiagonal spectral peaks are not necessary to represent mode–mode coupling
peaks, as in the case of 2D IR spectra.[38, 39]

In this study, we extend our analysis of 2D THz–Raman[40, 41, 42, 43] and 2D
IR–Raman spectroscopies[36, 37] by including an additional time duration for
laser excitation and calculate the fifth-order 2D IR–IR–Raman–Raman spectrum
defined as
$R_{\mathrm{IIRR}}^{(5)}(t_{3},t_{2},t_{1})=-i\langle[[\hat{\bm{\Pi}}(t_{3}+t_{2}+t_{1}),[[\hat{\bm{\Pi}}(t_{2}+t_{1}),\hat{\bm{\mu}}(t_{1})],\hat{\bm{\mu}}(0)]]\rangle/\hbar^{3}$,
in addition to the third-order 2D IR spectrum defined as
$R_{\mathrm{2DIR}}^{(3)}(t_{3},t_{2},t_{1})$. Through simulation, we attempt
to clarify the roles of anharmonic mode–mode coupling and vibrational
dephasing, as well as thermal relaxation, among the intermolecular and
intramolecular modes of balk water.

It should be noted that simulating the 2D spectrum of liquid water for a
mixture of intramolecular and intermolecular modes is challenging, especially
on account of the peak-splitting in third-order 2D IR spectra caused by the
transitions between the $|0\rangle$ $\rightarrow$ $|1\rangle$ $\rightarrow$
$|2\rangle$ and $|0\rangle$ $\rightarrow$ $|1\rangle$ $\rightarrow$
$|0\rangle$ states, where $|n\rangle$ is the $n$th vibrational energy
state.[44, 45, 46] The high-frequency intramolecular motions must be addressed
within the framework of quantum dynamics theory, and their dynamics involve
complex hydrogen bond networks for 2D spectroscopic measurements. These
motions have only been studied from classical MD simulations[47, 46] or
classical MD-based quantum simulations, such as using stochastic theory[48]
and wave function–based theory.[49, 50] Therefore, the peak profiles of
experimentally obtained 2D IR spectra–particularly for mode–mode coupling
peaks[25, 26, 27, 28]–have not been accurately predicted from a theoretical
approach.

As a practical approach, the Brownian oscillator (BO) model has been employed
for the analysis of both linear[51, 52, 53] and nonlinear spectra.[21, 20, 22]
In this approach, the vibrational modes representing the spectroscopic
properties of interest are described as functions of molecular coordinates,
whereas the environmental molecular motions are described using heat baths
that exert thermal fluctuations and dissipation on the vibrational modes.
Although analyses based on MD simulations have shown that vibrational
relaxation and dephasing are important mechanisms for characterizing molecular
motions,[54, 55] the inclusion of the nonlinear and non-Markovian interactions
with the heat baths,[56, 57, 58] as well as the anharmonic mode–mode
interactions,[59, 22] are also significant. Therefore, for the vibrational
motions of liquid water, we introduced the multimode nonlinear BO model, which
consists of interacting anharmonic oscillators coupled to their respective
heat baths with linear–linear (LL) and square–linear (SL) system–bath (SB)
interactions to describe vibrational relaxation and dephasing.[43, 36, 37] To
guarantee the equilibrium state of the system at finite temperature, the model
must be solved accurately such that the dephasing (fluctuations) and
relaxation (dissipation) satisfy the quantum version of the fluctuation-
dissipation theorem.[60, 61, 62]

The model parameters were then chosen to fit the classical MD simulation for
2D IR–Raman spectra conducted using the polarizable water model for inter- and
intramolecular vibrational spectroscopy (POLI2VS);[63] we then modified the
anharmonicity of intramolecular modes to fit experimentally obtained 1D Raman
and IR spectra,[37] which agree with quantum MD results obtained using the
POLI2VS model[64] and with the ab-initio MD results that have recently become
available. [65, 66, 67] To simulate 2D correlation spectra that detect the
relaxation time of vibrational dephasing, it is necessary to deal with non-
perturbative and non-Markovian SB interactions, especially for intramolecular
modes, in which vibrational motion and the environment are quantum-
mechanically entangled.[44, 45] Although various hierarchical equations of
motion (HEOM) approaches have been developed for investigating nonlinear
spectroscopies,[60, 61, 62] here we employ discretized hierarchical equations
of motion in mixed Liouville–Wigner space (DHEOM–MLWS) specifically developed
for the water BO model.[37]

This paper is organized as follows. In Sec. II, we explain our multimode LL+SL
BO model and the parameter set used for simulations. In Sec. III, we give the
definitions of the nonlinear optical observables to be computed. In Sec. IV,
we present the results for third-order 2D correlation IR and fifth-order 2D
correlation IR–IR–Raman–Raman spectra based on two-mode calculations. Section
V is devoted to our concluding remarks.

##  II Nonlinear BO model for vibrational modes of Water

The details of the model and the DHEOM–MLWS for the quantum case have been
explained in Refs. 36, 37. Here, we present an overview of the model because
it is helpful for the analysis of calculated 2D spectra.

Table 1:  Parameter values of multimode LL+SL BO model for (1) stretching, (2)
bending, (3) librational, and (4) translational modes obtained on the basis of
Refs. 36, 37. The normalized parameters are defined as
$\tilde{\zeta}_{s}\equiv(\omega_{0}/\omega_{s})^{2}\zeta_{s}$,
$\tilde{V}_{LL}^{(s)}\equiv(\omega_{s}/\omega_{0})V_{LL}^{(s)}$,
$\tilde{V}_{SL}^{(s)}\equiv V_{SL}^{(s)}$,
$\tilde{g}_{s^{3}}\equiv(\omega_{s}/\omega_{0})^{3}g_{s^{3}}$,$\tilde{\mu}_{s}\equiv(\omega_{0}/\omega_{s})\mu_{s}$,
$\tilde{\mu}_{ss}\equiv(\omega_{0}/\omega_{s})^{2}\mu_{ss}$,
$\tilde{\Pi}_{s}\equiv(\omega_{0}/\omega_{s})\Pi_{s}$, and
$\tilde{\Pi}_{ss}\equiv(\omega_{0}/\omega_{s})^{2}\Pi_{ss}$, where we set the
fundamental frequency as $\omega_{0}=4000$ cm-1.

s |  $\omega_{s}$ (cm-1) | $\gamma_{s}/\omega_{0}$ | $\tilde{\zeta}_{s}$ | $\tilde{V}_{LL}^{(s)}$ | $\tilde{V}_{SL}^{(s)}$ | $\tilde{g}_{s^{3}}$ | $\tilde{\mu}_{s}$ | $\tilde{\mu}_{ss}$ | $\tilde{\Pi}_{s}$ | $\tilde{\Pi}_{ss}$
---|---|---|---|---|---|---|---|---|---|---
1 | $3520$ | $5.0{\times}10^{-3}$ | $9$ | $0$ | $1.0$ | $-5.0{\times}10^{-1}$ | $3.3$ | $1.2{\times}10^{-2}$ | $3.3$ | $2.5{\times}10^{-2}$
2 | $1710$ | $2{\times}10^{-2}$ | $0.8$ | $0$ | $1.0$ | $-7{\times}10^{-1}$ | $1.8$ | $0$ | $0.47$ | $-3.9{\times}10^{-2}$
3 | $390$ | $8.5{\times}10^{-2}$ | $8.3$ | $3.4{\times}10^{-3}$ | $1.0$ | $7{\times}10^{-3}$ | $21$ | $0$ | $2.1$ | $-0.83$
4 | $125$ | $0.5$ | $2.8$ | $2.8{\times}10^{-3}$ | $1.0$ | $9.7{\times}10^{-2}$ | $26$ | $2.1$ | $9.0$ | $2.3$

The multimode nonlinear BO model for liquid water consists of four primary
oscillator modes representing (1) intramolecular OH stretching (“stretching”),
(2) intramolecular HOH bending (“bending”), (3) HB–intermolecular librational
(“librational”), and (4) HB–intermolecular translational (“translational”)
motions.[36, 37] Note that the OH stretching modes consist of symmetry and
anti-symmetry modes, but here we treat the stretching modes as a single mode
because we construct the model based on the broadened linear and nonlinear
spectra as benchmarks, as in previous studies. Although it is possible to
model the two modes separately by MD simulation with the machine learning
approach,[68] solving such a model is computationally demanding, so here we
treat them as one mode; thus, they are described by dimensionless vibrational
coordinates $\vec{q}=(\dots,q_{s},\dots)$ with $s=1,\dots,4$ indexing the four
vibrational modes. The total Hamiltonian is then expressed as

$\displaystyle{H}$
$\displaystyle=\sum_{s}\left({H}^{(s)}_{\mathrm{S}}+{H}^{(s)}_{\mathrm{B}}+{H}^{(s)}_{\mathrm{I}}\right)$
$\displaystyle+\frac{1}{2}\sum_{s>s^{\prime}}\left(g_{s^{2}s^{\prime}}q_{s}^{2}q_{s^{\prime}}+g_{s{s^{\prime}}^{2}}q_{s}q_{s^{\prime}}^{2}\right),$
(1)

where

$\displaystyle{H}^{(s)}_{\mathrm{S}}=\frac{{p}_{s}^{2}}{2m_{s}}+\frac{1}{2!}m_{s}\omega_{s}^{2}q_{s}^{2}+\frac{1}{3!}g_{s^{3}}q_{s}^{3}$
(2)

is the Hamiltonian for the $s$th mode with mass $m_{s}$, coordinate ${q_{s}}$,
momentum ${p_{s}}$, fundamental frequency $\omega_{s}$, and anharmonicity of
potential $g_{s^{3}}$, and $g_{s^{2}s^{\prime}}$ and $g_{s{s^{\prime}}^{2}}$
are the anharmonic coupling between modes $s$ and $s^{\prime}$,

$\displaystyle{H}^{(s)}_{\mathrm{B}}=\sum_{j_{s}}\left(\frac{{p}_{j_{s}}^{2}}{2m_{j_{s}}}+\frac{m_{j_{s}}\omega_{j_{s}}^{2}{x}_{j_{s}}^{2}}{2}\right)+\sum_{j_{s}}\left(\frac{\alpha_{j}^{2}V_{s}^{2}({q_{s}})}{2m_{j_{s}}\omega_{j_{s}}^{2}}\right)$
(3)

is the bath Hamiltonian for the $s$th mode with the momentum, coordinate,
mass, and frequency of the $j_{s}$th bath oscillator given by ${p}_{j_{s}}$,
${x}_{j_{s}}$, $m_{j_{s}}$ and $\omega_{{j_{s}}}$, respectively, and

$\displaystyle{H}^{(s)}_{\mathrm{I}}$
$\displaystyle=-V_{s}({q_{s}})\sum_{j_{s}}\alpha_{j_{s}}{x}_{j_{s}}$ (4)

is the SB interaction, which consists of the linear–linear (LL) and
square–linear (SL) SB interactions, and $V_{s}({q_{s}})\equiv
V^{(s)}_{\mathrm{LL}}{q_{s}}+V^{(s)}_{\mathrm{SL}}{q_{s}}^{2}/2$, with
coupling strengths of $V^{(s)}_{\mathrm{LL}}$, $V^{(s)}_{\mathrm{SL}}$, and
$\alpha_{j_{s}}$.[56, 57, 58, 44, 45, 43] The last term of the bath
Hamiltonian is a counter-term that maintains the translational symmetry of the
system.[61] While the LL interaction contributes mainly to energy relaxation,
the SL system–bath interaction leads to vibrational dephasing in the slow
modulation case due to the frequency fluctuation of system vibrations.[44, 45,
43] The sum of the bath coordinates,
${X}_{s}\equiv\sum_{j_{s}}\alpha_{{j_{s}}}{x}_{{j_{s}}}$, acts as a collective
coordinate that modulates mode $s$.[60, 61, 62] We introduce the spectral
distribution function (SDF)
$J_{s}(\omega)\equiv\sum_{j_{s}}\alpha_{j_{s}}^{2}\hbar\delta(\omega-\omega_{j_{s}})/2m_{j_{s}}\omega_{j_{s}}$
to characterize the environmental fluctuation. We assume that $J_{s}(\omega)$
has an Ohmic form with a Lorentzian cutoff:

$\displaystyle J_{s}(\omega)=\frac{\hbar
m_{s}\zeta_{s}}{2\pi}\frac{\omega\gamma_{s}^{2}}{\gamma_{s}^{2}+\omega^{2}}.$
(5)

Here, $\zeta_{s}$ is the system–bath coupling strength, and $\gamma_{s}$
represents the width of the spectral distribution for mode $s$ that relates to
the vibrational dephasing time, defined as $\tau_{s}=1/\gamma_{s}$.

Table 2: Parameter values of multimode LL+SL BO model for anharmonic mode–mode coupling and optical properties among (1) stretching, (2) bending, (3) librational, and (4) translational modes on the basis of Refs. 36, 37. The normalized parameters are defined as $\tilde{g}_{s^{2}s^{\prime}}\equiv(\omega_{0}^{3}/\omega_{s}^{2}\omega_{s}^{\prime})g_{s^{2}s^{\prime}}$, $\tilde{g}_{s{s^{\prime}}^{2}}\equiv(\omega_{0}^{3}/\omega_{s}\omega_{s^{\prime}}^{2})g_{s{s^{\prime}}^{2}}$, $\tilde{\mu}_{ss^{\prime}}\equiv(\omega_{0}^{2}/\omega_{s}\omega_{s}^{\prime})\mu_{ss^{\prime}}$, and $\tilde{\Pi}_{ss^{\prime}}\equiv(\omega_{0}^{2}/\omega_{s}\omega_{s}^{\prime})\Pi_{ss^{\prime}}$.  $\mathrm{s-s^{\prime}}$ | $\tilde{g}_{s^{2}s^{\prime}}$ | $\tilde{g}_{s{s^{\prime}}^{2}}$ | $\tilde{\mu}_{ss^{\prime}}$ | $\tilde{\Pi}_{ss^{\prime}}$
---|---|---|---|---
$\mathrm{1-2}$ | $0$ | $0.2$ | $2.0\times 10^{-3}$ | $2.6\times 10^{-3}$
$\mathrm{1-3}$ | $-3.9\times 10^{-2}$ | $-3.9\times 10^{-2}$ | $0.13$ | $0.19$
$\mathrm{1-4}$ | $-7.5\times 10^{-2}$ | $-7.5\times 10^{-2}$ | $0.43$ | $0.46$
$\mathrm{2-3}$ | $-1.5\times 10^{-2}$ | $-1.5\times 10^{-2}$ | $7.0$ | $4.0$
$\mathrm{2-4}$ | $-2.0\times 10^{-2}$ | $-2.0\times 10^{-2}$ | $3.1\times 10^{-2}$ | $3.1\times 10^{-2}$
$\mathrm{3-4}$ | $0.23$ | $0.23$ | $7.8\times 10^{-2}$ | $0.16$

The dipole operator and polarizability for the $s$th mode are defined as

$\hat{\mu}_{s}=\mu_{s}{q}_{s}+\sum_{s^{\prime}}\mu_{ss^{\prime}}{q}_{s}{q}_{s^{\prime}}$
(6)

and

$\hat{\Pi}_{s}=\Pi_{s}{q}_{s}+\sum_{s^{\prime}}\Pi_{ss^{\prime}}{q}_{s}{q}_{s^{\prime}},$
(7)

where $\mu_{s}$ and $\mu_{ss^{\prime}}$ are the linear and nonlinear elements
of the dipole moment, and $\Pi_{s}$ and $\Pi_{ss^{\prime}}$ are those of the
polarizability, respectively. The vibrational modes interact through the
mechanical anharmonic coupling (MAHC) described by ${g}_{s^{2}s^{\prime}}$ and
${g}_{ss^{\prime 2}}$ and the electric anharmonic coupling (EAHC) described by
${\mu}_{ss^{\prime}}$ and ${\Pi}_{ss^{\prime}}$. Although the effects of MAHC
and EAHC in 2D spectra have been discussed analytically,[38, 39, 69, 70, 70,
71, 72] no quantum theory exists for complex systems involving complex
system–bath interactions, such as those in this study.

In the present case, the difference between the fundamental frequencies of the
intermolecular and intramolecular modes is very large, and we can excite each
mode separately by choosing the light sources: a high-frequency intramolecular
mode is solely excited by an IR pulse, whereas a low-frequency intermolecular
mode is solely excited by a Raman or THz pulse. Note that, although the
contribution from EAHC arises from nonlinear dipole elements or nonlinear
polarizability elements are not important for a spectroscopic approach based
on the third-order response functions for $s=$ 1 and 2 and they can be
ignored,[38, 39] we used the same parameter values presented in Ref. 37 as
given in Tables 1 and 2 because the computational cost in our approach is the
same. Besides these, we used the same parameter values presented in Ref. 37.
The element $\mu_{23}$ is important even in 1D IR because this is the cause of
the combination band at the frequency $\omega=2100$ cm-1.

We chose a set of parameters for the model based on the MD results of 2D
spectra and experimentally obtained 1D spectra, taking advantage of the fact
that the 2D profiles differ significantly for different dynamic processes.[36,
37]

##  III Calculating 2D correlation spectra of third-order IR and fifth-order
IR–Raman measurements

Figure 1: Pulse sequence of third-order 2D IR and fifth-order 2D
IR–IR–Raman–Raman spectroscopies.

We consider 2D IR and IR–Raman spectroscopies whose pulse sequences are
depicted in Fig. 1. These spectra can be calculated from the three-body
nonlinear response function defined as[61, 62]

$R(t_{3},t_{2},t_{1})=\quantity(\frac{i}{\hbar})^{3}\mathrm{tr}\quantity{\hat{A}\mathcal{G}(t_{3})\hat{B}^{\times}\mathcal{G}(t_{2})\hat{C}^{\times}\mathcal{G}(t_{1})\hat{D}^{\times}\hat{\rho}^{\mathrm{eq}}},$
(8)

where the operators $\hat{A}$, $\hat{B}$, $\hat{C}$, and $\hat{D}$ are either
the dipole moment $\hat{\mu}$ or the polarizability $\hat{\Pi}$. Here, we have
employed the hyperoperator × defined as
$\hat{A}^{\times}\hat{{\hat{\rho}}}\equiv[\hat{{A}},{\hat{\rho}}]$,
$\mathcal{G}(t)$ is the Green’s function of the total Hamiltonian without a
laser interaction, and $\hat{\rho}_{\mathrm{eq}}$ is the equilibrium state.
Using this expression, the third-order IR signal is calculated from
$\hat{A}$$\sim$$\hat{D}=\hat{\mu}$ and the fifth-order IR–Raman signal is
calculated from $\hat{A}=\hat{B}=\hat{\Pi}$ and $\hat{C}=\hat{D}=\hat{\mu}$,
respectively.[57, 58, 59] Note that, in our model calculation, the role of IR
and Raman excitation for the intramolecular modes and that of THz and Raman
excitation for the intermolecular modes are similar because the nonlinear
elements of each mode, such as $\mu_{ss}$ and $\Pi_{ss}$, do not play a major
role in determining the profiles of the 2D spectrum, which is defined as Eq.
(8). This is because nonlinear spectroscopy based on a three-body response
function such as
$R_{\mathrm{IIR}}^{(3)}(t_{2},t_{1})=-\langle[[\hat{\bm{\Pi}}(t_{2}+t_{1}),\hat{\bm{\mu}}(t_{1})],\hat{\bm{\mu}}(0)]\rangle/\hbar^{2}$
is described as an odd number of optical operators, whereas that based on a
four-body response function such as
$R_{\mathrm{IIRR}}^{(5)}(t_{3},t_{2},t_{1})=-i\langle[[\hat{\bm{\Pi}}(t_{3}+t_{2}+t_{1}),[[\hat{\bm{\Pi}}(t_{2}+t_{1}),\hat{\bm{\mu}}(t_{1})],\hat{\bm{\mu}}(0)]]\rangle/\hbar^{3}$
is described as an even number of optical operators. When the anharmonicity of
the mode described by $g_{s^{3}}q_{s}^{3}$ is weak, the thermal statistical
average included in the calculation of the response function is performed as a
Gaussian integral of $q_{s}$. Therefore, $R_{\mathrm{IIR}}^{(3)}(t_{2},t_{1})$
must be of even order function of $q_{s}$ so that it does not disappear in the
Gaussian integral, and a contribution of ${\Pi}_{ss}$ or ${\mu}_{ss}$ is
required. In contrast, $R_{\mathrm{IIRR}}^{(5)}(t_{3},t_{2},t_{1})$ is fourth
order in $q_{s}$ even without ${\Pi}_{ss}$ or ${\Pi}_{ss}$, and, if they
contribute, it should be of the second order in ${\Pi}_{ss}$ and/or
${\mu}_{ss}$, such as ${\Pi}_{ss}{\mu}_{ss}$.[21] This indicates that the MAHC
contribution, which is often referred as the non-Condon effect,[73] in
$R_{\mathrm{IIR}}^{(3)}(t_{2},t_{1})$ and
$R_{\mathrm{IIRR}}^{(5)}(t_{3},t_{2},t_{1})$ measurements is minor. Thus,
although it would be difficult to conduct an actual measurement, the profile
of a seventh-order 2D Raman spectrum (i.e.,
$\hat{A}$$\sim$$\hat{D}=\hat{\Pi}$)[57, 58, 74] should be similar to that of
third-order 2D IR.

The calculated signals are plotted as the 2D correlation spectrum that is
obtained by adding the two terms corresponding to the rephasing and
nonrephasing parts of $R(t_{3},t_{2},t_{1})$ in equal weights. A commonly used
definition of a 2D correlation spectrum is expressed as[75, 76, 77]

$\displaystyle I_{\rm{C}}({\omega_{3}},{t_{2}},{\omega_{1}})$
$\displaystyle={\rm
Im}\left\\{\int_{0}^{\infty}\differential{t_{1}}\int_{0}^{\infty}\differential{t_{3}}e^{-i{\omega_{1}}{t_{1}}}e^{+i{\omega_{3}}{t_{3}}}R_{\rm{R}}({t_{3}},{t_{2}},{t_{1}})\right\\}$
$\displaystyle+{\rm
Im}\left\\{\int_{0}^{\infty}\differential{t_{1}}\int_{0}^{\infty}\differential{t_{3}}{e^{i{\omega_{1}}{t_{1}}}}{e^{+i{\omega_{3}}{t_{3}}}}R_{{\rm{NR}}}({t_{3}},{t_{2}},{t_{1}})\right\\},$

where $R_{{\rm{R}}}$ and $R_{{\rm{NR}}}$ are the rephasing and nonrephasing
responses that are experimentally detected separately using the phase-matching
conditions. Theoretically, although we can separate $R_{{\rm{R}}}$ and
$R_{{\rm{NR}}}$ by choosing the specific Liouville paths in the energy state
model, the process is not easy when we calculate the response functions in the
coordinate space. Therefore, here we eliminate the undesired rephasing
contribution by utilizing the Fourier transformation of $t_{2}$ for
$I_{\rm{C}}({\omega_{3}},{t_{2}},{\omega_{1}})$ with
$R_{\rm{R}}({t_{3}},{t_{2}},{t_{1}})=R_{{\rm{NR}}}({t_{3}},{t_{2}},{t_{1}})=R(t_{3},t_{2},t_{1})$
to remove the oscillating contribution in period $t_{2}$ period with frequency
$2\omega_{s}$, where $\omega_{s}$ is the frequency of the targeting mode
$s$.[78, 79] Note that we can also separate the plain response function in
terms of different phase-matching conditions using a phase-matching
transformation,[80] but the Fourier based method is simpler when only a 2D
correlation spectrum is calculated.

After removing the dispersive contribution of the response function, we obtain
the pure absorptive line shape. We observe a distorted line shape, called the
phase-twisted line, which exhibits the non-Markovian and nonperturbative
dynamics of the environment. [22, 23, 24, 75, 76, 77] To investigate not only
anharmonic mode–mode coupling but also vibrational dephasing and thermal
relaxation effects at finite temperature, we must employ nonperturbative and
non-Markovian open quantum dynamics theory, which enables us to treat the
quantum entanglement between the system and bath.[62] Therefore, we employ the
discretized hierarchical equations of motion in mixed Liouville–Wigner space
(DHEOM–MLWS): Lagrange–Hermite mesh discretization is employed in the
Liouville space of the intramolecular modes, and Lagrange–Hermite mesh
discretization and Hermite discretization are employed in the Wigner space of
the intermolecular modes.[37]

##  IV RESULTS AND DISCUSSION

Although the DHEOM–MLWS for the water BO model was created to treat the above
four modes simultaneously, the computational time required for 2D vibrational
spectra is about 1000 times longer than that for linear spectra, even with a
fixed $t_{2}$. Therefore, here we calculate 2D spectra for mode–mode coupling,
which can be investigated for two modes, as in combinations of
stretching–bending (1–2) and bending–librational (2–3) for 2D correlation IR
spectra, and stretching–librational (1–3) and stretching–translational (1–4)
for 2D correlation IR–Raman spectra. To calculate the spectrum, we integrate
the DHEOM–MWLS using the same numerical setup as in Ref. 37 with the parameter
set presented in Tables 1 and 2. Thus, the calculated 1D IR and 1D Raman
spectra (not shown) are identical to those presented in Ref. 37.

Note that the dipole moment and polarizability have similar functional forms
in this model calculation. Therefore, if the same modes are excited, 2D
Raman–Raman–IR–IR or 2D IR–Raman–IR–Raman exhibit a profile that is similar to
the present 2D IR–IR–Raman–Raman results.

###  IV.1 2D correlation IR spectra

####  IV.1.1 Stretching–bending modes

Figure 2: Third-order 2D correlation IR spectra for stretching motion (upper
panel) and stretching $\rightarrow$ bending motion (lower panel) for different
$t_{2}$. All spectral intensities are normalized with respect to the absolute
value of the maximum peak intensity of each diagonal peak. The direction of
the nodal lines (orange dashed lines) in the upper panel represents the extent
of correlation between the vibrational coherences of the $t_{1}$ and $t_{3}$
periods. For clarity, the data of off-diagonal peaks are multiplied by 3 in
the cases of $t_{2}=0$ fs and $50$ fs.

Figure 2 illustrates 2D correlation spectra calculated for the
stretching–bending (1–2) modes for which experimental results are
available.[26, 27, 28] Note that previous theoretical investigations for water
2D spectroscopy have been limited to stretching modes only for quantum
calculations based on classical MD[48, 49, 50] and individual stretching and
bending modes for classical MD.[47, 46] The present investigations provide the
first quantum simulation of 2D spectra for intramolecular–intramolecular and
intramolecular–intermolecular mode–mode couplings. To conduct the simulations,
we set $\hat{\bm{\mu}}=\hat{{\mu}}_{1}+\hat{{\mu}}_{2}$. Each peak near
$(\omega_{1},\omega_{3})=(3400\rm{cm}^{-1},3400cm^{-1})$ in the upper panel
represents the stretching mode; the red positive peaks and the blue negative
peaks correspond to the transitions between the
$|0\rangle_{s}$$\rightarrow_{s}$$|1\rangle_{s}$$\rightarrow$$|0\rangle_{s}$
and $|0\rangle_{s}$$\rightarrow$$|1\rangle_{s}$$\rightarrow$$|2\rangle_{s}$
states for $s=1$, respectively, where $|n\rangle_{s}$ is the $n$th vibrational
energy state of the $s$th mode.

In the lower panel of each figure in Fig. 2 present the cross peaks
representing the stretching$\rightarrow$bending transition (e.g.,
$|0\rangle_{1}|0\rangle_{2}$$\rightarrow$$|1\rangle_{1}|0\rangle_{2}$$\sim$$|0\rangle_{1}|1\rangle_{2}$$\rightarrow$$|0\rangle_{1}|0\rangle_{2}$).[22]
The coupling peaks that appear at $t_{2}=0$ are due to coherent transfer
between the two modes and quickly disappear at $t_{2}=50$fs, when coherence is
lost. The coupling peaks that appear after $t_{2}\geq$100 fs are due to
population transfer, and their intensity increases with $t_{2}$.

The direction of the orange nodal lines represents the extent of noise
correlation (non-Markovian effects) between the vibrational coherences of
periods $t_{1}$ and $t_{3}$. The direction parallel to the $\omega_{1}$ axis
is the uncorrelated case, whereas the direction parallel to the
$\omega_{1}=\omega_{3}$ line is the fully correlated case.[75, 76, 77] The
width of the peak parallel to the $\omega_{1}=\omega_{3}$ direction
corresponds to inhomogeneous broadening, whereas that perpendicular to
$\omega_{1}=\omega_{3}$ is homogeneous broadening.[81, 22]

The spectral width of the stretching peaks in the upper panel of Fig. 2 is
broadened due to the dephasing (spectral diffusion regime) that arises from
the non-perturbative and non-Markovian effects of the SL interaction.[44] When
$t_{2}$ becomes larger than the correlation time of the noise
$\tau_{1}=1/\gamma_{1}\approx 160$ fs, the correlation begins to weaken, and
at $t_{2}=500$ fs, the nodal line becomes horizontal. The SL interaction used
in this model causes the same effect as in stochastic theory, where $\tau_{1}$
determines the dephasing time.[61] We found that our result is consistent with
those of previous experiments [82, 25, 26, 27, 28] and theoretical
calculations.[48, 49, 50, 47, 46]

Upon increasing $t_{2}$, it was observed experimentally that the blue 0–1–2
peak extended in the small $\omega_{3}$ direction for $t_{2}\leq 250$fs, and a
peak, referred to as the hot ground state, appeared above the red 0–1–0 peak
for $t_{2}\geq 250$fs. [25, 26, 27, 28] In contrast, our results show that the
blue peak extends only slightly in the negative $\omega_{3}$ for $t_{2}\geq
500$fs direction due to the population relaxation arising from the SL
interaction,[44] while no hot ground state peak is observed. Note that the
inclusion of the LL interactions in the stretching[44, 45] and/or the presence
of non-radiative decay states[68] would enhance the elongation. To predict
such phenomena with a current handy model calculation, one must include
various hypothetical factors that are considered important in determining the
vibrational process of water, with reference to the experimentally obtained 2D
spectral profile. We have computed 2D spectra for the case of including
quartic Darling–Dennison (DD) coupling described as
$\tilde{g}_{1^{2}2^{2}}q_{1}^{2}g_{2}^{2}$.[72, 22] and for the case of
stronger cubic (Fermi) coupling ($\tilde{g}_{s{s^{\prime}}^{2}}=0.5$). While
inclusion of the DD coupling does not change the 2D profile much (not shown),
we found that, when the Fermi coupling is strong, the 2D profile is relatively
close to the experimental results (see Appendix A). However, some other
information, such as quantum MD results, is needed to explore such a
possibility in greater detail.

The cross peaks representing the stretching$\rightarrow$bending transition in
the lower panel of Fig. 2 are uncorrelated due to the process involved in the
bending mode having a very short noise correlation time, $1/\gamma_{2}\approx
40$fs, and the positive and negative peaks are elongated in the $\omega_{1}$
direction. [26, 27, 28] In Appendix B, we depict the 2D spectrum of the
bending mode. It can be seen that the time scale of the spectral diffusion of
the cross-peak is determined by the mode in which it is faster.

####  IV.1.2 Bending–librational modes

Figure 3: Third-order 2D correlation IR spectra for the
bending$\rightarrow$librational (2$\rightarrow$3) transition for different
$t_{2}$. All spectral intensities are normalized with respect to the absolute
value of the maximum peak intensity in the $t_{2}=0$ case.

Next, we present the calculated result for the bending$\rightarrow$librational
(2$\rightarrow$3) transition peak in Fig. 3. We chose the first two IR pulses
for the excitation of the bending mode and last two IR (or THz) pulses for the
detection of the librational modes. In the 2D correlation spectrum based on
the third-order (four body) response function, we can interpret the 2D
spectral profile as investigating the excitation of each vibrational mode in
the $\omega_{1}$ direction in the IR region and then detecting the dynamics of
that excited motion after $t_{2}$ as the absorption spectrum in the
$\omega_{3}$ direction in the THz region. In Appendix B, we depict the 2D
spectrum of the bending mode for the same conditions as in Fig. 3.

In the present case, the first two IR pulses excite the bending mode
($\omega_{1}=1600$ cm-1) and the combination band (around $\omega_{1}=2100$
cm-1) through $\mu_{2}$ and $\mu_{23}$, respectively. Because no population
transfer for 2$\rightarrow$3 occurs at time $t_{2}=0$, we cannot observe a
prominent peak in the $\omega_{3}$ direction along the line $\omega_{1}=1600$
cm-1. The positive and negative peaks of the combination band peak around
$\omega_{1}=2100$ cm-1 appear because the IR excitations described by
$\mu_{23}$ cause transitions of the librational vibrational states,
$|0\rangle_{2}|0\rangle_{3}$$\rightarrow$$|1\rangle_{2}|1\rangle_{3}$$\rightarrow$$|1\rangle_{2}|0\rangle_{3}$
and
$|0\rangle_{2}|0\rangle_{3}$$\rightarrow$$|1\rangle_{2}|1\rangle_{3}$$\rightarrow$$|1\rangle_{2}|2\rangle_{3}$,
which appear as the positive upper and negative lower peaks, respectively.

At $t_{2}=20$ fs, the elongated peak along line $\omega_{1}=1600$ cm-1 appears
as the result of the population transfer from 2$\rightarrow$3\. The elongation
of the peak along this line reflects the fact that the librational mode is
spread out. Because of this, the positive and negative combination band peaks
are also elongated in the $\omega_{3}$ direction. The nodal line of the
combination band peaks in the high-frequency region is twisted, indicating
that the contribution from the librational mode, which seems to dominate the
high-frequency part of the combination band, can be correlated with the IR
excitation of the librational mode. Because these combination band peaks
appear due to the coherent process
$|0\rangle_{2}$$\rightarrow$$|1\rangle_{2}$, they decay quickly and almost
disappear at $t_{2}=80$ fs, whereas the 2$\rightarrow$3 transition peak along
line $\omega_{1}=1600$cm-1 grows because this process arises from the
relaxation. The fast relaxation found between these two modes is also
consistent with other classical MD analyses.[55]

###  IV.2 Fifth-order 2D correlation IR–IR–Raman–Raman spectra

Next, we display our calculated results for the mode–mode coupling of the
high-frequency stretching mode to the low-frequency intramolecular modes.
Here, we chose IR pulses for the excitation of the stretching mode and Raman
pulses for the detection of the intermolecular modes because we believe that
such an experimental setup is feasible with current technology. However,
because the difference in the mathematical treatment of Raman and IR (or THz)
excitation is minor in our model calculation, the spectral profile would be
similar regardless of the laser configuration, for example, for a
configuration where Raman pulses excite the stretching mode and THz pulses
detect the intermolecular modes.

####  IV.2.1 Stretching–translational modes

Figure 4: Fifth-order 2D correlation IR–IR–Raman–Raman spectra for the
stretching$\rightarrow$translational (1$\rightarrow$4) coupling for different
$t_{2}$. All spectral intensities are normalized with respect to the absolute
value of the maximum peak intensity in the $t_{2}=0$ case.

We now present the simulation result of fifth-order 2D IR–IR–Raman–Raman
spectra for the stretching–translational coupling case. Note that this process
was investigated theoretically and experimentally through the use of 2D
THz–IR–visible (or 2D IR–IR–Raman) measurement.[36, 34, 35, 37] The cross
peaks predicted from the 2D THz–IR–visible spectra arise from both MAHC and
EAHC contributions of the stretching–translational modes, but, in the 2D
IR–IR–Raman–Raman spectra, the cross peaks always arise from MAHC, whereas
those from EHAC appear at different positions or do not appear, as can be seen
from the analytical expressions of the response functions.[69, 70, 71, 38, 39]

Moreover, although 2D THz–IR–visible measurement is designed to detect the
coherence between the two different modes through the mode–mode interactions,
the present 2D IR–IR–Raman–Raman spectroscopy is designed to detect the
population transfer among modes under dephasing and relaxation through the use
of period $t_{2}$. Thus, the information obtained from the two approaches is
different and complimentary.

The results of 2D IR–IR–Raman–Raman (or 2D correlation IR–IR–Raman–Raman)
spectra for different $t_{2}$ values are presented in Fig. 4. As in the
bending-librational case in Fig. 3, the first two IR pulses applied at time
duration $t_{1}$ excite not only the fast stretching mode but also the slow
translational mode due to the presence of the large $\mu_{14}$ element. As a
result, the phase of the stretching motion (with large $\omega_{1}$) is
twisted by the slow translational motion. Thus, the intensities of the
stretching–translational peaks vary as a function of $t_{2}$. Because the
vibrational motion of the translational mode decays quickly due to the large
LL interaction, the phase twist disappears around $t_{2}=160$ fs, and a single
positive peak is observed along line $\omega_{1}=3400$ cm-1 at $t_{2}=300$ fs.
This peak is elongated in the $\omega_{3}$ direction, reflecting the fact that
the translational mode extends into the high-frequency region.

####  IV.2.2 Stretching–librational modes

Figure 5: 2D correlation IR–IR–Raman–Raman spectra for
stretching$\rightarrow$librational ($1\rightarrow$3) coupling for different
$t_{2}$ values. All spectral intensities are normalized with respect to the
absolute value of the maximum peak intensity in the $t_{2}=0$ case.

Finally, 2D correlation IR–IR–Raman–Raman spectra for the
stretching–librational case are plotted in Fig. 5. As shown in Table 2, the
intensity of the $\mu_{13}$ element is comparable to $\mu_{1}$. Thus, the
first two infrared pulses excite not only the stretching mode but also the
librational mode, and, as in the stretching–translational case, the spectral
profile changes in time rapidly for small $t_{2}$ values due to the phase
twist. After librational oscillation decays around $t_{2}=100$ fs, we observe
a single positive peak at the stretching–librational coupling position, which
extends in the $\omega_{3}$ direction due to the large broadening of the
librational peak.

##  V CONCLUSION

We calculated third-order 2D correlation IR and fifth-order 2D correlation
IR–IR–Raman–Raman spectra for the stretching–bending (1-2),
bending–librational (2-3), stretching–librational (1-3), and
stretching–translational (1-4) modes of liquid water using the DHEOM–MLWS
approach based on the multimode LL BO model. Because the peaks that arise from
EAHC play a minor role in four-body response function–based spectroscopies, we
can easily discuss the role of mode–mode couplings that arise from MAHC
processes. As the peaks arising from MAHC in this order of spectroscopy are
associated with the population transfer process, they appear after the
population relaxation occurs, whereas the peaks arising from EAHC are
associated with the coherent transfer and appear at $t_{2}=0$ before decaying
rapidly.

To reproduce a 2D spectral profile, the model must capture the dynamical
properties of the vibrational motions accurately. Thus, the scientific
scenario underlying the model can be validated by calculating 2D spectra and
comparing them to those obtained from experiments or accurate simulations.
Taking advantage of the low computational cost and simplicity of the model, we
can easily examine the mechanism to determine the 2D profile of spectrum, as
briefly demonstrated in Appendix A. Although computationally expensive, within
the framework of HEOM formalism, it is possible to describe the OH–stretching
modes separately as symmetric and antisymmetric modes[68] or to introduce
another optically inactive mode that leads the excitation of the stretching
mode to the vibrational ground state.[83] However, the choice of models and
model parameters should be justified by comparing results obtained by advanced
experimental and simulation techniques.

With the limitation of CPU power, here we considered only the combination of
the two-mode cases. To investigate a pathway of energy relaxation from the
OH–Stretching mode to the low-frequency intramolecular modes, we should
consider at least three modes, such as the stretching $\rightarrow$ bending
$\rightarrow$ librational (1 $\rightarrow$ 2 $\rightarrow$ 3) modes or the
stretching $\rightarrow$ bending $\rightarrow$ translational (1 $\rightarrow$
2 $\rightarrow$ 4) modes including the anharmonicity of the three modes
described as
$\tilde{g}_{s{s^{\prime}}{s^{\prime\prime}}}q_{s}q_{s^{\prime}}q_{s^{\prime\prime}}$.
We leave such investigations to future studies.

Figure 6: Third-order 2D correlation IR spectra for stretching motion (upper
panel) and stretching $\rightarrow$ bending motion (lower panel) calculated
for a stronger cubic coupling ($\tilde{g}_{s{s^{\prime}}^{2}}=0.5$) than the
case in Fig. 2. All spectral intensities are normalized with respect to the
absolute value of the maximum peak intensity of each diagonal peak. The
direction of the nodal lines (orange dashed lines) in the upper panel
represents the extent of correlation between the vibrational coherences of the
$t_{1}$ and $t_{3}$ periods. The data of off-diagonal peaks are multiplied by
3 for clarity in the cases of $t_{2}=0$ fs and $50$ fs.  Figure 7: Third-order
2D correlation IR spectra of the bending mode ($s=2$) calculated with the
stretching–bending coupling for different $t_{2}$. All spectral intensities
are normalized with respect to the absolute value of the maximum peak
intensity in the $t_{2}=0$ case.

## Acknowledgments

The authors are thankful to Professor Shinji Saito and Professor Keisuke
Tominaga for helpful discussions. Y.T. was supported by JSPS KAKENHI (Grant
No. B21H01884). H.T. was supported by JST SPRING (Grant No. JPMJSP2110).

### Conflict of Interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A 2D IR spectra of stretching–bending modes in the enhanced mode-
mode coupling case

In the investigation of the stretching mode, it was observed experimentally
that as $t_{2}$ increased, the negative (blue) 0–1–2 peak extended in the
small $\omega_{3}$ direction, and a peak, referred to as the hot ground state,
appeared above the positive (red) 0–1–0 peak. [25, 26, 27, 28] To predict such
phenomena with the current approach, one must include various hypothetical
factors that are considered important in determining the vibrational process
of water, with reference to the 2D spectral profile. Here, we examine the
effects of the cubic coupling strength. Thus we made $\tilde{g}_{12^{2}}$ 2.5
times larger than the case in Fig. 2, leaving the other parameter values
unchanged. Figure 6 illustrates the calculated results. Although less
pronounced than the experimentally observed spectrum, the blue 0–1–2 peak
extends toward small $\omega_{3}$ direction from $t_{2}=0$, which indicates
that this peak is instantaneously generated by the coherent transfer rather
than a population transfer between the two modes. A small negative peak
appears for $t_{2}\geq 100$fs above the red 0–1–0 peak, although its peak
position is much higher than the experimentally observed hot grand state peak.
The appearance of this peak indicates that this peak arises from a population
transfer rather than coherent transfer as the stretching$\rightarrow$bending
peak for $t_{2}\geq 100$fs in the lower panels. The vibrational coherence,
represented by the angle of the orange nodal line, is slightly weaker than in
the case in Fig. 2.

##  Appendix B 2D IR spectra of bending mode

In this Appendix, we present third-order 2D IR spectra of the bending mode
calculated for the bending$\rightarrow$librational (2$\rightarrow$3)
transition case under the same conditions as in Fig. 3. This spectrum has been
measured experimentally.[29] In this case, we observe horizontally elongated
positive (red) and negative(blue) peaks along $\omega_{3}=1650$ cm-1. The
nodal line is always parallel to the $\omega_{1}$ axis because the noise
correlation time is very short ($1/\gamma_{2}\approx 40$ fs). Note that we
have also computed 2D spectra for the case of single bending mode (not shown)
and we found that the inclusion of the bending$\rightarrow$librational
(2$\rightarrow$3) coupling does not change the 2D profile much. This indicates
the effect of population and coherent transfer for the
bending$\rightarrow$librational transition is minor. As a result, the peak
profile changes only slightly as $t_{2}$ increases.

## references

## References

  * Ohmine and Saito [1999] I. Ohmine and S. Saito, “Water dynamics : fluctuation, relaxation, and chemical reactions in hydrogen bond network rearrangement,” Accounts of Chemical Research 32, 741–749 (1999).
  * Nibbering and Elsaesser [2004] E. T. J. Nibbering and T. Elsaesser, “Ultrafast vibrational dynamics of hydrogen bonds in the condensed phase.” Chemical reviews 104 4, 1887–1914 (2004).
  * Steinel  _et al._ [2004] T. Steinel, J. B. Asbury, S. Corcelli, C. Lawrence, J. Skinner, and M. Fayer, “Water dynamics: dependence on local structure probed with vibrational echo correlation spectroscopy,” Chemical Physics Letters 386, 295–300 (2004).
  * Bakker and Skinner [2010] H. J. Bakker and J. L. Skinner, “Vibrational spectroscopy as a probe of structure and dynamics in liquid water,” Chemical Reviews 110, 1498–1517 (2010).
  * Fecko  _et al._ [2003] C. Fecko, J. Eaves, J. Loparo, A. Tokmakoff, and P. Geissler, “Ultrafast hydrogen-bond dynamics in the infrared spectroscopy of water,” SCIENCE 301, 1698–1702 (2003).
  * Rey, Møller, and Hynes [2004] R. Rey, K. B. Møller, and J. T. Hynes, “Ultrafast vibrational population dynamics of water and related systems, a theoretical perspective,” Chemical Reviews 104, 1915–1928 (2004).
  * Liu  _et al._ [2011] J. Liu, W. H. Miller, G. S. Fanourgakis, S. S. Xantheas, S. Imoto, and S. Saito, “Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initio-based flexible and polarizable force field,” The Journal of Chemical Physics 135, 244503 (2011), https://doi.org/10.1063/1.3670960 .
  * Perakis  _et al._ [2016] F. Perakis, L. De Marco, A. Shalit, F. Tang, Z. R. Kann, T. D. Kühne, R. Torre, M. Bonn, and Y. Nagata, “Vibrational spectroscopy and dynamics of water,” Chemical Reviews 116, 7590–7607 (2016), https://doi.org/10.1021/acs.chemrev.5b00640 .
  * Ohmine and Tanaka [1993] I. Ohmine and H. Tanaka, “Fluctuation, relaxations, and hydration in liquid water. hydrogen-bond rearrangement dynamics,” Chemical Reviews 93, 2545–2566 (1993), https://doi.org/10.1021/cr00023a011 .
  * Yagasaki and Saito [2009] T. Yagasaki and S. Saito, “Molecular dynamics simulation of nonlinear spectroscopies of intermolecular motions in liquid water,” Accounts of Chemical Research 42, 1250–1258 (2009), https://doi.org/10.1021/ar900007s .
  * Yagasaki and Saito [2013] T. Yagasaki and S. Saito, “Fluctuations and relaxation dynamics of liquid water revealed by linear and nonlinear spectroscopy,” Annual Review of Physical Chemistry 64, 55–75 (2013), https://doi.org/10.1146/annurev-physchem-040412-110150 .
  * Imoto, Xantheas, and Saito [2015] S. Imoto, S. S. Xantheas, and S. Saito, “Ultrafast dynamics of liquid water: Energy relaxation and transfer processes of the OH stretch and the HOH bend,” The Journal of Physical Chemistry B 119, 11068–11078 (2015), pMID: 26042611, https://doi.org/10.1021/acs.jpcb.5b02589 .
  * Huse  _et al._ [2005] N. Huse, S. Ashihara, E. Nibbering, and T. Elsaesser, “Ultrafast vibrational relaxation of O-H bending and librational excitations in liquid H2O,” CHEMICAL PHYSICS LETTERS 404, 389–393 (2005).
  * Ramasesha  _et al._ [2013] K. Ramasesha, L. De Marco, A. Mandal, and A. Tokmakoff, “Water vibrations have strongly mixed intra- and intermolecular character,” Nature Chemistry 5, 935–940 (2013).
  * Bertie and Lan [1996] J. E. Bertie and Z. Lan, “Infrared intensities of liquids xx: The intensity of the OH stretching band of liquid water revisited, and the best current values of the optical constants of H2O(l) at 25°c between 15,000 and 1 cm-1,” Appl. Spectrosc. 50, 1047–1057 (1996).
  * Maréchal [2011] Y. Maréchal, “The molecular structure of liquid water delivered by absorption spectroscopy in the whole IR region completed with thermodynamics data,” Journal of Molecular Structure 1004, 146–155 (2011).
  * Brooker  _et al._ [1989] M. H. Brooker, G. W. Hancock, B. Rice, and J. G. Shapter, “Raman frequency and intensity studies of liquid H2O, H218O and D2O,” Journal of Raman Spectroscopy 20, 683–694 (1989).
  * Pattenaude, Streacker, and Ben-Amotz [2018] S. R. Pattenaude, L. M. Streacker, and D. Ben-Amotz, “Temperature and polarization dependent Raman spectra of liquid H2O and D2O,” Journal of Raman Spectroscopy 49, 1860–1866 (2018).
  * Heisler, Mazur, and Meech [2017] I. A. Heisler, K. Mazur, and S. R. Meech, “Raman vibrational dynamics of hydrated ions in the low-frequency spectral region,” Journal of Molecular Liquids 228, 45–53 (2017), from simple liquids to macromolecular solutions: recent experimental and theoretical developments. In Honor of the 70th birthday of Vojko Vlachy.
  * Mukamel [1999] S. Mukamel,  _Principles of nonlinear optical spectroscopy_ , 6 (Oxford University Press on Demand, 1999).
  * Tanimura and Mukamel [1993] Y. Tanimura and S. Mukamel, “Two-dimensional femtosecond vibrational spectroscopy of liquids,” The Journal of Chemical Physics 99, 9496–9511 (1993).
  * Tanimura and Ishizaki [2009] Y. Tanimura and A. Ishizaki, “Modeling, calculating, and analyzing multidimensional vibrational spectroscopies,” Accounts of Chemical Research 42, 1270–1279 (2009).
  * Cho [2009] M. Cho,  _Two-Dimensional Optical Spectroscopy_ (CRC Press, 2009).
  * Hamm and Zanni [2011] P. Hamm and M. T. Zanni,  _Concepts and Methods of 2D Infrared Spectroscopy_ (Cambridge University Press, 2011).
  * Kraemer  _et al._ [2008] D. Kraemer, M. L. Cowan, A. Paarmann, N. Huse, E. T. J. Nibbering, T. Elsaesser, and R. J. D. Miller, “Temperature dependence of the two-dimensional infrared spectrum of liquid H2O,” Proceedings of the National Academy of Sciences 105, 437–442 (2008).
  * De Marco  _et al._ [2016] L. De Marco, J. A. Fournier, M. Thämer, W. Carpenter, and A. Tokmakoff, “Anharmonic exciton dynamics and energy dissipation in liquid water from two-dimensional infrared spectroscopy,” The Journal of Chemical Physics 145, 094501 (2016), https://aip.scitation.org/doi/pdf/10.1063/1.4961752 .
  * Carpenter  _et al._ [2017] W. B. Carpenter, J. A. Fournier, R. Biswas, G. A. Voth, and A. Tokmakoff, “Delocalization and stretch-bend mixing of the HOH bend in liquid water,” The Journal of Chemical Physics 147, 084503 (2017), https://doi.org/10.1063/1.4987153 .
  * Lewis  _et al._ [2022] N. H. C. Lewis, B. Dereka, Y. Zhang, E. J. Maginn, and A. Tokmakoff, “From networked to isolated: Observing water hydrogen bonds in concentrated electrolytes with two-dimensional infrared spectroscopy,” The Journal of Physical Chemistry B 126, 5305–5319 (2022), pMID: 35829623, https://doi.org/10.1021/acs.jpcb.2c03341 .
  * Chuntonov, Kumar, and Kuroda [2014] L. Chuntonov, R. Kumar, and D. G. Kuroda, “Non-linear infrared spectroscopy of the water bending mode: direct experimental evidence of hydration shell reorganization?” Phys. Chem. Chem. Phys. 16, 13172–13181 (2014).
  * Hamm and Savolainen [2012] P. Hamm and J. Savolainen, “Two-dimensional-Raman-Terahertz spectroscopy of water: Theory,” The Journal of Chemical Physics 136, 094516 (2012), https://doi.org/10.1063/1.3691601 .
  * Savolainen, Ahmed, and Hamm [2013] J. Savolainen, S. Ahmed, and P. Hamm, “Two-dimensional Raman–Terahertz spectroscopy of water,” Proceedings of the National Academy of Sciences 110, 20402–20407 (2013), https://www.pnas.org/doi/pdf/10.1073/pnas.1317459110 .
  * Hamm [2014] P. Hamm, “2D-Raman-THz spectroscopy: A sensitive test of polarizable water models,” The Journal of Chemical Physics 141, 184201 (2014).
  * Hamm and Shalit [2017] P. Hamm and A. Shalit, “Perspective: Echoes in 2d-Raman-THz spectroscopy,” The Journal of Chemical Physics 146, 130901 (2017), https://doi.org/10.1063/1.4979288 .
  * Grechko  _et al._ [2018] M. Grechko, T. Hasegawa, F. D’Angelo, H. Ito, D. Turchinovich, Y. Nagata, and M. Bonn, “Coupling between intra- and intermolecular motions in liquid water revealed by two-dimensional terahertz-infrared-visible spectroscopy,” Nat Commun 9, 885 (2018).
  * Vietze  _et al._ [2021] L. Vietze, E. H. G. Backus, M. Bonn, and M. Grechko, “Distinguishing different excitation pathways in two-dimensional Terahertz-infrared-visible spectroscopy,” The Journal of Chemical Physics 154, 174201 (2021), https://doi.org/10.1063/5.0047918 .
  * Ito and Tanimura [2016] H. Ito and Y. Tanimura, “Simulating two-dimensional infrared-Raman and Raman spectroscopies for intermolecular and intramolecular modes of liquid water,” The Journal of Chemical Physics 144, 074201 (2016).
  * Takahashi and Tanimura [2023] H. Takahashi and Y. Tanimura, “Discretized hierarchical equations of motion in mixed liouville–wigner space for two-dimensional vibrational spectroscopies of liquid water,” The Journal of Chemical Physics 158, 044115 (2023), https://doi.org/10.1063/5.0135725 .
  * Okumura and Tanimura [2003] K. Okumura and Y. Tanimura, “Energy-level diagrams and their contribution to fifth-order Raman and second-order infrared responses: distinction between relaxation models by two-dimensional spectroscopy,” The Journal of Physical Chemistry A 107, 8092–8105 (2003).
  * Sidler and Hamm [2020] D. Sidler and P. Hamm, “A feynman diagram description of the 2d-raman-thz response of amorphous ice,” The Journal of Chemical Physics 153, 044502 (2020), https://doi.org/10.1063/5.0018485 .
  * Ito, Hasegawa, and Tanimura [2014] H. Ito, T. Hasegawa, and Y. Tanimura, “Calculating two-dimensional THz-Raman-THz and Raman-THz-THz signals for various molecular liquids: The samplers,” The Journal of Chemical Physics 141, 124503 (2014).
  * Ito, Jo, and Tanimura [2015] H. Ito, J.-Y. Jo, and Y. Tanimura, “Notes on simulating two-dimensional Raman and terahertz-Raman signals with a full molecular dynamics simulation approach,” Structural Dynamics 2, 054102 (2015).
  * Ito, Hasegawa, and Tanimura [2016] H. Ito, T. Hasegawa, and Y. Tanimura, “Effects of intermolecular charge transfer in liquid water on Raman spectra,” The Journal of Physical Chemistry Letters 7, 4147–4151 (2016).
  * Ikeda, Ito, and Tanimura [2015] T. Ikeda, H. Ito, and Y. Tanimura, “Analysis of 2d THz-Raman spectroscopy using a non-Markovian Brownian oscillator model with nonlinear system-bath interactions,” The Journal of Chemical Physics 142, 212421 (2015).
  * Ishizaki and Tanimura [2006] A. Ishizaki and Y. Tanimura, “Modeling vibrational dephasing and energy relaxation of intramolecular anharmonic modes for multidimensional infrared spectroscopies,” The Journal of Chemical Physics 125, 084501 (2006).
  * Sakurai and Tanimura [2011] A. Sakurai and Y. Tanimura, “Does $\hbar$ play a role in multidimensional spectroscopy? reduced hierarchy equations of motion approach to molecular vibrations,” The Journal of Physical Chemistry A 115, 4009–4022 (2011).
  * Jeon and Cho [2014] J. Jeon and M. Cho, “An accurate classical simulation of a two-dimensional vibrational spectrum: OD stretch spectrum of a hydrated HOD molecule,” The Journal of Physical Chemistry B 118, 8148–8161 (2014), pMID: 24601590, https://doi.org/10.1021/jp501182d .
  * Imoto, Xantheas, and Saito [2013] S. Imoto, S. S. Xantheas, and S. Saito, “Ultrafast dynamics of liquid water: Frequency fluctuations of the OH stretch and the HOH bend,” The Journal of Chemical Physics 139, 044503 (2013), https://doi.org/10.1063/1.4813071 .
  * Piryatinski, Lawrence, and Skinner [2003] A. Piryatinski, C. P. Lawrence, and J. L. Skinner, “Vibrational spectroscopy of HOD in liquid D2O. v. infrared three-pulse photon echoes,” The Journal of Chemical Physics 118, 9672–9679 (2003), https://doi.org/10.1063/1.1569474 .
  * Paarmann  _et al._ [2009] A. Paarmann, T. Hayashi, S. Mukamel, and R. J. D. Miller, “Nonlinear response of vibrational excitons: Simulating the two-dimensional infrared spectrum of liquid water,” The Journal of Chemical Physics 130, 204110 (2009), https://doi.org/10.1063/1.3139003 .
  * Jansen  _et al._ [2010] T. l. C. Jansen, B. M. Auer, M. Yang, and J. L. Skinner, “Two-dimensional infrared spectroscopy and ultrafast anisotropy decay of water,” The Journal of Chemical Physics 132, 224503 (2010).
  * Gottwald, Ivanov, and Kühn [2015] F. Gottwald, S. D. Ivanov, and O. Kühn, “Applicability of the Caldeira–Leggett model to vibrational spectroscopy in solution,” The Journal of Physical Chemistry Letters 6, 2722–2727 (2015).
  * Gottwald  _et al._ [2015] F. Gottwald, S. Karsten, S. D. Ivanov, and O. Kühn, “Parametrizing linear generalized Langevin dynamics from explicit molecular dynamics simulations,” The Journal of Chemical Physics 142, 244110 (2015).
  * Gottwald, Ivanov, and Kühn [2016] F. Gottwald, S. D. Ivanov, and O. Kühn, “Vibrational spectroscopy via the Caldeira–Leggett model with anharmonic system potentials,” The Journal of Chemical Physics 144, 164102 (2016).
  * Yagasaki, Ono, and Saito [2009] T. Yagasaki, J. Ono, and S. Saito, “Ultrafast energy relaxation and anisotropy decay of the librational motion in liquid water: A molecular dynamics study,” The Journal of Chemical Physics 131, 164511 (2009), https://doi.org/10.1063/1.3254518 .
  * Yagasaki and Saito [2011] T. Yagasaki and S. Saito, “A novel method for analyzing energy relaxation in condensed phases using nonequilibrium molecular dynamics simulations: Application to the energy relaxation of intermolecular motions in liquid water,” The Journal of Chemical Physics 134, 184503 (2011), https://doi.org/10.1063/1.3587105 .
  * Okumura and Tanimura [1997a] K. Okumura and Y. Tanimura, “Two-time correlation functions of a harmonic system nonbilinearly coupled to a heat bath: Spontaneous Raman spectroscopy,” Phys. Rev. E 56, 2747–2750 (1997a).
  * Tanimura and Steffen [2000] Y. Tanimura and T. Steffen, “Two-dimensional spectroscopy for harmonic vibrational modes with nonlinear system-bath interactions.ii. Gaussian-Markovian case,” Journal of the Physical Society of Japan 69, 4095–4106 (2000).
  * Kato and Tanimura [2004] T. Kato and Y. Tanimura, “Two-dimensional Raman and infrared vibrational spectroscopy for a harmonic oscillator system nonlinearly coupled with a colored noise bath,” The Journal of Chemical Physics 120, 260–271 (2004).
  * Ishizaki and Tanimura [2007] A. Ishizaki and Y. Tanimura, “Dynamics of a multimode system coupled to multiple heat baths probed by two-dimensional infrared spectroscopy,” The Journal of Physical Chemistry A 111, 9269–9276 (2007).
  * Tanimura and Kubo [1989] Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath,” Journal of the Physical Society of Japan 58, 101–114 (1989).
  * Tanimura [2006] Y. Tanimura, “Stochastic Liouville, Langevin, Fokker-Planck, and master equation qpproaches to quantum dissipative systems,” Journal of the Physical Society of Japan 75, 082001 (2006).
  * Tanimura [2020] Y. Tanimura, “Numerically ”exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM),” The Journal of Chemical Physics 153, 020901 (2020).
  * Hasegawa and Tanimura [2011] T. Hasegawa and Y. Tanimura, “A polarizable water model for intramolecular and intermolecular vibrational spectroscopies,” The Journal of Physical Chemistry B 115, 5545–5553 (2011).
  * Liu and Liu [2018] X. Liu and J. Liu, “Critical role of quantum dynamical effects in the Raman spectroscopy of liquid water,” Molecular Physics 116, 755–779 (2018), https://doi.org/10.1080/00268976.2018.1434907 .
  * Liu, Wang, and Bowman [2015] H. Liu, Y. Wang, and J. M. Bowman, “Quantum calculations of the ir spectrum of liquid water using ab initio and model potential and dipole moment surfaces and comparison with experiment,” The Journal of Chemical Physics 142, 194502 (2015), https://doi.org/10.1063/1.4921045 .
  * Hunter, Shakib, and Paesani [2018] K. M. Hunter, F. A. Shakib, and F. Paesani, “Disentangling coupling effects in the infrared spectra of liquid water,” The Journal of Physical Chemistry B 122, 10754–10761 (2018), pMID: 30403350, https://doi.org/10.1021/acs.jpcb.8b09910 .
  * Trenins, Willatt, and Althorpe [2019] G. Trenins, M. J. Willatt, and S. C. Althorpe, “Path-integral dynamics of water using curvilinear centroids,” The Journal of Chemical Physics 151, 054109 (2019), https://doi.org/10.1063/1.5100587 .
  * Ueno and Tanimura [2020] S. Ueno and Y. Tanimura, “Modeling intermolecular and intramolecular modes of liquid water using multiple heat baths: Machine learning approach,” Journal of Chemical Theory and Computation 16, 2099–2108 (2020).
  * Okumura and Tanimura [1997b] K. Okumura and Y. Tanimura, “The (2n+1)th-order off-resonant spectroscopy from the (n+1)th-order anharmonicities of molecular vibrational modes in the condensed phase,” The Journal of Chemical Physics 106, 1687–1698 (1997b).
  * Okumura and Tanimura [1997c] K. Okumura and Y. Tanimura, “Femtosecond two-dimensional spectroscopy from anharmonic vibrational modes of molecules in the condensed phase,” The Journal of Chemical Physics 107, 2267–2283 (1997c).
  * Okumura and Tanimura [1997d] K. Okumura and Y. Tanimura, “Sensitivity of two-dimensional fifth-order Raman response to the mechanism of vibrational mode-mode coupling in liquid molecules,” Chemical Physics Letters 278, 175–183 (1997d).
  * Okumura, Jonas, and Tanimura [2001] K. Okumura, D. M. Jonas, and Y. Tanimura, “Two-dimensional spectroscopy and harmonically coupled anharmonic oscillators,” Chemical Physics 266, 237–250 (2001).
  * Schmidt, Corcelli, and Skinner [2005] J. R. Schmidt, S. A. Corcelli, and J. L. Skinner, “Pronounced non-condon effects in the ultrafast infrared spectroscopy of water,” The Journal of Chemical Physics 123, 044513 (2005), https://doi.org/10.1063/1.1961472 .
  * Palese  _et al._ [1994] S. Palese, J. T. Buontempo, L. Schilling, W. T. Lotshaw, Y. Tanimura, S. Mukamel, and R. J. D. Miller, “Femtosecond two-dimensional Raman spectroscopy of liquid water,” The Journal of Physical Chemistry 98, 12466–12470 (1994).
  * Hybl, Albrecht Ferro, and Jonas [2001] J. D. Hybl, A. Albrecht Ferro, and D. M. Jonas, “Two-dimensional fourier transform electronic spectroscopy,” The Journal of Chemical Physics 115, 6606–6622 (2001), https://doi.org/10.1063/1.1398579 .
  * Ge, Zanni, and Hochstrasser [2002] N.-H. Ge, M. T. Zanni, and R. M. Hochstrasser, “Effects of vibrational frequency correlations on two-dimensional infrared spectra†,” Journal of Physical Chemistry A 106, 962–972 (2002).
  * Khalil, Demirdöven, and Tokmakoff [2003] M. Khalil, N. Demirdöven, and A. Tokmakoff, “Obtaining absorptive line shapes in two-dimensional infrared vibrational correlation spectra,” Phys. Rev. Lett. 90, 047401 (2003).
  * Hasegawa and Tanimura [2008] T. Hasegawa and Y. Tanimura, “Nonequilibrium molecular dynamics simulations with a backward-forward trajectories sampling for multidimensional infrared spectroscopy of molecular vibrational modes,” The Journal of Chemical Physics 128, 064511 (2008).
  * Yagasaki and Saito [2008] T. Yagasaki and S. Saito, “Ultrafast intermolecular dynamics of liquid water: A theoretical study on two-dimensional infrared spectroscopy,” The Journal of Chemical Physics 128, 154521 (2008), https://doi.org/10.1063/1.2903470 .
  * Kato and Tanimura [2001] T. Kato and Y. Tanimura, “Multi-dimensional vibrational spectroscopy measured from different phase-matching conditions,” Chemical Physics Letters 341, 329–337 (2001).
  * Lazonder, Pshenichnikov, and Wiersma [2006] K. Lazonder, M. S. Pshenichnikov, and D. A. Wiersma, “Easy interpretation of optical two-dimensional correlation spectra,” Opt. Lett. 31, 3354–3356 (2006).
  * Nienhuys  _et al._ [1999] H.-K. Nienhuys, S. Woutersen, R. A. van Santen, and H. J. Bakker, “Mechanism for vibrational relaxation in water investigated by femtosecond infrared spectroscopy,” The Journal of Chemical Physics 111, 1494–1500 (1999), https://doi.org/10.1063/1.479408 .
  * Ikeda and Tanimura [2018] T. Ikeda and Y. Tanimura, “Phase-space wavepacket dynamics of internal conversion via conical intersection: Multi-state quantum Fokker-Planck equation approach,” Chemical Physics 515, 203–213 (2018).


# The Effect of Model Size on LLM Post-hoc Explainability via LIME

Henning Heyen, Amy Widdicombe, Noah Y. Siegel, María Pérez-Ortiz, Philip
Treleaven
University College London
Contact<EMAIL_ADDRESS>
###### Abstract

Large language models (LLMs) are becoming bigger to boost performance.
However, little is known about how explainability is affected by this trend.
This work explores LIME explanations for DeBERTaV3 models of four different
sizes on natural language inference (NLI) and zero-shot classification (ZSC)
tasks. We evaluate the explanations based on their faithfulness to the models’
internal decision processes and their plausibility, i.e. their agreement with
human explanations. The key finding is that increased model size does not
correlate with plausibility despite improved model performance, suggesting a
misalignment between the LIME explanations and the models’ internal processes
as model size increases. Our results further suggest limitations regarding
faithfulness metrics in NLI contexts.

##  1 Introduction

Research has shown that performance in language models depends strongly on
scale and less on model shape (Kaplan et al., 2020), where scale refers to the
number of parameters, the training dataset size, and the amount of compute for
training. For instance, OpenAI’s series of Generative Pre-Trained Transformers
(GPT) has grown from 1.5 billion parameters for GPT-2 to 175 billion
parameters for GPT-3 which helped improve across various NLP tasks (Brown et
al., 2020). This trend seems likely to continue.

As LLMs grow in size and performance and are increasingly deployed in high-
stakes applications, the need to understand and explain their behaviour
becomes more crucial. Post-hoc explainability methods such as LIME (Ribeiro et
al., 2016) are one way of attempting to do this. Although these methods have
been widely applied to LLMs (Madsen et al., 2022), to the best of our
knowledge no research has been conducted on the impact of model size on the
quality of these kinds of explanations. Here we begin to fill this gap by
investigating the impact of model size on the quality of LIME explanations. We
apply two approaches to assess the quality of explanations, namely
faithfulness (Chan et al., 2022) and plausibility (DeYoung et al., 2020).
While faithfulness aims to measure the extent to which an explanation reflects
the true internal decision processes, plausibility assesses the quality of the
explanations based on their agreement with human-generated explanations.

We find that, even though model performance increases with model size, the
agreement between human-generated and LIME-generated explanations does not.
This indicates some extent of misalignment between the explanations and the
true internal decision processes. Our findings also imply possible flaws in
removal-based faithfulness metrics based on the NLP task which points to more
general limitations for highlight-based post-hoc explainability in NLP such as
lack of expressiveness. This study serves as a first attempt to understand how
post-hoc explainability is affected by model size. We hope that this research
encourages others to further explore this area and to that end we provide an
extensible code repository111https://github.com/henningheyen/Scalability-Of-
LLM-Posthoc-Explanations for others to build on.

##  2 Methodology

#### Models and Datasets

We use fine-tuned DeBERTaV3 models from Huggingface of four different sizes,
ranging from 22 to 304 million parameters222For architectural specifics see
Table 4 in the appendix.. Note that state-of-the-art models exhibit parameter
counts in the order of billions. Due to computational constraints, we could
not experiment with larger models. The models were fine-tuned on two standard
natural language inference (NLI) datasets, matched MNLI (Williams et al.,
2018) and SNLI (Bowman et al., 2015). Instead of SNLI, we use e-SNLI (Camburu
et al., 2018) which extends SNLI by human annotated highlights indicating the
most important tokens with respect to the label. Additionally, we apply the
models in a zero-shot classification (ZSC) setting using the CoS-e (Rajani et
al., 2019) dataset. CoS-e consists of commonsense questions with five
candidate labels where the candidate labels differ for each question333The
Huggingface API internally transforms zero-shot classification problems to an
NLI problem.. Similarly to e-SNLI, CoS-e contains human annotated
highlights444For examples from MNLI, e-SNLI and CoS-e refer to Table 5, 6 and
7 in the appendix..

#### Explainability Method

There exist different notions of explainability in NLP. One is in the form of
free-text natural language explanations and another is in the form of
highlight-based, typically post-hoc explanations. While the former approach
commonly requires human evaluation the latter can be measured more objectively
because post-hoc techniques in NLP are normally mappings from tokens to real-
valued importance scores. Common techniques can be categorised as gradient-
based, attention weight-based and perturbation-based. While gradient-based
techniques are highly vulnerable to adversarial perturbation (Wang et al.,
2020) several studies argue that explanations based on attention weights are
unreliable. The study “Attention is not Explanation“ (Jain & Wallace, 2019),
for instance, identified different attention distributions yielding equivalent
predictions. We believe perturbation-based methods avoid some of these
pitfalls - we use one such method, LIME (Ribeiro et al., 2016), in our
experiments. LIME is a well-known and widely used explainability technique
that is based on locally approximating the prediction with a simple
interpretable model. The weights of the simple model serve as feature
importance scores. Figure 1 illustrates a LIME explanation on a CoS-e
instance.

Figure 1: LIME example on a CoS-e instance using the xsmall DeBERTaV3 model.
LIME maps every token to a real-valued importance score.

#### Explainability Metrics

There is no single framework for evaluating post-hoc explanations due to a
lack of consensus on what constitutes a high-quality explanation. In this
work, we examine two different approaches, namely faithfulness and
plausibility:

Faithfulness as discussed in Chan et al. (2022) aims to measure to what extent
the explanation reflects the model’s internal decision process. Generally,
faithfulness metrics rely on removing tokens from the input sequence based on
the explanation and measuring the change in prediction. While several
faithfulness metrics exist, they are “not always consistent with each other
and even lead to contradictory conclusions“ (Chan et al., 2022) \- they
compared six faithfulness metrics and concluded that comprehensiveness is the
most diagnostic and the least complex. Based on this we use comprehensiveness
for our experiments. Proposed by DeYoung et al. (2020), comprehensiveness
suggests that an explanation is faithful if the prediction strongly deviates
when the most important tokens (as identified by the explanation method) are
removed from the input sequence555For a visualization of the comprehensiveness
metric see Figure 3 in the appendix.. More formally,

$\text{COMP}(\mathbf{x},c,k)=p(c\mid\mathbf{x};\theta)-p(c\mid\mathbf{x}\backslash\mathbf{x}_{k};\theta),$
(1)

where $p(c\mid\mathbf{x};\theta)$ denotes the model’s prediction for class $c$
on the entire input sequence and
$p(c\mid\mathbf{x}\backslash\mathbf{x}_{k};\theta)$ denotes the model’s
prediction when the top $k$ most important tokens $\mathbf{x}_{k}$ are removed
from the input string. Practically, we obtain the most important tokens by
taking the top $t$ percent tokens from the list of tokens with associated
real-valued importance scores generated by a post-hoc explainability method
such as LIME. We denote this explanation as $\mathbf{x}_{t}$. To enhance the
metric’s reliability, the original paper proposes aggregated comprehensiveness
which averages the comprehensiveness over different lengths of explanations.
In this work, we use bins $t\in T=\\{10\%,30\%,50\%\\}$ to vary the length of
explanations. The aggregated comprehensiveness can be defined as,

$\text{COMP}_{\text{agg}}(\mathbf{x},c)=\sum_{t\in
T}\text{COMP}(\mathbf{x},c,t)=\sum_{t\in
T}\left(p(c\mid\mathbf{x};\theta)-p(c\mid\mathbf{x}\setminus\mathbf{x}_{t};\theta)\right).$
(2)

Plausibility, as compared to faithfulness, defines the quality of an
explanation by the intersection between the highlights generated by a post-hoc
explainability technique and some human-annotated highlights. In other words,
plausibility measures the “agreement between extracted and human rationale“
(DeYoung et al., 2020). Plausibility fundamentally differs from faithfulness
in that plausible explanations do not reveal whether the model actually relied
on the explanation. Assessing plausibility typically requires human evaluation
(Strout et al., 2019). However, more recently some existing datasets have been
extended by human-annotated highlights (DeYoung et al., 2020) which allows for
more quantitative evaluation of plausibility. In this paper, we use two
datasets with human-annotated highlights, namely CoS-e (Rajani et al., 2019)
and e-SNLI (Camburu et al., 2018). Similarly to DeYoung et al. (2020) we
measure plausibility by the intersection over union (IOU). As proposed by
DeYoung et al. (2020), we take the number of most important tokens according
to the average explanation length provided by humans for each dataset666Mean
explanation-input-ratio e-SNLI: $0.19$ ($\pm 0.193$), CoS-e: $0.26$ ($\pm
0.137$).. Suppose $\mathbf{x}_{1}$ is the set of tokens from the human
explanation and $\mathbf{x}_{2}$ is the set of generated tokens. Then IOU can
be formalised by

$\text{IOU}(\mathbf{x}_{1},\mathbf{x}_{2})=\frac{\left|\mathbf{x}_{1}\cap\mathbf{x}_{2}\right|}{\left|\mathbf{x}_{1}\cup\mathbf{x}_{2}\right|}.$
(3)

##  3 Experiments and Results

The results for Experiment 1 and 2 are summarised in Table 1. Experiment 3
results are visualised in Figure 2.

Experiment 1 First, the four DeBERTaV3 models were evaluated in terms of
performance on the validation sets of all three datasets (MNLI, e-SNLI,
CoS-e). We report on model accuracy and 95% confidence intervals. We find
that, as expected, performance improves monotonically with increasing model
size for all three datasets. We can conclude that the models’ capabilities are
different enough to reason about the effect of model size on the LIME
explanations.

Experiment 2 We then computed LIME explanations for each model on a subset of
100 test samples from each dataset. We had to use a subset due to the
computational intensity of LIME. The explanations were always calculated with
respect to the predicted class, not necessarily the correct class. For each
explanation, the aggregated comprehensiveness and IOU scores were computed
according to equation 2 and 3 respectively. Note that IOU scores were only
feasible for e-SNLI and CoS-e instances since MNLI does not provide human-
annotated highlights. For each model, we report on the mean comprehensiveness
and mean IOU scores across all 100 explanations. We observe an overall
increase in comprehensiveness with the highest scores on the largest model for
all three datasets suggesting that faithfulness of LIME increases with model
size. IOU, on the other hand, stays almost constant across all model sizes for
both datasets indicating that the plausibility of the LIME explanations is
uncorrelated with model size.

Dataset | Model Size | Comprehensiveness | IOU | Accuracy | 95% C.I.
---|---|---|---|---|---
MNLI | xsmall | 0.785 ($\pm$ 0.022) | – | 0.878 | (0.871, 0.885)
small | 0.817 ($\pm$ 0.022) | – | 0.878 | (0.872, 0.884)
base | 0.796 ($\pm$ 0.027) | – | 0.900 | (0.894, 0.906)
large |  0.823 ($\pm$ 0.027) | – | 0.902 | (0.896, 0.908)
e-SNLI | xsmall | 0.726 ($\pm$ 0.022) |  0.282 ($\pm$ 0.017) | 0.920 | (0.915, 0.925)
small | 0.724 ($\pm$ 0.026) | 0.259 ($\pm$ 0.016) | 0.922 | (0.917, 0.927)
base | 0.764 ($\pm$ 0.025) | 0.254 ($\pm$ 0.016) | 0.931 | (0.926, 0.936)
large |  0.778 ($\pm$ 0.025) | 0.256 ($\pm$ 0.017) | 0.932 | (0.927, 0.937)
CoS-e | xsmall | 0.304 ($\pm$ 0.018) | 0.233 ($\pm$ 0.013) | 0.331 | (0.305, 0.355)
small | 0.316 ($\pm$ 0.019) | 0.231 ($\pm$ 0.014) | 0.336 | (0.306, 0.362)
base | 0.356 ($\pm$ 0.020) |  0.235 ($\pm$ 0.012) | 0.359 | (0.330, 0.383)
large |  0.391 ($\pm$ 0.022) | 0.230 ($\pm$ 0.012) | 0.378 | (0.349, 0.406)
Table 1: Mean comprehensiveness and IOU scores with mean standard errors on
100 test samples for each dataset across all model sizes and accuracy scores
on full validation sets with 95% confidence intervals. IOU could not be
computed on MNLI as this dataset does not provide human annotated highlights
as ground truth explanations.

Experiment 3 Lastly, we investigated both metrics with respect to the labels
(entailment, neutral, contradiction) for MNLI and e-SNLI777Exact numbers are
shown in Table 2 in the appendix, the labels are balanced, see Table 8.. The
goal was to see how the metrics change with model size when we condition on
the label. We observe that comprehensiveness improves for contradictory
sentence pairs with larger model sizes in MNLI, while no consistent pattern
emerges in e-SNLI. Generally, we find that neutral sentence pairs achieved
lower comprehensiveness scores than contradiction pairs. IOU stays almost
constant across different model sizes regardless of the label.

Figure 2: Mean comprehensiveness scores by labels for MNLI (left), e-SNLI
(middle) and mean IOU scores by labels for e-SNLI (right) with mean standard
errors on 100 test samples across all model sizes. Note how neutral sentence
pairs achieve generally lower comprehensiveness scores than contradictive
sentence pairs and how IOU scores are almost constant as the model size
increases regardless of the label. Exact numbers are displayed in Table 2 in
the appendix.

##  4 Discussion

Overall, we find that for all three datasets, the largest model achieved the
highest comprehensiveness score suggesting that with LIME larger models yield
more faithful explanations. However, the IOU score stays constant across
different model sizes suggesting that the plausibility of the explanations is
uncorrelated with model size and performance. Interestingly, this would imply
that the agreement with human-annotated highlights does not improve with model
performance which indicates an inherent misalignment between the generated
explanations and the true internal decision process. This finding seems
contradictory to our previous result that with LIME larger models yield more
faithful explanations. Splitting the metrics by labels could reveal potential
flaws with the comprehensiveness metric in the NLI setting as we found
significantly lower scores for neutral sentence pairs. The problem with
comprehensiveness in an NLI setting could be that removing highlighted tokens
from a neutral pair might very well result in another neutral prediction which
limits the applicability of comprehensiveness in this case. More generally,
this shows how post-hoc explanations in NLP lack expressiveness. Highlighting
may not be sufficient to fully explain LLMs. While highlighting can provide
some intuition about which tokens LLMs pay attention to, some higher-level
reasoning concepts such as token dependencies or logical relations can hardly
be expressed by highlighted tokens. This observation limits our finding that
LIME explanations are more faithful for larger models. We conclude that the
applicability of comprehensiveness is task-dependent and more coherent and
expressive explainability metrics and techniques are needed.

##  5 Future Work

The exclusive focus on LIME cannot capture the full range of challenges and
opportunities that other techniques might present for explainability. Future
research could repeat our experiments using other perturbation-based post-hoc
techniques such as Anchors (Ribeiro et al., 2018) or SHAP (Lundberg & Lee,
2017) to validate our observations. Additionally, other tasks such as
sentiment analysis, text summarisation or language modelling could be
explored. Furthermore, none of our models have gone through Reinforcement
Learning from Human Feedback (RLHF) (Ouyang et al., 2022). We believe that the
results on plausibility might significantly change with RLHF. While this work
explored explainability primarily quantitatively, future investigations may
benefit from qualitative analysis to deepen understanding of why discrepancies
between model explanations and human understanding occur. Note that our large
models had 304 million parameters. The current state-of-the-art models,
however, are much bigger, consisting of billions of parameters. While
computational feasibility with perturbation-based techniques is one concern,
the effect on the explainability of very large models is still unexplored and
should be the subject of future investigation. More broadly, this work aims to
emphasize the urgent need for an LLM explainability framework that is human-
interpretable, objective, scalable and expressive. A promising direction from
a recent study showed empirically that the faithfulness of post hoc
explanations could improve when the models are optimized for robustness
against adversarial attacks (Li et al., 2023). The method incorporates the
explanations into the model training process.

##  6 Conclusion

Our work investigated LIME explanations on fine-tuned DeBERTaV3 models of
different sizes in an NLI and ZSC setting. We applied two approaches to
capture the quality of explanations, namely faithfulness and plausibility. We
measured faithfulness with comprehensiveness and plausibility with IOU. We
identified limitations of the comprehensiveness metric in the NLI setting
which suggests that removal-based faithfulness metrics are task-dependent and
that their general applicability is questionable. Given that performance
increased with model size raises questions on why agreement with human
annotations does not increase. We suggest that there is some extent of
misalignment between the model’s internal decision process and its LIME
explanation. This work aims to serve as an initial step towards understanding
the effect of model size on LLM post-hoc explainability. We believe there is
an urgent need for further investigations on LLM explainability and a more
coherent explainability framework for LLMs. If we fail to faithfully explain
the decisions of increasingly large language models we risk that those models
pursue unexpected objectives rather than agreeing with human values and
intentions.

## References

  * Bowman et al. (2015) Samuel R. Bowman, Gabor Angeli, Christopher Potts, and Christopher D. Manning.  A large annotated corpus for learning natural language inference.  In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing_ , pp. 632–642, Lisbon, Portugal, September 2015. Association for Computational Linguistics.  doi: 10.18653/v1/D15-1075.  URL https://aclanthology.org/D15-1075.
  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901, 2020.
  * Camburu et al. (2018) Oana-Maria Camburu, Tim Rocktäschel, Thomas Lukasiewicz, and Phil Blunsom.  e-snli: Natural language inference with natural language explanations.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Chan et al. (2022) Chun Sik Chan, Huanqi Kong, and Liang Guanqing.  A comparative study of faithfulness metrics for model interpretability methods.  In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pp. 5029–5038, Dublin, Ireland, May 2022. Association for Computational Linguistics.  doi: 10.18653/v1/2022.acl-long.345.  URL https://aclanthology.org/2022.acl-long.345.
  * DeYoung et al. (2020) Jay DeYoung, Sarthak Jain, Nazneen Fatema Rajani, Eric Lehman, Caiming Xiong, Richard Socher, and Byron C. Wallace.  ERASER: A benchmark to evaluate rationalized NLP models.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pp. 4443–4458, Online, July 2020. Association for Computational Linguistics.  doi: 10.18653/v1/2020.acl-main.408.  URL https://aclanthology.org/2020.acl-main.408.
  * Jain & Wallace (2019) Sarthak Jain and Byron C. Wallace.  Attention is not Explanation.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pp. 3543–3556, Minneapolis, Minnesota, June 2019. Association for Computational Linguistics.  doi: 10.18653/v1/N19-1357.  URL https://aclanthology.org/N19-1357.
  * Kaplan et al. (2020) Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei.  Scaling laws for neural language models.  _arXiv preprint arXiv:2001.08361_ , 2020.
  * Li et al. (2023) Dongfang Li, Baotian Hu, Qingcai Chen, and Shan He.  Towards faithful explanations for text classification with robustness improvement and explanation guided training.  In Anaelia Ovalle, Kai-Wei Chang, Ninareh Mehrabi, Yada Pruksachatkun, Aram Galystan, Jwala Dhamala, Apurv Verma, Trista Cao, Anoop Kumar, and Rahul Gupta (eds.), _Proceedings of the 3rd Workshop on Trustworthy Natural Language Processing (TrustNLP 2023)_ , pp. 1–14, Toronto, Canada, July 2023. Association for Computational Linguistics.  doi: 10.18653/v1/2023.trustnlp-1.1.  URL https://aclanthology.org/2023.trustnlp-1.1.
  * Lundberg & Lee (2017) Scott M. Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17, pp. 4768–4777, Red Hook, NY, USA, 2017. Curran Associates Inc.  ISBN 9781510860964.
  * Madsen et al. (2022) Andreas Madsen, Siva Reddy, and Sarath Chandar.  Post-hoc interpretability for neural nlp: A survey.  _ACM Computing Surveys_ , 55(8):1–42, 2022.
  * Ouyang et al. (2022) Long Ouyang, Jeffrey Wu, Xu Jiang, Diogo Almeida, Carroll Wainwright, Pamela Mishkin, Chong Zhang, Sandhini Agarwal, Katarina Slama, Alex Ray, et al.  Training language models to follow instructions with human feedback.  _Advances in Neural Information Processing Systems_ , 35:27730–27744, 2022.
  * Rajani et al. (2019) Nazneen Fatema Rajani, Bryan McCann, Caiming Xiong, and Richard Socher.  Explain yourself! leveraging language models for commonsense reasoning.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pp. 4932–4942, Florence, Italy, July 2019. Association for Computational Linguistics.  doi: 10.18653/v1/P19-1487.  URL https://aclanthology.org/P19-1487.
  * Ribeiro et al. (2016) 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_ , pp. 1135–1144, 2016.
  * Ribeiro et al. (2018) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  Anchors: High-precision model-agnostic explanations.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 32, 2018.
  * Strout et al. (2019) Julia Strout, Ye Zhang, and Raymond Mooney.  Do human rationales improve machine explanations?  In _Proceedings of the 2019 ACL Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP_ , pp. 56–62, Florence, Italy, August 2019. Association for Computational Linguistics.  doi: 10.18653/v1/W19-4807.  URL https://aclanthology.org/W19-4807.
  * Wang et al. (2020) Junlin Wang, Jens Tuyls, Eric Wallace, and Sameer Singh.  Gradient-based analysis of NLP models is manipulable.  In _Findings of the Association for Computational Linguistics: EMNLP 2020_ , pp. 247–258, Online, November 2020. Association for Computational Linguistics.  doi: 10.18653/v1/2020.findings-emnlp.24.  URL https://aclanthology.org/2020.findings-emnlp.24.
  * Williams et al. (2018) Adina Williams, Nikita Nangia, and Samuel Bowman.  A broad-coverage challenge corpus for sentence understanding through inference.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pp. 1112–1122, New Orleans, Louisiana, June 2018. Association for Computational Linguistics.  doi: 10.18653/v1/N18-1101.  URL https://aclanthology.org/N18-1101.

##  Appendix A Appendix

###  A.1 Additional Tables

Dataset | Model Size | Label | Comprehensiveness | IOU
---|---|---|---|---
MNLI | xsmall | contradiction | 0.810 ($\pm$ 0.038) | -
entailment | 0.759 ($\pm$ 0.042) | -
neutral | 0.788 ($\pm$ 0.034) | -
small | contradiction | 0.853 ($\pm$ 0.041) | -
entailment | 0.805 ($\pm$ 0.039) | -
neutral | 0.794 ($\pm$ 0.035) | -
base | contradiction | 0.895 ($\pm$ 0.034) | -
entailment | 0.768 ($\pm$ 0.050) | -
neutral | 0.728 ($\pm$ 0.049) | -
large | contradiction | 0.939 ($\pm$ 0.018) | -
entailment | 0.750 ($\pm$ 0.057) | -
neutral | 0.789 ($\pm$ 0.046) | -
e-SNLI | xsmall | contradiction | 0.805 ($\pm$ 0.034) | 0.289 ($\pm$ 0.031)
entailment | 0.713 ($\pm$ 0.039) | 0.315 ($\pm$ 0.025)
neutral | 0.663 ($\pm$ 0.035) | 0.244 ($\pm$ 0.028)
small | contradiction | 0.744 ($\pm$ 0.042) | 0.286 ($\pm$ 0.029)
entailment | 0.783 ($\pm$ 0.051) | 0.249 ($\pm$ 0.025)
neutral | 0.652 ($\pm$ 0.040) | 0.242 ($\pm$ 0.030)
base | contradiction | 0.808 ($\pm$ 0.038) | 0.264 ($\pm$ 0.027)
entailment | 0.786 ($\pm$ 0.043) | 0.291 ($\pm$ 0.025)
neutral | 0.701 ($\pm$ 0.045) | 0.211 ($\pm$ 0.028)
large | contradiction | 0.759 ($\pm$ 0.046) | 0.259 ($\pm$ 0.024)
entailment | 0.809 ($\pm$ 0.038) | 0.292 ($\pm$ 0.026)
neutral | 0.768 ($\pm$ 0.042) | 0.220 ($\pm$ 0.036)
Table 2: Mean comprehensiveness and IOU scores with mean standard errors on 100 test samples split by label for both NLI datasets across all model sizes. IOU could not be computed on MNLI as this dataset does not provide human-annotated highlights as ground truth explanations. | MNLI | e-SNLI | CoS-e
---|---|---|---
xsmall | 2min 3s | 1min 8s | 34min 35s
small | 2min 40s | 1min 40s | 44min 28s
base | 5min 20s | 3min 35s | 1h 27min 7s
large | 15min 38s | 12min 45s | 4h 35min 50s
Table 3: Compute time for all LIME explanations of 100 test instances from each dataset across all model sizes on Nvidia’s T4 GPU, 51GB RAM. | Parameters | Layers | Hidden | Attention
---|---|---|---|---
| (in millions) |  | Size | Heads
large | 304 | 24 | 1024 | 12
base | 86 | 12 | 768 | 12
small | 44 | 6 | 768 | 12
xsmall | 22 | 12 | 384 | 6
Table 4: Architecture comparison for DeBERTaV3 models. Premise |  Hypothesis |  Label
---|---|---
Look, there’s a legend here. |  See, there is a well-known hero here. |  Entailment
Yeah, I know, and I did that all through college and it worked too. |  I did that all through college but it never worked. |  Contradiction
Boats in daily use lie within feet of the fashionable bars and restaurants. |  Bars and restaurants are interesting places. |  Neutral
Table 5: Natural language inference examples from the MNLI dataset. Premise |  Hypothesis |  Label
---|---|---
An adult dressed in black holds a stick. |  An adult is walking away, empty-handed. |  Contradiction
A child in a yellow plastic safety swing is laughing as a dark-haired woman in pink and coral pants stands behind her. |  A young mother is playing with her daughter in a swing. |  Neutral
A man in an orange vest leans over a pickup truck. |  A man is touching a truck. |  Entailment
Table 6: Natural language inference examples from the e-SNLI dataset. Highlighted tokens indicate human-annotated explanations. Question |  Candidate Labels |  Label
---|---|---
He was a sloppy eater, so where did he leave a mess? |  sailboat, desk, closet, table, apartment |  table
Where can someone get a new saw? |  hardware store, toolbox, logging camp, tool kit, auger |  hardware store
Many homes in this country are built around a courtyard. Where is it? |  hospital, park, spain, office complex, office |  spain
Table 7: Zero shot classification examples from the CoS-e dataset. Highlighted tokens indicate human-annotated explanations. Dataset | Contradiction | Entailment | Neutral
---|---|---|---
MNLI | 32 | 36 | 32
e-SNLI | 33 | 32 | 35
Table 8: Number of observations by label for MNLI and e-SNLI for 100 explained
test samples.

###  A.2 Additional Figures

Figure 3: Visualistion of comprehensiveness metric on CoS-e instance from
(DeYoung et al., 2020). Comprehensiveness suggests that an explanation is
faithful if the prediction strongly deviates when the most important tokens
(as identified by the explanation method) are removed from the input sequence.
Figure 4: Mean comprehensiveness (left) and IOU (right) scores with mean
standard errors on 100 test samples for each dataset across all model sizes.
IOU could not be computed on MNLI as this dataset does not provide human-
annotated highlights as ground truth explanations.