#  3D Grid-Based Monte Carlo Code for Radiative Transfer through Raman and
Rayleigh Scattering with Atomic Hydrogen – STaRS

Seok-Jun Chang  Department of Physics and Astronomy, Sejong University, Seoul,
Korea;  Korea Astronomy and Space Science Institute, Daejeon, Korea  Hee-Won
Lee  Department of Physics and Astronomy, Sejong University, Seoul, Korea;
<EMAIL_ADDRESS>
(Received September 30, 2020; accepted , 2020)

##  1 Introduction

Scattered radiation conveys special information regarding both the emission
and scattering regions. In particular, linear polarization may develop in an
anisotropic scattering geometry and the relative motion between the emission
and scattering regions gives rise to complicated line profiles including
multiple-peak structures and broad wings. An excellent example is provided by
spectropolarimetric observations of Seyfert 2 galaxies exhibiting broad lines
in the linearly polarized fluxes, lending strong support to the unification
model of active galactic nuclei (Miller & Goodrich, 1990; Tran, 2010).

Interaction of electromagnetic radiation with an atomic electron may be
classified into Rayleigh and Raman scattering (e.g. Sakurai, 1967). Rayleigh
scattering refers to the elastic process where the scattered photon has the
same wavelength as the incident one. Otherwise, we have Raman scattering in
which the initial and final electronic states differ so that the scattered
photon emerges with the energy difference, which is enforced by the law of
energy conservation. Raman spectroscopy is particularly useful in revealing
the complicated energy level structures of a molecule in chemistry.

In astrophysics, Nussbaumer et al. (1989) provided a pioneering discussion on
Raman scattering with atomic hydrogen, introducing a new spectroscopic tool to
diagnose gaseous emission nebulae including symbiotic stars and active
galactic nuclei. Relevant Raman scattering processes start with a far UV
photon more energetic than Ly$\alpha$ incident on an hydrogen atom in the
ground state and end with an outgoing photon with energy less than that of
Ly$\alpha$ leaving behind the hydrogen atom in the $2s$ state. They introduced
basic atomic physics of Rayleigh and Raman scattering illustrating the cross
sections and presented a number of candidate far UV spectral lines that may
result in detectable Raman-scattered features (e.g., Saslow & Mills, 1969).
One notable point is that the branching ratio of Raman and Rayleigh scattering
of far UV radiation near Ly$\beta$ is approximately 0.15 so that most
Ly$\beta$ photons are Rayleigh (or resonantly) scattered several times before
they are converted into H$\alpha$ photons to escape from the thick neutral
region.

Schmid (1989) identified the broad emission features at 6825 Å and 7082 Å
found in about half of symbiotic stars as Raman scattered features of O VI
$1032$ and $1038$ Å emission lines, respectively (Akras et al., 2019).
Symbiotic stars are wide binary systems composed of a hot white dwarf and a
mass-losing red giant. Hydrodynamical studies suggest that some fraction of
slow stellar wind is gravitationally captured to form an accretion disk (de
Val-Borro et al., 2017; Chen et al., 2017; Saladino et al., 2018), where O
VI$\lambda\lambda$1032 and 1038 are important coolants.

Considering very small Raman scattering cross sections $\sim 10^{-23}{\rm\
cm^{2}}$ for O VI$\lambda\lambda$1032 and 1038 doublet lines, the operation of
Raman scattering requires a special condition that a very thick neutral region
is present in the vicinity of a strong O VI emission region. This special
condition is ideally met in symbiotic stars, where a thick neutral region
surrounding the giant component is illuminated by strong far UV radiation
originating from the nebular region that may be identified with the accretion
flow onto the hot component.

Raman-scattered O VI features in symbiotic stars exhibit complicated profiles
with multiple peaks separated by $\sim 30-50{\rm\ km\ s^{-1}}$ indicative of
the O VI emission regions with physical dimension of $\sim 1{\rm\ au}$ (e.g.
Shore et al., 2010; Heo & Lee, 2015; Lee et al., 2019). Because Raman and
Rayleigh scattering sufficiently off resonance shares the same scattering
phase function as Thomson scattering(Schmid, 1995; Yoo et al., 2002; Chang et
al., 2017), strong linear polarization may develop in an anisotropic
scattering geometry. Harries & Howarth (1996) conducted spectropolarimetric
observations of many symbiotic stars to show that Raman-scattered O VI
features are strongly polarized. They also found that most Raman-scattered O
VI features show polarization flip in the red wing part, where the
polarization develops nearly perpendicularly to the direction along which the
main part is polarized. Lee & Park (1999) proposed that the polarization flip
is closely associated with the bipolar structure of symbiotic stars (see also
Heo et al., 2016).

Raman scattering plays an interesting role of redistributing far UV radiation
near Ly$\beta$ and Ly$\gamma$ into near H$\alpha$ and H$\beta$, respectively.
He II being a single electron ion with $Z=2$, the transition to $n=2$ from an
energy level from $n=2k,k>1$ gives rise to emission lines with wavelengths
slightly shorter than those of H I Lyman series $k\to 1$, for which the cross
sections for Rayleigh and Raman scattering are conspicuously large. van
Groningen (1993) found Raman-scattered He II features near H$\beta$ in the
symbiotic nova RR Telescopii. Pequignot el al. (1997) discovered the same
spectral feature in the young planetary nebula NGC 7027, which constitutes the
first discovery of a spectral feature formed through Raman scattering with
atomic hydrogen in planetary nebulae. Subsequently, Raman-scattered He II at
$6545$ Å has been detected in several symbiotic stars (Birriel, 2004; Jung &
Lee, 2004; Sekeráš & Skopal, 2015) and in the young planetary nebulae NGC
6302, IC 5117, NGC 6790, NGC 6886, and NGC 6881 (Groves et al., 2002; Lee et
al, 2001; Lee et al., 2006; Kang et al., 2009; Choi & Lee, 2020).

It is particularly notable that the case B recombination theory allows one to
deduce the strengths of incident far UV He II lines (Storey & Hummer, 1995),
yielding faithful estimates of Raman conversion efficiencies. This is
extremely useful in the measurement of H I content in symbiotic stars and
young planetary nebulae, which makes Raman spectroscopy a totally new approach
to probing the mass loss processes occurring in the late stage of stellar
evolution (Lee et al., 2006; Choi et al., 2020). Furthermore, H$\alpha$ and
H$\beta$ in symbiotic stars and young planetary nebulae often display fairly
extended wings that may indicate the presence of fast tenuous stellar wind
(e.g., Arrieta & Torres-Peimbert, 2002). Broad wings around Balmer lines may
also arise via Raman scattering far-UV continuum around Lyman lines, which
requires further investigation (Lee, 2000; Yoo et al., 2002; Chang et al.,
2018).

Additional examples include Raman scattering of C II$\lambda\lambda$1036 and
1037 forming optical features at 7023 Å and 7054 Å, which were reported in the
symbiotic nova V1016 Cygni by Schild & Schmid (1996). Dopita et al. (2016)
investigated the H II regions in the Orion Nebula (M42) and five H II regions
in the Large and Small Magellanic Clouds to discover Raman-scattered features
at 6565 Å and 6480 Å formed through Raman scattering of O I $\lambda$1025.76
and Si II $\lambda$1023.70, respectively.

The Monte Carlo approach is an efficient numerical technique to describe
radiative transfer in various regions having a dust component(e.g. Celnikier &
Lefévre, 1974; Seon, 2015) and a molecular component(e.g. Brinch &
Hogerheijde, 2010). A similar approach has been applied to radiative transfer
of an electron scattering (e.g. Angel, 1969; Seon et al., 1994). Ly$\alpha$
deserves special attention being characterized by large scattering optical
depth (e.g. Eide et al., 2018; Seon & Kim, 2020).

Schmid (1992) performed Monte-Carlo simulations to investigate the formation
of Raman-scattered O VI features in an expanding H I region with an assumption
that a given line photon has an invariant scattering cross section as it
propagates through the H I medium. Similar studies were presented by Lee & Lee
(1997), who adopted a density matrix formalism to determine the physical
information of scattered radiation including polarization. Chang et al. (2015)
investigated the formation of broad Balmer wings near H$\alpha$ and H$\beta$
in the unification scheme of active galactic nuclei and presented
quantitatively the asymmetry of the wings formed in neutral regions with
extremely high H I column density $\sim 10^{23}\rm\ cm^{-2}$.

In this paper, we introduce a new grid-based Monte Carlo code entitled ”Sejong
Radiative Transfer for Rayleigh and Raman Scattering (STaRS), We also present
our test of the code by revisiting the formation of Balmer wings and Raman O
VI features through Raman scattering with atomic hydrogen.

##  2 Grid Based Monte Carlo Simulation

In this section, we describe ’STaRS’ and discuss the basic atomic physics of
Rayleigh and Raman scattering with atomic hydrogen. Fig. 1 is a schematic
illustration of a few representative transitions pertaining to the two types
of scattering. Thus far, detected Raman-scattered features are limited to
those associated with the final de-excitation to $2s$ state. The second order
time-dependent perturbation theory is used to compute the scattering cross
sections known as the Kramers-Heisenberg formula (e.g., Bethe & Salpeter,
1967; Sakurai, 1967; Saslow & Mills, 1969).

In our grid based Monte Carlo code, we divide the region of interest into a
large number of small cubes or cells in the Cartesian coordinate system, where
each cell is characterized by uniform physical properties in the three
dimensional space. Here, the uniform physical properties include H I number
density $n_{{\rm HI},G}$ and the velocity ${\bf v}_{G}$. We assign the
emissivity $j_{e}$ to each cell and generate an initial photon using
$j_{e}(\lambda_{i},x,y,z)$, where $\lambda_{i}$ is the wavelength of the
initial photon. STaRS is mainly written in Fortran with Message Passing
Interface implemented for parallel computing. We also adopt the shared memory
technique by intel MPI.

