# Correct-by-Construction Navigation Functions
with Application to Sensor Based Robot Navigation

Savvas G. Loizou1 and Elon D. Rimon2 1Dept. of Mechanical and Materials
Science Engineering, Cyprus University of Technology, Limassol, Cyprus.
<EMAIL_ADDRESS>of Mechanical Engineering, Technion Israel
Institute of Technology, Haifa, Israel<EMAIL_ADDRESS>
###### Abstract

This paper brings together the concepts of navigation transformation and
harmonic functions to form navigation functions that are correct-by-
construction in the sense that no tuning is required. The form of the
navigation function is explicitly related to the number of obstacles in the
environment. This enables application of navigation functions for autonomous
robot navigation in partially or fully unknown environments, with the
capability of on-the-fly adjustment of the navigation function when new
obstacles are discovered by the robot. Appropriate navigation controllers,
applicable to robots with local, sector bounded sensing, are presented and
analyzed for a kinematic point-mass robot and then for the dynamic point-mass
robot system. The closed form nature of the proposed navigation scheme
provides for online, fast-feedback based navigation. In addition to the
analytic guarantees, simulation studies are presented to verify the
effectiveness of the methodology.

##  I Introduction

Navigation Functions [6] is an active research topic in the field of robotic
navigation. There are numerous successful examples in the literature of
Navigation Function applications to navigation problems of increasing
complexity [5]. One of the major challenges when setting up a Navigation
Functions has to do with its tuning. A construction can be a candidate
Navigation Function, but to actually have the Navigation Function properties,
an appropriate selection of the $k$ parameter has to be made. Such a $k$ is
always guaranteed to exist, however calculation of its lower bound is not
straightforward. Several solutions have appeared in the literature to tackle
the issue of sensor based reactive planning, see e.g. [1], [9] however setting
up a correct-by-construction Navigation Function on an arbitrary workspace is
still an open issue.

Recent results on the Navigation Transformation [7] and Harmonic Function
based Navigation Functions [8] indicated the feasibility of a tuning
controller. However the analysis provided in the current work demonstrates
that a tuning controller is not required and that for Harmonic Function based
Navigation Functions utilizing the Navigation Transformation, parameter $k$,
when explicitly provided as a function of the number of workspace obstacles is
sufficient to generate a correct behavior. Kinematic and dynamic controllers
were produced to handle environments with only partial workspace knowledge and
robot sensing limited in a bounded sensing sector. The contributions of the
current work are summarized in the list below:

  * •

A correct-by-construction Navigation Function

  * •

On-the-fly addition of obstacles in the workspace

  * •

A relaxed version of the Morse Property of the Navigation Functions that
provides equivalent performance characteristics

  * •

Navigation capability with arbitrary, non-zero, sector-bounded local sensing

  * •

Kinematic controller whose vector field does not vanish at the vicinity of
saddle points

  * •

Dynamic controller with critical damping and bounded maximum velocity

The rest of the paper is organized as follows: Section II presents preliminary
notions, sections III and IV present the construction and analysis of the
navigation function, section V presents the controller design, section VI
presents simulation results and section VII concludes the paper.

##  II Preliminaries

This section introduces navigation functions terminology, the two mobile robot
types, the robot sensor model, then the navigation function and controller
design problem.

Basic terminology: If $K$ is a set, then The $n$-dimensional _sphere world_ as
defined in [6] is a compact and connected set in $\mathbb{R}^{n}$ that is
bounded by an outer $(n-1)$-dimensional sphere and contains $M$ disjoint
internal spheres that represent obstacles. It serves as a topological model
for robot workspaces defined below.

Given a smooth function $f(\cdot):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, we
denote the Jacobian matrix of this function as $J_{f}(\cdot)$ and the Jacobian
determinant by $\left\|J_{f}(\cdot)\right\|$. Given a smooth function
$\phi:\mathbb{R}^{n}\rightarrow\mathbb{R}$ we denote the Hessian matrix of
$\phi$ as $\mathcal{H}_{\phi}\left(\cdot\right)$.

Robot workspaces: The navigation functions will be constructed in the
following point worlds [7].

###### Definition 1

Let $P_{i}\in\mathbb{R}^{n}$ for $i=1\ldots M$ be $M$ discrete point obstacles
in $\mathbb{R}^{n}$. Then a _point world_ is defined as the set
$\mathcal{P}=\mathbb{R}^{n}-\\{P_{1},\ldots,P_{M}\\}$.

We will restrict our attention to robot workspaces that are topologically
equivalent to sphere worlds as follows.

###### Definition 2

The _robot workspace_ $\mathcal{W}$ is a compact and connected set in
$\mathbb{R}^{n}$ topologically equivalent to a sphere world.

A valid robot workspace has an outer boundary denoted ${\cal O}_{0}$ and
disjoint internal obstacles denoted ${\cal O}_{i}$ for $i=1\ldots M$. By the
nature of the underlying sphere world, the _interior_ of a robot workspace is
diffeomorphic to a point world $\mathcal{P}$. The mapping of complex robot
workspaces to their point-world models will use the following coordinate
transformation.

######  Definition 3 ([7])

A _navigation transformation_ is a smooth coordinate transformation (a
diffeomorphism)
$\Phi:\stackrel{{\scriptstyle\circ}}{{\mathcal{W}}}\rightarrow\mathcal{P}$,
that maps the interior of the robot workspace to a point world.

Robot models: This paper considers two types of mobile robots that navigate in
planar environments populated by obstacles. The first is the trivial kinematic
integrator that models a kinematic point-mass robot:

$\dot{x}(t)=u(t)$ (1)

where $x\in\mathbb{R}^{2}$ is the robot position and $u(t)\in\mathbb{R}^{2}$
is a piecewise continuous control input. The more demanding model is the
second-order dynamic integrator for a point robot of mass $m$:

$m\cdot\ddot{x}(t)=f(t)$ (2)

where $(x,\dot{x})$ is the robot state and $f(t)\in\mathbb{R}^{2}$ is a
piecewise continuous control input representing force applied to the robot.

Sensor model: Instead of just considering a uniform sensing radius for the
robot, this paper will use a more general robot centered _symmetric sensing
sector_ defined as

$S_{r,\theta}(x,\dot{x})=\left\\{(r^{\prime},\theta^{\prime})\left|\;\right.0\leq
r^{\prime}\leq r,\;\left|\theta^{\prime}\right|\leq\frac{\theta}{2}\right\\}$
(3)

where the sensing sector’s pole is at $x$ and the sensing sector’s axis is
aligned with the robot velocity $\dot{x}$.

Navigation function: Navigation functions were originally defined as analytic
Morse functions [6]. This paper introduces a smooth version that relaxes the
Morse property as follows. Denote by $x_{0}$ and $x_{d}$ the robot initial and
destination positions.

###### Definition 4

Let $\mathcal{W}\subset\mathbb{R}^{n}$ be a compact connected smooth manifold
with boundary representing a valid robot workspace. A map
$\varphi:\mathcal{W}\rightarrow[0,1]$, forms a  _navigation function_ if it
is:

  1. 1.

Smooth on $\mathcal{W}$.

  2. 2.

Polar on $\mathcal{W}$, $\varphi$ has a unique minimum at
$x_{d}\in\stackrel{{\scriptstyle\circ}}{{\mathcal{W}}}$ such that
$\varphi(x_{d})=0$.

  3. 3.

The critical points of $\varphi$ except $x_{d}$ are isolated saddle points
with basins of attraction of zero measure.

  4. 4.

Admissible on $\mathcal{W}$, $\varphi$ attains maximal value of unity on
$\partial\mathcal{W}$.

Problem Description: In the problem considered in this paper the workspace is
assumed to contain static obstacles. Each obstacle has a neighborhood that is
disjoint from the neighborhoods of the other obstacles. In each neighborhood a
local _navigation transformation_ is provided that maps the obstacle
neighborhood to a point-obstacle neighborhood in the corresponding point
world. However, the obstacle locations are _not_ known to the robot. As the
robot navigates to the destination it detects new obstacles. The position and
shape of the obstacles are determined at _discrete time instants_ and
accumulated by the robot during navigation. We can now state the problem
considered in this paper.

Problem statement: Construct correct-by-construction navigation functions and
provide the corresponding feedback control law so that the robot systems (1)
and (2) can navigate from any initial position to a specified destination
avoiding collisions with obstacles in the environment.

Note that navigation functions guarantee arrival to the destination from any
initial position except for a set of zero measure, which is unstable in the
sense that small disturbances will set the robot on a path that reaches the
destination.

##  III Construction of Navigation Function based on Harmonic Potentials

Assume a valid robot workspace $\mathcal{W}$ and a navigation transformation
$\Phi:\mathcal{W}\rightarrow\mathcal{P}$ having the additional property that
$\lim\limits_{x\rightarrow\mathcal{O}_{0}}\left\|\Phi(x)\right\|=\infty$,
where $\mathcal{O}_{0}$ is the workspace outer boundary (see [8] for such an
instance). A _point obstacle harmonic potential_ is defined as

$\phi_{i}(h)=\ln\left(\left\|h-P_{i}\right\|^{2}\right),$

and the _destination harmonic potential_ is defined as

$\phi_{d}(h)=\ln\left(\left\|h-P_{d}\right\|^{2}\right)$

where $P_{d}=\Phi(x_{d})$. Note that $-\phi_{i}(P_{i})=+\infty$ while
$\phi_{d}(P_{d})=-\infty$. The _harmonic point world potential function,_
$\phi_{k}:\mathcal{P}\rightarrow R$, is defined as

