# A gap-sharing planet pair shaping the crescent in HD 163296: a disk sculpted
by a resonant chain

Juan Garrido-Deutelmoser Instituto de Astrofísica, Pontificia Universidad
Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile
Núcleo Milenio de Formación Planetaria (NPF), Chile  Cristobal Petrovich
Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña
Mackenna 4860, 782-0436 Macul, Santiago, Chile  Millennium Institute for
Astrophysics, Chile  Carolina Charalambous naXys, Department of Mathematics,
University of Namur, Rue de Bruxelles 61, 5000 Namur, Belgium  Viviana V.
Guzmán Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av.
Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile  Núcleo Milenio de
Formación Planetaria (NPF), Chile  Ke Zhang Department of Astronomy,
University of Wisconsin-Madison, 475 N. Charter Street, Madison, WI 53706, USA

###### Abstract

ALMA observations of the disk around HD 163296 have resolved a crescent-shape
substructure at around 55 au, inside and off-center from a gap in the dust
that extends from 38 au to 62 au. In this work we propose that both the
crescent and the dust rings are caused by a compact pair (period ratio $\simeq
4:3$) of sub-Saturn-mass planets inside the gap, with the crescent
corresponding to dust trapped at the $L_{5}$ Lagrange point of the outer
planet. This interpretation also reproduces well the gap in the gas recently
measured from the CO observations, which is shallower than what is expected in
a model where the gap is carved by a single planet. Building on previous works
arguing for outer planets at $\approx 86$ and $\approx 137$ au, we provide
with a global model of the disk that best reproduces the data and show that
all four planets may fall into a long resonant chain, with the outer three
planets in a 1:2:4 Laplace resonance. We show that this configuration is not
only an expected outcome from disk-planet interaction in this system, but it
can also help constraining the radial and angular position of the planet
candidates using three-body resonances.

protoplanetary disks — planet–disk interactions — hydrodynamics — planets and
satellites: dynamical evolution and stability — radiative transfer

††software: Fargo3D (Benítez-Llambay et al., 2019), Numpy (van der Walt et
al., 2011), Matplotlib (Hunter, 2007). Rebound (Rein & Liu, 2012), ReboundX
(Tamayo et al., 2019), RADMC3D (Dullemond et al., 2012).

##  1 Introduction

Substructures are ubiquitous in protoplanetary disks, particularly in the dust
density distribution exhibited by high angular resolution observations
(Andrews, 2020; Bae et al., 2022). The Atacama Large Millimeter Array (ALMA)
has revealed a variety of substructures, whereas a large population of rings
and gaps are shown in continuum observations, to a lesser extent, in molecular
line emissions (e.g., van der Marel et al. 2019).

The advances in spatial resolution have been able to resolve non-axisymmetric
substructures within gaps, including systems such as PDS 70 (Benisty et al.,
2021), HD 163296 (Isella et al., 2018), HD 100546 (Pérez et al., 2020), HD
97048 (Pinte et al., 2019), and LkCa 15 (Long et al., 2022). These
substructures may be due to embedded planets induced by gravitational
interactions (e.g., Bae et al. 2022), which vary depending on the emission
shape. Point-like emissions are generally associated with an accreting planet
surrounded by a circumplanetary disk (CPD) (Perez et al., 2015; Szulágyi et
al., 2018), while crescent shapes may be related to the stable Lagrange points
$L_{4}$ and $L_{5}$ of a star-planet system (Rodenkirch et al., 2021; Long et
al., 2022) or vortices. These hypotheses frequently assume that a single
substructure is caused by a single planet, often in a Jovian mass regime.
However, we have recently shown in Garrido-Deutelmoser et al. (2022), that a
pair of lower-mass and gap-sharing planets can sculpt compact and/or elongated
vortices within the gap that last for several thousands orbits.

The disk surrounding HD 163296 contains ringed structures in the mm-continuum
(Isella et al., 2018) and several molecular tracers (Law et al., 2021; Zhang
et al., 2021). In particular, inside the dust density gap that extends from 38
to 62 au, a crescent-shaped substructure resides at around 55 au (Teague et
al., 2021). Recently, it was suggested that the emission comes from the dust
trapped around the stable point $L_{5}$ of a Jupiter mass planet orbiting at
48 au (Rodenkirch et al., 2021). Even though this method seems to reproduce
the broad features of the dust continuum distribution, two aspects remain
unclear:

  1. 1.

the crescent feature resides at 55 au, off-center from the dust gap, while the
$L_{5}$ point is co-orbital to the planet at 48 au. Varying the planet’s
eccentricity to account for this shift is unlikely to help as the stability of
the crescent is damaged, leading to its prompt disruption.

  2. 2.

the Jupiter-mass planet needed to retain observable amounts of dust at $L_{5}$
and open a wide enough gap in the dust, is expected to open a deep gap111An
increase in the local viscosity to produce a much shallower gap comes at the
expense of reducing the lifetime of the $L_{5}$ crescent, or even prevent its
formation in the first place (Rodenkirch et al., 2021). (Duffell, 2020). This
prediction disagrees with the recent results provided by Zhang et al. (2021),
who found that the dust density gap has a corresponding CO gap $\sim$ 10 times
shallower than the predictions involving a Jupiter. The local gas depletion
depends on the planetary mass ($\propto m_{\rm p}^{-2}$) (Kanagawa et al.,
2015), whereby opting for a lower mass planet to carve a shallower gap, may
not produce a sufficient gravitational interaction to enforce the dust
trapping at $L_{5}$.

In this work, we propose a that a compact pair of sub-Saturn-mass planets can
solve these issues, simultaneously accounting for the dust emission (the
shifted crescent and dust rings) and shallow gap in the CO. This scenario is
largely motivated by our recent work in Garrido-Deutelmoser et al. (2022)
where we showed that a compact pair of gap-sharing planets generally lead to
nonaxisymmetric substructures like that observed in HD 163296.

##  2 Setup

The hydrodynamics simulations and radiative transfer calculations in this work
largely follow the scheme in Rodenkirch et al. (2021) and are only briefly
summarized here. We carried out 2D hydrodynamic simulations using the FARGO3D
multifluid code (Benítez-Llambay et al., 2019; Masset, 2000; Weber et al.,
2019) to produce gas and dust density distribution for a fiducial disk model.
The resulting density maps are read into the RADMC3D code (Dullemond et al.,
2012) to calculate the radiative transfer image at $\lambda\sim$ 1.25 mm. We
use this image and the HD 163296 template (that contain all the technical
properties of the observation) in the
SIMIO222https://www.nicolaskurtovic.com/simio package to achieve the synthetic
ALMA observation comparable to the dust continuum observation provided by
Isella et al. (2018).

###  2.1 Hydrodynamic Simulations

The initial surface density profiles for the gas (sub-index g) and dust (sub-
index d) are given by

$\Sigma_{\rm g/d}=\Sigma_{{\rm g/d,0}}\left(\frac{r}{r_{0}}\right)^{-0.8}{\rm
exp}\left[-\left(\frac{r}{r_{\rm c,g/d}}\right)^{\gamma_{\rm g/d}}\right],$
(1)

where we set $r_{0}=48$ au, the initial surface density $\Sigma_{\rm
g,0}=37.4\;\rm gr\;cm^{-2}$, the cut-off radius $r_{\rm c,g}=165$ au, and the
exponent $\gamma_{\rm g}=1$. Similarly, for the dust we set $\Sigma_{\rm
d,0}=\Sigma_{\rm g,0}/100=0.374\;\rm gr\;cm^{-2}$, $r_{\rm c,d}=90$ au, and
$\gamma_{\rm d}=2$. We include the evolution for five independent dust
species. These grains have sizes in cm of $0.02$, $0.071$, $0.13$, $0.26$, and
$1.92$. We use an aspect ratio of $h(r)=h_{0}(r/r_{0})^{f}$ with $h_{0}=0.05$
and a flaring index $f=0.25$, which implies a mid-plane temperature profile
described by $T=25(r/r_{0})^{-0.5}$ K. This setup coincides with that from
Rodenkirch et al. (2021).

The disk extends from $r_{\rm in}=5$ au to $r_{\rm out}=197$ au, implying an
initial disk mass of $\sim 0.15\;M_{\rm\odot}$. The computational domain is
composed of $n_{\rm r}=512$ logarithmically spaced radial cells and
$n_{\rm\theta}=768$ equally spaced cells in the azimuthal $[0,2\pi]$ domain.
We include a radially variable viscosity of the standard parameter $\alpha$
(Shakura & Sunyaev, 1973) as

$\alpha(r)=\alpha_{\rm in}-\frac{\alpha_{\rm in}-\alpha_{\rm
out}}{2}\left[1+{\rm tanh}\left(\frac{r-\xi}{\sigma r_{0}}\right)\right],$ (2)

where the inner and outer viscosities are $\alpha_{\rm in}=1\times 10^{-4}$
and $\alpha_{\rm out}=5\times 10^{-3}$, $\xi=144$ au indicate the midpoint of
the transition and $\sigma=1.25$ defines the slope. Similar to the values
provided by Liu et al. (2018).

A system of 4 planets was embedded. The location of the two outer ones is
indicated in Teague et al. (2018) through kinematic detections. The third
planet’s position (i.e., the inner planet of the outer pair) is strongly
associated with the potential velocity kink reported by (Pinte et al., 2020).
The two inner planets are tightly packed and their parameters (masses and
orbits) were derived numerically guided by the results in Garrido-Deutelmoser
et al. (2022), where it was found that the planets should be gap-sharing with
forming vortices at their Lagrange points, implying a condition in the
planetary separation333This expression has been tested for planets with masses
near the thermal mass of $M_{\rm th}=M_{\star}h^{3}\sim 0.25M_{\rm J}\simeq
80M_{\oplus}$. of:

$\displaystyle\Delta a\lesssim 4.6H\simeq 11.5\mbox{ au},$ (3)

where $H$ the scale high of the disk. In turn, the masses are constrained by
the width of the gap. After a few dozen simulations attempting to match disk
morphology in continuum observations at the $\sim 48$ au region, we choose to
place the planets at $a_{1}=46$, $a_{2}=54$, $a_{3}=84.5$, and $a_{4}=137$ au
with their respective masses of $M_{1}=85M_{\oplus}$, $M_{2}=60M_{\oplus}$,
$M_{3}=0.4M_{\rm Jup}$, and $M_{4}=1M_{\rm Jup}$. The four bodies can
gravitationally interact between them, but they do not feel the disk. We ran
an extra model to compare against previous works, substituting a Jupiter at 48
au instead of the inner package of planets. Both cases were evolved for $0.48$
Myrs, equivalent to 2000 orbits of the innermost planet.

###  2.2 Radiative Transfer

We convert the 2D dust surface density into a 3D volume density assuming the
vertical approximation given by

$\rho_{{\rm d}_{j}}(r,\phi,\theta)=\frac{\Sigma(r,\phi)}{\sqrt{2\pi}H_{{\rm
d}_{j}}}{\rm exp}\left(-\frac{z^{2}}{2H_{{\rm d}_{j}}^{2}}\right),$ (4)

where $z=r\cos(\theta)$ and the dust settling follows the diffusion model
$H_{{\rm d}_{j}}=\sqrt{D_{z}/(D_{z}+{\rm St}_{j})}H$, with $D_{z}=0.6\alpha$
the vertical diffusion coefficient, and ${\rm St}_{j}$ the Stokes number of
the species $j$ (Weber et al., 2022). We assume an intrinsic volume density
for the particles $\rho_{s}=2\;\rm gr\;cm^{-3}$ and a power law for the grain
size distribution, such that $n(a)\propto a^{-3.5}$. We assumed a dust
composition of 20$\%$ amorphous carbon, 20$\%$ water ices and 60$\%$
silicates, where the corresponding dust opacities were computed with the code
provided by (Bohren & Huffman, 1983). The polar direction is distributed in 64
equally spaced cells and extended in [$80.6^{\circ},99.4^{\circ}$] inclination
domain. We use $n_{\rm phot}=10^{8}$ photon packages to calculate the dust
temperature, and $n_{\rm scatt}=10^{7}$ photon packages to trace the thermal
emission. We use a full anisotropic scattering with polarization treatment.
The system is assumed to be at a distance of $101$ pc with a central star of
mass $1.9$ $M_{\odot}$ and effective temperature $T_{\rm eff}=9,330$ K. The
inclination is taken to be $i=46^{\circ}$ and position angle PA$=133^{\circ}$.