Fig. 2 shows a flow chart for radiative transfer simulations using STaRS. A
simulation starts with the setup of the scattering geometry by assigning to
each cell appropriate physical conditions. Initial far UV photons are
generated in accordance with our prescription of $j_{e}$. Each photon is
tagged with the information of the unit wavevector $\bf{\hat{k}}$, the
position vector $\bf r$, the wavelength $\lambda$, and the density matrix
$\rho$ composed of the Stokes parameters, $I$, $Q$, $U$, and $V$. The free
path $d$ for next scattering position $\bf r^{\prime}$ is computed by
transforming the scattering optical depth into the physical depth.

Decision is made whether the photon escapes from the region or is scattered
into another direction. If the next scattering position $\bf r^{\prime}$ is
outside the scattering region, the photon is assumed to reach the observer as
a far UV photon. Otherwise, we generate a new photon using the scattering
phase function for an electric dipole process and determine the scattering
type. We assume that the region is transparent to Raman-scattered photons. The
scattered photon escapes from the scattering region if the scattering is
Raman. Otherwise, the Rayleigh scattered photon is regarded as an incident
photon propagating to a new scattering position. The procedure is repeated
until escape.

In Sec. 2.1 and 2.2, we provide more detailed descriptions on the scattering
geometry and generation of initial far UV photon prescribed by $j_{e}$. In
Sec. 2.3, we describe the computation of a free path $d$. We discuss the basic
properties of the scattered photons in Sec. 2.4.

Figure 1: Schematic illustrations of energy levels and electronic transitions
associated with Raman and Rayleigh scattering with a hydrogen atom of far UV
electromagnetic radiation.  Figure 2: A flow chart of STaRS.

###  2.1 Geometry : Scattering Region

The medium for radiative transfer through Raman and Rayleigh scattering
corresponds to a thick H I region that is easily found in the slow stellar
wind from a red giant (Lee & Lee, 1997; Lee et al., 2019). No consideration on
the thermal motion of neutral hydrogen is given because the variation of the
cross section and the branching ratio in the scale of thermal speed is
negligible (Chang et al., 2015, 2018).

In our simulation, we divide the scattering region into a large number of
cells having a fixed size. Each cell is identified by a set of three indices
$(iX,iY,iZ)$ with $iX,iY$ and $iZ$ running from 1 to $N_{X},N_{Y}$ and
$N_{Z}$. If we denote by $X_{min}$ and $X_{max}$ the range of $x$ coordinates
for the scattering region, the boundary of the $iX$th cell on the $x$ axis is
given by

$X_{G}(iX)=X_{min}+\frac{X_{max}-X_{min}}{N_{X}}(iX-1).$ (1)

We also define $Y_{G}(iY)$ and $Z_{G}(iZ)$ in a similar way.

Therefore, for any point $P(x,y,z)$ in the cell identified with $(iX,iY,iZ)$,
we have

$\displaystyle X_{G}(iX)<x<X_{G}(iX+1)$ (2) $\displaystyle
Y_{G}(iY)<y<Y_{G}(iY+1)$ $\displaystyle Z_{G}(iZ)<z<Z_{G}(iZ+1).$

Each cell is also identified by its center point, whose coordinates are given
by the following relations

$\displaystyle X_{C}(iX)=\frac{X_{G}(iX)+X_{G}(iX+1)}{2}$ (3) $\displaystyle
Y_{C}(iY)=\frac{Y_{G}(iY)+Y_{G}(iY+1)}{2}$ $\displaystyle
Z_{C}(iZ)=\frac{Z_{G}(iZ)+Z_{G}(iZ+1)}{2}.$

###  2.2 Emission Source : Initial Photon

The operation of Raman scattering with atomic hydrogen requires the
coexistence of a strong far UV source and a thick neutral region. In the case
of symbiotic stars and young planetary nebulae, a strong far UV emission
region is formed near an accreting white dwarf or a hot central star and a
thick neutral region is also present in association with the mass loss of a
giant star. In the simulation, the emissivity $j_{e}(\lambda,x,y,z)$ is
prescribed as a function of wavelength and position. Regarding $j_{e}$ as the
normalized probability density function, we pick a wavelength $\lambda_{i}$
and a starting position $(x_{i},y_{i},z_{i})$ of an initial photon in
accordance with $j_{e}(\lambda,x,y,z)$.

We find the spatial index $(iX,iY,iZ)$ from $(x_{i},y_{i},z_{i})$ to determine
the cell that contains the starting position. For simplicity, we assume that
the initial photon is completely unpolarized and that the unit wavevector of
the initial photon is selected from an isotropic distribution. Thus, the
initial unit wavevector ${\bf\hat{k}}=(k_{x},k_{y},k_{z})$ is obtained using
two uniform random numbers $r_{1}$ and $r_{2}$ between 0 and 1 with the
following prescription

$\displaystyle\mu$ $\displaystyle=$ $\displaystyle\cos\theta=2r_{1}-1$ (4)
$\displaystyle\phi$ $\displaystyle=$ $\displaystyle 2\pi r_{2}$ $\displaystyle
k_{x}$ $\displaystyle=$ $\displaystyle\sin\theta\cos\phi$ $\displaystyle
k_{y}$ $\displaystyle=$ $\displaystyle\sin\theta\sin\phi$ $\displaystyle
k_{z}$ $\displaystyle=$ $\displaystyle\cos\theta.$

Because scattered radiation is polarized in an anisotropic geometry, it is
important to carry the polarization information. The four Stokes parameters
$(I,Q,U,V)$ are necessary to describe the polarization state. An equivalent
way is provided by considering the $2\times 2$ density matrix defined by

$\rho=\begin{bmatrix}{(I+Q)/2}&{(U+iV)/2}\\\ {(U-iV)/2}&{(I-Q)/2}\\\
\end{bmatrix}.$ (5)

.

In our simulation, the Stokes parameter $V$, representing circular
polarization, is always set to zero because no circular polarization develops
from initially unpolarizd photons in electric dipole processes associated with
Raman and Rayleigh scattering. Initially unpolarized photons are described by
a simple $\rho$ is given by

$\displaystyle\rho_{11}$ $\displaystyle=$ $\displaystyle 0.5$ (6)
$\displaystyle\rho_{22}$ $\displaystyle=$ $\displaystyle 0.5$
$\displaystyle\rho_{12}$ $\displaystyle=$ $\displaystyle\rho_{21}=0\ .$

The wavelength measured by an observer lying on the photon path is kept in the
simulation. We assume that the emitters are in random motion and also subject
to the bulk motion associated with the grid. If we let the velocity of the
emitter be ${\bf v}_{\rm emit}$ including bulk and random velocities, the
wavelength in the grid frame $\lambda_{g}$ is given by

$\lambda_{g}=\lambda_{i}\left(1-{\frac{{\bf v}_{\rm
emit}\cdot{\bf\hat{k}}}{c}}\right)\left(1+{\frac{{\bf v}_{\rm
G}\cdot{\bf\hat{k}}}{c}}\right),$ (7)

where ${\bf v}_{\rm G}$ is the velocity of the cell $(iX,iY,iZ)$.

###  2.3 Journey of Photon: Optical Depth and Free Path

Figure 3: The total scattering cross section(black), the branching ratio to
$n=2$ (blue) and $n=3$ (red) computed in Chang et al. (2015).

The generation of an initial photon is followed by the estimate of the free
optical depth $\tau$ given by

$\tau=-\ln r,$ (8)

where $r$ is a uniform random number between $0$ and $1$. We adopt the method
in Seon (2009) to compute the free path and the next scattering position in
grid-based geometry. In this section, we provide a brief description to
compute the free path and determine the scattering position.

If the starting position ${\bf r}$ of the photon is in cell $A$, we measure
the distance to the boundary of the cell $d_{A}$ along the photon ray from the
starting position. With $d_{A}$ we define the scattering optical depth
$\tau_{A}$ to the cell boundary by

$\tau_{A}=\sigma_{\rm tot}(\lambda_{A})n_{{\rm HI},A}d_{A},$ (9)

where $\lambda_{A}$ is the wavelength in cell $A$. Fig. 3 shows the total
scattering cross section $\sigma_{\rm tot}(\lambda)$ as a function of
wavelength. If $\tau_{A}>\tau$, then the next scattering position ${\bf
r}^{\prime}$ is found in cell $A$ as follows

${\bf r^{\prime}}={\bf
r}+\left(\frac{\tau}{\tau_{A}}\right)d_{A}{\bf\hat{k}}.$ (10)

In the opposite case where $\tau_{A}<\tau$, the photon enters the neighboring
cell $B$. In this case, the same problem is obtained if we regard the entry
point ${\bf r}_{e}={\bf r}+d_{A}{\bf\hat{k}}$ as the new starting point of the
photon with a new free optical path $\tau^{\prime}$ reduced by $\tau_{A}$, or

${\tau}^{\prime}={\tau}-{\tau_{A}}.$ (11)

It should be noted that we are dealing with radiative transfer in a medium in
motion. Therefore, cell $B$ may move with a velocity different from that of
cell $A$, in which case the photon wavelength along its propagation direction
may change on entering cell $B$ from cell $A$. Denoting by ${\bf v}_{\rm G,A}$
and ${\bf v}_{\rm G,B}$ the velocities of cells $A$ and $B$, respectively, we
have

$\lambda_{B}=\lambda_{A}\left(1-{\frac{{\bf v}_{\rm
G,A}\cdot{\bf\hat{k}}}{c}}\right)\left(1+\frac{{\bf v}_{\rm
G,B}\cdot{\bf\hat{k}}}{c}\right).$ (12)

Iterations from Eq. 9 are made with necessary updates ${\bf r}={\bf r}_{e}$,
$\tau=\tau^{\prime}$, and $\lambda_{A}=\lambda_{B}$ and new naming of cell $B$
as cell $A$ until we have $\tau_{A}>\tau$. Fig. 4 shows a schematic
illustration of this procedure. In cases when a neighboring cell $B$ may not
exist and $\bf r$ is outside the geometry, the photon is regarded as Rayleigh-
escaped.

