# Ro-vibrational energy and thermodynamic properties of molecules subjected to
Deng-Fan potential through an improved approximation

Debraj Nath , Amlan K. Roy Department of Mathematics, Vivekananda College,
Thakurpukur, Kolkata-700063, India. Email<EMAIL_ADDRESS>of
Chemical Sciences, Indian Institute of Science Education and Research (IISER)
Kolkata, Nadia, Mohanpur-741246, WB, India. Email<EMAIL_ADDRESS><EMAIL_ADDRESS>dedicates this article to the loving memory of his kind-
hearted teacher, late Prof. Hrishikesh Pradhan, P. K. College, Contai. Through
his relentless, selfless service, he has influenced and inspired many
chemists, over the years.

###### Abstract

Accurate solution of the Schrödinger equation with Deng-Fan potential is
presented by means of Nikiforov-Uvarov method. A modified Pekeris-type
approximation is proposed for the centrifugal term, from a linear combination
of the $r\to 0$ and $r\to r_{e}$ limits. It can potentially offer a series of
approximations (depending on an adjustable parameter $\lambda$). The existing
approximations in the literature can then be recovered in certain special
cases. Its efficiency and feasibility is demonstrated by a critical comparison
of eigenvalues produced at various $\lambda$’s for four molecules, _viz._ ,
H2, LiH, HCl and CO. Analytical expressions are derived for energies,
eigenfunctions and the thermodynamic properties such as vibrational mean free
energy, vibrational free energy, vibrational entropy and vibrational specific
heat. The effect of quantum correction on partition function and thermodynamic
properties is discussed by including the correction up to 10th-order, for H2
and LiH. The effect of $\lambda$ parameter on these properties is also
studied.

Key words: Deng-Fan potential, Nikiforov-Uvarov Method, Greene-type
approximation, Pekeris-type approximation, ro-vibrational energy, vibrational
partition function.

##  1 Introduction

A correct and proper representation of potential energy as a function of
internuclear distance is an essential and indispensable tool for accurate
understanding of structure and dynamics of a molecular system. An ideal
potential energy function should show proper limiting behavior, i.e.,
$V(0)=\infty$ and approach a constant at infinite distance. Analytical
eigenvalues and eigenfunctions are desirable, as they facilitate accurate
estimation of transition frequencies and matrix elements. The oldest and
arguably the most celebrated empirical Morse potential [1] was suggested
almost ninety years ago. This widely used simple analytical functional has
found valuable applications [2] in several branches in physics and chemistry
including ro-vibrational states of diatomic molecules, adsorption of
atoms/molecules on solid surface, deformation of cubic metals etc. The three-
parameter exponential function is represented by,

$U(r)=D_{e}\left[1-\exp^{-\alpha(r-r_{e})}\right]^{2},$ (1)

where $D_{e}$ denotes dissociation energy, $r_{e}$ corresponds to equilibrium
internuclear distance and the Morse parameter $\alpha$ is defined in terms of
force constant $k_{e}$ as $\alpha=\sqrt{k_{e}/2D_{e}}$. Despite its remarkable
success, the potential has a finite value, as internuclear distance approaches
zero, and does not lead to the expected infinity nature. Thus when employed
for the vibrational levels of diatomic molecule, the Morse oscillator shows
considerable deviation from experimentally observed values. There has been
many attempts to address these issues leading to a variety of proposals in the
literature over the years. Some of the most prominent and notable ones are
listed as follows: Kratzer [3, 4, 5], Rosen-Morse [6], Manning-Rosen [7, 8, 9,
10], Pöschl-Teller [11, 12, 13], Linnett [14], Hulthén [15, 16, 17, 18],
Lippincott [19], shifted/Deng-Fan (DF) [20, 2, 21, 22, 23, 24, 25], Tietz-Hua
[26, 27, 28, 29], Schlöberg [30], Zavitsas [31], deformed Rosen-Morse [32, 33,
34], Woods-Saxon [35, 36, 37], pseudoharmonic [38, 39], Hajigeorgiou [40].

The current communication deals with the DF potential for diatomic molecules
put forth about 70 years ago. This three-parameter function given in the
following expression

$V(r)=D_{e}\left(1-\frac{b}{e^{ar}-1}\right)^{2},\ \ \ \ b=e^{ar_{e}}-1,\ \ \
\ r\in(0,\infty),$ (2)

has aroused considerable interest in last two decades. Here the three positive
parameters $D_{e},r_{e},a$ represent dissociation energy, equilibrium
internuclear distance and radius of the potential well respectively. As
internuclear distance approaches the limiting values zero and infinity, it
offers a qualitatively correct asymptotic nature. This is also often referred
as Generalized Morse potential due to its qualitative resemblance to Morse
potential [41, 2]. This is also connected to another important diatomic
potential due to Manning-Rosen [7].

Like the case of many other practical and realistic potentials, the relevant
Schrödinger equation with centrifugal term can not be solved analytically for
eigenvalues and eigenfunctions. This has lead to the development of a number
of theoretical methodologies of approximation to the centrifugal term in both
non-relativistic and relativistic case. Some of the prominent ones may be
mentioned. The authors in [41] have discussed the exact solvability problem
via an $SO(2,2)$ symmetry algebra. The eigenvalues and eigenfunctions for
$\ell=0$ states were studied by an algebraic method [42]. Approximate
analytical solutions for rotating DF potential for arbitrary $n,\ell$ quantum
numbers were given in terms of generalized hyper-geometric functions
${}_{2}F_{1}(a,b;c;z)$ [43]. An improved approximation scheme for the
centrifugal term, along with a super-symmetric shape invariance approach was
employed in [44]; analytical solution for Dirac equation has been presented
with latter approach [45]. Solutions of Klein-Gordon equation with spinless
particle for $\ell\neq 0$ states were considered in [46]. A Pekeris
approximation for centrifugal term and a Nikiforov-Uvarov (NU) method was
employed for Dirac, Klein-Gordon and Schrödinger equation [47, 48]. Some other
methods include a numerical generalized pseudospectral (GPS) scheme [22],
exact quantization rule [49, 23], a Feynman path integral formalism for
D-dimensional problem [25], etc. A comparative analysis for the performance of
Morse, Manning-Rosen, Schlöberg and DF potential in the context of diatomic
molecules has been offered in [50].

Our objective in this communication is two-fold. At first, we present accurate
eigenvalues and eigenfunctions of the DF potential by means of a simple novel
approximation for the centrifugal term, within the framework of NU method.
This is accomplished by taking a combination of approximations at $r\to 0$ and
$r\to r_{e}$ limits. Its performance is critically analyzed by presenting the
eigenvalues at various values of tunable parameter. These are compared with
existing approximations available in the literature. The calculated ro-
vibrational eigenvalues are offered for two representative diatomic molecules
(H2, LiH). Both circular ($s$) as well as non-zero-$\ell$ states are obtained
accurately. Then in the second stage, analytical expressions for vibrational
partition function and related thermodynamic properties such as internal
energy, free energy, entropy, specific heat capacity are derived. They are
examined in terms of (i) temperature and (ii) the impact of quantum correction
on these. The article is organized as follows. Section 2.1 provides the
relevant details pertaining to the new approximation. The necessary
expressions for eigenvalues and eigenfunctions within the NU method are
detailed in Sec. 2.2. The associated energy spectrum is analyzed in Sec. 2.3.
The expressions for thermodynamical quantities are given and discussed in Sec.
3 for a representative set of four diatomic molecules (H2, LiH, HCl and NO).
Additionally the effect of quantum correction on the thermal quantities are
discussed as well. Finally a few comments are made in Sec. 4.

##  2 The methodology

###  2.1 Analytical solutions of DF potential with centrifugal term within a
new approximation

The time-independent non-relativistic Schrödinger equation for a diatomic
molecule in presence of DF potential in Eq. (2), can be written as,

$-\frac{\hbar^{2}}{2\mu}\nabla^{2}\psi({\bf r})+V({\bf r})\psi({\bf
r})=E\psi({\bf r}),$ (3)

where

$\nabla^{2}=\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial}{\partial
r}\right)+\frac{1}{r^{2}\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right)+\frac{1}{r^{2}\sin^{2}\theta}\frac{\partial^{2}}{\partial\phi^{2}},$
(4)

with $\mu$ and $\hbar$ denoting the reduced mass and Planck constant
respectively. Let us define,

$\psi({\bf
r})=\frac{1}{\sqrt{2\pi}}\frac{R(r)}{r}Y_{l}^{m}(\theta)\,e^{im_{l}\phi},$ (5)

where $R(r)$ and $Y_{l}^{m_{l}}$ signify the usual radial and angular
functions, while $l,m_{l}$ represent the angular quantum numbers. Then one
obtains the following radial Schrödinger equation,

$\frac{d^{2}R(r)}{dr^{2}}+\left[\frac{2\mu}{\hbar^{2}}E-\frac{2\mu}{\hbar^{2}}V(r)-\frac{l(l+1)}{r^{2}}\right]R(r)=0,$
(6)

where $l(l+1)$ is the separation constant. Now one can expand (Pekeris-type
approximation) $\frac{1}{r^{2}}$ as below,

$\begin{array}[]{ll}\frac{1}{r^{2}}&=\frac{f(x)}{r_{e}^{2}}=\frac{1}{r_{e}^{2}}\displaystyle\sum\limits_{n=0}^{\infty}\frac{x^{n}}{n!}\left[\frac{d^{n}f(x)}{dx^{n}}\right]_{x=0}\approx\frac{1}{r_{e}^{2}}\displaystyle\sum\limits_{n=0}^{2}\frac{x^{n}}{n!}\left[\frac{d^{n}f(x)}{dx^{n}}\right]_{x=0}+O(x^{3}),\end{array}$
(7)

where

$f(x)=c_{0}+\frac{c_{1}\,s}{1-s}+\frac{c_{2}\,s^{2}}{\left(1-s\right)^{2}},$
(8)

and

$s=e^{-\alpha r},r=r_{e}(1+x),$ (9)

where $x$ is a dimensionless variable. The coefficients $c_{0},c_{1}$ and
$c_{2}$ are dimensionless parameters. In order to find these coefficients for
non-zero $l$, one may express $\frac{1}{r^{2}}$ as [51, 27, 52],

$\begin{array}[]{ll}\frac{1}{r^{2}}&=\frac{1}{r_{e}^{2}}\displaystyle\sum\limits_{n=0}^{\infty}(-1)^{n}(n+1)x^{n},\\\
&\approx\frac{1}{r_{e}^{2}}\left(1-2x+3x^{2}+O(x^{3})\right).\end{array}$ (10)

Comparing the coefficients of Eq. (7) and Eq. (10), [52, 53, 25], one obtains,

$\begin{array}[]{ll}c_{0}&=\frac{1}{u^{2}}\left[3-3u+u^{2}+\left(2u-6\right)s_{e}+\left(u+3\right)s_{e}^{2}\right],\\\
c_{1}&=\frac{2}{s_{e}u^{2}}\left(1-s_{e}\right)^{2}\left((3+u)s_{e}+2u-3\right),\\\
c_{2}&=-\frac{1}{s_{e}^{2}u^{2}}\left(1-s_{e}\right)^{3}\left((3+u)s_{e}+u-3\right),\\\
s_{e}&=e^{-u},~{}u=\alpha r_{e}.\end{array}$ (11)

Let us now consider another approximation (Greene-type) [54, 23, 55, 21, 48]
of the form,