$\phi_{k}(h)=\phi_{d}(h)-\frac{\mbox{\small$1$}}{k}\sum\limits_{i=1}^{M}\phi_{i}(h)$

where $k$ a positive parameter and $M$ is the total number of internal
obstacles in $\mathcal{W}$. The function $\phi_{k}$ attains values in the
extended real line while navigation functions are required to attain values in
$[0,1]$. The mapping to the unit interval is achieved by composition with the
function

$\sigma\left(x\right)\triangleq\frac{e^{x}}{1+e^{x}}\,.$ (4)

that maps the extended real line to the unit interval $[0,1]$, such that
$\sigma(-\infty)=0$ and $\sigma(+\infty)=1$. The candidate _workspace
navigation function,_ $\varphi_{k}:\mathcal{W}\rightarrow R$, is obtained by
composition with the workspace to point-world coordinate transformation:

$\varphi_{k}(x)=\sigma\circ\phi_{k}\circ\Phi(x).$ (5)

##  IV Navigation Function Analysis

This section builds up several properties of the point world harmonic
navigation function, $\phi_{k}(h)$, that will lead to the conclusion that
$\varphi_{k}(x)=\sigma\circ\phi_{k}\circ\Phi(x)$ is a navigation function. We
start with the following result.

###### Lemma 1

The point world harmonic potential $\phi_{k}:\mathcal{P}\rightarrow R$ is free
of any local minima.

###### Proof:

Since $\phi_{k}:\mathcal{P}\rightarrow R$ is a superposition of the harmonic
potentials $\phi_{i}(\cdot)$, $i\in\left\\{1,\ldots M\right\\}$ and
$\phi_{d}(\cdot)$, it will also be harmonic, which implies by definition that

$\Delta\phi_{k}(h)\triangleq 0.$

Since $\Delta\phi_{k}(h)=\mathrm{tr}(\mathcal{H}_{\phi_{k}}(h))$, we will have
that $\mathrm{tr}(\mathcal{H}_{\phi_{k}}(h))=0$. The trace of a matrix is the
first principal invariant of a matrix, hence
$I_{1}=\mathrm{tr}(\mathcal{H}_{\phi_{k}}(h))=\lambda_{1}+\lambda_{2}=0$ This
implies that at the critical points we cannot have both eigenvalues positive,
hence a local minimum is not possible. ∎

Next consider the _attractivity_ of $\phi_{k}$, i.e. that the negated gradient
trajectories of $\phi_{k}$ do not escape to infinity in $\mathcal{P}$.

###### Lemma 2

For $k>M$ we have that
$\lim\limits_{\left\|h\right\|\rightarrow\infty}\phi_{k}(h)=+\infty$.

###### Proof:

By construction

$\phi_{k}(h)=\ln\left(\left\|h-P_{d}\right\|^{2}\right)-\frac{1}{k}\sum_{i=1}^{M}\ln\left(\left\|h-P_{i}\right\|^{2}\right).$

Hence

$\lim\limits_{\left\|h\right\|\rightarrow\infty}\phi_{k}(h)=\ln\left(\lim\limits_{\left\|h\right\|\rightarrow\infty}\left\|h\right\|^{2(k-M)}\right)=+\infty$

for $k>M$. ∎

Let us next verify that $\phi_{k}$ has well behaved critical points in
$\mathcal{P}$. The following result concerns special obstacle arrangements for
which $\phi_{k}$ has  _degenerate_ critical points.

###### Proposition 1

For every choice of $k\geq 0$, there exist workspace arrangements for which
$\phi_{k}$ is non-Morse in $\mathcal{P}$.

###### Proof:

A smooth function is called a Morse function iff all its critical points are
non-degenerate. Degenerate critical points, are critical points where the
Hessian has zero eigenvalue(s). Let $\mathcal{C}(\phi_{k})$ denote the set of
critical points of $\phi_{k}$. We will denote with $h_{d}=h-q_{d}$ and
$h_{i}=h-P_{i}$. Taking the gradient of $\phi_{k}$ we have:

$\nabla\phi_{k}(h)=2\frac{h_{d}}{\left\|h_{d}\right\|^{2}}-\frac{2}{k}\sum\limits_{i=1}^{M}\frac{h_{i}}{\left\|h_{i}\right\|^{2}}$
(6)

Now since the Hessian is a linear operator, the Hessian of $\phi_{k}$ can be
written as the superposition of the Hessians of $\phi_{i}$ and $\phi_{d}$:

$\mathcal{H}_{\phi_{k}}(h)=\mathcal{H}_{\phi_{d}}(h)+\sum\limits_{i=1}^{M}\mathcal{H}_{\phi_{i}}(h)$

Using eigendecomposition on the Hessians of $\phi_{d}$ and $\phi_{i}$, the
Hessian of $\phi_{k}$ assumes the following form:

$\mathcal{H}_{\phi_{k}}=Q_{d}\cdot\Lambda_{d}\cdot
Q_{d}^{T}-\frac{1}{k}\sum\limits_{i=1}^{M}Q_{i}\cdot\Lambda_{i}\cdot
Q_{i}^{T}$ (7)

where the eigenvalues of the Hessians of $\phi_{d}$ and $\phi_{i}$ form:

$Q_{x}=\left[\hat{h}_{x}\;\vdots\;J\cdot\hat{h}_{x}\right],\quad
x\in\left\\{d,i\right\\}$

with $J=\left[\begin{array}[]{cc}0&-1\\\ 1&0\end{array}\right]$. Observe that
$Q_{x}^{-1}=Q_{x}^{T}$. Now the eigenvalues of the Hessians of $\phi_{d}$ and
$\phi_{i}$ form the following diagonal matrix:

$\Lambda_{x}=\frac{2}{\left\|h_{x}\right\|^{2}}L$

where $L=\begin{bmatrix}-1&0\\\ 0&1\end{bmatrix}$. Since $\phi_{k}$ is a
harmonic function it will have two eigenvalues of equal magnitude and opposite
signs. To establish that $\mathcal{H}_{\phi_{k}}$ is degenerate we need to
show that there is a test direction $v$ and its perpendicular $J\cdot v$, such
that:

$v^{T}\cdot\mathcal{H}_{\phi_{k}}\cdot v=v^{T}\cdot
J^{T}\cdot\mathcal{H}_{\phi_{k}}\cdot J\cdot v=0.$

We will use as a test direction $v=\hat{h}_{d}$. Then:

$\hat{h}_{d}^{T}\cdot\mathcal{H}_{\phi_{k}}\cdot\hat{h}_{d}=\hat{h}_{d}^{T}\cdot
Q_{d}\cdot\Lambda_{d}\cdot
Q_{d}^{T}\cdot\hat{h}_{d}^{T}-\frac{1}{k}\sum\limits_{i=1}^{M}\hat{h}_{d}^{T}\cdot
Q_{i}\cdot\Lambda_{i}\cdot Q_{i}^{T}\cdot\hat{h}_{d}^{T}$

expanding:

$\hat{h}_{d}^{T}\cdot\mathcal{H}_{\phi_{k}}\cdot\hat{h}_{d}=-\frac{2}{\left\|h_{d}\right\|^{2}}-\frac{1}{k}\sum\limits_{i=1}^{M}\left[c_{id}\;\vdots\;s_{id}\right]\cdot\Lambda_{i}\left[c_{id}\;\vdots\;s_{id}\right]^{T}$

where:

$c_{id}=\cos(\theta_{id})=\hat{h}_{d}^{T}\cdot\hat{h}_{i}$

and

$s_{id}=\sin(\theta_{id})=\hat{h}_{d}^{T}\cdot
J\cdot\hat{h}_{i}=-\hat{h}_{d}^{T}\cdot J^{T}\cdot\hat{h}_{i}$

where $\theta_{id}$ is the angle between $h_{i}$ and $\hat{h}_{d}$.

We can now write:

$\hat{h}_{d}^{T}\cdot\mathcal{H}_{\phi_{k}}\cdot\hat{h}_{d}=-\frac{2}{k}\lambda$

and following similar analysis for the perpendicular test direction:

$\hat{h}_{d}^{T}\cdot J^{T}\cdot\mathcal{H}_{\phi_{k}}\cdot
J\cdot\hat{h}_{d}=\frac{2}{k}\lambda.$

where:

$\lambda=\frac{k}{\left\|h_{d}\right\|^{2}}-\sum\limits_{i=1}^{M}\frac{c_{id}^{2}}{\left\|h_{i}\right\|^{2}}+\sum\limits_{i=1}^{M}\frac{s_{id}^{2}}{\left\|h_{i}\right\|^{2}}.$
(8)

So the problem now becomes one of identifying at least one example of a
symmetry that can cause $\lambda$ to vanish. One such example is the
following: Assume $M=2\mu+1$ obstacles where $\mu\in\mathbb{Z}_{+}$, placed
axisymmetrically and $h$ is placed such that $\theta_{id}=135^{\circ}$,
$i\in\left\\{2,\ldots\mu+1\right\\}$, $\theta_{jd}=-135^{\circ}$,
$j\in\left\\{\mu+2,\ldots 2\mu+1\right\\}$ and $\theta_{1d}=0^{\circ}$. Then
the sum
$-\sum\limits_{i=2}^{M}\frac{c_{id}^{2}}{\left\|h_{i}\right\|^{2}}+\sum\limits_{i=2}^{M}\frac{s_{id}^{2}}{\left\|h_{i}\right\|^{2}}=0$
and