Figure 4: Schematic illustration of computing the next scattering position in
the grid-based simulation.

###  2.4 Rayleigh and Raman Scattering

Raman scattering with atomic hydrogen with the initial and final electronic
states being $1s$ and $2s$ shares the same the scattering phase function with
Rayleigh scattering with atomic hydrogen in $1s$ (e.g. Chang et al., 2015). It
is also notable that the scattering is sufficiently far from resonance, the
scattering phase function is the same as that of Thomson scattering. In this
section, we describe the density matrix formalism with which the polarization
and the unit wavevector $\bf\hat{k}^{\prime}$ of the scattered photon is
determined (Ahn & Lee, 2015; Chang et al., 2017). We also describe the
wavelength conversion and the line broadening associated with Raman
scattering.

According to the density matrix formalism, the probability density of the
scattered wavevector
${\bf\hat{k}^{\prime}}=(\sin\theta^{\prime}\cos\phi^{\prime},\sin\theta^{\prime}\sin\phi^{\prime},\cos\theta^{\prime})$
is given by

$I^{\prime}({\bf\hat{k}^{\prime}})=\rho^{\prime}_{11}+\rho^{\prime}_{22},$
(13)

where $\rho^{\prime}_{ij}$ is defined by

$\rho^{\prime}_{ij}=\sum_{kl=1,2}(\hat{\epsilon}^{\prime}_{i}\cdot\hat{\epsilon_{k}})\rho_{kl}(\hat{\epsilon_{l}}\cdot\hat{\epsilon_{l}}^{\prime}).$
(14)

Here, $\hat{\epsilon}_{1,2}$ and $\hat{\epsilon}_{1,2}^{\prime}$ are the
polarization basis vector associated with $\bf\hat{k}$ and
$\bf\hat{k}^{\prime}$, respectively. Specifically,
$\hat{\epsilon}_{1}=(-\sin\phi,\cos\phi,0)$ and
$\hat{\epsilon}_{2}=(\cos\theta\cos\phi,\cos\theta\sin\phi,-\sin\theta)$ so
that $\hat{\epsilon}_{1}$ represents polarization in the direction
perpendicular to the plane spanned by the photon wavevector and the $z-$axis.

The components of the density matrix associated with the scattered radiation
are related to those of incident radiation by

$\displaystyle\rho^{\prime}_{11}$ $\displaystyle=$
$\displaystyle(\cos^{2}\Delta\phi)\rho_{11}$ $\displaystyle-$
$\displaystyle(\cos\theta\sin 2\Delta\phi)\rho_{12}$ $\displaystyle+$
$\displaystyle(\sin^{2}\Delta\phi\cos^{2}\theta)\rho_{22}$
$\displaystyle\rho^{\prime}_{12}$ $\displaystyle=$
$\displaystyle(\frac{1}{2}\cos\theta^{\prime}\sin 2\Delta\phi)\rho_{11}$
$\displaystyle+$ $\displaystyle(\cos\theta\cos\theta^{\prime}\cos
2\Delta\phi+\sin\theta\sin\theta^{\prime}\cos\Delta\phi)\rho_{12}$
$\displaystyle-$
$\displaystyle\cos\theta(\sin\theta\sin\theta^{\prime}\sin\Delta\phi+\frac{1}{2}\cos\theta\cos\theta^{\prime}\sin
2\Delta\phi)\rho_{22}$ $\displaystyle\rho^{\prime}_{22}$ $\displaystyle=$
$\displaystyle(\cos^{2}\theta^{\prime}\sin^{2}\Delta\phi)\rho_{11}$ (15)
$\displaystyle+$
$\displaystyle\cos\theta^{\prime}(2\sin\theta\sin\theta^{\prime}\sin\Delta\phi+\cos\theta\cos\theta^{\prime}\sin
2\Delta\phi)\rho_{12}$ $\displaystyle+$
$\displaystyle(\cos\theta\cos\theta^{\prime}\cos\Delta\phi+\sin\theta\sin\theta^{\prime})^{2}\rho_{22},$

where $\Delta\phi=\phi^{\prime}-\phi$. In the code, a selection is made for
${\bf\hat{k}}^{\prime}$ from an isotropic distribution to compute
$I^{\prime}$. A new random deviate $r^{\prime}$ is compared with $I^{\prime}$.
If $r^{\prime}<I^{\prime}$, then the selection of ${\bf\hat{k}}^{\prime}$ is
accepted. Otherwise, the process is iterated until the acceptance is obtained.
The scattered photon now becomes a new incident photon propagating to a new
scattering site. This is followed by necessary updates for the wavevector and
the density matrix given by

$\displaystyle{\bf\hat{k}}$ $\displaystyle=$
$\displaystyle{\bf\hat{k}^{\prime}}$ (16) $\displaystyle\rho_{ij}$
$\displaystyle=$ $\displaystyle\rho^{\prime}_{ij}.$

The scattering type is determined by the branching ratio. If an incident far
UV photon is more energetic than Ly$\gamma$, the final states available for
Raman scattering include $2s$, $3s$ and $3d$ states. For photons near
Ly$\gamma$, the total cross section is given by

$\sigma_{\rm tot}=\sigma_{1s}+\sigma_{2s}+\sigma_{3s+3d},$ (17)

where $\sigma_{1s}$, $\sigma_{2s}$, and $\sigma_{3s+3d}$ are the scattering
cross sections corresponding to the final states $1s$, $2s$, and $3s+3d$,
respectively. Rayleigh branching ratio $BR1$ is given by

$BR1=\frac{\sigma_{1s}}{\sigma_{\rm tot}}.$ (18)

In a similar way, Raman branching ratios corresponding to the final energy
levels $n=2$ and $3$ are

$BR2=\frac{\sigma_{2s}}{\sigma_{\rm tot}},\quad
BR3=\frac{\sigma_{3s+3d}}{\sigma_{\rm tot}},.$ (19)

In Fig. 3, the blue and red solid lines represent $BR2$ and $BR3$,
respectively.

The energy difference between the incident and Raman-scattered photons is the
same as that between the initial and final atomic states. This is translated
into the relation between the wavelengths $\lambda$ and $\lambda^{\prime}$ of
the Raman-scattered and incident photons, respectively, which is given by

$\frac{1}{\lambda}=\frac{1}{\lambda^{\prime}}+\frac{1}{\lambda_{\rm res}},$
(20)

where $\lambda_{\rm res}$ is the wavelength corresponding to the energy
difference between the initial and final states.

One very important aspect in Raman scattering can be found in the conspicuous
change in line width. Differentiating Eq. (20), we have

$\frac{d\lambda^{\prime}}{\lambda^{\prime}}=\left(\frac{\lambda^{\prime}}{\lambda}\right)\left(\frac{d\lambda}{\lambda}\right),$
(21)

from which we see immediately that the line width of Raman-scattered feature
is broadened by the factor $(\lambda^{\prime}/\lambda)^{2}$. For example, a
typically observed line width of Raman-scattered O VI at 6825 Åamounts to
$\sim 30$ Å whereas the far UV parent line O VI$\lambda\lambda$1032 exhibits a
line width $\sim 1$ Å in many symbiotic stars.

Due to the line broadening effect, far UV radiation near Lyman series of
hydrogen will be considerably diluted and redistributed around Balmer emission
lines to appear as broad wings (e.g. Yoo et al., 2002; Chang et al., 2015,
2018). Another important consequence of the line broadening effect is that the
line profiles of Raman-scattered features mainly reflect the relative motion
between the far UV emitters and the neutral scatterers and quite independent
of the observer’s line of sight (Heo et al., 2016; Choi et al., 2020).

##  3 Code Test

As a check of our code, we present our simulation results for two exemplary
cases. The first example is a static spherical H I region surrounding a far UV
continuum source located at the center, in which Balmer wings are formed
through Raman scattering. Analytic solutions are available for this case,
against which our result obtained from STaRS is compared. The second case is
reproduction of the result of Lee & Lee (1997), who investigated Raman
scattering of O VI in symbiotic stars. In this case, the H I region is an
expanding spherical wind around the giant component. Schematic illustrations
of the two cases are shown in Fig. 5.

###  3.1 Formation of Balmer Wings in a Static Spherical H I Region

The central point-like far UV source surrounded by a spherical H I region with
radius $R$ is characterized by a flat continuum. The static neutral region is
assumed to be of uniform H I density $n_{\rm HI}$. We fix the radial column
density $N_{\rm HI}$ defined by

$N_{\rm HI}=n_{\rm HI}R=10^{23}{\rm\ cm^{-2}},$ (22)

and vary the number of cells. We set $N_{xyz}=N_{x}=N_{y}=N_{z}$, so that the
total number of cells is given by $N_{xyz}^{3}$. We generate $10^{7}$ photons
for each simulation. The initial photons are generated at the center of the H
I sphere in accordance with the emisivity $j_{e}$ given by the three
dimensional Dirac delta function

$j_{e}(\lambda,x,y,z)={\lambda\over{\lambda_{max}-\lambda_{min}}}\delta^{(3)}({\bf
r}),$ (23)

where $\lambda_{max}$ and $\lambda_{min}$ are the maximum and minmum
wavelengths of initial photons.

Fig. 6 shows optical spectra formed through Raman scattering for the cases of
$N_{xyz}=3,10$ and 100. The vertical axis shows the Raman conversion
efficiency ($RCE$) defined as the number ratio per unit wavelength of incident
far UV photons and Raman-scattered optical photons. The left and right panels
show Balmer wings formed through Raman scattering around H$\alpha$ and
H$\beta$, respectively. The solid lines represent the analytic solutions and
open circles show simulation results obtained using STaRS. The analytic
solutions are obtained from the non-grid based simulation in Chang et al.
(2015, 2018).