###  2.3 Synthetic Observations

We use SIMIO that contains a suit of functions for CASA 5.6.2. We select the
template designed for HD 163296 to create images with the same $uv$-coverage
as observation from Isella et al. (2018). We set the rescale flux option in
0.4 to get similar intensities. In addition, we add simple thermal
noise444https://simio-
continuum.readthedocs.io/en/main/tutorials/tutorial_3.html of level $12{\rm
mJy}$ to finally get RMS noise of $0.022\;{\rm mJy\;beam^{-1}}$.

##  3 Local gap depletion

Figure 1 shows the surface density after $\sim 0.5$ Myr (2,000 orbits) for the
single-Jovian case (panel a) and the two sub-Saturns with masses
$85M_{\oplus}$ and $60M_{\oplus}$ (panel b). The corresponding azimuthally-
averaged profiles in panel (c) show that $\sim 95\%$ of gas is depleted for
the single-Jovian, while only $\sim 55\%$ is depleted for the two-planet case.
Despite of their lower masses the planets pair creates the same gap width as
the Jovian. This shallower gaps for fixed gap width are expected in compact
multi-planet systems due the planet lower masses (depth $\Sigma_{\rm
gap}/\Sigma_{0}\propto M_{p}^{-2}$) and angular flux transferred by the
neighbouring planets (Duffell & Dong, 2015; Garrido-Deutelmoser et al., 2022).

As argued by Zhang et al. (2021), if a Jupiter-mass planet had opened the
corresponding CO gap, it would be $10$ times deeper than what is actually
observed. Instead, by embedding two planets we can alleviate these differences
and largely reduce this discrepancy as shown in panel (c). Therefore, the
depletion values would be closer to the results from observations. In order to
better quantify this, we compare the $\rm CO$ column density gaps with the
surface density from both models555The $\Sigma$ profiles were processed under
smooth function methodology described described in Appendix A.. As shown in
Figure 2, this approach reproduces the gas gap reasonably well in the two-
planet case, largely matching the depths and widths with a small offset of the
peaks by $\sim 4$ au. In turn, the single-Jovian case is too deep compared to
the observations as expected.

Beyond the third planet at $\sim 85$ au, neither of the two models (single
Jovian or compact pair) is able to reproduce the gas gaps. This was already
noted in Zhang et al. (2021) when comparing with the models from Teague et al.
(2021) and it may partly be explained by the presence of a CO snowline at 65
au. Outside the mid-plane CO snowline, CO freezes out at the disk mid-plane
and therefore the CO gap properties (e.g., width and depth) may deviate from
that of gas gaps due to vertical temperature and CO abundance variations. We
recall that our work mostly focuses on the gap and crescent at $\sim 50$ au
where these issues can be more securely avoided.

Figure 1: Panels show a time evolution after 2000 orbits at 48 au ($\sim
4.8\times 10^{5}$ yrs). (a) Surface density $\Sigma$ maps in log-scale for a
single Jupiter planet at 48 au. (b) The same as (a), but for a two-planet
system with $85M_{\oplus}$ and $60M_{\oplus}$ instead the Jupiter. The crosses
denote the position of the planets and the white lines indicate their orbits.
The disk rotates in a clockwise direction. (c) Azimuthally averaged
$\Sigma/\Sigma_{0}$ profiles for both models.

##  4 The Crescent Feature

A closely-packed planet pair can directly affect the gas around each other’s
coorbital regions by depositing angular momentum from wave steepening and
subsequent shocks. These shocks may strengthen the vortencity around the
stable $L_{4}$ and $L_{5}$ Lagrange points, thus enabling the effective gas
and dust trapping for at least thousands of orbital periods (Garrido-
Deutelmoser et al., 2022).

The overdensity around either $L_{4}$ or $L_{5}$ for the inner and outer
planet is a highly dynamic problem, with the most prominent structure
chaotically alternating location. This said, we observe that for several
combinations of parameters (surface density, aspect ratio, planetary masses,
and so on), the outer planet often retains large amounts of material around
$L_{5}$. The depicted behavior would leave an off-center substructure inside
the gap that greatly reproduces the distinctive emission in the disk as shown
in Figure 3. More quantitatively, our models matches the observed azimuthal
extent of $\sim 45^{\circ}$ and the radial intensity profiles passing through
the crescent (peaks and troughs inside the ring at $\sim 68$ au; see Appendix
B for more details).

Figure 2: Comparison of column density gap in $\rm C^{18}O\;(2-1)$ line
observation found by Zhang et al. (2021) with gas gaps derived from surface
density profiles in this work.

###  4.1 Dynamical behavior

Figure 4 compares the gas and dust density distribution for different dust
grain sizes. All the dust species show a clear shared gap between the two
inner planets. In the largest size, the co-orbital regions of each one display
an overdensity at $L_{5}$, but the more prominent substructure belongs to the
second planet. This output has taken into account two conditions. The choice
of planetary masses in the inner system must be lower for the body orbiting
the substructure, and the presence of the two outer planets. If we neglect
either, the Lagrange points can still trap dust, but the mass distribution may
change to make some $L_{4}$ or the inner $L_{5}$ the most prominent. As shown
in Garrido-Deutelmoser et al. (2022), this evolution is highly dynamic and
irregular so that the overdensities around $L_{4}$ and $L_{5}$ constantly
change. However, the final morphology in our configuration comes from the
early stages of evolution.

Theoretically, the stable Lagrange point $L_{5}$ is located at $60^{\circ}$
from the planet at its trailing position. However, our model shows that the
center of the crescent is slightly shifted and can vary between $65^{\circ}$
and $85^{\circ}$ for reason we still do not understand and deserve further
investigation. Accordingly, the position of the proposed planet will be also
slightly shifted from the Lagrange point. Finally, we observe that azimuthal
extent of the crescent remains roughly constant and equal to $\sim 45^{\circ}$
similar to the observations.

###  4.2 Behavior for different dust species

In our fiducial model, the vortensity for $L_{5}$ of the outer planet is
stronger than that of the inner one (not shown). Therefore, the dust
accumulation is generally expected to be greater around the orbit of the outer
planet for all sizes. This is especially true for small sizes that are well
coupled to the gas and librate with larger amplitudes around $L_{5}$ (panels b
and c in Fig. 4). As the size of the grains increases (Stokes numbers approach
unity, panels d to f), the dust distribution becomes compact toward the center
of the Lagrange point (Montesinos et al., 2020) and we can even see some
tenuous accumulation around the inner planet’s $L_{5}$ point at 1.9 cm (panel
f).

Figure 3: Band 6 ($\lambda\sim$1.25 mm) comparison between dust continuum
image from ALMA observation (Isella et al., 2018) and our synthetic
observation after $\approx 4.8\times 10^{5}$ yrs. The synthesized beam are the
same for both images (0.038′′ $\times$ 0.048′′, 82.5∘), represented by white
ellipse at the bottom left corner for each image. The synthetic image is
projected with an inclination $i=46^{\circ}$ and a position angle
PA$=133^{\circ}$. The enumerated white dots indicate the position of potential
planets associated with: (1) and (4) resonance angles in Laplace chains (see
§5.1), (2) Lagrange point $L_{5}$ from the simulation, (3) velocity kink
reported by Pinte et al. (2020). The white dashed lines denote the orbits of
planets in the simulation. The right bottom Cartesian coordinate describe the
prescription to estimate the azimuthal angles. The disk rotates in a clockwise
direction.

##  5 A Laplace resonance chain

Our fiducial simulation has 4 planets with period ratios
$P_{2}/P_{1}=(54/46)^{3/2}=1.27$, $(84.5/54)^{3/2}=1.96$, and
$(137/84.5)^{3/2}=2.06$. Therefore, the three outer planets lie near a $1:2:4$
commensurability, which becomes nearly exact (within $1\%$) if planet 3
changes from 84.5 au to 86 au. This fact begs the question of whether disk-
driven migration may have placed the planets in their current, near resonant,
orbits666We note that a disk-driven migration may also lead to offsets from
the exact commensurabilities by either wake-planet interactions (Baruteau &
Papaloizou, 2013), disk-driven precession (e.g., Tamayo et al., 2015) or
resonant repulsion (e.g., Papaloizou, 2011). .

From Figure 1, the planets carved relatively shallow gaps, so we may estimate
the rate of orbital migration following Kanagawa & Szuszkiewicz 2020 as:

$\displaystyle\tau_{a}\equiv\left|\frac{a}{\dot{a}}\right|$
$\displaystyle\simeq$ $\displaystyle\left(\frac{M_{\star}}{M_{\rm
p}}\right)\left(\frac{M_{\star}}{\Sigma_{\rm min}a_{\rm
p}^{2}}\right)\frac{h_{\rm p}^{2}}{\Omega_{\rm K,p}},$ (5)
$\displaystyle\simeq$ $\displaystyle 0.4\;\mbox{Myr}\left(\frac{10\;{\rm
gr\;cm^{-2}}}{\Sigma_{\rm min}}\right)\left(\frac{100\;{\rm au}}{a_{\rm
p}}\right)^{1/2}$
$\displaystyle\left(\frac{2M_{\odot}}{M_{\star}}\right)\left(\frac{1M_{\rm
J}}{M_{\rm p}}\right)\left(\frac{h_{\rm p}}{0.1}\right)^{2}.$

where $\Sigma_{\rm min}$ corresponds to the local density at the base of the
density gap. Using the fiducial planetary parameters and
$M_{\star}=1.9\,M_{\odot}$, a fix aspect ratio $h=0.1$ and the surface density
constraints from Zhang et al. (2021, Table 5 therein), we observe that the
migration timescales are all comparable to the age of the system making
migration a plausible scenario.

As a proof of concept, in Figure 5 we show an N-body integration using REBOUND
(Rein & Liu, 2012) and prescribing the damping timescales $\tau_{a}$ and
$\tau_{e}\equiv e/|\dot{e}|$ for each planet in order to mimic planet-disk
interactions in the REBOUNDx library (Tamayo et al., 2019). We set
$\tau_{a}/\tau_{e}=100$ with $\tau_{a}$ computed using Eq. 5, see values for
the orbital decay timescales in Table 1. We begin the simulations with the
planets further away from their current positions and let them migrate due to
their interaction with the gaseous component of the disk, where kinematic
evidence has been detected in the gas at $\sim 260$ au (Pinte et al., 2020;
Teague et al., 2021), and extensions in CO (2-1) up to $\sim 500$ au (Zhang et
al., 2021). The evolution shows that all planet pairs are captured into a long
resonant chain after $\sim 0.5\times 10^{6}$ yrs. These correspond to a two-
body 4:3 mean-motion resonance (MMR) for the innermost planets, and a double
2:1 - 2:1 MMR for the two outer pairs, finally leading to the libration of the
following three-body angles:

$\displaystyle\phi_{123}$ $\displaystyle=$ $\displaystyle
3\lambda_{1}-5\lambda_{2}+2\lambda_{3},\mbox{ and}$ (6)
$\displaystyle\phi_{234}$ $\displaystyle=$ $\displaystyle
2\lambda_{2}-6\lambda_{3}+4\lambda_{4}.$ (7)

Both $\phi_{123}$ and $\phi_{234}$ have small-amplitude libration ($\sim
2^{\circ}$) and their libration centers are $187^{\circ}$ and $218^{\circ}$,
respectively. Note that when the four planets reach the reported semi-major
axis at approximately the same time $\sim 10^{6}$ yrs (gray vertical line in
panel b), they are already captured in the two- and three-planet resonances,
showing that the proposed configuration with our hydrodynamical simulations
presented in the previous sections is possible.

Figure 4: Face-on gas and dust surface density $\Sigma$ from the hydrodynamic
model after $\approx 4.8\times 10^{5}$ yrs (2000 orbits at 48 au). The panels
correspond to different fluids. The crosses denote the position of the planets
and the white dashed lines indicate their orbits. The disk rotates in a
clockwise direction.

###  5.1 Predicting the position angle (PA) of the planet candidates using
3-body resonances