$\begin{array}[]{ll}\frac{1}{r^{2}}&\approx\alpha^{2}\left(\frac{1}{12}+\frac{s}{1-s}+\frac{s^{2}}{(1-s)^{2}}\right),\\\
\end{array}$ (12)

If we express
$\alpha^{2}\left(\frac{1}{12}+\frac{s}{1-s}+\frac{s^{2}}{(1-s)^{2}}\right)$ in
the explicit form of $r$, then we obtain
$\alpha^{2}\left(\frac{1}{12}+\frac{s}{1-s}+\frac{s^{2}}{(1-s)^{2}}\right)\approx\frac{1}{r^{2}}+\alpha^{2}\left(\frac{(\alpha
r)^{2}}{240}-\frac{(\alpha r)^{4}}{6048}+\frac{(\alpha
r)^{6}}{172800}+O((\alpha r)^{8})\right)$. Hence it is better than Eq. (7),
when $\alpha r\ll 1$. In the limiting case,

$\begin{array}[]{ll}\frac{1}{r^{2}}&=\lim\limits_{\alpha\rightarrow
0}\alpha^{2}\left(\frac{1}{12}+\frac{s}{1-s}+\frac{s^{2}}{(1-s)^{2}}\right).\end{array}$
(13)

Therefore, for a fixed $\alpha$, Eq. (12) is a good approximation near $x=-1$
i.e., around $r=0$ and Eq. (7) is good near $x=0$ i.e., in region close to
$r=r_{e}$. This prompts us to consider a linear combination of above two
approximations, as in the following,

$\frac{1}{r^{2}}\approx\alpha^{2}\left(d_{0}(\lambda)+\frac{d_{1}(\lambda)\,s}{1-s}+\frac{d_{2}(\lambda)\,s^{2}}{(1-s)^{2}}\right),$
(14)

where

$d_{0}(\lambda)=\frac{\lambda}{12}+\frac{(1-\lambda)c_{0}}{u^{2}},~{}d_{1}(\lambda)=\lambda+\frac{(1-\lambda)c_{1}}{u^{2}},~{}d_{2}(\lambda)=\lambda+\frac{(1-\lambda)c_{2}}{u^{2}}.$
(15)

From the above equation, it is observed that, $\lambda=0$ implies Eq. (7) and
$\lambda=1$ signifies Eq. (12). Therefore, one can use Eq. (14) as a more
general form of approximation for the centrifugal term $\frac{1}{r^{2}}$. As
$r$ increases, we can consider Eq. (14) for negative value of $\lambda$. In
Fig. 1, the centrifugal term $\frac{1}{r^{2}}$ is plotted along with three
approximations defined in Eqs. (7), (12) and (14), for H2. The panels (A), (B)
refer to the corresponding functions in $(-)$ ve and $(+)$ ve domains of $x$.
For $0\leq\lambda\leq 1$, the approximation of Eq. (14) is a convex
combination of Eqs. (7) and (12), and the corresponding centrifugal term is a
good approximation for $0\leq r\leq r_{e}$. For negative values of $\lambda$,
Eq. (14) is a good approximation to the centrifugal term for $r>r_{e}$, which
is clear from Fig. 1(B). It is noticed that the approximation defined in Eq.
(14) is the best amongst them. This is generally found to be true for other
molecules as well.

Figure 1:  Plot of the centrifugal term $\frac{1}{r^{2}}$ and the three
approximations given in Eqs. (7), (12) and (14) for H2 molecule. The two
situations $\lambda=0,1$ correspond to Eqs. (7) and (12) respectively. Panels
(A), (B) refer to $(-)$ ve and $(+)$ ve regions of $x$.

Using the transformation in Eq. (9) and approximation of Eq. (14) in Eq. (6)
we obtain,

$\begin{array}[]{l}\frac{d^{2}R}{ds^{2}}+\frac{1-s}{s(1-s)}\frac{dR}{ds}+\frac{-As^{2}+Bs-C}{s^{2}(1-s)^{2}}R=0,\end{array}$
(16)

where

$\begin{array}[]{ll}A&=\bar{\varepsilon}^{2}+2b\beta^{2}+b^{2}\beta^{2}+\gamma(d_{2}-d_{1}),\\\
B&=2\bar{\varepsilon}^{2}+2b\beta^{2}-\gamma d_{1},\\\
C&=\bar{\varepsilon}^{2},\\\ \bar{\varepsilon}^{2}&=\beta^{2}+\gamma
d_{0}-\varepsilon,\\\ \beta^{2}&=\frac{2\mu D_{e}}{\alpha^{2}\hbar^{2}},\\\
\varepsilon&=\frac{2\mu E}{\alpha^{2}\hbar^{2}},\\\
\gamma&=l(l+1).\end{array}$ (17)

Therefore, it reduces to an energy-dependent problem of a second-order
differential equation. In the next subsection, we provide the solution of Eq.
(16) by making use of NU method.

###  2.2 Nikiforov-Uvarov Method

One can compare Eq. (16) with the following Eq. [56, 57, 58, 21, 59],

$\frac{d^{2}R}{ds^{2}}+\frac{\widetilde{\tau}(s)}{\sigma(s)}\frac{dR}{ds}+\frac{\widetilde{\sigma}(s)}{\sigma^{2}(s)}R=0,$
(18)

where

$\begin{array}[]{ll}\widetilde{\tau}(s)&=1-s,~{}\sigma(s)=s(1-s),~{}\widetilde{\sigma}(s)=-A\,s^{2}+B\,s-C.\end{array}$
(19)

Let us consider the solution of Eq. (18) in following form,

$R(s)=\phi(s)y(s).$ (20)

Then one may obtain [56],

$\sigma(s)y^{{}^{\prime\prime}}(s)+\tau(s)y^{\prime}(s)+\nu\,y(s)=0,$ (21)

and

$\begin{array}[]{l}\phi(s)=\displaystyle
e^{\int\frac{\pi(s)}{\sigma(s)}ds},\end{array}$ (22)

where

$\begin{array}[]{ll}\pi(s)&=\frac{\sigma^{\prime}(s)-\widetilde{\tau}(s)}{2}\pm\sqrt{\left(\frac{\sigma^{\prime}(s)-\widetilde{\tau}(s)}{2}\right)^{2}-\widetilde{\sigma}(s)+k\,\sigma(s)},\\\
\tau(s)&=\widetilde{\tau}(s)+2\pi(s),\\\ \nu&=k+\pi^{\prime}(s).\end{array}$
(23)

Here $\nu$ and $k$ are real constants independent of $s$. Since $\pi(s)$ is a
polynomial in $s$,
$\left(\frac{\sigma^{\prime}(s)-\widetilde{\tau}(s)}{2}\right)^{2}-\widetilde{\sigma}(s)+k\,\sigma(s)$
is a square of a polynomial in $s$, and the solutions of Eq. (21) are given
by,

$y_{n}(s)=\frac{1}{\rho(s)}\frac{d^{n}}{ds^{n}}\left[\sigma^{n}(s)\rho(s)\right],$
(24)

where

$\begin{array}[]{l}\rho(s)=\displaystyle\left[\sigma(s)\right]^{-1}\,e^{\int\frac{\tau(s)}{\sigma(s)}ds}.\end{array}$
(25)

Therefore, the corresponding eigenvalues are obtained from the relation,

$\nu_{n}=-n\,\tau^{\prime}(s)-\frac{n(n-1)}{2}\sigma^{{}^{\prime\prime}}(s)=-n\,\tau^{\prime}(s)+n(n-1),~{}n=0,1,2,\dots.$
(26)

According to NU method, we obtain a pair of possible values of $k$, which are
given by,

$\begin{array}[]{ll}k_{\pm}&=-\gamma
d_{1}+2b\beta^{2}\pm\bar{\varepsilon}\sqrt{4b^{2}\beta^{2}+4\gamma
d_{2}+1}.\end{array}$ (27)

Now using Eq. (27), we obtain a set of possible values of $\pi(s)$ as,

$\begin{array}[]{ll}\pi(s)&=-\frac{s}{2}\pm\left\\{\begin{array}[]{lll}\left(\bar{\varepsilon}-\frac{1}{2}\sqrt{4b^{2}\beta^{2}+4\gamma
d_{2}+1}\right)s+\sqrt{C},&k=k_{+},&k_{+}-B>0\\\
\left(\bar{\varepsilon}+\frac{1}{2}\sqrt{4b^{2}\beta^{2}+4\gamma
d_{2}+1}\right)s+\sqrt{C},&k=k_{-},&k_{-}-B>0\\\
\left(\bar{\varepsilon}-\frac{1}{2}\sqrt{4b^{2}\beta^{2}+4\gamma
d_{2}+1}\right)s-\sqrt{C},&k=k_{+},&k_{+}-B<0\\\
\left(\bar{\varepsilon}+\frac{1}{2}\sqrt{4b^{2}\beta^{2}+4\gamma
d_{2}+1}\right)s-\sqrt{C},&k=k_{-},&k_{-}-B<0\end{array}\right\\}.\end{array}$
(28)

Since $\tau^{\prime}(s)<0$, in this work, we have considered $k=k_{-}$, where
$k_{-}-B<0$ and select the following $\pi(s)$,

$\begin{array}[]{ll}\pi(s)&=-\frac{s}{2}-\left(\bar{\varepsilon}+\frac{1}{2}\sqrt{4b^{2}\beta^{2}+4\gamma
d_{2}+1}\right)s+\sqrt{C}.\end{array}$ (29)

In that case, one gets,

$\begin{array}[]{ll}\rho(s)&=s^{2\bar{\varepsilon}}(1-s)^{2L-1},\\\
\phi(s)&=s^{\bar{\varepsilon}}(1-s)^{L},\end{array}$ (30)

and

$\begin{array}[]{ll}y_{n}&=s^{-2\bar{\varepsilon}}(1-s)^{-2L+1}\frac{d^{n}}{dx^{n}}\left[s^{n+2\bar{\varepsilon}}(1-s)^{n+2L-1}\right]=n!\,P_{n}^{(2\bar{\varepsilon},2L-1)}(1-2s),\\\
&=\frac{\Gamma\left(n+2\bar{\varepsilon}+1\right)}{\Gamma\left(2\bar{\varepsilon}+1\right)}{}_{2}F_{1}\left(-n,n+2\bar{\varepsilon}+2L;2\bar{\varepsilon}+1;s\right),\end{array}$
(31)

where

$L=\frac{1}{2}+\sqrt{\frac{1}{4}+b^{2}\beta^{2}+\gamma d_{2}}.$ (32)

The ro-vibrational energies $E_{nl}$ of the DF potential can then be derived
from the relation,

$k+(2n+1)\pi^{\prime}(s)=n^{2}.$ (33)

Finally, the ro-vibrational energy spectrum can be secured from,

$\begin{array}[]{ll}E_{nl}&=d_{4}-\frac{\alpha^{2}\hbar^{2}}{2\mu}\left(\frac{d_{3}^{2}}{(n+L)^{2}}+\frac{(n+L)^{2}}{4}\right),\end{array}$
(34)

where

$\begin{array}[]{ll}d_{3}&=\frac{1}{2}(2+b)b\beta^{2}+\frac{1}{2}\gamma(d_{2}-d_{1}),\\\
d_{4}&=D_{e}+\frac{l(l+1)\hbar^{2}\alpha^{2}}{2\mu}d_{0}+\frac{\alpha^{2}\hbar^{2}d_{3}}{2\mu}.\end{array}$
(35)