$\lambda=\frac{k}{\left\|h_{d}\right\|^{2}}-\frac{1}{\left\|h_{1}\right\|^{2}}$

Requiring this to be zero we get:

$\frac{1}{\left\|h_{d}\right\|^{2}}=\frac{1}{k\left\|h_{1}\right\|^{2}}$ (9)

Now that we have identified a symmetry, we need to investigate whether this
symmetry can render the point $h$ under investigation to be a critical point.
To this end, we will project $\nabla\phi_{k}(h)$ across two perpendicular
directions, $\hat{h}_{d}$ and $\hat{h}_{d}^{\bot}=J\cdot\hat{h}_{d}$ and
investigate whether the projection is zero for both. Define the sets
$R=\left\\{2,\ldots\mu+1\right\\}$ and $L=\left\\{\mu+2,\ldots M\right\\}$. We
have that:

$\hat{h}_{d}^{T}\cdot
J^{T}\cdot\nabla\phi_{k}(h)=\frac{2}{k\sqrt{2}}\left(\sum\limits_{i\in
R}\frac{1}{\left\|h_{i}\right\|}-\sum\limits_{i\in
L}\frac{1}{\left\|h_{i}\right\|}\right)$

where we have substituted $s_{id}=-s_{jd}$ for $i\in R$ and $j\in L$. Since
$h_{i}$ and $h_{j}$ for $i\in R$ and $j\in L$ are symmetrically placed the
sums in the above equation will cancel each other and $\hat{h}_{d}^{T}\cdot
J^{T}\cdot\nabla\phi_{k}(h)=0$. Now taking the projection across $h_{d}$, we
have:

$\hat{h}_{d}^{T}\cdot\nabla\phi_{k}(h)=\frac{2}{k}\left(2\sqrt{k}-\sum\limits_{i=2}^{M}\frac{1}{\sqrt{2}}\frac{1}{\left\|h_{i}\right\|}\right)$
(10)

where we have made use of Eq. (9) and that $c_{id}=-\frac{1}{\sqrt{2}}$ for
$i\in\left\\{2,\ldots M\right\\}$. As can be seen from Eq. (10) for every
choice of $k$ there is an appropriate choice of obstacle positions in the
symmetry such that $\hat{h}_{d}^{T}\cdot\nabla\phi_{k}(h):=0$. Hence
$h\in\mathcal{C}(\phi_{k})$ and the selected configuration will be a
degenerate critical point. ∎

This tells us that there is no universal $k$ that depends only upon abstract
information about the workspace such as the number of obstacles used in Lemma
2 that would enable us to establish on-the-fly that for a workspace with a
given number of obstacles, $\phi_{k}$ will be Morse. Hence, ensuring on-the-
fly that $\phi_{k}$ has the Morse property seems unrealistic. However, the
following proposition and two lemmas establish the ”nice” behaviour of
$\phi_{k}$ in $\mathcal{P}$.

###### Proposition 2

The critical points of $\phi_{k}$ are isolated

###### Proof:

Since the logarithmic function is an analytic function, the sum of logarithmic
functions will also be analytic. Hence $\phi_{k}$ is an analytic function.
This implies that its Taylor expansion converges on any open set of its
domain. Let $E$ be an open set containing a critical point $h_{c}$ of
$\phi_{k}$. Assume that the critical point is not isolated. This implies that
there is either (a) an open neighborhood $Z$ around $h_{c}$ where $\forall
h_{z}\in Z:\phi_{k}(h_{z})=\phi_{k}(h_{c})$ or (b) there is a 1-D manifold
$\mathcal{L}$, where $\forall
h_{L}\in\mathcal{L}:\phi_{k}(h_{L})=\phi_{k}(h_{c})$.

Let us now examine case (a) where ”there is an open neighborhood $Z$ around
$h_{c}$ where $\forall h_{z}\in Z:\phi_{k}(h_{z})=\phi_{k}(h_{c})$ ”. Then the
Taylor expansion of $\phi_{k}$ around $h_{c}$ will be:

$\phi_{k}(h)=\phi_{k}(h_{c})+\sum\limits_{n=1}^{\infty}\left\\{\frac{1}{n!}\sum\limits_{m=0}^{n}\left(\begin{array}[]{c}n\\\
m\end{array}\right)D_{n,m}\left.\right|_{h_{c}}h_{cx}^{n-m}h_{cy}^{m}\right\\}$

where $D_{n,m}=\frac{\partial^{n}\phi_{k}}{\partial x^{n-m}\partial y^{m}}$
and $h-h_{c}=\left[h_{cx}\;\;h_{cy}\right]^{T}$. Since
$\phi_{k}(h_{z})=\phi_{k}(h_{c}),\forall h_{z}\in Z$, this necessarily implies
that $\forall n>0,m\geq 0:D_{n,m}=0$. Since $\phi_{k}$ is analytic, for any
open connected component $G$ of the domain of $\phi_{k}$ containing $Z$ as a
subset, it will hold $\phi_{k}(h_{g})=\phi_{k}(h_{c}),\forall h_{g}\in G$,
which according to the Principle of Permanence implies that $\phi_{k}$ will be
identically constant. However this is not true and by contradiction case (a)
is not possible.

Let us now proceed to examining case (b) where ”there is a 1-D manifold
$\mathcal{L}$, where $\forall
h_{L}\in\mathcal{L}:\phi_{k}(h_{L})=\phi_{k}(h_{c})$”. We can now construct an
analytic function $\lambda:\mathbb{R}_{\geq 0}\rightarrow\mathcal{L}$ with
$\lambda(0)=h_{c}$ and such that $\lambda(\ell)=h_{\ell}$ where $\ell$ is the
length of the curve from $h_{c}$ to $h_{\ell}$. Define the function
$\phi^{L}(\cdot)\triangleq(\phi_{k}\circ\lambda)(\cdot)$. Since
$\phi^{L}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}$ is a composition of
analytic functions, it will also be analytic. However by construction
$\phi^{L}(\ell)=\phi_{k}(h_{c})$, and this will be true for its whole domain
$\ell\in\mathbb{R}_{\geq 0}$ according to the Principle of Permanence. This
implies that $\mathcal{L}$ will have an infinite length. Then $\mathcal{L}$
will have infinitely large elements. Pick one of them, lets say
$\ell_{i}\in\mathcal{L}$ such that $\left\|\ell_{i}\right\|\rightarrow\infty$.
Then we have that $\phi_{k}(\ell_{i})=\phi_{k}(h_{c})$. However from Lemma 2
we have that $\phi_{k}(\ell_{i})\rightarrow\infty$ which is a contradiction
since for a valid workspace $\phi_{k}(h_{c})$ is finite. Hence the assumed
case (b) is not possible.

Since cases (a) and (b) were shown to not be possible, then by contradiction
the assumption that the critical point is not isolated is false and
consequently all the critical points of $\phi_{k}$ are isolated. ∎

Degenerate critical points are atypical for analytical functions in the sense
that the set of Morse functions is dense and any non-Morse smooth function can
be approximated by a Morse function [2]. Nevertheless the nature of a critical
point is not clear when the Hessian becomes degenerate (in our case its rank
may become zero). The following establishes the nature of the critical points.

###### Lemma 3

The critical points of $\phi_{k}$ are saddle points.

###### Proof:

For the case of non-degenerate critical points this is trivial since from the
proof of Lemma 1 we have that $\lambda_{1}=-\lambda_{2}\neq 0$ at a critical
points. For the case of degenerate critical points, let
$h_{d}\in\mathcal{C}(\phi_{k})$ be a degenerate critical point and
$\theta_{d}=\phi_{k}(h_{d})$. Then by the Mean Value Property of Harmonic
Functions [3] we have that:

$\theta_{d}=\int\limits_{S}\phi_{k}(h_{d}+r\zeta)d\sigma(\zeta)$

where $S$ is the boundary of a ball centered at $\theta_{d}$ with radius $r$
and $\sigma$ the normalized perimeter measure. Since $h_{d}$ is isolated,
$\phi_{k}$ on $S$ cannot be identically $\theta_{d}$. In particular for the
mean value on $S$ to be $\theta_{d}$, $\phi_{k}$ needs to take both higher and
lower values than $\theta_{d}$ which implies that $h_{d}$ is a saddle point. ∎

The next Lemma establishes the nature of the basins of attraction of the
critical points:

###### Lemma 4

The basin of attraction of the saddle points of $\phi_{k}$ is a set of measure
zero.

###### Proof:

Since $\Phi(\cdot)$ is a diffeomorphism, it is sufficient to prove this for
$\phi_{k}$. Assume a basin of attraction $\Omega$ of the critical point
$h_{c}\in\mathcal{C}(\phi)$ such that $\mu(\Omega)>0$ and
$h_{c}\in\partial\Omega$, where $\mu(\cdot)$ is the Lebesque measure of a set.
Then $\Omega$ is a positively invariant set with respect to the negated
gradient flows of $\phi_{k}$, which will be pointing inward along
$\partial\Omega$, i.e. $\nabla\phi_{k}(h_{\omega})\cdot\hat{n}\geq 0$,
$\forall h_{\omega}\in\partial\Omega$, where $\hat{n}$ the outward pointing
normal to set’s boundary. Take $\Omega_{n}\subseteq\Omega$ to be the largest
subset of $\Omega$ that contains $h_{c}$ and for which
$\nabla\phi_{k}(h_{\omega})\cdot\hat{n}>0$, $\forall
h_{\omega}\in\partial\Omega_{n}\setminus h_{c}$. Such a set is guaranteed to
exist as long as $\phi_{k}$ is not flat in $\Omega$. Hence