Far from the line centers, the simulation results for $N_{xyz}=3$ are slightly
higher than the analytic solutions. Other than this, the agreement is fairly
good, indicating little dependence on $N_{xyz}$ in the case of a static H I
region. The volume of the H I region with $N_{xyz}=3$ is larger than the
sphere with the radius $R$ as the H I density associated with a cell is
determined by the central position $(X_{C},Y_{C},Z_{C})$ of the cell. When
$(X_{C}^{2}+Y_{C}^{2}+Z_{C}^{2})^{1/2}$ is smaller than $R$, the H I density
of the cell is assigned to be $n_{HI}$. Otherwise, the density is assigned to
be zero. The number ratio between Raman scattered photons near H$\alpha$ and
the total initial photons near Ly$\beta$ is 20.77 % for the analytic solution,
whereas the simulations give 21.17, 20.91, and 20.77 % for $N_{xyz}=3$, 10,
and 100, respectively.

In Fig. 7, we present the polarimetric data of Raman scattered H$\alpha$
projected to the celestial sphere. The surface brightness, the degree of
polarization $p$, and the direction of polarization are shown in the left,
middle and the right panels, respectively. Here, in terms of the Stokes
parameters, the degree of polarization $p$ and the position angle $\phi_{p}$
are given by

$p=\frac{\sqrt{Q^{2}+U^{2}}}{I},\quad\phi_{p}=\frac{1}{2}\tan^{-1}\left(\frac{U}{Q}\right).$
(24)

In the left panels, the surface brightness is shown in logarithmic scale. With
our choice of a rather large value of $N_{\rm HI}$, a considerable fraction of
Raman-scattered photons are formed near the source, which leads to excellent
agreement between the analytic results and those obtained using STaRS. One may
notice increase in smoothness of the surface brightness as $N_{xyz}$
increases. It is also noticeable that $p$ becomes large with increasing
distance from the center. The concentric polarization pattern reflects the
spherically symmetric scattering geometry. In this particular case,
$N_{xyz}>10$ appears to be sufficient to describe the the analytic result of
the static medium.

Figure 5: Schematic illustrations of two cases to test the code.  Figure 6:
Raman conversion efficiency of H$\alpha$ (left panel) and H$\beta$ (right
panel) broad wings of the static spherical case. The black solid lines are the
spectra by the analytic method. The $x$ axis is the wavelength. The red, blue,
and green open circles are the spectra by STaRS for $N_{xyz}=3$, 10, and 100.

###  3.2 Raman Scattering of O VI in Expanding H I Region

Lee & Lee (1997) presented their basic study of line formation of Raman O VI
in a symbiotic star consisting of a white dwarf and a mass losing giant. The O
VI emission region near the white dwarf component was assumed to be a point
source broadened thermally with $T=10^{4}\,\rm K$. For simplicity, the slow
stellar wind from the red giant component is assumed to be entirely neutral
ignoring the photoionization by the white dwarf. The orbital separation is
$10R_{*}$ of which $R_{*}$ is the radius of the red giant. The positions of
the white dwarf and the red giant are $(0,0,10R_{*})$ and $(0,0,0)$,
respectively. The initial photons are generated at the position of the white
dwarf according to the emissivity $j_{e}$ given by

$\displaystyle j_{e}(\lambda,x,y,z)=$ (25)
$\displaystyle{1\over{\sigma_{th}\sqrt{2\pi}}}\exp\left[{(\lambda-\lambda_{1032})^{2}\over{2\sigma_{th}^{2}}}\right]\delta^{(3)}(x,y,z-10R_{*})$

where $\lambda_{1032}$ is the center wavelength of O VI $\lambda$ 1032 and
$\sigma_{th}\sim 0.008{\rm\ \AA}$ is the thermal width of O VI $\lambda$ 1032
with $T=10^{4}K$.

The velocity $\bf v(r)$ and H I number density $n({\bf r})$ are the functions
of the distance from the red giant $r=|\bf r|$ are given by

$\displaystyle{\bf v}({\bf r})$ $\displaystyle=$ $\displaystyle
v_{\infty}(1-R_{*}/r)\frac{\bf r}{r}$ (26) $\displaystyle\quad n({\bf r})$
$\displaystyle=$ $\displaystyle n_{0}(R_{*}/r)^{2}(1-R_{*}/r)^{-1},$

where $n_{0}$ is the characteristic number density defined in Eq. 2.11 of Lee
& Lee (1997) and $v_{\infty}$ is the terminal velocity. As a second code test
case, we revisit Raman O VI formation illustrated in Figs. 1 and 6 of Lee &
Lee (1997). Fixing $v_{\infty}=20\,\rm km\,s^{-1}$, we consider the three
values of $\tau_{0}=0.5$, 1, and 10, where $\tau_{0}=n_{0}R_{*}\sigma_{\rm
tot}$. We generate $10^{9}$ photons for each simulation.

Fig. 8 shows the spectra and the Stokes parameter $P=Q/I$ for Raman scattered
O VI at 6825 Å. It should be noted that the integrated $U$ vanishes due to the
axial symmetry about the $z$ axis. Therefore, the signed ratio $Q/I$
represents the degree and direction of polarization, where a positive and a
negative $Q$ correspond to the polarization in direction perpendicular and
parallel to the $z$ axis, respectively. The observer’s line of sight lies in
$x$-$y$ plane. Noting that it is perpendicular to the symmetry axis, maxium
polarization can be developed along this direction or the symmetry $z-$axis.
We collecte the photons escaping toward the observer. The fraction of the
collected Raman photons near 6825 Å is 0.11, 0.21, and 1.82 % for
$\tau_{0}=0.5$, 1, and 10, respectively.

The line profiles obtained using STaRS differ slightly from those Lee & Lee
(1997) presented. In particular, the red peaks obtained from STaRS are more
enhanced than those presented by Lee & Lee (1997), which is attributed to the
geometrical truncation adopted by them. The difference is in the range of
10-15 per cent with respect to the red peak. We find that the agreement gets
better where the full range of the scattering region is taken into account. In
the bottom panels, we find overall agreement in the polarization behaviors.
The noisy features with $\Delta V>200\rm\,km\,s^{-1}$ are attributed to the
small number statistics of collected photons. Blue photons are scattered
mostly in the compact and dense region between the red giant and white dwarf,
resulting in development of strong polarization in the direction perpendicular
to the $z$ axis. In contrast, red photons are scattered mainly in a quite
extended region near the red giant, leading to weak polarization and enhanced
line flux.

Fair agreement shown in Fig. 6 and 8 demonstrates that STaRS has been well-
tested. Furthermore, as illustrated in Fig. 7, STaRS is capable of study of
radiative transfer for spectropolarimetric imaging observations.

Figure 7: The surface brightness (left), the degree of polarization (center),
and the direction of polarization (right) of the projected H$\alpha$ photons.
The panels in the first, second, and third rows represent the results obtained
using STaRS for $N_{xyz}=3$, 10, and 100. The bottom panels represent the
result obtained using an analytic method.  Figure 8: Line formation of Raman O
VI 6825 features in an expanding H I region. The blue solid lines are computed
by ’STaRS’. The red dashed lines are the results of Fig. 6 in Lee & Lee
(1997).

##  4 Summary

We have developed a 3D grid-based Monte Carlo code ’STaRS’ for radiative
transfer through Raman and Rayleigh scattering, which can be mainly used to
investigate line formation of Raman-scattered features in a thick neutral
region illuminated by a strong far UV emission source. Favorable conditions
for Raman scattering with atomic hydrogen are easily met in symbiotic stars,
young planetary nebulae and active galactic nuclei. Through a couple of tests,
we have successfully demonstrated that ’STaRS’ is a flexible code to deal with
radiative transfer in a thick neutral media yielding multidimensional
spectropolarimetric and imaging data. ’STaRS’ is easily accessed in Github
’https://github.com/csj607/STaRS’.

###### Acknowledgements.

This research was supported by the Korea Astronomy and Space Science Institute
under the R&D program (Project No. 2018-1-860-00) super-vised by the Ministry
of Science, ICT and Future Planning. This work was also supported by a
National Research Foundation of Korea (NRF) grant funded by the Korea
government (MSIT; No. NRF-2018R1D1A1B07043944). Seok-Jun is very grateful to
Dr. Kwang-Il Seon for his help in adoption of grid-based technique.

