In order to distinguish policies that prescribe
good from bad actions in transient states,
we need to evaluate the so-called bias of a policy from transient states.
However, we observe that most (if not all) works in
approximate discounting-free policy evaluation thus far are developed
for estimating the bias solely from recurrent states.
We therefore propose a system of approximators for the bias
(specifically, its relative value) from transient and recurrent states.
Its key ingredient is a seminorm LSTD (least-squares temporal difference),
for which we derive its minimizer expression that enables approximation by
sampling required in model-free reinforcement learning.
This seminorm LSTD also facilitates the formulation of a general unifying procedure
for LSTD-based policy value approximators.
Experimental results validate the effectiveness of our proposed method.

§ INTRODUCTION

Consider an environment where there are two types of states:
those that are visited infinitely many times by an agent, and
those that are not (even though the agent is modelled to operate up to infinity).
The members of the former group are called recurrent states, whereas
those of the latter are called transient states.
If all recurrent states form a single closed irreducible set, then we have
the so-called unichain Markov chain (MC).
It is closed in that once the agent is in any member of the set,
the agent cannot go outside to any non-member state.
It is irreducible because from any member of the set, the agent can
visit any other member.
In reinforcement learning (RL), such a unichain MC is induced by (at least)
one of the stationary policies of the Markov decision process (MDP) model,
for which we call such a model a unichain MDP.

The work in this  is concerned with evaluating a stationary policy $\pi$
from both transient and recurrent states in terms of
discounting-free evaluation functions.
Particularly, the policy value function of interest is the bias (denoted by $v_b$)
of $\pi$ as follows,

\begin{equation}
v_b(\pi, s)
\eqdef \lim_{t_{\mathrm{max}} \to \infty}
\E{A_t \sim \pi(\cdot|s_t), S_{t+1} \sim p(\cdot| s_t, a_t)}{
\sum_{t=0}^{t_{\mathrm{max}} - 1}
\Big( r(S_t, A_t) - v_g(\pi) \Big) \Big| S_0 = s, \pi},
\quad \forall s \in \setname{S},
\label{equ:bias_def}
\end{equation}

where $S_t$ and $A_t$ are discrete state and action random variables on
the state set $\setname{S}$ and action set $\setname{A}$
of an infinite-horizon MDP with one-step state transition distribution $p$,
and reward function $r$.
Here, $v_g$ denotes the gain (the average-reward) value function, which is
state-invariant whenever the induced MC is unichain.
Both bias and gain do not involve any discount factor (hence, they are said to
be discounting-free; the discounted reward value function).

Evaluating the bias from both state types is essential for carrying out
further policy selection on gain-optimal policies that induce unichain MCs.
This is because the gain value function only concerns with the long-run rewards
(which are earned in recurrent states).
It ignores rewards earned at the outset in transient states, in which
gain-optimal policies therefore cannot distinguish “good” from “bad” actions.
In other words, they are suboptimal in transient states with respect to
the finest (the most selective) optimality criterion, the Blackwell optimality.

Despite the aforementioned importance, we observe that most works in RL are
designed for approximately evaluating a policy from recurrent states,
specifically for MDPs whose all induced MCs have only recurrent states.
The stationary state distribution is used for weighting the state-wise error terms
in the error function.
This applies to both discounted and discounting-free policy evaluation,
<cit.>, <cit.>.
There are a few works that estimate the policy values from transient states.
However, they are applicable merely for MCs with a single recurrent state
whose reward is zero and known (thus no estimation is needed).
For instance, <cit.> proposed
a discounted-reward estimator for environments with
multiple transient states and a single 0-reward absorbing terminal state.
The error terms are weighted by the visitation probabilities of transient states
from the initial time until absorption, which is known to happen at
the last timestep of a trial (as the agent reaches the absorbing terminal state).

In this , we propose techniques that approximate the bias value
(of any stationary policy) from multiple transient and recurrent states
in model-free RL.
This requires value approximation from both state types,
instead of either one out of two types as the above-mentioned existing works.
Moreover, the state classification is unknown since
the agent does not know the state transition distribution
and does not attempt to estimate it.
This also implies that the agent does not know when it is absorbed into
the closed, irreducible recurrent state set (the absorption time).
If the state classification was known, the state set could be sliced
and two approximators could be built:
one for the recurrent states and one for the transient states;
taking advantage of the existing works.
However in that case, some additional work would still be needed for two reasons.
First is that all recurrent states cannot simply be separated from the state set
(for the sake of the transient state value approximator)
because there must be transitions from some transient states to recurrent states,
which may have non-zero rewards affecting transient state values.
Second is because those two individual bias approximators have different offsets
from the true bias (due to the nature of the error function that they minimize).
Therefore, their approximation results need to be calibrated before being used
simultaneously in a formula that involves bias values of
both transient and recurrent states.
More exposition about these two issues are provided later
in sec:seminorm_lstdsec:approxbias.

We present a system of approximators for the bias (relative) value.
Each approximator is based on least-squares temporal difference (LSTD).
In one extreme (where the computation cost is put aside), the system
instantiates stepwise LSTD approximators that use stepwise state distributions,
denoted as $p_\pi^t$, to weight the state-wise error terms.
Here, $p_\pi^t(s) \eqdef \E{S_0 \sim \isd}{\prob{S_t = s| s_0, \pi}},
\forall s \in \setname{S}$, which indicates the probability of visiting a state $s$
in $t$ timesteps when the agent begins at an initial state $S_0 \sim \isd$
then follows a policy $\pi$.
Such a system dismisses the need for state classification, but
poses at least three challenges, for which we contribute some solutions.

First, the stepwise state distribution $p_\pi^t$ may not have
the whole state set as its support.
For example, when the agent can only begin in transient states,
$p_\pi^{t=0}$ has zero probability for any recurrent state.
The same goes for transient states with respect to $p_\pi^t$ after
the stationary state distribution is reached
(as there is no chance of visiting transient states
once the recurrent class is entered).
Consequently, the diagonal matrix derived from $p_\pi^t$ may be
positive semidefinite (PSD).
This necessitates an LSTD approximator that involves a seminorm,
hence a generalized pseudoinverse.
The main difficulty comes from the fact that the reverse product law
does not apply to pseudoinverse.
We derive the optimal (minimizing) parameter of the seminorm LSTD, and
its sampling-based estimator in sec:seminorm_lstd.

Second, a system of stepwise approximators
requires at least one parameter per timestep.
This implies an infinite number of parameters whenever there is
an infinite number of timesteps (as in infinite-horizon MDPs).
To be practical, we propose a procedure that accommodates the specification of
the desired number of approximators (hence, the desired number of parameters in a system).
It determines a number of timestep neighborhoods according to
sampling-based estimated distances among stepwise state distributions.
For every neighborhood, it then applies a seminorm LSTD that weights
the state-wise error terms using the average state distributions on that neighborhood.
This procedure is general and unifies existing LSTD-based methods,
as explained in sec:unify_lstd_based.

Third, the approximators of the proposed system estimate the bias
only up to some offset; generally one unique offset for each approximator.
This is inherently due to the limitation of the bias error-function
that those approximators individually minimize.
In order to use the resulting approximations (of all approximators)
in a formula that jointly involves transient and recurrent state values,
we need to calibrate those offsets so that all approximations have
one common offset with respect to the true bias, which is time-invariant.
In sec:approxbias, we describe such offset calibration along with
the pseudocode of the proposed system of seminorm LSTD approximators
for the bias (its relative value).

We provide experimental results in nbwpval:xprmtresult, which is
preceded by experimental setup in nbwpval:xprmtsetup.
Finally, we conclude this work in nbwpval:conclude, where
we also describe its limitations as well as some avenues for future research.
The next nbwpval:prelim presents some necessary prerequisites
for this work.

§ PRELIMINARIES

Given a stationary policy $\pi$, we are interested in computing
its bias value $v_b(\pi, s), \forall s \in \setname{S}$
as in (<ref>).
Since the policy and the value function type are fixed,
we often simplify their notations and write $v(s) \eqdef v_b(\pi, s)$.
This $v(s)$ is then called a state value of $s$.

A parametric state-value approximator is parameterized by a parameter vector
$\vecb{w} \in \setname{W} = \real{\dim(\vecb{w})}$.
Linear parameterization (which is the focus of this work) gives

\begin{equation} \label{equ:vhat_linear}
\hat{v}(s, \vecb{w}) \eqdef \vecb{w}^\intercal \cdot \vecb{f}(s) \approx v(s),
\quad \forall s \in \setname{S}, \forall \vecb{w} \in \setname{W},
\quad \text{equivalently,} \quad
\hat{\vecb{v}}(\vecb{w}) = \mat{F} \vecb{w},
\end{equation}

where $v(s)$ denotes the true (ground-truth) state value of $s$,
$\vecb{f}(s) \in \real{\dim(\vecb{w})}$ the state feature vector of $s$,
and $\mat{F} \in \real{\setsize{S} \times \dim(\vecb{w})}$
the corresponding state feature matrix
(whose $s$-th row contains $\vecb{f}^\intercal(s)$).
Here, the approximate state value vector $\hat{\vecb{v}} \in \real{\setsize{S}}$
is obtained by stacking all scalar approximations
$\hat{v}(s) \in \real{} , \forall s \in \setname{S}$ on top of each other.

One way to learn $\vecb{w}$ is by minimizing
the weighted mean squared projected Bellman error (MSPBE),
denoted as $\epb(\vecb{w})$ in (<ref>).
This error function is derived based on the identity in
the average-reward Bellman equation, namely

\begin{equation}
v(s_t) = r(s_t) - g + \E{p(\cdot|s_t)}{v(S_{t+1})},
\quad \forall s_t \in \setname{S},\ \text{equivalently,} \quad
\vecb{v} = \vecb{r} - \vecb{g} + \mat{P}\vecb{v} \eqdefr \bo[\vecb{v}],
\label{equ:poisson_avgrew}
\end{equation}

