# Unveiling a hidden bar-like structure in NGC 1087: kinematic and photometric
evidence using MUSE/VLT, ALMA and JWST

Carlos López-Cobá Institute of Astronomy and Astrophysics, Academia Sinica,
No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan  Lihwai Lin Institute
of Astronomy and Astrophysics, Academia Sinica, No. 1, Section 4, Roosevelt
Road, Taipei 10617, Taiwan  Sebastián F. Sánchez Instituto de Astronomía,
Universidad Nacional Autonoma de México, Circuito Exterior, Ciudad
Universitaria, Ciudad de México 04510, Mexico

###### Abstract

We report a faint non-axisymmetric structure in NGC 1087 through the use of
JWST Near Infrared Camera (NIRCam), with an associated kinematic counterpart
observed as an oval distortion in the stellar velocity map, H$\alpha$ and CO
$J=2\rightarrow 1$ velocity fields. This structure is not evident in the MUSE
optical continuum images but only revealed in the near-IR with the F200W and
F300M band filters at $2\mu$m and $3\mu$m respectively. Due to its elongation,
this structure resembles a stellar bar although with remarkable differences
with respect to conventional stellar bars. Most of the near-IR emission is
concentrated within $6\arcsec~{}\sim 500$ pc with a maximum extension up to
1.2 kpc. The spatial extension of the large-scale non-circular motions is
coincident with the bar, which undoubtedly confirms the presence of a non-
axisymmetric perturbation in the potential of NGC 1087. The oval distortion is
enhanced in CO due to its dynamically cold nature rather than in H$\alpha$. We
found that the kinematics in all phases including stellar, ionized and
molecular, can be described simultaneously by a model containing a bisymmetric
perturbation; however, we find that an inflow model of gas along the bar major
axis is also likely. Furthermore the molecular mass inflow rate associated can
explain the observed star formation rate in the bar. This reinforces the idea
that bars are mechanisms for transporting gas and triggering star formation.
This work contributes to our understanding of non-axisymmetry in galaxies
using the most sophisticated data so far.

##  1 Introduction

Stellar bars or “bars” are one of the most visual signs of non-axisymmetry in
galaxies (de Vaucouleurs et al., 1991). Like other morphological structures in
galaxies, bars are not well defined, although they are relatively easy to
identity in composite images because of their bar-shape structure (e.g.,
Masters et al., 2011; Cheung et al., 2013), and their prominent bar dust lanes
(Athanassoula, 1992). NGC 1087 is a clear example where the combination of
high spatial resolution data allowed us to reveal a hidden bar that is only
observed in the near-infrared (NIR). The James Webb Space Telescope (JWST) is
allowing us to reveal structures not seen before due to the lack of resolution
and sensitivity. Their infrared bands are prone to detect the emission from
old-stellar structures like bars $\sim 10$ Gyr (Sánchez-Blázquez et al.,
2011). In the optical, dust obscuration can prevent their detection, affecting
the global statistics of galaxies hosting bars in the nearby Universe (e.g.,
Sellwood & Wilkinson, 1993). Apart from their optical and NIR characteristics,
bars leave imprints of non-axisymmetry in the kinematics of gas and stars
(e.g., Wong et al., 2004; Fathi et al., 2005; López-Cobá et al., 2022); thus
studying the dynamical effects of bars at local scales is crucial for
understanding the role they play in galaxy evolution. A comprehensive study of
bars could be addressed if kinematics and wide-range photometric observations
were accessible for a considerably large sample of galaxies. Yet, spatial
resolution plays an important role in separating the different structural
components of galaxies. For example, the large integral field spectroscopic
surveys like CALIFA (Sánchez et al., 2012) or MaNGA (Bundy et al., 2015),
although providing a large statistical sample of galaxies, their nominal
resolution of $2\farcs 5\sim 1$ kpc inhibits a detailed study of individual
galaxies, while biasing the sample towards bar-lengths larger than the FWHM
resolution. In a similar manner NIR photometric and optical catalogs provide
resolution of few arcseconds which again limits the detection of small bars
(Sheth et al., 2010).

In this work, we make use of the most sophisticated data to unravel the
detection of a bar in NGC 1087. This paper is structured as follows: in
Section 2 we describe the data and data analysis; in Section 3 we address the
detection of a faint stellar bar in the infrared while their ionized and
molecular counterpart are addressed in Section 4; Section 5 describes the oval
distortion and in Section 6 we present the discussion and conclusions.

##  2 Data and Data analysis

NGC 1087 is an intermediate Sc galaxy located at 14 Mpc (e.g., Kourkchi &
Tully, 2017); at this distance the physical spatial resolution is 68
pc$\mathrm{~{}arcsec^{-1}}$.

This work is based on public data from the Multi-Unit Spectroscopic Explorer
(MUSE, Bacon et al., 2010); data from the JWST Near Infrared Camera (NIRCam);
and ALMA CO $J=2\rightarrow$1 observations, in particular we used the data
products from the PHANGS–ALMA survey (e.g., Leroy et al., 2021a, b), namely
moments 0 and 1. The MUSE-VLT data was obtained as part of the recent release
of the MUSE-PHANGS datacubes (e.g., Schinnerer, 2021; Emsellem et al., 2022).
MUSE is an integral field spectrograph which provides spatially resolved
spectra over a $1\arcmin\times 1\arcmin$ field of view (FoV), covering the
optical spectrum, with a spectral resolution at full width at half maximum
(FWHM) of $\sim 2.6$Å. The estimated spatial resolution of the $2\arcmin\times
3\arcmin$ MUSE mosaic covering NGC 1087 is not worse than $0\farcs 9$/FWHM
(e.g., Emsellem et al., 2022). We used fully calibrated data products from
JWST from the JWST Science Calibration Pipeline version 1.10.1 (e.g., Bushouse
et al., 2022). Specifically, we made use of the PHANGS-JWST first result
products (e.g., Rosolowsky, Erik, 2022), namely, the NIRCam F200W filter (FWHM
$\sim 0\farcs 05$ and $0\farcs 031/\mathrm{pixel}$ sampling), which has its
nominal wavelength at $1.99\mu$m, in addition to the F360M and F335M band
filters to trace the stellar continuum. Finally the spatial resolution of the
ALMA data is $1\farcs 60$/FWHM following Leroy et al. (2021a). The CO moment 0
map was transformed to surface density ($\mathrm{\Sigma_{mol}}$), assuming a
standard Milky Way
$\mathrm{\alpha_{\mathrm{CO}}^{1-0}=4.35\,M_{\odot}/pc^{2}}$ conversion factor
and $\mathrm{\Sigma_{mol}[M_{\odot}pc^{2}]=6.7\,I_{CO(2-1)}[K\,km/s]}\cos i$
(e.g., Eq. 11 from Leroy et al., 2021a), with $i$ being the disk inclination
angle estimated in $44.5^{\circ}$.

The MUSE data analysis was made with the pypipe3d tool (e.g., Lacerda et al.,
2022). In short, pypipe3d performs a decomposition of the observed stellar
spectra into multiple simple stellar populations (SSPs) each with different
age and metalicities. For concordance with Emsellem et al. (2022), we use the
same stellar libraries based on the extended MILES (E-MILES, Vazdekis et al.,
2016), FWHM $=2.5$Å; We convolved the spectral resolution of the SSPs to that
of the MUSE line-spread function (e.g., Bacon et al., 2017). To increase the
signal-to-noise (SN) of the stellar continuum, we performed a Voronoi binning-
segmentation on a 2D-slice centered around the V-band, ensuring bin-sizes of
the order of the MUSE-FWHM resolution and SN $\sim 50$. This binned map will
serve to compute different properties of the underlying stellar continuum,
while the ionized gas properties are analyzed in a spaxel-wise sense. The
result of pypipe3d is a set of maps comprising information about the stellar
populations (stellar velocity, stellar mass, among others products), and the
ionized gas (emission-line fluxes including: H$\alpha$, H$\beta$, [S
ii]$\lambda 6717,6731$, [N ii]$\lambda 6584$, [O iii]$\lambda 5007$, [O
i]$\lambda 6300$, emission-line velocities, equivalent widths etc); see
Lacerda et al. (2022) for a thorough description of the analysis and
dataproducts. From this analysis we estimate the stellar mass of this object
in $\mathrm{\log M_{\star}/M_{\odot}=9.9}$.

Figure 1: Central panel: False color image showing the distribution of the
ionized gas and NIR continuum in NGC 1087 with fluxes taken from MUSE-VLT and
NIRCam imaging respectively (red: [S ii], yellow: H$\alpha$, blue: [O iii],
brown: F360M, orange: F335M, white: F200W). Right panels: A zoom-in of the
$30\arcsec\times 20\arcsec$ innermost region is shown to highlight the bar-
structure. The upper panel shows a true color image from the MUSE cube (R:
i-band, G: r, B: g). The black cross in the middle represents the NIR nucleus.
The middle panel shows the JWST-F200W image together with a set of isophotes
describing the bar NIR light distribution; the isophotes show a constant
orientation in the sky at $\phi^{\prime}_{bar}=304^{\circ}\pm 3^{\circ}$. Two
isophotes located at $6\arcsec$ and $14\arcsec$ are highlighted with thicker
lines for reference. In each case the FWHM resolution is shown with a yellow
circle. Bottom inset: Orientation of the object in the sky.

##  3 FAINT stellar bar

The optical continuum image of NGC 1087 exhibits a bright, fuzzy and
featureless nucleus with several dust lanes as revealed by the $gri$ MUSE
image in Figure 1, because of that it is unclear whether this object shows a
stellar bar in the optical. Two bright spots are observed in the central
region where the near-IR (NIR) nucleus is found. Unlike the majority of
galaxies hosting stellar bars, NGC 1087 does not show a well defined bar,
nucleus neither a star forming ring (e.g., de Vaucouleurs et al., 1991).
Conversely, the high spatial resolution from the F200W imaging filter allows
to resolve a clear elongated and clumpy structure of presumably stellar
clusters aligned along a preferential position angle (P.A.), which differs
from the disk orientation111From now on primed variables make reference to
values measured on the sky plane, otherwise in the disk plane.
$\phi_{disk}^{\prime}=359^{\circ}$, as observed in Figure 1. At $2\mu$m this
structure resembles a faint stellar bar and it is obscured in the optical due
to dust absorption. F200W traces mostly stellar continuum, hence, it is
expected that most of the continuum emission along this structure comes from
old stars (e.g., Leitherer et al., 1999). Elliptical isophotes with variable
position angle and ellipticity ($\varepsilon^{\prime}$) show a preferential
alignment of the $2\mu$m emission as observed in the top-right panel from
Figure 1.