Because the orbits of planets $i$ are coplanar and nearly circular
($e_{i}\lesssim 0.07$), the mean longitudes $\lambda_{i}$ are close to the
true longitudes $\varpi_{i}+f_{i}$ and will likely librate around the same
angles. Defining an arbitrary reference frame, rotated by ${\rm PA}_{0}$ we
write the position angles (PAs) as ${\rm PA}_{i}=\varpi_{i}+f_{i}+{\rm
PA}_{0}$, and define the following three-body angle, frame-independent777Since
the critical angles satisfy the D’Alembert property, these combinations are
independent of the reference frame., combinations:

$\displaystyle{\rm PA}_{123}$ $\displaystyle=$ $\displaystyle 3{\rm
PA}_{1}-5{\rm PA}_{2}+2{\rm PA}_{3},~{}~{}\mbox{ and}$ (8) $\displaystyle{\rm
PA}_{234}$ $\displaystyle=$ $\displaystyle 2{\rm PA}_{2}-6{\rm PA}_{3}+4{\rm
PA}_{4}.$ (9)

Figure 5: Potential migratory history of the four-planet system locking the
planets in a long orbital resonance chain leaving the outer three planets near
a consecutive 2:1 commensurability and the innermost pair near 4:3 (panels a
and b). The eccentricities remain small after the capture (panel c) and the
three-body resonant angles $\phi_{123}=3\lambda_{1}-5\lambda_{2}+2\lambda_{3}$
and $\phi_{234}=2\lambda_{2}-6\lambda_{3}+4\lambda_{4}$ undergo small-
amplitude libration (panel d). The bottom panel exhibits the corresponding
combinations of the position angles in the system: ${\rm PA}_{123}=3{\rm
PA}_{1}-5{\rm PA}_{2}+2{\rm PA}_{3}$ and ${\rm PA}_{234}=2{\rm PA}_{2}-6{\rm
PA}_{3}+4{\rm PA}_{4}$.

From panel (e) in Figure 5, we observe that these angles librate around ${\rm
PA}_{123}\sim 185^{\circ}$ and ${\rm PA}_{234}\sim 210^{\circ}$ similar to
$\phi_{123}$ and $\phi_{234}$, but with larger amplitudes (near $100^{\circ}$
in both cases). This is expected due to the non-zero eccentricities.

Assuming that we know two angles, say ${\rm PA}_{2}$ due to the crescent and
${\rm PA}_{3}$ due to a velocity kink, we can use the above relations to
constrain ${\rm PA}_{1}$ and ${\rm PA}_{4}$ as:

$\displaystyle{\rm PA}_{1}$ $\displaystyle\sim$
$\displaystyle{\textstyle\frac{1}{3}{\rm PA_{123}}}+{\textstyle\frac{5}{3}{\rm
PA_{2}}}-{\textstyle\frac{2}{3}{\rm PA_{3}}},~{}~{}\mbox{ and}$ (10)
$\displaystyle{\rm PA}_{4}$ $\displaystyle\sim$
$\displaystyle{\textstyle\frac{1}{4}{\rm PA_{234}}}-{\textstyle\frac{1}{2}{\rm
PA_{2}}}+{\textstyle\frac{3}{2}{\rm PA_{3}}}.$ (11)

Our model indicates that planet 2 (54 au) is $\sim 75^{\circ}\pm 10^{\circ}$
ahead of the crescent center888Due to the clockwise rotation of the disk, our
angle convention PA is given the coordinate axes at the bottom of Figure 3.,
which corresponds to $\rm PA_{2}\simeq 32^{\circ}$. In addition, as reported
by Pinte et al. (2020), the planet 3 at 86 au has $\rm PA_{3}\simeq
357^{\circ}$ associated with a velocity kink. Thus, our model predicts that
the planets at $46$ au and 137 au should have position angles of $\rm
PA_{1}\sim 237^{\circ}$ and $\rm PA_{4}\sim 212^{\circ}$ respectively (see the
left panel in Figure 3).

Quite recently999After the submission of our manuscript to the journal.,
Alarcón et al. (2022) localized strong kinematic deviation in C I line
emission. The position of this structure lies inside the gap at 48 au, which
azimuthally coincides with our predicted planet 1 at $\rm PA_{1}\sim
237^{\circ}$. We note that this predicted planet differs from the one proposed
by Isella et al. (2018) and incorrectly quoted by Alarcón et al. (2022) as
coincident with the outflow. The reason is that the disk rotates in a
clockwise direction so the proposed planet invoked to explain the crescent as
a $L_{5}$ feature, similar to Rodenkirch et al. (2021), will actually show
ahead of the crescent. In this way, the C I deviation cannot be explained by a
co-rotational planet that is also responsible for the crescent, unless the
dust accumulation around $L_{4}$ becomes more prominent than that of $L_{5}$,
which is unlikely (Rodenkirch et al., 2021; Garrido-Deutelmoser et al., 2022).

Table 1: Migration rate estimates

— | — | — | —
---|---|---|---
$a_{\textrm{p}}$ | $M_{\textrm{p}}$ | $\Sigma_{\textrm{min}}\textrm{ [gr/cm/cm]}$ |  $\tau_{\textrm{a}}$ [Myr]
46 au | 85 $M_{\oplus}$ | 12 | 1.8
54 au | 60 $M_{\oplus}$ | 19 | 1.5
84.5 au | 127 $M_{\oplus}$ | 9.3 | 1.1
137 au | 317 $M_{\oplus}$ | 4.2 | 0.8

Note. — The values of $\Sigma_{\rm min}$ are taken from Table 5 in Zhang et
al. (2021), except for the innermost one provided by our model.

###  5.2 Resonances in other systems

We remark that resonances may be a common outcome in these young systems,
including the embedded planets in PDS 70 (Bae et al., 2019), as well as young,
but disk-free systems, like HR 8799 also in a long resonance chain involving
four planets (Goździewski & Migaszewski, 2020).

Similar to our work, a compact multi-planet system has been proposed using the
axisymmetric dust gaps and rings of HL Tau (ALMA Partnership et al., 2015),
where a resonant configuration may promote the system’s dynamical stability
(Tamayo et al., 2015). In our case, we use not only the system’s migration
history and dust rings and gaps, but also add the constraints from the
crescent shape structure and the CO gas emission.

##  6 Conclusions

We have provided a global model for HD 163296 with four planets (semi-major
axes in the range of $40-140$ au) that can reproduce the rings and gaps in the
dust continuum and the shallow gaps in the gas constrained by the CO emission.
A key ingredient in our model is the presence of two sub-Saturn-mass planets
near the 4:3 resonance opening the gap at $\sim 48$ au, where the crescent
corresponds to the $L_{5}$ Lagrange point of the outer planet at 54 au.

We show that the four-planet system may be part of a long resonance chain with
the inner two in a 4:3 MMR and the outer three in a 1:2:4 Laplace resonance
chain, consistent with a history of convergent migration within the disk. Our
proposed three-body resonances allow to relate the planetary radial and
angular positions, and based on the crescent location at 55 au and the
proposed location by Pinte et al. (2020) for the planet at $\simeq 86$ au, our
model predicts two planets: i) a sub-Saturn at 46 au and $\rm PA\sim
237^{\circ}$; ii) a Jovian at 137 au and $\rm PA\sim 212^{\circ}$(Figure 3).

Overall, our work shows that tightly-spaced planetary systems, often found at
small orbital distances in transiting surveys, may leave detectable imprints
in protoplanetary disks at much larger separations.

_Acknowledgements_ The authors would like to thank Andrew Youdin, Kaitlin
Kratter, Diego Muñoz, Matt Russo, Pablo Benítez-Llambay, Simón Cassasus, and
Ximena S. Ramos for helpful discussions that improved the quality of this work
and Juan Veliz for his support with the cluster logistics. Finally we thank
the anonymous reviewer for the thorough and useful report. J.G. acknowledge
support by ANID, – Millennium Science Initiative Program – NCN19_171 and
FONDECYT Regular grant 1210425. The Geryon cluster at the Centro de Astro-
Ingenieria UC was extensively used for the calculations performed in this
paper. BASAL CATA PFB-06, the Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL
130008 provided funding for several improvements to the Geryon cluster. C.P.
acknowledges support from ANID Millennium Science Initiative-ICN12_009, CATA-
Basal AFB-170002, ANID BASAL project FB210003, FONDECYT Regular grant 1210425,
CASSACA grant CCJRF2105, and ANID+REC Convocatoria Nacional subvencion a la
instalacion en la Academia convocatoria 2020 PAI77200076. C.C. acknowledges
FNRS Grant No. F.4523.20 (DYNAMITE MIS-project). V.V.G. acknowledges support
from FONDECYT Regular 1221352, ANID project Basal AFB-170002, and ANID, –
Millennium Science Initiative Program – NCN19_171. K.Z. acknowledges the
support of the Office of the Vice Chancellor for Research and Graduate
Education at the University of Wisconsin – Madison with funding from the
Wisconsin Alumni Research Foundation.

## References

  * Alarcón et al. (2022) Alarcón, F., Bergin, E. A., & Teague, R. 2022, ApJ, 941, L24, doi: 10.3847/2041-8213/aca6e6
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3, doi: 10.1088/2041-8205/808/1/L3
  * Andrews (2020) Andrews, S. M. 2020, ARA&A, 58, 483, doi: 10.1146/annurev-astro-031220-010302
  * Bae et al. (2022) Bae, J., Isella, A., Zhu, Z., et al. 2022, arXiv e-prints, arXiv:2210.13314.  https://arxiv.org/abs/2210.13314
  * Bae et al. (2019) Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJ, 884, L41, doi: 10.3847/2041-8213/ab46b0
  * Baruteau & Papaloizou (2013) Baruteau, C., & Papaloizou, J. C. B. 2013, ApJ, 778, 7, doi: 10.1088/0004-637X/778/1/7
  * Benisty et al. (2021) Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJ, 916, L2, doi: 10.3847/2041-8213/ac0f83
  * Benítez-Llambay et al. (2019) Benítez-Llambay, P., Krapp, L., & Pessah, M. E. 2019, ApJS, 241, 25, doi: 10.3847/1538-4365/ab0a0e
  * Bohren & Huffman (1983) Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles
  * Duffell (2020) Duffell, P. C. 2020, ApJ, 889, 16, doi: 10.3847/1538-4357/ab5b0f
  * Duffell & Dong (2015) Duffell, P. C., & Dong, R. 2015, ApJ, 802, 42, doi: 10.1088/0004-637X/802/1/42
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool.  http://ascl.net/1202.015
  * Garrido-Deutelmoser et al. (2022) Garrido-Deutelmoser, J., Petrovich, C., Krapp, L., Kratter, K. M., & Dong, R. 2022, arXiv e-prints, arXiv:2204.09074.  https://arxiv.org/abs/2204.09074
  * Goździewski & Migaszewski (2020) Goździewski, K., & Migaszewski, C. 2020, ApJ, 902, L40, doi: 10.3847/2041-8213/abb881
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Isella et al. (2018) Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49, doi: 10.3847/2041-8213/aaf747
  * Kanagawa & Szuszkiewicz (2020) Kanagawa, K. D., & Szuszkiewicz, E. 2020, ApJ, 894, 59, doi: 10.3847/1538-4357/ab862f
  * Kanagawa et al. (2015) Kanagawa, K. D., Tanaka, H., Muto, T., Tanigawa, T., & Takeuchi, T. 2015, MNRAS, 448, 994, doi: 10.1093/mnras/stv025
  * Law et al. (2021) Law, C. J., Loomis, R. A., Teague, R., et al. 2021, ApJS, in press.  https://arxiv.org/abs/2109.06210
  * Liu et al. (2018) Liu, S.-F., Jin, S., Li, S., Isella, A., & Li, H. 2018, ApJ, 857, 87, doi: 10.3847/1538-4357/aab718
  * Long et al. (2022) Long, F., Andrews, S. M., Zhang, S., et al. 2022, ApJ, 937, L1, doi: 10.3847/2041-8213/ac8b10
  * Masset (2000) Masset, F. 2000, A&AS, 141, 165, doi: 10.1051/aas:2000116
  * Montesinos et al. (2020) Montesinos, M., Garrido-Deutelmoser, J., Olofsson, J., et al. 2020, A&A, 642, A224, doi: 10.1051/0004-6361/202038758
  * Papaloizou (2011) Papaloizou, J. C. B. 2011, Celestial Mechanics and Dynamical Astronomy, 111, 83, doi: 10.1007/s10569-011-9344-4
  * Perez et al. (2015) Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJ, 811, L5, doi: 10.1088/2041-8205/811/1/L5
  * Pérez et al. (2020) Pérez, S., Casassus, S., Hales, A., et al. 2020, ApJ, 889, L24, doi: 10.3847/2041-8213/ab6b2b
  * Pinte et al. (2019) Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nature Astronomy, 3, 1109, doi: 10.1038/s41550-019-0852-6
  * Pinte et al. (2020) Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9, doi: 10.3847/2041-8213/ab6dda
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rodenkirch et al. (2021) Rodenkirch, P. J., Rometsch, T., Dullemond, C. P., Weber, P., & Kley, W. 2021, A&A, 647, A174, doi: 10.1051/0004-6361/202038484
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Szulágyi et al. (2018) Szulágyi, J., Plas, G. v. d., Meyer, M. R., et al. 2018, MNRAS, 473, 3573, doi: 10.1093/mnras/stx2602
  * Tamayo et al. (2019) Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. 2019, Monthly Notices of the Royal Astronomical Society, 491, 2885–2901, doi: 10.1093/mnras/stz2870
  * Tamayo et al. (2015) Tamayo, D., Triaud, A. H. M. J., Menou, K., & Rein, H. 2015, ApJ, 805, 100, doi: 10.1088/0004-637X/805/2/100
  * Teague et al. (2018) Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12, doi: 10.3847/2041-8213/aac6d7
  * Teague et al. (2021) Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18, doi: 10.3847/1538-4365/ac1438
  * van der Marel et al. (2019) van der Marel, N., Dong, R., di Francesco, J., Williams, J. P., & Tobin, J. 2019, ApJ, 872, 112, doi: 10.3847/1538-4357/aafd31
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Weber et al. (2022) Weber, P., Casassus, S., & Pérez, S. 2022, MNRAS, 510, 1612, doi: 10.1093/mnras/stab3438
  * Weber et al. (2019) Weber, P., Pérez, S., Benítez-Llambay, P., et al. 2019, ApJ, 884, 178, doi: 10.3847/1538-4357/ab412f
  * Zhang et al. (2021) Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5, doi: 10.3847/1538-4365/ac1580