and some projection to obtain the representation of $\bo[\vecb{v}]$
in the parameter space.
Here, the reward function $r(s_t) = \E{\pi}{r(s_t, A_t)}$ corresponds to
the reward vector $\vecb{r} \in \real{\setsize{S}}$,
the gain $g \eqdef v_g(\pi)$ corresponds to the gain vector
$\vecb{g} \eqdef g \vecb{1} \in \real{\setsize{S}}$
(using the vector $\vecb{1}$, whose entries are all 1's), and
$\mat{P} \in \real{\setsize{S} \times \setsize{S}}$ is
the one-step state transition stochastic matrix of an induced MC, whose
$s_t$-th row represents a next-state conditional distribution
$p_\pi(S_{t+1}|s_t) = \sum_{a_t \in \setname{A}} \pi(a_t|s_t) p(S_{t+1}| s_t, a_t)$.
The operator $\bo: \real{\setsize{S}} \mapsto \real{\setsize{S}}$
is termed as the Bellman policy-evaluation operator on
$\vecb{v} \in \real{\setsize{S}}$.

The MSPBE is defined as follows,

\begin{equation}
\epb(\vecb{w})
\eqdef \norm{
\underbrace{
\hat{\vecb{v}}(\vecb{w}) - \tilde{\po} \bo \hat{\vecb{v}}(\vecb{w})
}_{ \Delta_{\hat{\vecb{v}}} }
= \Delta_{\hat{\vecb{v}}}^{\!\!\intercal}\ \diagptilde\ \Delta_{\hat{\vecb{v}}}
= \sum_{s \in \setname{S}} \tilde{p}(s)\ \Delta_{\hat{v}[s]}^{\!2},
\label{equ:epb}
\end{equation}

where $\tilde{\vecb{p}} \in \real{\setsize{S}}$ is a vector of probability values
of some state distribution $\tilde{p}(s) = \prob{S = s}, \forall s \in \setname{S}$,
and $\diagptilde$ is an $\setsize{S}$-by-$\setsize{S}$ diagonal matrix with
$\tilde{\vecb{p}}$ along its diagonal.
Here, $\tilde{\po}$ denotes a projection operator such that

\begin{equation} \label{equ:epb_projector}
\tilde{\po} \vecb{v} = \mat{F} \vecb{w}^\diamond,\ \text{where}\
\vecb{w}^\diamond = \argmin_{\vecb{w} \in \setname{W}}
\Big\{ \norm{\mat{F} \vecb{w} - \vecb{v}}_{\tilde{\vecb{p}}}^2
= \sum_{s \in \setname{S}} \tilde{p}(s)\
( \vecb{f}^{\intercal}(s)\ \vecb{w} - v(s) )^{2} \Big\}.
\end{equation}

At this stage, what is left to fully define $\epb$ is
the state distribution $\tilde{p}$, whose probability values serve as weights
in (<ref>) and (<ref>).
For recurrent MCs, one natural choice for $\tilde{p}$ is
the stationary state distribution $p^\star$ that indicates the state visitation
frequency in the long-run.
More precisely,

\begin{equation} \label{equ:pstar_lim}
p^\star(s) = \mathbb{E}_{S_0 \sim \isd}
\Big[ p^\star(s|s_0)
\eqdef \lim_{\tmax \to \infty} \frac{1}{\tmax} \sum_{t=0}^{\tmax - 1} p^t(s | s_0)
= \underbrace{
\lim_{\tmax \to \infty} p^{\tmax}(s | s_0)
}_\text{when the MC is aperiodic}
\Big],
\quad \forall s \in \setname{S},
\end{equation}

where $p^\star(s|s_0)$ is the limiting distribution of the stepwise $p^t(s|s_0)$
as $t$ goes to infinity
(nonetheless, $p^\star$ may be achieved in finite time).
Since all states of a recurrent MC are recurrent, its $p^\star$ has
the whole state set as its support, $p^\star(s) > 0, \forall s \in \setname{S}$.
This is advantageous because $\diagpstar$ is positive definite (PD) so that
(<ref>) and (<ref>) involve
a (weighted Euclidean) norm, and the inverse $\diagpstar^{-1}$ exists.

The state feature matrix $\mat{F}$ has a full column rank.
This is equivalent to saying that all state feature vectors are linearly independent.
(Remark: this assumption is not required by our proposed seminorm LSTD
in sec:seminorm_lstd.)

In fact, setting $\tilde{\vecb{p}} \gets \vecb{p}^\star$
leads to the LSTD method for recurrent MDPs [31].
Whenever assume:indep_fea is satisfied, the projection operator
in (<ref>) is defined as
$\po \eqdef
\mat{F} (\mat{F}^\intercal \diagpstar \mat{F})^{-1}
\mat{F}^\intercal \diagpstar$.
Then, the optimal parameter value (which minimizes $\epb$) is given by

\begin{equation}
\vecb{w}^* = \mat{X}^{-1} \vecb{y},
\qquad \text{where $\vecb{w}^* \in \setname{W}$ and}
\label{equ:td_fixedpoint}
\end{equation}

\begin{align}
\mat{X}
& = \sum_{s \in \setname{S}} p^\star(s) \sum_{s' \in \setname{S}} p(s'| s)
\Big[ \vecb{f}(s)  \Big( \vecb{f}(s) - \vecb{f}(s') \Big)^{\!\!\intercal} \Big]
= \sum_{s \in \setname{S}} p^\star(s) \vecb{f}(s)
\Big[ \Big(
\vecb{f}(s) - \sum_{s' \in \setname{S}} p(s'| s) \vecb{f}(s')
\Big)^{\!\!\intercal} \Big] \notag \\
& = \mat{F}^\intercal \diagpstar (\mat{I} - \mat{P}) \mat{F}
\qquad \in \real{\dim(\vecb{w}) \times \dim(\vecb{w})},\
\text{with an identity matrix
$\mat{I} \in \real{\setsize{S} \times \setsize{S}}$, and}
\label{equ:x_sampling} \\
\vecb{y}
& = \sum_{s \in \setname{S}} p^\star(s) \sum_{a \in \setname{A}} \pi(a|s)
\Big[ (r(s, a) - g) \vecb{f}(s) \Big]
= \sum_{s \in \setname{S}} p^\star(s) \Big[ (r(s) - g) \vecb{f}(s) \Big]
\notag \\
& = \mat{F}^\intercal \diagpstar (\vecb{r} - \vecb{g})
\qquad \in \real{\dim(\vecb{w})}. \label{equ:y_sampling}
\end{align}

The minimizer $\vecb{w}^*$ in (<ref>) involves $\mat{X}^{-1}$,
which exists whenever assume:indep_fea is satisfied, and
the singularity of $(\mat{I} - \mat{P})$ is remedied.
For example, by introducing an eligibility factor $\lambda \in (0, 1)$ such that
$\mat{X}$ then involves $(\mat{I} - \mat{P}^{(\lambda)})$, where
\eqdef (1 - \lambda) \sum_{\tau=0}^\infty \lambda^\tau \mat{P}^{\tau+1}$.
Another technique is replacing $\mat{X}$ altogether with its non-singular approximation
by some perturbation <cit.>.
Note that $(\mat{I} - \mat{P})$ is not invertible <cit.>.

An LSTD-based method approximately computes the minimizer $\vecb{w}^*$
in (<ref>) by the sample means of $\mat{X}$ and $\vecb{y}$
according to (<ref>) and (<ref>), respectively.
That is,

\begin{equation}
\hat{\vecb{w}}^*
= \hat{\mat{X}}^{-1} \hat{\vecb{y}}
= \Big( \frac{1}{\nsample} \sum_{i=1}^{\nsample}
\vecb{f}(s_i) ( \vecb{f}(s_i) - \vecb{f}(s_i'))^{\intercal} \Big)^{-1}
\Big( \frac{1}{\nsample} \sum_{i=1}^{\nsample} (r_i - g) \vecb{f}(s_i) \Big),
\label{equ:norm_lstd_sampling}
\end{equation}

where $\nsample$ denotes the number of
state $s_i$, next state $s_i'$, and reward $r_i$ samples,
which are collected by the agent through interaction with its environment.
Typically in practice, $\hat{\mat{X}} \gets \hat{\mat{X}} + \epsilon \mat{I}$ for
some small positive $\epsilon > 0$ in order to ensure
the approximation matrix $\hat{\mat{X}}$ is invertible.

One interesting property of LSTD based on $p^\star$ is
that its minimizer (<ref>) is also the solution of
the semi-gradient TD method for recurrent MDPs [28].
This method minimizes the weighted mean squared error (MSE) as follows,

\begin{equation}
\ems(\vecb{w})
\eqdef \sum_{s \in \setname{S}} \tilde{p}(s) [ v(s) - \hat{v}(s; \vecb{w}) ]^2
= \norm{\vecb{v} - \hat{\vecb{v}}(\vecb{w})}_{\tilde{\vecb{p}}}^2
= [\vecb{v} - \hat{\vecb{v}}(\vecb{w})]^\intercal\ \diagptilde\
[\vecb{v} - \hat{\vecb{v}}(\vecb{w})],
\label{equ:ems}
\end{equation}

where $\vecb{v}$ denotes the true (ground-truth) value, while
$\tilde{\vecb{p}} \gets \vecb{p}^\star$ and $\diagptilde$ are
the state distribution and the corresponding diagonal matrix, respectively.
The semi-gradient TD method follows the stochastic gradient descent (SGD) for
updating its parameter $\vecb{w}$.
For linear pameterization (<ref>) such that
$\nabla \hat{v}(S; \vecb{w}) = \vecb{f}(S)$, the SGD update rule for
an approximation iterate $\hat{\vecb{w}}^* \approx \vecb{w}^*$
is given by

\begin{equation}
\underbrace{
\hat{\vecb{w}}^* \gets  \hat{\vecb{w}}^* - \alpha \hat{\nabla} \ems(\hat{\vecb{w}}^*),
}_\text{\cf the LSTD estimator in \eqref{equ:norm_lstd_sampling}}
\quad \text{with}\
\hat{\nabla} \ems(\hat{\vecb{w}}^*)
\eqdef - \Big(
\underbrace{
r(s) - g + \hat{v}(s'; \hat{\vecb{w}}^*)
}_\text{$\approx v(s)$ based on \eqref{equ:poisson_avgrew}}
- \hat{v}(s; \hat{\vecb{w}}^*) \Big) \vecb{f}(s),
\label{equ:semigrad_update}
\end{equation}

where $\alpha$ is some positive learning rate and $\hat{\nabla} \ems$ is
the stochastic estimate of the gradient of $\ems$ (<ref>) by
one current state $s$ and one next state $s'$ sampled from
$p^\star$ and $p(\cdot|s)$, respectively.
Since the true $v(s)$ is unknown in RL, an approximation is substituted for it
in (<ref>) only after taking the gradient
(hence, the term semi-gradient[
In contrast, LSTD methods are based on (true) gradients of
the error function $\epb$.
This is possible since $\nabla \epb$ does not involve the true value $v(s)$,
see (<ref>).
Such approximation is based on the Bellman equation (<ref>).
It can be shown that the approximation iterate
$\hat{\vecb{w}}^*$ (<ref>) converges to the LSTD's minimizer
$\vecb{w}^*$ in (<ref>),
for which $\vecb{w}^*$ is called the TD fixed point <cit.>.
Note that a semi-gradient TD method needs the specification of
the learning rate $\alpha$ and the initial value for $\hat{\vecb{w}}^*$.

§ SEMINORM LSTD APPROXIMATORS

In this section, we present an LSTD approximator that minimizes $\epb$
(<ref>), whose state distribution $\tilde{p}$ does not necessarily
have the whole state set as its support,
$\tilde{p}(s) \ge 0, \forall s \in \setname{S}$.
This $\tilde{p}$ induces a positive semidefinite (PSD) diagonal matrix $\diagptilde$,
hence a seminorm $\epb$.
Consequently, minimizing $\epb$ and deriving its projector $\tilde{\po}$
(<ref>) require solving $\diagptilde$-seminorm LS problems.
We call the corresponding state-value approximator based on such $\epb$
as a seminorm LSTD.

A seminorm LSTD is useful for unichain MCs with multiple transient and
recurrent states (and with certain reward structures).
For example, since each state type has different timing
(transient states are visited at the outset before absorption, whereas
recurrent states in the long-run),
a proper $\tilde{p}$ is different for each type so that
the support of a proper type-specific $\tilde{p}$ only contains a subset of
the state set; inducing a seminorm $\epb$.
It is proper in that it provides reasonable weighting for the state-wise error terms
in (<ref>) and (<ref>), and
that it enables state sampling in (<ref>) and (<ref>).
More importantly, a seminorm LSTD facilitates the derivation of
a general approximation procedure (see sec:unify_lstd_based).

The main result of this Section is a sampling-enabler expression for the minimizer
of $\epb$ of a seminorm LSTD.
It is presented in thm:epb_minimizer (sec:samplingenabler_minimizer).
For that, the preceding sec:seminorm_part contains
the projection operator for the seminorm $\epb$ and two necessary lemmas for the minimizer.

§ A GENERAL PROCEDURE FOR LSTD-BASED POLICY EVALUATION

Equipped with seminorm LSTD (sec:seminorm_lstd), we are now ready to
devise a general unifying procedure for LSTD-based policy evaluation,
which leads to a system of LSTD approximators,
as illustrated in fig:unify_diagram.
The proposed procedure is formally presented in def:general_lstd,
for which we need the definitions of its main components as follows.

A timestep neighborhood, denoted as $\setname{N}$, is an ordered set of
consecutive timesteps from an anchor timestep $t$ to $(t + \setsize{N} - 1)$.
Every neighborhood $\setname{N}$ has a unique anchor $t$
(the earliest timestep in $\setname{N}$).
Hence, the notation $\setname{N}_t$ denotes a neighboorhod anchored at $t$.
The non-anchor member of $\setname{N}_t$, if any, is called a neighbor.
Hence, every anchor has $(\setsizesubt{N} - 1)$ neighbors.

A state probability distribution of a neighborhood $\setname{N}_t$,
denoted as $\bar{p}_t$, is a lumpsum of stepwise state probabilities $p^\tau$
from $\tau=t$ to $(t + \setsize{N} - 1)$.
That is,

\begin{equation}
\bar{p}_t(s) \eqdef
\frac{1}{\setsizesubt{N}} \sum_{\tau = t}^{t + \setsizesubt{N} - 1} p^\tau(s),
\quad \text{with}\ p^\tau(s) = \E{S_0 \sim \isd}{p^\tau(s| s_0)},
\qquad \forall s \in \setname{S},
\label{equ:p_lumpsum}
\end{equation}

where $p^\tau(s| s_0) \eqdef \prob{S_\tau = s| S_0 = s_0}$, which indicates
the probability of visiting the state $s$ in $\tau$ timesteps from
an initial state $s_0$.
This $p^\tau(s| s_0)$ is equivalent to the $[s_0, s]$-entry of $\mat{P}^\tau$,
which is the one-step transition matrix $\mat{P}$ raised to the power of $\tau$.
That is, $p^\tau(s| s_0) = \vecb{e}_{s_0}^\intercal \mat{P}^\tau \vecb{e}_s$,
where $\vecb{e}_i \in \real{\setsize{S}}$ denotes the $i$-th standard basis vector.
The $s_0$-th row of $\mat{P}^\tau$ therefore contains the probability values
of the stepwise conditional state distribution $p^\tau(\cdot| s_0)$.

The support of a neighborhood $\setname{N}_t$ is defined as
the support of its state distribution $\bar{p}_t$,
denoted as $\setname{S}(\bar{p}_t)$.
That is,
\eqdef \{ s : \bar{p}_t(s) > 0, \forall s \in \setname{S} \}
= \bigcup_{\tau=t}^{t + \setsizesubt{N} - 1} \setname{S}(p^\tau)
\subseteq \setname{S}$, where
$\setname{S}(p^\tau)$ is the support of a stepwise state distribution $p^\tau$.
Note that $\setname{S}(\bar{p}_t)$ may be a proper subset of $\setname{S}$.

A general procedure for LSTD-based policy evaluation has three steps as follows.

* Specify a number of timestep neighborhoods (def:neighborhood)
over the whole time-horizon.

* Train a seminorm LSTD approximator (sec:samplingenabler_minimizer)
for every neighborhood $\setname{N}_t$.
This approximator minimizes $\epb$ that is based on
the neighborhood state distribution $\tilde{p} \gets \bar{p}_t$ (def:pbar)
and a desired type of policy value functions.
It is termed as a seminorm LSTD-$\bar{p}_t$.

* Predict the state values at timestep $\tau$ using
the approximator of a neighborhood $\setname{N}_t$ where $\tau$ belongs
(either as an anchor or a neighbor member of $\setname{N}_t$).

The number of timestep neighborhoods (equivalently, the number of anchors or approximators)
is denoted as $\nanchor$.
This procedure forms a system of $\nanchor$ seminorm LSTD-$\bar{p}_t$
(as linear approximators).

This general procedure unifies two existing approaches to
approximate policy evaluation (sec:unify_existing),
as summarized in tab:unify_existing.
It also gives a spectrum of benefits by controlling the number of neighborhoods
More importantly, it enables value approximation for both
transient and recurrent states in unichain MDPs, which is
the main motivation for this work and is presented in sec:approxbias.

§ APPROXIMATING THE BIAS OF TRANSIENT AND RECURRENT STATES

In this section, we describe our proposed approach to approximating
the bias values of unichain MDPs in model-free RL settings.
It is devised from the general procedure (sec:unify_lstd_based) with
two additional components specific to bias computation.
They are about reference states and offsets, presented in sec:sref_offset.
Subsequently, we explain our proposed pseudocode in sec:approxbias_pcode.
Its entry point is alg:training_nbwpval, which is about
training (learning) the estimator in model-free RL settings.

§ EXPERIMENTAL SETUP

In this section, we describe the setup of our experiments, whose
results are presented in nbwpval:xprmtresult.
We begin with the environment specifications in sec:env_nbwpval,
followed by state features and state kernels (sec:state_fea_ker).
Then, various experiment schemes are described in sec:xprmt_scheme.
Lastly, we explain the evaluation metrics and protocols
in nbwpval:xprmt_metric.

§ EXPERIMENTAL RESULTS

In this section, we present the experimental results, whose setup is
described in the previous nbwpval:xprmtsetup.
There are two groups of results, namely non-sampling and sampling experiments,
as explained in sec:xprmt_scheme.
Each is evaluated with two error metrics, namely $\epbtot$ and $\emstot$

§ CONCLUSIONS, LIMITATIONS, AND FUTURE WORKS

We propose a system of seminorm LSTD approximators for estimating
the relative bias value from multiple transient and multiple recurrent states
in unichain MDPs.
To this end, we derive an expression for the minimizer of the seminorm LSTD
that enables approximation through sampling (as for model-free RL).
We also devise a general procedure for LSTD-based policy evaluation,
from which the relative value approximator emerges as a special case;
so do the other existing LSTD-based approximators for recurrent MDPs and
unichain MDPs with one 0-reward recurrent state.
Experimental results validate that a system with more approximators yields
the lower projected Bellman errors.
It is also empirically shown that timestep-neighborhoods can reasonably
be specified based on estimating the squared MMD among stepwise state distributions
(using their state samples).

The proposed method addresses the problem of minimizing MSBPE that is based on
the one-step Bellman operator $\bo$ in (<ref>).
It is known that $\bo$ can be extended to its multi-step variant.
That is,

\begin{equation*}
\vecb{v} = \bo^m [\vecb{v}]
\eqdef \left\{ \sum_{\tau=0}^{m - 1} \mat{P}^\tau (\vecb{r} - \vecb{g}) \right\}
+ \mat{P}^m \vecb{v},
\quad \text{for $m \ge 1$ (where $m=1$ gives the one-step variant)}.
\end{equation*}

Another known extension is to utilize a weighted average over all $\bo^m$
for $m=1, 2, \ldots$, which gives

\begin{equation*}
\vecb{v} = \bo_{\!\lambda} [\vecb{v}],
\quad \text{where}\
\bo_{\!\lambda} \eqdef (1 - \lambda) \sum_{m=1}^{\infty} \lambda^{m-1} \bo^m,\
\text{for a trace-decay factor $\lambda \in [0, 1)$}.
\end{equation*}

These two extensions essentially provide a way to control the bias-variance trade-off
in the value approximation.[
One way to see this bias-variance trade-off is from
the semi-gradient TD viewpoint as follows.
Setting $m$ to $1$ leads to the one-step TD algorithm.
As in (<ref>), it approximates the true value
$v(s_t) \approx r_{t+1} - g + \hat{v}(s_{t+1})$, which has low variance
(as it involves only one-step next state and reward samples)
but is biased towards the estimator $\hat{v}$.
On the other hand, $m$-step TD with $m$ approaches infinity approximates
$v(s_t) \approx (r_{t+1} - g) + (r_{t+2} - g) + \ldots$, which is
unbiased (due to no involvement of $\hat{v}$)
but has high variance (due to an infinitely long sequence of reward samples).
Note that in practice, the gain $g$ should also be approximated.
They lead to $m$-step TD, TD($\lambda$) and LSTD($\lambda$) algorithms
[25, 2].
It is interesting therefore to extend our proposed method to a system of
multiple $m$-step seminorm LSTD (or seminorm LSTD($\lambda$)) approximators.
This includes examination about how the size of a neighborhood affects
the suitable values for $m$ and $\lambda$.

This work has not taken the advantage of iterative techniques for calculating
the Moore-Penrose pseudoinverse (for computing the minimizer of the seminorm LSTD).
It also has not exploited the fact that $(1, 3)$-pseudoinverse is sufficient
for an LS solution that is not necessarily a minimum-norm solution
Note that the pseudoinverse used throughout this work,
$\mat{M}^\dagger$ for a matrix $\mat{M}$, is the full $(1, 2, 3, 4)$-pseudoinverse.

Additionally, the (finite) sample complexity of the proposed method deserves
a careful study.
This includes the relationship between the number of neighborhoods and
the number of samples to achieve a certain level of errors.
We anticipate an intricate interplay because a fewer number of neighborhoods
(equivalently, more neighbors per neighborhood) yields more violation to
the i.i.d sample condition in computing the sample means for
the minimizing parameter of the seminorm LSTD (thm:epb_minimizer).

Lastly, our proposed method can be modified to become
a system of semi-gradient seminorm TD approximators.
This is worth investigating because semi-gradient algorithms involve
an initial parameter value, which is paradoxically beneficial for
policy iteration RL methods in that it can be set to the parameter of
the last policy's value approximator.
The modification should leverage the fact that
both semi-gradient TD and LSTD algorithms converge to the same TD fixed point

[1]
Ben-Israel, A. and Greville, T. (2003).
Generalized Inverses: Theory and Applications.
CMS Books in Mathematics. Springer.

[2]
Boyan, J. A. (2002).
Technical update: Least-squares temporal difference learning.
Machine Learning, 49.

[3]
Bradtke, S. J. and Barto, A. G. (1996).
Linear least-squares algorithms for temporal difference learning.
Machine Learning, 22(1–3).

[4]
Campbell, S. L. and Meyer, C. D. (2009).
Generalized Inverses of Linear Transformations.
Society for Industrial and Applied Mathematics.

[5]
Castro, P. S. (2020).
Scalable methods for computing state similarity in deterministic
Markov decision processes.
Proceedings of the AAAI Conference on Artificial Intelligence.

[6]
Dann, C., Neumann, G., and Peters, J. (2014).
Policy evaluation with temporal differences: A survey and comparison.
Journal of Machine Learning Research, 15(24).

[7]
Ferns, N., Castro, P. S., Precup, D., and Panangaden, P. (2006).
Methods for computing state similarity in Markov decision
In Conference in Uncertainty in Artificial Intelligence.

[8]
Flamary, R., Courty, N., Gramfort, A., Alaya, M. Z., Boisbunon, A., Chambon,
S., Chapel, L., Corenflos, A., Fatras, K., Fournier, N., Gautheron, L.,
Gayraud, N. T., Janati, H., Rakotomamonjy, A., Redko, I., Rolet, A., Schutz,
A., Seguy, V., Sutherland, D. J., Tavenard, R., Tong, A., and Vayer, T.
POT: Python optimal transport.
Journal of Machine Learning Research, 22(78).

[9]
Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A.
A kernel two-sample test.
Journal of Machine Learning Research, 13(25).

[10]
Grinstead, C. and Snell, J. (2012).
Introduction to Probability.
American Mathematical Society.

[11]
Grünewälder, S., Lever, G., Baldassarre, L., Pontil, M., and Gretton,
A. (2012).
Modelling transition dynamics in MDPs with RKHS embeddings.
In Proceedings of the 29th International Conference on Machine

[12]
Haasdonk, B. and Bahlmann, C. (2004).
Learning with distance substitution kernels.
In Pattern Recognition, volume 3175. Springer Berlin

[13]
Hartwig, R. (1986).
The reverse order law revisited.
Linear Algebra and its Applications, 76.

[14]
Harville, D. A. (1997).
Matrix Algebra From a Statistician’s Perspective.

[15]
Hordijk, A., Dekker, R., and Kallenberg, L. C. M. (1985).
Sensitivity analysis in discounted Markovian decision problems.
Operations Research Spektrum, 7.

[16]
Lan, C. L., Bellemare, M. G., and Castro, P. S. (2021).
Metrics and continuity in reinforcement learning.
In Proceedings of the AAAI Conference on Artificial

[17]
Lewis, T. O. and Newman, T. G. (1968).
Pseudoinverses of positive semidefinite matrices.
SIAM Journal on Applied Mathematics, 16(4).

[18]
Liu, R. and Olshevsky, A. (2021).
Temporal difference learning as gradient splitting.
In Proceedings of the 38th International Conference on Machine

[19]
Muandet, K., Fukumizu, K., Sriperumbudur, B., and Schölkopf, B. (2017).
Kernel mean embedding of distributions: A review and beyond.
Foundations and Trends® in Machine Learning, 10.

[20]
Proszynski, W. and Sosnowski, A. (1995).
Matrix expressions for minimum seminorm least squares solutions of
linear systems and the related inverse.
Demonstratio Mathematica, 28(3).

[21]
Puterman, M. L. (1994).
Markov Decision Processes: Discrete Stochastic Dynamic
John Wiley & Sons, Inc., 1st edition.

[22]
Song, T., Li, D., Cao, L., and Hirasawa, K. (2016).
Kernel-based least squares temporal difference with gradient
IEEE Transactions on Neural Networks and Learning Systems,

[23]
Strens, M. J. A. (2000).
A Bayesian framework for reinforcement learning.
In Proceedings of the 27th International Conference on Machine

[24]
Sugiyama, M., Suzuki, T., and Kanamori, T. (2012).
Density Ratio Estimation in Machine Learning.
Cambridge University Press.

[25]
Sutton, R. S. (1988).
Learning to predict by the methods of temporal differences.
Machine Learning, 3(1).

[26]
Sutton, R. S. and Barto, A. G. (2018).
Introduction to Reinforcement Learning.
MIT Press.

[27]
Tian, Y. (2019).
Reverse order laws for generalized inverses of products of two or
three matrices with applications.

[28]
Tsitsiklis, J. N. and Roy, B. V. (1999).
Average cost temporal-difference learning.
Automatica, 35(11).

[29]
Ueno, T., Kawanabe, M., Mori, T., Maeda, S.-i., and Ishii, S. (2008).
A semiparametric statistical approach to model-free policy
In Proceedings of the 25th International Conference on Machine

[30]
Xu, X., Xie, T., Hu, D., and Lu, X. (2005).
Kernel least-squares temporal difference learning.
Information Technology - IT, 11.

[31]
Yu, H. and Bertsekas, D. P. (2009).
Convergence results for some temporal difference methods based on
least squares.
IEEE Transactions on Automatic Control, 54(7).


# Generalization of Legendre functions applied to Rosen-Morse scattering
states

F. L. Freitas<EMAIL_ADDRESS>
###### Abstract

A generalization of associated Legendre functions is proposed and used to
describe the scattering states of the Rosen-Morse potential. The functions are
then given explicit formulas in terms of the hypergeometric function, their
asymptotic behavior is examined and shown to match the requirements for states
in the regions of total and partial reflection. Elementary expressions are
given for reflection and transmission coefficients, and an integral identity
for the generalized Legendre functions is proven, allowing the calculation of
the spectral measure of the induced integral transform for the scattering
states. These methods provide a complete classical solution to the potential,
without need of path integral techniques.

##  I Introduction

It is well knownNatanzon (1979) that the hypergeometric function provides
closed form solutions to the Schrödinger equation for a large number of
potentials. Nevertheless, a full classification of these potentials and their
solutions remains difficult, in part because the orthogonality of their states
cannot be described in terms of classical Sturm-Liouville theory, since their
differential equations possess singularities.

Perhaps the most famous instance of this theoretical phenomenon is in the
radial differential equation for the hydrogen atom. See, for example, exercise
13.2.11 of Arfken and WeberArfken and Weber (2005), where although Laguerre
polynomials appear in the potential wave functions, the orthogonality
relations they obey are unrelated to the one actually derived in the
construction of said polynomials.

For this reason, generaliations of the classical theory are often proposed,
leading to new insights about the mathematics behind these analytically
solvable systems. One example is the rediscovery of Romanovski polynomials and
their finite orthogonalityWeber (2007) in the bound states of a hypergeometric
potentialCompean and Kirchbach (2005).

Of all these solvable potentials, one of the most useful and insteresting is
the one proposed by Rosen and MorseRosen and Morse (1932), given by the
expression

$V(x)=-\frac{\hbar^{2}\alpha(\alpha+1)}{2m\Delta^{2}}\operatorname{sech}^{2}\left(\frac{x}{\Delta}\right)+\frac{\hbar^{2}\beta}{m\Delta^{2}}\tanh\left(\frac{x}{\Delta}\right).$
(1)

Although initially proposed to study the properties of molecules, the
potential soon drew attention to its mathematical properties. Not only it
describes a well trapping a finite number of bound states, but it’s also like
a smooth step potential, possessing a zone where particles are totally
reflected, and another where they have a chance of being transmitted with
changes in frequency.

The Rosen-Morse potential can be studied in the context of path
integralsKleinert and Mustapic (1992); Lin (1997); Grosche (1989), where its
propagator is derived by relating the problem to motion in a manifold and
applying transformations. It is also a generalization of the modified Pöschl-
Teller potentialFlügge (2012), the Woods-Saxon potentialWoods and Saxon (1954)
and the Morse-Feshbach potentialMorse and Feshbach (1953).

It has been noticed recentlyFreitas (2018) that the operator linking adjacent
states in the Rosen-Morse potential is nonlocal, and needs tools from
fractional calculus for its construction. This is in stark contrast to its
symmetric form, where a differential operator does exist and is obtainable
from properties of associated Legendre functionsDong and Lemus (2002),
suggesting that the orthogonality of the Rosen-Morse eigenstates is in a sense
deeper than the modified Pöschl-Teller and that of the Coulomb potential.

In this work, we propose that the most natural description of these states is
not in terms of the plain hypergeometric function, but of closely related
special functions. These functions are the solution of a differential equation
very similar to that of Legendre functions, and elegantly describe both the
scattering and bound states of the Rosen-Morse potential.

Since the normalization constants of the bound states have already been
derived previouslyFreitas (2018), we focus mostly on scattering states,
calculating reflection and transmission coefficients by analyzing the
asymptotic behavior of the functions, and deriving the normalization of the
continuum spectrum with a generalization of an integral identity for
associated Legendre functions reported by Magnus, Oberhettinger and SoniMagnus
_et al._ (1966) and applied to quantum mechanical problems by GroscheGrosche
(1990). This normalization is then used to evaluate the spectral measure of
the integral transform induced by these functions, and we arrive at much
simpler expressions than those obtained with path integral methodsGrosche and
Steiner (1998).

This paper is organized as follows. In Section II the generalized Legendre
functions are defined, shown to be solutions to the Schrödinger equation for
the Rosen-Morse potential, and expressed in terms of the hypergeometric
function. In Section III their asymptotic behavior is elucidated, and it’s
shown that they can oscillate at different frequencies or even decay or grow
exponentially at different rates at $x\to+\infty$ and at $x\to-\infty$. In
section IV, we exploit our asymptotic analysis to map the parameters of the
generalized functions to the energy of Rosen-Morse scattering states, while
discussing briefly how one could also use them to find bound states. In
section V, we calculate reflection and transmission coefficients for states
with energy above the potential barrier, using several identities to simplify
the final expression. In section VI, we prove a useful integral identity for
the generalized Legendre functions. In section VII, the identity is used to
evaluate the continuum normalization of the scattering states and the spectral
measure for their integral transform. Finally, in section VIII, the
conclusions are outlined.

##  II Generalized Legendre functions

Before we start, let us briefly state the conventions for usual associated
Legendre functions, taken from Gradshteyn and RyzhikGradshteyn and Ryzhik
(2014). An associated Legendre function is a particular solution to the
differential equation

$(1-z^{2})\frac{d^{2}u}{dz^{2}}-2z\frac{du}{dz}+\left[\nu(\nu+1)-\frac{\mu^{2}}{1-z^{2}}\right]u=0,$
(2)

where $\mu$ and $\nu$ are understood as arbitrary complex parameters.

These functions are defined in the complex plane with a branch cut in the
interval $(-\infty,1)$. In our case, though, we are interested in solutions in
the interval $(-1,1)$, and these are constructed by averaging the previously
mentioned complex function above and below the cut, which also solve (2).
These functions _on the cut_ are the ones we are interested in generalizing.

With this in place, let us consider functions $D^{\mu,\eta}_{\nu}(x)$, which
satisfy a more general equation than (2), given by

$(1-x^{2})\frac{d^{2}D^{\mu,\eta}_{\nu}}{dx^{2}}-2x\frac{dD^{\mu,\eta}_{\nu}}{dx}+\left[\nu(\nu+1)-\frac{\mu^{2}+2\mu\eta
x+\eta^{2}}{1-x^{2}}\right]D^{\mu,\eta}_{\nu}=0.$ (3)

The motivation for considering this equation becomes clearer when we perform
the change of variables $x=\tanh(v)$. Writing $\psi=D^{\mu,\eta}_{\nu}$, we
can do the transformation

$\frac{d\psi}{dx}=\cosh^{2}v\frac{d\psi}{dv}$ (4)

$\frac{d^{2}\psi}{dx^{2}}=\cosh^{2}v\frac{d}{dv}\left(\cosh^{2}v\frac{d\psi}{dv}\right)=\cosh^{4}v\frac{d^{2}\psi}{dv^{2}}+2\sinh
v\cosh^{3}v\frac{d\psi}{dv}$ (5)

Substituting on (3), we obtain

$\cosh^{2}v\frac{d^{2}\psi}{dv^{2}}+\left[\nu(\nu+1)-2\mu\eta\sinh v\cosh
v-(\mu^{2}+\eta^{2})\cosh^{2}v\right]\psi=0.$ (6)

Another way of writing this differential equation is

$\psi^{\prime\prime}+[\nu(\nu+1)\operatorname{sech}^{2}v-2\mu\eta\tanh
v-\mu^{2}-\eta^{2}]\psi=0,$ (7)

and comparing it with the adimensional form of (1), obtained with
$\hbar=2m=\Delta=1$:

$\psi^{\prime\prime}+[\alpha(\alpha+1)\operatorname{sech}^{2}x-2\beta\tanh
x+E]\psi=0,$ (8)

it is readily seen that they are equivalent by means of the substitutions
$E=-\mu^{2}-\eta^{2}$, $\mu\eta=\beta$ and $\nu=\alpha$.

Equation (7) can be transformed into the hypergeometric differential equation.
We will do this by closely following the method used in a previous workFreitas
(2018), but there are important differences, since we want to reduce to
associated Legendre functions when $\eta=0$. The basic idea is to employ the
_ansatz_

$\psi(v)=e^{-av}\cosh^{-b}vF(v).$ (9)

The equation changes toFreitas (2018):

$F^{\prime\prime}-2(a+b\tanh
v)F^{\prime}+\left\\{\left[\nu(\nu+1)-b(b+1)\right]\operatorname{sech}^{2}v+(2ab-2\mu\eta)\tanh
v+a^{2}+b^{2}-\mu^{2}-\eta^{2}\right\\}F=0.$ (10)

The behavior at $v\to+\infty$ is

$F^{\prime\prime}-2(a+b)F^{\prime}+\left[(a+b)^{2}-2\mu\eta-\mu^{2}-\eta^{2}\right]F=0,$
(11)

and at $v\to-\infty$:

$F^{\prime\prime}-2(a-b)F^{\prime}+\left[(a-b)^{2}+2\mu\eta-\mu^{2}-\eta^{2}\right]F=0.$
(12)

Parameters $\mu$ and $\eta$ are chosen so that the factors multiplying $F$ at
(11) and (12) vanish. In order to obtain the correct Legendre functions at the
$\eta=0$ limit, we must use

$a=-\mu,\qquad b=-\eta.$ (13)

Using these values for $a$ and $b$ we get many cancellations on (10), and
changing to $u=(1+\tanh v)/2$, we findFreitas (2018):

$u(1-u)\frac{d^{2}F}{du^{2}}+[\mu-\eta+1-2(1-\eta)u]\frac{dF}{du}+[\nu(\nu+1)+\eta(1-\eta)]F=0.$
(14)

This is the hypergeometric equation, typically written asAbramowitz and Stegun
(1965)

$u(1-u)F^{\prime\prime}+[t-(r+s+1)u]F^{\prime}-rsF=0.$ (15)

For our problem, the parameters are

$r=-\eta-\nu,\qquad s=-\eta+\nu+1,\qquad t=-\eta+\mu+1.$ (16)

Equation (15) is a second order differential equation, and we’re interested in
solutions in the interval $[0,1]$. Its set of solutions forms a two-
dimensional vector space, and there are several ways to express the general
form. Our strategy will be to pick a solution of (15) and define our
generalized Legendre function $D^{\mu,\eta}_{\nu}(x)$ in terms of it, so that
when $\eta=0$ we reduce to the usual associated Legendre function. It turns
out that, to achieve this goal, the solution of (15) we need to choose
isAbramowitz and Stegun (1965):

$F(r,s;r+s+1-t;1-u)=F\left(-\nu-\eta,\nu+1-\eta;1-\mu-\eta;(1-x)/2\right)$
(17)

In view of the hyperbolic relations

$x=\tanh v=\frac{e^{2v}-1}{e^{2v}+1}\Rightarrow e^{2v}(x-1)=-1-x\Rightarrow
e^{2v}=\frac{1+x}{1-x},$ (18)

and

$\cosh v=\sqrt{\frac{1}{1-\tanh^{2}v}}=\sqrt{\frac{1}{1-x^{2}}},$ (19)

we can finally write the definition for our generalized Legendre function:

$D^{\mu,\eta}_{\nu}(x)=(1-x^{2})^{-\frac{\eta}{2}}\left(\frac{1+x}{1-x}\right)^{\frac{\mu}{2}}F\left(-\nu-\eta,\nu+1-\eta;1-\mu-\eta;(1-x)/2\right).$
(20)

If we compare this with the definition of usual associated Legendre
functionsGradshteyn and Ryzhik (2014):

$P^{\mu}_{\nu}(x)=\frac{1}{\Gamma(1-\mu)}\left(\frac{1+x}{1-x}\right)^{\frac{\mu}{2}}F(-\nu,\nu+1;1-\mu;(1-x)/2)$
(21)

Notice that, when $\beta=0$, (20) reduces to

$D^{\mu,0}_{\nu}(x)=\Gamma(1-\mu)P^{\mu}_{\nu}(x).$ (22)

##  III Asymptotic analysis

Since our wavefunctions will have the form $D^{\mu,\eta}_{\nu}(\tanh x)$,
let’s analyze the asymptotics of this function in $x\to+\infty$.

When we send $x\to\tanh x$ in (20), we get

$\frac{1+\tanh x}{1-\tanh x}=e^{2x},\qquad
1-\tanh^{2}x=\operatorname{sech}^{2}x=\cosh^{-2}x.$ (23)

Since $\cosh x\to e^{x}/2$ at $x\to+\infty$ and the hypergeometric factor goes
to $1$ for an argument of $0$, we obtain

$\lim_{x\to+\infty}D^{\mu,\eta}_{\nu}(\tanh x)=2^{-\eta}e^{(\mu+\eta)x}.$ (24)

Because of this, $\mu+\eta$ cannot have a positive real part, or else we don’t
get a normalizable solution.

To find the asymptotics at $x\to-\infty$, use the hypergeometric linear
transformation formulaAbramowitz and Stegun (1965)

$\displaystyle F(a,b;c;z)=$
$\displaystyle\frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}F(a,b;a+b-c+1;1-z)$
(25)
$\displaystyle+(1-z)^{c-a-b}\frac{\Gamma(c)\Gamma(a+b-c)}{\Gamma(a)\Gamma(b)}F(c-a,c-b;c-a-b+1;1-z).$

For our generalized Legendre function:

$\displaystyle D^{\mu,\eta}_{\nu}(x)=$
$\displaystyle(1-x^{2})^{-\frac{\eta}{2}}\left(\frac{1+x}{1-x}\right)^{\frac{\mu}{2}}\bigg{[}\frac{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}{\Gamma(1-\mu+\nu)\Gamma(-\mu-\nu)}F(-\nu-\eta,\nu+1-\eta;1+\mu-\eta;(1+x)/2)$
(26)
$\displaystyle+\left(\frac{1+x}{2}\right)^{-\mu+\eta}\frac{\Gamma(1-\mu-\eta)\Gamma(\mu-\eta)}{\Gamma(-\nu-\eta)\Gamma(\nu+1-\eta)}F(1-\mu-\nu,\nu-\mu;1-\mu+\eta;(1+x)/2)\bigg{]}.$

Notice that, when we send $x\to\tanh x$ in (26), in the limit $x\to-\infty$,
the hypergeometric factors go to $1$, and the factor $(1+\tanh x)/2\to
e^{2x}$, because we have the limit

$\lim_{x\to-\infty}\frac{1+\tanh x}{e^{2x}}=2.$ (27)

With $\cosh x\to e^{-x}/2$ at $x\to-\infty$, we have

$\lim_{x\to-\infty}D^{\mu,\eta}_{\nu}(x)=2^{-\eta}\left[\frac{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}{\Gamma(1-\mu+\nu)\Gamma(-\mu-\nu)}e^{(\mu-\eta)x}+\frac{\Gamma(1-\mu-\eta)\Gamma(\mu-\eta)}{\Gamma(-\nu-\eta)\Gamma(\nu+1-\eta)}e^{-(\mu-\eta)x}\right].$
(28)

##  IV Potential eigenstates

The goal of this section is to make formulas (24) and (28) more physical by
relating the parameters in the generalized Legendre function with those of the
Rosen-Morse potential. As stated before, we must perform the following
substitution on (8):

$\mu^{2}+\eta^{2}=-E,\qquad\mu\eta=\beta.$ (29)

If we are looking for bound states, then we must have $E<-2\beta$. In this
case, we can write $E=-2\beta-b^{2}$. The parameters $\mu$ and $\eta$ are both
real and are given by

$\mu=\frac{b-\sqrt{4\beta+b^{2}}}{2},\qquad\eta=\frac{-b-\sqrt{4\beta+b^{2}}}{2}.$
(30)

Notice that we must choose $\mu$ and $\eta$ to be both negative, or else the
asymptotic behavior at $+\infty$ from (24) would not lead to a normalizable
state. Besides that, the limit at $-\infty$ from (28) has two contributions
from $e^{(\mu-\eta)x}$ and $e^{-(\mu-\eta)x}$. Since in our case $\mu-\eta>0$,
the $e^{-(\mu-\eta)x}$ term needs to vanish, and we can only do this by
placing the arguments of the gamma functions at the poles at nonpositive
integers.

From this, the discrete spectrum can be obtained and the wavefunctions,
together with their normalization constants, obtained from properties of
Jacobi polynomials. Details can be found in previous workFreitas (2018).
Currently, we are interested in the continuum spectrum, so we discuss
scattering states next.

There are two scattering regimes, the first with $-2\beta<E<2\beta$, where we
scatter below the barrier and get total reflection, and the regime with
$E>2\beta$, where we scatter above the barrier and can have nonzero
transmission.

For the first regime, let’s set $E=-2\beta+k^{2}$. The parameters assume the
values

$\mu=-\frac{\sqrt{4\beta-k^{2}}}{2}+\frac{ik}{2},\qquad\eta=-\frac{\sqrt{4\beta-k^{2}}}{2}-\frac{ik}{2}.$
(31)

We see that both $\mu$ and $\eta$ have a real and an imaginary part, and are
complex conjugates of each other. There are two degenerate states for this
energy, one for $k$ and the other for $-k$, and they will contain the
information regarding the phase shift caused by the reflection depending on
the polarization of the incoming wave. Since $\mu+\eta<0$, all states decay to
zero at $x\to+\infty$, as expected. The qualitative behavior of solutions in
this regime is plotted in Figure 1.

In the second regime, above the barrier but also with $E=-2\beta+k^{2}$, the
square root in (31) gives us an imaginary number, and we have

$\mu=\frac{i\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}}{2}+\frac{ik}{2},\qquad\eta=\frac{i\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}}{2}-\frac{ik}{2}.$
(32)

Now, both $\mu$ and $\eta$ are purely imaginary. At $x\to+\infty$, our
solution oscillates with a single mode
$e^{i\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}}$. At $x\to-\infty$, there are
two modes, $e^{ikx}$ and $e^{-ikx}$. The qualitative behavior is plotted in
Figure 2. Notice that our choice of the sign of the square root is dictated by
the sign of $k$, as emphasized by the $\operatorname{sgn}$ function. This is
important to guarantee that the states corresponding to $k$ and $-k$ are
orthogonal.

Figure 1: Rosen-Morse scattering state $D^{\mu,\eta}_{\nu}(\tanh x)$ below the
energy barrier. The solution is a distorted wave close to the well, decays
into the classically forbidden region as $x\to+\infty$ and shows a wave
pattern typical of a wave interfering with its reflection as $x\to-\infty$.

Figure 2: Rosen-Morse scattering state $D^{\mu,\eta}_{\nu}(\tanh x)$ above the
energy barrier. The solution is a distorted wave close to the well. The
frequency at $x>0$ is lower than that at $x<0$. A large portion of the wave is
transmitted, so the amplitude increases at $x>0$ to conserve probability
current at lower frequency.

##  V Reflection and transmission coefficients

With the asymptotic analysis done, we are ready to evaluate reflection and
transmission coefficients. We need a setup like

$\psi_{k}(x)=\begin{cases}Ae^{ikx}+Be^{-ikx}&x\to-\infty\\\
Ce^{i\sqrt{k^{2}-4\beta}x}&x\to+\infty\end{cases}.$ (33)

Notice that this is exactly the asymptotics of our eigenstate, and it
corresponds to particles going from left to right, and being either reflected
or transmitted through the potential. The coefficients are given by

$R=\frac{|B|^{2}}{|A|^{2}},\qquad
T=\frac{\sqrt{k^{2}-4\beta}}{k}\frac{|C|^{2}}{|A|^{2}}=\frac{(\mu+\eta)}{(\mu-\eta)}\frac{|C|^{2}}{|A|^{2}}.$
(34)

Although it is not intuitive, if we instead set up a wave coming from the
right, the reflection and transmission coefficients would be the same, even
though the particle would see a drop in the potential. This follows from very
general properties of the Schrödinger equationLandau and Lifshitz (2013).

Let’s first consider the regime below the barrier, where transmission is
trivially $T_{below}=0$. We apply (28), with $\mu-\eta=ik$ from (31). To
evaluate the reflection coefficient we first divide the amplitude multiplying
$e^{-ikx}$ by the amplitude of $e^{ikx}$.

$\frac{B}{A}=\frac{\Gamma(\mu-\eta)\Gamma(1-\mu+\nu)\Gamma(-\mu-\nu)}{\Gamma(-\nu-\eta)\Gamma(\nu+1-\eta)\Gamma(-\mu+\eta)}$
(35)

In view of (31), we can write $\mu=a+bi$ and $\eta=a-bi$ for real $a$ and $b$.
We also substitute $\nu=\alpha$, which is also real.

$\frac{B}{A}=\frac{\Gamma(2bi)\Gamma(1-a-bi+\alpha)\Gamma(-\alpha-a-
bi)}{\Gamma(-\alpha-a+bi)\Gamma(\alpha+1-a+bi)\Gamma(-2bi)}$ (36)

Notice that $B/A=1/(B/A)^{*}$, and since we have
$R_{below}=\frac{|B|^{2}}{|A|^{2}}$, we have

$R_{below}=1,\qquad T_{below}=0,$ (37)

as expected.

Now, for the more interesting process, where we scatter above the barrier. In
this situation, the reflection coefficient is a formula identical to (35), but
instead we substitute (32). When we calculate the magnitude squared, we must
remember that both $\mu$ and $\eta$ are purely imaginary.

$R=\frac{\Gamma(\mu-\eta)\Gamma(1-\mu+\alpha)\Gamma(-\mu-\alpha)}{\Gamma(-\alpha-\eta)\Gamma(\alpha+1-\eta)\Gamma(-\mu+\eta)}\frac{\Gamma(-\mu+\eta)\Gamma(1+\mu+\alpha)\Gamma(\mu-\alpha)}{\Gamma(-\alpha+\eta)\Gamma(\alpha+1+\eta)\Gamma(\mu-\eta)}.$
(38)

Cancelling several terms and applying the gamma reflection formula

$\Gamma(z)\Gamma(1-z)=\frac{\pi}{\sin(\pi z)},$ (39)

we write

$R=\frac{\sin(\pi(-\alpha+\eta))\sin(\pi(-\alpha-\eta))}{\sin(\pi(-\alpha+\mu))\sin(\pi(-\alpha-\mu))}.$
(40)

Since $\alpha$ is real, and $\mu$ and $\eta$ are imaginary, we can apply the
trigonometric identity

$\sin(x+iy)\sin(x-iy)=\sin^{2}(x)+\sinh^{2}(y).$ (41)

Therefore,

$R=\frac{\sin^{2}(\pi\alpha)+\sinh^{2}\left(\frac{\pi}{2}(k-\sqrt{k^{2}-4\beta})\right)}{\sin^{2}(\pi\alpha)+\sinh^{2}\left(\frac{\pi}{2}(k+\sqrt{k^{2}-4\beta})\right)}.$
(42)

To evaluate transmission we divide the coefficient multiplying the
$e^{i\sqrt{k^{2}-4\beta}x}$ in (24) by the coefficient of $e^{ikx}$ in (28):

$\frac{C}{A}=\frac{\Gamma(1-\mu+\alpha)\Gamma(-\mu-\alpha)}{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}.$
(43)

The transmission coefficient is then

$T=\frac{\sqrt{k^{2}-4\beta}}{k}\frac{|C|^{2}}{|A^{2}|}=\frac{(\mu+\eta)}{(\mu-\eta)}\frac{\Gamma(1-\mu+\alpha)\Gamma(-\mu-\alpha)}{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}\frac{\Gamma(1+\mu+\alpha)\Gamma(\mu-\alpha)}{\Gamma(1+\mu+\eta)\Gamma(\mu-\eta)}.$
(44)

To simplify it we need (39) and the gamma shift identity

$\Gamma(z+1)=z\Gamma(z),$ (45)

we have

$T=\frac{\sin(\pi(\mu-\eta))\sin(\pi(-\mu-\eta))}{\sin(\pi(\mu-\alpha))\sin(\pi(-\mu-\alpha))}.$
(46)

Using (41) again:

$T=\frac{\sinh(\pi
k)\sinh\left(\pi\sqrt{k^{2}-4\beta}\right)}{\sin^{2}(\pi\alpha)+\sinh^{2}\left(\frac{\pi}{2}(k+\sqrt{k^{2}-4\beta})\right)}.$
(47)

As expected, we have $R+T=1$, in view of the hyperbolic identity

$\sinh(a)\sinh(b)=\sinh^{2}\left(\frac{a+b}{2}\right)-\sinh^{2}\left(\frac{a-b}{2}\right).$
(48)

Our results agree in the suitable limits with those for the modified Pöschl-
Teller potentialFlügge (2012); Landau and Lifshitz (2013), the Woods-Saxon
potentialLandau and Lifshitz (2013) and the Morse-Feshbach potentialMorse and
Feshbach (1953).

##  VI Integral identity

The continuum eigenstates states of the Rosen-Morse potential should form an
orthogonal basis for the space of scattering states. However, in order to
create an integral transform that implements this change of basis we need to
evaluate their normalization. Before doing that, we prove a very useful
integral identity, which highlights that the generalized Legendre functions
have interesting properties of their own, and are not just a simple
reparametrization of the hypergeometric function.

We start by defining the function

$y(x)=D^{\mu,\eta}_{\nu}\frac{dD^{\rho,\gamma}_{\sigma}}{dx}-D^{\rho,\gamma}_{\sigma}\frac{dD^{\mu,\eta}_{\nu}}{dx}.$
(49)

If we differentiate this expression:

$\frac{dy}{dx}=D^{\mu,\eta}_{\nu}\frac{d^{2}D^{\rho,\gamma}_{\sigma}}{dx^{2}}-D^{\rho,\gamma}_{\sigma}\frac{d^{2}D^{\mu,\eta}_{\nu}}{dx^{2}}.$
(50)

Notice we can rearrange (3) as

$\frac{d^{2}D^{\mu,\eta}_{\nu}}{dx^{2}}=\frac{2x}{1-x^{2}}\frac{dD^{\mu,\eta}_{\nu}}{dx}+\left[\frac{\mu^{2}+2\mu\eta
x+\eta^{2}}{(1-x^{2})^{2}}-\frac{\nu(\nu+1)}{1-x^{2}}\right]D^{\mu,\eta}_{\nu},$
(51)

and substituting (51) into (50) we get

$\frac{dy}{dx}-\frac{2x}{1-x^{2}}y=\left[\frac{\rho^{2}+\gamma^{2}-\mu^{2}-\eta^{2}+2x(\rho\gamma-\mu\eta)}{(1-x^{2})^{2}}-\frac{\sigma(\sigma+1)-\nu(\nu+1)}{1-x^{2}}\right]D^{\mu,\eta}_{\nu}D^{\rho,\gamma}_{\sigma}.$
(52)

This is a first order differential equation which can be solved by multiplying
by the integrating factor
$\exp\left(-\int\frac{2xdx}{1-x^{2}}\right)=(1-x^{2})$. At the end, we obtain

$\displaystyle\frac{d}{dx}$
$\displaystyle\left[(1-x^{2})\left(D^{\mu,\eta}_{\nu}\frac{dD^{\rho,\gamma}_{\sigma}}{dx}-D^{\rho,\gamma}_{\sigma}\frac{dD^{\mu,\eta}_{\nu}}{dx}\right)\right]$
(53)
$\displaystyle=\left[\frac{\rho^{2}+\gamma^{2}-\mu^{2}-\eta^{2}+2x(\rho\gamma-\mu\eta)}{1-x^{2}}-\sigma(\sigma+1)+\nu(\nu+1)\right]D^{\mu,\eta}_{\nu}D^{\rho,\gamma}_{\sigma}.$

From all of this, we have proved

$\displaystyle\int_{a}^{b}$
$\displaystyle\left[\frac{\rho^{2}+\gamma^{2}-\mu^{2}-\eta^{2}+2x(\rho\gamma-\mu\eta)}{1-x^{2}}-\sigma(\sigma+1)+\nu(\nu+1)\right]D^{\mu,\eta}_{\nu}D^{\rho,\gamma}_{\sigma}dx$
(54)
$\displaystyle=\left[(1-x^{2})\left(D^{\mu,\eta}_{\nu}\frac{dD^{\rho,\gamma}_{\sigma}}{dx}-D^{\rho,\gamma}_{\sigma}\frac{dD^{\mu,\eta}_{\nu}}{dx}\right)\right]\bigg{|}_{a}^{b}.$

##  VII Normalization and spectral measure

Recall the mapping between the wave number $k$ and our generalized Legendre
functions

$\psi_{k}(x)=\begin{cases}D^{-\frac{\sqrt{4\beta-k^{2}}}{2}+\frac{ik}{2},-\frac{\sqrt{4\beta-k^{2}}}{2}-\frac{ik}{2}}_{\alpha}(\tanh
x),&|k|\leq 2\sqrt{\beta}\\\
D^{\frac{i\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}}{2}+\frac{ik}{2},\frac{i\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}}{2}-\frac{ik}{2}}_{\alpha}(\tanh
x),&|k|>2\sqrt{\beta}\end{cases}.$ (55)

Our goal is to evaluate the integral

$I=\int_{-\infty}^{\infty}\psi_{p}^{*}(x)\psi_{k}(x)dx.$ (56)

For this, we use (54). We set $\nu=\sigma=\alpha$ and notice that, for the
choice of parameters in (55), we have $\mu\eta=\rho\gamma=\beta$,
$\mu^{2}+\eta^{2}=2\beta-k^{2}$ and $\rho^{2}+\gamma^{2}=2\beta-p^{2}$. In
this setting, lots of simplifications occur in (54) and we obtain

$(k^{2}-p^{2})\int_{a}^{b}\frac{\bar{D}^{ip}(u)^{*}\bar{D}^{ik}(u)}{1-u^{2}}du=\left[(1-u^{2})\left(\bar{D}^{ik}\frac{d\bar{D}^{ip*}}{du}-\bar{D}^{ip*}\frac{d\bar{D}^{ik}}{du}\right)\right]\bigg{|}_{a}^{b},$
(57)

where we have introduced $\bar{D}^{ik}$ as an abreviation for the generalized
Legendre functions in (55).

After doing the change of variables $u=\tanh x$:

$\displaystyle(k^{2}-p^{2})\int_{a}^{b}$ $\displaystyle\bar{D}^{ip*}(\tanh
x)\bar{D}^{ik}(\tanh x)dx$ (58) $\displaystyle=\left(\bar{D}^{ik}(\tanh
x)\frac{d}{dx}\bar{D}^{ip*}(\tanh x)-\bar{D}^{ip*}(\tanh
x)\frac{d}{dx}\bar{D}^{ik}(\tanh x)\right)\bigg{|}_{a}^{b}.$

And we may finally express our orthogonality integral as

$I=\frac{1}{k^{2}-p^{2}}\left(\bar{D}^{ik}(\tanh
x)\frac{d}{dx}\bar{D}^{ip*}(\tanh x)-\bar{D}^{ip*}(\tanh
x)\frac{d}{dx}\bar{D}^{ik}(\tanh x)\right)\bigg{|}_{-\Lambda}^{\Lambda},$ (59)

where we have introduced a cutoff $\Lambda\to\infty$ to regularize the
integral.

We need to subtract two terms which are very similar. Let’s focus on the first
one. The asymptotics at $\Lambda$ is only nonzero for $|k|>2\sqrt{\beta}$,
when it’s equal to

$\bar{D}^{ik}(\tanh\Lambda)=2^{-\eta}e^{i\sqrt{k^{2}-4\beta}\Lambda},$ (60)

and the derivative is

$\frac{d}{dx}\bar{D}^{ip*}(\tanh\Lambda)=-i2^{-\gamma}\operatorname{sgn}(p)\sqrt{p^{2}-4\beta}e^{-i\operatorname{sgn}(p)\sqrt{p^{2}-4\beta}\Lambda},$
(61)

so that the total contribution at $\Lambda$ is given by

$T(\Lambda)=-i2^{-\eta-\gamma}\left(\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}+\operatorname{sgn}(p)\sqrt{p^{2}-4\beta}\right)e^{i\left(\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}-\operatorname{sgn}(p)\sqrt{p^{2}-4\beta}\right)\Lambda}.$
(62)

To make further progress, we need the distributional relations

$\lim_{\Lambda\to\infty}\frac{\sin(\Lambda
x)}{x}=\pi\delta(x),\qquad\lim_{\Lambda\to\infty}\frac{\cos(\Lambda x)}{x}=0.$
(63)

They allow us to get rid of the regulator with the substitution

$\frac{e^{-i(f(k)-f(p))\Lambda}}{f(k)-f(p)}\to-i\pi\delta(f(p)-f(k)).$ (64)

Using some algebra, when we divide by $k^{2}-p^{2}$ we obtain

$\frac{T(\Lambda)}{k^{2}-p^{2}}=-i2^{-\eta-\gamma}\frac{e^{i\left(\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}-\operatorname{sgn}(p)\sqrt{p^{2}-4\beta}\right)\Lambda}}{\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}-\operatorname{sgn}(p)\sqrt{p^{2}-4\beta}}\to\pi
2^{-\eta-\gamma}\delta\left(\operatorname{sgn}(k)\sqrt{k^{2}-4\beta}-\operatorname{sgn}(p)\sqrt{p^{2}-4\beta}\right).$
(65)

We now use the identity

$\delta(f(x))=\sum_{i}\frac{\delta(x-x_{i})}{|f^{\prime}(x_{i})|},$ (66)

where $x_{i}$ denotes the zeros of $f(x)$. This is where the
$\operatorname{sgn}$ function is seen to be important, because the only zero
is when $k=p$. We get, after substituting $\gamma=-\eta$:

$I(\Lambda)=\frac{T(\Lambda)}{k^{2}-p^{2}}=\pi\frac{\sqrt{k^{2}-4\beta}}{|k|}\delta\left(k-p\right).$
(67)

We now look at the more interesting $x=-\Lambda$ asymptotics. It’s given by
(28):

$\bar{D}^{ik}(\tanh-\Lambda)=2^{-\eta}\left[\frac{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}{\Gamma(1-\mu+\alpha)\Gamma(-\mu-\alpha)}e^{-(\mu-\eta)\Lambda}+\frac{\Gamma(1-\mu-\eta)\Gamma(\mu-\eta)}{\Gamma(-\alpha-\eta)\Gamma(\alpha+1-\eta)}e^{(\mu-\eta)\Lambda}\right],$
(68)

and the derivative is

$\displaystyle\frac{d}{dx}$
$\displaystyle\bar{D}^{ip*}(\tanh-\Lambda)=2^{-\gamma}\bigg{[}(\rho-\gamma)\frac{\Gamma(1-\rho-\gamma)\Gamma(-\rho+\gamma)}{\Gamma(1-\rho+\alpha)\Gamma(-\rho-\alpha)}e^{-(\rho-\gamma)\Lambda}$
(69)
$\displaystyle+(\gamma-\rho)\frac{\Gamma(1-\rho-\gamma)\Gamma(\rho-\gamma)}{\Gamma(-\alpha-\gamma)\Gamma(\alpha+1-\gamma)}e^{(\rho-\gamma)\Lambda}\bigg{]}.$

When we multiply the terms in (69), the crossed terms will have exponential
factors like

$e^{(\rho-\mu-\gamma+\eta)\Lambda}=e^{-i(p+k)\Lambda}.$ (70)

We get rid of the regulator by doing something analogous to (64). We can
gobble a $p+k$ term from the denominator of (59) and write

$\frac{e^{-i(p+k)\Lambda}}{p+k}\to-i\pi\delta(p+k).$ (71)

Notice, though, that once we’ve gathered all contibutions from the first term
of (59), the function multiplying it has to be symmetric with respect to the
interchange $k\to-p$ and $p\to-k$. However, we would get the same contribution
from the second term of (59), because it’s the same calculation under the same
relabeling. Since we need to subtract these two, the factors multiplying
$\delta(p+k)$ necessarily vanish, which means we can completely ignore the
crossed terms in the multiplication of (68) and (69), and focus only on the
uncrossed ones.

The terms are

$T_{1}=-i2^{-\eta-\gamma}p\frac{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}{\Gamma(1-\mu+\alpha)\Gamma(-\mu-\alpha)}\frac{\Gamma(1-\rho-\gamma)\Gamma(-\rho+\gamma)}{\Gamma(1-\rho+\alpha)\Gamma(-\rho-\alpha)}e^{-i(k-p)\Lambda},$
(72)

$T_{2}=i2^{-\eta-\gamma}p\frac{\Gamma(1-\mu-\eta)\Gamma(\mu-\eta)}{\Gamma(-\alpha-\eta)\Gamma(\alpha+1-\eta)}\frac{\Gamma(1-\rho-\gamma)\Gamma(\rho-\gamma)}{\Gamma(-\alpha-\gamma)\Gamma(\alpha+1-\gamma)}e^{i(k-p)\Lambda}$
(73)

To compute the second term of (59) we simply send $k\to-p$ and $p\to-k$. With
this transformation, $(k-p)$ stays unchanged, and we have $\mu\to\rho$,
$\eta\to\gamma$, and vice-versa. Since the gamma functions are all symmetric
with respect to parameter changes, the only thing that changes is the $p$
factor at the beginning, which goes to $-k$. However, since we need to
subtract, this gives a $(k+p)$ factor in both $T_{1}$ and $T_{2}$, which
cancels the one in the denominator in (59).

Finally, we gobble the $(k-p)$ factor in the denominator of (59) to perform a
substitution like the one in (71), and at the end we have the total terms at
$-\Lambda$:

$I_{1}(-\Lambda)=-\pi
2^{-\eta-\gamma}\frac{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}{\Gamma(1-\mu+\alpha)\Gamma(-\mu-\alpha)}\frac{\Gamma(1-\rho-\gamma)\Gamma(-\rho+\gamma)}{\Gamma(1-\rho+\alpha)\Gamma(-\rho-\alpha)}\delta(k-p),$
(74)

$I_{2}(-\Lambda)=-\pi
2^{-\eta-\gamma}\frac{\Gamma(1-\mu-\eta)\Gamma(\mu-\eta)}{\Gamma(-\alpha-\eta)\Gamma(\alpha+1-\eta)}\frac{\Gamma(1-\rho-\gamma)\Gamma(\rho-\gamma)}{\Gamma(-\alpha-\gamma)\Gamma(\alpha+1-\gamma)}\delta(k-p),$
(75)

where

$I=I(\Lambda)-I_{1}(-\Lambda)-I_{2}(-\Lambda).$ (76)

The parameters $\rho$ and $\gamma$ are obtained by substituting $p$ in either
(31) or (32) and complex conjugating. The delta function forces $p=k$. If we
have $k^{2}<4\beta$, then we should replace $\rho=\eta$ and $\gamma=\mu$. If
we have $k^{2}>4\beta$, we should replace $\rho=-\mu$ and $\gamma=-\eta$.
Let’s do $k^{2}>4\beta$ first.

$I_{1}(-\Lambda)=-\pi\frac{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}{\Gamma(1-\mu+\alpha)\Gamma(-\mu-\alpha)}\frac{\Gamma(1+\mu+\eta)\Gamma(\mu-\eta)}{\Gamma(1+\mu+\alpha)\Gamma(\mu-\alpha)}\delta(k-p),$
(77)

$I_{2}(-\Lambda)=-\pi\frac{\Gamma(1-\mu-\eta)\Gamma(\mu-\eta)}{\Gamma(-\alpha-\eta)\Gamma(\alpha+1-\eta)}\frac{\Gamma(1+\mu+\eta)\Gamma(-\mu+\eta)}{\Gamma(-\alpha+\eta)\Gamma(\alpha+1+\eta)}\delta(k-p).$
(78)

Applying the shift and reflection identities of the gamma function, we
simplify:

$I_{1}(-\Lambda)=-\pi\frac{\mu+\eta}{-\mu+\eta}\frac{\sin(\pi(\mu-\alpha))\sin(\pi(-\mu-\alpha))}{\sin(\pi(\mu-\eta))\sin(\pi(\mu+\eta))}\delta(k-p),$
(79)

$I_{2}(-\Lambda)=-\pi\frac{\mu+\eta}{-\mu+\eta}\frac{\sin(\pi(\eta-\alpha))\sin(\pi(-\eta-\alpha))}{\sin(\pi(\mu-\eta))\sin(\pi(\mu+\eta))}\delta(k-p).$
(80)

And trigonometric identities:

$I_{1}(-\Lambda)=-\pi\frac{\sqrt{k^{2}-4\beta}}{|k|}\frac{\sin^{2}(\pi\alpha)+\sinh^{2}\left(\frac{\pi}{2}(\sqrt{k^{2}-4\beta}+k)\right)}{\sinh(\pi
k)\sinh(\pi\sqrt{k^{2}-4\beta})}\delta(k-p),$ (81)

$I_{2}(-\Lambda)=-\pi\frac{\sqrt{k^{2}-4\beta}}{|k|}\frac{\sin^{2}(\pi\alpha)+\sinh^{2}\left(\frac{\pi}{2}(\sqrt{k^{2}-4\beta}-k)\right)}{\sinh(\pi
k)\sinh(\pi\sqrt{k^{2}-4\beta})}\delta(k-p),$ (82)

When evaluating $I(\Lambda)-I_{1}(-\Lambda)-I_{2}(-\Lambda)$ for this region,
we can simplify the sum using (48).

Now the region $k^{2}<4\beta$:

$I_{1}=I_{2}=-\pi
2^{-\eta-\mu}\frac{\Gamma(1-\mu-\eta)\Gamma(-\mu+\eta)}{\Gamma(1-\mu+\alpha)\Gamma(-\mu-\alpha)}\frac{\Gamma(1-\eta-\mu)\Gamma(-\eta+\mu)}{\Gamma(1-\eta+\alpha)\Gamma(-\eta-\alpha)}\delta(k-p),$
(83)

Substituting $\mu$ and $\eta$ from (31):

$I_{1}=-\frac{\pi^{2}2^{-\sqrt{4\beta-k^{2}}}}{k\sinh(\pi
k)}\frac{\Gamma\left(1+\sqrt{4\beta-k^{2}}\right)^{2}}{\left|\Gamma\left(1+\alpha+\frac{\sqrt{4\beta-k^{2}}}{2}+\frac{ik}{2}\right)\right|^{2}\left|\Gamma\left(-\alpha+\frac{\sqrt{4\beta-k^{2}}}{2}+\frac{ik}{2}\right)\right|^{2}}\delta(k-p)$
(84)

Therefore, if we let

$w_{\alpha,\beta}(k)=\begin{cases}2\pi\frac{\sqrt{k^{2}-4\beta}}{|k|}\frac{\sin^{2}(\pi\alpha)+\sinh^{2}\left(\frac{\pi}{2}(\sqrt{k^{2}-4\beta}+k)\right)}{\sinh(\pi
k)\sinh(\pi\sqrt{k^{2}-4\beta})},&|k|>2\sqrt{\beta}\\\
\frac{2\pi^{2}2^{-\sqrt{4\beta-k^{2}}}}{k\sinh(\pi
k)}\frac{\Gamma\left(1+\sqrt{4\beta-k^{2}}\right)^{2}}{\left|\Gamma\left(1+\alpha+\frac{\sqrt{4\beta-k^{2}}}{2}+\frac{ik}{2}\right)\right|^{2}\left|\Gamma\left(-\alpha+\frac{\sqrt{4\beta-k^{2}}}{2}+\frac{ik}{2}\right)\right|^{2}},&|k|\leq
2\sqrt{\beta}\end{cases}$ (85)

we can write the normalization of the generalized Legendre functions as

$\int_{-\infty}^{\infty}\bar{D}_{\alpha,\beta}^{ip*}(\tanh
x)\bar{D}_{\alpha,\beta}^{ik}(\tanh x)dx=w_{\alpha,\beta}(k)\delta(k-p).$ (86)

We can setup something analogous to the Fourier transform. If we suppose an
arbitrary scattering state can be expanded into Legendre functions, we write

$f(x)=\int_{-\infty}^{\infty}F(k)\bar{D}^{ik}(\tanh x)dk$ (87)

In view of (86), we can find $F(k)$ by multiplying by another Legendre
function and integrating:

$\int_{-\infty}^{\infty}f(x)\bar{D}_{\alpha,\beta}^{ik*}(\tanh
x)dx=F(k)w_{\alpha,\beta}(k)$ (88)

But substituting (88) on (87) we have

$f(x)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{\bar{D}_{\alpha,\beta}^{ik*}(\tanh
y)\bar{D}_{\alpha,\beta}^{ik}(\tanh x)}{w_{\alpha,\beta}(k)}f(y)dydk$ (89)

But, since we know that

$f(x)=\int_{-\infty}^{\infty}\delta(x-y)f(y)dy,$ (90)

we must identify the transform on (89) as a delta function, and we obtain the
resolution of the identity for the scattering states:

$\int_{-\infty}^{\infty}\frac{\bar{D}_{\alpha,\beta}^{ik*}(\tanh
y)\bar{D}_{\alpha,\beta}^{ik}(\tanh x)}{w_{\alpha,\beta}(k)}dk=\delta(x-y)$
(91)

We also note that, when $\beta=0$, (85) reduces to

$w_{\alpha,0}(k)=2\pi\left(1+\frac{\sin^{2}(\pi\alpha)}{\sinh^{2}(\pi
k)}\right),$ (92)

and in this limit our functions reduce to regular associated Legendre
functions like

$\bar{D}_{\alpha,0}^{ik}(\tanh x)=\Gamma(1-ik)P^{ik}_{\alpha}(\tanh x),$ (93)

so that when we substitute this on (91), we have the resolution of the
identity for scattering states of the symmetric potential with respect to
Legendre functions:

$\int_{-\infty}^{\infty}\frac{kdk\sinh(\pi k)P_{\alpha}^{-ik}(\tanh
y)P_{\alpha}^{ik}(\tanh x)}{2(\sinh^{2}(\pi
k)+\sin^{2}(\pi\alpha))}=\delta(x-y),$ (94)

which matches with the one obtained by means of path integralsGrosche and
Steiner (1998).

##  VIII Conclusions

In this article, we have generalized associated Legendre functions and shown
how they naturally describe states of the Rosen-Morse potential, allowing for
a completely classical solution of the system, from which we can obtain
reasonably simple expressions for the reflection and transmission coefficients
and the normalization of the continuum spectrum, all of which reduced in the
suitable limit to the ones reported in the literature. The integral identity
obtained in this work shows that these functions are interesting in their own
right, and more investigation of their properties should give more insight
into exactly solvable potentials.

## References

  * Natanzon (1979) G. A. Natanzon, Theor. Math. Phys. 38, 146 (1979).
  * Arfken and Weber (2005) G. Arfken and H. Weber,  _Mathematical Methods For Physicists_ (Elsevier Science, 2005).
  * Weber (2007) H. Weber, Centr. Eur. J. Math. 5, 415 (2007).
  * Compean and Kirchbach (2005) C. B. Compean and M. Kirchbach, J. Phys. A: Math. Gen. 39, 547 (2005).
  * Rosen and Morse (1932) N. Rosen and P. M. Morse, Phys. Rev. 42, 210 (1932).
  * Kleinert and Mustapic (1992) H. Kleinert and I. Mustapic, J. Math. Phys. 33, 643 (1992).
  * Lin (1997) D.-H. Lin, . Phys. A: Math. Gen. 30, 3201 (1997).
  * Grosche (1989) C. Grosche, J. Phys. A 22, 5073 (1989).
  * Flügge (2012) S. Flügge,  _Practical Quantum Mechanics_ (Springer Berlin Heidelberg, 2012).
  * Woods and Saxon (1954) R. D. Woods and D. S. Saxon, Phys. Rev. 95, 577 (1954).
  * Morse and Feshbach (1953) P. Morse and H. Feshbach,  _Methods of Theoretical Physics_ (McGraw-Hill, 1953).
  * Freitas (2018) F. L. Freitas,  (2018), arXiv:1802.06617 [math-ph] .
  * Dong and Lemus (2002) S.-H. Dong and R. Lemus, Int. J. Quantum Chem. 86, 265 (2002).
  * Magnus  _et al._ (1966) W. Magnus, F. Oberhettinger, and R. Soni,  _Formulas and Theorems for the Special Functions of Mathematical Physics_ (Springer Berlin Heidelberg, 1966).
  * Grosche (1990) C. Grosche, Fortschr. Phys. 38, 531 (1990).
  * Grosche and Steiner (1998) C. Grosche and F. Steiner,  _Handbook of Feynman Path Integrals_ (Springer Berlin Heidelberg, 1998).
  * Gradshteyn and Ryzhik (2014) I. Gradshteyn and I. Ryzhik,  _Table of Integrals, Series, and Products_ (Elsevier Science, 2014).
  * Abramowitz and Stegun (1965) M. Abramowitz and I. Stegun,  _Handbook of Mathematical Functions_ (Dover Publications, 1965).
  * Landau and Lifshitz (2013) L. Landau and E. Lifshitz,  _Quantum Mechanics: Non-Relativistic Theory_ (Elsevier Science, 2013).