Figure 2: NGC 1087 as observed with the F200W and F300M NIRCam band filters.
First panel: zoom around a $30\arcsec\times 20\arcsec$ size window containing
the bar-like structure. second panel: Best non-parametric model for the
$2\mu$m light distribution using the photometric version of the $\mathtt{XS}$
code (in prep.). The coordinates of the photometric center was estimated in
$\alpha_{2000}=2^{h}46^{m}25.170^{s}$,
$\delta_{2000}=-00^{\circ}29^{\prime}55.728\arcsec$, while the position angle
of the disk describing the bar light distribution was estimated in
$\phi_{bar}^{\prime}=305^{\circ}$ with a $0.5$ ellipticity. Third panel:
residuals of the modeling (observed-model). Fourth panel: the black line
represents the light distribution profile of the model shown in the second
panel, while the blue dashed line represents the cumulative light distribution
of the model. Distances are measured following these projection angles. Figure
3: Diagnostic diagrams for the ionized gas in the bar structure. Top panel: [O
i]6300/H$\alpha$ vs. [O iii]/H$\beta$ diagram. The black dashed line
represents the star forming demarcation curve from Kewley et al. (2006). Shock
grid models from the updated MAPPINGS-V (Sutherland et al., 2018) computed by
Alarie & Morisset (2019) are shown for $\mathrm{n_{e}=10\,cm^{-3}}$,
$\mathrm{Z_{ISM}=0.008~{}(equivalent~{}to~{}12+\log O/H=8.55)}$, and a wide
range of magnetic fields (blue lines) and pre-shock velocities (red lines).
The colors of the points map the F200W light distribution shown in the inset
figure. Bottom panel: EW(H$\alpha$) in absolute value vs. [N ii]/H$\alpha$
line ratio. Old and evolved stars are expected to produce EW(H$\alpha$)
$<3$Å..  Figure 4: Oxygen abundance in the bar-region computed with the
Pilyugin & Grebel (2016) calibrator. As in Figure 1 the white ellipses
describe the NIR bar. The innermost ellipse shows an average metallicity of
$\mathrm{12+\log O/H=8.5}$.  Figure 5: Top: Specific star formation rate
computed through the H$\alpha$ based SFR and the stellar mass density maps,
$\mathrm{sSFR=SFR/M_{\star}}$. Bottom: Molecular surface density around the
bar structure. The H$\alpha$ flux is overlaid with reddish colors for
comparison. The yellow circles represent the FWHM spatial resolution in each
case. The white ellipses delineate the NIR bar extension.  Figure 6: First
column: Line-of-sight velocities of CO, H$\alpha$ and the stellar velocity in
NGC 1087. Second column: Zoom-in around the bar-structure, with $\pm 25$ km
s-1 spaced iso-velocity shown on-top; the black straight line shows the P.A.
of the faint bar at $\phi_{bar}^{\prime}=305^{\circ}$. Third column: residual
velocities after subtracting the circular rotation in each velocity map.
Fourth column: Bisymmetric model after fixing the oval orientation at
$\phi_{bar}^{\prime}$. Fifth column: inflow model assuming the gas streams
along $\phi_{bar}^{\prime}$; the individual radial profiles and expressions
are shown in Figure 7. The root mean square of the models is shown at the
bottom right of each panel. All maps share the same color-bar except for the
residual maps which cover $\pm 50$ km s-1 following the same colors scheme.

###  3.1 Bar light modeling

The method adopted here is similar to the non-parametric method from Reese et
al. (2007), which is the same algorithm used by $\mathtt{XS}$ (e.g., Lopez-
Coba et al., 2021) for extracting the different velocities in the kinematic
models, as we explain in the following sections. Our non-parametric model
minimizes the function $\chi^{2}=(I_{obs}-I_{k}(r)w_{k}(r))^{2}$, where
$I_{obs}$ is the observed intensity map, $I_{k}$ represents a set of
intensities that will be estimated at different annuli, and $w_{k}$ is a set
of weighting factors that will serve for performing a linear interpolation
between the estimated $I_{k}$ intensities. The implementation on the F200W and
F300M band filters is shown in Figure 2.

The NIR light from the bar-structure can be successfully modeled by an
elliptical light distribution with with a constant orientation in the sky. The
P.A. of such ellipse is $\phi_{bar}^{\prime}=305^{\circ}$222Sky angles
relative the disk major axis can be translated to angles measured in the
galaxy plane and vise versa through, $\tan\phi_{bar}=\tan\Delta\phi/\cos i$,
with $\Delta\phi=\phi_{bar}^{\prime}-\phi_{disk}^{\prime}$. and
$\varepsilon^{\prime}_{bar}=0.5$, as shown in Figure 2. From this analysis the
estimated semi-major axis length of the NIR bar is
$a_{bar}^{\prime}=14\arcsec$ or $18\arcsec\sim 1.2$ kpc333Distances measured
at different P.A. in the sky plane can be translated to distances in the
galaxy plane following:
$r_{disk}^{2}=r_{sky}^{2}(\cos^{2}\Delta\phi+\sin^{2}\Delta\phi/\cos^{2}i)$.
on the disk at the considered distance. However, this structure is very
diffuse, with most of the 2$\mu$m emission $\sim 70\%$, concentrated within
the inner $6\arcsec\sim 500$ pc. The latter region encloses the two bright
spots observed in the MUSE optical continuum images.

##  4 Ionized and molecular gas

The central panel of Figure 1 shows the ionized gas distribution traced by the
H$\alpha$, [S ii] and [O iii] emission-line fluxes. In general, galaxies
hosting stellar bars do not frequently present cold or warm gas along the bar.
This is often explained as a lost of angular momentum after the gas encounters
the offset ridges, resulting in an infall of gas towards a central ring (e.g.,
Athanassoula, 1992). NGC 1087 shows plenty of ionized gas throughout the disk
and within the bar region. Furthermore, as we will see later, multiple
enhanced CO bloops are spatially coincident with the F200W emission. To
investigate the dominant ionizing source along the bar-like structure we made
use of line-ratios sensitive to the ionization, in particular the [O
i]6300/H$\alpha$ and [O iii]/H$\beta$. The [O i]6300/H$\alpha$ ratio has the
advantage of being a good indicator of the presence of shocks. The top panel
from Figure 3 shows these line-ratios color-coded with the $2\mu$m emission.
The observed ionized gas is compatible with being produced by star-formation
(SF) according to the Kewley et al. (2006) demarcation line. However, dust
lanes in bars are often associated with shocks as a result of the complex gas
dynamics (Athanassoula, 1992). The bar-region in NGC 1087 shows several dust
lanes, thus we adopt shock models to investigate whether shock ionization can
reproduce the observed line ratios. For this purpose we used the
photoionization grids from MAPPINGS V (e.g., Sutherland et al., 2018) computed
by Alarie & Morisset (2019); we adopt the Gutkin et al. (2016) metallicities.
The latter parameter was set from the oxygen gas-phase abundance, adopting the
Pilyugin & Grebel (2016) calibrator based on strong emission-lines (see Figure
4). For the bar-region we find $\mathrm{12+\log O/H\sim 8.5}$ corresponding to
$\mathrm{Z_{ISM}=0.008}$ in the Gutkin et al. (2016) models. For this
metallicity, the predicted emission-line intensities from shock models can not
explain the observed [O i]/H$\alpha$ line ratios as observed in the top panel
of Figure 3. Thus, if shocks are happening in the bar-region, then their
optical emission is not expected to be the dominant ionizing source. Although
there is still a possibility that a combination of SF plus shocks could
contribute to the observed line-ratios (e.g., Davies et al., 2016).

Additionally, Figure 3 shows the equivalent width of H$\alpha$
(EW($\rm{H}\alpha$)) for the ionized gas in the bar-region. Hot, old, and low-
mass evolved stars are characterized for producing
EW($\rm{H}\alpha$)$~{}\lesssim 3$Å (e.g., Stasińska et al., 2008; Cid
Fernandes et al., 2010; Lacerda et al., 2018). The large values of
EW($\rm{H}\alpha$) $\gtrsim 100$ Å found in the bar, can not be explained by
the ionizing continuum emitted by this population of stars. Therefore, if most
of the $2\mu$m emission observed along the bar is due to old stars, then their
ionizing continuum is not sufficient to explain the observed line-ratios.

The star formation rate (SFR) was estimated from the H$\alpha$ luminosity
(Kennicutt, 1998), after correcting the H$\alpha$ flux from dust extinction
adopting the Cardelli et al. (1989) extinction law with $\mathrm{R_{V}=3.1}$,
and assuming an intrinsic flux ratio of H$\alpha$/H$\beta$= 2.86 and ionized
gas temperature of T $\sim 10^{4}$ K corresponding to case B recombination
(e.g., Osterbrock, 1989). The integrated SFR within the NIR bar (i.e., Figure
2) is $\mathrm{SFR(H\alpha)_{bar}}=0.08~{}\mathrm{M_{\odot}/yr}$, while the
total $\mathrm{SFR(H\alpha)_{total}}=0.4~{}\mathrm{M_{\odot}/yr}$. The
specific star formation rate (sSFR) was obtained from dividing the SFR surface
density (=SFR/pixel area) by the stellar mass surface density as shown in
Figure 5. We note a clear enhancement in SF along the bar main axis;
furthermore, the molecular surface density shows an enhancement in the same
region.

Overall, our analysis suggests that the ionized gas in the bar structure is
mostly associated to SF processes. Following we investigate whether the CO and
SF concentration is induced by the bar potential.

##  5 Inner Oval distortion

If a perturbation in the gravitational potential, such as that induced by a
bar potential, is causing the $2\mu$m light distribution to elongate along a
preferred direction; then the particles in such structure are expected to
follow quasi-elliptical orbits (e.g., Athanassoula, 1992). Hence kinematics of
gas and stars should reflect this perturbation (Pence & Blackman, 1984).

Figure 6 shows the CO, ionized and stellar kinematics around the bar region.
The gas, being collisional, is more sensitive to non-axisymmetric
perturbations. This is clearly reflected in the CO moment 1 map where a strong
distortion of the semi-minor axis is observed. The former misalignment is a
signature of an oval distortion produced by a bar (e.g., Pence et al., 1988;
López-Cobá et al., 2022), here observed at high spatial resolution. The
distortion is also observed in the H$\alpha$ velocity map, however, since the
ionized gas traced by H$\alpha$ is hotter ($\mathrm{T\sim 10^{4}}$ K), and
therefore it presents a larger intrinsic velocity dispersion than the one of
the molecular gas traced with CO ($\mathrm{T\sim 10}$ K), it is less
pronounced in the ionized gas and is probably affected by local SF. The
stellar kinematics although affected by the pixel-coadding during the SSP
analysis, it still reveal a clear distortion in the iso-velocities near the
bar as noted in Figure 6.

So far, it is clear that the stars, the molecular gas traced with CO and the
ionized gas traced with H$\alpha$ respond in a similar way to the presence of
this bar-like structure detected in NGC 1087.

###  5.1 Kinematic interpretation of the oval distortion

The distortion in the velocity field such as the ones observed before could be
caused by an elongated potential. However, only a few kinematic models in the
literature attempt to describe the flow caused by an oval distortion (e.g.,
Spekkens & Sellwood, 2007; Maciejewski et al., 2012). A bisymmetric distortion
induced by a second order perturbation to the potential, such as that induced
by an elongated potential, has been shown to successfully reproduce the
velocity field in bars (e.g., Spekkens & Sellwood, 2007). Since the F200W
image evidences the presence of a faint bar-like structure oriented at
$305^{\circ}$ in the sky, we performed bisymmetric models over the stellar,
H$\alpha$ and CO velocity maps, with a fixed orientation of the oval
distortion.

The bisymmetric model from Spekkens & Sellwood (2007) is described by the
following expression:

$V_{\mathrm{LoS}}=\sin i\Big{(}V_{t}(r)\cos\theta-V_{2,t}(r)\cos
2\theta_{\mathrm{bar}}\cos\theta\\\ -V_{2,r}(r)\sin
2\theta_{\mathrm{bar}}\sin\theta\Big{)}+V_{\mathrm{sys}}$ (1)

where $\theta_{bar}=\theta-\phi_{bar}$; with $\theta$ being the azimuthal
angle measured on the disk plane from the line of nodes and $\phi_{bar}$ the
position angle of the bar-like structure on the disk plane. $V_{t}$ is the
tangential or circular rotation and $V_{2r}$ and $V_{2t}$ represent the radial
and tangential velocities that result from a bisymmetric distortion to the
gravitational potential.

Before this analysis we performed a circular rotation model to obtain the disk
projection angles and the rotational curves. In this case and in the
subsequent, we used the $\mathtt{XS}$ code for generating the kinematic models
(e.g., Lopez-Coba et al., 2021). This code derives an interpolated model over
a set of concentric rings evenly spaced $r_{k}$, minimizing the function
$\chi^{2}=(V_{obs}-V_{k}(r)w_{k}(r))^{2}/\sigma^{2}$, where $V_{obs}$ is the
observed velocity map; $V_{k}$ are the set of velocities inferred at $r_{k}$;
$w_{k}$ is a set of weights that depend on the specific kinematic model
adopted and will serve to create a 2D model; $\sigma$ is the error velocity
map. We refer to Lopez-Coba et al. (2021) for a thorough description of this
analysis. It is important to mention that the modeling does not assume any
parametric function on the velocity profiles.