##  Appendix A GAS GAP CALCULATION

Figure 6: Black line denote the convolved surface density profile of models
and grey dashed line the respective smooth function. Panel (a) and (b)
represent the single-planet case and two-planet case respectively. The crosses
indicate the position of the planets.

In Zhang et al. (2021) a smooth function was subtracted from $\rm N_{CO}$
column density profiles to better characterize substructures in the residual
values. To compare these results with our models, we follow the same
procedure. First, the $\Sigma$ maps from hydrodynamic simulations were
convolved with a circular Gaussian beam of $0.15^{\prime\prime}$, which has
the same size as MAPS CO (2-1) line observations. Then, their azimuthally
averaged profiles were interpolated every 2 au. In addition, the radial region
$\\{r\;{\rm[au]}:0<r_{0}<35,59<r_{1}<72,98<r_{2}<110,r_{3}>170\\}$ was
selected to describe the gaps. Both were taken as input for the
smoothfit101010https://pypi.org/project/smoothfit/ module. The Figure 6 shows
the outputs of smoothed profile represented by the grey dashed line and the
convolved surface density profile in black lines. The cyan dots denote the
regions in which the function acts. The Figure 2 shows the residual between
lines to provide a reasonable comparison with CO gaps observations.

##  Appendix B Radial intensity

We quantify the intensity around the substructure region of our synthetic
model with the ALMA observation. First, we deproject the images obtaining a
face-on view to convert them to polar coordinates and then generate a radial
profile by taking the azimuthal average between PA of $300^{\circ}$ and
$350^{\circ}$. This extension fully covers the emission from the crescent. The
results are shown in the Figure 7, which is accompanied by a diagram showing
the angular slice.

Figure 7 show that radial intensity through the crescent region reaches
amplitudes higher than those observed by a factor of 1.2 at 55 au. The
emission from the substructure is clearly off-centered on the gap and resolved
in spatial resolution, showing a gap in intensity between it and the ring. The
first ring reproduces the intensities in a good way, while the second is
noticeably 1.7 times fainter.

Figure 7: Azimuthally averaged intensity profiles for synthetic and ALMA
observations after $\approx 4.8\times 10^{5}$ yrs. The crosses denote the
semi-major axes of the planets. The insert shows the ALMA observation in Band
6 with contours that reproduce the crescent and rings as well as the angular
slice used for the azimuthal average denoted by $\hat{\phi}$.

##  Appendix C Effect from disk gravity acting on planets

We briefly test whether turning on the full disk-planet interaction may lead
to morphological changes in the structure of the crescent. We recall that in
our fiducial simulation (see §2.1), while the disk do feel the planets’
gravity, the planets do not feel the disk.

We perform two-planet simulations considering only the inner planet pair near
the 4:3 commensurability (46 au and 55 au) for up to 2000 orbits of the inner
planet. The initial density has been reduced by a factor of 100 to avoid
significant migration. In figure 8 we show the density distribution for two
dust fluids of 0.02 cm and 1.9 cm grain sizes in two cases: the full disk-
planet interaction is considered (left panels, displaying a slight inward
migration at the $\sim 10\%$ level), and the disk gravity acting on the
planets is ignored (right panels, with no migration). Despite of the slight
orbital migration, we do not observe any significant changes regarding the
amount and distribution of captured material at the $L_{4}$ and $L_{5}$
Lagrange points.

Figure 8: Dust surface density $\Sigma$ for 0.02 cm and 1.9 cm grain sizes.
The crosses denote the position of the planets. The upper and bottom white
rectangles, indicate the Lagrange points $L_{4}$ and $L_{5}$ respect to the
outer planet.


# Measuring Exploration in Reinforcement Learning
via Optimal Transport in Policy Space

Reabetswe M. Nkhumise Computer Science Dept.
The University of Sheffield
Sheffield, South Yorkshire, UK  Debabrota Basu  Eq́uipe Scool
Univ. Lille, Inria, CNRS, Centrale Lille, UMR 9189- CRIStAL, France  Tony J.
Prescott  Computer Science Dept.
The University of Sheffield
Sheffield, South Yorkshire, UK  Aditya Gilra  Centrum Wiskunde & Informatica
Amsterdam, Netherlands
Computer Science Dept.
The University of Sheffield
Sheffield, South Yorkshire, UK

###### Abstract

Exploration is the key ingredient of reinforcement learning (RL) that
determines the speed and success of learning. Here, we quantify and compare
the amount of exploration and learning accomplished by a Reinforcement
Learning (RL) algorithm. Specifically, we propose a novel measure, named
Exploration Index, that quantifies the relative effort of knowledge transfer
(transferability) by an RL algorithm in comparison to supervised learning (SL)
that transforms the initial data distribution of RL to the corresponding final
data distribution. The comparison is established by formulating learning in RL
as a sequence of SL tasks, and using optimal transport based metrics to
compare the total path traversed by the RL and SL algorithms in the data
distribution space. We perform extensive empirical analysis on various
environments and with multiple algorithms to demonstrate that the exploration
index yields insights about the exploration behaviour of any RL algorithm, and
also allows us to compare the exploratory behaviours of different RL
algorithms.

##  1 Introduction

In Reinforcement Learning (RL), exploration impacts the overall performance of
the agent by determining the quality of experiences from which the agent
learns [Zhang et al., 2019]. Exploration is the activity of searching to
maximise the reward signal, and to overcome local minima (or maxima), or areas
of flat rewards [Ladosz et al., 2022]. Quantifying exploration is important as
it would allow evaluation of exploration qualities of different RL algorithms,
and exploration strategies used in the same algorithm. However, a standardised
metric to compare exploratory efforts that is algorithm-agnostic is currently
absent [Hollenstein et al., 2020, Ladosz et al., 2022].

In recent years, significant advancements have been achieved in developing
exploration techniques that improve learning [Bellemare et al., 2016, Osband
et al., 2016, Burda et al., 2019, Eysenbach et al., 2019]. Yet, there remains
a lack of consensus over approaches that can quantitatively compare these
techniques across RL algorithms [Amin et al., 2021, Ladosz et al., 2022]. This
is attributed to some methods being algorithm-specific [Tang et al., 2017],
while others provide theoretical guarantees for tabular settings which are not
extensible to non-tabular settings [Zhang et al., 2019].

The goal of exploration is to maximise the agent’s knowledge to facilitate
identification of (near-)optimal policies, i.e. discovering highly rewarding
sequences of state-action pairs. However, the environment dictates the
difficulty of visiting all the states to identify the rewarding state-actions,
which is also known as the visitation complexity [Conserva and Rauber, 2022].
An exploration strategy has to overcome this complexity to gather meaningful
samples of state-actions. We quantify the relative effort of RL algorithms to
overcome the visitation complexity and name it Exploration Index.

We consider that the crux of exploration effort lies in transferability, a
property indicating the ease or difficulty of transferring knowledge learned
across tasks [Nguyen et al., 2020, Tan et al., 2021]. During training in RL,
each policy model update encountered by the agent results in a policy, and
each policy generates an associated dataset of state-action pairs (aka
trajectories). Therefore, a classification model learning sequentially on
these datasets can be equated to the learning process of the RL agent
(illustrated in Figure 1). Since the exploration strategy is responsible for
policies encountered during training, we can use the relative effort of
knowledge transfer induced by exploration to capture the exploratory behaviour
of RL algorithm.

In Supervised Learning (SL), Alvarez-Melis and Fusi [2020] showed that
transferability can be predicted based on the similarity (alternatively,
distance) between task datasets. Here, we conceptualise learning in RL as a
sequence of SL tasks and further leverage the optimal transport dataset
distance (OTDD) developed by Alvarez-Melis and Fusi [2020]. This allows us to
compute the relative exploration effort as the ratio of the overall
transferability in RL over the transferability in SL.

Figure 1: Illustration of RL as a sequence of SL tasks (seq. SL). Policies
$\pi_{n}$ encountered during learning in RL have corresponding datasets
$D_{n}$ of state-action pairs. Training a classifier successively on these
datasets, while transferring knowledge (knwl), is analogous to training the RL
agent.

Contributions. We introduce an Exploration Index to measure the relative
effort of knowledge transfer (transferability) involved while learning in RL.
We compute it by comparing similarities of datasets of state-action pairs of
subsequent policies yielded during training using the optimal transport
dataset distance (OTDD). Further, we perform an extensive analysis of
Exploration Index across RL environments and algorithms to show its universal
applicability to compare different exploration strategies used in the same RL
algorithm, exploratory efforts of different RL algorithms, and hardness of
different tasks.

##  2 Preliminaries

We introduce the basic setup of RL and Markov decision processes assumed in
this paper. Then, we discuss relevant aspects of transfer learning,
Wasserstein distance, and OTDD. Finally, we propose an extension of OTDD to
the space of probability distributions.

RL setup. We consider an RL setup consisting of an agent interacting with an
environment in discrete timesteps. At each timestep ${t}\in\mathbb{N}$, the
agent observes a state ${s_{t}}$, executes an action ${a_{t}}$, and receives a
scalar reward ${r_{t}}$. The behaviour of the agent is defined by a policy
${\pi}$, which maps the observed states to actions. The environment is
modelled as a Markov Decision Process (MDP) with a state space $\mathcal{S}$,
action space $\mathcal{A}$, transition dynamics
$\mathcal{T}:\mathcal{S}\times\mathcal{A}\to\mathcal{S}$, and reward function
${r:\mathcal{S}\times\mathcal{A}\to\mathbb{R}}$. During task execution, the
agent issues actions in response to states visited. Thus, we obtain a sequence
of states and actions ${\tau=(s_{0},a_{0},s_{1},a_{1},...,s_{T},a_{T})}$,
called a trajectory. The return, over a finite a time horizon $T$, is the
discounted cumulative reward over the trajectory
${R(\tau)=\sum_{t=0}^{T}\gamma^{t}r(s_{t},a_{t})}$, with a discounting factor
${\gamma\in[0,1]}$. Note that the trajectory and return depend on the actions
selected, and therefore on the policy ${\pi}$. The goal in RL is to learn a
policy that maximises the expected return
${J=\mathbb{E}_{\tau\sim\pi}[R(\tau)]}$, i.e. learn a policy that generates a
trajectory class $C(\pi)=\\{\tau_{i}\\}_{i=0}^{n\to\infty}$ corresponding to
the highest possible expected return.