## References

  * Ahn & Lee (2015) Ahn, S.-H., Lee, H.-W., 2015, Polarization of Lyman $\alpha$ Emergent from a Thick Slab of Neutral Hydrogen, JKAS, 48, 195
  * Akras et al. (2019) Akras, S., Guzman-Ramirez, L., Leal-Ferreira, M. L., Ramos-Larios, G., 2019, A Census of Symbiotic Stars in the 2MASS, WISE, and Gaia Surveys, ApJS, 240, 21
  * Angel (1969) Angel, J. R. P., 1969, Polarization of Thermal X-Ray Sources, ApJ, 158, 219
  * Arrieta & Torres-Peimbert (2002) Arrieta, A., Torres-Peimbert, S., 2002, Broad H$\alpha$ Wings in Young Planetary Nebulae RMxAC, 12, 154
  * Bethe & Salpeter (1967) Bethe, H. A., Salpeter, E. E., 1967, Quantum Mechanics of One and Two Electron Atoms, Academic Press, New York
  * Brinch & Hogerheijde (2010) Brinch, C., Hogerheijde, M. R., 2010, LIME – a flexible, non-LTE line excitation and radiation transfer method for millimeter and far-infrared wavelengths, A&A, 523, A25
  * Birriel (2004) Birriel, J. J., 2004, Raman-scattered He II at 4851 Åin the symbiotic stars HM Sagittae and V1016 Cygni, ApJ, 612, 1136
  * Celnikier & Lefévre (1974) Celnikier, L. M., Lefévre, J., 1974, Radiation transport in circumstellar dust – A Monte Carlo approach, A&A, 36, 429
  * Chang et al. (2015) Chang, S.-J., Heo, J.-E., Di Mille, F., Angeloni, R., Palma, T., Lee, H.-W., 2015, Formation of Raman scattering wings around H$\alpha$, H$\beta$, and Pa$\alpha$ in active galactic nuclei, ApJ, 814, 98
  * Chang et al. (2017) Chang, S.-J., Lee, H.-W., Yang, Y., Polarization of Rayleigh scattered Ly$\alpha$ in active galactic nuclei, 2017, MNRAS, 464, 5018
  * Chang et al. (2018) Chang, S.-J., Lee, H.-W., Lee, H.-G., Hwang, N., Ahn, S.-H., Park, B.-G., 2018, Broad Wings around H$\alpha$ and H$\beta$ in the Two S-type Symbiotic Stars Z Andromedae and AG Draconis, ApJ, 866, 129
  * Chen et al. (2017) Chen, Z., Frank, A., Blackman, E. G., Nordhaus, J., Carroll-Nellenback, J., 2017, Mass transfer and disc formation in AGB binary systems, MNRAS, 468, 4465
  * Choi et al. (2020) Choi, B.-E., Chang, S.-J., Lee, H.-G., Lee, H.-W., 2020, Line Formation of Raman-scattered He II $\lambda$4851 in an Expanding Spherical H I Shell in Young Planetary Nebulae, ApJ, 889, 2
  * Choi & Lee (2020) Choi, B.-E., Lee, H.-W., 2020, Discovery of Raman-scattered He II $\lambda$6545 in the Planetary Nebulae NGC 6886 and NGC 6881, ApJL, 903L, 39
  * de Val-Borro et al. (2017) de Val-Borro, M., Karovska, M., Sasselov, D. D., Stone, J. M., 2017, Three-dimensional hydrodynamical models of wind and outburst-related accretion in symbiotic systems, MNRAS, 468, 3408
  * Dopita et al. (2016) Dopita, M. A., Nicholls, D. C., Sutherland, R. S., Kewley, L. J., Groves, B. A., 2016, The discovery of Raman scattering in H II regions, ApJL, 824, L13
  * Eide et al. (2018) Eide, M. B., Gronke, M., Dijkstra, M., Hayes, M., 2018, ApJ, Unlocking the Full Potential of Extragalactic Ly$\alpha$ through Its Polarization Properties, 856, 156
  * Groves et al. (2002) Groves, B., Dopita, M. A., Williams, R. E., Hua, C.-T., 2002, The internal extinction curve of NGC 6302 and its extraordinary spectrum, PASA, 19, 425
  * Harries & Howarth (1996) Harries, T. J., Howarth, I. D., 1996, Raman scattering in symbiotic stars. I. Spectropolarimetric observations, A&AS, 119, 61
  * Heo & Lee (2015) Heo, J.-E., Lee, H.-W., 2015, Accretion flow and disparate profiles of Raman scattered O VI $\lambda\lambda$1032, 1038 in the symbitic star V1016 Cygni, JKAS, 48, 105
  * Heo et al. (2016) Heo, J.-E., Angeloni, R., Di Mille, F., Palma, T., Lee, H.-W., 2016, A Profile Analysis of Raman-scattered O VI Bands at 6825 Å and 7082 Å in Sanduleak’s Star, ApJ, 833, 286
  * Jung & Lee (2004) Jung, Y.-C., Lee, H.-W., 2004, Centre shift of the Raman scattered He II $\lambda$4850 in the symbiotic star V1016 Cygni, MNRAS, 355, 221
  * Kang et al. (2009) Kang, E.-H., Lee, B.-C., Lee, H.-W., 2009, Raman-Scattered He II $\lambda$6545 in the Young and Compact Planetary Nebula NGC 6790, ApJ, 695, 542
  * Lee & Park (1999) Lee, H.-W., Park. M.-G., 1999, Toward the Evidence of the Accretion Disk Emission in the Symbiotic Star RR Telescopii, ApJL, 515, L89
  * Lee (2000) Lee, H.-W., 2000, Raman-Scattering Wings of H$\alpha$ in Symbiotic Stars, ApJ, 541, 25
  * Lee et al (2001) Lee, H.-W., Kang, Y.-W., Byun, Y.-I., 2001, Raman-scattered He II Line in the Planetary Nebula M2-9 and in the Symbiotic Stars RR Telescopii and He 2-106, ApJ, 551, 121
  * Lee et al. (2006) Lee, H.-W., Jung, Y.-C., Song, I.-O., Ahn, S.-H., 2006, Raman-scattered He II $\lambda\lambda$4850, 6545 in the Young and Compact Planetary Nebula IC 5117, ApJ, 636, 1045
  * Lee & Lee (1997) Lee, K.-W., Lee, H.-W., 1997, On the profiles and the polarization of Raman-scattered emission lines in symbiotic stars – II. Numerical simulations, MNRAS, 292, 573
  * Lee et al. (2019) Lee, Y.-M., Lee, H.-W., Lee, H.-G., Angeloni, R., 2019, Stellar-wind accretion and Raman-scattered O VI features in the symbiotic star AG Draconis, MNRAS, 487, 2166
  * Miller & Goodrich (1990) Miller, J. S., Goodrich, R. W., 1990, Spectropolarimetry of high-polarization Seyfert 2 galaxies and unified Seyfert theories, ApJ, 355, 456
  * Nussbaumer et al. (1989) Nussbaumer, H., Schmid, H. M., Vegel, M., 1989, Raman scattering as a diagnostic possibility in astrophysics, A&A, 211L, 27
  * Pequignot el al. (1997) Pequignot, D., Baluteau, J.-P., Morisset, C., Boisson, C., 1997, NGC7027: Discovery of a Raman line in a planetary nebula, A&A, 323, 217
  * Saladino et al. (2018) Saladino, M. I., Pols, O. R., van der Helm, E., Pelupessy, I., Portegies Zwart, S., 2018, Gone with the wind: the impact of wind mass transfer on the orbital evolution of AGB binary systems, A&A, 618, A50
  * Saslow & Mills (1969) Saslow, W. M., Mills, D. L., 1969, Raman scattering by hydrogenic systems, Phys. Rev., 187, 1025
  * Sakurai (1967) Sakurai, J. J., 1967, Advanced Quantum Mechanics, Addison-Wesley, Reading, MA.
  * Schild & Schmid (1996) Schild, H., Schmid, H. M., 1996, Spectropolarimetry of symbiotic stars. On the binary orbit and the geometric structure of V1016 Cygni, A&A, 310, 211
  * Schmid (1989) Schmid, H. M., 1989, Identification of the emission bands at 6830, 7088 Å, A&A, 211L, 31
  * Schmid (1992) Schmid, H. M., 1992, Monte-Carlo simulations of Raman scattered O VI emission lines in symbiotic stars, A&A, 254, 224
  * Schmid (1995) Schmid, H. M., 1995, Monte Carlo simulations of the Rayleigh scattering effects in symbiotic stars, MNRAS, 275, 227
  * Seon et al. (1994) Seon, K.-I., Min, K. W., Choi, C. S., Nam, U. W., 1994, Monte Carlo simulation of comptonization in a spherical shell geometry, JKAS, 27, 45
  * Seon (2009) Seon, K.-I., 2009, Monte-Carlo Simulation of the Dust Scattering, PKAS, 24, 43
  * Seon (2015) Seon, K.-I., 2015, Monte-Carlo Radiative Transfer Model of the Diffuse Galactic Light, JKAS, 47, 57
  * Seon & Kim (2020) Seon, K.-I., Kim, C.-G., 2020, Ly-alpha Radiative Transfer: Monte-Carlo Simulation of the Wouthuysen-Field Effect, ApJS, 250, 9
  * Sekeráš & Skopal (2015) Sekeráš, M., Skopal, A., 2015, Mass-loss Rate by the Mira in the Symbiotic Binary V1016 Cygni from Raman Scattering, ApJ, 812, 162
  * Shore et al. (2010) Shore, S. N., Wahlgren, G. M., Genovali1, K., Bernabei, S., Koubsky, P., 2010, The spectroscopic evolution of the symbiotic star AG Draconis: I. The O VI Raman, Balmer, and helium emission line variations during the outburst of 2006-2008, A&A, 510, A70
  * Storey & Hummer (1995) Storey, P. J., Hummer, D. G., 1995, Recombination line intensities for hydrogenic ions – IV. Total recombination coefficients and machine-readable tables for $Z=1$ to 8, MNRAS, 272, 41
  * Tran (2010) Tran, H. D., 2010, Hidden Double-peaked Emitters in Seyfert 2 Galaxies, ApJ, 711, 1174
  * van Groningen (1993) van Gronigen, E., 1993, Further evidence for Raman scattering in RR Tel, MNRAS, 264, 975
  * Yoo et al. (2002) Yoo, J. J., Bak, J.-Y., Lee, H.-W., 2002, Polarization of the broad H$\alpha$ wing in symbiotic stars, MNRAS, 336, 467


# Textual Explanations for Automated Commentary Driving
††thanks: This work was partially supported by a fellowship within the IFI
programme of the German Academic Exchange Service (DAAD) and by the EPSRC
project RAILS (grant reference: EP/W011344/1).

Marc Alexander Kühn§ UVC Partners
Munich, Germany
<EMAIL_ADDRESS>Daniel Omeiza  Oxford Robotics Institute
Department of Engineering Science
University of Oxford
Oxford, United Kingdom
<EMAIL_ADDRESS>Lars Kunze  Oxford Robotics Institute
Department of Engineering Science
University of Oxford
Oxford, United Kingdom
<EMAIL_ADDRESS>
###### Abstract