A pure circular rotation, without non-circular motions, is described by the
first term on the right from Equation 1. The disk orientation of NGC 1087 was
estimated from modeling the H$\alpha$ and stellar velocity maps with circular
rotation only, obtaining $\phi_{disk}^{\prime}=358.9^{\circ}$,
$i=44.5^{\circ}$ and $V_{sys}=1526$ km s-1. The residual velocity maps from
this model i.e., $V_{obs}-V_{circular}$, are shown in the third column from
Figure 6. The circular rotation models leave large-scale residual velocities
of the order of $30$ km s-1 with a mirror symmetry about the nucleus; this
corresponds to de-projected amplitudes of $\sim 50$ km s-1 for the non-
circular velocities on the disk plane. These residual patterns have been
observed in larger scales in galaxies with long bars in H$\alpha$ (e.g., Lang
et al., 2020; López-Cobá et al., 2022) and molecular gas (e.g., Pence et al.,
1988; Mazzalay et al., 2014).

The bisymmetric model with a fixed oval orientation at $\phi_{bar}^{\prime}$
is shown in the fourth column in Figure 6. This model reproduces
simultaneously the twisted iso-velocities in the three velocity maps. The root
mean square (rms) of the models decrease compared with the circular rotation
one, therefore, in terms of the residuals, the LoS-velocities observed around
the bar region can be reproduced by a kinematic model that considers a bar-
like perturbation as the main source of non-circular motions. However, at such
high spatial resolution $\sim 70$ pc, there are still residual velocities that
the bisymmetric model can not account for, this is reflected in the 15 km s-1
rms in the models, which is similar or larger than the level of turbulence of
the ISM (e.g., Moiseev et al., 2015).

An alternative interpretation to the observed non-circular motions in the
H$\alpha$ and CO velocity maps is the presence of gas inflow induced by the
NIR-bar (Mundell & Shone, 1999); in fact, hydro-dynamical models simulating
bars predict inflow of gas along the offset ridges, or bar dust lanes, (e.g.,
Athanassoula, 1992). As observed in Figure 1, the MUSE image of NGC 1087 shows
several dust lanes in the central region making difficult the identification
of those associated with the bar. Assuming the spiral arms are trailing (see
bottom-right panel from Figure 1), NGC 1087 rotates counterclockwise in the
sky. Thus, positive (negative) residuals in the near (far) side represent
inflow; hence, a radial flow ($V_{rad}$) along the bar major axis is inflowing
if $V_{rad}<0$. We implemented a non-axisymmetric model, with a flow streaming
along $\phi_{bar}^{\prime}$. The non-axisymmetric inflow model is described by
the following expression:

$V_{LoS}=(V_{t}(r)\cos\theta+V_{rad}(r)\cos\theta_{bar}\sin\phi_{bar})\sin
i\\\ +V_{sys}$ (2)

where $V_{rad}$ is the radial velocity of the flow, and $\phi_{bar}$ is the
misalignment between the disk and the bar position angles on the disk plane.
This model assumes the gas is flowing along $\phi_{bar}$. This expression is
similar to Hirota et al. (2009) and Wu et al. (2021) assuming that the gas
flows parallel to the bar major axis.

In order to consider pixels likely affected by this motion, we consider the
elliptical region that better describes the bar-like light distribution and
defined in Figure 2. The 2D representation of this model is shown in Figure 6,
while the kinematic radial profiles of all models are shown in Figure 7. As
observed, the radial velocity $V_{rad}$ results negative in the model, with
maximum inflow velocities of the order of 50 km s-1.

Figure 7: Radial distribution of the different velocity components for the
considered kinematic models. $V_{t}$ is the disk pure circular rotation, or
rotational curve; $V_{2r}$ and $V_{2t}$ are the bisymmetric velocities;
$V_{rad}$ the radial flow from the non-axisymmetric flow model. In all non-
circular models the position angle of the oval distortion was fixed to
$\phi_{bar}^{\prime}=305^{\circ}$ (namely, $297^{\circ}$ in the galaxy plane).
Lines in red colors represent results of models computed on the H$\alpha$
velocity map, blue lines on the CO moment 1 and black lines on the stellar
velocity map. Shaded regions represent $1\sigma$ errors. The non-circular
velocities were estimated up to $18\arcsec$, which covers the de-projected
length of the bar. A $2\arcsec$ sampling step was adopted, i.e., larger than
the FWHM spatial resolution of the data.

The molecular mass flow rate ($\dot{M}_{mol}$) associated can be computed
following the expression (e..g, Di Teodoro & Peek, 2021):

$\dot{M}_{mol}(r)=2\pi r\Sigma_{mol}(r)V_{inflow}(r)$ (3)

with $\mathrm{\Sigma_{mol}}$ being the de-projected molecular mass surface
density, $V_{inflow}$ the CO inflow velocity and $r$ the galactocentric
distance. Figure 8 shows the spatially resolved $\dot{M}_{mol}$ and its
average radial profile. We find an average
$\dot{M}_{mol}\sim-20\mathrm{M_{\odot}/yr}$. For comparison, spiral arms
induce radial flows and radial velocities of the order of
$1\mathrm{M_{\odot}/yr}$ and $10$ km s-1 respectively, (e.g., Di Teodoro &
Peek, 2021).

Figure 8: Molecular mass inflow induced by the bar-like structure. Top panel:
2D distribution of $\dot{M}_{mol}$. Bottom panel: Average radial distribution
of $\dot{M}_{mol}$ in black and radial velocity profile of the inflow velocity
in red colors. Shadow black regions represent the standard deviation of
$\dot{M}_{mol}$ at each radial bin. A $\mathrm{SN>3}$ was applied to the
moment 0 map to exclude spurious CO detection.

##  6 Discussion and Conclusions

The combination of spatially resolved spectra provided by MUSE-VLT, with the
high spatial resolution from JWST and the ALMA CO $J=2\rightarrow 1$ data
allowed us to reveal a central kinematic oval distortion, as well as a small
scale elongated structure in the NIR continuum in NGC 1087. This elongated
structure is not evident in optical images, but only revealed with the NIRCam
F200W and F360M filters at $2\mu$m and $3\mu$m respectively, thanks to its
exquisite resolution. The optical counterpart of this bar is likely to be
affected by dust absorption and by the intense radiation from young stars. Its
preferential elongation suggests the presence of a non-axisymmetric
perturbation in the gravitational potential, in particular a bar-like type. At
first order, the 2$\mu$m light distribution of this structure can be described
with an exponential disk profile oriented at
$\phi_{bar}^{\prime}=305^{\circ}$, and the ellipticity of the elongated
structure is $\varepsilon^{\prime}=0.5$, (see Figure 2), with a half-major
axis length of $a_{bar}^{\prime}\sim 14\arcsec$ or $\sim$ 1.2 kpc in the
galaxy plane. Overall the shape of this structure resembles to a faint stellar
bar, although with remarkable differences with respect to conventional stellar
bars, in terms of bar-length, the lack of a SF nuclear ring and not clear
associated dust lanes along the bar axis. The same argument however applies
for an elongated bulge. The main reason this is not considered is due to the
large reservoirs of molecular gas observed, since bulges tend to be quiescent
structures with null or low SFRs (Shapiro et al., 2010). Furthermore, the [O
i]$/$H$\alpha$ and [O iii]$/$H$\beta$ line-ratios and large EWH$\alpha$ from
the bar (see Figure 1) indicate that the ionized gas is product from a recent
or ongoing SF event given its spatial coincidence with the CO emission.

The ionized and molecular gas kinematics and the stellar velocity, reveal
simultaneously the presence of an oval distortion in the velocity maps in
concordance with the extension of the faint bar. This evidences undoubtedly
the imprints of a non-axisymmetric perturbation in the potential, here
captured in different phases and enhanced in dynamically cold gas.

Kinematic models for a second order potential perturbation successfully
described the LoS velocities in the three velocity maps. We showed that this
is achieved just by fixing the orientation of the stream flow to the bar major
axis in the models, and the radial profiles of the noncircular velocities are
consistent with the residuals from circular rotation, with possible
differences in amplitudes due to assymmetric drift and other random motions
not considered in the model as we showed in Figure 7. Apart from the
elliptical motions considered in the bisymmetric models, radial inflows are
expected in bars, and these are observed to occur along the bar-dust lanes (Wu
et al., 2021; Sormani et al., 2023). Our implementation of an inflow model is
therefore physically motivated. An inflow of cold gas caused by a lose of
angular momentum could trigger the SF observed along the bar. Although a
fraction of ionized gas could arise by the shock with the dust lanes, the
contribution of shocks is expected to be minor. Moreover, the sign of the
radial flow suggests the gas is inflowing to the center at a maximum speed of
$50$ km s-1. This kinematic model yields the lowest rms despite being the
model with the lowest number of free variables. Statistically speaking, the
inflow scenario is favored over the bisymmetric model, although given the
complex behavior of gas in bars, the contribution of gas oval orbits can not
be ruled out. In fact both type of motions are expected to happen
simultaneously in bars (Athanassoula, 1992; Regan et al., 1999). To date,
there is no a 2D-kinematic model that includes such dynamics of bars. The gas
inflow model is in concordance with the observed enhancement in H$\alpha$ flux
along the bar. The inferred inflow velocities translate in molecular mass
inflow rates of $\dot{M}_{mol}=20~{}\mathrm{M_{\odot}/yr}$. This is larger
than the measured SFR(H$\alpha$) along the bar, thus we argue that a major
inflow of gas might be feeding the SF in the bar. In fact SF prevents the so-
called continuity problem (e.g., Simon et al., 2003; Maciejewski et al., 2012)
when often considering radial flow models. This large inflow rate is partially
due to the large inflow velocities; although this amplitude of non-circular
motions is observed in bars, the global inflow rate expected by hydrodynamical
models is lower (Athanassoula, 1992; Mundell & Shone, 1999).

To summarize, in this work we have taken the advantage of multiple archival
data to reveal a hidden non-axisymmetric structure in the optical, thanks to
the high angular resolution of the data. We were able to confirm the presence
of an oval distortion and successfully model it with a bar-like flow model. We
also show that the SFR in the bar could be explained by an inflow of gas along
the bar major-axis. Our results contribute to understand the overall picture
of non-axisymmetry in galaxies with the most sophisticated data so far.
Additionally, it highlights the importance of using the IR instead of optical
bands to detect stellar bars, which could increase the fraction of barred
galaxies detected in the local Universe.

## Acknowledgment

We thank the anonymous referee for their comments and suggestions which help
to improve the quality of this manuscript. CLC acknowledges support from
Academia Sinica Institute of Astronomy and Astrophysics. LL acknowledges the
Ministry of Science & Technology of Taiwan under the grant NSTC
112-2112-M-001-062 -.