#### Transfer Learning and Task Similarity in Supervised Learning.

Transfer learning studies how the performance on a target task can be improved
by knowledge learned from a source task [Tan et al., 2021]. It is known that
transferring knowledge across similar tasks should be easier than across
distant ones [Alvarez-Melis and Fusi, 2020]. Hence, transferability between
tasks is often quantified by their similarity. High transferability implies
low efforts of knowledge transfer as the model pre-trained on source data can
generalise well on target task [Gao and Chaudhari, 2021].

Task similarity in SL has been studied extensively to capture transferability.
For instance, various methods evaluate task similarity according to model
performance [Zamir et al., 2018], embedded vector representations [Achille et
al., 2019], discrepancy measures between domains [Ben-David et al., 2006],
information theory [Tran et al., 2019], and optimal transport theory [Alvarez-
Melis and Fusi, 2020, Gao and Chaudhari, 2021, Tan et al., 2021]. We employ a
method called optimal transport dataset distance (OTDD) [Alvarez-Melis and
Fusi, 2020] to determine the effort of knowledge transfer in SL tasks because
it is simple and flexible to implement in practice.

#### Optimal Transport: Wasserstein Distance.

Optimal transport studies how to optimally shift a probability measure
$u\in\mathcal{P}(\mathcal{Z})$ on a metric space $\mathcal{Z}$ to a
probability measure $v\in\mathcal{P}(\mathcal{Z})$ on the same metric space.
The optimal shift is quantified by Wasserstein distance [Peyré and Cuturi,
2019]. Specifically, the 2-Wasserstein distance of $u$ and $v$ is

$W_{2}^{2}(u,v)\triangleq\min_{\pi\in\Pi(u,v)}\int_{\mathcal{Z}\times\mathcal{Z}}c(z,z^{\prime})d\pi(z,z^{\prime}).$
(1)

Here, $\Pi(u,v)$ is the set of all transport plans between $u$ and $v$.
$c(z,z\textquoteright)\in\mathbb{R}^{+}$ is the cost function quantifying the
cost of transporting unit mass from $z$ to $z^{\prime}$. Note that when
$\mathcal{Z}$ is endowed with a metric $d_{\mathcal{Z}}$, the cost function is
$c(z,z\textquoteright)=d_{\mathcal{Z}}^{2}(z,z^{\prime})$. We use Wasserstein
distance to compute OTDD.

#### Optimal Transport Dataset Distance (OTDD).

Let dataset $\mathcal{D}$ be a set of feature-label pairs
$(x,y)\in\mathcal{X,Y}$ over a feature space $\mathcal{X}$ and label set
$\mathcal{Y}$. Suppose there are two datasets
$\mathcal{D_{A}}=\\{(x_{A}^{i},y_{A}^{i})\\}_{i=1}^{m}$ and
$\mathcal{D_{B}}=\\{(x_{B}^{i},y_{B}^{i})\\}_{i=1}^{n}$ with
$x_{A},x_{B}\in\mathcal{X}$ and $y_{A},y_{B}\in\mathcal{Y}$. The samples in
these datasets are thought to be drawn from joint distributions $P_{A}(x,y)$
and $P_{B}(x,y)$. The goal of SL task is to estimate the conditional
distribution model $P(y|x)$ for each dataset.

OTDD [Alvarez-Melis and Fusi, 2020] is a statistical divergence between joint
distributions $P_{A}(x,y)$ and $P_{B}(x,y)$ that captures transferability from
dataset $A$ to $B$. Here, transferability means the effort of learning the
conditional distribution $P(y|x)$ for $\mathcal{D_{B}}$ using a model pre-
trained on $\mathcal{D_{A}}$.

Following the OTDD framework, we consider the feature-label space as a metric
space, i.e. $\mathcal{Z}=\mathcal{X}\times\mathcal{Y}$, and define the cost
function in the 2-Wasserstein distance as

$d_{\mathcal{Z}}((x,y),(x^{\prime},y^{\prime}))^{2}=d_{\mathcal{X}}(x,x^{\prime})^{2}+W^{2}_{2}(\alpha_{y},\alpha_{y^{\prime}}),$
(2)

where labels are encoded as distributions over the feature space by the map
$y\to\alpha_{y}(X)\doteq P(X\mid Y=y)$.

OTDD between $P_{A}(x,y)$ and $P_{B}(x,y)$, i.e. $d_{OT}(P_{A},P_{B})$, is the
2-Wasserstein distance (Equation (1)) between $P_{A}$ and $P_{B}$ computed
with the cost function of Equation (2). Hereafter, we use either
$d_{OT}(P_{A},P_{B})$ or $d_{OT}(\mathcal{D_{A}},\mathcal{D_{B}})$
interchangeably to denote the OTDD between joint distributions $P_{A}(x,y)$
and $P_{B}(x,y)$.

In summary, OTDD is a hierarchical Wasserstein distance that first measures
the transportation cost between labels by reflecting them back to feature
space, and then sums this with the cost between the input features. Note that
OTDD is a metric [Alvarez-Melis and Fusi, 2020]. We use OTDD to compute
overall effort of knowledge transfer (i.e. transferability) during a learning
process.

###  2.1 Formulation: OTDD on the manifold of probability distributions

Let $P_{\theta}$ be a parameterised joint probability distribution from which
samples in the dataset $\mathcal{D_{\theta}}$ are drawn, i.e.
$\mathcal{D_{\theta}}=\\{(x_{A}^{i},y_{A}^{i})\\}_{i=i}^{m}\sim P_{\theta}$.
When the parameters $\theta$ undergo a slight modification, through some
update rule, the resulting distribution is denoted as $P_{\theta+d\theta}$ and
samples in the dataset $\mathcal{D}_{\theta+d\theta}$ are drawn from it.

We consider a Riemannian manifold of these parameterised joint probability
distributions $\mathcal{M}=\\{P_{\theta}(x,y)\mid\theta\in\mathbb{R}^{N}\\}$.
These types of manifolds are studied in information geometry [Amari, 2010],
and extended to the context of RL [Basu et al., 2020]. We endow the manifold
with the OTDD metric to compute the distance between points
$\theta,\theta+d\theta\in\mathcal{M}$ as
$ds=d_{OT}(P_{\theta},P_{\theta+d\theta}).$ Now, suppose throughout some
learning process, the parameters undergo continuous updates until a predefined
convergence criterion is satisfied. This process yields a continuously
differentiable curve on the manifold $\mathcal{M}$ [Gao and Chaudhari, 2021].
The length of the curve is approximated by summation of the finite distances
along it [Lott, 2008].

$C\triangleq\sum_{k=0}^{K-1}d_{OT}(P_{\theta_{k}},P_{\theta_{k+1}}).$ (3)

$\theta_{0}$ and $\theta_{K}$ represent the initial and final parameter values
before and after the learning process, respectively. Henceforth, we refer to
the space of joint probability distributions on $\mathcal{M}$ equipped with
the OTDD metric as the OTDD space. To define distances between the
distributions induced by policies computed during RL training, we use such an
OTDD space, interchangeably terming it ’policy space’.

##  3 Measuring Exploration in RL

In this section, we study the trajectories generated by policies as datasets,
and leverage the representation of RL as sequential SLs (Figure 1) to present
our exploration metric.

###  3.1 Episodic trajectories as datasets

We define the dataset of a policy $\pi$ as the state-action pairs originating
from the set of possible trajectories obtained by playing the policy, i.e.
$C(\pi)$. That is, $\mathcal{D}_{\pi}=\\{(s_{(i)},a_{(i)})\\}_{i=0}^{m}\in
C(\pi)\subseteq\mathcal{S}\times\mathcal{A}$. These samples in the dataset are
drawn from a joint distribution $P_{\pi}(s,a)$.

Proposition 3.1 Given an MDP, the joint distribution $P_{\pi}(s,a)$ is unique
to policy ${\pi}$. (Proof in Appendix A)

If we are given $\mathcal{D}_{\pi}$, generated by an expert policy, we can
train a policy model on it in a supervised manner, via behaviour cloning
[Hussein et al., 2017]. Thus, knowing $\mathcal{D}_{\pi}$ can reduce an RL
problem to an SL task.

###  3.2 Exploration Index

A model that learned a policy $\pi_{i}$ during pre-training can be used to
generate $\mathcal{D}_{\pi_{i}}$. Suppose the objective is to learn $\pi_{j}$
using this trained model, and we are given $\mathcal{D}_{\pi_{j}}$, i.e.
state-action pairs generated by $\pi_{j}$. The model can simply be retrained
in a supervised manner on $\mathcal{D}_{\pi_{j}}$. Samples in
$\mathcal{D}_{\pi_{i}}$ and $\mathcal{D}_{\pi_{j}}$ can be treated as support
points of empirical joint distributions $\hat{P}_{\pi_{i}}$ and
$\hat{P}_{\pi_{j}}$, respectively. Then the effort required for knowledge
transfer in the model is $d_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{j}})$.

During training in RL, the policy model updates through
$\pi_{0},\dots,\pi_{N}$. $\pi_{0}$ is usually given by randomly sampling the
model parameters, and $\pi_{N}$ is achieved after convergence. Initially,
$\pi_{0}$ is evaluated to generate experiences, which are then used to update
the model to $\pi_{1}$. This iterative process continues until reaching
$\pi_{N}$. This process is analogous to transfer learning sequentially over
datasets $\mathcal{D}_{\pi_{0}},\dots,\mathcal{D}_{\pi_{N}}$ in a supervised
manner.

Proposition 3.2 In the OTDD space, learning sequentially over datasets
$\mathcal{D}_{\pi_{0}},\dots,\mathcal{D}_{\pi_{N}}$ using supervised learning,
is a close approximation of a policy model updating through
$\pi_{0},\dots,\pi_{N}$ during training in RL. (Proof in Appendix A)

If $\mathcal{D}_{\pi_{N}}$ was known a priori, then the effort of learning
$\pi_{N}$ given the model for $\pi_{0}$ would be
$L=d_{OT}(\mathcal{D}_{\pi_{0}},\mathcal{D}_{\pi_{N}})$, i.e. the geodesic
distance on the OTDD space. However, due to exploration the policy model moves
through $\pi_{0},\dots,\pi_{N}$, thus tracing a curve in the OTDD space and
the overall effort of learning $\pi_{N}$ is given by Equation (3). Figure 2
illustrates learning $\mathcal{D}_{\pi_{N}}$ in a supervised manner yielding
$L$ versus using RL which traces a curve $C$ mainly due to exploration.

Figure 2:  Depiction of the curve $C$ traced on the OTDD space during learning
in RL. The geodesic $L$ is the shortest path between initial and final points
on the curve, i.e. ${P}_{\pi_{0}}$ and ${P}_{\pi_{N}}$, respectively.

By definition, $L\leq C$ since $L$ is the geodesic, i.e. shortest distance
between between $\pi_{0}$ and $\pi_{N}$. We assess the effort of knowledge
transfer involved in the model learning $\pi_{N}$ from initial $\pi_{0}$,
using RL, in relation to simply using SL as:

$\mathrm{EI}\triangleq\frac{C}{L}=\frac{\sum_{k=0}^{N-1}d_{OT}(\mathcal{D}_{\pi_{k}},\mathcal{D}_{\pi_{k+1}})}{d_{OT}(\mathcal{D}_{\pi_{0}},\mathcal{D}_{\pi_{N}})}$
(4)

We define $C/L$ as the Exploration Index (EI) since it signifies the relative
effort induced by exploration during learning in RL compared to learning in a
supervised manner. This can be interpreted as the efficiency of exploration.
An EI close to 1 indicates highly efficient exploration i.e. the exploration
scheme gathers high quality experiences successfully, and hence the overall
effort of knowledge transfer (transferability) is minimised.

Let $\Phi$ be the set of all possible exploration schemes for a given MDP. An
optimal scheme can be defined as follows:

Definition 3.1 An optimal exploration scheme $\phi^{*}\in\Phi$ is the one that
minimises, the Exploration Index, $\frac{C}{L}$.