Therefore, the radial solution can be expressed as,

$R_{n}(r)=N_{nl}\,s^{\bar{\varepsilon}}(1-s)^{L}\,{}_{2}F_{1}\left(-n,n+2\bar{\varepsilon}+2L;2\bar{\varepsilon}+1;s\right),$
(36)

where

$N_{nl}=\left[\frac{2\bar{\varepsilon}\alpha(n+\bar{\varepsilon}+L)\Gamma(n+2\bar{\varepsilon}+1)\Gamma(n+2\bar{\varepsilon}+2L)}{n!(n+\bar{\varepsilon})\Gamma(n+2L)\left[\Gamma(2\bar{\varepsilon}+1)\right]^{2}}\right]^{\frac{1}{2}},$
(37)

is the normalization constant, derived from,

$\displaystyle\int|\psi_{n,l,m_{l}}({\bf r})|^{2}\,d{\bf r}=1,~{}d{\bf
r}=r^{2}\sin\theta\,dr\,d\theta\,d\phi.$ (38)

Finally, the explicit form of eigenfunctions of DF potential can be written
as,

$\psi_{n,l,m_{l}}({\bf
r})=\frac{N_{nl}}{\sqrt{2\pi}}\frac{R_{n}(r)}{r}Y_{l}^{m_{l}}(\theta)\,e^{im_{l}\phi},$
(39)

where

$Y_{l}^{m_{l}}(\theta)=\left[\frac{(2l+1)(l-|m_{l}|)!}{2(l+|m_{l}|)!}\right]^{\frac{1}{2}}\,P_{l}^{|m_{l}|}(\cos\theta)\,e^{im_{l}\phi}.$
(40)

###  2.3 Energy spectrum

To illustrate the effectiveness of our approximation, selected energies of DF
potential are tabulated for two diatomic molecules, _viz._ , H2 (homo-nuclear)
and LiH (hetero-nuclear). In each case, 10 states are considered to cover a
broad range. The parameters used are as follows: Bohr radius = 0.52917721092Å,
Hartree energy = 27.21138505 eV, electron rest mass = 5.48577990946$\times
10^{-4}$ amu. Spectroscopic parameters for H2 are: $D_{e}=4.7446$ eV,
$r_{e}=0.7416$Å, $\alpha=1.9426$Å-1, $\mu=0.50391$ amu, while for LiH, these
are: $D_{e}=2.515283695$ eV, $r_{e}=1.5956$Å, $\alpha=1.1280$Å-1,
$\mu=0.8801221$ amu. In columns 3-7 of Table 1, our estimated energies from
Eq. (34) are presented for five $\lambda$, namely $-$0.25, 0, 0.75, 1 and 1.25
respectively. A total of 10 states are given for each molecule for arbitrary
$n,l$ quantum numbers. It may be noted that, the energies were considered by
means of super-symmetric quantum mechanics [53] and path integral approach
[25]; these correspond to the case of $\lambda=0$. On the other hand, the
energies of NU method [48, 21], asymptotic iteration [55] and Hellmann-Feynman
theorem [23] represent the $\lambda=1$ case. For the so-called circular or
$l=0$ states, all the columns produce identical energies, which is expected
from Eq. (14). But for $l\neq 0$ states, the energies in five columns differ
from each other. It is noticed that the energies in 6th column (relating to
$\lambda=1$) completely reproduce those reported in [21].

Table 1:  Selected ro-vibrational energies (in eV) of H2 and LiH, represented
by DF potential. These correspond to five representative $\lambda$’s,
indicated in parentheses. For more details, see text.

$n$ | $l$ | $E_{nl}(-0.25)$ | $E_{nl}(0)$ | $E_{nl}(0.75)$ |  $E_{nl}(1)$11footnotemark: 1 | $E_{nl}(1.25)$
---|---|---|---|---|---|---
|  | H2 |
0 | 0 | 0.349980221 | 0.349980221 | 0.349980221 | 0.349980221 | 0.349980221
1 | 0 | 0.996777053 | 0.996777053 | 0.996777053 | 0.996777053 | 0.996777053
1 | 1 | 1.009941717 | 1.010018022 | 1.010246935 | 1.010323238 | 1.010399541
2 | 0 | 1.580248366 | 1.580248366 | 1.580248366 | 1.580248366 | 1.580248366
2 | 1 | 1.592275484 | 1.592360547 | 1.592615732 | 1.592700793 | 1.592785853
2 | 2 | 1.616260599 | 1.616516418 | 1.617283848 | 1.617539648 | 1.617795443
5 | 0 | 2.986148433 | 2.986148433 | 2.986148433 | 2.986148433 | 2.986148433
5 | 5 | 3.118869192 | 3.120534519 | 3.12552957 | 3.127194276 | 3.128858827
5 | 10 | 3.452477485 | 3.458819802 | 3.477833081 | 3.484166293 | 3.490497236
5 | 15 | 3.935487614 | 3.950184048 | 3.994199652 | 4.008847068 | 4.023482314
|  | LiH |
0 | 0 | 0.103334650 | 0.103334650 | 0.103334650 | 0.103334650 | 0.103334650
1 | 0 | 0.302005955 | 0.302005955 | 0.302005955 | 0.302005955 | 0.302005955
1 | 1 | 0.303739903 | 0.303759653 | 0.303818903 | 0.303838653 | 0.303858402
2 | 0 | 0.490685861 | 0.490685861 | 0.490685861 | 0.490685861 | 0.490685861
2 | 1 | 0.492346421 | 0.4923672887 | 0.492429891 | 0.492450759 | 0.492471626
2 | 2 | 0.495665829 | 0.495728464 | 0.495916365 | 0.495978997 | 0.496041629
5 | 0 | 0.999006401 | 0.999006401 | 0.999006401 | 0.999006401 | 0.999006401
5 | 5 | 1.020697858 | 1.021060955 | 1.022150163 | 1.022513206 | 1.022876235
5 | 10 | 1.077958316 | 1.079301598 | 1.083330322 | 1.084672857 | 1.086015205
5 | 15 | 1.169212076 | 1.1721857 | 1.181101026 | 1.184070957 | 1.187039967

11footnotemark: 1These energies coincide with those from [21].
---

These observations of Table 1 are graphically shown in Fig. 2, where the ro-
vibrational energies of H2 molecule are depicted for three representative
states. The left and right panels give $D_{e}-E_{nl}$ and $E_{nl}$
respectively; the state labels are provided in parentheses. Shifted energies
are plotted as well, to ease the facilitation of comparison with select
literature values. Five values of $\lambda$ ranging from $-0.5$ to 1.5
including 0, are considered. In each case, a line is generated. In the right
panel, only reference energies of [21] are available, which are marked in blue
color. For $\lambda=1$, they fall on the straight line. In left panels,
reference energies are available in scattered manner. But in any case, the
available ones related to $\lambda=1$ directly fall on the line as expected.
The GPS energies of [22] remain slightly left off $\lambda=0$, in $l\neq 0$
states in panels (B) and (F). However, this can happen as it is a purely
numerical approach. A similar plot is also found for LiH molecule in Fig. 3.
In this case, the reference results of [23] falls slightly right off
$\lambda=1$. This is presumably because of the differences in implementation
of their method from ours; we use a NU scheme while they employed a Hellmann-
Feynman approach. Note that in panel (C), the path-integral [25] and GPS [22]
energies overlap. The latter results, like H2, hover in the neighborhood of
$\lambda=0$. For a few other molecules, quite similar observations are made;
hence they are omitted here.

Figure 2:  Plot of $(D_{e}-E_{nl})$ (left) and $E_{nl}$ (right) of H2 molecule
in eV, against $\lambda$. The state indices are given inside the panels. The
colored points show literature energies from various references as: red [22],
magenta [23], green [25], yellow [55] and blue [21]. Figure 3:  Plot of
$(D_{e}-E_{nl})$ (left) and $E_{nl}$ (right) of LiH molecule in eV, against
$\lambda$. The state indices are given inside the panels. The colored points
show literature energies from various references as: red [22], magenta [23],
green [25], yellow [55] and blue [21].

##  3 Thermodynamical properties in DF potential

Now we move on to the derivation of thermodynamic quantities using the
approximate solution of Sec. 2. The starting point for such calculations is
the vibrational partition function. This can be constructed from a direct
summation over all possible vibrational energy levels available to the system,
namely,

$Q(\beta)=\sum\limits_{n=0}^{N_{max}}e^{-\beta E_{nl}}$ (41)

where $\beta=\frac{1}{kT}$, $k$ denotes the Boltzmann constant and $T$ is the
temperature. For a finite summation with an upper bound $N$, the Poisson
summation formula can be written as [60],

$\sum\limits_{n=0}^{N}f(n)=\frac{1}{2}\left[f(0)-f(N+1)\right]+\sum\limits_{m=-\infty}^{\infty}\displaystyle\int_{0}^{N+1}f(x)e^{-i2\pi
mx}dx.$ (42)

Under the lowest-order approximation $(m=0)$, the above summation formula
leads to [60],

$\sum\limits_{n=0}^{N}f(n)=\frac{1}{2}\left[f(0)-f(N+1)\right]+\displaystyle\int_{0}^{N+1}f(x)dx.$
(43)

With the help of above relation, one can rewrite Eq. (41) in following form
(‘cl’ stands for classical),

$\begin{array}[]{lr}Z=Q_{cl}(\beta)&=e^{-\beta\,d_{4}}\left\\{\frac{1}{2}e^{\frac{\alpha^{2}\hbar^{2}\beta}{2\mu}\left(\frac{d_{3}^{2}}{L^{2}}+\frac{L^{2}}{4}\right)}-\frac{1}{2}e^{\frac{\alpha^{2}\hbar^{2}\beta}{2\mu}\left(\frac{d_{3}^{2}}{(N_{\max}+L)^{2}}+\frac{(N_{\max}+L)^{2}}{4}\right)}+J^{(0)}_{0}(\beta)\right\\}\end{array}$
(44)

where

$\begin{array}[]{ll}J^{(j)}_{m}(\beta)&=\displaystyle\int_{L}^{N_{max}+L}\left(\frac{d_{3}^{2}}{x^{2}}+\frac{x^{2}}{4}\right)^{j}\exp\left[\frac{\alpha^{2}\hbar^{2}\beta}{2\mu}\left(\frac{d_{3}^{2}}{x^{2}}+\frac{x^{2}}{4}\right)-i2\pi
mx\right]\,dx\\\ N_{max}&=\left[\sqrt{2\,|d_{3}|}-L\right]+1.\end{array}$ (45)

$N_{max}$ is a (+)ve integer, $\left[x\right]$ denotes the integral part of
$x$, and $\left[x\right]\leq x$. The energy is maximum at level $n=N_{max}$.
The total partition function is a sum of classical $(m=0)$ and quantum
corrected ($Q_{QC}$) contribution corresponding to $(m\neq 0)$, i.e.,

$Q_{T}(\beta)=Z(\beta)+Q_{QC}(\beta)$ (46)

where

$Q_{QC}(\beta)=\displaystyle\exp\left[-\beta\,d_{4}\right]\sum\limits_{\begin{subarray}{c}m=-N\\\
m\neq 0\end{subarray}}^{N}\,\exp\left[i2\pi mL\right]J^{(0)}_{m}(\beta)$ (47)