Based on data obtained from the ESO Science Archive Facility with DOI(s):
https://doi.org/10.18727/archive/47 (catalog
https://doi.org/10.18727/archive/47). The specific PHANGS–JWST observations
analyzed can be can be found in MAST: https://doi.org/10.17909/9bdf-jn24
(catalog https://doi.org/10.17909/9bdf-jn24). This paper makes use of the
following ALMA data: ADS/JAO.ALMA#2018.1.01651.S.

## References

  * Alarie & Morisset (2019) Alarie, A., & Morisset, C. 2019, Rev. Mexicana Astron. Astrofis., 55, 377, doi: 10.22201/ia.01851101p.2019.55.02.21
  * Athanassoula (1992) Athanassoula, E. 1992, MNRAS, 259, 345, doi: 10.1093/mnras/259.2.345
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Proc. SPIE, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, 773508, doi: 10.1117/12.856027
  * Bacon et al. (2017) Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1, doi: 10.1051/0004-6361/201730833
  * Bundy et al. (2015) Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7, doi: 10.1088/0004-637X/798/1/7
  * Bushouse et al. (2022) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022, JWST Calibration Pipeline, 1.8.2, Zenodo, doi: 10.5281/zenodo.7229890
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Cheung et al. (2013) Cheung, E., Athanassoula, E., Masters, K. L., et al. 2013, ApJ, 779, 162, doi: 10.1088/0004-637X/779/2/162
  * Cid Fernandes et al. (2010) Cid Fernandes, R., Stasińska, G., Schlickmann, M. S., et al. 2010, MNRAS, 403, 1036, doi: 10.1111/j.1365-2966.2009.16185.x
  * Davies et al. (2016) Davies, R. L., Groves, B., Kewley, L. J., et al. 2016, MNRAS, 462, 1616, doi: 10.1093/mnras/stw1754
  * de Vaucouleurs et al. (1991) de Vaucouleurs, G., de Vaucouleurs, A., Corwin, Herold G., J., et al. 1991, Third Reference Catalogue of Bright Galaxies
  * Di Teodoro & Peek (2021) Di Teodoro, E. M., & Peek, J. E. G. 2021, ApJ, 923, 220, doi: 10.3847/1538-4357/ac2cbd
  * Emsellem et al. (2022) Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191, doi: 10.1051/0004-6361/202141727
  * Fathi et al. (2005) Fathi, K., van de Ven, G., Peletier, R. F., et al. 2005, MNRAS, 364, 773, doi: 10.1111/j.1365-2966.2005.09648.x
  * Gutkin et al. (2016) Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757, doi: 10.1093/mnras/stw1716
  * Hirota et al. (2009) Hirota, A., Kuno, N., Sato, N., et al. 2009, PASJ, 61, 441, doi: 10.1093/pasj/61.3.441
  * Kennicutt (1998) Kennicutt, Jr., R. C. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Kewley et al. (2006) Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961, doi: 10.1111/j.1365-2966.2006.10859.x
  * Kourkchi & Tully (2017) Kourkchi, E., & Tully, R. B. 2017, ApJ, 843, 16, doi: 10.3847/1538-4357/aa76db
  * Lacerda et al. (2022) Lacerda, E. A. D., Sánchez, S. F., Mejía-Narváez, A., et al. 2022, New A, 97, 101895, doi: 10.1016/j.newast.2022.101895
  * Lacerda et al. (2018) Lacerda, E. A. D., Cid Fernandes, R., Couto, G. S., et al. 2018, MNRAS, 474, 3727, doi: 10.1093/mnras/stx3022
  * Lang et al. (2020) Lang, P., Meidt, S. E., Rosolowsky, E., et al. 2020, ApJ, 897, 122, doi: 10.3847/1538-4357/ab9953
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3, doi: 10.1086/313233
  * Leroy et al. (2021a) Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021a, ApJS, 257, 43, doi: 10.3847/1538-4365/ac17f3
  * Leroy et al. (2021b) Leroy, A. K., Hughes, A., Liu, D., et al. 2021b, ApJS, 255, 19, doi: 10.3847/1538-4365/abec80
  * Lopez-Coba et al. (2021) Lopez-Coba, C., Lin, L.-H., & Sanchez Sanchez, S. F. 2021, arXiv e-prints, arXiv:2110.05095.  https://arxiv.org/abs/2110.05095
  * López-Cobá et al. (2022) López-Cobá, C., Sánchez, S. F., Lin, L., et al. 2022, ApJ, 939, 40, doi: 10.3847/1538-4357/ac937b
  * Maciejewski et al. (2012) Maciejewski, W., Emsellem, E., & Krajnović, D. 2012, MNRAS, 427, 3427, doi: 10.1111/j.1365-2966.2012.22040.x
  * Masters et al. (2011) Masters, K. L., Nichol, R. C., Hoyle, B., et al. 2011, MNRAS, 411, 2026, doi: 10.1111/j.1365-2966.2010.17834.x
  * Mazzalay et al. (2014) Mazzalay, X., Maciejewski, W., Erwin, P., et al. 2014, MNRAS, 438, 2036, doi: 10.1093/mnras/stt2319
  * Moiseev et al. (2015) Moiseev, A. V., Tikhonov, A. V., & Klypin, A. 2015, MNRAS, 449, 3568, doi: 10.1093/mnras/stv489
  * Mundell & Shone (1999) Mundell, C. G., & Shone, D. L. 1999, MNRAS, 304, 475, doi: 10.1046/j.1365-8711.1999.02330.x
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae and active galactic nuclei
  * Pence & Blackman (1984) Pence, W. D., & Blackman, C. P. 1984, MNRAS, 210, 547, doi: 10.1093/mnras/210.3.547
  * Pence et al. (1988) Pence, W. D., Taylor, K., Freeman, K. C., de Vaucouleurs, G., & Atherton, P. 1988, ApJ, 326, 564, doi: 10.1086/166117
  * Pilyugin & Grebel (2016) Pilyugin, L. S., & Grebel, E. K. 2016, MNRAS, 457, 3678, doi: 10.1093/mnras/stw238
  * Reese et al. (2007) Reese, A. S., Williams, T. B., Sellwood, J. A., Barnes, E. I., & Powell, B. A. 2007, AJ, 133, 2846, doi: 10.1086/516826
  * Regan et al. (1999) Regan, M. W., Sheth, K., & Vogel, S. N. 1999, ApJ, 526, 97, doi: 10.1086/307960
  * Rosolowsky, Erik (2022) Rosolowsky, Erik. 2022, PHANGS-JWST First Results, STScI/MAST, doi: 10.17909/9BDF-JN24
  * Sánchez et al. (2012) Sánchez, S. F., Kennicutt, R. C., Gil de Paz, A., et al. 2012, A&A, 538, A8, doi: 10.1051/0004-6361/201117353
  * Sánchez-Blázquez et al. (2011) Sánchez-Blázquez, P., Ocvirk, P., Gibson, B. K., Pérez, I., & Peletier, R. F. 2011, MNRAS, 415, 709, doi: 10.1111/j.1365-2966.2011.18749.x
  * Schinnerer (2021) Schinnerer, E. 2021, 3D Mapping of the Dynamic Feedback and Chemistry in the Star Formation Process across Galaxy Disks, European Southern Observatory (ESO), doi: 10.18727/ARCHIVE/47
  * Sellwood & Wilkinson (1993) Sellwood, J. A., & Wilkinson, A. 1993, Reports on Progress in Physics, 56, 173, doi: 10.1088/0034-4885/56/2/001
  * Shapiro et al. (2010) Shapiro, K. L., Falcón-Barroso, J., van de Ven, G., et al. 2010, MNRAS, 402, 2140, doi: 10.1111/j.1365-2966.2009.16111.x
  * Sheth et al. (2010) Sheth, K., Regan, M., Hinz, J. L., et al. 2010, PASP, 122, 1397, doi: 10.1086/657638
  * Simon et al. (2003) Simon, J. D., Bolatto, A. D., Leroy, A., & Blitz, L. 2003, ApJ, 596, 957, doi: 10.1086/378200
  * Sormani et al. (2023) Sormani, M. C., Barnes, A. T., Sun, J., et al. 2023, MNRAS, 523, 2918, doi: 10.1093/mnras/stad1554
  * Spekkens & Sellwood (2007) Spekkens, K., & Sellwood, J. A. 2007, ApJ, 664, 204, doi: 10.1086/518471
  * Stasińska et al. (2008) Stasińska, G., Vale Asari, N., Cid Fernandes, R., et al. 2008, MNRAS, 391, L29, doi: 10.1111/j.1745-3933.2008.00550.x
  * Sutherland et al. (2018) Sutherland, R., Dopita, M., Binette, L., & Groves, B. 2018, MAPPINGS V: Astrophysical plasma modeling code, Astrophysics Source Code Library, record ascl:1807.005.  http://ascl.net/1807.005
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409, doi: 10.1093/mnras/stw2231
  * Wong et al. (2004) Wong, T., Blitz, L., & Bosma, A. 2004, ApJ, 605, 183, doi: 10.1086/382215
  * Wu et al. (2021) Wu, Y.-T., Trejo, A., Espada, D., & Miyamoto, Y. 2021, MNRAS, 504, 3111, doi: 10.1093/mnras/stab1087


#  Draft
Hitting time for Markov decision process

Ruihcao Jiang, Javad Tavakoli, Yiqiang Zhao

###### Abstract

We define the hitting time for a Markov decision process (MDP). We do not use
the hitting time of the Markov process induced by the MDP because the induced
chain may not have a stationary distribution. Even it has a stationary
distribution, the stationary distribution may not coincide with the
(normalized) occupancy measure of the MDP. We observe a relationship between
the MDP and the PageRank. Using this observation, we construct an MP whose
stationary distribution coincides with the normalized occupancy measure of the
MDP and we define the hitting time of the MDP as the hitting time of the
associated MP.

##  1 Introduction

††footnotetext: The first version of this paper pointed out some issues in the
old version of [3]. The authors of [3] have then addressed these issues
according to our suggestions in a new version. We therefore updated our paper,
in which we removed contents related to these issues.

In [3], Gromov-Wasserstein distance was used to compare different MDPs. In
their approach, a distance function is needed on the space of state-action
pairs of the underlying MDP. They used Euclidean distance. However, it is not
entirely clear how to define the distance between actions. Moreover, it is
desirable that the distance reflects some useful information about the MDP.
For this reason, we propose to use the first-hitting time quasi-distance in
places of the distance function. The rationale is as follows.

  * •

The occupancy measure to the MDP is whats the stationary distribution to an
Markov process (MP);

  * •

The stationary distribution of an MP is simply the inverse of the mean of the
first-return time;

  * •

The first-return time is simply the diagonal element for the first-hitting
time matrix.

##  2 Hitting time of MDP

To define the first-hitting time of MDP, we first look at the occupancy
measure of MDP.

######  Definition 1 (Occupancy measure).

Let $(S,A,R,P,\gamma)$ be an MDP with the initial distribution $\rho_{0}(s)$
on $S$ and the stationary policy $\pi$, the occupancy measure $\rho:S\times
A\to\mathbb{R}$ is defined as follows.

$\rho_{\pi}(s,a)\coloneqq\sum_{t=0}^{\infty}\gamma^{t}\mathbb{P}\left(s_{t}=s,a_{t}=a|s_{0}\sim\rho_{0},\pi\right).$
(1)

We give a new interpretation of the occupancy measure. First, Tulcea’s theorem
(Proposition C.10 in [4]) says that a policy $\pi$ induces a Markov Process
(MP) on $S\times A$ with the transition probability

$P_{\pi}(s^{\prime},a^{\prime}|s,a)\coloneqq
P(s^{\prime}|s,a)\pi(a^{\prime}|s^{\prime})$ (2)

and the initial distribution $\rho_{0}(s,a)\coloneqq\rho_{0}(s)\pi(a|s)$. This
describes a process as follows: Given the state-action pair $(s,a)$, transit
to the next state $s^{\prime}$ according to the transition probability of the
original MDP; then immediately choose an action $a^{\prime}$ according to the
policy $\pi$.

Note that $\rho_{\pi}$ is not the stationary distribution (if it exists at
all) of $(S\times A,P_{\pi})$. First of all, $\rho_{\pi}$ is not a probability
measure as it sums to $\frac{1}{1-\gamma}$. Even after a normalization,
$\left(1-\gamma\right)\rho_{\pi}$ may differ from the stationary distribution.
Consider a scenario where a car moves straight on a straight road and if it
turns left or right, it will go off the road. Then none of the state-action
pairs related to steering will be occupied and hence $\rho_{\pi}$ is zero on
those pairs. However, for any stationary distribution, each component must be
strictly positive. We will see later that it does not matter whether $(S\times
A,P_{\pi})$ is ergodic or not. Treating $\rho_{\pi}$ and $\rho_{0}$ as
$1\times|S||A|$ vectors and $P_{\pi}$ as an $|S||A|\times|S||A|$ matrix, the
following equation is well-known.

$\rho_{\pi}=\left(I-\gamma P_{\pi}\right)^{-1}\rho_{0}.$ (3)

From Equation (3), we make the following key observation.

######  Observation (PageRank).

$(1-\gamma)\rho_{\pi}$ coincides with the personalized PageRank vector on
$S\times A$ with transition matrix $P_{\pi}$ and with probability $1-\gamma$
to restart with the initial distribution $\rho_{0}$.

Second, we look at the first-hitting time of an MP.

######  Definition 2 (First-hitting time of MP).

Let $\\{X_{t}\\}$ be an MP. The first-hitting time $T$ is the random variable
defined as follows.

$T_{ij}=\begin{cases}\inf_{t\in\mathbb{N}}\\{t\ |\ X_{0}=j,X_{1}\neq
i,\cdots,X_{t-1}\neq i,X_{t}=i\\}&\text{if}\ i\neq j\\\
0,&\text{otherwise.}\end{cases}$ (4)

If $i$ is never reached from $j$, $T_{ij}=+\infty$, where we have adopted the
convention that $\inf\emptyset=+\infty$.

By abuse of the notation, we call the expectation $\mathbb{E}T_{ij}$ of the
random variable $T_{ij}$ also the first-hitting time and denote it by
$T_{ij}$.

$T_{ij}\coloneqq\mathbb{E}T_{ij}=\begin{cases}\sum_{t=1}^{\infty}t\mathbb{P}(T_{ij}=t)&\text{if}\
i\neq j\\\ 0&\text{otherwise}.\end{cases}$ (5)

If the MP has transition matrix $P_{\pi}$, a one-step analysis shows that
$T_{ij}$ satisfies the following recursive relation.

$T_{ij}=\begin{cases}1+\sum_{k}T_{ik}P_{\pi}(k|j)&\text{if}\ i\neq j\\\
0&\text{otherwise}.\end{cases}$ (6)

By our Observation, we reduce the task of defining a first-hitting time for an
MDP $(S,A,R,P,\gamma)$ to that of defining a first-hitting time for a PageRank
$(S\times A,P_{\pi},1-\gamma)$. A PageRank can be described equivalently by a
new process $(S\times A,\tilde{P}_{\pi})$ with the transition matrix

$\tilde{P}_{\pi}(x^{\prime}|x)\coloneqq(1-\gamma)\rho_{0}(x^{\prime})+\gamma
P_{\pi}(x^{\prime}|x).$ (7)

$(S\times A,\tilde{P}_{\pi})$ is obviously Markovian. The initial condition
for this new process can be arbitrary, in particular, not necessarily
$\rho_{0}$ since it turns out that $(S\times A,\tilde{P}_{\pi})$ is ergodic
and its stationary distribution is given by $(1-\gamma)\rho_{\pi}$ (See
Proposition 1 and Corollary 1 of [1]). The ergodicity of $(S\times
A,\tilde{P}_{\pi})$ implies that the first-hitting time is finite on
supp$(\rho_{\pi})$. In this case, Equation (6) reads

$T_{ij}=\begin{cases}1+(1-\gamma)\sum_{k}T_{ik}\rho_{0}(k)+\gamma\sum_{k}T_{ik}P_{\pi}(k|j)&\text{if}\
i\neq j\\\ 0&\text{otherwise}.\end{cases}$ (8)

There is also a related quantity $L_{ij}$ satisfying the following recursive
relation.

$L_{ij}=\begin{cases}1+\gamma\sum_{k}L_{ik}P_{\pi}(k|j)&\text{if}\ i\neq j\\\
0&\text{otherwise}.\end{cases}$ (9)

The probabilistic interpretation of $L_{ij}$ is the expectation of the
discounted path length from $j$ to $i$ for the first time. $L_{ij}$ and
$T_{ij}$ are related by the following formula (Theorem 1(b) of [1]).

$T_{ij}=\frac{L_{ij}}{1-(1-\gamma)\sum_{k}L_{ik}\rho_{0}(k)}.$ (10)

######  Remark.

$L_{ij}$ is called the (1-$\gamma$)-discounted hitting time in [6] although
$L_{ij}$ is not the first-hitting time of $(S\times A,P_{\pi})$, whose first-
hitting time is given by Equation (6), nor of $(S\times A,\tilde{P}_{\pi})$,
whose first-hitting time is given by Equation (8). However, $L_{ij}$ is easier
to estimate empirically.

We have shown that

  * •

The occupancy measure of an MDP $(S,A,R,P,\gamma)$ with the initial
distribution $\rho_{0}$ coincides with the PageRank vector of the personalized
PageRank $(S\times A,P_{\pi},1-\gamma)$ with the initial distribution
$\rho_{0}$.

  * •

The PageRank vector of the personalized PageRank $(S\times
A,P_{\pi},1-\gamma)$ with the initial distribution $\rho_{0}$ coincides with
the stationary distribution of the ergodic MP $(S\times A,\tilde{P}_{\pi})$
with arbitrary initial distribution.

Based on the above facts, we propose the following definition.

######  Definition 3 (First-hitting time of MDP).

Let $(S,A,R,P,\gamma)$ be an MDP with a stationary policy $\pi$. Let $P_{\pi}$
be the induced transition matrix in Equation (2) and $\tilde{P}_{\pi}$ be the
modified transition matrix in Equation (7). The first-hitting time of the MDP
is defined as the first-hitting time of $(S\times A,\tilde{P}_{\pi})$ with its
expectation given by Equation (8).

Since we define the first-hitting time of MDP as the first-hitting time of
some MP, we immediately have the following.

######  Proposition 1.

The first-hitting time of an MDP gives rise to a quasi-distance
$T_{\pi}(i,j)\coloneqq T_{ij}$, i.e. a distance function without satisfying
the symmetry condition.

##  3 Hitting time quasidistance induces distance

It might be a concern that we use a quasi-distance, not a distance. In this
section, we show that applying $\mathcal{GW}$ to a quasi-metric measure spaces
still yield a distance (our Theorem 2).

We have constructed a triple $(S\times A,\rho_{\pi},T_{\pi})$ from the MDP
$(S,A,R,P,\gamma)$. Now we can apply the $\mathcal{GW}$ construction [5] to
define a quantity between two MDPs

######  Definition 4.

Let $(S_{1},A_{1},R_{1},P_{1},\gamma_{1})$ and
$(S_{2},A_{2},R_{2},P_{2},\gamma_{2})$ be MDPs. Denote by
$(\mathcal{X},\rho_{\mathcal{X}},T_{\mathcal{X}})$ and
$(\mathcal{Y},\rho_{\mathcal{Y}},T_{\mathcal{Y}})$ their induced triples,
respectively. Then

$\begin{split}&\mathcal{GW}((S_{1},A_{1},R_{1},P_{1},\gamma_{1}),(S_{2},A_{2},R_{2},P_{2},\gamma_{2}))\\\
\coloneqq&\mathcal{GW}((\mathcal{X},(1-\gamma_{1})\rho_{\mathcal{X}},T_{\mathcal{X}}),(\mathcal{Y},(1-\gamma_{2})\rho_{\mathcal{Y}},T_{\mathcal{Y}}))\\\
=&\min_{\mu}\frac{1}{2}\left[\sum_{\mathcal{X}\times\mathcal{Y}}\sum_{\mathcal{X}\times\mathcal{Y}}|T_{\mathcal{X}}(x,x^{\prime})-T_{\mathcal{Y}}(y,y^{\prime})|^{2}\mu(x,y)\mu(x^{\prime}.y^{\prime})\right]^{\frac{1}{2}},\end{split}$
(11)

where
$\mu\in\mathcal{M}((1-\gamma_{1})\rho_{\mathcal{X}},(1-\gamma_{2})\rho_{\mathcal{Y}})$
is the set of couplings, i.e.

$\mathcal{M}((1-\gamma_{1})\rho_{\mathcal{X}},(1-\gamma_{2})\rho_{\mathcal{Y}})=\left\\{\mu\left|\sum_{y\in\mathcal{Y}}\mu=(1-\gamma_{1})\rho_{\mathcal{X}}\
\text{and}\
\sum_{x\in\mathcal{X}}\mu=(1-\gamma_{2})\rho_{\mathcal{Y}}\right\\}\right.$
(12)

######  Remark.

It is necessary to normalize the occupancy measures. Otherwise, if
$\gamma_{1}\neq\gamma_{2}$, then no coupling exists since

$\sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}}\mu=\sum_{x\in\mathcal{X}}\rho_{\mathcal{X}}=\frac{1}{1-\gamma_{1}}\neq\frac{1}{1-\gamma_{2}}=\sum_{y\in\mathcal{Y}}\rho_{\mathcal{Y}}=\sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}}\mu.$