The provision of natural language explanations for the predictions of deep-
learning-based vehicle controllers is critical as it enhances transparency and
easy audit. In this work, a state-of-the-art (SOTA) prediction and explanation
model is thoroughly evaluated and validated (as a benchmark) on the new
Sense–Assess–eXplain (SAX). Additionally, we developed a new explainer model
that improved over the baseline architecture in two ways: (i) an integration
of part of speech prediction and (ii) an introduction of special token
penalties. On the BLEU metric, our explanation generation technique
outperformed SOTA by a factor of 7.7 when applied on the BDD-X dataset. The
description generation technique is also improved by a factor of 1.3. Hence,
our work contributes to the realisation of future explainable autonomous
vehicles.

###### Index Terms:

advanced driver assistance systems, automated vehicles, natural language
explanations, deep learning

§§footnotetext: During the project, the author was with the Department of
Informatics, Technical University of Munich.

©2022 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

##  I Introduction

Deep neural networks typically operate as black-boxes without providing
adequate insights into their decision-making process. This makes the
justification of their usage in high-stakes scenarios difficult [1].
Explainable models are handy as they provide better insights into the internal
workings of the underlying opaque models. While the end-to-end learned model
in [1] reports state-of-the-art results, it has only been tested on a single
dataset that incorporates only post-hoc explanations in free text form. Hence,
a thorough evaluation of the language generation model on additional datasets
is essential.

Figure 1: On a high level, the approach takes a series of image frames as
input into an opaque neural network and outputs a natural language
explanation.

In this paper, we first use a novel _Explanation Driving Dataset_ 111Contact
authors for access. which was collected during the Sense–Assess–eXplain (SAX)
project [2]. The dataset incorporates real-time driver audio commentary into
its annotation corpus which allows learning a model on more realistic
explanations. Moreover, the natural language annotations were generated and
presented in a more structured form which was guided by a multi-label visual
annotation process. We evaluated the SOTA model [1] on the new dataset
(Section IV). Also, we implemented and evaluated a new architecture similar to
[1] but with two modifications; (i) an integration of part of speech
prediction and (ii) an introduction of special token penalties. Overall, this
paper makes the following contributions:

  1. 1.

an evaluation of the state-of-the-art pipeline on the new SAX data.

  2. 2.

integration of part of speech prediction in the language generator.

  3. 3.

an improved sentence structure via the inclusion of special token penalties
into the loss function.

Fig. 1 shows the high-level process overview. Our code is accessible via
GitHub222https://github.com/cognitive-robots/causal-natural-language-
explanations.

##  II Related Work

###  II-A Explainable AI for Self-Driving Vehicles

It is argued that explainable AI in autonomous driving holds many benefits,
among which are the enhancement of accountability and facilitation of trust
[3, 4]. One commonly used categorisation of explainable AI techniques is based
on the explanation provision timing—that is, model intrinsic (inherently
interpretable), and post-hoc (built upon a black-box model to explain
predictions). Interpretable driver controllers have been proposed in previous
works. For example, real-time visual attention maps which highlight causally
influential regions in driving image frames were proposed in [5], [6]. The
architecture pipeline consists of a visual attention-based CNN with a
subsequent causal filtering technique. Kim et al. [1] imported a human-created
knowledge corpus consisting of free-text action descriptions and explanations
into an attention-based deep learning vehicle controller to provide natural
language action descriptions and explanations based on a series of image
frames. This approach can also be referred to as an image captioning method.
Image captioning is the task to generate a sentence explaining an input image
[7]. Most recent methods are based on a CNN- and RNN-based encoder-decoder
model like presented by [8] in an unsupervised set-up, by [9] incorporating
scene graphs or by [10] using a hierarchical attention network. While those
publications use domain-unspecific general image data, [1] and [7] showed that
image captioning methods trained on specific driving data can create textual
explanations for driving decisions. In a related work [11], a knowledge corpus
consisting of natural language human advice to the driver/AI controller was
added. A tree-based approach for automated driving commentary generation has
also been explored in [12] for select driving actions, such as stop, move, and
lane changes.

###  II-B Berkeley DeepDrive eXplanation Dataset

The Berkeley DeepDrive eXplanation (BDD-X) dataset [1] consists of videos
captured by a dashcam. In addition to the video data, the set also contains
timestamped sensor measurements like the vehicle’s velocity, course and GPS
location. For each scene, the textual annotation consists of an action
description and an explanation part with the corresponding timestamps. An
example annotation is: "The car is moving into the right lane because it is
safe to do so.". The annotator writes the description and the explanation part
separately in two individual free-text boxes. The dataset consists of over 77
hours of driving scenes within 6,984 videos with an average of 3.8 high-level
driving actions. The vocabulary of training action descriptions is 906 words
and of explanations 1,668 words. All annotations are conducted post-hoc such
that no insights into the driver’s real decision-making process are provided.

###  II-C SAX Explanation Driving Dataset

One key contribution of this work is the evaluation of the state-of-the-art
natural language explanation model with the new Sense-Assess-eXplain (SAX)
explanation driving dataset. The dataset was obtained during the SAX project
[2] and it extends existing explanation datasets like BDD-X to provide more
structured information through real-time commentary driving. It contains 9.5
hours of driving scenes from London roads with real-time driver audio
commentary. The explanation subset contains camera, GPS and CAN-bus data and
is visually annotated with the agents’ type, action and position. In addition,
agents are labelled to indicate how much of an influence they were on the
driver’s action. This was jointly based on the annotators’ judgments and the
commentary provided by the driver. Annotators were familiar with the UK road
rules. After the completion of the visual annotation, the labels are parsed
and agents that had influences are selected and structured to make ground
truth explanations for the corresponding ego vehicle’s actions. This creates
an annotation corpus with standardised sentences and allows to have insights
into the real decision-making process of the driver. An example annotation is:
"Car is stopping because pedestrian is crossing on ego’s lane". Most action
sequences are circa 3 seconds long. In total, the dataset consists of 491
driving sequences and the used vocabulary consists of 34 words.

##  III Prediction and Explanation Model

We implemented the SOTA attention-based deep learning model proposed in [1].
The model predicts high-level driving actions, such as moving, stopping,
changing lanes, and turning. The explanations consist of attention maps
accompanied by generated natural language texts. A description of the model is
provided in Fig. 2. The model is designed to obtain video frames as input and
generates natural language action descriptions and explanations as output. In
the training, vehicle acceleration and vehicle course values are fed to the
model along with corresponding ground truth natural language annotations.
Structurally, the pipeline consists of a CNN-based visual feature encoder that
extracts visual information which represent high-level object descriptions.
Afterwards, an LSTM-based vehicle controller with implemented spatial
attention predicts vehicle acceleration and change of course for a series of
image frames. In the end, an LSTM-based language generator uses the
controller’s spatial attention together with temporal attention over multiple
video frames. The textual action descriptions and explanations are generated
by outputting per-word softmax probabilities.

Figure 2: The pipeline consists of a visual feature encoder, a vehicle
controller that uses spatial attention and a language generator that uses
strongly aligned spatial attention plus temporal attention. We extended the
standard model with PoS prediction and special token penalties.

In order to achieve higher-quality textual action descriptions and
explanations, we modified the language generator in [1] with the aim of
generating better sentence grammar. Further, we integrated part of speech
prediction in the generator and included special token penalties in the loss
function.

###  III-A Integration of Part of Speech Prediction

The idea is to include an explicit understanding of part of speech (PoS) into
the model. In theory, having constraints on the part of speech could benefit
the sentence grammar and improve the similarity score between the generated
and ground truth sentences. First, the PoS needs to be added to the natural
language annotation corpus for both datasets. For that, each caption is fed
into the pre-trained part-of-speech-tagger by [13] using the universal tagset.
The tagger processes sequences of words and adds PoS tags to each word. The
used "universal tagset" categorizes adjectives, adpositions, adverbs,
conjunctions, determiners, nouns, numerals, particles, pronouns, verbs and
punctuation marks. The tokens <sep>, <START>, <END> and <NULL> are replaced
with ";" tokens such that they get tagged uniformly as punctuation marks.

During training, the deep learning model needs to be adjusted to predict the
PoS for each word. This involves the implementation of an additional decoder
layer that takes as input the LSTM hidden state $\bm{h}_{k}$ and temporal
context vector $\bm{z}_{k}$ and produces part-of-speech predictions
$\bm{o}_{k}^{pos}$ for each word according to Equation 1. W refers to
trainable weight matrices. These predictions are an additional input to the
subsequent decoder layer that produces the word predictions $\bm{o}_{k}$.

$\bm{o}_{k}^{pos}=\operatorname{softmax}\left(\left(\mathbf{h}_{k}\cdot\mathbf{W}^{h\\_pos}+\mathbf{z}_{k}\mathbf{W}^{z\\_pos}\right)\\\
\cdot\mathbf{W}^{gen\\_out\\_pos}\right)$ (1)

The loss function of the language generator is additively extended with a new
cross entropy loss that takes as input the ground truth one-hot encoded PoS
data $p_{k}^{pos}(x)$ together with the generated PoS outputs
$\bm{o}_{k}^{pos}$ according to Equation 2. $k$ is an iterator over each word
of the sentence and $T$ describes the tagset containing all possible PoS tags.

$\mathcal{L}^{pos}_{\mathrm{g}}=-\sum_{k}\sum_{x\in
T}p_{k}^{pos}\left(x\right)\log\left(\bm{o}_{k}^{pos}\left(x\right)\right)$
(2)

The overall generator loss with PoS prediction is constructed out of the new
PoS loss $\mathcal{L}^{pos}_{\mathrm{g}}$ weighted by $\lambda_{pos}$ and the
original generator loss $\mathcal{L}_{\mathrm{g}}$ weighted by
$(1-\lambda_{pos})$.

###  III-B Integration of Special Token Penalties

The raw caption data used in this project has a unique syntax consisting of a
<START> token in the beginning, a <sep> token between action description and
explanation, an <END> token in the sentence end and possibly several <NULL>
tokens as padding towards the maximum sentence length. An example caption
looks like this: "<START> the car accelerates to a constant speed <sep>
because the light has turned green <END> <NULL> <NULL> <NULL> <NULL> <NULL>
<NULL>". To steer the model towards generating the correct syntax, deviations
from an ideal sentence structure should be penalised.