This approach effectively links exploration and knowledge, recognizing that
the objective of exploration is to maximise the agent’s knowledge such that it
facilitates (near)optimal policy identification. The ability to collect high-
quality experiences enhances the learning process [Zhang et al., 2019] as
these data points are profoundly informative, thereby reducing unnecessary
knowledge transfers.

Computing Exploration Index. The model parameters $\theta\in\mathbb{R}^{N}$
are initialised using some prior distribution $\theta_{0}\sim Q(\theta)$ (see
Appendix B). The model is evaluated for $T$ number of episodes to generate
state-action pairs that are collected in $\mathcal{D}_{\theta_{0}}$. The
parameters are then updated using the RL algorithm’s update scheme. After
every update, the model is evaluated to generate the corresponding policy’s
dataset. The datasets are collected for every model update until convergence.
Then we use the datasets to compute the $C$ and $L$ (as in Algorithm 1).

Algorithm 1 Computing Exploration Index (EI)

1: Initialise model parameters $\theta\sim p(\theta)$

2: {Perform standard RL training}

3: for $k=0,1,\ldots,K$ model updates do

4: for $t=0,1,\ldots,T$ episodes do

5: Evaluate policy $\theta_{k}$ on $s_{init}$ and $s_{goal}$

6: Store $\\{(s_{(i)},a_{(i)})\\}_{i=0}^{m}$ in $\mathcal{D}_{k}$

7: end for

8: end for

9: $L=d_{OT}(\mathcal{D}_{0},\mathcal{D}_{K})$

10: $C=0$

11: for $k=0,1,\ldots,K-1$ do

12: $C\leftarrow C+d_{OT}(\mathcal{D}_{k},\mathcal{D}_{k+1})$

13: end for

14: Output: $C/L$

##  4 Experimental Analysis

Now, we conduct extensive experiments to evaluate the usefulness of EI in
capturing exploration effort and its flexibility in comparing different RL
algorithms.

Experimental Setup. Our experiments consist of five deterministic 2D-Gridworld
tasks with actions {up, right, down, left}. The start and goal states always
located at top-left and bottom-right, respectively, except for a single task.
The tasks are illustrated in Figure 3. They are selected for their simplicity,
to enable clear interpretation of EI across tasks.

Figure 3: Five grid-world tasks with the same action space, but different
rewards, state space and location of the goal state. Tasks (a) and (d) have
rewards described by the L1 norm between the current and goal states $\lVert
s_{c}-s_{g}\rVert_{1}$. The rest issue a reward of -0.04 at every state and 1
at the goal.

Our experiments aim to address the following questions:
1\. Does EI capture expected behaviours within a simple deterministic setup
employing Q-learning with $\epsilon$-greedy and decaying $\epsilon$-greedy
strategies?
2\. Does EI scale proportionally with the task difficulty for any specific RL
algorithm?
3\. How do various exploration strategies within both identical and distinct
RL algorithms compare?
4\. How do convergence criteria impact EIs of RL algorithms?
Summary of Results. In Section 4.1, we validate EI and show that, as expected,
the EI values reduce with higher $\epsilon$ for $\epsilon$-greedy exploration.
Further, EIs are lower for Q-learning with decaying $\epsilon$-greedy strategy
compared to $\epsilon$-greedy strategy. Furthermore, the EI values are higher
for settings with sparse rewards than with dense rewards. In Section 4.2, we
show that EI scales proportionally with task difficulty, and thus, is able to
reflect task difficulty on exploration. In Section 4.3, we learn that studying
the variation and magnitude of segment lengths in policy space between
consecutive policies, as learning progresses, gives us insight into the
exploration behaviour of RL algorithms. This further enables us to decide if
the exploration behaviour is suitable or unsuitable for the task. We also
observe (Section 4.4) that when the convergence criterion is not strict, the
EI provides complementary diagnostic information, which is not provided by
convergence times of RL algorithms.

###  4.1 Validation of Exploration Index: Q-learning with $\epsilon$-greedy
strategies

The tabular Q-learning algorithm [Watkins and Dayan, 1992] is innately
deterministic. Its primary source of stochasticity is the exploration
strategy, apart from the random initialization of Q-values. This attribute
makes tabular Q-learning a testbed to observe the effects of exploration
without any effect of external noise. $\epsilon$-greedy is a simple and
undirected exploration scheme [Amin et al., 2021] that picks exploratory
actions at random. Here, we validate our index employing Q-learning with
$\epsilon$-greedy strategies.

We compare $\epsilon$-greedy and decaying $\epsilon$-greedy strategies,
deployed in tasks (a) and (b), by examining both their EI values and
convergence times (shown in Figure 4). The results are for
$\epsilon=\\{0.05,0.2,0.4,0.9\\}$, representing a wide range from 0 to 1. We
aim to evaluate the capability of the exploration strategy to stably settle at
an optimal policy. The corresponding convergence criterion is satisfied when
the policy model consistently yields maximum returns over 50 consecutive
updates.

Figure 4: Plots of EI $C/L$ vs $\epsilon$ and convergence time (in steps)
$T_{c}$ vs $\epsilon$ for a 5x5 grid-world. Top panels (a) and (b) depict
plots in the dense rewards setting, while bottom panels (c) and (d) show plots
in the sparse rewards setting. The results were averaged over 40 runs per
$\epsilon$-value in each task with Q-values initialised randomly between -1
and 1. Solid lines represents the median of the data, while the darker and
lighter shaded regions show the interquartile range and limits, respectively.

Figure 4 shows that the EI values are lower for the decaying $\epsilon$-greedy
strategy compared to the $\epsilon$-greedy strategy. This is consistent with
convergence time in describing decaying $\epsilon$-greedy as a more efficient
exploration scheme which aligns with our intuition and literature [Wei, 2020].
Once the agent possesses information about the environment, it needs to
exploit it and avoid recollecting known experiences [Maroti, 2019]enabling it
to reach a stable and optimal policy quicker.

EI values, along with convergence times, are smaller in the dense rewards
setting than the sparse rewards setting. This is consistent with literature
[Ladosz et al., 2022] which highlights that employing dense rewards rather
than sparse ones in reward shaping is a technique for improving exploration.

Larger $\epsilon$’s yield lower EI values for both algorithms in both tasks.
This shows that in this simple task, higher exploration as opposed to
exploitation is good. However, exploring extensively initially and gradually
exploiting with decaying $\epsilon$ is better.

###  4.2 Exploration Index increases with task difficulty

Figure 5 illustrates the EI values for Q-learning with decaying
$\epsilon$-greedy strategy across tasks (a) to (e) (ref. Figure 3). We chose
to assess this algorithm because it is simple and yet managed to complete all
the tasks. As in Section 4.1, convergence is reached when a policy model
yields maximum returns over 50 consecutive updates.

We observe that the EI is lowest for task (a) and highest for task (e), as
anticipated. Looking closely we can see that the medians of both EI and
convergence time for task (b) are higher than for task (c). This aligns with
our intuition as task (c) is a simpler version of task (b).

The results demonstrate that EI scales proportionally with task difficulty,
much like convergence time. This is sensible as more challenging tasks demand
greater exploration compared to simpler ones.

Figure 5: Plots depict Q-learning with decaying $\epsilon$-greedy exploration
scheme across various tasks. Left: EI $C/L$ vs tasks. Right: convergence time
$T_{c}$ vs tasks. The tasks (illustrated in Figure 3) are placed in order of
ascending difficulty. The results were averaged over 40 runs per task with
$\epsilon=0.9$.

###  4.3 Exploration Index for comparing exploration strategies

We utilised tasks (a), (b), and (d) of Figure 3 to evaluate the following RL
algorithms: Tabular Q-learning with 1) decaying $\epsilon$-greedy, 2)
$\epsilon$-greedy, 3) softmax (Boltzmann), and 4) greedy strategies, and also,
5) UCRL [Auer and Ortner, 2006], 6) SAC [Haarnoja et al., 2018, Christodoulou,
2019], and 7) DQN [Mnih et al., 2013] with decaying $\epsilon$-greedy.

The results of the experiment, depicted in Figure 6, are for the best-
performing parameters for each algorithm. In Q-learning, values in the Q-table
were randomly initialised between -1 and 1. Over 80 runs per algorithm in each
task were performed. The considered convergence criterion is met when the
policy model’s returns were consistently maximal over 5 consecutive updates,
which is typically employed, but different from the 50 updates required
earlier. Even with this reduced criterion, some algorithms had difficulty
succeeding (see Table 1). In these cases, we set a cut-off training time of
500 episodes (with 60 steps per episode) and still computed the EI (see
Methods in Appendix B) to differentiate exploration in the various algorithms
even without convergence.

Figure 6: Violin plots of EI $C/L$ and convergence time $T_{c}$ vs RL
algorithms in gridworld tasks. Panels (a) and (b) are for a 5x5 dense rewards
task, (c) and (d) for 5x5 sparse rewards task, and (e) and (f) for 15x15 dense
rewards task. The first four RL algorithms are tabular Q-learning based with
different exploration schemes, while the rest are assorted RL algorithms. We
plot mean values and distributions of the data collected over 80 runs per
algorithm for each task.

Task (a) 5x5 grid and dense rewards. We see that the decaying
$\epsilon$-greedy, $\epsilon$-greedy and softmax algorithms equivalently have
the lowest EI values. The trio is followed by SAC, greedy and DQN,
respectively. UCRL has the highest EI. However, greedy has the best
convergence time and SAC the worst in comparison to the rest. This highlights
that the relative effort of knowledge transfer does not necessarily correlate
with the convergence time across algorithms. This can be well explained by
contrasting exploration behaviours of UCRL, SAC, and Q-learning greedy as in
Figure 7.

The high values of ‘segment lengths’, i.e. OTDD distances in policy space
between consecutive policy updates, at the start of learning (Figure 7 (a)),
reveal that UCRL initially explores aggressively by traversing broader within
the OTDD space, then gradually reduces coverage distance, and then settles
into a rather high segment length between successive policies. In this manner,
UCRL expends notable effort to reach optimality quicker. However, the high
segment length ($\sim 1.0$) even after settling, compared to other algorithms
in Figure 7 suggests that the exploration strategy is less inclined to
stabilise reliably at an optimal policy, i.e. it is likely to escape from the
optimal policy region, wander into other areas, before eventually returning or
encountering another optimal policy (see Figure 9 as well). Such an
exploration strategy can be called diverse since it explores a wide range of
policies without honing in on the most promising one.

Figure 7: The behaviour of the ’segment length’ (distance between consecutive
policies on the OTDD space) over time for RL algorithms UCRL, SAC, and
Q-learning greedy on the 5x5 dense rewards task. The cumulative length plots
demonstrate how the segment length accrues across time to determine the length
of the curve $C$ on the OTDD space. The vertical dotted line represents the
mean convergence time averaged over 80 runs. For segment length plots, the
solid line represents a smoothed rendition of the noisy mean curve, shown
shaded. For cumulative length plots, the solid line depicts the mean curve and
the shaded areas indicate the range of standard deviation values.

By contrast, Figure 7 (b) shows that while the greedy strategy at first
explores widely due to random Q-value initialization, it rapidly narrows down
its coverage distance as the Q-values are updated following a greedy strategy.
Once it finds an optimal policy, it ceases further exploration. This strategy
prioritises precision in finding a promising policy and overlooks the
exploration of alternative policies. The drawback is that it might miss out on
potentially superior policies in unexplored regions. This type of exploration
can be described as precise. It conserves effort of knowledge transfer and
time by honing in on the most promising policy. Such a behaviour was also
observed in the other exploration strategies in Q-learning (shown in Appendix
C).

Rather than starting with wide traversals, SAC maintains a low search distance
fairly consistently (shown in Figure 7 (c)), and thus, conserving effort but
sacrificing time. The exploration strategy in SAC can hone in on a promising
policy, search within its region for other policies, and only deviate when
nearby policies are similar to the promising one. It is therefore likely to
alternate between equally promising policies provided they are clustered
closely. The exploration strategy in SAC is somewhere between diverse and
precise.

Task (b) 5x5 grid and sparse rewards. We increased the task difficulty by
discarding reward shaping. As seen in panels (c) and (d) of Figure 6, the EI
values, and also convergence times, are an order of magnitude or two higher
compared to the 5x5 dense rewards task.