Different MDPs may give rise to the same triple. To eliminate this redundancy,
we introduce an equivalent relation on $\mathcal{MDP}$, the set of all MDPs.

######  Definition 5 (Equivalent relation on $\mathcal{MDP}$).

$\textup{MDP}_{1}$ and $\textup{MDP}_{2}$ are said to be equivalent, denoted
by $\textup{MDP}_{1}\sim\textup{MDP}_{2}$ if they give rise to the same
triple, i.e. there exists a bijection $\varphi:S_{1}\times A_{1}\to
S_{2}\times A_{2}$ s.t.
$T_{1}(x,x^{\prime})=T_{2}(\varphi(x),\varphi(x^{\prime}))$ and
$(1-\gamma_{1})\rho_{1}(x)=(1-\gamma_{2})\rho_{2}(\varphi(x))$ for all
$x,x^{\prime}\in S_{1}\times A_{1}$.

We have not yet called $\mathcal{GW}(\textup{MDP}_{1},\textup{MDP}_{2})$ a
distance. We know that $\mathcal{GW}$ is a distance function on the space of
metric-measure spaces after quotient out redundancy [5]. However, the first-
hitting time is only a quasi-distance, hence $(S\times
A,(1-\gamma)\rho_{\pi},T_{\pi})$ is not necessarily a metric-measure space. We
prove that $\mathcal{GW}$ is indeed a distance function on
$\mathcal{MDP}/\sim$.

######  Theorem 1.

$\mathcal{GW}$ is a distance on $\mathcal{MDP}/\sim$.

###### Proof.

Theorem 16 in [2] says that if $T_{\mathcal{X}}$ and $T_{\mathcal{Y}}$ are
measurable, then $\mathcal{GW}$ is a pseudo-distance. Hence,to show that
$\mathcal{GW}$ is a distance, we need to show that
$\mathcal{GW}(\textup{MDP}_{1},\textup{MDP}_{2})=0$ iff
$\textup{MDP}_{1}\sim\textup{MDP}_{2}$. The “if" part is trivial. We only need
to show
$\mathcal{GW}(\textup{MDP}_{1},\textup{MDP}_{2})=0\implies\textup{MDP}_{1}\sim\textup{MDP}_{2}$.
Let $\mu$ be a coupling s.t. Equation (11) evaluates to $0$. Then
$T_{\mathcal{X}}(x,x^{\prime})=T_{\mathcal{Y}}(y,y^{\prime})$ for all pairs of
$(x,y)$ and $(x^{\prime},y^{\prime})$ s.t. $\mu(x,y)>0$,
$\mu(x^{\prime},y^{\prime})>0$. Let $x=x^{\prime}$. Then
$T_{\mathcal{X}}(x,x)=T_{\mathcal{Y}}(y,y^{\prime})=0$. Since
$T_{\mathcal{Y}}$ is a quasi-distance, it implies that $y=y^{\prime}$. Hence
$\mu$ does not split mass i.e. $\mu$ induces a bijection $\varphi$ s.t.
$T_{\mathcal{X}}(x,x^{\prime})=T_{\mathcal{Y}}(\varphi(x),\varphi(x^{\prime}))$.
It remains to check measures. By the boundary condition Equation (12),

$\begin{split}(1-\gamma_{2})\rho_{\mathcal{Y}}(\varphi(x))&=\sum_{x\in\mathcal{X}}\mu(x,\varphi(x))\\\
&=\sum_{\varphi(x)\in\mathcal{Y}}\mu(x,\varphi(x))\\\
&=(1-\gamma_{1})\rho_{\mathcal{X}}(x)\end{split}$

as desired. ∎

######  Remark.

Intuitively, the definition of $\mathcal{GW}$ involves a minimization over all
couplings. If $T_{\pi}(x,x^{\prime})\neq T_{\pi}(x^{\prime},x)$, the
minimization will choose the smaller one.

######  Remark.

Theorem 16 in [2] only says that their $\mathcal{GW}$ is a pseudo-distance.
They explicitly constructs an example in their Remark 15 where their
$\mathcal{GW}$ fails to be a distance. This is because they assume nothing
about their network weight function $\omega$, except measurability. In
particular, $\omega(x,x)$ in their Remark 15 does not evaluate to zero.
However, for the first-hitting time, by definition $T(x,x)=0$.

##  4 Conclusion

We defined the hitting time for an MDP via a relationship between the MDP and
the PageRank. The hitting time is a quasi-distance. By applying $\mathcal{GW}$
construction, the hitting time quasi-distance induces a distance on the set of
discrete MDPs.