$\oint_{\partial\Omega_{n}}\nabla\phi_{k}\cdot\hat{n}\,dl>0$ (11)

Using the 2-dimensional Divergence Theorem we have that:

$\oint_{\partial\Omega_{n}}\nabla\phi_{k}\cdot\hat{n}\,dl=\int_{\Omega_{n}}\nabla\cdot\nabla\phi_{k}\,dA$
(12)

where $dA$ the area differential. Since $\phi_{k}$ is a Harmonic function,
then by definition $\nabla\cdot\nabla\phi_{k}=\Delta\phi_{k}=0$ and
consequently $\oint_{\partial\Omega_{n}}\nabla\phi_{k}\cdot\hat{n}\,dl=0$.
This is in contradiction with eq. (11) and since $\phi_{k}$ in $\Omega$ is
non-flat, $\Omega$ reduces to a set of measure zero. ∎

###### Theorem 1

Given a smooth navigation transformation $\Phi(\cdot)$, the function
$\varphi_{k}(x)=\sigma\circ\phi_{k}\circ\Phi(x)$ is a navigation function on
$\mathcal{W}$ as per Definition 4 for $k>M$.

###### Proof:

Let us examine the properties of Definition 4.

#### Property 1

The function $\varphi_{k}(\cdot)$ is smooth since it is a composition of
smooth functions.

#### Property 2

From Lemma 1 we have that $\phi_{k}$ is free of local minima. From (5), using
the chain rule, we have that
$\nabla\varphi_{k}=(\frac{d\sigma}{d\phi}\circ\phi_{k}\circ\Phi)J^{T}_{\Phi}(\nabla\phi\circ\Phi)$,
where $J_{\Phi}$ is the Jacobian of $\Phi$. Since $\Phi$ is a diffeomorphism,
$J_{\Phi}$ is non-singular and since $\sigma$ is strictly increasing, away
from the destination $x_{d}$ the set of critical points in $\mathcal{W}$ and
the corresponding point world match,
$\mathcal{C}(\phi_{k})=\Phi(\mathcal{C}(\varphi_{k})-\\{x_{d}\\})$. Also at
the critical points away from $x_{d}$ we have that:
$\nabla^{2}\varphi_{k_{\left|\mathcal{C}_{\phi_{k}}\right.}}=(\frac{d\sigma}{d\phi}\circ\phi_{k}\circ\Phi)J^{T}_{\Phi}(\nabla^{2}\phi\circ\Phi)J_{\Phi}$.
Hence, using the same arguments as in [6], away from $x_{d}$ we have
$index(\varphi)_{\left|\mathcal{C}_{\varphi_{k}}\right.}=index(\phi)_{\left|\mathcal{C}_{\phi_{k}}\right.}$.
Since $\phi_{k}$ has no local minima in $\mathcal{P}$, the same will be true
for $\varphi$ in $\mathcal{W}$. Now let us examine $x_{d}$. From Eq. (5) and
using the notation $h=\Phi(x)$ we have that:

$\varphi_{k}(h)=\frac{\left\|h-P_{d}\right\|^{2}}{\left\|h-P_{d}\right\|^{2}+\prod\limits_{i=1}^{M}\left\|h-P_{i}\right\|^{\frac{2}{k}}}$

Hence $\nabla\varphi_{k}(h)=\frac{2(h-P_{d})D-\left\|h-P_{d}\right\|^{2}\nabla
D}{D^{2}}$ where
$D=\left\|h-P_{d}\right\|^{2}+\prod\limits_{i=1}^{M}\left\|h-P_{i}\right\|^{\frac{2}{k}}$.
The Hessian of $\varphi_{k}$ at $x_{d}$ is thus

$\mathcal{H}_{\varphi_{k}}(x_{d})=\frac{2}{\prod\limits_{i=1}^{M}\left\|P_{d}-P_{i}\right\|^{\frac{4}{k}}}I.$
(13)

The destination $x_{d}$ is thus the unique minimum of $\varphi_{k}$ in
$\mathcal{W}$.

#### Property 3

By Proposition 2 the critical points of $\phi_{k}$ are isolated. As discussed
in the previous,
$\mathcal{C}(\phi_{k})=\Phi(\mathcal{C}(\varphi_{k})-\\{x_{d}\\})$. Hence
$\varphi_{k}$ has no additional critical points except $x_{d}$, which was
shown to be a non-degenerate global minimum. Hence all critical points of
$\varphi_{k}$ are isolated. By Lemma 3 the isolated critical points except
$x_{d}$ are saddle points and by Lemma 4 the basin of attraction of these
critical points has zero measure. Since $\Phi(\cdot)$ is a diffeomorphism, the
above properties will also hold for $\varphi_{k}(x)$.

#### Property 4

By construction $\lim\limits_{h\rightarrow P_{i}}\phi_{k}(h)\\!=\\!\infty$ for
$i=1\ldots M$. Moreover, we have the requirement that
$\lim\limits_{x\rightarrow\mathcal{O}_{0}}\left\|\Phi(x)\right\|=+\infty$
which in combination with Lemma 2 implies that

$\lim\limits_{x\rightarrow\partial\mathcal{W}}\phi_{k}\circ\Phi(x)=+\infty.$
(14)

From Eq. (4) we have that

$\lim\limits_{x\rightarrow\infty}\sigma(x)=1.$ (15)

Combining eqs. (14) and (15)) we obtain

$\lim\limits_{x\rightarrow\partial\mathcal{W}}\varphi=1$

which implies the admissibility property. ∎

##  V Navigation Function Controller Design

###  V-A Setup

In the previous section we have shown that $\varphi_{k}$ where $k\geq M+1$ is
a Navigation Function where $M$ is the number of workspace obstacles. As the
robot moves in the environment, the number $M$ of obstacles encountered is
increasing. Let $d$ be the detection radius of the robot. Any obstacle
entering this radius is registered as a new obstacle and the number of
obstacles $M$ increases when fist detected. Let
$T=\left\\{T_{1},\ldots,T_{N}\right\\}$ be the time instants where new
obstacles are discovered during the navigation of the robot and assume that
$N$ is the total number of obstacles that are discovered during a navigation
instance. Let $n(t)$ denote the number of obstacles that were discovered up to
the time instant $t$ and define:

$\Theta(x,t)\triangleq\varphi_{k(t)}(x)$ (16)

where $k(t)>n(t)$ an integer111It is sufficient to set $k(t)=n(t)+1$. Eq.
(16), is a continuous (smooth) function in the state but discontinuous in
time, representing in every time interval $T_{n(t)}\leq t<T_{n(t)+1}$ the
Navigation Function that is applied to the discovered environment at time $t$.

To ensure that the Navigation Transformation $\Phi$ is well defined at any
time instant, we will have the following:

###### Assumption 1

Sensed obstacles’ shapes and orientations are known.

In addition we need to have an obstacle ”memory” in our system as follows:

###### Assumption 2

Obstacles become part of the known workspace after being sensed.

Considering the definition of the symmetric sensing sector $S_{r,\theta}$, it
can be shown that the minimum distance that a robot can approach to an
obstacle before detecting it, is given by:

$d_{\min}(\theta)=\left\\{\begin{array}[]{ll}\min\left\\{r\sin\frac{\theta}{2},\frac{\rho_{\min}}{\cos\frac{\theta}{2}}\right\\},&\theta<\pi\\\
r,&\pi\leq\theta\leq 2\pi\end{array}\right.$ (17)

where $\rho_{\min}$ is the minimum radius of curvature of the workspace
obstacles.

The following assumption is about knowledge of the starting neighborhood:

###### Assumption 3

For a symmetric sensing sector $S_{r,\theta}$, any obstacles that are within a
radius of $d_{\min}(\theta)$ of the initial robot configuration, are assumed
known.

###  V-B Kinematic Controller Design

Let us now consider the first-order kinematic system model of Eq. (1). Our
motivation in setting up the control law for this system, is to create a
vector field that provides kinetic energy to the system that matches (or is
proportional to) the level of the Navigation Function. We have the following
result:

###### Proposition 3

System (1) with a symmetric sensing sector $S_{r,\theta}$ under the control
law:

$u=-K\sqrt{2\Theta}\cdot\widehat{\nabla\Theta}$ (18)

with $K>0$ a tuning parameter, is globally asymptotically stable almost
everywhere as long as $r>0$ and $\theta>0$

###### Proof:

Case 1: Global Sensing, i.e. $r\rightarrow\infty$ and $\theta=2\pi$:
For the case of global sensing in an environment with $M$ obstacles we have
that $n(t)=M,\forall t$ and $\Theta(x,t)=\varphi_{k}(x)$, with $k\geq M+1$, is
a Navigation Function by construction. Choosing $V=\Theta$ as a Lyapunov
function candidate we have that:

$\dot{V}=-K\sqrt{2\Theta}\left\|\nabla\Theta(x,t)\right\|\leq 0$

The set where $\dot{V}=0$, consists only of the critical points and the
destination configuration (global minimum). Since by Proposition 1, function
$\Theta$ is a Navigation Function, it has no local minima and its critical
points are isolated saddle points with attractive basins that are sets of
measure zero. Hence by LaSalle’s Invariance Principle the system will converge
to the largest invariant set that includes the destination configuration and
the saddle points that are the $\omega$ limit set of a set of measure zero of
initial conditions. Hence $\dot{V}\overset{a.e.}{<}0$ and the control law of
Eq. (18) is globally asymptotically stable, almost everywhere.

Case 2: Local sensing i.e. $r$ is finite and $0<\theta\leq 2\pi$:
In this case, the robot senses an obstacle’s existence before the minimum
distance to the obstacle becomes $d_{\min}(\theta)$.

Let $t=T_{\nu}$ be the time instant that the robot discovered the $\nu$’th
obstacle. At this time instant the control law (18) will have a discontinuity
in time. This discontinuity will be an isolated event due to Assumption 2.
Hence the trajectories of System (1) under control law (18), will flow along
Carathéodory solutions, and following the same Lyapunov analysis as the
previous case, we get:

$\underset{t}{\rm{ess}\,\sup}\left\\{-K\sqrt{2\Theta}\left\|\nabla\Theta(x,t)\right\|\right\\}\leq
0$

which implies $\dot{V}\overset{a.e.}{<}0$ and the control law of Eq. (18) is
globally asymptotically stable, almost everywhere. ∎

###  V-C Dynamic Controller Design

Now let us consider the case of the dynamical System (2). Let

$V(t)=\mu\Theta+\frac{1}{2}m\dot{x}^{T}\dot{x}$ (19)

be the total energy of the system. In this case our system will be converting
potential energy to kinetic energy and vice-versa while dissipating energy
through an appropriate dissipation term. We have the following result:

###### Proposition 4

System (2) with global sensing, under the control law:

$f=-\mu\nabla\Theta-\lambda\dot{x}$ (20)

with $\mu,\lambda$ positive gains, is globally asymptotically stable almost
everywhere.

###### Proof:

Using the same reasoning as in the Case 1 of the proof of Proposition 3 and
choosing $V$ as a Lyapunov function candidate we have that:

$\dot{V}=\mu\nabla^{T}\Theta\cdot\dot{x}+\dot{x}^{T}\cdot(-\mu\nabla\Theta-\lambda\dot{x})=-\lambda\left\|\dot{x}\right\|^{2}\leq
0$

Invariant sets where $\dot{V}=0$, consists only of the critical points and the
destination configuration (global minimum). Since by Proposition 1, function
$\Theta$ is a Navigation Function, it has no local minima and its critical
points are isolated saddle points with attractive basins that are sets of
measure zero. Hence by LaSalle’s Invariance Principle the system will converge
to the largest invariant set that includes the destination configuration and
the saddle points that are the $\omega$ limit set of a set of measure zero of
initial conditions. Hence $\dot{V}\overset{a.e.}{<}0$ and the control law of
Eq. (20) is globally asymptotically stable, almost everywhere. ∎

The following result enable us to choose the minimum damping to avoid
oscillations at the destination configuration.

###### Corollary 1

Control law (20) with

$\lambda=\lambda_{c}=2\sqrt{2\mu
m}\prod\limits_{i=1}^{M}\left\|P_{d}-P_{i}\right\|^{\frac{-1}{k}}$ (21)

provides critical damping at the neighborhood of the destination
configuration.

###### Proof:

From Eq. (13), the Navigation Function $\mu\Theta$ will behave locally at the
destination neighborhood as an elastic potential, i.e.
$\mu\Theta\approx\frac{1}{2}k_{sp}x^{T}x$, with spring constant:

$k_{sp}=2\mu\prod\limits_{i=1}^{M}\left\|P_{d}-P_{i}\right\|^{\frac{-2}{k}}$

Since $\lambda$ acts as the damping coefficient for system (2) under control
law (20), the damping ratio of the second order system will be:

$\zeta=\frac{\lambda}{2\sqrt{mk_{sp}}}.$

Setting $\zeta=1$ for critical damping, we get the result in Eq. (21) ∎

We can now state the following result in the case of local sensing:

###### Proposition 5

System (2), with a symmetric sensing sector $S_{r,\theta}$ under the control
law:

$f=-\mu\nabla\Theta-\lambda(t)\dot{x}$ (22)

with

$\lambda(t)=\left\\{\begin{array}[]{ll}2\sqrt{2\mu
m}\prod\limits_{i=1}^{n(t)}\left\|P_{d}-P_{i}\right\|^{\frac{-1}{k}},&V<\mu\\\
\frac{m}{d_{\min(\theta)}}\left\|\dot{x}(T_{n(t)})\right\|,&V\geq\mu\end{array}\right.$
(23)

and $\mu>0,r>0$ and $0<\theta\leq 2\pi$ ensures collision avoidance
everywhere, critical damping at the neighborhood of the destination
configuration and, almost everywhere global asymptotic stability.

###### Proof:

System (2) under control law (22) can be considered as a mechanical system
with dissipation under a conservative force field. This implies that if the
system starts from rest, then collisions are not possible as long as the total
energy of the system is bounded below $\mu$. This is a direct result of the
admissibility Property of the Navigation Function.

However, it is possible that the jump in the potential level of $\Theta$,
right after a new obstacle is encountered, can cause the total energy to go
above $\mu$. In this case we need to ensure that, for the worst case scenario,
there is sufficient dissipation in the system to absorb the energy added due
to introduction of a new obstacle.

Now assume that the robot is heading towards the newly discovered obstacle. To
simplify our analysis let us not take into account the repulsive force
generated by the navigation function and assume that the robot is moving
towards the obstacle with only the dissipation force slowing it down. Then the
differential equation describing the motion of the robot is given by
$m\ddot{x}+\lambda\dot{x}=0$ where $x$ is the distance across the straight
line connecting the robot position $x(T_{n(t)})$ when the obstacle was
detected and the closest obstacle point. The solution to the differential
equation is

$x(t)=\|\dot{x}(T_{n(t)})\|\frac{m}{\lambda}\left(1-\exp^{-\frac{\lambda}{m}t}\right)$
(24)

and its time derivative is

$\dot{x}(t)=\|\dot{x}(T_{n(t)})\|\exp^{-\frac{\lambda}{m}t}$ (25)

Combining equations (24) and (25) we get

$x(t)=\frac{m}{\lambda}\left(\dot{x}(T_{n(t)})-\dot{x}(t)\right)$

Setting the final velocity to zero and requiring that the travelled distance
is $d_{\min(\theta)}$, we get requirement on the damping coefficient:
$\lambda(t)=\frac{m}{d_{\min(\theta)}}\left\|\dot{x}(T_{n(t)})\right\|$, that
ensures that there is no collision due to the potential jump caused by new
obstacle discovery.

Since there are finite obstacles in the workspace, and since no collision is
possible, $\dot{V}(t)=-\lambda(t)\|\dot{x}\|^{2}$ for every time interval
$T_{n(t)}\leq t<T_{n(t)+1}$, and following the same reasoning as in the proof
of Proposition 4, we have that $\dot{V}(t)\overset{a.e.}{<}0$ in every time
interval. Now after the last workspace obstacle is discovered, the proof
follows as the proof of Proposition 4 and we recover global asymptotic
stability, almost everywhere. ∎

The following result is useful for determining the mechanical limits of the
robot actuation system:

###### Corollary 2

For system (2) under control laws (20) or (22), it holds that:

$\left\|\dot{x}\right\|<\sqrt{\frac{2\mu}{m}}$

###### Proof:

Since system (2) under control laws (20) or (22) represents a mechanical
system with dissipation under a conservative force field the maximum velocity
will occur at the lower potential level in the case that there is no
dissipation and for a system that has started at the highest potential level.
Due to the conservation of energy and from Eq. (19) we get the result. Please
note that this represents a conservative bound. ∎

##  VI Simulation Results

The effectiveness of the methodology was verified through a set of non-trivial
simulations. A symmetric sensing sector with $\theta=60^{o}$ and $r=1m$ was
used. The Navigation Transformation was built based on [7] and [4].

In the first simulation study, the kinematic controller (18) was applied to
the system (1). The robot started at its initial position in a ”seemingly”
empty workspace. As the robot moved, new obstacles were discovered through the
symmetric sensing sector as shown in Fig. 1, starting from $O_{1}$ up to
$O_{6}$. As can be seen the proposed controller performed successfully,
maintaining the robot in the workspace, handling the occurrence of new
obstacles, avoiding collisions and stabilizing at the destination
configuration.

Figure 1: [Kinematic Controller] Trajectory of system (1) under control law
(18). Symmetric sensing sector shown at destination and every instance a new
obstacle is discovered.

In the second simulation study, the dynamic controller (22) was applied to the
system (2). The mass of the system was set as $m=1Kg$. The controller
parameters were chosen as $mu=10$, $\rho_{\min}=0.05m$ whereas the dissipation
parameter $\lambda(t)$ was chosen as per Eq. (23). The robot started at the
same position as in the first simulation study with the same workspace. The
robot moved along the trajectory as shown in Fig. 2, sequentially discovering
obstacles $O_{1}$ up to $O_{6}$. The robot trajectory in this case follows a
similar - but not identical - path as the kinematic system. We can observe
limited oscillatory behavior as the robot approaches $O_{5}$ and $O_{6}$.
However, as predicted by our analysis, the robot demonstrated critical damping
without oscillations at the destination configuration. The maximum velocity of
the robot in this case study was $\left\|\dot{x}\right\|_{\max}=0.93m/s$, much
lower than the upper bound of $4.47m/s$ provided by Corollary 2. As can be
seen the proposed controller performed successfully, maintaining the robot in
the workspace, handling the appearance of new obstacles, avoiding collisions
and stabilizing without oscillations at the destination configuration.

Figure 2: [Dynamic Controller] Trajectory of system (2) under control law
(22). Symmetric sensing sector shown at destination and every instance a new
obstacle is discovered.

##  VII Conclusions

This paper presents a correct-by-construction methodology to create Navigation
Functions, without the need for tuning. The methodology builds on the concepts
of Navigation Transformation [7], Harmonic Function Based Navigation Functions
[8] and of course the Navigation Functions [6]. Without the need for tuning,
on-the-fly addition of new obstacles to the workspace is possible, while
maintaining the convergence and stability properties. Considering issues from
sensor based robot navigation, the concept of symmetric sensing sector is
introduced, modelling typical robot sensors like the laser scanner. It is
shown that navigation is possible for any non-zero symmetric sensing sector. A
kinematic controller that creates a vector field with kinetic energy that
matches the level of the Navigation Function is proposed, in order to avoid
the slowdowns close to the saddle points. The proposed dynamic controller, in
addition to the stability and collision avoidance guarantees, provides
critical damping at the destination configuration and bounds on the maximum
velocity of the system. In addition to the analytical guarantees, simulation
results verify the performance of the system.

Future work includes extending the methodology to handle moving obstacles,
obstacle and workspace mapping, multi-agent scenarios and 3D navigation.

## ACKNOWLEDGMENT

The first author would like to acknowledge the contribution of European
Union’s Horizon 2020 research and innovation program under grant agreements
824990 (RIMA), 767642 (L4MS) and by the Cyprus Research and Innovation
Foundation grants EXCELLENCE/1216/0365 (HOD-ICCCS) and EXCELLENCE/1216/0296
(RETuNE).

## References

  * [1] Omur Arslan and Daniel E Koditschek.  Sensor-based reactive navigation in unknown convex sphere worlds.  The International Journal of Robotics Research, 38(2-3):196–223, 2019.
  * [2] M. Audin and M. Damian.  Morse Theory and Floer Homology.  Springer-Verlag, London, 2014.
  * [3] S. Axler, P. Bourdon, and W. Ramey.  Harmonic Function Theory.  Springer-Verlag New York, Inc., 2nd edition, 2001.
  * [4] N. Constantinou and S. G. Loizou.  Robot navigation on star worlds using a single-step navigation transformation.  59th IEEE Conference on Decision and Control (CDC), pages 1537–1542, 2020.
  * [5] D. V. Dimarogonas.  Sufficient conditions for decentralized potential functions based controllers using canonical vector fields.  IEEE Transactions on Automatic Control, 57(10):2621–2626, 2012\.
  * [6] D. E. Koditschek and E. Rimon.  Robot navigation functions on manifolds with boundary.  Advances Appl. Math., 11:412–442, 1990.
  * [7] S. G. Loizou.  The navigation transformation.  IEEE Transactions on Robotics, 33(6):1516–1523, Dec 2017.
  * [8] S.G. Loizou.  Closed Form Navigation Functions Based on Harmonic Potentials.  50th IEEE Conference on Decision and Control and European Control Conference, 2011.
  * [9] V. Vasilopoulos, G. Pavlakos, K. Schmeckpeper, K. Daniilidis, and D. E. Koditschek.  Reactive navigation in partially familiar planar environments using semantic perceptual feedback.  CoRR, abs/2002.08946, 2020.


# Zero-energy photoelectric effect

Sajad Azizi, Ulf Saalmann and Jan M. Rost  Max-Planck-Institut für Physik
komplexer Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany

(August 28, 2024)

###### Abstract

We predict a near-threshold (“zero energy”) peak in multi-photon ionization
for a dynamical regime where the photon frequency is large compared to the
binding energy of the electron. The peak position does not depend on the laser
frequency, but on the binding energy and the pulse duration. The effect
originates from the fact that bound-continuum dipole transitions are stronger
than continuum-continuum ones. To clearly observe this zero-energy
photoelectric effect, the spectral width of the laser pulse should be
comparable to the binding energy of the ionized orbital, and the second
ionization potential should be larger than the photon energy. This suggests
negative ions as ideal candidates for corresponding experiments.

A threshold for an observable $A$ indicates the transition of the system from
one regime to another upon change of the relevant parameter $\varepsilon$ and
provides therefore important information about the system. Often, the
observable changes near threshold at $\varepsilon_{0}$ with a certain power of
the parameter, i.e.,
$A(\varepsilon{\to}{+}\varepsilon_{0})\propto\left(\varepsilon/\varepsilon_{0}{-}1\right)^{\alpha}$.
Probably best known are thermodynamical variables near phase transitions [1],
but also quantum critical points in condensed matter [2] or
fragmentation/ionization thresholds in atomic and molecular physics [3] are
examples. In the latter case the so-called Wigner threshold law [4] for
fragmentation cross sections of particles under short-range forces is a
universal property, similar as its counterpart for long-range (Coulomb)
forces, the Wannier law [5, 6].

The Wigner law has been verified in fragmentation scenarios including multi-
photon detachment of a negative ions [7], which were a popular target for
theoretical considerations [8] regarding above-threshold ionization (ATI).
However, so far it has gone unnoticed that the combination of the Wigner
threshold behavior with _short_ intense pulses can give rise to a peculiar
zero-energy photoelectric effect (ZEPE) with a characteristic maximum, whose
position at very low photo-electron energy does not depend on the photon
energy, but on the duration of the (short) laser pulse. That short intense
pulses can lead to unusual electron dynamics has been pointed out in the
context of non-adiabatic photo-ionization [9, 10] where short pulses can even
be used for coherent control [11].

The ZEPE effect requires pulses short enough such that their spectral
bandwidth $\Delta E$ is larger than the binding energy $E_{\rm EA}$ of the
detached electron, typically of the order of a few femtoseconds. This was
clearly not the case in the ATI detachment experiments with 800 nm light. ZEPE
is enabled by a two-photon process, where a bound electron absorbs a photon
and emits a photon, ending at the same (bound) energy as it started from.
Therefore, this mechanism does typically not contribute to the photo-electron
spectrum. Only, if the laser pulse is short enough that the spectral peak at
the binding energy “leaks” into the continuum, see Fig. 1, this process
becomes visible in form of photo-electrons and can lead, together with the
Wigner power law, to a pronounced ZEPE maximum.

Negative ions are an ideal target to clearly identify the ZEPE, since they
combine a low electron affinity (EA), i.e., the ionization potential of the
negative ion, with a typically large gap to the next ionization potential
(IP), i.e., the ionization potential of the (neutral) atom. Therefore, a
relatively large photon frequency, still fulfilling $E_{\rm EA}<\omega<E_{\rm
IP}$, can be chosen, which on the one hand energetically prevents ionization
of more deeply bound electrons and on the other hand lets the first ATI peak
appear at relatively high energy, keeping an energy interval just above
threshold pristine for a clean signature of ZEPE.

Figure 1: Multi-photon detachment of a an electron, weakly bound by electron
affinity $E_{\rm EA}$ to an ion, will give rise to peaks in the photo-electron
spectrum at energies $E\,{=}\,{-}E_{\rm EA}{+}n\omega$, shown here for
$n{=}1,2$. However, following absorption of the 1st photon it is more likely
that a 2nd photon is not absorbed but emitted (${-}\hbar\omega$) since the
dipole for this transition is much stronger. This process with net-zero energy
absorption becomes only visible in the spectrum for very short pulses (blue)
leading to a broad peak at $-E_{\rm EA}$ whose tail reaches into the continuum
(lightblue-shaded), but not for longer pulses (red).

We will demonstrate the effect in the following with negative ions of hydrogen
and oxygen in a simple and transparent fashion. To this end, we use an
effective potential for the bound electron of the negative ion which is
designed to reproduce the EA well [12, 13],

$V(r)=-\frac{Z}{r}\frac{b}{1+c\,[\exp(r/r_{0}){-}1]}\,,$ (1)

where $Z$ is the nuclear charge, $Z\,{=}\,1$ and $Z\,{=}\,8$ for hydrogen and
oxygen, respectively. We use atomic units unless stated otherwise. With
parameter values $b\,{=}\,1.1$, $c\,{=}\,1$, and $r_{0}\,{=}\,0.5292$ Å the
computed $E_{\rm 1s}\,{=}\,{-}0.75$ eV for H- closely matches the experimental
value [14]. Similarly, with $b\,{=}\,1$, $c\,{=}\,1.9607$, and
$r_{0}\,{=}\,0.4689$ Å, the computed $E_{\rm 2p}\,{=}\,{-}1.464$ eV for O- is
in good agreement with the experimental value of $-1.461$ eV [15].

With the potential (1) the Hamilton operator reads

$\widehat{H}(t)=\widehat{\bm{p}}{\,}^{2}\\!/2+V(r)+\widehat{\bm{p}}\cdot\bm{e}_{z}\,A_{0}g(t)\cos(\omega
t),$ (2)

where $g(t)\,{=}\,{\exp}({-}t^{2}/T^{2})$ is the laser-pulse envelope with the
pulse duration $T\,{=}\,T_{\rm fwhm}/\sqrt{2\ln\\!2}$, linearly polarized
along the $z$-axis and dipole coupled in velocity gauge. In single-active-
electron approximation (2), the propagation of the 3D time-dependent
Schrödinger equation (TDSE), for processes discussed here, is straightforward,
since the number of photons involved is moderate. What is challenging,
however, is the numerically accurate description of the _very_ small
detachment probabilities close to threshold. To this end we calculate
eigenstates and dipole matrix elements for the potentials given in (1) for
$\ell\,{=}\,0\ldots\ell_{\rm max}$ in a box $r\,{=}\,0\ldots R$ up to a cutoff
energy $E_{\rm cut}$. The parameters used are $\ell_{\rm max}\,{=}\,4$,
$R\,{=}\,3{\times}10^{3}$ Å and $E_{\rm cut}\,{=}\,3$ keV. The TDSE is
propagated numerically with this field-free basis. From the final amplitudes
obtained we calculate the photo-electron spectrum. Attaching to every
eigenstate a normalized Gaussian, whose height is the absolute square of the
amplitude and whose width is given by the level spacing at the corresponding
eigenenergy, produces a continuous spectrum. Thereby, the density of states
and the detachment process are correctly included.

Figure 2: Spectra from H- for $T_{\rm fwhm}\,{=}\,1$ fs,
$I\,{=}\,10^{14}$W/cm2 and 4 different photon frequencies
($\omega\,{=}\,3,6,9,12$ eV). Note, while higher-order peaks (at energies
$E\,{\gtrsim}\,6$ eV) shift with increasing $\omega$ to higher energies, the
zero-peak position ($E_{\rm s}\,{\approx}\,0.3$ eV and $E_{\rm
d}\,{\approx}\,1.4$ eV) is independent of $\omega$. Spectra are scaled such
that they agree for $E\,{=}\,0.1$ eV.

Figure 2 shows the resulting photo-electron spectra for H- exposed to a laser
pulse of 1 fs duration and a peak intensity of
$I\,{=}\,10^{14}\,\mathrm{W/cm^{2}}$ for different photon frequencies. The
unique feature of the zero-energy photo-effect is evident: There is a
pronounced maximum close to threshold ($E_{\rm s}{\approx}0.3$ eV and $E_{\rm
d}{\approx}1.4$ eV) which is independent of the photon frequency, clearly
different in shape from the two-photon ATI peaks which can be identified at
the respective energies $E\,{=}\,2\omega\,{-}\,E_{\rm EA}$. If the ZEPE
maximum is much smaller than and too close to the ATI peaks, typically for
higher final partial waves $\ell>0$, then it may get buried under the rise to
the two-photon ATI peak, which is the case for H${}^{-}{\to}$ d detachment
with 3 eV photons, as visible in Fig. 2b.

Figure 3: Photo-electron spectrum (colored solid lines) of H- for a Gaussian
pulse with different pulse duration $T_{\rm fwhm}=1,2,3$ fs, carrier frequency
$\omega=9\,$eV, and intensity $I=10^{14}\,\mathrm{W/cm^{2}}$ for the s (a) and
the d-channel (b), respectively, and fits (black dashed lines) according to
Eq. (4).

To understand how ZEPE comes about qualitatively and quantitatively, we take a
closer look at the near-threshold energy range with the spectra of Fig. 3 on a
linear scale for different pulse durations $T$. One immediately sees that the
spectra depend on $T$, as already anticipated. The Wigner threshold law [4]
states that for break-up of two fragments under short-range forces (which is
the case for electron detachment), the ionization probability near threshold
is given by the available continuum states in momentum space, i.e.,

$P(E{\to}0)\propto\mbox{\large$\int$}\\!\mathrm{d}p\>p^{2(\ell+1)}\delta(p^{2}\\!/2{-}E)\propto
E^{\ell+1/2},$ (3)

where $\ell$ is the angular momentum of the fragment pair. Following the
intuition that ZEPE originates in the combination of the Wigner law and a
Gaussian energy distribution induced by the short laser pulse, located at the
binding energy $-E_{\rm EA}$ due to the two-photon zero-energy process, the
detachment probability should be given by

$\displaystyle P_{\rm ZEPE}^{\ell,\beta}(E)$
$\displaystyle=P_{*}(E)\,s_{\ell,\beta}(E/E_{\rm EA}),$ (4) $\displaystyle
P_{*}(E)$ $\displaystyle=[1+E/E_{*}]^{-1}$ (4a) $\displaystyle
s_{\ell,\beta}(x)$
$\displaystyle=\beta^{4}x^{\ell+\frac{1}{2}}\exp(-\beta^{2}[x{+}1]^{2}).$ (4b)

The universal shape $s_{\ell,\beta}(x)$ of the zero-energy photo-electron
spectrum is fully determined by the angular momentum $\ell$ of the photo-
electron and $\beta\,{=}\,E_{\rm EA}/\Delta E$ the ratio of electron affinity
$E_{\rm EA}$ to spectral pulse width $\Delta E$. From the spectral
representation of the pulse and the fact that ZEPE is a two-photon process
follows $\Delta E\,{=}\,2/T$, which is confirmed by 2nd-order perturbation
theory discussed below. Furthermore, $P_{*}(E)$ takes care of the slowly-
varying background, being different for each ion.

Figure 4: Photo-electron spectra for $H^{-}$ (a,b) as in Fig. 3 but also for
the two final channels of O- (c,d) in double-logarithmic scale. Note the
different vertical scale for each panel. The gray line in panel (c) is for a
lower intensity of $I{=}10^{12}$W/cm2 and therefore scaled by a factor of
$10^{4}$. The dotted lines in panel (d) are computed with a basis of states up
to $E_{\rm cut}{=}30$ eV which would allow to describe 3-photon processes, the
solid curves are for $E_{\rm cut}{=}3$ keV. Arrows mark the maxima according
to Eq. (5).

The dashed lines in Fig. 3 show $P^{\ell,\beta}_{\rm ZEPE}(E)$ according to
Eq. (4). The only fit parameter $E_{\ast}$ assumes the values
$E_{\ast}{=}0.896$ eV and $E_{\ast}{=}0.684$ eV for the “s” and “p” photo-
electrons of hydrogen and oxygen, respectively, while $E_{\ast}{=}13.4$ eV and
$E_{\ast}{>}1000$ eV for the electrons detached with angular momentum “d” and
“f” reveal that the background is nearly constant in these cases over the
energy interval considered. Obviously, Eq. (4) represents the numerical
results very well, giving confidence in the interpretation and description of
ZEPE provided. Deviations for O- at “longer” pulses can be attributed to
higher-order effects occurring for intense pulses, since the agreement with
the analytic description for a weaker pulse with $I{=}10^{12}$W/cm2, as shown
by the gray line in Fig. 4c, is excellent.

Peculiar at second glance is the variation of the maximal detachment with the
pulse length $T$ or $\Delta E$, respectively. Note, that the highest maximum
is achieved in Fig.3 with $T_{\rm fwhm}{=}\,2$ fs for the $\ell{=}0$ case,
while for all other $\ell$ the maxima seem to increase monotonically. To
elucidate their behavior systematically, we determine the maximum of $P_{\rm
ZEPE}^{\ell,\beta}$ with respect to the two-dimensional parameter space
$\\{\beta,x\\}$, i. e., the (scaled) pulse width $\beta\,{=}\,E_{\rm
EA}/\Delta E$ and the (scaled) excess energy $x\,{=}\,E/E_{\rm EA}$. The
corresponding optimal parameters are given by analytical but lengthy
expressions. Inserting $\beta(x)\,{=}\,\sqrt{2}/(1{+}x)$, the solution of
$\partial\\!P_{\rm ZEPE}^{\ell,\beta}/\partial\\!\beta{=}0$ into
$\partial\\!P_{\rm ZEPE}^{\ell,\beta}/\partial\\!x{=}0$ gives an analytical
but lengthy expression for $\beta_{\rm Max}(x^{*},\ell)$, i. e. the pulse
width that gives the highest peak of the detachment probability. Ignoring the
background modification (4a) by taking the limit $x^{*}{\to}\infty$ one
obtains with the simpler expressions $\beta_{\rm
Max}\,{=}\,(7{-}2\ell)/{4\sqrt{2}}$ and $x_{\rm
Max}\,{=}\,(1{+}2\ell)/(7{-}2\ell)$ directly the optimal pulse duration
$T_{\rm fwhm}{=}\sqrt{8\log 2}\,\beta_{\rm Max}/E_{\rm EA}$. They values are
2.55 fs and 1.09 fs for H- (s- and d-channel) and 0.92 fs and 0.18 fs for O-
(p- and f-channel), consistent with the data shown in Fig. 4. Within the same
approximation one can determine the location of the maxima of $P_{\rm
ZEPE}^{\ell,\beta}$ for a given pulse durations $\beta$. They read

$E_{\rm max}=\frac{E_{\rm
EA}}{2}\big{(}\sqrt{1{+}(2\ell{+}1)/\beta^{2}}-1\big{)}$ (5)

and are shown in Fig. 4 for the respective pulses with vertical arrows.
Although for $\ell{=}0$ (H-) and $\ell{=}1$ (O-) the background is relevant,
since $E_{*}$ is in the energy range of interest, the values (5) explain the
maxima there quite well.

Hence, the overview of ZEPE for H- and O- in Fig. 4, highlighting the
influence of different partial waves $\ell=0,...,3$ of the fragments reveal
that essentially all features discussed so far can be identified and are
confirmed. For completeness we estimate the number of ZEPE electrons per shot
one could detect in an experiment,

$N_{\rm exp}=N_{\rm ion}\big{[}I_{\rm exp}\big{/}10^{14}{\rm
W/cm}^{2}\big{]}^{2}P_{T},$ (6)

where $N_{\rm ion}$ is the number of ions in the target volume, $I_{\rm exp}$
the experimental laser intensity and $P_{T}$ is given in Tab. 1.

Table 1: Total ionization probability $P_{T}$ for H- and O- for accessible channels and three pulse durations $T_{\rm fwhm}$. |  | H${}^{-}{\to}\,$s |  | H${}^{-}{\to}\,$d |  | O${}^{-}{\to}\,$s |  | O${}^{-}{\to}\,$d
---|---|---|---|---|---|---|---|---
$T_{\rm fwhm}$=1fs |  | $8.9{\times}10^{-4}$ |  | $1.5{\times}10^{-6}$ |  | $1.8{\times}10^{-4}$ |  | $7.1{\times}10^{-8}$
2fs |  | $4.4{\times}10^{-4}$ |  | $9.1{\times}10^{-8}$ |  | $3.1{\times}10^{-6}$ |  | $9.7{\times}10^{-11}$
3fs |  | $1.1{\times}10^{-4}$ |  | $6.4{\times}10^{-9}$ |  | $3.9{\times}10^{-8}$ |  | $6.8{\times}10^{-13}$

The double-logarithmic scale in Fig. 4 was chosen to emphasize the threshold
behavior and one can see that the quasi-analytical formula (4) performs in all
cases extremely well compared to the numerical TDSE calculations, particularly
with regard to the large dynamic ranges considered. In fact, the TDSE results
in panel 4d contain a surprise: They appear not to be converged for the
$\ell\,{=}\,3$ spectrum of oxygen. Indeed, despite very small photo-electron
energies of the order of 1 meV, numerically one has to include continuum
electrons up to 3 keV to achieve convergence (see also Fig. 5). An explanation
is provided by 2nd-order perturbation theory for the two-photon zero-energy
process comprising the two different events of absorption followed by emission
($\eta{=}``+"$), typically much stronger than emission followed by absorption
($\eta{=}``-"$). The ionization amplitude to energy $E$ for these processes
reads [16]

$\displaystyle a(E)$ $\displaystyle=\sum_{k}d_{Ek}d_{k\,\rm
EA}\int_{-\infty}^{+\infty}\\!\\!\\!\mathrm{d}t\,A(t)\mathrm{e}^{{\mathrm{i}[E-E_{k}]t}}\int_{-\infty}^{t}\\!\\!\\!\mathrm{d}t^{\prime}\,A(t^{\prime})\mathrm{e}^{{\mathrm{i}[E_{k}+E_{\rm
EA}]t^{\prime}}}$
$\displaystyle=\frac{\pi}{8}T^{2}A_{0}{\\!}^{2}\sum_{k}d_{Ek}d_{k\,\rm
EA}\,\sum_{\eta=\pm}\Big{[}\mathrm{e}^{{-[[\Delta^{\eta}_{k}+E_{\rm
EA}]^{2}+[E-\Delta^{\eta}_{k}]^{2}]T^{2}\\!/4}}-\tfrac{2\mathrm{i}}{\sqrt{\pi}}\mathrm{e}^{{-[E+E_{\rm
EA}]^{2}T^{2}\\!/8}}\,F\big{(}[2\Delta^{\eta}_{k}{+}E_{\rm
EA}{-}E]T/\sqrt{8}\big{)}\Big{]}$ (7)

with the dipole matrix elements
$d_{jk}\,{\equiv}\,\bm{e}_{z}\,{\cdot}\,\langle\phi_{j}|\bm{d}|\phi_{k}\rangle$
and the detunings $\Delta^{\pm}_{k}\,{\equiv}\,E_{k}{\mp}\omega$. We have not
explicitly specified that $|\phi_{0}\rangle$, $|\phi_{k}\rangle$,
$|\phi_{E}\rangle$ have different angular momenta $\ell$, dictated by
selection rules in the dipole matrix elements, and for simplicity we have
assumed that all virtual states $|\phi_{k}\rangle$ are discrete in accordance
with our numerical treatment of the continuum. Finally, due to the necessary
cutoff at $E_{\rm cut}$, the sum over virtual states $k$ is finite. From the
real part of (7) one sees that continuum states, which are resonant with
initial one-photon absorption $\Delta_{k}\,{=}\,0$, are dominant as well as
the final energy equal to the initial energy, $E\,{=}\,{-}E_{\rm EA}$, as
expected for the zero-energy photoelectric effect.

Figure 5: Photo-electron spectra in linear (left panel) and double-logarithmic
(right) scale for the f-channel of O- from TDSE calculations ($\omega$ and $I$
as in Fig. 3, $T_{\rm fwhm}{=}1$ fs) with different cut-off energies $E_{\rm
cut}$, see text.

Coming back to the requirement of including very high energies $E_{k}$ in the
calculation to reach convergence for small finite energy $E$, we note that the
Dawson function $F$ [17] in the imaginary part of Eq. (7) falls off very
slowly $F(x)\,{\sim}\,1/(2x)$ for large $x$. That implies that non-resonant
states $E_{k}\gg E$ can contribute significantly. Since the detachment
probability ${\propto}\,(E/E_{\rm EA})^{\ell+1/2}$ is very small near
threshold, it must be determined with high absolute accuracy which explains
why unexpectedly high virtual energies need to be taken into account. Note,
that for the very same reason the two-photon process of first emitting and
then absorbing a photon, in standard situations never considered, has a non-
negligible contribution.

In summary, we have discussed a universal two-photon process optimally
observable in negative ions, where the weakly-bound electron interacts with a
short laser pulse of spectral width comparable to the electronic binding
energy. Together with the typical Wigner power law of the electron-detachment
cross section near threshold, a characteristic maximum forms at very-low
electron energies, which does not dependent on the photon frequency, as
accurate numerical calculations confirm. They require the inclusion of
surprisingly energetic electronic continuum states to converge, highlighting
subtle features of this unconventional process which also surface in 2nd-order
time-dependent perturbation theory. Following physical intuition, we have
derived an analytical form of the spectrum. It describes the relative peak
height and location as well as the dependence on pulse length and electron
affinity of this peculiar zero-energy photoelectric effect very well which
will facilitate its experimental realization.

S. A. acknowledges discussion with Jonathan Dubois in the early stage of this
work.

## References

  * [1] H.E. Stanley,  Introduction to phase transitions and critical phenomena.  International series of monographs on physics. Oxford University Press 1971.
  * [2] S. Sachdev,  Quantum phase transitions.  Cambridge University Press 2000.
  * [3] J. M. Rost, _Critical phenomena in atomic physics_ Physica E 9, 467 (2001).
  * [4] E. P. Wigner, _On the behavior of cross sections near thresholds_ Phys. Rev. 73, 1002 (1948).
  * [5] G. H. Wannier, _The threshold law for single ionization of atoms or ions by electrons_ Phys. Rev. 90, 817 (1953).
  * [6] J. M. Rost, _Two-electron escape near threshold: A classical process?_ Phys. Rev. Lett. 72, 1998 (1994).
  * [7] R. Reichle, H. Helm, and I Yu. Kiyan, _Photodetachment of H - in a strong infrared laser field_ Phys. Rev. Lett. 87, 243001 (2001).
  * [8] K. Krajewska, I. I. Fabrikant, and A. F. Starace, _Threshold effects in strong-field detachment of H - and F-: Plateau enhancements and angular distribution variations_ Phys. Rev. A 74, 053407 (2006).
  * [9] K. Toyota, O. I. Tolstikhin, T. Morishita, and S. Watanabe, _Slow electrons generated by intense high-frequency laser pulses_ Phys. Rev. Lett. 103, 153003 (2009).
  * [10] Q.-C. Ning, U. Saalmann, and J. M. Rost, _Electron dynamics driven by light-pulse derivatives_ Phys. Rev. Lett. 120, 033203 (2018).
  * [11] S. Azizi, U. Saalmann, and Jan M. Rost, _Non-adiabatic ionization with tailored laser pulses_ J. Phys. B 54, 134001 (2021).
  * [12] P. S. Ganas, J. D. Talman, and A. E. S. Green, _Independent-particle models for light negative atomic ions_ Phys. Rev. A 22, 336 (1980).
  * [13] R. Shakeshaft and X. Tang, _Integral-equation approach to multiphoton ionization by intense fields. II. Application to H and H -_ Phys. Rev. A 36, 3193 (1987).
  * [14] K. R. Lykke, K. K. Murray, and W. C. Lineberger, _Threshold photodetachment of H -_ Phys. Rev. A 43, 6104 (1991).
  * [15] M. K. Kristiansson, K. Chartkunchand, G. Eklund, O. M. Hole, E. K. Anderson, N. de Ruette, M. Kamińska, N. Punnakayathil, J. E. Navarro-Navarrete, St. Sigurdsson, J. Grumer, A. Simonsson, M. Björkhage, St. Rosén, P. Reinhed, M. Blom, A. Källberg, J. D. Alexander, H. Cederquist, H. Zettergren, H. T. Schmidt, and D. Hanstorp, _High-precision electron affinity of oxygen_ Nat. Commun. 13, 5906 (2022).
  * [16] S. Azizi,  Three aspects of photo-ionization in ultrashort pulses PhD thesis, Technische Universität Dresden 2023.
  * [17] The Dawson function $F(z)$ [18] is related to the (complex) error function, $F(z)=\frac{\sqrt{\pi}}{2\mathrm{i}}\mathrm{e}^{{-z^{2}}}\mbox{erf}(\mathrm{i}z)$.
  * [18] H. G. Dawson, _On the numerical value of $\int_{0}^{\infty}\mathrm{d}x\,\mathrm{e}^{x^{2}}$_ Proc. London Math. Soc. XXIX, 519 (1897).