All the exploration strategies used in the Q-learning algorithm had low
success rates here as seen in Table 1. Since the strategies did not traverse
very far in the policy space, often not reaching optimal policies, they had
low geodesic distances $L$, leading to a large $C/L$. The convergence times in
these cases reached the cut-off limit of learning time. However, the EI $C/L$
continues to be informative about the exploration effort, even when the
algorithm does not succeed, enabling it to be used as a diagnostic tool.

SAC had the highest success rate followed closely by DQN, while the others
have poor success rates. We observe that the convergence time can be
misleading in providing information about exploration efficiency, while the EI
is not. In the case of DQN, its success rate aligns with the EI value while
the convergence time would have suggested it is as poor as others.

As expected, deep learning algorithms SAC and DQN outperformed Q-learning
algorithms, and thus, are consistent with literature Mnih et al. [2013]. The
exploration strategies in Q-learning seem to have the same EI median. However,
decaying $\epsilon$-greedy has data points concentrated around the median,
while others seem to be widely distributed. The wide distribution of the data
indicates that there are instances where optimal policies were successfully
reached as well as not, contrasting with distributions exhibiting low
dispersion.

Algorithm | [5x5] D | [5x5] S | [15x15] D
---|---|---|---
Decaying $\epsilon$-greedy | 100.0% | 5.0% | 88.7%
$\epsilon$-greedy | 100.0% | 15.0% | 96.2%
Greedy | 100.0% | 13.0% | 100.0%
UCRL | 100.0% | - | 92.5%
Softmax | 100.0% | 15.0% | 92.5%
DQN | 100.0% | 83.0% | 77.4%
SAC | 100.0% | 98.0% | 20.8%
Table 1: The success rates of exploration schemes on tasks (a) 5x5 dense
rewards, (b) 5x5 grid-world sparse rewards, and (d) 15x15 grid-world with
dense rewards. The success criterion was formulated to be satisfied when the
model policy consistently yields maximum returns for the task across a span of
50 successive updates. UCRL was not tested on task (b) since it cannot handle
transient states Jaksch et al. [2010], often found in sparse reward settings.

Task (d) 15x15 grid and dense rewards. Here, we increased the task difficulty
by increasing the state space. We see in Figure 6 (e) and (f), that SAC is the
worst-performing algorithm. Although its EI has a small median, the wide
distribution highlights instances of failure to reach optimal policies. Its
poor performance was unexpected, and in this task, it seems that starting off
with exploring widely has an advantage as algorithms with this exploration
strategy performed well (see Table 1 for the success rates).

In this task, Q-learning with greedy strategy seems to be a superior
algorithm, followed by other Q-learning algorithms. Although UCRL had better
success rate than DQN, it scored worse in EI and convergence time. This means,
UCRL takes more time and effort to consistently attain optimal policies. We
notice that both EI and convergence time cannot always predict success.
However, they are useful in understanding exploration behaviour, particularly
EI.

In summary, we see that convergence time lacks the gradation to concisely
convey the algorithms’ ability to reach the optimal policy, while EI does
(Figure 6 (c)-(d)). Also, with segment length plots, we can analyse the
inherent exploratory behaviour of an RL algorithm for a task.

###  4.4 Effects of Convergence Criterion on Exploratory Index

We set a stricter convergence criterion of 50 consecutive policy updates with
maximal returns in Section 4.2 compared to the more typical convergence over 5
policy updates (Section 4.3), yielding EI values for 5x5 dense rewards task
for the decaying $\epsilon$-greedy that do not match across the two cases as
seen in Figures 5 and 6. Here, we test all the aforementioned RL algorithms on
the 5x5 dense rewards task (a) across these two convergence criteria.

Figure 8: Plots demonstrate the relationship between EI and convergence time
of various algorithms across convergence criteria of varying strictness. (a)
and (b) are for a span of 5 consecutive policy model updates yielding maximum
returns, while (c) and (d) are for 50. The task is a 5x5 gridworld with dense
rewards.

In Figure 8(c) and (d), the EI strongly correlates with convergence time when
the convergence criterion is stricter, but doesn’t correlate as much under the
less restrictive typical convergence criterion (Figure 8(a) and (b)).

With a weaker criterion, a diverse exploration strategy with high EI may reach
an optimal policy fast i.e. with low convergence time. With a stricter
criterion, it may need to meander in and around the space of optimal policies
till it settles stably within the optimal region, leading to a higher
convergence time and greater correlation between EI and convergence time. A
more precise strategy with lower EI will in any case settle quicker as
discussed in Section 4.3. These behaviours are illustrated in Figure 9, where
$C_{1}$ and $C_{2}$ denote trajectories in the OTDD space taken by a diverse
exploration strategy (e.g. UCRL) and a precise exploration scheme (e.g. decay
$\epsilon$-greedy), respectively.

Figure 9: Illustration of exemplary exploration behaviours of exploration
strategies with diverse and precise characteristics. Curve $C_{1}$ shows how
typically diverse exploration behaves in comparison to precise exploration
represented by curve $C_{2}$.

##  5 Discussion

Our work is an attempt to theoretically and quantitatively understand
exploration in RL. We tackle this by introducing Exploration Index (EI), which
contrasts learning in RL against learning in SL. This helps us identify
exploration efficiency as the relative effort of knowledge transfer involved
during learning in RL. With EI, we are able to compare exploration strategies
from both tabular and non-tabular RL algorithms. Using this concept, we were
able to study and extract insights about the exploration behaviour of various
schemes and also witness the effect of the strictness of convergence criterion
on exploration. Not only is that EI is based on sound theory, but we also
provide several empirical evaluations to demonstrate its usefulness in
practice.

It should be noted that EI presupposes that for any RL algorithm, exploration
is the primary source of stochasticity contributing to its path within the
OTDD space. This assumption is true for deterministic RL algorithms but may
require careful consideration for algorithms that include other sources of
noise, such as stochastic gradient descent (e.g. SAC and DQN). Future work
could assess when exploration is the dominant source of stochasticity in non-
deterministic RL algorithms, or various mitigations otherwise.

Briefly list author contributions. This is a nice way of making clear who did
what and to give proper credit. This section is optional.

H. Q. Bovik conceived the idea and wrote the paper. Coauthor One created the
code. Coauthor Two created the figures.

###### Acknowledgements.

R. Nkhumise thanks Pawel Pukowski for helpful discussions on the theoretical
aspects. R. Nkhumise was supported by the EPSRC Doctoral Training Partnership
(DTP) - Early Career Researcher Scholarship. A. Gilra and D. Basu were
supported under the CHIST-ERA grant (CHIST-ERA-19-XAI-002), by the Engineering
and Physical Sciences Research Council, United Kingdom (grant reference
EP/V055720/1), and L’Agence Nationale de la Recherche, France (grant reference
ANR-21-CHR4-0007) respectively, as part of the Causal Explanations in
Reinforcement Learning (CausalXRL) project.

## References

  * Achille et al. [2019] A. Achille, M. Lam, R. Tewari, A. Ravichandran, S. Maji, C. C. Fowlkes, S. Soatto, and P. Perona.  Task2vec: Task embedding for meta-learning.  _Proceedings of the IEEE/CVF International conference on computer vision_ , pages 6430––6439, 2019.
  * Alvarez-Melis and Fusi [2020] D. Alvarez-Melis and N. Fusi.  Geometric dataset distances via optimal transport.  In _Proceedings of the 34th International Conference on Neural Information Processing Systems_ , Red Hook, NY, USA, 2020. Curran Associates Inc.  ISBN 9781713829546.
  * Amari [2010] S. Amari.  Information geometry in optimization, machine learning and statistical inference.  _Frontiers of Electrical and Electronic Engineering in China_ , 5:241––260, 2010.
  * Amin et al. [2021] S. Amin, M. Gomrokchi, H. Satija, H. van Hoof, and D. Precup.  A survey of exploration methods in reinforcement learning, 2021.
  * Auer and Ortner [2006] P. Auer and R. Ortner.  Logarithmic online regret bounds for undiscounted reinforcement learning.  In B. Schölkopf, J. Platt, and T. Hoffman, editors, _Advances in Neural Information Processing Systems_ , volume 19. MIT Press, 2006.  URL https://proceedings.neurips.cc/paper_files/paper/2006/file/c1b70d965ca504aa751ddb62ad69c63f-Paper.pdf.
  * Basu et al. [2020] D. Basu, P. Senellart, and S. Bressan.  Belman: An information-geometric approach to stochastic bandits.  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2019, Würzburg, Germany, September 16–20, 2019, Proceedings, Part III_ , pages 167–183. Springer, 2020.
  * Bellemare et al. [2016] M. Bellemare, S. Srinivasan, G. Ostrovski, T. Schaul, D. Saxton, and R. Munos.  Unifying count-based exploration and intrinsic motivation.  In _Advances in Neural Information Processing Systems_ , pages 1471–1479. Curran Associates, Inc., 2016.
  * Ben-David et al. [2006] S. Ben-David, J. Blitzer, K. Crammer, and F. Pereira.  Analysis of representations for domain adaptation.  _Advances in neural information processing systems_ , 2006.
  * Burda et al. [2019] Y. Burda, H. Edwards, A. Storkey, and O. Klimov.  Exploration by random network distillation.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=H1lJJnR5Ym.
  * Christodoulou [2019] P. Christodoulou.  Soft actor-critic for discrete action settings, 2019.
  * Conserva and Rauber [2022] M. Conserva and P. Rauber.  Hardness in markov decision processes: Theory and practice.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, _Advances in Neural Information Processing Systems_ , volume 35, pages 14824–14838. Curran Associates, Inc., 2022.
  * Eysenbach et al. [2019] B. Eysenbach, A. Gupta, J. Ibarz, and S. Levine.  Diversity is all you need: Learning skills without a reward function.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=SJx63jRqFm.
  * Flamary et al. [2021] R. Flamary, N. Courty, A. Gramfort, M. Z. Alaya, A. Boisbunon, S. Chambon, L. Chapel, A. Corenflos, K. Fatras, N. Fournier, L. Gautheron, N. T.H. Gayraud, H. Janati, A. Rakotomamonjy, I. Redko, A. Rolet, A. Schutz, V.Seguy, D. J. Sutherland, R.Tavenard, A. Tong, and T. Vayer.  Pot: Python optimal transport.  _Journal of Machine Learning Research_ , 22(78):1–8, 2021.  URL http://jmlr.org/papers/v22/20-451.html.
  * Gao and Chaudhari [2021] Y. Gao and P. Chaudhari.  An information-geometric distance on the space of tasks.  _Proceedings of the 38th International Conference on Machine Learning_ , 139, 2021.
  * Haarnoja et al. [2018] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine.  Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor.  In _Proceedings of the 35th International Conference on Machine Learning_ , 2018.
  * Haarnoja et al. [2019] T. Haarnoja, A. Zhou, K. Hartikainen, G. Tucker, S. Ha, J. Tan, V. Kumar, H. Zhu, A. Gupta, P. Abbeel, and S. Levine.  Soft actor-critic algorithms and applications, 2019.
  * Hollenstein et al. [2020] J. J. Hollenstein, S. Auddy, M. Saveriano, E. Renaudo, and J. Piater.  How do offline measures for exploration in reinforcement learning behave?, 2020.
  * Hussein et al. [2017] A. Hussein, M. M. Gaber, E. Elyan, and C. Jayne.  Imitation learning: A survey of learning methods.  _ACM Computing Surveys_ , 50:1––35, 2017.
  * Jaksch et al. [2010] T. Jaksch, R. Ortner, and P. Auer.  Near-optimal regret bounds for reinforcement learning.  _Journal of Machine Learning Research_ , 11(51):1563–1600, 2010.  URL http://jmlr.org/papers/v11/jaksch10a.html.
  * Kingma and Ba [2017] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization, 2017.
  * Ladosz et al. [2022] P. Ladosz, L. Weng, M. Kim, and H. Oh.  Exploration in deep reinforcement learning: A survey.  _Information Fusion_ , 85:1–22, 2022.
  * Lott [2008] J. Lott.  Some geometric calculations on wasserstein space.  _Communications in Mathematical Physics_ , 277:423––437, 2008.
  * Maroti [2019] A. Maroti.  RBED: Reward based epsilon decay, 2019.
  * Mnih et al. [2013] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller.  Playing atari with deep reinforcement learning, 2013.
  * Nguyen et al. [2020] C. V. Nguyen, T. Hassner, M. Seeger, and C. Archambeau.  Leep: a new measure to evaluate transferability of learned representations.  In Hal Daume III and Aarti Singh, editors, _Proceedings of the 37th International Conference on Machine Learning (ICML 2020)_ , pages 7294–7305. PMLR, 2020.
  * Osband et al. [2016] I. Osband, B. V. Roy, and Z. Wen.  Generalization and exploration via randomized value functions.  In _Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48_ , ICML’16, pages 2377––2386. JMLR.org, 2016.
  * Peyré and Cuturi [2019] G. Peyré and M. Cuturi.  Computational optimal transport.  _Foundations and Trends in Machine Learning_ , 11(5-6):355–607, 2019.
  * Tan et al. [2021] Y. Tan, Y. Li, and S. L. Huang.  OTCE: A transferability metric for cross-domain cross-task representations.  _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021.
  * Tang et al. [2017] H. Tang, R. Houthooft, D. Foote, X. Chen A. Stooke, Y. Duan, J. Schulman, F. De Turck, and P. Abbeel.  $\\#$exploration: A study of count-based exploration for deep reinforcement learning.  In _Advances in Neural Information Processing Systems_ , pages 2753–2762. Curran Associates, Inc., 2017.
  * Tran et al. [2019] A. T. Tran, C. V. Nguyen, and T. Hassner.  Trasnferability and hardness of supervised classification tasks.  _Proceedings of the IEEE/CVF International conference on computer vision_ , pages 1395––1405, 2019.
  * Watkins and Dayan [1992] Christopher JCH Watkins and Peter Dayan.  Q-learning.  _Machine learning_ , 8:279–292, 1992.
  * Wei [2020] P. Wei.  Exploration-exploitation strategies in deep q-networks applied to route-finding problems.  _Journal of Physics Conference Series_ , 1684:012073, 2020.
  * Zamir et al. [2018] A. R. Zamir, A. Sax, W. Shen, L. J. Guibas, J. Malik, and S. Savarese.  Taskonomy: Disentangling task transfer learning.  _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 3712––3722, 2018.
  * Zhang et al. [2019] L. Zhang, K. Tang, and X. Yao.  Explicit planning for efficient exploration in reinforcement learning.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32, pages 7488–7497. Curran Associates, Inc., 2019.