## References

  * [1] K. Avrachenkov, A. Piunovskiy and Y. Zhang  “Hitting Times in Markov Chains with Restart and their Application to Network Centrality”  In _Methodology and Computing in Applied Probability_ 20, 2018, pp. 1173–1188  URL: https://doi.org/10.1007/s11009-017-9600-5
  * [2] S. Chowdhury and F. Mémoli  “The Gromov–Wasserstein distance between networks and stable network invariants”  In _Information and Inference: A Journal of the IMA_ 8, 2019, pp. 757–787  URL: https://doi.org/10.1093/imaiai/iaz026
  * [3] Arnaud Fickinger, Samuel Cohen, Stuart Russell and Brandon Amos  “Cross-Domain Imitation Learning via Optimal Transport”  In _International Conference on Learning Representations_ , 2022  URL: https://openreview.net/forum?id=xP3cPq2hQC
  * [4] O. Hernández-Lerma and J. Lasserre  “Discrete-Time Markov Control Processes”, Applications of Mathematics  Springer, 1996  URL: https://doi.org/10.1007/978-1-4612-0729-0
  * [5] F. Mémoli  “Gromov–Wasserstein distances and the metric approach to object matching”  In _Foundations of computational mathematics_ 11, 2011, pp. 417–487  URL: https://doi.org/10.1007/s10208-011-9093-5
  * [6] P. Sarkar and G. Gordon  “Random Walks with Random Projections”  In _Workshop on Analyzing Networks and Learning with Graphs NIPS_ , 2009  URL: http://snap.stanford.edu/nipsgraphs2009/papers/sarkar-paper.pdf


# Robust Environment Perception for Automated Driving: A Unified Learning
Pipeline for Visual-Infrared Object Detection