As a first step, each generated <NULL> token will be penalised with the
penalty $\gamma_{null}$. This should steer the model towards trying to predict
a word instead of a <NULL> token in case of uncertainty towards the end of the
sentence. Therefore, the sentence lengths of generated and ground truth
sentences should be further aligned. The penalties are additively included in
the model loss function. Secondly, each deviation from having the tokens
<START>, <sep>, <END> exactly once per sentence should be penalised as well.
In this process, the mentioned tokens are counted in each generated caption.
In the case of $\\#token\neq 1$, a penalty $\gamma_{other}$ once per
respective token is added to the loss function.

##  IV Experiments and Evaluation

###  IV-A Training Procedure and Evaluation Set-Up

For the training of the language generator, an exponential learning rate decay
with a decay rate of 0.96 was implemented together with an Adam optimizer
[14]. For the BDD-X dataset, the decay was applied every 11,600 steps and for
the SAX dataset every 3,500 steps. For the BDD-X dataset, a start learning
rate of $1e^{-5}$, a batch size of 32 and a dictionary with 1,300 words were
used. For the SAX dataset, a start learning rate of $6e^{-7}$, a batch size of
16 and a dictionary with 34 words were used. Furthermore, TensorFlow version
1.15 was used and the graph-level random seed for the default graph was fixed
for better reproducibility [15]. Moreover, an NVIDIA GeForce RTX 2080 Ti with
11GB GDDR6 memory was used. The datasets were split between training,
validation and test sets with ratios of 80%/10%/10% for BDD-X and
75%/12.5%/12.5% for SAX, respectively.

Two quantitative performance metrics were used: METEOR and BLEU. The METEOR
metric was implemented following the algorithm in [16] and using the Natural
Language Toolkit from [13]. METEOR functions by using word-to-word matches
between a reference and candidate string which create word alignment. The
METEOR score is then calculated based on a parameterized harmonic mean of
unigram precision and unigram recall taking also into account a penalty
referring to the correct word order. For the BLEU metric, SacreBleu [17] which
wraps the original reference implementation in [18] with additional features
to provide comparable corpus-level scores was used. The main idea of BLEU is
to compare n-grams of the candidate and reference string by counting the
number of position-independent matches. For that, the modified n-gram
precision metric is introduced with an additional brevity penalty referring to
the sentence length. In addition, we performed a qualitative analysis of
selected generated texts. The samples were categorised as good or bad samples
according to their assigned METEOR score. The best-performing model parameters
from our experiments were used: $\lambda_{pos}=0.3$, $\gamma_{null}=4$ and
$\gamma_{other}=50$.

###  IV-B Standard Model Performance on SAX and BDD-X

####  IV-B1 Quantitative Results

TABLE I: This table shows the corpus-level METEOR and BLEU scores as
percentages.

| METEOR [%] | BLEU [%]
---|---|---
Description | Explanation | Description | Explanation
BDD-X (70%) | 45.92 | 10.47 | 14.82 | 0.33
SAX | 62.82 | 68.30 | 29.57 | 41.11

Table I shows METEOR and BLEU scores as percentage values for our model
implementation trained on 70% of the BDD-X dataset (due to unavailability of
the remaining 30%) and on the SAX dataset. The scores are individually
calculated for the description (e.g. "Car is stopping") and explanation part
(e.g. "because traffic light is not green on ego’s lane") of each generated
sentence.

On the BDD-X dataset, the model generated more accurate phrases for the action
description part than for the explanation part: The METEOR score for the
description is $4.4\times$ higher while the BLEU score is $44.9\times$ higher.
For the SAX dataset, such a significant difference was not observed. The
scores for the explanation were slightly higher than the description part.
Explanation scores were 41.11% against 29.57% (BLEU) and 68.30% against 62.82%
(METEOR). Overall, the language generator creates significantly higher scoring
samples for the SAX dataset than for the BDD-X dataset, especially, for the
explanation part.

####  IV-B2 Qualitative Analysis

From the qualitative analysis performed using the BDD-X data, results indicate
that no sample with a completely correct sentence was generated. The most
common errors were word repetitions and single missing words in the generated
sentences. For the good scoring samples, the semantics were captured correctly
by the language generator. The image frames had clearer visuals. The bad
scoring samples did not contain correct grammar (e.g. "the car is forward" or
"a to turn") and did not express the correct meaning. The ground truth
sentences are longer and more complicated compared to the good-scoring
samples. Additionally, the video frames had unclear visuals like sun glare or
bad illumination.

For the SAX dataset, clear differences in quality between action description
and explanation were observed: The description part is within the good scoring
samples mostly correctly generated. The explanation part is not once perfectly
generated. A few samples come close, for example by substituting the word
"ego" with "car". Although the bad sentences have mostly valid grammar,
semantics are wrongly captured.

####  IV-B3 Discussion

One reason why the BDD-X trained model scores higher on the description than
the explanation part is because describing a scenario is semantically easier
than finding a causal explanation. Another possible reason is based on
sentence grammar: The free text annotation process of the BDD-X set results in
largely varying samples, therefore, learning the sentence syntax and length
could be especially difficult for the model. Probably, this has a stronger
impact on the explanation part due to more variant expressions for post-hoc
explanations than for descriptions.

The model trained on the SAX dataset generated higher scores because the
ground truth sentences were annotated in a more structured way. This results
in less variety of used words and sentence lengths which allows the model to
generate high-scoring samples easily. One reason for this might be that the
generation of the right amount of NULL tokens towards the end of the sentence
was comparably easy using the uniform SAX samples and positively influences
the explanation scores. Another reason might be due to more realistic ground
truth explanations, because of the incorporated real-time driver commentary of
the SAX set. This could increase data quality and therefore ease the model’s
learning process.

On the qualitative analysis, one can conclude that further improvements in
language generation are needed. In the case of the BDD-X dataset, especially
word repetitions, grammar and sentence lengths need to be improved. The
generation of shorter sentences is probably due to the tendency of the network
to generate <NULL> tokens towards the end of the sentence instead of trying to
predict the correct word as a result of uncertainty. For the SAX dataset,
grammar is mostly correctly generated, but wrong semantics are often
expressed. A reason for that could be the small dataset size nature of SAX
compared to the BDD-X set.

###  IV-C Evaluation of Language Generator Modifications

####  IV-C1 Quantitative Results

Table II shows the METEOR and BLEU scores for the new model variants compared
to the standard model using the BDD-X dataset. For both description and
explanation generation, the model variations generated higher-scoring samples.
The best explanations were generated using the combination of PoS prediction
and special token penalties. The improvement is significant as it yields
METEOR scores that were $1.5\times$ better and BLEU scores that were
$7.7\times$ higher. Although, absolute scores were still significantly lower
compared to the description scores. The best descriptions were generated using
only token penalties. Here, the improvement is less significant: The METEOR
score improves only slightly by $1.6\%$ whereas the BLEU score improves more
by $5.1\%$ which equals $1.3\times$ better performance.

Table III shows the METEOR and BLEU scores for the model variants using the
SAX dataset. Here, a more significant improvement can be seen for description
generation: The combination of PoS prediction and token penalties performs
best with a $10.9\%$ higher METEOR score and $5,8\%$ higher BLEU score.
Regarding explanation generation, the best-performing model consists only of
PoS prediction with marginal improvement: $0.8\%$ higher according to METEOR
and $0.7\%$ higher according to BLEU. Both models that included special token
penalties generated worse scoring explanations compared to the standard model
on SAX: 5.9-8.2% lower METEOR scores and 3.3-8.1% lower BLEU scores.

TABLE II: This table shows the corpus-level METEOR and BLEU scores as
percentages for the variations of the language generator on the BDD-X dataset.

| METEOR [%] | BLEU [%]
---|---|---
Description | Explanation | Description | Explanation
Standard Model | 45.92 | 10.47 | 14.82 | 0.33
Standard + PoS Prediction | 46.06 | 12.49 | 17.54 | 0.42
Standard + Token Penalties | 47.56 | 12.47 | 19.95 | 1.12
Standard + PoS Predicton \+ Token Penalties | 44.64 | 15.95 | 17.44 | 2.53

TABLE III: This table shows the corpus-level METEOR and BLEU scores as
percentages for the variations of the language generator on the SAX dataset.

| METEOR [%] | BLEU [%]
---|---|---
Description | Explanation | Description | Explanation
Standard Model | 62.82 | 68.30 | 29.57 | 41.11
Standard + PoS Prediction | 66.27 | 69.10 | 29.81 | 41.84
Standard + Token Penalties | 64.69 | 60.12 | 34.75 | 32.95
Standard + PoS Predicton \+ Token Penalties | 73.76 | 62.40 | 35.37 | 37.84

Table IV shows variations of $\lambda_{pos}$ and $\gamma_{null}$ and their
impact on the METEOR and BLEU scores. Varying the part of speech weight
$\lambda_{pos}$ does not significantly change the performance. The differences
between $\lambda_{pos}=0.15$ and $\lambda_{pos}=0.3$ according to METEOR
(0.2-0.5%) and BLEU (0.1-0.7%) are minimal. Varying the NULL penalty
$\gamma_{null}$ led to more significant performance differences: Here, a lower
penalty led to higher scoring sample generation - especially for the
description part: 47.23% (METEOR) and 20.38% (BLEU) for $\gamma_{null}=4$
compared to 44.15% (METEOR) and 17.09% (BLEU) for $\gamma_{null}=12$.

TABLE IV: This table shows the corpus-level METEOR and BLEU scores as
percentages on the BDD-X dataset for different $\lambda_{pos}$ and
$\gamma_{null}$. For the shown results: $\gamma_{other}=0$.