Now one can proceed for evaluation of quantities like internal energy
$(U_{T}(\beta,N),U_{cl}(\beta))$, specific heat
$(C_{T}(\beta,N),C_{cl}(\beta))$, free energy $(F_{T}(\beta,N),F_{cl}(\beta))$
and entropy $(S_{T}(\beta,N),S_{cl}(\beta))$ by using the usual expressions.
The total partition function and corresponding thermal properties can be
conveniently expressed in following simplified manner,

$\begin{array}[]{lll}Q_{T}(\beta,N)&=H_{0}(\beta,N),&Z=Q_{cl}=H_{0}(\beta,0)\\\
U_{T}(\beta,N)&=-\displaystyle\frac{1}{Q_{T}}\frac{\partial
Q_{T}}{\partial\beta}=d_{4}-\frac{H_{1}}{Q_{T}},&U_{cl}(\beta)=U_{T}(\beta,0)\\\
C_{T}(\beta,N)&=-k_{\beta}\beta^{2}\frac{\partial
U}{\partial\beta}=\displaystyle
k_{\beta}\beta^{2}\left(\frac{H_{2}}{Q_{T}}-\frac{H_{1}^{2}}{Q_{T}^{2}}\right),&C_{cl}(\beta)=C_{T}(\beta,0)\\\
F_{T}(\beta,N)&=-k_{\beta}T\ln Q_{T}=-\frac{1}{\beta}\ln
Q_{T},&F_{cl}(\beta)=F_{T}(\beta,0)\\\ S_{T}(\beta,N)&=k_{\beta}\ln
Q_{T}+k_{\beta}T\frac{\partial\ln Q_{T}}{\partial T}=k_{\beta}\ln
Q_{T}+k_{\beta}T\left(\frac{H_{1}}{Q_{T}}-d_{4}\right),&S_{cl}(\beta)=S_{T}(\beta,0)\end{array}$
(48)

where

$\begin{array}[]{ll}H_{j}(\beta,N)&=\left(\frac{\alpha^{2}\hbar^{2}}{2\mu}\right)^{j}e^{-\beta\,d_{4}}\left[\frac{1}{2}\left(\frac{L^{2}}{4}+\frac{d_{3}^{2}}{L^{2}}\right)^{j}e^{\frac{\alpha^{2}\hbar^{2}\beta}{2\mu}\left(\frac{L^{2}}{4}+\frac{d_{3}^{2}}{L^{2}}\right)}-\frac{1}{2}\left(\frac{(N_{max}+L)^{2}}{4}+\frac{d_{3}^{2}}{(N_{max}+L)^{2}}\right)^{j}\right.\\\
&e^{\frac{\alpha^{2}\hbar^{2}\beta}{2\mu}\left(\frac{(N_{max}+L)^{2}}{4}+\frac{d_{3}^{2}}{(N_{max}+L)^{2}}\right)}\left.+\sum\limits_{\begin{subarray}{c}m=-N\end{subarray}}^{N}\,\exp\left[i2\pi
mL\right]J^{(j)}_{m}(\beta)\right],\\\ \frac{\partial
H_{j}(\beta,N)}{\partial\beta}&=H_{j+1}(\beta,N)-d_{4}\,H_{j}(\beta,N).\end{array}$
(49)

The classical partition function and the corresponding classical thermal
quantities are given in the right side.

Figure 4:  Plot of classical and $N$-corrected thermal properties as functions
of $\beta$, for H2 (left) and LiH (right), in atomic unit. The classical
quantities are shown in black, while the effect of quantum correction on
thermal properties are shown for $N=5$ (red) and $N=10$ (blue). Inset of each
figure shows the difference between classical and quantum corrected thermal
properties. Insets of (A), (B) represent $(Q_{cl}(\beta)-Q_{T}(\beta,5))$
(red) and $(Q_{cl}(\beta)-Q_{T}(\beta,10))$ (blue) for H2 and LiH. In a
similar manner, other quantities are depicted as follows: (C),(D): internal
energies, (E),(F): specific heats, (G),(H): free energies, and (I),(J):
entropies.

Figure 4 depicts various calculated quantities, such as $Z,U,C,F,S$, for H2
and LiH in left and right panels respectively, with respect to $\beta$, in
atomic unit. The classical partition function $Q_{cl}(\beta)$ along with its
5th-order quantum correction, $Q_{T}(\beta,5)$ and 10th order,
$Q_{T}(\beta,10)$. All these calculations are performed with energies
corresponding to $\lambda=-0.25$ in Eq. (34), while $N_{max}$ for a particular
molecule is given in Eq. (45). They all monotonically increase (decrease) with
temperature ($\beta$). The differences in $Z$ from the classical and quantum-
corrected counterparts are shown in the onset panels in red and blue colors
for 5th and 10th orders respectively. It is seen that, the two differences
tend to converge in the low-$\beta$ region and show some deviations for large
$\beta$. In other words, the $m\neq 0$ terms in Eq. (42) should be taken in to
consideration for vibrational partition function in the low temperature
region. Next, panels (C), (D) illustrate the changes in vibrational mean
energies in H2 and LiH, along with the effects of quantum correction inside.
In both cases, $U$ fall off with $\beta$ starting from a certain positive
value at $\beta\to$0, i.e., it grows monotonically with T. There is a slight
difference in the nature of plots for two molecules; for LiH it is almost
straight, but for H2, one finds a minor bent. However, like the plots in (A),
(B), the differences in this occasion also remain small; and as $\beta$
progresses, the 5th and 10th order corrections show some variations (somehow
rather more pronounced for LiH). Slight fluctuation is observed in case of
LiH. For H2, the curve shows some flattish character at low T. The dependence
of vibrational specific heat $U$ for the two molecules are now recorded in
segments (E) and (F). At first it gradually increases with $\beta$ to reach a
maximum and thereafter tends to decline. In both molecules, the quantum
corrections make significant impacts in the behavior of $C$ in low T region,
with vigorous fluctuations in the corresponding differences. The vibrational
free energies $F$ are produced in (G) and (H), which show a sharp rise at low
$\beta$ to reach a maximum value, and then approaching a constant with further
growth in $\beta$. The differences in classical and quantum corrections in
this case appear to be rather small, while the 5th and 10th order results
remain in harmony with each other. The last two panels (I) and (J) depict the
analogous plots for vibrational entropy. After attaining a sharp maximum at
lower $\beta$, for both molecules it monotonically decay. Similar plots are
given in Fig. (5) for HCl and CO in left and right sides. The vibrational
partition functions drop sharply in low $\beta$ region in contrast to the
previous figure. The differences in $Z$ from classical to $N$th corrected
quantity apparently become more significant in CO. The features in $U$ plots
remain very similar. One notices much stronger deviations from classically and
corrected $C$ in case of HCl amongst all the four molecules considered in
these two figures. Also for CO, the variation of $C$ differs from the rest of
molecules. The $F$ plots in these molecules resemble those of previous figure.
The vibrational entropies very sharply decay at low $\beta$ and thereafter
remains practically constant throughout the whole range of $\beta$ studied.

Figure 5:  Plots of $Z,U,C,F,S$ as functions of $\beta$, for HCl (left) and CO
(right), in atomic unit, as in Fig. (4).

Finally, to estimate the role of parameter $\lambda$ on the thermodynamical
quantities, Fig. (6) portrays the relevant differences in $Z,U,C,F,S$ for H2
and LiH molecules in left and right panels respectively. In this case we
employ classical partition functions for all the calculations. In each case
the deviations of corresponding quantities are calculated with respect to
$\lambda=0$. Five $\lambda$ are chosen for this, _viz._ , $-$0.25, 0.25, 0.75,
1 and 1.25. It is observed that, the differences are in general, not very
significant. One finds the impact to be small through out the whole range of
$\beta$ undertaken here.

Figure 6:  The effect of $\lambda$ on thermal properties are plotted for H2
(left) and LiH (right) in atomic unit. Following differences are given:
classical partition functions in (A)-(B), internal energies in (C)-(D),
specific heats in (E)-(F), free energies in (G)-(H), and entropies in (I)-(J).
They are defined as: $\Delta
Z(H_{2})(\lambda_{i})=\left[Z(\beta)\right]_{\lambda=0}-\left[Z(\beta)\right]_{\lambda=\lambda_{i}}$,
where $\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5}$ correspond
to $-$0.25, 0.25, 0.75, 1 and 1.25. Other quantities follow accordingly.

##  4 Conclusion

A simple useful approximation for the centrifugal term is proposed for
solution of Schrödinger equation in certain potentials through an amalgamation
of $r\to 0$ and $r\to r_{e}$ limits. This is applied to derive expressions for
eigenvalues and eigenfunctions in Deng-Fan molecular potential. The method
offers accurate results ro-vibrational energies for $l=0$ as well as $l\neq 0$
states. This is tested for four molecules (H2, LiH, HCl, CO) through a
comparison with other existing results in the literature. The effect of
$\lambda$ on energies is analyzed. These are now utilized to calculate the
vibrational partition function and thermodynamic properties such as internal
energy, specific heat, free energy and entropy in two representative diatomic
molecules (H2, LiH). The impact of quantum correction as well as $\lambda$ on
these properties is discussed by considering the quantum correction up to 10th
order, which suffices our purpose. The method can be easily extended to other
potentials (both central and non-central) including molecular potentials, such
as Manning-Rosen, Pöschl-Teller, Hülthen, etc. It would be worthwhile to
verify its applicability and performance in these cases, some of which may be
taken up in future.

##  5 Acknowledgement

AKR gratefully acknowledges financial support from MATRICS, DST-SERB, New
Delhi (sanction order: MTR/2019/000012).