Mohsen Vadidar1, Ali Kariminezhad1, Christian Mayr1, Laurent Kloeker2 and Lutz
Eckstein2 1The authors are with the Elektronische Fahrwerksysteme GmbH, 85080
Gaimersheim, Germany{mohsen.vadidar, ali.kariminezhad, christian5.mayr}@efs-
auto.de2The authors are with the research area Vehicle Intelligence &
Automated Driving, Institute for Automotive Engineering, RWTH Aachen
University, 52074 Aachen, Germany {laurent.kloeker<EMAIL_ADDRESS>aachen.de

###### Abstract

The RGB complementary metal-oxide-semiconductor (CMOS) sensor works within the
visible light spectrum. Therefore it is very sensitive to environmental light
conditions. On the contrary, a long-wave infrared (LWIR) sensor operating in
8-14 µm spectral band, functions independent of visible light.

In this paper, we exploit both visual and thermal perception units for robust
object detection purposes. After delicate synchronization and (cross-)
labeling of the FLIR [1] dataset, this multi-modal perception data passes
through a convolutional neural network (CNN) to detect three critical objects
on the road, namely pedestrians, bicycles, and cars. After evaluation of RGB
and infrared (thermal and infrared are often used interchangeably) sensors
separately, various network structures are compared to fuse the data at the
feature level effectively. Our RGB-thermal (RGBT) fusion network, which takes
advantage of a novel entropy-block attention module (EBAM), outperforms the
state-of-the-art network [2] by 10% with 82.9% mAP.

##  I INTRODUCTION

A statistical projection of traffic fatalities in the United States for the
first half of 2021 shows that an estimated 20,160 people died in motor vehicle
traffic crashes. This represents an increase of about 18.4 percent as compared
to 17,020 fatalities that were reported in the first half of 2020 [3]. Looking
at the fatal accidents of 2019 based on the time, one can see that there are
1,000 more fatal accidents during the night-time compared to the day-time [4].
Given less average traffic during the night-time, the importance of visibility
in dark is inevitable.

The number of publications on RGB-IR sensor fusion for multi-spectral object
detection in the automotive sector has increased within the past two years.
However, the lack of data in this research area is still noticeable. There are
two main sources of data, namely FLIR thermal dataset [1] and KAIST multi-
spectral pedestrian detection benchmark [5], which provide a dataset
containing IR and RGB pair images. FLIR mainly provides three classes car,
pedestrian, and bicycle, whereas KAIST only contains pedestrians.

(a) RGB Frame

(b) IR Frame

Figure 1: A sample of IR and RGB frames in heavy rain from the same scene: one
can notice that the pedestrians are not detected on the RGB frame and the cold
vehicle coming out of a parking lot is not recognized on the IR frame.
Realizing the complementary data in this multispectral setup motivated us to
fuse the information.

Nevertheless, the FLIR dataset comes only with IR labels. That introduces the
first challenge to researchers. Previously published papers [2], [6] and [7]
have made various objections to the dataset. For instance, the usage of
different camera resolutions at multiple instances ranging from 0.3 to 3.1 MP,
misalignment of resolution/aspect ratio, a vast difference between the field
of views and resolutions between the sensors, and no possibility to align the
RGB and IR frames without having a correlation matrix. These issues currently
made fusing the data on this dataset an unsolved problem.

Within this work, we address the FLIR dataset challenges and propose a
solution to properly utilize it. Further, we introduce a CNN architecture to
fuse RGB-IR data. The proposed network outperforms the state-of-the-art
networks for multi-spectral object detection purposes [2]. In section II, the
latest advancements and publications on RGB-IR sensor fusion networks will be
presented. Later in section III, a cross-labeling algorithm and pre-processing
pipeline will be introduced to provide suitable labels for RGB frames. In
section IV, we start with an RGB network and compare its detection performance
with a monospectral IR network. The best monospectral detector will be
considered as our baseline. Afterward, we proceed with the simplest form of
fusion and develop our proposed RGBT network step by step. Finally, sections V
and VI are dedicated to results, analysis and conclusion of this work.

##  II Related Works

In July 2020 Ravi Yadav et al. attempted fusing color and thermal images to
detect objects for self-driving applications [6]. They have used VGG16 [8] and
Faster-RCNN [9] as their encoder and detector, respectively. They achieved a
miss rate of 29% on the KAIST benchmark. However, they could not perform well
on the FLIR benchmark due to its complications. Later that month, Chaitanya
Devaguptapu et al. [10] have used the GAN framework to use the FLIR dataset
and produce RGB images from the IR frames. The idea was to borrow knowledge
from the data-rich pre-trained RGB domain and use it on IR images. Heng Zhang
et. al. introduced cyclic fuse and refine blocks for object detection [7] in
September 2020. The FSDD (Feature Fusion Single Shot Multibox Detector)
network [11] is used as their object detector and two VGG16 backbones are
implemented to work independently from each other for each sensor. No
information is given regarding the number of parameters or complexity of the
network. However, looking at the designed cyclic module, the number of its
occurrences, and considering the used VGG16, which is one of the most
expensive feature extractors with 138 million parameters, one can realize that
this network has more than 276 million parameters. In this work, a so-called
"well-aligned" version [7] of FLIR is released. However, as mentioned in the
paper, the misalignment problem made exploiting FLIR dataset for multispectral
object detection impractical. Therefore, the pre-processing part, could not
tackle all the issues from the dataset. We will have a look at a failed
attempt in the next section. Nevertheless, Cyclic Fuse-and-Refine Network
(CFR_3) could reduce the miss rate on the KAIST dataset and improve the mAP on
the FLIR benchmark to 72.39%. In November 2020 A. Sai Charan et al. attempted
to fuse IR and RGB data by giving weight to each sensor using attention
mechanism [12]. The attention module is placed in every residual block [13]
after the RELU activation function and first batch normalization. They have
tried both the late and mid-fusion approaches. Late fusion in this case is
done by pooling features of the thermal and RGB camera separately from their
respective region proposals and then concatenating them together. Comparing
the performance of the two mid and late fusions, mid-fusion has shown a better
result with 62.12% mAP. In January 2021 Heng Zhang et. al. proposed guided
attentive feature fusion (GAFF) [2]. This paper claims to be the first work
that is considering multi-spectral feature fusion as a sub-task in network
optimization. While bringing accuracy gains, GAFF has a lower computational
cost compared to common addition or concatenation methods. They reduced the
miss rate on the KAIST dataset by 2% and achieved an mAP of 72.9% on the FLIR
dataset. To the best of our knowledge, this is the highest mAP achieved using
the FLIR dataset. Thus, this work will be taken as the state-of-the-art
network on this benchmark. ZHANG et. al. in a comprehensive review studied
different methods and levels of fusion between visible and IR cameras for the
task of tracking [14]. Since the feature extraction process is the same in
tracking and object detection, this study can provide a good insight into
feature-level fusion using deep learning-based methods.

##  III Dataset and Data Pre-processing

This section is an important part of our work, without which the rest of our
research would not be possible.

###  III-A Dataset

Table I summarizes open-source datasets including IR frames. Most of the
datasets are concerned with pedestrian class only. Those which cover a wider
range of objects like the Dense dataset [15] confront other downsides such as
non-availability of RGB-IR data pairs, annotation issues, strong misalignment
between IR and RGB fields of view (FOV), and lack of extrinsic calibration
information (e.g. correlation matrix).

TABLE I: Available Datasets including IR Frames (NG stands for not given) Name | Year | Size (Frames) | Downsides
---|---|---|---
FLIR | 2020 | 13813 | IR label only
Dense | 2020 | 121500 | Narrow IR FOV
ZUT-FIR | 2020 | NG | IR data available only
SCUT-CV | 2019 | 211011 | IR data available only
|  |  | Pedestrian class only
CAMEL | 2018 | NG | Not on-board frames
Multispec. | 2017 | 7512 | Different Ref. Sys.
CVC-14 | 2016 | 3695 | Gray RGB
|  |  | Pedestrian class only
KAIST | 2015 | 95000 | Crowd label as human
|  |  | Pedestrian class only
CVC-09 | 2014 | 5990 | Pedestrian class only

###  III-B Data Pre-processing

To include the Vulnerable Road Users (VRU) as well as motorized users in our
study, we have decided to work with FLIR dataset [1]. The non-availability of
a dataset for fusing the IR and RGB frames, led us to put a great deal of
effort into developing a reliable cross-labeling algorithm, additional manual
labeling of the missing objects, and checking the data one by one, to make
sure the entire dataset is reliable and usable.

Sensor data fusion between an IR and an RGB camera was a challenge to
researchers since there is no open-source dataset providing synchronized and
labeled RGB-IR image pairs. Given that the ground truth in FLIR dataset is
provided for IR frames only, we introduce a method to cross-label the RGB
frames using IR annotations semi-automatically.

Figure 2: IR-RGB overlay: Union of the filed of views in green and the
transferred labels from IR to RGB frame in red box

Algorithm 1 demonstrates our approach for pre-processing the data. We first
convert the RGB frames to gray-scale images and extract the edges. IR and RGB
frames are matrices of size $640\times 512$ and $1800\times 1600$,
respectively. Thus, we introduce a scaling-parameter to resize the IR frame
gradually in a loop, until the IR extracted edges fall on RGB edges. The
cross-correlation value between the two matrices is a good metric to
understand the accuracy of the alignment. The highest value is obtained, when
the edges of two frames are completely aligned. Since the two frames are
different in size, the cross-correlation between them takes place in a
sliding-window manner. As figure 2 illustrates, vertical and horizontal
offsets exist between the RGB and IR frames. Figure 3 represents a comparison
between our result and the "well-aligned" version [7]. Ultimately, we utilize
6924, 1982, and 3960 pair images for training, validation, and test,
respectively.

Result: X_Offset, Y_Offset, Scaling_Factor

1 load both frames

2 Smooth both frames by a Gaussian filter

3 Detect the edges using Canny edge detector

4 Use binary threshold to set the edge pixels to one and all the other pixels
to zero

5

6Scaling_Factor = 2 # A constant to resize the IR frame

7 while  _Width_IR < Width_RGB_ do

8 Width_IR_new := Width_IR $\times$ Scaling_Factor

9 Resize IR binary frame with the Width_IR_new and the corresponding
Height_IR_new

10 Cross correlate (CC) the IR and RGB binary frames in a sliding window
manner and save the pixel coordinate of the maximum value from the correlation
matrix

11

12 if  _CC_Value > Best_CC_Value_ then

13 Update Best_CC_Value, Pixel_coordinates, Best_Scaling_Factor

14

15 end if

16 Scaling_Factor += 0.001

17

18 end while

19Best_X_Offset, Best_Y_Offset := Pixel_coordinates

20 return Best_X_Offset, Best_Y_Offset, Best_Scaling_Factor

Algorithm 1 Cross Labeling Algorithm

(a) "well-aligned" version

(b) our version

Figure 3: Comparison of two cross-labeling methods: the left image represents
the previous attempt [7] to process the FLIR dataset [1] and transfer the
labels from IR to RGB coordinate system. The right picture illustrates the
result of our proposed cross-labeling algorithm.

##  IV Network Structures

One of the most accurate networks on MSCOCO benchmark [16], namely Scaled-
Yolov4 [17] is selected as a starting point for this work. The performance of
each sensor will be evaluated separately and the best monospectral network
will be chosen as our baseline. The CSP (Cross-Stage-Partial-Connections) [18]
version of the network is used for both sensors. That is possible since the IR
sensor is delivering image-like data and the encoder is very good at picking
up the features and adapting itself to the IR sensors’ data texture (i.e.
domain adaptation).

To take the most of each monospectral backbone, we start with our vanilla
fusion. The designed network fuses the feature maps across scales as shown in
figures 4 and 5. The smallest scale is designed for large objects and the
biggest one detects the smaller objects. After concatenation, channels are
halved using $1\times 1$ convolution layers to prepare the feature maps for
the detector head.

The most crucial part of fusing visible and IR sensors is finding a way to
weight each backbone, depending on how informative the extracted features are.
As each sensor receives a different wavelength, we want to ensure that
complementary data from both spectrums increase the chance of detecting an
object as well as improving the robustness of the network. The goal is to
detect the objects, even in lack of strong emissions within visible-light
frequency. To do that, we integrate two Convolutional Block Attention Modules
(CBAM) [19], such that more weight is given to more informative features
depending on the input data (see figure 5).

Figure 4: RGBT Network Structure: Two identical encoders are used for feature
extraction from both RGB and IR frames. The feature-level fusion block
receives features in three scales from each encoder and fuses them. Then, the
data is passed to the detector, which is the same as CSP version of
Scaled_Yolov4. At last, three tensor outputs refer to three scales. Elements
of the fusion block are demonstrated in figure 5.

In what follows, we will discuss our proposed entropy-block attention module
(EBAM) and how it is a better substitute for the conventional CBAM [19]. In
this regard, we will proceed with channel attention, which can further be
generalized to spatial attention. Suppose, average and max pooling for feature
map (channel) $k$ are represented by $a_{k}$ and $m_{k}$, as:

$\displaystyle a_{k}$
$\displaystyle=\frac{1}{MN}\sum_{i}^{M}\sum_{j}^{N}X^{<k>}_{ij},$ (1)
$\displaystyle m_{k}$ $\displaystyle=\max_{i,j}\\{X^{<k>}_{ij}\\},$ (2)

respectively, where $\mathbf{X}^{<k>}$ is the matrix of $k$th feature map and
$X^{<k>}_{ij}$ is the $k$th feature map grid value for the grid located in
$i$th row and $j$th column. Further, we define the vector of average and max
poolings as

$\displaystyle\mathbf{a}$
$\displaystyle=[a_{1},\cdots,a_{k},\cdots,a_{K}]^{T},$ (3)
$\displaystyle\mathbf{m}$
$\displaystyle=[m_{1},\cdots,m_{k},\cdots,m_{K}]^{T}.$ (4)

According to the CBAM, the vectors $\mathbf{a}$ and $\mathbf{m}$ pass through
a shared multi-layer perceptron, for which the weights are optimized in the
backpropagation. Here, we assume a single layer for convenience in further
mathematical presentation. This can later be generalized to multi-layer
networks. Thus, the channel attention portion of the CBAM with a single-layer
perceptron delivers the following weight vector:

$\displaystyle\mathbf{c}_{\textrm{CBAM}}=\mathbf{W}(\mathbf{m}+\mathbf{a})+\mathbf{b},$
(5)

where the matrix of shared weights and the bias are represented by
$\mathbf{W}\in R^{K\times K}$ and $\mathbf{b}\in R^{K}$, respectively.
Therefore, the attention weight corresponding to the $k$th feature map can
readily be formulated as

$\displaystyle
c^{<k>}_{\textrm{CBAM}}=\mathbf{w}_{k}(\mathbf{m}+\mathbf{a})+b_{k},$ (6)

where $\mathbf{w}_{k}$ is the $k$th row of matrix $\mathbf{W}$. Interestingly,
one could observe that the average and max poolings are equally weighted due
to the shared layers among them. Consequently, the contribution of the average
and max poolings are not captured and only the sum of the values plays a role.
For instance, assume a particular feature map consists of multiple very high
values, while the majority of the grids contain low values. This way, the
average pooling downgrades the importance of that particular feature map.
However, due to the shared weights, this phenomenon can not be remedied in the
backpropagation by weight optimization. This motivates us to develop a more
robust and flexible attention mechanism as a replacement for CBAM.
Uncertainty of a data-driven feature is a valuable metric delivering a notion
for assessing how much extra information can still be mined from that feature.

$\displaystyle H_{k}=-\sum_{i=1}^{M}\sum_{j=1}^{N}P(X^{<k>}_{ij})\log
P(X^{<k>}_{ij}),$ (7)

where $P(X^{<k>}_{ij})$ are the probabilities delivered by the softmax
operation as,

$\displaystyle
P(X^{<k>}_{ij})=\frac{e^{X^{<k>}_{ij}}}{\sum_{i=1}^{M}\sum_{j=1}^{N}e^{X^{<k>}_{ij}}}.$
(8)

Notice that, equal probabilities maximizes the entropy.

$\displaystyle P(X^{<k>}_{ij})=\frac{1}{MN},\quad\forall i,j.$ (9)

It is important to discern that, the max entropy is reached if the max pooling
equals the average pooling of the feature map $k$, i.e., $a_{k}=m_{k}$, which
holds if:

$\displaystyle
X^{<k>}_{ij}=\frac{1}{\sum_{i}^{M}\sum_{j}^{N}X^{<k>}_{ij}},\forall i,j.$ (10)

This reflects the fact that in an extreme case the sum of max and average,
i.e., $m_{k}+a_{k}$, corresponds with the entropy upper-bound, i.e.,
$H^{\textrm{max}}_{k}$. This extreme case, i.e., $a_{k}=m_{k}$ can not be
differentiated by CBAM, since it considers only the sum. However, entropy
captures the interaction between max and average as follows:

$\displaystyle\lim_{m_{k}-a_{k}\rightarrow 0}H_{k}$
$\displaystyle=H^{\textrm{max}}_{k},$ (11)
$\displaystyle\lim_{m_{k}-a_{k}\rightarrow\infty}H_{k}$ $\displaystyle=0.$
(12)

That means, as the max pool gets closer to the average pool, the entropy
maximizes and as it gets higher value than the average pool, the entropy
approaches zero. It is important to mention that, the entropy of a random
process does not reduce except by providing information. Furthermore, notice
that more data does not necessarily mean more information. This observation
motivated us to design a network, which focuses on informative features in the
fusion process. This focus is represented by strongly weighting the highly-
informative feature maps. We call this new attention mechanism entropy-block
attention module (EBAM), for which the $k$th channel attention weight is given
by:

$\displaystyle c^{<k>}_{\textrm{EBAM}}=\mathbf{v}_{k}\mathbf{h}+b_{k},$ (13)

where $\mathbf{h}=[H_{1},\cdots,H_{k},\cdots,H_{K}]^{T}$. Note that
$\mathbf{v}_{k}$ is the channel attention weight vector designated for the
$k$th feature map optimized in the backpropagation. In what follows we
generalize the EBAM to capture both channel and spatial attention.
Let’s define the intermediate feature map $F\in R^{C\times H\times W}$ as a
tensor consisting of $\mathbf{X}^{<k>},\ \forall k$. That means,
$F=[\mathbf{X}^{<1>},\cdots,\mathbf{X}^{<k>},\cdots,\mathbf{X}^{<K>}]$. This
feature map tensor is the input to our attention module, from which attention
weights are derived. Hence, the generalized self-attention can be written as:

$W(F)=\alpha(F)\>\otimes\>F,$ (14)

where $\otimes$ denotes element-wise multiplication. $C,H$ and $W$ refer to
channel, height and width, respectively and $\alpha(F):R^{C\times H\times
W}\rightarrow R^{C\times H\times W}$ is an attention function of the feature
map $F$.

Figure 5: Upper block: Vanilla fusion is referred to a simplistic approach to
fuse the data by just concatenating the extracted features in each scale and
selecting the most valuable information by a convolution layer. Middle block:
CBAM is added as an attention module. Lower block: EBAM is designed to select
informative features more efficiently.

EBAM sequentially infers a 1D channel attention map $A_{c}\in R^{C\times
1\times 1}$ and a 2D spatial attention $A_{s}\in R^{1\times W\times H}$ as
illustrated in figure 6. The overall attention process can be summarized as:

$\displaystyle F^{\prime}=A_{c}(F)\otimes F,$ (15) $\displaystyle
F^{\prime\prime}=A_{s}(F^{\prime})\otimes F^{\prime},$

Figure 6: EBAM: detailed illustration of entropy-block attention module, which
contains channel and spatial domains. MLP refers to multilayer perceptron.

In the channel attention by giving more weight to more uncertain feature maps
(assuming each feature map is represented by a random variable), the network
is signaled to extract more information by prioritizing the gradient direction
in backpropagation. Therefore, the network is forced to gain more information
(i.e. kernel activation) from different kernels. As a matter of fact, the
proposed attention module operating only in the channel domain, showed a
relatively better performance than a CBAM including both spatial and channel
attention. One can calculate the explained channel attention as follows:

$\displaystyle A_{c}(F)$ $\displaystyle=\sigma(MLP(H(F))),$ (16)
$\displaystyle=\sigma(W_{0}(F^{C}_{H})),$

where $\sigma$ denotes the sigmoid function, $H$ is the entropy matrix, $F$
stands for feature map and $W_{0}$ represents the MLP weights.

The more the entropy of a random process, the more uncertain we are about the
outcome. In the case of the attention module in the spatial domain, each grid
pattern can contribute to detecting an object. As the training process goes
on, the grids, which refer to an object will have lower entropy (i.e. lower
uncertainty). Hence, if we use the entropy values directly, more weight will
be given to grids with high entropy, which would be the background information
in this case. Thus, more weights are intended to be allocated for the grids
with less entropy values through the feature maps.

Therefore, the entropy-based spatial attention can be denoted as follows:

$\displaystyle A_{s}(F^{\prime})$
$\displaystyle=\sigma(Conv.(1-\dfrac{H_{ij}}{\max_{\\{i,j\\}}H_{ij}})),$ (17)
$\displaystyle=\sigma(W_{1}(F^{S}_{H})),$

The networks are trained to the point, where an overfitting effect appeared.
The hyperparameters are kept unchanged as the original Scaled-Yolov4 [17]
implementation.

##  V Results and Analysis

In this chapter, the results of our pre-processing, baseline, and fusion
networks are presented. The processed data and the code can be found in the
following GitHub link: https://github.com/SamVadidar/RGBT.

###  V-A Data Pre-processing

Exploiting our proposed pre-processing step (i.e. algorithm 1), not only the
accuracy of the bounding boxes is improved, but more data (i.e. 13,848 images)
could be recovered from the original dataset. That means 34.6% more data with
respect to what Heng Zhang’s version [7] could use in training. Thus, a
dataset is produced, where the labels can be referred to both frames, which
are $640\times 512$ pixels. After running the pre-processing algorithm on the
original dataset, one receives 6924, 1982, and 3960 paired images for train,
validation, and test sets, respectively. Since the intention is to compare a
CMOS stand-alone sensor with a CMOS-IR sensor pair, the dataset is split into
night and day subsets. Therefore, one can test and evaluate a trained network
based on these two scenarios.

###  V-B Studied Networks

To reduce the computational costs, all the network trainings and their
comparisons will be done using $320\times 320$ image size and only the final
network will be fine-tuned with the original images ($640\times 512$ pixels).
Generally speaking, car class being the largest object among the three, shows
the best mAP. On the contrary, the person class being the smallest in size is
the most challenging object to detect.

TABLE II: Summary of studied networks: Fusion+H_C refers to a fusion network with entropy-block attention modules in the channel domain only. RGBT contains the full network structure as discussed in chapter IV. CFR [7] and GAFF [2] illustrate the two best available mAPs (NG stands for not given) on the FLIR dataset [1]. RGBT* is trained and tested using previously published aligned version [7] of the FLIR dataset. Network |<EMAIL_ADDRESS>| Param.
---|---|---
| Person | Bicycle | Car | Overall |
RGB Baseline | 39.6% | 50.4% | 79.4% | 56.6% | 52.5 M
IR Baseline | 49.6% | 54.9% | 84.4% | 63.0% | 52.5 M
Vanilla Fusion | 56.9% | 56.7% | 82.0% | 65.2% | 81.8 M
Fusion+CBAM | 57.6% | 60.5% | 83.6% | 67.2% | 82.7 M
Fusion+H_C | 62.6% | 65.9% | 86.0% | 71.5% | 82.7 M
RGBT | 63.7% | 67.1% | 86.4% | 72.4% | 82.7 M
CFR(640) | 74.4% | 57.7% | 84.9% | 72.3% |  $>$276 M
GAFF(640) | NG | NG | NG | 72.9% | NG
RGBT*(640) | 85.3% | 64.0% | 88.6% | 79.3% | 82.7 M
RGBT(640) | 80.1% | 76.7% | 91.8% | 82.9% | 82.7 M

A very interesting finding is that the infrared camera outperforms the CMOS
sensor even as a stand-alone sensor. Its performance in table II shows a 6.4%
overall improvement in the mAP metric using the same training pipeline, same
network, and an equal number of epochs. One can observe the improvement in all
classes, especially for the pedestrian class with a 10% jump in mAP. High
contrast in the frames between the objects and surroundings, very high sensor
sensitivity of less than 50 mK as well as an operating temperature of -40°C to
+85°C, delivers a data-rich infrared frame for feature extraction. Later in
this chapter, some qualitative and quantitative analyses will be illustrated
to support this statement.

The first attempt for fusing the two backbones (i.e. Vanilla Fusion) has shown
a 2.2% improvement in the overall mAP compared to the IR network (See table
II). The boldest improvement among the classes belongs to the pedestrian
class, which has experienced a drastic improvement of 7.3% with respect to the
IR network and 17.3% to the RGB network considering the mAP metric. Overall
8.6% mAP rise compared to the RGB network by adding no additional module,
motivated us to use the CBAM attention module [19]. By the way the CBAM is
implemented, each beneficial sensor data is enhanced and downsides are
compensated by focusing on the right and reliable features automatically.
Therefore, an additional performance gain of 2.2% compared to the vanilla
network is achieved and the bicycle class mAP is improved by 3.8%. After
designing the entropy-block, despite spatial attention absence, the Fusion+H_C
network could outperform the Fusion+CBAM network. One can see 4.3% mAP
improvement overall, and 5.0% for pedestrian class.

The proposed RGBT Network contains two complete (i.e. channel and spatial)
attention modules. The spatial attention added 0.9% on top of the channel-only
attention network and pushed the overall mAP up to 72.4%. The network shows a
noticeable overall mAP improvement to RGB and IR baselines with 15.8% and
9.4%, respectively. This shows that a fusion network can have an immense
impact on improving the mAP. Since the RGBT network performs better than both
IR and RGB object detectors, one can observe the importance of an efficient
feature-fusion by selecting the best features from each encoder through the
entropy-block attention module. This work could outperform the results of the
CFR network [7], by only utilizing half of the pixel resolution. We also have
trained and tested our proposed network with the same aligned version of the
FLIR dataset [7], which CFR and GAFF used. This helps to differentiate the
contribution of our preprocessing method from the proposed network. As the
table II illustrates, 3.6% mAP improvement is solely achieved by better
preprocessing the data. The test-set of the "well-aligned" dataset has almost
46% less number of pedestrian instances. Fewer test data and relatively easier
scenarios might be the underlying reason behind the higher mAP (i.e. 85.3%)
for this class. Considering the number of parameters and what the CFR paper
[7] reports, two VGG-16 networks [8] are used as separate feature extractors
of each sensor. Hence, knowing that each VGG-16 has around 138 million
parameters, one can find out that the number of the parameters in the state-
of-the-art network is greater than 276 million. Looking at the table II, it is
demonstrated that the RGBT network could outperform the CFR network by using
at least 70% less number of parameters. Therefore, not only the mAP is
improved in this work, but the complexity of the network is reduced
drastically. This will cause a considerably faster inference speed compared to
the CFR network. The GAFF, which showed slightly (i.e. 0.6%) better mAP than
CFR network, does not report the per class mAPs and the total number of
parameters of the complete network. Due to limited resources and long training
time, only our proposed network is trained with the original image size. Table
II shows that using the full resolution, one can reach 82.9% mAP with 10%
performance gain compared to the state-of-the-art network.

The results using the original images are compared to other networks in figure
7. 10.5% raise in the mAP due to a higher input resolution, which indicates
the importance of the pixel per object value. Working with FLIR dataset [1],
one can observe that it contains very small and challenging objects,
especially in the pedestrian class. Therefore, when the number of input
pixels, which are fed to the network is doubled (i.e. $640\times 640$), the
overall outcome elevates significantly.

Figure 7: Qualitative comparison of three scenarios: Columns from left to
right: RGB, IR, Fusion with CBAM, and Fusion with EBAM. Note that the fusion
results are represented using TGB (i.e. thermal, green, blue) instead of RGB
colormap.

Looking at figure 7, it is to be noted that the detected objects on the
sidewalk in the first two scenarios (i.e. first two rows) were not labeled.
That means the network is using the extracted features from other instances of
the pedestrian class to classify and localize the objects. Since in the
training process, detected objects, which do not have a corresponding label
are punished by the loss function, the performance of the network in this
scene is admirable. The third scenario (i.e. third-row) illustrates a group of
heavily occluded pedestrians sitting around a table. The night scenarios are
the most challenging scenarios for the RGB cameras and in this case, the
target object was neither labeled nor provided a sufficient amount of pixel
information. Despite these challenges, the RGBT network having only $320\times
320$ resolution could detect the pedestrians as well as the bicycles in front
of them.

##  VI Conclusion and Outlook

In general, cameras are one of the most informative sensors on the vehicle.
However, they are highly sensitive to light and weather conditions. Such
disturbances affect their performance severely. In order to introduce
robustness to the system, the functionality of an infrared camera as a
complementary sensor is investigated. Two groups of researchers [2], [6] and
[7] have considered data pre-processing and preparation of FLIR dataset [1] as
one of the major challenges and referred to it as an unsolved problem. High-
quality data being the core of machine learning approaches motivated us to
address this issue with our cross-labeling algorithm. Hence, 34.6% more data
is retrieved with respect to the previous approach [7]. Considering the trend
in recent publications, the rally was to gain a single additional mAP
improvement, whereas thanks to the pre-processing pipeline and the entropy-
block attention module, the proposed RGBT network could outperform the state-
of-the-art network [2] by 10% mAP. Undoubtedly, the availability of higher
quantity and quality labeled data is essential for further investigation in
this domain. In our next step, we are planning to combine our RGB-Radar
network with the proposed RGBT network to extract depth information from the
detected objects. We also have seen that Radar can improve the mAP
considerably. At last, the network will be optimized by TensorRT for our
target hardware and run on our demonstrator vehicle to perform real-world
tests.

## ACKNOWLEDGMENT

The authors would like to thank EFS MLOps team, specially Mr. Sebastian Balz
for his unceasing support, without which the training and evaluation of the
networks would not be possible.

## References

  * [1] T. F. LLC, “Free flir thermal dataset for algorithm training,” https://www.flir.com/oem/adas/adas-dataset-form/, 2020, accessed: 01.05.2021.
  * [2] H. Zhang, E. Fromont, S. Lefevre, and B. Avignon, “Guided attentive feature fusion for multispectral pedestrian detection,” in _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV)_ , January 2021, pp. 72–80, accessed: 13.06.2021.
  * [3] N. N. H. T. S. Administration), “Early estimate of motor vehicle traffic fatalities for the first half (jan–june) of 2021,” https://crashstats.nhtsa.dot.gov/Api/Public/ViewPublication/813199, 2021, accessed: 08.11.2021.
  * [4] ——, “Fatality analysis reporting system (fars) auxiliary datasets analytical user’s manual 1982-2019,” https://crashstats.nhtsa.dot.gov/Api/Public/ViewPublication/813071, 2021, accessed: 20.08.2021.
  * [5] S. Hwang, J. Park, N. Kim, Y. Choi, and I. S. Kweon, “Multispectral pedestrian detection: Benchmark dataset and baselines,” in _Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015, accessed: 15.05.2021.
  * [6] R. Yadav, A. Samir, H. Rashed, S. Yogamani, and R. Dahyot, “Cnn based color and thermal image fusion for object detection in automated driving,” 07 2020, accessed: 05.05.2021.
  * [7] H. Zhang, E. Fromont, S. Lefèvre, and B. Avignon, “Multispectral fusion for object detection with cyclic fuse-and-refine blocks - icip 2020 - ieee international conference on image processing,” Abou Dabi, United Arab Emirates, Oct. 2020, pp. 1–5, accessed: 12.05.2021. [Online]. Available: https://hal.archives-ouvertes.fr/hal-02872132
  * [8] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv 1409.1556_ , 09 2014, accessed: 05.05.2021.
  * [9] S. Ren, K. He, R. B. Girshick, and J. Sun, “Faster R-CNN: towards real-time object detection with region proposal networks,” _CoRR_ , vol. abs/1506.01497, 2015, accessed: 21.05.2021. [Online]. Available: http://arxiv.org/abs/1506.01497
  * [10] C. Devaguptapu, N. Akolekar, M. M Sharma, and V. N Balasubramanian, “Borrow from anywhere: Pseudo multi-modal object detection in thermal imagery,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR) Workshops_ , June 2019, accessed: 13.05.2021.
  * [11] Z. Li and F. Zhou, “FSSD: feature fusion single shot multibox detector,” _CoRR_ , vol. abs/1712.00960, 2017, accessed: 18.06.2021. [Online]. Available: http://arxiv.org/abs/1712.00960
  * [12] A. S. Charan, M. Jitesh, M. Chowdhury, and H. Venkataraman, “Abifn: Attention-based bi-modal fusion network for object detection at night time,” _Electronics Letters_ , vol. 56, pp. 1309–1311(2), November 2020, accessed: 13.06.2021. [Online]. Available: https://digital-library.theiet.org/content/journals/10.1049/el.2020.1952
  * [13] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” _CoRR_ , vol. abs/1512.03385, 2015, accessed: 05.08.2021. [Online]. Available: http://arxiv.org/abs/1512.03385
  * [14] X. Zhang, P. Ye, H. Leung, K. Gong, and G. Xiao, “Object fusion tracking based on visible and infrared images: A comprehensive review,” _Information Fusion_ , vol. 63, pp. 166–187, 2020, accessed: 21.08.2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1566253520302657
  * [15] M. Bijelic, T. Gruber, F. Mannan, F. Kraus, W. Ritter, K. Dietmayer, and F. Heide, “Seeing through fog without seeing fog: Deep multimodal sensor fusion in unseen adverse weather,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2020, accessed: 30.06.2021.
  * [16] T.-Y. Lin, M. Maire, S. Belongie, L. Bourdev, R. Girshick, J. Hays, P. Perona, D. Ramanan, C. L. Zitnick, and P. Dollár, “Microsoft coco: Common objects in context,” 2014, cite arxiv:1405.0312Comment: 1) updated annotation pipeline description and figures; 2) added new section describing datasets splits; 3) updated author list 4) Accessed: 18.06.2021. [Online]. Available: http://arxiv.org/abs/1405.0312
  * [17] C.-Y. Wang, A. Bochkovskiy, and H.-Y. M. Liao, “Scaled-YOLOv4: Scaling cross stage partial network,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2021, pp. 13 029–13 038, accessed: 23.06.2021.
  * [18] C.-Y. Wang, H.-Y. M. Liao, Y.-H. Wu, P.-Y. Chen, J.-W. Hsieh, and I.-H. Yeh, “Cspnet: A new backbone that can enhance learning capability of cnn,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops_ , June 2020, accessed: 13.06.2021.
  * [19] S. Woo, J. Park, J.-Y. Lee, and I. S. Kweon, “Cbam: Convolutional block attention module,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , September 2018, accessed: 11.07.2021.