| METEOR [%] | BLEU [%]
---|---|---
Description | Explanation | Description | Explanation
$\lambda_{pos}=0$ \+ $\gamma_{null}=0$ | 45.92 | 10.47 | 14.82 | 0.33
$\lambda_{pos}=0.15$ \+ $\gamma_{null}=0$ | 45.60 | 12.69 | 16.76 | 0.52
$\lambda_{pos}=0.3$ \+ $\gamma_{null}=0$ | 46.06 | 12.49 | 17.54 | 0.42
$\lambda_{pos}=0$ \+ $\gamma_{null}=4$ | 47.23 | 11.57 | 20.38 | 0.54
$\lambda_{pos}=0$ \+ $\gamma_{null}=12$ | 44.15 | 10.64 | 17.09 | 0.43

####  IV-C2 Qualitative Analysis

Qualitatively, it can be seen that the good-scoring samples capture the
correct semantics. Additionally, the sentences were almost identical to the
ground truth. For the bad scoring samples, the explanation part represents no
complete and syntactically correct sentence. The description parts had valid
syntax, but the wrong meaning. On top of that, the ground truth explanations
had complex sentence structures.

For the SAX dataset, Fig. 3 shows a good and a bad scoring sample. For the
good, the action description and explanation are correctly generated. The bad
sentence has description parts with correct grammar, whereas the explanation
sentences had few cases of word repetitions, with the semantics mostly
captured wrongly.

Figure 3: This figure shows a good (top) and a bad (bottom) scoring generated
sample using the model with PoS prediction trained on the SAX dataset. A
corresponding video frame with the spatial attention map as an overlay is
shown next to each sample.

In Fig. 4 one sample of the BDD-X dataset is shown (top) as ground truth
together with three generated sentences according to the standard model, the
special token penalty model and the combination model of PoS prediction and
token penalties. These two variations were chosen, because they were the
quantitatively best-performing ones. In the standard model’s generation, the
last word of each sentence part is missing. Additionally, instead of the term
"driving down", the term "driving forward" with a similar meaning was
generated. In contrast, the combination model generates a nearly perfect
sentence. The model using only token penalties has the worst qualitative
performance.

In Fig. 4 equivalently one sample of the SAX dataset is shown (bottom) with
its ground truth annotation. Here, the generated sentences are produced by the
standard model, the model with PoS prediction and the combination model.
Again, these two variations were chosen, because they were the quantitatively
best-performing ones. In the standard models’ sentence, the explanation was
correctly generated, but the description portrays the wrong meaning, even
though the grammar was correct. The PoS prediction model achieves to improve
the description by capturing the correct meaning and therefore generating the
correct complete sentence. The combination model consisting of PoS prediction
and token penalties generated a faulty sentence: Wrong semantics for the
description and explanation and wrong sentence structure for the explanation
part.

Figure 4: This figure shows one sample of the BDD-X (top) and of the SAX
dataset (bottom) with its ground truth annotation and three generated texts by
the standard model, the model with special token penalties (BDD-X) / with PoS
prediction (SAX) and the combination model with PoS prediction and token
penalties.

####  IV-C3 Discussion

On the BDD-X dataset, the model adjustments led to a relevant performance
increase. Especially, the token penalties had a significant impact according
to the BLEU scores. The highest performance increase can be observed for the
explanation part with the combination of PoS prediction and token penalties. A
possible reason could be that the quality of the explanation part usually is
dependent on the correct prediction of the separator token and the correct
amount of NULL tokens in the end. Therefore, the explanation generation is
strongly influenced by the token penalties.

For the SAX dataset, only the description part could be significantly
improved. Again, the combination of PoS prediction and token penalties yielded
the best performance. For the explanation part, token penalties had a negative
influence. This could be due to the uniform annotation corpus of the dataset
that makes it already manageable to learn the correct syntax and sentence
lengths with the original model. Adding additional NULL penalties might have
led to predicting wrong random words instead of the correct NULL tokens.

Testing different strengths of $\lambda_{pos}$ and $\gamma_{null}$ led to the
conclusion that the weight of the PoS loss $\lambda_{pos}$ has no significant
influence on the model performance, but the strength of the null penalty
$\gamma_{null}$ has. An optimal $\gamma_{null}$ value exists.

Qualitatively, the proposed modifications also increased the performance. The
sentence structure was mostly correctly generated and the remaining challenges
refer to the model capturing the wrong semantics of a scene. This is
especially true for examples with complex ground truth captions. The
performance of the individual models varies depending on the used dataset -
probably due to different complex ground truth annotations. On an individual
sample level, the qualitative assessment of the generated sample is not always
aligned with the overall quantitative score on the complete dataset as can be
seen in Fig. 4.

##  V Conclusion

This work validates the performance of the SOTA prediction and explanation
pipeline for high-level driving actions on the new SAX dataset with improved
performance. The unique features of the SAX dataset, e.g., uniformly
structured annotation corpus and real-time driver commentary allow the model
to score higher and generate sentences with better sentence grammar compared
to the BDD-X set. To further address the faulty sentence grammar in many BDD-
X-based generated samples, we modify the language generator by adding part of
speech prediction and special token penalties in the loss function. This
resulted in a significant improvement of the explanation generation by a
factor of $7.7$ and $1.5$ with the BLEU and METEOR metrics, respectively. The
description generation was also improved by a factor of $1.3$. Our modified
pipeline outperforms the SOTA end-to-end model.

Nevertheless, our work has limitations: First, the text evaluation metrics
BLEU and METEOR mainly rely on sentence similarity and, thus, are limited in
comparing semantics. In future work, we would validate our results using human
evaluation. Second, optimal weight and penalty values would benefit from a
more rigorous estimation. Third, the model’s ability to learn correct
semantics could be further improved by using a higher performance model for
feature processing, e.g., a completely attention-based transformer [19] or
using more sophisticated image representations outputted by an advanced
encoder like CLIP [20].

## References

  * [1] J. Kim, A. Rohrbach, T. Darrell, J. Canny, and Z. Akata, “Textual explanations for self-driving vehicles,” in Computer Vision – ECCV 2018, pp. 577–593, Springer International Publishing, 2018.
  * [2] M. Gadd, D. D. Martini, L. Marchegiani, P. Newman, and L. Kunze, “Sense-assess-explain (SAX): building trust in autonomous vehicles in challenging real-world driving scenarios,” CoRR, vol. abs/2005.02031, 2020\.
  * [3] A. B. Arrieta, N. Díaz-Rodríguez, J. Del Ser, A. Bennetot, S. Tabik, A. Barbado, S. García, S. Gil-López, D. Molina, R. Benjamins, et al., “Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai,” Information Fusion, vol. 58, pp. 82–115, 2020.
  * [4] D. Omeiza, H. Webb, M. Jirotka, and L. Kunze, “Explanations in autonomous driving: A survey,” IEEE Transactions on Intelligent Transportation Systems, pp. 1–21, 2021.
  * [5] X.-H. Li, C. C. Cao, Y. Shi, W. Bai, H. Gao, L. Qiu, C. Wang, Y. Gao, S. Zhang, X. Xue, and L. Chen, “A survey of data-driven and knowledge-aware explainable ai,” IEEE Transactions on Knowledge and Data Engineering, vol. 34, no. 1, pp. 29–49, 2022.
  * [6] J. Kim and J. Canny, “Interpretable learning for self-driving cars by visualizing causal attention,” in Proceedings of the IEEE International Conference on Computer Vision (ICCV), Oct 2017.
  * [7] Y. Mori, T. Hirakawa, T. Yamashita, and H. Fujiyoshi, “Image captioning for near-future events from vehicle camera images and motion information,” in 2021 IEEE Intelligent Vehicles Symposium (IV), pp. 1378–1384, 2021.
  * [8] Y. Feng, L. Ma, W. Liu, and J. Luo, “Unsupervised image captioning,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [9] X. Yang, K. Tang, H. Zhang, and J. Cai, “Auto-encoding scene graphs for image captioning,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [10] W. Wang, Z. Chen, and H. Hu, “Hierarchical attention network for image captioning,” Proceedings of the AAAI Conference on Artificial Intelligence, vol. 33, pp. 8957–8964, Jul. 2019.
  * [11] J. Kim, T. Misu, Y.-T. Chen, A. Tawari, and J. Canny, “Grounding human-to-vehicle advice for self-driving vehicles,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2019.
  * [12] D. Omeiza, S. Anjomshoae, H. Webb, M. Jirotka, and L. Kunze, “From spoken thoughts to automated driving commentary: Predicting and explaining intelligent vehicles’ actions,” in 33rd IEEE Intelligent Vehicles Symposium (IV), IEEE Press, 2022.
  * [13] S. Bird, E. Klein, and E. Loper, Natural Language Processing with Python.  O’Reilly Media Inc, 2009.
  * [14] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014.
  * [15] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015.  Software available from tensorflow.org.
  * [16] A. Lavie and A. Agarwal, “Meteor: An automatic metric for mt evaluation with high levels of correlation with human judgments,” in Proceedings of the second workshop on statistical machine translation, pp. 228–231, Association for Computational Linguistics, 2007.
  * [17] M. Post, “A call for clarity in reporting BLEU scores,” in Proceedings of the Third Conference on Machine Translation: Research Papers, (Belgium, Brussels), pp. 186–191, Association for Computational Linguistics, Oct. 2018\.
  * [18] K. Papineni, S. Roukos, T. Ward, and W.-J. Zhu, “Bleu: a method for automatic evaluation of machine translation,” in Proceedings of the 40th annual meeting of the Association for Computational Linguistics, pp. 311–318, 2002\.
  * [19] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, “Attention is all you need,” in Advances in Neural Information Processing Systems (I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, eds.), vol. 30, Curran Associates, Inc., 2017.
  * [20] A. Radford, J. W. Kim, C. Hallacy, A. Ramesh, G. Goh, S. Agarwal, G. Sastry, A. Askell, P. Mishkin, J. Clark, G. Krueger, and I. Sutskever, “Learning transferable visual models from natural language supervision,” 2021.