## References

  * 1. Morse, P.M. Phys. Rev. 1929, 34, 57.
  * 2. Rong, Z., Kjaergaard, G., Sage, M.L. Mol. Phys. 2003, 20, 2285.
  * 3. Kratzer, A. Z. Phys. 1920, 3, 289.
  * 4. Bayrak, O., Boztosun, I., Ciftci, H. Int. J. Quant. Chem. 2007, 107, 540.
  * 5. Ikhdair, S.M., Sever, R. J. Math. Chem. 2009, 45, 1137.
  * 6. Rosen, N., Morse, P.M. Phys. Rev. 1932, 42, 210.
  * 7. Manning, M.F., Rosen, N. Phys. Rev. 1933, 44, 953.
  * 8. Chen, Z.Y., Li, M., Jia, C.S. Mod. Phys. Lett. A 2009, 24, 1863.
  * 9. Nasser, I., Abdelmonem. M.S., Abdel-Hady, A. Mol. Phys. 2013, 111, 1.
  * 10. Roy, A.K. Mod. Phys. Lett. 2014, 9, 1450042.
  * 11. Pöschl, G., Teller, E. Z. Phys. 1933, 83, 143.
  * 12. Dong, S.H., Lemus, R. Int. J. Quant. Chem. 2002, 86, 265.
  * 13. Atkas, M., Sever, R. J. Mol. Struct. THEOCHEM 2006, 710, 223.
  * 14. Linnett, J.W. Trans. Faraday Soc. 1940, 36, 1123.
  * 15. Hulthén, L. Ark. Mat. Astron. Fys. 1942, 28A, 5.
  * 16. Roy, A.K. Pramana-J. Phys. 2005, 65, 1.
  * 17. Jia, C.S., Diao, Y.F., Yi, L.Z., Chen, T. Int. J. Mod. Phys. A 2009, 24, 4519.
  * 18. Gu, X.Y., Sun, J.Q. J. Math. Phys. 2010, 51, 022106.
  * 19. Lippincott, E.R. J. Chem. Phys. 1953, 21, 2070.
  * 20. Deng, Z.H., Fan, Y.P. Shandong Univ. J. 1957, 7, 162.
  * 21. Oyewumi, K.J., Oluwadare, O.J., Sen, K.D., Babalola, O.A. J. Math. Chem. 2013, 51, 976.
  * 22. Roy, A.K. Int. J. Quant. Chem. 2014, 114, 383.
  * 23. Oluwadare, O.J., Oweyumi, K.J. Eur. Phys. J. Plus 2018, 133, 422.
  * 24. Vogt, E., Sage, D.L., Kjaergaard, H.G. Mol. Phys. 2018, 117, 1629.
  * 25. Boukabcha, H., Hachama, H., Diaf, A. Appl. Math. Comput. 2018, 321, 121.
  * 26. Tietz, T. J. Chem. Phys. 1963, 38, 3036.
  * 27. Hua, W. Phys. Rev. A 1990, 42, 2524.
  * 28. Hamzavi, M., Rajabi, A.A., Thylwe, K.E. Int. J. Quant. Chem. 2012, 112, 2701.
  * 29. Roy, A.K. J. Math. Chem. 2014, 52, 1405.
  * 30. Schiöberg, D. Mol. Phys. 1986, 59, 1123.
  * 31. Zavitsas, A.A. J. Am. Chem. Soc. 1991, 113, 4755.
  * 32. Eğrifes, H., Demirhan, D., Büyükkiliç, F. Phys. Scr. 1999, 60, 195.
  * 33. Sun, Y., He, S., Jia, C.S. Phys. Scr. 2013, 87, 025301.
  * 34. Jia, C.S., Chen, T., Yi, L.Z., Lin, R.S. J. Math. Chem. 2013, 51, 2165.
  * 35. Fakhri, H., Sadeghi, J. Mod. Phys. Lett. A 2004, 19, 615.
  * 36. Guo, J.Y., Sheng, Z.Q. Phys. Lett. A 2005, 338, 90.
  * 37. Ikhdair, S.M., Sever, R. Int. J. Mod. Phys. A 2010, 25, 3941.
  * 38. Ikhdair, S.M., Sever, R. J. Mol. Struct. THEOCHEM 2006, 806, 155.
  * 39. Oyewumi, K.J., Akinpelu, F.O., Agboola, A.D. Int. J. Theor. Phys. 2008, 47, 1039.
  * 40. Hajigeorgiou, P.G. J. Mol. Spect.THEOCHEM 2010, 263, 101.
  * 41. Mesa, A.D.S., Quesne, C., Smirnov, Y.F. J. Phys. A 1998, 31, 321.
  * 42. Codriansky, S., Cordero, P., Salamó, S. J. Phys. A 1999, 32, 6287.
  * 43. Dong, S.H., Gu, X.Y. J. Phys. Conf. Ser. 2008, 96, 012109.
  * 44. Zhang, L.H., Li, X.P., Jia, C.S. Int. J. Quant. Chem. 2011, 111, 1870.
  * 45. Zhang, L.H., Li, X.P., Jia, C.S. Phys. Scr. 2009, 80, 035003.
  * 46. Dong, S.H. Commun. Theor. Phys. 2011, 55, 969.
  * 47. Oluwadare, O.J., Oyewumi, K.J., Akoshile, C.O., Babalola, O.A. Phys. Scr. 2012, 86, 035002.
  * 48. Hamzavi, M., Ikhdair, S.M., Thylwe, K.E. J. Math. Chem. 2013, 51, 227.
  * 49. Falaye, B.J., Ikhdair, S.M., Hamzavi, M., Z. Naturforsch. 2015, 70 85\.
  * 50. Wang, P.Q., Zhang, L.H., Jia, C.S., Liu, J.Y. J. Mol. Spect. 2012, 274, 5.
  * 51. Pekeris, C.L. Phys. Rev. 1934, 45, 98.
  * 52. Badawi, M., Bessis, N., Bessis, G. J. Phys. B 1972, 5, L157.
  * 53. Mustafa, O. J. Phys. B 2015, 48, 065101.
  * 54. Greene, R.L., Aldrich, C. Phys. Rev. A 1976, 14, 2363.
  * 55. Oyewumi, K.J., Falaye, B.J., Onate, C.A., Oluwadare, O.J., Yahya, W.A. Mol. Phys. 2013, 112, 127.
  * 56. Nikiforov, A.F., Uvarov, V.B. Special Functions of Mathematical Physics (Birkhäuser, Basel) (1988).
  * 57. Berkdemir, C., Han, J. Chem. Phys. Lett. 2005, 409, 203.
  * 58. Tezcan, C., Sever, R. Int. J. Theor. Phys. 2009, 48, 337.
  * 59. Udoh, M.E., Okorie, U.S., Ngwueke, M.I., Ituen, E.E., Ikot, A.N. J. Mol. Model. 2019, 25, 170.
  * 60. Strekalov, M.L. Chem. Phys. Lett. 2007, 439, 209.


# Feeling Optimistic? Ambiguity Attitudes for Online Decision Making

Jared J. Beard1, R. Michael Butts2, and Yu Gu3 1Jared J. Beard is a Ruby
Distinguished Doctoral Fellow with the Department of Mechanical and Aerospace
engineering at<EMAIL_ADDRESS>R. Michael Butts is a Master’s student
with the Lane Department of Computer Science and Electrical Engineering at
<EMAIL_ADDRESS>Gu is faculty with the Department of Mechanical and
Aerospace engineering at<EMAIL_ADDRESS>West Virginia University,
Morgantown, WV 26506, USA

###### Abstract

As autonomous agents enter complex environments, it becomes more difficult to
adequately model the interactions between the two. Agents must therefore cope
with greater ambiguity (e.g., unknown environments, underdefined models, and
vague problem definitions). Despite the consequences of ignoring ambiguity,
tools for decision making under ambiguity are understudied. The general
approach has been to avoid ambiguity (exploit known information) using robust
methods. This work contributes ambiguity attitude graph search (AAGS),
generalizing robust methods with ambiguity attitudes—the ability to trade-off
between seeking and avoiding ambiguity in the problem. AAGS solves online
decision making problems with limited budget to learn about their environment.
To evaluate this approach AAGS is tasked with path planning in static and
dynamic environments. Results demonstrate that appropriate ambiguity attitudes
are dependent on the quality of information from the environment. In
relatively certain environments, AAGS can readily exploit information with
robust policies. Conversely, model complexity reduces the information conveyed
by individual samples; this allows the risks taken by optimistic policies to
achieve better performance.

##  I Introduction

Ambiguity is a growing problem as robots are increasingly used in real-world
environments. Robots face ambiguities in the form of incomplete information
such as unknown terrain, conflicting sensory information, as well as vague
definitions of problems and concepts. Failure to consider ambiguities in the
context of decision making under uncertainty leaves robots oblivious to
potential fault conditions. Despite a rich theoretical body of literature,
tools for decision making under ambiguity are notably understudied [1]. Among
existing solutions, robust planners are designed to recognize information gaps
and ensure safe handling of worst case scenarios [2, 3]. Such a problem
description aims to avoid ambiguity and will act conservatively. Robust
behavior may be desirable for low uncertainty systems and safety critical
systems, however, conservatism limits use of such tools in many other
applications. The work here focuses on giving robots the option to trade-off
between robust policies and more exploratory or optimistic policies when faced
with ambiguity (Fig. 1).

Figure 1: In the tunnel world environment above, the agent is given a
relatively small planning horizon, leaving the the world largely ambiguous.
The agent starts near a region with small reward, but is expected to reach a
region with better reward beyond its horizon. Robust policies avoid the
unknown areas of the map and become stuck in the local minima. To overcome
this, ambiguity attitude graph search (AAGS) employs ambiguity attitudes
(specified by the tuning parameter $\alpha$) to generalize how robots consider
model ambiguities. Using optimistic attitudes ($\alpha=1$), agents can
actively explore and achieve a higher reward than robust attitudes
($\alpha=0$).

The definitions surrounding ambiguity in literature are unfortunately used
loosely and are often times conflicting. We do not aim to precisely define
these terms here, rather to contextualize them and better communicate the
motivation of our problem. Our focus will center on two discussions: the first
being the relationship between uncertainty and ambiguity in probability
distributions and the second on ambiguity in the definition of models (Fig.
2).

Cuzzolin argues uncertainty is of two kinds: predictable and unpredictable
[4]. Predictable uncertainty (sometimes called first-order uncertainty) is
likely most familiar to us. When sampling a given probability distribution, it
is expected that each outcome occurs in proportion to their likelihood
function as the number of samples tends to infinity. Unpredictable
uncertainties (sometimes called second-order or Knightian uncertainty [4]) are
ambiguous and represent gaps in knowledge about a first-order distribution
(e.g. a coin with an unknown probability of getting heads). Here the
unpredictability arises when the likelihood of distinct outcomes cannot be
distinguished. Such a scenario has a set of probability distributions, but we
cannot define which is more “correct” due to the lack of information. It is
the authors’ opinion to avoid use of the term second-order uncertainty as this
is often obfuscated with parameter uncertainty. In this case, the ambiguity is
implicit in the parameter estimation process [5], not the randomness of the
underlying parameter distribution. In practice, this distinction may be
irrelevant, but the implication is that if the parameter distribution is
known, then there is no ambiguity, hence no information gap implied by the use
of second-order uncertainty. Some economists brand unpredictable uncertainty
and ambiguity simply as uncertainty and first-order uncertainty as risk [6,
7]. However, the robotics community tends to overload uncertainty to mean both
concepts. Compare the concept of optimism in the face of uncertainty [8] with
decision making under uncertainty [9]. In the prior, ambiguous distributions
must be uncovered. In the latter, there is a sequential decision making
problem which may have a known stochastic model. Both types of uncertainty
have the concept of entropy to measure how uncertain a quantity is; in the
case of unpredictable uncertainties, entropy measures ambiguity. This
generalization quantifies not only how much mass is distributed to each
outcome, but also to what extent mass could be distributed [4]. This is of
interest in estimation processes because relatively few samples from a low
entropy source are more likely to reflect the underlying distribution that a
high entropy one.

Model uncertainty arises when an agent must decide from a set of possibly
unknown models to describe a problem [7]. The concept of model uncertainty
goes further than unpredictable uncertainties in probability distributions; it
also questions the very structure and assumptions of the problem. For the
purpose of this work, we restrict ourselves to the ambiguity present from
estimating Markov decision process (MDP) transition probabilities with sampled
data. Without loss of generality we refer to our corresponding ambiguous mass
functions as models as is commonly done in literature [10, 11]. Lastly, we
stress that model uncertainty is not necessarily ambiguous if the models are
drawn from a known generator. While this may not appear in many problems, it
is of conceptual importance as we aim to better understand the challenges in
our field.

Figure 2: Ambiguity can be expressed as a continuum. An unambiguous concept is
known or deterministic. In the case of randomness or stochastic processes,
events may be unknown beforehand, but they occur with some known frequency. As
gaps in knowledge grow, a model or distribution becomes more ambiguous or may
be completely unknown. Estimation processes are attempts (often indirectly) to
fill such gaps. This work focuses on estimating probability distributions
representing state transition models in MDP’s. These models are represented as
the belief functions described in Sec. II.