Measuring Exploration in Reinforcement Learning
via Optimal Transport in Policy Space (Supplementary Material)

This Supplementary Material should be submitted together with the main paper.

##  Appendix A Proof Details

Proof of Proposition 3.1. Suppose that the environment transitions and the
policy are stochastic, then the trajectory distribution is

$p(\Gamma|\pi)=\rho_{0}(s_{0})\Pi_{t=0}^{T}p(s_{t+1}|s_{t},a_{t})\pi(a_{t}|s_{t}),$
(5)

where $\Gamma$ is the space of possible trajectories, $\tau_{i}\sim
p(\Gamma|\pi)$ and $C(\pi)=\\{\tau_{i}\\}_{i=0}^{n\to\infty}$. If two
trajectory distributions $p(\Gamma|\pi_{1})$ and $p(\Gamma|\pi_{2})$ produce
the same dataset $C(\pi)$, then
$p(\Gamma=\tau|\pi_{1})=p(\Gamma=\tau|\pi_{2})$ for all $\tau$. However,
$p(\Gamma|\pi_{1})$ and $p(\Gamma|\pi_{2})$ are only truly distinct
distributions if $\pi_{1}\neq\pi_{2}$. Therefore, there must exist at least
one outcome $\tau$ for which $p(\Gamma=\tau|\pi_{1})\neq
p(\Gamma=\tau|\pi_{2})$. This contradiction proves that two distinct policies
cannot produce the same trajectory class $C(\pi)$, and if $C(\pi)$ is unique
to the policy ${\pi}$, then so is $P_{\pi}(s,a)$.

Proof of Proposition 3.2. On the OTDD space, the true, i.e. unapproximated
distance between consecutive policies $\pi_{i}$ and $\pi_{i+1}$ after an
update is $d_{OT}(P_{\pi_{i}},P_{\pi_{i+1}})$. This is induced by $i^{th}$ RL
policy update. While in supervised learning, we estimate the distance using
datasets of policies, i.e. approximated joint distributions, using
$d_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{i}+1})$. Therefore,

$|d_{OT}(P_{\pi_{i}},P_{\pi_{i+1}})-d_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{i+1}})|\leq\epsilon$
(6)

where $\epsilon$ is some small value. By summing sequentially through policies
encountered during RL training:

$\sum_{i}|d_{OT}(P_{\pi_{i}},P_{\pi_{i+1}})-d_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{i+1}})|\leq
k\epsilon$ (7)

Since from reverse triangle inequality we know,

$\left|\sum_{i}d_{OT}(P_{\pi_{i}},P_{\pi_{i+1}})-\sum_{i}d_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{i+1}})\right|\leq\sum_{i}|d_{OT}(P_{\pi_{i}},P_{\pi_{i+1}})-d_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{i+1}})|$
(8)

Therefore

$\left|\sum_{i}d_{OT}(P_{\pi_{i}},P_{\pi_{i+1}})-\sum_{i}d_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{i+1}})\right|\leq
k\epsilon$ (9)

This means as the size of the dataset increases $\hat{P}_{\pi}\to P_{\pi}$,
then
$\sum_{i}d_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{i+1}})\to\sum_{i}d_{OT}(P_{\pi_{i}},P_{\pi_{i+1}})$.

##  Appendix B Methods

###  B.1 Results in Figures

Model parameter initialisation. We initialised model parameters for deep
learning RL algorithms like DQN and SAC by uniformly sampling weight values
between $-3\cdot 10^{-4}$ and $3\cdot 10^{-4}$ and the biases at $0$. For
tabular Q-learning algorithms, we randomly initialized the Q-values between
$-1.0$ and $1.0$.

Results in Figure 4. To determine results in Figure 4, we worked with
$\epsilon=\\{0.05,0.2,0.4,0.9\\}$. For each $\epsilon$ value we ran 40
training trials and quantified the numerical data as displayed in the figures.
With decaying $\epsilon$-greedy, we begin with the given $\epsilon$-value and
exponentially using Equation 10 below. The convergence criterion was met when
the policy model consistently produced maximal returns over 50 updates.

$\epsilon=max\\{0.9999\cdot\epsilon,0.0001\\}$ (10)

Results in Figure 5. Q-learning with decaying $\epsilon$-greedy was employed
on the grid-world tasks (Figure 3) for a convergence criterion of 50 model
updates with maximum returns. We selected $\epsilon=0.9$ and aggregated the
result on over 40 training trials.

Results in Figure 6 and Table 1 Q-learning with greedy has its Q-values
updated based on a greedy (exploitative) strategy. For panels (a) and (b), all
algorithms achieved $100\%$ success rate unlike in other panels. The success
rate was computed based on the number of training trials each algorithm met
the convergence criterion of producing maximum returns over 5 consecutive
updates. For $\epsilon$-greedy and decaying $\epsilon$-greedy $\epsilon=0.9$
was utilized, and $\tau=2.0$ and $\delta=0.1$ for Softmax(Boltzmann) and UCRL,
respectively. The $\alpha$ parameter for SAC was autotuned using the approach
in [Haarnoja et al., 2019] along with hyperparameters described in Table 2.
While DQN began with $\epsilon=1.0$ and the value decayed using Equation 10.
Table 13 shows hyperparameters for DQN.

For panels (c) - (f), there were instances when the geodesic length between
the initial and final policies
$L=d_{OT}(\mathcal{D}_{\pi_{0}},\mathcal{D}_{\pi_{N}})$ was equal to 0. This
happens usually occurs when the algorithm could not reach an optimal policy
and randomly landed on the initial policy. To avoid infinite EI values we
added $10^{-3}$ on all $L$ values in our implementation. This provides
reasonable large numbers.

Table 2: SAC Hyperparameters. Parameter | Value
---|---
optimizer |  Adam Kingma and Ba [2017]
learning rate | $5\cdot 10^{-4}$
discount($\gamma$) | 0.99
replay buffer size | $10^{4}$
number of hidden layers (all networks) | 1
number of hidden units per layer | 32
number of samples per minibatch | 64
nonlinearity | ReLU
entropy target | -4
target smoothing coefficient ($\tau$) | 0.01
target update interval | 1
gradient steps | 1
initial exploration steps |
before model starts updating | 500
Table 3: DQN Hyperparameters. Parameter | Value
---|---
optimizer |  Adam Kingma and Ba [2017]
learning rate | $5\cdot 10^{-2}$
discount($\gamma$) | 0.99
replay buffer size | $10^{4}$
number of hidden layers (all networks) | 1
number of hidden units per layer | 32
number of samples per minibatch | 64
nonlinearity | ReLU
target smoothing coefficient ($\tau$) | 0.001
target update interval | 1
gradient steps | 1
initial exploration steps |
before $\epsilon$ decays | 1500

Results in Figure 7 The plots show segments length
$ds=d_{OT}(P_{\theta},P_{\theta+d\theta})$ vs time-steps. While the cumulative
plots show accumulation of the segment lengths over time. The results in
segment plots were smoothed using exponential smoothing with smoothing factor
of 0.9. The plots are were made, similarly to Figure 6, with over 80 training
trials for the 5x5 dense rewards task.

Results in Figure 8 We used the same approach as in Figure 6, however the
convergence criteria were 5 and 50, respectively. The task was a 5x5 grid-
world with dense rewards.

###  B.2 Calculations

OTDD. To compute the OTDD given that we have collected datasets
$\mathcal{D}_{\pi}=\\{(s_{(i)},a_{(i)})\\}_{i=0}^{m}$, we define a set of
states with action $a$ in $\mathcal{D}_{\pi}$ as
$\mathcal{N}_{\mathcal{D}_{\pi}}(a)=\\{s\in\mathcal{S}\mid(s,a)\in\mathcal{D}\\}$
with cardinality $n_{a}$. This is used as a finite sample to define
distribution that encodes action $a$ as

$\alpha_{a}=\frac{1}{n_{a}}\sum_{s\in\mathcal{N}_{\mathcal{D}_{\pi}}(a)}\delta_{s}$
(11)

We then precompute and store in memory action-to-action pairwise Wasserstein
distances $W_{2}(\alpha_{a},\alpha_{a^{\prime}})$ using an off-the-shelf
optimal transport solver toolbox POT [Flamary et al., 2021]. The cost function
used was L1-norm (cityblock or Manhattan distance) since the tasks are grid-
worlds restricted to vertical and horizontal movements between states.
Following Equation (2), OTDD between two policy distributions is

$d^{2}_{OT}(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{j}})=\min_{\gamma\in\Gamma(\hat{P}_{\pi_{i}},\hat{P}_{\pi_{j}})}\int_{\mathcal{S}\times\mathcal{A}}\left(\lVert
s-s^{\prime}\rVert^{2}_{1}+W^{2}_{2}(\alpha_{a},\alpha_{a^{\prime}})\right)d\gamma((s,a),(s^{\prime},a^{\prime}))$
(12)

Equation (12) is also solved using the toolbox POT, while action-to-action
pairwise distances are retrieved from memory on-demand. Due to the simplicity
of our tasks, an Earth Movers distance algorithm [Peyré and Cuturi, 2019] in
OT is used.

##  Appendix C Additional experimental results

This section encompasses segment length and cumulative length plots of
Q-learning with $\epsilon$-greedy (Figure 10), decaying $\epsilon$-greedy
(Figure 11), and Softmax (Figure 12) strategies, along with DQN (13).

Figure 10:  These plots show the behaviour of the segment length over time for
Q-learning with $\epsilon$-greedy strategy. The vertical dotted line
represents the mean convergence time averaged over 80 runs. For segment length
plots, the solid line represents a smoothed rendition of the noisy mean curve,
which is shaded. While for cumulative length plots, the solid line depicts the
mean curve and the shaded areas indicate the range of standard deviation
values. Figure 11:  These plots show the behaviour of the segment length over
time for Q-learning with decaying $\epsilon$-greedy strategy. Figure 12:
These plots show the behaviour of the segment length over time for Q-learning
with softmax strategy. Figure 13:  These plots show the behaviour of the
segment length over time for DQN.