General purpose decision making tools are usually ambiguity neutral. Upper
confidence trees (UCT) [12] is a good example. While its search process
attempts to offset ambiguities using the number of samples, UCT arrives at its
policy with the expected value of each action (an expression for which
ambiguity is undefined). As mentioned above, the selection of tools which
explicitly consider ambiguity are relatively limited. Among the most powerful
tools for decision making under ambiguity are robust solvers such as the
graph-based optimistic planning (GBOP) algorithm [8]. While robust planning
offers a deeper look into the structure of ambiguity, their conservatism
imposes a limit on the scope of problems they are intended to solve. This
motivates the need for ambiguity attitudes, wherein an agent’s behavior can be
designed to exist in the space between actively seeking out ambiguity
(optimistic policies) and avoiding it (robust policies) [10, 13]. The work
here extends online planning techniques to leverage ambiguous distributions
and integrates ambiguity attitudes to solve a broader class of problems (Fig.
1).

This paper contributes the following:

  * •

Ambiguity attitudes are introduced to online decision making through the novel
ambiguity attitude graph search (AAGS) algorithm. This is demonstrated to
generalize robust decision making and incorporates exploratory behaviors
directly into an agent’s decision making process.

  * •

A technique for evaluating sets of transition models from confidence intervals
is developed. This technique maps confidence intervals to ambiguous
probability mass functions known as belief functions.

Additionally, all code for environments is open source at the Interactive
Robotics Laboratory gym (https://github.com/wvu-irl/irl-gym) and algorithms
can be found in our decision making toolbox (https://github.com/wvu-
irl/decision-making).

The remainder of the paper is as follows. Section II contains details
regarding decision making under ambiguity: defining ambiguous MDP’s and belief
functions. Section III describes the approach used to generate belief
functions, as well as the AAGS algorithm. Section IV discusses results and
identifies limitations of the proposed approach. Lastly, in Section V, the
reader will find conclusions and future work.

##  II Ambiguity MDP

MDP’s have broad use for formulating decision problems [9, 14]. As such, they
have a substantial base of literature from which to build tools. Here we
extend MDP’s to ambiguous Markov decision processes (AMDP) and define belief
functions from the theory of evidence. Unlike the “belief” in POMDP
literature, these do not capture uncertainty in a standard mass function.
Rather they generalize the concept of a mass function to support multiple
outcomes [4].

While ambiguous MDP’s have been defined before [10, 11], we suggest a
competing formulation. Existing frameworks are best suited to solutions with
parametric distributions and problems which can be defined analytically.
Conversely, the proposed AMDP formulation better supports advances from
sampling based decision making tools. Define the AMDP as
$\mathcal{M}\coloneqq\langle\mathcal{S},\mathcal{A},\mathcal{T},R,\gamma\rangle$.
Let $s\in\mathcal{S}$ be a state in the set of states and $a\in\mathcal{A}$ be
an action in the action space. The reward $R(s,a,s^{\prime})$ is a function of
a given state transition. The temporal discount factor is $\gamma$. Lastly,
$\mathcal{T}(s^{\prime}|s,a)$ is the set the transition models, represented
here as belief functions.

To define a model with a belief function, let
$\mathcal{S}^{\prime}\subseteq\mathcal{S}$ be the space of outcomes from some
state-action pair $(s,a)$. Then distribute mass to the power set of outcomes
$2^{\mathcal{S}^{\prime}}$, such that the mass of some proposition
$m(s^{\prime}):2^{\mathcal{S}^{\prime}}\rightarrow[0,1]$ and
$\sum_{s^{\prime}}m(s^{\prime})=1$, given $s^{\prime}\in\mathcal{S}^{\prime}$
[4]. As an example let our belief function over
$\mathcal{S}^{\prime}=\\{s_{1},s_{2}\\}$ be such that
$m(\\{s_{1},s_{2}\\})=0.6$, $m(\\{s_{1}\\})=0.1$, and $m(\\{s_{2}\\})=0.3$.
For the latter two, mass acts as is typical of a standard probability. For
$m(\\{s_{1},s_{2}\\})$, this posits that, given our limited knowledge, 60% of
the time we would be unable to predict the outcome. These belief functions are
typically mapped to lower and upper probabilities rather than using the
undistributed mass directly. Let there be two sets
$\mathcal{B},\mathcal{C}\in\mathcal{S}^{\prime}$ The lower probability or
belief represents the minimum mass supporting a set of outcomes

$Bel(\mathcal{B})=\sum_{\mathcal{C}\subseteq\mathcal{B}}m(\mathcal{C}),$ (1)

while the upper probability or plausibility represents corresponding maximum

$Pl(\mathcal{B})=\sum_{\mathcal{B}\cap\mathcal{C}\neq\emptyset}m(\mathcal{C}).$
(2)

These bounds on probability are demonstrated in Figure 3 (a) where each mass
is known to $\pm\epsilon$. Figure 3 (b) captures the distribution as a belief
function where mass can now support multiple hypotheses.

With a set of transition models, it is no longer possible to use Bellman’s
equation and the expected value to solve the AMDP. The $\alpha$-Hurwicz
function is a common choice as it allows for models with varying ambiguity
attitudes [13]. More recent work sought to integrate expectation value, risk,
and Knightian uncertainty into a single metric termed the additive model [6].
The added generality is unnecessary for the work here, and $\alpha$-Hurwicz
has some the desirable properties explored at greater length in Sec. III-A.
According to $\alpha$-Hurwicz, the value is defined as follows

$V(s)=\max_{a}[(1-\alpha)\underline{\mathbb{E}}(s,a)+\alpha\overline{\mathbb{E}}(s,a)].$
(3)

Here $\alpha\in[0,1]$ is a tuning parameter representing the ambiguity
attitude. The lower expectation is

$\underline{\mathbb{E}}(s,a)=min_{T\in\mathcal{T}}\sum_{s^{\prime}}T(s^{\prime}|s,a)[(r(s,a,s^{\prime})+\gamma\underline{\mathbb{E}}(s^{\prime},a^{*})].$
(4)

The upper expectation similarly follows suit, with the $max$ function. Note
that $a^{*}$ indicates the optimal action for $s^{\prime}$.

##  III Methods

The goal of this work is to generalize robust decision making to include
ambiguity attitudes. The AAGS algorithm used to do so approximates a decision
making problem as a set of models using samples from a blackbox simulator. As
with several other sampling based methods [8, 9, 12, 14], agents supply an
$(s,a)$-pair and the model returns observations in the form of a state
transition and reward $\omega\coloneqq(s^{\prime},r)$. It is further assumed
that the underlying model is an MDP with discrete state and actions spaces.
AAGS is a graph-based solver, thus building up transition models for each
$(s,a)$-pair. Whereas prior methods, such as GBOP, assume ambiguity is
adequately described by confidence intervals, we contribute a novel way to
generate belief functions from these intervals to more accurately represent
the distribution of mass. Based on our mass functions upper and lower
expectations are computed over the models. The online decision making
algorithm AAGS then uses the specified ambiguity attitude to select the best
action.

###  III-A Computing Belief Functions

Approximate decision making algorithms are often evaluated in part by the
accuracy of estimated models. Here we concern ourselves with fixed confidence
methods [15]. Such algorithms specify a confidence interval parameterized by
accuracy $\epsilon$ and confidence $\delta$. The idea is that a relevant
quantity $z$ is approximated with $\hat{z}$ such that
$P(|z-\hat{z}|<\epsilon)\geq 1-\delta$. Whereas these bounds are generally
applied to the reward or value estimates, AAGS directly bounds the mass terms
instead, similar to the “knows what it knows” framework [16, 17]. In doing so,
AAGS adopts a distributional perspective (shown to be more expressive than
value based ones [18]). Treating the accuracy bounds as constraints, we
evaluate the belief function as a linear optimization. Subsequent discounting
of the model against the solution space allows us to distribute previously
unused mass in the confidence term, improving the precision of our model.

(a)

(b)

Figure 3: Our approach accepts multinomial distributions estimated by
collecting samples from a blackbox simulator. (a) Using confidence intervals
of size $\epsilon$ for each mass term, we constrain the model. This yields
upper $Pl$ and lower $Bel$ estimates on the mass of each outcome. (b) Using
these bounds, we distribute the mass to a belief function. Such a formulation
assigns mass to multiple hypotheses when information is ambiguous.

Begin by assuming we have a multinomial model, specified to some accuracy and
confidence $(\epsilon,\delta)$. Our approach is agnostic to how these are
obtained (implementation details in Appendix). We wish to redistribute the
$\epsilon$ mass from each outcome to second-order mass terms, generating a set
of transition models. Define $\Theta$ as the set of outcomes $(s^{\prime},r)$
from our state transition, and let $2^{\Theta}$ be the power set of all such
outcomes. Let each proposition of size one be $\omega_{i}\in\Theta$ and
propositions of size two or more be
$\theta_{j}\in\\{2^{\Theta}\setminus\Theta\\}$. Given the accuracy is a
confidence interval, let it constrain our belief functions such that the lower
and upper bounds are, respectively, $Bel(\omega)=p_{\omega}-\epsilon$ and
$Pl(\omega)=p_{\omega}+\epsilon$, $\forall\omega\in\Theta$. Let 0 and 1 bound
these values, as necessary.

Using these constraints, the problem is structured as a linear optimization
$Ax=B$. Note, for each proposition $\omega$,

$Pl(\omega)-Bel(\omega)=\sum_{\omega\in\theta}m(\theta).$ (5)

Let the columns $j$ of $A$ correspond to each $\theta_{j}$ and each row $i$ to
the singleton outcomes $\omega_{i}$. Then compose the $A$ matrix using the
indicator function $\mathds{1}(\omega_{i})\coloneqq\omega_{i}\in\theta_{j}$.
From Eq. 5, each row of $B$ is equal to $Pl(\omega_{i})-Bel(\omega_{i})$. Next
append a row of ones to the end of $A$ and $1-\sum Bel(\omega)$ to $B$ to
constrain the total mass by that available. The variable $x$ is a column of
all $\theta_{j}$, as these are the elements we wish to assign our mass. This
optimization problem is as follows:

$\begin{bmatrix}\lx@intercol\hfil\mathds{1}(\omega_{i})\hfil\lx@intercol\\\
1&1&...&1\end{bmatrix}\begin{bmatrix}2^{\Theta}\setminus\Theta\\\
\end{bmatrix}=\begin{bmatrix}Pl(\omega_{i})-Bel(\omega_{i})\\\ 1-\sum
Bel(\omega_{i})\end{bmatrix}.$ (6)

Solving for $2^{\Theta}\setminus\Theta$ yields the distribution of mass to the
ambiguous terms. In practice the number of terms in $2^{\Theta}$ can become
intractable. We overcome this by limiting $|\Theta|\leq 12$ and precomputing
$A^{-1}$. This offered a reasonable balance between accuracy and computational
constraints: storing power sets is memory intensive and online matrix
inversions incur compute expenses. Outcomes were binned when the number
exceeded this limit.

Next, we discount our model by the confidence to account for the chance that
the true model is not within the accuracy specification. Selection of the
$\alpha$-Hurwicz decision criterion has benefits here. Theory of evidence
requires us to know the entire outcome space, while in our scenario, we are
provided only sampled outcomes and bounds. The upper and lower expectations
used in $\alpha$-Hurwicz criterion permit us to neglect these terms.

###### Lemma 1

Let $U\cup L\cup\Theta$ be the entire outcome space of the decision, where $U$
is the upper bound of the outcome space and $L$ the lower bound. For a given
belief function $b$ over $U\cup L\cup\Theta$, inclusion of any other element
$v:L\leq v\leq U$ does not change the upper $\overline{\mathbb{E}}(b)$ or
lower $\underline{\mathbb{E}}(b)$ expectations.

Proof. We begin by looking at the upper expectation; the lower expectation
follows trivially. By definition of the upper expectation, the mass goes to
the element in each set having the greatest magnitude. All $v$ less than $U$
result in all mass allocated to $U$. If $v=U$, the terms are the same and we
are left with the original set $\\{U\cup L\cup\Theta\\}$ of outcomes.
$\square$

Let $m_{b}$ be the mass of elements in our belief function. The discounted
distribution is then evaluated using the approach of [19] with the confidence
as the discount parameter. Thus, the original model to $(\epsilon,\delta)$
accuracy/confidence is equivalent to the following mass function
$m(\theta)=(1-\delta)m_{b}(\theta)$ for $\theta\in 2^{\Theta}$, $m(\\{U\cup
L\cup\Theta\\})=\delta$. For simplicity, let $\theta$ now include elements of
size one. This is reflected in the Algorithm 1; here $P$ is our sampled
distribution. Notice in line two the accuracy is computed from the model
confidence and number of samples used to make the distribution $P.N$. This is
used to generate our constraints $Bel,Pl$ on the model. Notice also, that in
the for loop, we are not omitting all elements in the set of outcomes
$P.\Theta$, rather the proposition $\\{P.\Theta\\}$ itself.

$b\leftarrow\\{\\}$, $n=size(P)$

$Bel,Pl\leftarrow f(\delta,P.N)$

$b\leftarrow b\cup(1-\delta)Bel$

$m=A^{-1}B$

for  _$\theta\in\\{2^{P.\Theta}\setminus\\{P.\Theta\\}\\}$_ do

$b\leftarrow b\cup(\theta,(1-\delta)m(\theta))$

return $b\cup(\\{U\cup L\cup P.\Theta\\},\delta)$

Algorithm 1 dist2belief($P$,$\delta$,$U$,$L$)

###  III-B Ambiguity Attitude Graph Search

Our algorithm ambiguity attitude graph search (Alg. 2), generalizes robust
planners to permit ambiguity attitudes. These attitudes are specified by
$\alpha\in[0,1]$, where $\alpha=0$ is a conservative or a robust policy.
Conversely, $\alpha=1$ is an optimistic or ambiguity seeking policy. We assume
the user knows the reward bounds $R_{min}$, $R_{max}$. For an infinite horizon
MDP, the bounds on the value follow as $R_{max}/(1-\gamma)$ and
$R_{min}/(1-\gamma)$ [8]. Users also specify accuracy and confidence
requirements $\epsilon$, $\delta$.

AAGS is a graph based planner, thus using samples more efficiently than tree
methods. Assuming the simulator is much more expensive than the algorithm
itself, the cost of any added computation are far outweighed by the reduced
sample complexity. We allow for cyclic graphs with no constraints on
connectivity. In discrete static environments, users may wish to reuse the
graph, as nodes are more likely to connect to existing neighbors. Multiple
independent subgraphs may occur in dynamic environments, making it preferable
to throw out the graph after each iteration if reconnection is unlikely.

AAGS begins from state $s_{0}$ with an empty graph $\mathcal{G}$ and searches
the solution space with $N$ trajectories up to horizon $D$. For each time
step, the action having the greatest upper expectation is selected (in line
with the optimism in the face of uncertainty paradigm [8]). Both the upper and
lower expectations compute the belief functions from Alg. 1 internally. If a
state has not been sampled enough to be considered known (i.e., satisfying the
constraints $\epsilon$, $\delta$ for all actions), the generative model $G$ is
sampled. Otherwise, we borrow from the “knows what it knows” framework [16]
and sample the transition estimated model $T$ to reduce computational
overhead. New outcomes are added to the graph and the current state is added
to the list of parents. Next, our bounds are backpropagated following the
approach of Leurent, et al. [8]. Therein, upper and lower bounds are updated
through an approximate value iteration over the set of parents. Parents are
appended to the list for those states, whose bounds have changed by more than
$\beta=\dfrac{1-\gamma}{\gamma}\Delta$. Here $\Delta$ is the regret or
difference in value between the second best and optimal actions for our model
at a state. The key difference with respect to Leurent, et al. being we
compute our upper $s.U$ and lower $s.L$ value bounds at a state through belief
functions rather than confidence intervals. To conclude, the ambiguity
attitude $\alpha$ and the expectations are fed into Eq. 3 to select the best
action.

input: _$\alpha,R_{max},R_{min},\epsilon,\delta,N,D$_

initialize: _$\mathcal{G}\leftarrow\\{\\}$_

Procedure  _search( $s_{0}$)_

if  _$s_{0}\not\in\mathcal{G}$_ then

$\mathcal{G}\leftarrow\mathcal{G}\cup s_{0}$,

$n\leftarrow 0$

while  _$n <N$_ do

$s\leftarrow s_{0}$

$parents\leftarrow s$

while  _$s\neq terminal$ and $depth<D$_ do

$a\leftarrow\arg\max_{b}\overline{\mathbb{E}}(s_{0},b)$

if  _$t(s,a) <t(\epsilon,\delta)$_ then

$s^{\prime},r\sim G(s,a)$

else

$s^{\prime},r\sim T(s,a)$

$\mathcal{G}(s,a).update(s^{\prime},r)$

$parents\leftarrow parents\cup s^{\prime}$, $s\leftarrow s^{\prime}$

$backpropagate(parents)$

$return\arg\max_{a}\alpha\underline{\mathbb{E}}(s_{0},a)+(1-\alpha)\overline{\mathbb{E}}(s_{0},a)$

Procedure  _backpropagate( $P$)_

for  _$s\in parents$_ do

$parents.pop(s)$

$L\leftarrow s.L$, $U\leftarrow s.U$

$s.L\leftarrow\underline{\mathbb{E}}(s)$, $s.U\leftarrow\hat{\mathbb{E}}(s)$

if  _$s.L-L >\beta$ or $U-s.U>\beta$_ then

$parents\leftarrow parents\cup s.parents$

Algorithm 2 Ambiguity Attitude Graph Search

##  IV Results

Three environments were used to evaluate AAGS. Path planning was conducted in
a grid world (Fig. 4 (a)), sailing world (Fig. 4 (b)), and tunnel world (Fig.
1). These scenarios are selected to highlight interactions of ambiguity
attitudes with the complexity of the environment. By analyzing grid world and
sailing world, we can control the entropy of the random models to illustrate
how optimistic and robust policies can both mislead agents. Tunnel world
focuses on ambiguity beyond the planning horizon and illustrates how robust
policies can limit an agent’s potential. Environment parameters can be found
from the default parameters of the gyms at https://github.com/wvu-irl/irl-gym.

UCT [12] is used to validate the performance of AAGS. UCT was to tuned to have
a rollout horizon of 25 and its exploration parameter $c$ was set to 0.5 and 8
for grid world and sailing, respectively. We draw attention to the fact that
UCT is ambiguity neutral so this exploration parameter affects its search, not
the final policy. In a similar way, AAGS and GBOP use optimism in the face of
uncertainty to select actions during search. When it is time to select an
action, however, GBOP selects a robust policy and AAGS selects one in line
with its ambiguity attitude.

Because GBOP [8] is a robust planner, it is considering ambiguity to some
extent. However, this does not provide a suitable foundation to compare the
implications of ambiguity attitudes. To illustrate the differences of
structuring ambiguity in the value estimate rather than the transition models,
we extended GBOP to include ambiguity attitudes. Such an extension was
trivial; vanilla GBOP provides upper and lower bounds on the value estimate so
these were used with Equation 3. The new algorithm GBOP+ has vanilla GBOP as
the limiting case with $\alpha=0$.

All algorithms presented here were allocated 500 samples per times step. We
collected data for other values, however they were omitted due to similarity
of their results. Each environment limited the number of time steps to 100.
Grid world and sailing world had 325 sampled state-goal pairs for each trial
considered here.

(a)

(b)

Figure 4: Example environments for (a) grid world and (b) sailing world. In
grid world, the agent (blue circle) aims to reach the goal in green and is
increasingly rewarded (darker regions) as it gets closer to the final goal. In
sailing, the agent (blue triangle) is similarly rewarded. Though now it must
try to avoid boundaries and balance following the wind (red triangles) with
moving towards the goal.

###  IV-A For Robustness

Robust policies attempt to improve the lower bound of potential outcomes [2,
3]. As a consequence, they emphasize exploiting knowledge. Exploitation can be
beneficial in relatively low entropy environments, as results are likely to
occur again. The grid world environment is a prime example. We take as example
a 50 x 50 grid world with stochastic transitions. The robot can move to each
of its four cardinal neighbors with a 10% chance of remaining in place.

Figure 5 compares the three algorithms’ expected reward as a function of the
initial distance. Conservative AAGS ($\alpha=0$) and UCT perform comparably,
yielding better rewards than both GBOP+ attitudes and optimistic AAGS
($\alpha=1$). Paired with Figure 6, offers a better idea of what is happening.
First, we see that the reward achieved by AAGS decreases roughly in proportion
with an increase in optimism. Meanwhile, GBOP+ is largely unaffected by
changes in the ambiguity attitude.

Both results are a byproduct of the relative lack of complexity in this
environment. Every time the algorithms sample their environment, there is a
significant possibility that result will be replicated in the actual
environment. Robust policies are more inclined to follow this information,
leading directly to the goal. On the other hand, optimistic policies attempt
to explore the environment to ensure there is no better outcome. this leads to
extra efforts despite clear signals from the environment where to go. GBOP+
has much tighter bounds than AAGS; thus it can readily detect the simplicity
of the environment and converge quickly to a solution. This leaves the
ambiguity attitudes with little room to change the policy.

(a)

(b)

Figure 5: Comparison of the reward achieved by each algorithm with respect to
distance in grid world. In (a), AAGS and GBOP+ have conservative policies
($\alpha=0$), whereas in (b), they are optimistic ($\alpha=1$). Conservative
AAGS and UCT achieve the highest rewards over the greatest distances. Both
GBOP+ policies and the optimistic AAGS policy perform worse. Optimism is
penalized here by unnecessary exploration because the relatively simple
environment provides clear signals to the goal.

(a)

(b)

Figure 6: Comparison of the effect $\alpha$ on the reward earned by (a) AAGS
and (b) GBOP+ with respect to starting distance from the goal in grid world.
Due to the relative tightness of its bounds and static nature of this
environment, GBOP+ converged rather quickly. Robust AAGS policies, however,
were rewarded for their exploitation of the simple model behavior.

###  IV-B For Optimism

As we look to use optimistic attitudes, our study leads to more complex
environments. Sailing world offers complexity on two fronts. First, inclusion
of the wind directions in the state, leads to a large state space. Second, the
high entropy associated with transitions in these large state spaces means the
likelihood of experiencing a sampled state transition again is nearly 0. As
such, the exploitation of robust policies can leave an agent chasing down
fruitless outcomes. A 40 x 40 sailing world is considered, where wind is
dynamic with a 10% chance of direction change. In sailing world, the agent can
move forward or change direction by $\pm 45^{\circ}$; the agent is rewarded
for moving towards the goal and penalized for moving against the wind or near
the map borders.

Again we begin by comparing all three algorithms (Fig. 7). In this scenario,
the optimistic AAGS policies achieve the highest rewards over the greatest
distances. UCT and optimistic GBOP+ policies perform well following closely
behind, though they suffer greater losses as distance increases. All taper off
roughly linearly with the distance. At the middle of the pack is the
conservative AAGS policy. The conservative GBOP+ policy mostly avoids
penalties, but does not achieve much in the way of reaching the goal. These
trends are reflected in Figure 8. For both AAGS and GBOP+ conservative
policies becomes mired, while even moderately optimstic policies are able to
reach the goal reliably.

These results can be traced back to the increased complexity of the
environment. Due to the high entropy of a large state space and dynamic
transitions, states are rarely visited. Robust policies quite literally “go
where the winds take them”, as the saying goes. The constant fluctuations of
the wind give them conflicting signals. This leaves the agent with a noisy
lower bound causing it to get lost. Conversely, optimistic policies can
capture the consistent positive reward of moving toward the goal reflected in
the upper bound.

(a)

(b)

Figure 7: Comparison of the reward achieved by each algorithm with respect to
distance in the sailing world. In (a), AAGS and GBOP+ have conservative
policies ($\alpha=0$), whereas in (b), they are optimistic ($\alpha=1$). In
the optimistic scenario, AAGS has the highest and most consistent rewards. UCT
and the optimistic GBOP+ follow, achieving similar results to each other.
Conservative AAGS policies are met with moderate success, while conservative
GBOP+ policies can do little other than avoid penalties. This demonstrates how
optimistic policies can capture the consistent positive signals, but
conservative policies become lost in conflicting perturbations from small
penalties.

(a)

(b)

Figure 8: Comparison of the effect $\alpha$ on the reward earned by (a) AAGS
and (b) GBOP+ with respect to starting distance from the goal in the sailing
world. Even moderately optimistic policies are able to reach the goal
consistently. Meanwhile, conservative policies are dominated by the small
negative fluctuations of the wind. They are left following the wind to avoid
penalties.

###  IV-C Optimism for Exploration

Lastly, we turn to optimism as a means of exploration. Often information will
be beyond the horizon of a planner. Some works are dedicated solely to
information gathering [20], but what happens when we wish to do so and perform
other objectives? Optimistic ambiguity attitudes naturally induce exploratory
behavior, offering a simple way to integrate exploration with other mission
objectives.

For this scenario, we introduce the tunnel world (Fig. 1). In tunnel world,
the agent is situated in a narrow passage near a reward that is locally
substantial but globally undesirable. Some distance away there is a larger
reward that is the true goal. The horizons of AAGS and GBOP+ are set to 10
steps so the world is largely unknown and the optimal goal can not be seen by
the agents Given the small scale of this problem, UCT is not included as its
dynamic horizon would find the goal for the number of samples given; in a
future work this is a feature we would like to include with AAGS. As this is
illustrative we only conduct 25 trials for each ambiguity attitude and
distance (to the large reward) pair.

Figure 1 demonstrates AAGS attempting to reach the large goal. The search area
progressively increases as ambiguity attitudes become more optimistic, until
the agent can efficiently search and find the ultimate goal. Figure 9
represents this behavior as a function of the number of steps it takes to
reach the final goal. With AAGS, the level of optimism need to reach the goal
in a short time increases relatively quickly. GBOP+, on the other hand, has a
more gradual response to the increase in distance, giving a user finer control
over the exploration.

(a)

(b)

Figure 9: Controlling the distance and ambiguity attitudes of (a) AAGS and (b)
GBOP+, the average number of steps required to reach the high reward goal are
shown. Notice that AAGS responds rather abruptly to increases in distance.
GBOP+, conversely has a more gradual response offering more control over the
amount of exploration.

###  IV-D General Discussion

As demonstrated in the results, AAGS can solve path planning tasks with
comparable performance to UCT with respect to accumulated reward and distance.
More importantly, these results show the impact ambiguity attitudes can have
on performance. Given that robust methods are already in use it is important
to understand their limitations. Similarly, the generalization of robustness
to ambiguity attitudes, motivates a need to understand optimistic behaviors as
well. This is accomplished by contextualizing the result with respect to the
complexity of environments and their reward structure.

###  IV-E Limitations

There are a few limitations the reader should be aware of. First, whereas
recent approaches have sought to tie accuracy to estimates of the value and
the structure of the algorithm, we have here implemented a generic mapping to
estimate accuracy of a probability mass function (Appendix). Thus, the bounds
on distributions in AAGS converge relatively slowly, leaving greater ambiguity
in the model. AAGS, however, is designed to more precisely distribute
ambiguity in the models and factors confidence into the models. While it
requires more thorough study, this helps compensate for the slower
convergence. This also introduces the question of whether methods can be
designed to reduce ambiguity by better understanding how information gaps are
distributed.

We performed analysis in discrete environments. Connectivity of graphs can be
tricky in continuous domains, so this planner would behave more like a tree
search without a method to abstract states. Sailing world approximates this
with its large state space, however this does not appear to have negatively
impacted performance. Lastly, as noted before, we artificially limit the
number of outcomes to compute our bounds. This is due to storage reasons as
the the $A$ matrix scales with complexity $\mathcal{O}(2^{\Theta})$. In the
future, we would like to develop an error constrained method for distributing
the mass.

##  V Conclusions

The increasing complexity of tasks we assign to robots entails increasing
ambiguity. Robustness is already a consideration when developing autonomous
systems. However, as demonstrated here, it is not always the best avenue when
ambiguity is present. This work supports a more complete conceptualization of
ambiguities in the decision making process through the introduction of
ambiguity attitudes. Whereas robust attitudes avoid ambiguity and exploit
their knowledge, optimistic attitudes seek ambiguity and are naturally
exploratory. In simple domains with clear signals or safety critical
applications, robust methods can be used to safely and efficiently reach an
objective. However, as signals become more muddied (e.g., in the sailing
world), taking risks can have a significant payout. Ambiguity attitudes
demonstrates the flexibility to make this trade-off. The results for this work
were supported by the ambiguity attitude graph search (AAGS) algorithm
developed here. AAGS was shown to perform comparably to existing techniques in
both general decision making under uncertainty (via UCT) and in robust
decision making (via GBOP).

Looking to the future there are a few key directions this work should be
taken. Among the most important are extensions to continuous domains and with
partially observable MDP’s. Some work has begun in this direction [21], though
they do not consider the implications of ambiguity attitudes. Furthermore, as
ambiguity increases, the model assumptions come into questions, not just their
probability distributions [5, 7]. Lastly, as a more exploratory measure, we
consider the option of making ambiguity attitudes a function of the state.
Agents could then assess for themselves the danger or payoff associated with a
more or less optimistic attitude.

## References

  * [1] A. J. Keith and D. K. Ahner, “A survey of decision making and optimization under uncertainty,” _Annals of Operations Research_ , vol. 300, no. 2, pp. 319–353, 2021.
  * [2] G. N. Iyengar, “Robust dynamic programming,” _Mathematics of Operations Research_ , vol. 30, no. 2, pp. 257–280, 2005.
  * [3] W. Wiesemann, D. Kuhn, and B. Rustem, “Robust markov decision processes,” _Mathematics of Operations Research_ , vol. 38, no. 1, pp. 153–183, 2013\.
  * [4] F. Cuzzolin, _The Geometry of Uncertainty: The Geometry of Imprecise Probabilities_. Springer Nature, 2020.
  * [5] A. H. Briggs, M. C. Weinstein, E. A. Fenwick, J. Karnon, M. J. Sculpher, and A. D. Paltiel, “Model parameter estimation and uncertainty analysis: a report of the ispor-smdm modeling good research practices task force working group–6,” _Medical decision making_ , vol. 32, no. 5, pp. 722–732, 2012\.
  * [6] Y. He, J. S. Dyer, J. C. Butler, and J. Jia, “An additive model of decision making under risk and ambiguity,” _Journal of Mathematical Economics_ , vol. 85, pp. 78–92, 2019.
  * [7] M. Marinacci, “Model uncertainty,” _Journal of the European Economic Association_ , vol. 13, no. 6, pp. 1022–1100, 2015.
  * [8] E. Leurent and O.-A. Maillard, “Monte-carlo graph search: the value of merging similar states,” in _Asian Conference on Machine Learning_. PMLR, 2020, pp. 577–592.
  * [9] M. J. Kochenderfer, _Decision making under uncertainty: theory and application_. MIT press, 2015.
  * [10] S. Saghafian, “Ambiguous partially observable markov decision processes: Structural results and applications,” _Journal of Economic Theory_ , vol. 178, pp. 1–35, 2018.
  * [11] N. Bäuerle and U. Rieder, “Markov decision processes under ambiguity,” _arXiv preprint arXiv:1907.02347_ , 2019.
  * [12] L. Kocsis and C. Szepesvári, “Bandit based monte-carlo planning,” in _European conference on machine learning_. Springer, 2006, pp. 282–293.
  * [13] K. J. Arrow and L. Hurwicz, “An optimality criterion for decision-making under ignorance,” _Uncertainty and expectations in economics_ , vol. 1, 1972.
  * [14] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [15] A. Jonsson, E. Kaufmann, P. Ménard, O. Darwiche Domingues, E. Leurent, and M. Valko, “Planning in markov decision processes with gap-dependent sample complexity,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 1253–1263, 2020.
  * [16] L. Li, M. L. Littman, T. J. Walsh, and A. L. Strehl, “Knows what it knows: a framework for self-aware learning,” _Machine learning_ , vol. 82, no. 3, pp. 399–443, 2011.
  * [17] S. M. Kakade, _On the sample complexity of reinforcement learning_. University of London, University College London (United Kingdom), 2003.
  * [18] M. G. Bellemare, W. Dabney, and R. Munos, “A distributional perspective on reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 449–458.
  * [19] G. Shafer, _A mathematical theory of evidence_. Princeton university press, 1976, vol. 42.
  * [20] J. J. Beard, _Environment Search Planning Subject to High Robot Localization Uncertainty_. West Virginia University, 2020.
  * [21] S. Ross, B. Chaib-draa, and J. Pineau, “Bayes-adaptive pomdps,” _Advances in neural information processing systems_ , vol. 20, 2007.

## APPENDIX

In this paper we use confidence intervals in a few different ways. First, we
wish to compute the accuracy given a confidence and the number of samples. We
also use the accuracy and confidence requirements to specify a stopping rule
for resampling the model. As such, it is convenient to have a relation that
affords this flexibility. To arrive at one, we performed Monte Carlo trials to
sample multinomial distributions and estimate the likelihood a model would
meet a given accuracy. Following, we arrived at this relation

$t(\epsilon,\delta)=\dfrac{ln\dfrac{1}{1.25(1-\delta)-1/6}}{\epsilon(-ln\dfrac{1}{1.5(1-\epsilon)+1/3})^{2}}.$
(7)

As intuition would suggest, the following constraints are also satisfied.
First, $\delta(\epsilon=1,t)=0$: no requirement on accuracy implies it is met
trivially. Next, $\epsilon(\delta=0,t)=0$: no confidence in a model implies it
is known to perfect precision. We found this approximation performed well and
better matched the distribution than existing (algorithm agnostic) bounds [16,
17], while keeping $\epsilon,\delta\in[0,1]$. It is an open question to better
approximate the relation between these values. Furthermore, algorithm specific
relations may improve performance by exploiting structural components of their
solution. Given the nonlinearity in $\epsilon$, we used gradient descent to
find the solution.