# Doppler Signature in Electrodynamic Retarded Potentials

Giovanni Perosa1,2, Simone Di Mitri1,2, William A Barletta3 and Fulvio
Parmigiani1,2,4<EMAIL_ADDRESS>1 Dipartimento di Fisica,
Università di Trieste, via A. Valerio 2, 34127 Trieste, Italy  2 Elettra
Sincrotrone Trieste, I-34149 Basovizza, Trieste, Italy  3 Dept. of Physics,
Massachusetts Institute of Technology Bldg. 24-402, 77 Massachusetts Avenue
Cambridge, MA 02139  4 International Faculty, University of Cologne, Albertus-
Magnus-Platz, 50923 Cologne, Germany

###### Abstract

The presence of the term $\left(\frac{v}{c}\right)$ that characterizes the
electrodynamic retarded potentials, also known as Liénard-Wiechert potentials,
is thought to be reminiscent of a Doppler effect. Here, we show that these
potentials are consistent with the Doppler shifted electromagnetic (e.m.)
field generated by a charge moving along a generic trajectory with respect to
the laboratory reference frame. Of course, the retarded potentials derived
here are formally the same as those reported in the current literature.
Nonetheless, this work sheds a new light on the origin of the electrodynamics
retarded fields, while offering a direct physical interpretation of the term
$\left(\frac{v}{c}\right)$ characterizing the related potentials.

##  I Introduction

The retarded potentials, also known as Liénard-Wiechert potentials, here after
L-W potentials, are the potentials associated with a charge moving along a
generic trajectory, hence of paramount importance in classical and quantum
electrodynamics. Griffith
Assuming that an electromagnetic (e.m.) signal generated by an accelerated
charge is traveling at a constant speed, $c$, it will reach a point
$\mathbf{r}$, at rest in the laboratory frame, at a later time, $t$, with
respect to the starting signal time $t_{r}$ (see figure. 1). That is why the
potentials related to the e.m. field, $\psi(\mathbf{r},t)$ and
$\mathbf{A}(\mathbf{r},t)$, are called ”retarded potentials”. Interestingly,
these potentials depend on the normalized velocity of the charge,
$\left(\frac{v}{c}\right)\equiv[\beta\mathbf{]}$, at the retarded time
($t_{r}$) and retarded distance respect to the time $t$ and distance
$\mathbf{r}$ at which the signal is detected, see point in $P$ in figure. 1.

Figure 1: Point charge $q$ moving in the reference frame $S$, following the
trajectory $\mathbf{r}_{O}(t)$. The retarded potential at $P$ is computed as
the superposition of outgoing spherical wave coming from $q$, schematically
represented as circular wavefronts in red. The relative motion of the charge
with respect to $P$ is the origin of the Doppler shift of the frequencies.

Most textbooks Jackson ; Stratton ; Zangwill ; Feynman ; Panofsky ; Orahilly ;
Reitz ; Tommasini  derive the L-W potentials following the original works of
Liénard Lienard  and Wiechert Wiechert , which are based on the superposition
principle of the fields produced by different volume elements of a charge
particle assumed to have a finite size volume. Otherwise, L-W potentials are
derived via Lorentz transformations, following the hints first given by
Minkowski.Minkowski ; Padmanabhan ; Streltsov  Sommerfeld also provided two
derivations of the retarded potential by starting from the Green’s function.
Haus ; Frank
After Sommerfeld, other authors have derived the L-W potentials using the
Green function formalism, see for example the text books of J. Schwinger et
al. Schwinger , and Zangwill Zangwill .
The Green’s functions formalism is a powerful and elegant mathematical tool,
whenever applicable, to describe a physical phenomenon that evolves from an
initial state to a final state. However, it leaves the physics hidden, whether
it is not properly disclosed.
This work has been inspired by a note reported in the widespread textbook of
D. J. Griffith where the author, referring to the $\left(\frac{v}{c}\right)$
term, writes, The argument is somewhat reminiscent of the Doppler effect.G2
Hence, our aim is to derive the L-W potentials starting from the Doppler
effect associated to an e.m. signal generated by an accelerated charge as
detected in the laboratory reference frame.
In this context, it is remarkable to note that Haus Haus  derived the e.m.
field generated by a single charge moving along a generic trajectory by
considering the Fourier-transformed components of the associated current
density. He showed that at constant velocity $v<c$, the charge does not
radiate because the light-like current density Fourier components are absent.
However, while this work helps to elucidate the close relationship between
current density and the radiation associated with a moving charge, thus paving
the way for an elegant interpretation of the Cherenkov effect, it does not
shed light on the Doppler mechanism characterizing the delayed potentials.
In this scenario, the difficulty in understanding the nature of the L-W
potentials is related to the physical origin of the Doppler-like factor
$\left(\frac{v}{c}\right)$ that characterises their mathematical
representation.

##  II Derivation of the retarded potentials

###  II.1 Retarded Green’s function

Let’s start from the retarded potentials in the Lorenz gauge Griffith ;
Jackson

$\displaystyle\left(\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}-\nabla^{2}\right)\psi(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{\rho(\mathbf{r},t)}{\epsilon_{0}}$ (1)
$\displaystyle\left(\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}-\nabla^{2}\right)\mathbf{A}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\mu_{0}\mathbf{J}(\mathbf{r},t)$ (2)

where $\psi(\mathbf{r},t)$ is the electric scalar potential,
$\mathbf{A}(\mathbf{r},t)$ is the magnetic vector potential,
$\rho(\mathbf{r},t)$ is the electric charge density,
$\mathbf{J}(\mathbf{r},t)$ is the current density, and $\nabla^{2}$ is the
Laplacian ($\nabla^{2}=\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial y^{2}}+\frac{\partial^{2}}{\partial
z^{2}}$). These differential equations are often referred to as wave
equations. Hassani  The inhomogeneous part, i.e. the non trivial right-hand
side (RHS), is regarded as the source the electromagnetic waves. In a matter
free space, $c$ is the e.m. wave phase velocity. Thereafter, we focus on the
case of e.m. fields generated by a moving charge $q$, having charge and
current density,

$\rho(\mathbf{r},t)=q\delta(\mathbf{r}-\mathbf{r}_{O}(t))\quad\text{a}nd\quad{\mathbf{J}(\mathbf{r},t)=q\mathbf{v}(t)\delta(\mathbf{r}-\mathbf{r}_{O}(t))},$
(3)

respectively, being $\mathbf{r}_{O}(t)$ the charge trajectory in the
laboratory reference frame (as shown in figure. (1)) and $\mathbf{v}(t)$ its
velocity. These equations can be solved by starting from the retarded Green’s
function for the wave equation Zangwill ; Hassani

$G_{-}(\mathbf{r},t)=\frac{1}{4\pi r}\delta(t-r/c).$ (4)

Assuming that the observer is affected only by the presence of a source that
acts at an earlier time, we can discard the advanced solution, proportional to
$(t+r/c)$, and use directly the retarded one.
The Fourier expansion of $G_{-}(\mathbf{r},t)$ in the frequency domain
describes a spherical wave propagating outward from the e.m. field
sourceZangwill ,

$G_{-}(\mathbf{r},t)=\frac{1}{2\pi}\int
d\omega\tilde{G}(\mathbf{r},\omega)e^{i\omega
t},\quad\text{with}\quad\tilde{G}(\mathbf{r},\omega)=\frac{e^{ikr}}{4\pi r}.$
(5)

###  II.2 Consequences of the relative motion

Figure 1 shows that, since the signal has a finite velocity ($c=\omega/k$),
the information carried by an e.m. wave will take some time to reach the point
$P$. The time at which the wavefront registered in $P$ at the time $t$ was
originated is called retarded time $t_{r}$. This can be retrieved by the
identity of the wavefront phase in $P$ at time $t$ and in the origin at time
$t_{r}$, where the source particle’s trajectory is properly taken into
account. Let us indicate with $\phi(\mathbf{r},t)$, the phase of a spherical
wave

$\phi(\mathbf{r},t)=kr-\omega t,$ (6)

and the condition for the retarded time translates into

$\phi(\mathbf{r},t)=\phi(\mathbf{r}_{O}(t_{r}),t_{r}).$ (7)

Inserting the expression for the phase of a spherical phase (6), we obtain:

$\phi(\mathbf{r},t)-\phi(\mathbf{r}_{O}(t_{r}),t_{r})=k|\mathbf{r}-\mathbf{r}_{O}(t_{r})|-\omega(t-t_{r})=-\omega\left(t-t_{r}-\frac{|\mathbf{r}-\mathbf{r}_{O}(t_{r})|}{c}\right)=0,$

hence,

$t_{r}=t-\frac{|\mathbf{r}-\mathbf{r}_{O}(t_{r})|}{c}.$ (8)

The second relevant consequence of the motion of the source in the reference
frame $S$ is the Doppler shift of the wave frequency

$\displaystyle\omega^{\prime}=-\frac{d\phi}{dt}$ $\displaystyle=$
$\displaystyle-\mathbf{k}\cdot\frac{d\mathbf{r}}{dt}+\omega=\omega\left(1-\frac{\hat{\mathbf{k}}\cdot\mathbf{v}}{c}\right)=\frac{\omega}{D_{\hat{\mathbf{k}}}(\mathbf{v})}$
(9) $\displaystyle D_{\hat{\mathbf{k}}}(\mathbf{v})$ $\displaystyle=$
$\displaystyle\left(\frac{c}{c-\hat{\mathbf{k}}\cdot\mathbf{v}}\right)\quad\text{and}\quad\hat{\mathbf{k}}=\hat{\mathbf{r}}$
(10)

Equation (10) can also be interpreted as the change of the phase velocity due
to the relative motion of the source in the laboratory reference frame, in a
consistent analogy with the apparent pitch change of the sound coming from a
moving emitter.
A single reference frame has been considered so far, allowing us to avoid any
reference to special relativity and relativistic Doppler effect.

###  II.3 Computation of the retarded potentials

According to the superposition principle, the solution to equation (1) is
given by the convolution of the Green’s function and the point charge
distribution

$\psi(\mathbf{r},t)=\frac{1}{\epsilon_{0}}\int dt^{\prime}\int
d\tau^{\prime}G_{-}(\mathbf{r},t|\mathbf{r}^{\prime},t^{\prime})\rho(\mathbf{r^{\prime}},t^{\prime}),$
(11)

where $d\tau^{\prime}$ is the volume three dimensional element and

$G_{-}(\mathbf{r},t|\mathbf{r}^{\prime},t^{\prime})=\frac{1}{4\pi|\mathbf{r}-\mathbf{r}^{\prime}|}\delta(t-t^{\prime}-|\mathbf{r}-\mathbf{r}^{\prime}|/c).$
(12)

In view of the interpretation that we give of the Green’s function in the
frequency domain, equation (11) represents the superposition of outgoing
spherical waves.
We can substitute the expression (5) and the first of equations (3) in (11) to
get

$\psi(\mathbf{r},t)=\frac{q}{2\pi\epsilon_{0}}\int dt^{\prime}\int
d\tau^{\prime}\delta(\mathbf{r^{\prime}}-\mathbf{r}_{O}(t^{\prime}))\int
d\omega\frac{1}{4\pi|\mathbf{r}-\mathbf{r}^{\prime}|}\exp{\big{\\{}i\left[k|\mathbf{r}-\mathbf{r}^{\prime}|-\omega(t-t^{\prime})\right]\big{\\}}}.$
(13)

The first step is to simplify the inner couple of integrals (the one over
$\tau^{\prime}$ and the one over $\omega$):

$\displaystyle\psi(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{q}{2\pi\epsilon_{0}}\int dt^{\prime}\int
d\tau^{\prime}\delta(\mathbf{r}^{\prime}-\mathbf{r}_{O}(t^{\prime}))\int
d\omega\frac{1}{4\pi|\mathbf{r}-\mathbf{r}^{\prime}|}\exp{\big{\\{}i\left[k|\mathbf{r}-\mathbf{r}^{\prime}|-\omega(t-t^{\prime})\right]\big{\\}}}=$
(14) $\displaystyle=$ $\displaystyle\frac{q}{2\pi\epsilon_{0}}\int
dt^{\prime}\int
d\omega\frac{1}{4\pi|\mathbf{r}-\mathbf{r}_{O}(t^{\prime})|}\exp{\big{\\{}i\left[k|\mathbf{r}-\mathbf{r}_{O}(t^{\prime})|-\omega(t-t^{\prime})\right]\big{\\}}}=$
$\displaystyle=$ $\displaystyle\frac{q}{2\pi\epsilon_{0}}\int dt^{\prime}\int
d\omega\frac{1}{4\pi|\mathbf{r}-\mathbf{r}_{O}(t^{\prime})|}\exp{\big{\\{}i[\phi(\mathbf{r},t)-\phi(\mathbf{r}_{O}(t^{\prime}),t^{\prime})]\big{\\}}}$

We expand the argument of the exponential, i.e. the phase difference, around
$t_{r}$. Actually, the only relevant contributions in the integration over
$\omega$ are the ones for which this difference is zero. Clearly, the
expansion of the phase difference around $t_{r}$ is proportional to the phase
velocity, which appears to be Doppler shifted, as shown in equation (10):

$\phi(\mathbf{r},t)-\phi(\mathbf{r}_{O}(t^{\prime}),t^{\prime})\approx\frac{\omega}{D_{\hat{\mathbf{k}}}[\mathbf{v}(t_{r})]}(t^{\prime}-t_{r})\quad\text{with}\quad\hat{\mathbf{k}}=\frac{\mathbf{r}-\mathbf{r}_{O}(t_{r})}{|\mathbf{r}-\mathbf{r}_{O}(t_{r})|},.$
(15)

We will omit the dependence of $\mathbf{v}$ from $t_{r}$ hereafter.
We substitute this new expression in (14),

$\displaystyle\psi(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{q}{2\pi\epsilon_{0}}\int dt^{\prime}\int
d\omega\frac{e^{i\omega(t^{\prime}-t_{r})/D_{\hat{\mathbf{k}}}(\mathbf{v})}}{4\pi|\mathbf{r}-\mathbf{r}_{O}(t^{\prime})|}$
(16) $\displaystyle=$ $\displaystyle\frac{q}{2\pi\epsilon_{0}}\int
dt^{\prime}\int
d\omega^{\prime}\frac{D_{\hat{\mathbf{k}}}(\mathbf{v})e^{i\omega^{\prime}(t^{\prime}-t_{r})}}{4\pi|\mathbf{r}-\mathbf{r}_{O}(t^{\prime})|}=$
$\displaystyle=$ $\displaystyle\frac{q}{4\pi\epsilon_{0}}\int
dt^{\prime}\frac{D_{\hat{\mathbf{k}}}(\mathbf{v})\delta(t^{\prime}-t_{r})}{|\mathbf{r}-\mathbf{r}_{O}(t^{\prime})|},$

where, the first step is a change of integration variable, from $\omega$ to
$\omega^{\prime}$, while the second step exploits the following identity

$\delta(x-x_{0})=\frac{1}{2\pi}\int dpe^{ip(x-x_{0})}.$ (17)

As a result of the integration over $t^{\prime}$ in equation (16), we find
that

$\psi(\mathbf{r},t)=\frac{q}{4\pi\epsilon_{0}}\frac{D_{\hat{\mathbf{k}}}(\mathbf{v})}{|\mathbf{r}-\mathbf{r}_{O}(t_{r})|}=\frac{1}{4\pi\epsilon_{0}}\frac{qc}{(c-\hat{\mathbf{k}}\cdot\mathbf{v})|\mathbf{r}-\mathbf{r}_{O}(t_{r})|}$
(18)

An analogous strategy can be applied to obtain $\mathbf{A}$, leading to the
expression

$\mathbf{A}(\mathbf{r},t)=\frac{\mathbf{v}}{c^{2}}\psi(\mathbf{r},t).$ (19)

We can add a few comments to what we have derived so far. Firstly, it is
possible to link our derivation with the one proposed in literature. The
evaluation of the charge density at the retarded time is usually expressed as
an integration over time of the Dirac delta

$\delta(t^{\prime}-t_{r}^{\prime})=\frac{\delta(t^{\prime}-t_{r})}{|g^{\prime}(t_{r})|}=\frac{1}{2\pi}\int
d\nu\frac{e^{i\nu(t^{\prime}-t_{r})}}{|g^{\prime}(t_{r})|}=\frac{1}{2\pi}\int
d\omega e^{i\omega(t^{\prime}-t_{r})|g^{\prime}(t_{r})|}$

where we used the fact that

$\delta[g(x)]=\sum_{n}\frac{\delta(x-x_{n})}{|g^{\prime}(x_{n})|},\quad\text{with}\quad
g(x_{n})=0,g^{\prime}(x_{n})\neq 0.$ (20)

In this case, $g(t^{\prime})=t^{\prime}-t+\mathcalligra{r}\,(t^{\prime})/c$,
whose zero is, by definition, the retarded time $t_{r}$. We use the notation
$\mathcalligra{r}\,(t^{\prime})=|\mathbf{r}-\mathbf{r}^{\prime}(t^{\prime})|$
and find that

$g^{\prime}(t^{\prime})=\frac{d}{dt^{\prime}}[t^{\prime}-t+\mathcalligra{r}\,(t^{\prime})/c]=1-\mathbf{v}(t^{\prime})\cdot\hat{\boldsymbol{\mathcalligra{r}}\,}/c=\frac{c-\mathbf{v}(t^{\prime})\cdot\hat{\boldsymbol{\mathcalligra{r}}\,}}{c}=\frac{1}{D_{\hat{\boldsymbol{\mathcalligra{r}}\,}}(\mathbf{v})}.$
(21)

which coincides with our results.
Other authors Zangwill ; Panofsky ; Orahilly  introduce the Doppler term in
different ways and suggest that it can be absorbed in the infinitesimal volume
element, stressing with this argument its geometrical origin. In our
procedure, instead, that factor is naturally embedded in the treatment and
explicitly derived from a Doppler mechanism.
More expert scholars could have spotted the presence of a Dirac delta in
expression (14) and, using equation (20), obtaining the final result without
passing through all the steps. Although it is undoubtedly correct and more
common, we think this haste becomes a lost opportunity to bring out the
physics.

##  III Conclusions

L-W potentials are a fundamental tool in classical and quantum
electrodynamics, although, sometimes the way they are derived hides the
physics underneath their origin.
Of course, the L-W as reported in the actual literature and textbooks of
electrodynamics are formally correct and in many cases exploring the physical
mechanism remains an unneeded subtext for their practical use.
Nevertheless, there could be cases where the physics underlying L-W potentials
could come in handy for the interpretation of processes closely related to L-W
potentials.
Here we have shown that the L-W potentials can be derived on the basis of a
Doppler effect, as suggested by the characteristic correction factor
$\left(\frac{v}{c}\right)$ presents in their analytical formulation.
Overall, the route followed here, besides being formally correct, makes
evident the physical phenomenology from which the electrodynamics L-W
potentials originate.

## Authors’ Declaration

The authors have no conflicts to disclose.

##  Appendix A Retarded fields

In this appendix we apply the same procedure presented above to find the
equations for the retarded fields, or Jefimenko-Feynman equations.Feynman
We start from Maxwell’s equations in the form of wave equations

$\displaystyle\left(\nabla^{2}-\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\right)\mathbf{E}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon_{0}}\nabla\rho+\mu_{0}\frac{\partial\mathbf{J}}{\partial
t}$ (22)
$\displaystyle\left(\nabla^{2}-\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\right)\mathbf{B}$ $\displaystyle=$
$\displaystyle-\nabla\times\mu_{0}\mathbf{J}$ (23)

whose solutions are

$\displaystyle\mathbf{E}$ $\displaystyle=$ $\displaystyle-\int dt^{\prime}\int
d\tau^{\prime}G_{-}(\mathbf{r},t|\mathbf{r}^{\prime},t^{\prime})\left[\frac{1}{\epsilon_{0}}\nabla\rho(\mathbf{r}^{\prime},t^{\prime})+\mu_{0}\frac{\partial\mathbf{J}(\mathbf{r}^{\prime},t^{\prime})}{\partial
t}\right]$ (24) $\displaystyle\mathbf{B}$ $\displaystyle=$ $\displaystyle\int
dt^{\prime}\int
d\tau^{\prime}G_{-}(\mathbf{r},t|\mathbf{r}^{\prime},t^{\prime})\left[\nabla\times\mu_{0}\mathbf{J}(\mathbf{r}^{\prime},t^{\prime})\right].$
(25)

Since the operators $\nabla$ and $\partial/\partial t$ act on the unprimed
coordinates, for the electric field we have that

$\mathbf{E}=-\frac{1}{4\pi\epsilon_{0}}\nabla\left(\int dt^{\prime}\int
d\tau^{\prime}G_{-}(\mathbf{r},t|\mathbf{r}^{\prime},t^{\prime})\rho(\mathbf{r}^{\prime},t^{\prime})\right)-\frac{\mu_{0}}{4\pi}\frac{\partial}{\partial
t}\left(\int dt^{\prime}\int
d\tau^{\prime}G_{-}(\mathbf{r},t|\mathbf{r}^{\prime},t^{\prime})\mathbf{J}(\mathbf{r}^{\prime},t^{\prime})\right),$

while for the magnetic field,

$\mathbf{B}=\mu_{0}\nabla\times\left(\int dt^{\prime}\int
d\tau^{\prime}G_{-}(\mathbf{r},t|\mathbf{r}^{\prime},t^{\prime})\mathbf{J}(\mathbf{r}^{\prime},t^{\prime})\right).$

The terms in parenthesis, called $I$ and $\mathbf{\Pi}$ hereafter, are exactly
the same integrals already encountered, when the definition (3) is inserted.
The remaining step is the evaluation of the derivatives with respect to space
and time

$\mathbf{E}=-\frac{1}{4\pi\epsilon_{0}}\nabla
I(\mathbf{r},t)-\frac{\mu_{0}}{4\pi}\frac{\partial}{\partial
t}\mathbf{\Pi}(\mathbf{r},t),\quad\mathbf{B}=\mu_{0}\nabla\times\mathbf{\Pi}(\mathbf{r},t)$
(26)

Following the steps proposed by Griffith Griffith , we need some useful
identities to carry out the derivation. Writing equation (8) as
$\mathcalligra{r}\,(t_{r})=c(t-t_{r})$, it is easy to show that

$\nabla\mathcalligra{r}\,=-c\nabla
t_{r}=\nabla(\sqrt{\boldsymbol{\mathcalligra{r}}\,\cdot\boldsymbol{\mathcalligra{r}}\,})=\frac{1}{2\sqrt{\boldsymbol{\mathcalligra{r}}\,\cdot\boldsymbol{\mathcalligra{r}}\,}}\nabla(\boldsymbol{\mathcalligra{r}}\,\cdot\boldsymbol{\mathcalligra{r}}\,)=\frac{1}{\mathcalligra{r}\,}\boldsymbol{\mathcalligra{r}}\,\cdot(\nabla\boldsymbol{\mathcalligra{r}}\,),$

multiplying both sides for $\mathcalligra{r}\,$ and considering the $i$-th
component

$-c\mathcalligra{r}\,\frac{\partial t_{r}}{\partial
x_{i}}=\mathcalligra{r}\,_{j}\frac{\partial\mathcalligra{r}\,^{j}}{\partial
x_{i}}=\mathcalligra{r}\,_{j}\left(\delta_{i}^{j}-v^{j}\frac{\partial
t_{r}}{\partial
x_{i}}\right)=\mathcalligra{r}\,_{i}-\mathcalligra{r}\,_{j}v^{j}\frac{\partial
t_{r}}{\partial x_{i}}$

or, in vectorial form,

$-c\mathcalligra{r}\,\nabla
t_{r}=\boldsymbol{\mathcalligra{r}}\,-(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v})\nabla
t_{r},$

hence

$\nabla\mathcalligra{r}\,=-c\nabla
t_{r}=\frac{c}{c-\hat{\boldsymbol{\mathcalligra{r}}\,}\cdot\mathbf{v}}\hat{\boldsymbol{\mathcalligra{r}}\,}=D_{\hat{\boldsymbol{\mathcalligra{r}}\,}}(\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}.$
(27)

Moreover,

$\displaystyle c\left(1-\frac{\partial t_{r}}{\partial
t}\right)=\frac{\partial}{\partial
t}c(t-t_{r})=\frac{\partial\mathcalligra{r}\,}{\partial
t}=\frac{\partial(\sqrt{\boldsymbol{\mathcalligra{r}}\,\cdot\boldsymbol{\mathcalligra{r}}\,})}{\partial
t}=\frac{1}{\mathcalligra{r}\,}\boldsymbol{\mathcalligra{r}}\,\cdot\frac{\partial\boldsymbol{\mathcalligra{r}}\,}{\partial
t}=-\frac{1}{\mathcalligra{r}\,}\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}\frac{\partial
t_{r}}{\partial t},$

and rearranging the first and last terms, we get

$\mathcalligra{r}\,c=(\mathcalligra{r}\,c-\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v})\frac{\partial
t_{r}}{\partial t},\quad\text{hence}\quad\frac{\partial t_{r}}{\partial
t}=D_{\hat{\boldsymbol{\mathcalligra{r}}\,}}(\mathbf{v}).$ (28)

Given a vectorial quantity $\mathbf{w}(t_{r})$, evaluated at the retarded
time, the directional derivative along a vector $\mathbf{d}$ and its curl are

$(\mathbf{d}\cdot\nabla)\mathbf{w}(t_{r})=\frac{d\mathbf{w}(t_{r})}{dt_{r}}\left(\mathbf{d}\cdot\nabla
t_{r}\right),\quad\nabla\times\mathbf{w}(t_{r})=-\frac{d\mathbf{w}(t_{r})}{dt_{r}}\times\nabla
t_{r}.$ (29)

and using (27),

$(\mathbf{d}\cdot\nabla)\mathbf{w}(t_{r})=-\frac{D_{\hat{\boldsymbol{\mathcalligra{r}}\,}}(\mathbf{v})}{c}\frac{d\mathbf{w}(t_{r})}{dt_{r}}\left(\mathbf{d}\cdot\hat{\boldsymbol{\mathcalligra{r}}\,}\right),\quad\nabla\times\mathbf{w}(t_{r})=\frac{D_{\hat{\boldsymbol{\mathcalligra{r}}\,}}(\mathbf{v})}{c}\frac{d\mathbf{w}(t_{r})}{dt_{r}}\times\hat{\boldsymbol{\mathcalligra{r}}\,}.$
(30)

We start from the gradient of $I$:

$\displaystyle\nabla I(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\nabla\left[\frac{D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v}(t_{r}))}{\mathcalligra{r}\,}\right]=D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v}(t_{r}))\nabla\left(\frac{1}{\mathcalligra{r}\,}\right)+\frac{\nabla
D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v}(t_{r}))}{\mathcalligra{r}\,}=$
$\displaystyle=$ $\displaystyle-
D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v}(t_{r}))\frac{\nabla\mathcalligra{r}\,}{\mathcalligra{r}\,^{2}}+\frac{\nabla
D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v}(t_{r}))}{\mathcalligra{r}\,}=$
$\displaystyle=$
$\displaystyle-\frac{D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v}(t_{r}))\hat{\boldsymbol{\mathcalligra{r}}\,}}{\mathcalligra{r}\,^{2}}+\frac{\nabla
D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v}(t_{r}))}{\mathcalligra{r}\,}$

We have to simplify the term $\nabla
D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})$, which leads to

$\nabla
D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})=\nabla\left(\frac{c}{c-\hat{\mathbf{\mathcalligra{r}\,}}\cdot\mathbf{v}}\right)=\frac{1}{c}\left(\frac{c}{c-\hat{\mathbf{\mathcalligra{r}\,}}\cdot\mathbf{v}}\right)^{2}\nabla(\hat{\mathbf{\mathcalligra{r}\,}}\cdot\mathbf{v})=\frac{1}{c}D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\nabla(\hat{\mathbf{\mathcalligra{r}\,}}\cdot\mathbf{v}).$

Now,

$\displaystyle\nabla(\hat{\mathbf{\mathcalligra{r}\,}}\cdot\mathbf{v})$
$\displaystyle=$
$\displaystyle\nabla\left(\frac{\boldsymbol{\mathcalligra{r}}\,}{\mathcalligra{r}\,}\cdot\mathbf{v}\right)=\nabla\left(\frac{1}{\mathcalligra{r}\,}\right)\left(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}\right)+\frac{1}{\mathcalligra{r}\,}(\nabla(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}))=$
$\displaystyle=$
$\displaystyle\frac{1}{\mathcalligra{r}\,}\left[-\frac{\nabla\mathcalligra{r}\,}{\mathcalligra{r}\,}\left(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}\right)+(\nabla(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}))\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\mathcalligra{r}\,}\left[-\frac{D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}}{\mathcalligra{r}\,}\left(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}\right)+(\nabla(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}))\right]$

We can exploit the vectorial identities (30) to compute
$\nabla(\mathbf{\mathcalligra{r}\,}\cdot\mathbf{v})$ and get

$\nabla(\mathbf{\mathcalligra{r}\,}\cdot\mathbf{v})=\mathbf{v}-\frac{1}{c}(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{a}-v^{2})D_{\hat{\boldsymbol{\mathcalligra{r}}\,}}(\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}.$
(31)

Collecting all together the steps, we arrive at

$\nabla
I(\mathbf{r},t)=-\frac{D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\mathcalligra{r}\,^{2}}\left\\{\hat{\boldsymbol{\mathcalligra{r}}\,}-\frac{1}{c}\left[\mathbf{v}-\frac{1}{c}(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{a}-v^{2})D_{\hat{\boldsymbol{\mathcalligra{r}}\,}}(\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}-\frac{D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}}{\mathcalligra{r}\,}\left(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}\right)\right]\right\\}.$

The first and last terms in the parenthesis simplify into

$1+\frac{\left(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{v}\right)}{c\mathcalligra{r}\,}D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})=\left[1+\frac{(\hat{\boldsymbol{\mathcalligra{r}}\,}\cdot\mathbf{v})}{c}\frac{c}{c-\hat{\boldsymbol{\mathcalligra{r}}\,}\cdot\mathbf{v}}\right]=D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})$

and finally,

$\displaystyle\nabla I(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle-\frac{D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\mathcalligra{r}\,^{2}}\left\\{D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}-\frac{1}{c}\left[\mathbf{v}-\frac{1}{c}(\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{a}-v^{2})D_{\hat{\boldsymbol{\mathcalligra{r}}\,}}(\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}\right]\right\\}=$
(32) $\displaystyle=$
$\displaystyle\frac{D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\mathcalligra{r}\,^{2}}\left[\frac{\mathbf{v}}{c}-\frac{D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{c^{2}}(c^{2}-v^{2}+\boldsymbol{\mathcalligra{r}}\,\cdot\mathbf{a})\hat{\boldsymbol{\mathcalligra{r}}\,}\right].$

Similarly, one can determine $\partial_{t}\mathbf{\Pi}$,

$\displaystyle\frac{\partial\mathbf{\Pi}(\mathbf{r},t)}{\partial t}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
t}[\mathbf{v}I(\mathbf{r},t)]=\left[I(\mathbf{r},t)\frac{\partial
t_{r}}{\partial t}\frac{\partial\mathbf{v}}{\partial
t_{r}}+\mathbf{v}\frac{\partial I(\mathbf{r},t)}{\partial
t}\right]=\left[I(\mathbf{r},t)D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\mathbf{a}+\mathbf{v}\frac{\partial
I(\mathbf{r},t)}{\partial t}\right]=$ $\displaystyle=$
$\displaystyle\left[\frac{D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\mathcalligra{r}\,}\mathbf{a}+\mathbf{v}\frac{\partial}{\partial
t}\left(\frac{D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\mathcalligra{r}\,}\right)\right]=\left[\frac{D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\mathcalligra{r}\,}\mathbf{a}+\frac{\mathbf{v}}{\mathcalligra{r}\,}\frac{\partial
D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\partial
t}+\mathbf{v}D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\frac{\partial}{\partial
t}\left(\frac{1}{\mathcalligra{r}\,}\right)\right].$

For the second term in RHS,

$\displaystyle\frac{\partial
D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\partial t}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
t}\left(\frac{c}{c-\mathbf{v}\cdot\hat{\boldsymbol{\mathcalligra{r}}\,}}\right)=\frac{1}{c}D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\frac{\partial}{\partial
t}(\mathbf{v}\cdot\hat{\boldsymbol{\mathcalligra{r}}\,})=\frac{1}{c}D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\frac{\partial}{\partial
t}\left(\mathbf{v}\cdot\frac{\boldsymbol{\mathcalligra{r}}\,}{\mathcalligra{r}\,}\right)=$
$\displaystyle=$
$\displaystyle\frac{1}{c}D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})(\mathbf{v}\cdot\boldsymbol{\mathcalligra{r}}\,)\frac{\partial}{\partial
t}\left(\frac{1}{\mathcalligra{r}\,}\right)+\frac{1}{c\mathcalligra{r}\,}D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\frac{\partial}{\partial
t}(\mathbf{v}\cdot\boldsymbol{\mathcalligra{r}}\,),$

in which

$\displaystyle\frac{\partial}{\partial
t}\left(\frac{1}{\mathcalligra{r}\,}\right)=-\frac{c}{\mathcalligra{r}\,^{2}}[1-D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})],\quad\text{and}\quad\frac{\partial}{\partial
t}(\mathbf{v}\cdot\boldsymbol{\mathcalligra{r}}\,)=\left(\mathbf{a}\cdot\boldsymbol{\mathcalligra{r}}\,-v^{2}\right)D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v}).$
(33)

The time derivative is

$\frac{\partial\mathbf{\Pi}(\mathbf{r},t)}{\partial
t}=\frac{D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{\mathcalligra{r}\,^{2}}\left[\mathcalligra{r}\,\mathbf{a}-c\mathbf{v}+\frac{\mathbf{v}D_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{c}\left(c^{2}-v^{2}+\mathbf{a}\cdot\boldsymbol{\mathcalligra{r}}\,\right)\right]$
(34)

As explained in Griffith , the calculation of $\nabla\times\mathbf{\Pi}$ can
be simplified using the expressions already presented.
The final results are

$\displaystyle\mathbf{E}$ $\displaystyle=$
$\displaystyle-\frac{q}{4\pi\epsilon_{0}}\left(\nabla
I(\mathbf{r},t)+\frac{1}{c^{2}}\frac{\partial\mathbf{\Pi}(\mathbf{r},t)}{\partial
t}\right)=-\frac{qD^{3}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})}{4\pi\epsilon_{0}\mathcalligra{r}\,^{2}}\left[(c^{2}-v^{2})\mathbf{u}+\boldsymbol{\mathcalligra{r}}\,\times(\mathbf{u}\times\mathbf{a})\right]$
(35) $\displaystyle\mathbf{B}$ $\displaystyle=$
$\displaystyle\frac{1}{c}\hat{\boldsymbol{\mathcalligra{r}}\,}\times\mathbf{E}(\mathbf{r},t),$
(36)

where $\mathbf{u}=\hat{\boldsymbol{\mathcalligra{r}}\,}-\mathbf{v}/c$. Notice
that, if we define the new directions

$\displaystyle\mathbf{u}_{\parallel}$ $\displaystyle=$
$\displaystyle(\mathbf{u}\cdot\hat{\boldsymbol{\mathcalligra{r}}\,})\hat{\boldsymbol{\mathcalligra{r}}\,}=D^{-1}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}$
(37) $\displaystyle\mathbf{u}_{\perp}$ $\displaystyle=$
$\displaystyle\mathbf{u}-\mathbf{u}_{\parallel}=\frac{1}{c}\left[(\hat{\boldsymbol{\mathcalligra{r}}\,}\cdot\mathbf{v})\hat{\boldsymbol{\mathcalligra{r}}\,}-\mathbf{v}\right],$
(38)

the field along $\mathbf{u}_{\parallel}$ scales as
$D^{2}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})/\mathcalligra{r}\,^{2}$,
while in the other direction, namely $\mathbf{u}_{\perp}$, is proportional to
$D^{3}_{\hat{\mathbf{\mathcalligra{r}\,}}}(\mathbf{v})/\mathcalligra{r}\,^{2}$.
The latter corresponds to the inverse direction of the velocity component
perpendicular to the source position.
As suggested in Heras1 , it is possible to retrieve the retarded potentials
from equations (26), invoking the retarded Helmholtz theorem Heras2 .

## References

## References

  * (1) D. J. Griffith, Introduction to Electrodynamics (Prentice Hall, Upper Saddle River, NJ, 1999), chap. 10.
  * (2) ”Please note that this (e.g. the L-W potentials) has nothing whatever to do with special relativity or Lorentz contraction, L is the length of the moving train, and its rest length is not at issue. The argument is somewhat reminiscent of the Doppler effect.” David J. Griffith, Introduction to electrodynamics, (Cambridge University, UK, 2013) pag. 452, note 12.
  * (3) J. D. Jackson, Classical Electrodynamics (Wiley, New York, NY, 1999), 3rd ed.
  * (4) J. A. Stratton, Electromagnetic Theory (Wiley, New York, NY, 2015).
  * (5) A. Zangwill, Modern Electrodynamics (Cambridge, New York, 2013).
  * (6) R. P. Feynman, R. B. Leighton, and M. Sands The Feynman Lecture on Physics (Addison-Wesley, Reading, MA, 1964), chap. 21.
  * (7) W. K. H. Panofsky and M. Phillips, Classical Electricity and Magnetism (Addison-Wesley, Reading, MA, 1975), 2nd ed.
  * (8) A. O’Rahilly, Electromagnetic Theory, (Dover, New York, NY, 1965), 213 and 214.
  * (9) J. R. Reitz, F. J. Milford, R. W. Christy, Foundations of Electromagnetic Theory, (Addison-Wesley, Reading, MA 1979) 3rd ed, sec. 21.1.
  * (10) F. Tommasini and A. Morgante, Correnti, Radiazioni e Quanti (Edizioni Goliardiche, Trieste, Italy, 2003).
  * (11) A. Liénard, ”Champ électrique e magnétique” L’Èclairage Èlect., Tome XVI, Annèe 5, 27, 5-14 (1898).
  * (12) E. Wiechert, ”Elektrodynamische Elementargesetze” Arch. Neerl. Sci. Exact Natur. 5, 549 (1900); ”Elektrodynamische Elementargesetze”, Ann. Phys. 309, 667 (1901).
  * (13) H. Minkoski, ”Raum un Zeit”, Jahresb. Deutschen Math.-Verein 18, 75-88 (1909).
  * (14) H. Padmanabhan, ”A simple derivation of the electromagnetic field of an arbitrarily moving charge”, Am. J. Phys. 77, 151 (2009).
  * (15) V. N. Strel’tsov, ”Liénard-Wiechert Potential as a Consequence of Lorentz Transformation of Coulomb Potential”, Technical Report JINR-D-2-93-437 (1993).
  * (16) H. A. Haus, ”On the radiation of point charge”, Am. J. Phys. 54, 1126 (1986).
  * (17) P. Frank and R. V. Mises, Die Differential-and Integralgleichungen der Mechanik und Physik, Zweiter phylikalischer Teil (Wieveg und Sohn, Braunschweig, 1935).
  * (18) J. Schwinger, L. L. DeRaad, Jr., K. A. Milton, W.Tsai, Classical Electrodynamics (CRC Press Taylor & Francis Group, 2018), pag. 346-347.
  * (19) S. Hassani, Mathematical Physics (Springer, 2013) 2nd ed.
  * (20) R. Heras, ”Alternative routes to the retarded potentials”, Eur. J. Phys. 38, 055203 (2017).
  * (21) R. Heras, ”The Helmholtz theorem and retarded fields”, Eur. J. Phys. 37, 065204 (2016).


# Robotic Manifold Tracking of Coherent Structures in Flows

M. Ani Hsieh, Eric Forgoston, T. William Mather, and Ira B. Schwartz This work
was supported by the Office of Naval Research (ONR). MAH was an ONR Summer
Faculty Fellow and supported by ONR Contract No. N0001411WX20079. EF is
supported by the Naval Research Laboratory (Award No. N0017310-2-C007) TWM is
supported by the National Science Foundation under Grant No. DGE-0947936.M.
Ani Hsieh and T. William Mather are with the SAS Laboratory, Mechanical
Engineering & Mechanics Department, Drexel University, Philadelphia, PA 19104,
USA<EMAIL_ADDRESS>Forgoston is with the Department of
Mathematical Sciences, Montclair State University, Montclair, New Jersey
07043, USA<EMAIL_ADDRESS>B. Schwartz is with the Nonlinear
Systems Dynamics Section, Plasma Physics Division, Code 6792, U.S. Naval
Research Laboratory, Washington, DC 20375, USA<EMAIL_ADDRESS>
###### Abstract

Tracking Lagrangian coherent structures in dynamical systems is important for
many applications such as oceanography and weather prediction. In this paper,
we present a collaborative robotic control strategy designed to track stable
and unstable manifolds. The technique does not require global information
about the fluid dynamics, and is based on local sensing, prediction, and
correction. The collaborative control strategy is implemented on a team of
three robots to track coherent structures and manifolds on static flows as
well as a noisy time-dependent model of a wind-driven double-gyre often seen
in the ocean. We present simulation and experimental results and discuss
theoretical guarantees of the collaborative tracking strategy.

##  I INTRODUCTION

In this paper, we present a collaborative control strategy for a class of
autonomous underwater vehicles (AUVs) to track the coherent structures and
manifolds on flows. In realistic ocean flows, these time-dependent coherent
structures, or Lagrangian coherent structures (LCS), are similar to
separatrices that divide the flow into dynamically distinct regions. LCS are
extensions of stable and unstable manifolds to general time-dependent flows
[1], and they carry a great deal of global information about the dynamics of
the flows. For two-dimensional (2D) flows, LCS are analogous to ridges defined
by local maximum instability, and quantified by local measures of Finite-Time
Lyapunov Exponents (FTLE) [2].

Recently, LCS have been shown to coincide with optimal trajectories in the
ocean which minimize the energy and the time needed to traverse from one point
to another [3, 4]. Furthermore, to improve weather and climate forecasting,
and to better understand various physical, chemical, and geophysical processes
in the ocean, there has been significant interest in the deployment of
autonomous sensors to measure a variety of quantities of interest. One
drawback to operating sensors in time-dependent and stochastic environments
like the ocean is that the sensors will tend to escape from their monitoring
region of interest. Since the LCS are inherently unstable and denote regions
of the flow where more escape events may occur [5], knowledge of the LCS are
of paramount importance in maintaining a sensor in a particular monitoring
region.

Existing work in cooperative boundary tracking for robotic teams that relies
on one-dimensional (1D) parameterizations include [6, 7] and [8, 9] for static
and time-dependent cases respectively. Formation control strategies for
distributed estimation of level surfaces and scalar fields in the ocean are
presented in [10, 11, 12] and pattern formation for surveillance and
monitoring by robot teams is discussed in [13, 14, 15]. Our work is
distinguished from existing work in that we use cooperative robots to find
coherent structures without requiring a global picture of the ocean dynamics.
We take inspiration from [16] and design a strategy to enable a team of robots
to track the stable/unstable manifolds of general 2D conservative flows
through local sensing alone. We verify the feasibility of our method through
simulations and experiments and show how the proposed strategy can be extended
to track coherent structures in time-dependent conservative flows with
measurement noise. To our knowledge, this is the first attempt in the
development of tracking strategies for mapping LCS in the ocean using AUVs.

The novelty of this work lies in the use of nonlinear dynamical and chaotic
system analysis techniques to derive a tracking strategy for a team of robots.
The cooperative control strategy leverages the spatio-temporal sensing
capabilities of a team of networked robots to track the boundaries separating
the regions in phase space that support distinct dynamical behavior.
Additionally, our boundary tracking relies solely on local measurements of the
velocity field. Our technique is quite general, and may be applied to any
conservative flow.

The paper is structured as follows: We formulate the problem and outline key
assumptions in Section II. The cooperative control strategy is presented in
Section III and its theoretical properties analyzed in Section IV. Section V
presents our simulation and experimental results. The extension of the
proposed strategy to a noisy time-dependent model of a wind-driven double-gyre
is presented in VI. We conclude with a discussion of our results and
directions for future work in Sections VII and VIII respectively.

##  II PROBLEM FORMULATION

We consider the problem of controlling a team of $N$ planar AUVs to
collaboratively track the material lines that separate regions of flow with
distinct fluid dynamics. This is similar to the problem of tracking the stable
(and unstable) manifolds of a general nonlinear dynamical system where the
manifolds separate regions in phase space with distinct dynamical behaviors.
We assume the following 2D kinematic model for each of the AUVs:

$\dot{x}_{i}=V_{i}\cos\theta_{i}+u_{i},$ (1a)
$\dot{y}_{i}=V_{i}\sin\theta_{i}+v_{i},$ (1b)

where $\mathbf{x}_{i}=[x_{i},\,y_{i}]^{T}$ is the vehicle’s planar position,
$V_{i}$ and $\theta_{i}$ are the vehicle’s linear speed and heading, and
$\mathbf{u}_{i}=[u_{i},\,v_{i}]^{T}$ is the velocity of the fluid current
experienced/measured by the $i^{th}$ vehicle. Additionally, we assume each
agent can be circumscribed by a circle of radius $r$, i.e., each vehicle can
be equivalently described as a disk of radius $r$.

In this work, $\mathbf{u}_{i}$ is provided by a 2D planar conservative vector
field described by a differential equation of the form

$\mathbf{\dot{x}}=F(\mathbf{x}).$ (2)

In essence, $u_{i}=F_{x}(\mathbf{x}_{i})$ and $v_{i}=F_{y}(\mathbf{x}_{i})$.
Let $B_{S}$ and $B_{U}$ denote the stable and unstable manifolds of (2). In
general, $B_{S}$ and $B_{U}$ are the separating boundaries between regions in
phase space with distinct dynamics. For 2D flows, $B_{*}$ are simply one-
dimensional curves where $*$ denotes either stable ($S$) or unstable ($U$)
boundaries. For a small region centered about a point on $B_{*}$, the system
is unstable in one dimension. Finally, let $\rho(B_{*})$ denote the radius of
curvature of $B_{*}$ and assume that the minimum of the radius of curvature
$\rho_{min}(B_{*})>r$. This last assumption is needed to ensure the robots do
not lose track of the $B_{*}$ due to sharp turns.

The objective is to develop a collaborative strategy to enable a team of
robots to track $B_{*}$ in general 2D planar conservative flow fields through
local sampling of the velocity field. While the focus is on the development of
a tracking strategy for $B_{S}$, the proposed method can be easily extended to
track $B_{U}$ since $B_{U}$ are simply stable manifolds of (2) for $t<0$. We
present our methodology in the following section.

##  III METHODOLOGY

Our methodology is inspired by the Proper Interior Maximum (PIM) Triple
Procedure [16] – a numerical technique designed to find stationary
trajectories in chaotic regions with no attractors. While the original
procedure was developed for chaotic dynamical systems, the approach can be
employed to reveal the stable set of a saddle point of a general nonlinear
dynamical system. The procedure consists of iteratively finding an appropriate
PIM Triple on a saddle straddling line segment and propagating the triple
forward in time. We briefly summarize the procedure in the following section
and refer the interested reader to [16] for further details.

###  III-A The PIM Triple Procedure

Given the dynamical system described by (2), let ${\cal D}\in\mathbb{R}^{2}$
be a closed and bounded set such that ${\cal D}$ does not contain any
attractors of (2). Given a point $\mathbf{x}\in{\cal D}$, the escape time of
$\mathbf{x}$, denoted by $T_{E}(\mathbf{x})$, is the time $\mathbf{x}$ takes
to leave the region ${\cal D}$ under the differential map (2).

Let $J$ be a line segment that crosses the stable set $B_{S}$ in ${\cal D}$,
i.e., the endpoints of the $J$ are on opposite sides of $B_{S}$. Let
$\\{\mathbf{x}_{L},\mathbf{x}_{C},\mathbf{x}_{R}\\}$ denote a set of three
points in $J$ such that $\mathbf{x}_{C}$ denotes the interior point. Then
$\\{\mathbf{x}_{L},\mathbf{x}_{C},\mathbf{x}_{R}\\}$ is an Interior Maximum
triple if
$T_{E}(\mathbf{x}_{C})>\max\\{T_{E}(\mathbf{x}_{L}),T_{E}(\mathbf{x}_{R})\\}$.
Furthermore, $\\{\mathbf{x}_{L},\mathbf{x}_{C},\mathbf{x}_{R}\\}$ is a Proper
Interior Maximum (PIM) triple if it is an Interior Maximum triple and the
interval $[\mathbf{x}_{L},\mathbf{x}_{R}]$ in $J$ is a proper subset of $J$.

Then the numerical computation of any PIM triple can be obtained iteratively
starting with an initial saddle straddle line segment $J_{0}$. Let
$\mathbf{x}_{L_{0}}$ and $\mathbf{x}_{R_{0}}$ denote the endpoints of $J_{0}$
and apply an $\epsilon_{0}>0$ discretization of $J_{0}$ such that
$\mathbf{x}_{L_{0}}=\mathbf{q}_{0}<\mathbf{q}_{1}<\ldots<\mathbf{q}_{M}=\mathbf{x}_{R_{0}}$.
For every point $\mathbf{q}_{i}$, determine $T_{E}(\mathbf{q}_{i})$ by
propagating $\mathbf{q}_{i}$ forward in time using (2). Then the PIM triple in
$J_{0}$ is given by the the points
$\\{\mathbf{q}_{k-1},\mathbf{q}_{k},\mathbf{q}_{k+1}\\}$ where
$\mathbf{q}_{k}=\arg\max\limits_{i=1,\ldots,M}T_{E}(\mathbf{q}_{i})$. This PIM
triple can then be further refined by choosing $J_{1}$ to be the line segment
containing $\\{\mathbf{q}_{k-1},\mathbf{q}_{k},\mathbf{q}_{k+1}\\}$ and
reapplying the procedure with another $\epsilon_{1}>0$ discretization where
$\epsilon_{1}<\epsilon_{0}$.

Given an initial saddle straddling line segment $J_{0}$, it has been shown
that the line segment given by any subsequent PIM triple on $J_{0}$ is also a
saddle straddling line segment [16]. Furthermore, if we use a PIM triple
$\mathbf{x}(t)=\\{\mathbf{x}_{L},\mathbf{x}_{C},\mathbf{x}_{R}\\}$ as the
initial conditions for the dynamical system given by (2) and propagate the
system forward in time by $\Delta t$, then the line segment containing the set
$\mathbf{x}(t+\Delta t)$, $J_{t+\Delta t}$, remains a saddle straddle line
segment. As such, the same numerical procedure can be employed to determine an
appropriate PIM trip on $J_{t+\Delta t}$. This procedure can be repeated to
eventually reveal the entire stable set $B_{S}$ and unstable set $B_{U}$
within ${\cal D}$ if time was propagated forwards and backwards respectively.
Furthermore, since the procedure always begins with a valid saddle straddling
line segment, by construction, the procedure always results in a non-empty
set.

Inspired by the PIM Triple Procedure, we propose a cooperative saddle straddle
control strategy for a team of $N\geq 3$ robots to track the stable (and
unstable) manifolds of a general conservative time-independent flow field
$F(\mathbf{x})$. Different from the procedure, our robots will solely rely on
information that can be gathered via local sensing and shared through the
network. In contrast, a straight implementation of the PIM Triple Procedure
would require global knowledge of the structure of the system dynamics
throughout a given region given its reliance on computing escape times. We
describe our approach in the following section.

###  III-B Controller Synthesis

Consider a team of three robots and identify them as robots $\\{L,C,R\\}$.
While the robots may be equipped with similar sensing and actuation
capabilities, we propose a heterogeneous cooperative control strategy.

Let
$\mathbf{x}(0)=[\mathbf{x}_{L}^{T}(0),\,\mathbf{x}_{C}^{T}(0),\,\mathbf{x}_{R}^{T}(0)]^{T}$
be the initial conditions for the three robots. Assume that $\mathbf{x}(0)$
lies on the line segment $J_{0}$ where $J_{0}$ is a saddle straddle line
segment and $\\{\mathbf{x}_{L}(0),\mathbf{x}_{C}(0),\mathbf{x}_{R}(0)\\}$
constitutes a PIM triple. Similar to the PIM Triple Procedure, the objective
is to enable the robots to maintain a formation such that a valid saddle
straddle line segment can be maintained between robots $L$ and $R$. Instead of
computing the escape times for points on $J_{0}$ as proposed by the PIM Triple
Procedure, robot $C$ must remain close to $B_{S}$ using only local
measurements of the velocity field provided by the rest of the team. As such,
we refer to robot $C$ as the tracker of the team while robots $L$ and $R$
maintains a straddle formation across the boundary at all times. Robots $L$
and $R$ may be thought of herding robots, since they control and determine the
actions of the tracking robot.

####  III-B1 Straddling Formation Control

The controller for the straddling robots consists of two discrete states: a
passive control state, $U_{P}$, and an active control state, $U_{A}$. The
robots initialize in the passive state $U_{P}$ where the objective is to
follow the flow of the ambient vector field. Therefore, $V_{i}=0$ for $i=L,R$.
Robots execute $U_{P}$ until they reach the maximum allowable separation
distance $d_{Max}$ from robot $C$. When
$\|\mathbf{x}_{i}-\mathbf{x}_{C}\|>d_{Max}$, robot $i$ switches to the active
control state, $U_{A}$, where the objective is to navigate to a point
$\mathbf{p}_{i}$ on the current projected saddle straddle line segment
$\hat{J}_{t}$ such that $\|\mathbf{p}_{i}-\mathbf{p}_{c}\|=d_{Min}$ and
$\mathbf{p}_{C}$ denotes the midpoint of $\hat{J}_{t}$. When robots execute
$U_{A}$ , $V_{i}=\|(\mathbf{p}_{i}-\mathbf{x}_{i})-\mathbf{u}_{i}\|$ and
$\theta_{i}(t)=\alpha_{i}(t)$ where $\alpha_{i}$ is the angle between the
desired, $(\mathbf{p}_{i}-\mathbf{x}_{i})$, and current heading,
$\mathbf{u}_{i}$, of robot $i$ as shown in Fig. 1. In summary, the straddling
control strategy for robots $L$ and $R$ is given by

$\displaystyle V_{i}$ $\displaystyle=\left\\{\begin{array}[]{ll}0&\textrm{if
}d_{Min}<\|\mathbf{x}_{i}-\mathbf{x}_{C}\|<d_{Max}\\\
\|(\mathbf{p}_{i}-\mathbf{x}_{i})-\mathbf{u}_{i}\|&\textrm{otherwise}\end{array}\right.,$
(3c) $\displaystyle\theta_{i}$
$\displaystyle=\left\\{\begin{array}[]{ll}0&\textrm{if
}d_{Min}<\|\mathbf{x}_{i}-\mathbf{x}_{C}\|<d_{Max}\\\
\alpha_{i}&\textrm{otherwise}\end{array}\right.~{}.$ (3f)

Figure 1: Three robots tracking $B_{S}$ in a conservative vector field. The
blue dash-dot lines represent the robot trajectories, the green dashed line
represents the saddle straddle line segment $J$, and $\mathbf{p}_{L}$ and
$\mathbf{p}_{R}$ denotes the target positions for $L$ and $R$ respectively
when executing $U_{P}$ and $U_{A}$.

We note that while the primary control objective for robots $L$ and $R$ is to
maintain a straddle formation across $B_{S}$, robots $L$ and $R$ are also
constantly sampling the velocity of the local vector field and communicating
these measurements and their relative positions to robot $C$. Robot $C$ is
then tasked to use these measurements to track the position of $B_{S}$.

####  III-B2 Manifold Tracking Control

Let $\mathbf{\hat{u}}_{L}(t)$, $\mathbf{\hat{u}}_{C}(t)$, and
$\mathbf{\hat{u}}_{C}(t)$ denote the current velocity measurements obtained by
robots $L$, $C$, and $R$ at their respective positions. Let $d(\cdot,\cdot)$
denote the Euclidean distance function and assume that
$d(\mathbf{x}_{C},B_{S})<\epsilon$ such that $\epsilon>0$ is small. Given the
straddle line segment $J_{t}$ such that $\mathbf{x}_{L}(k)$ and
$\mathbf{x}_{R}(k)$ are the endpoints of $J_{t}$, we consider an
$\epsilon_{t}<\epsilon$ discretization of $J_{t}$ such that
$\mathbf{x}_{L}=\mathbf{q}_{1}<\mathbf{q}_{2}<\ldots<\mathbf{q}_{M}=\mathbf{x}_{R}$.
The objective is to use the velocity measurements provided by the team to
interpolate the vector field at the points
$\mathbf{q}_{1},\ldots,\mathbf{q}_{M}$. Since (2) has ${\cal C}^{1}$
continuity and if $\mathbf{x}_{C}$ is $\epsilon$-close to $B_{S}$, then the
point
$\mathbf{q}_{B}=\arg\max\limits_{k=1,\ldots,M}\mathbf{u}(\mathbf{q}_{k})^{T}\mathbf{\hat{u}}_{C}(t)$
should be $\delta$-close to $B_{S}$ where $\epsilon<\delta<A$ and $A$ is a
small enough positive constant.

While there are numerous vector field interpolation techniques available [17,
18, 19], we employ the inverse distance weighting method described in [17].
For a given set of velocity measurements $\mathbf{\hat{u}}_{i}(t)$ and
corresponding position estimates $\mathbf{\hat{x}}_{i}(t)$, the velocity
vector at some point $\mathbf{q}_{k}$ is given by

$\displaystyle\mathbf{u}(\mathbf{q}_{k})=\sum_{j}\sum_{i=1}^{N}\frac{w_{ij}\mathbf{\hat{u}}_{i}(j)}{\sum_{j}\sum_{i=1}^{N}w_{ij}}$

where $w_{ij}=\|\mathbf{\hat{x}}_{i}(j)-\mathbf{q}_{i}\|^{-2}$. Rather than
rely solely on the current measurements provided by the three robots, it is
possible to include the recent history of $\mathbf{\hat{u}}_{i}(t)$ to improve
the estimate of $\mathbf{u}(\mathbf{q}_{k})$, i.e.,
$\mathbf{\hat{u}}_{i}(t-\Delta T)$, $\mathbf{\hat{u}}_{i}(t-2\Delta T)$, and
so on, where $\Delta T$ is the sampling period and $i=\\{L,C,R\\}$. Thus, the
control strategy for the tracking robot $C$ is given by

$\displaystyle V_{C}$
$\displaystyle=\|[(\mathbf{q}_{B}+b\mathbf{\hat{u}}_{B})-\mathbf{x}_{C}]-\mathbf{u}_{C}\|$
(4a) $\displaystyle\theta_{C}$ $\displaystyle=\beta_{C}$ (4b)

where $\beta_{C}$ denotes the difference in the heading of robot $C$ and the
vector $(\mathbf{q}_{B}-\mathbf{\hat{u}}_{B})$ and $b>r$ is a small number.
The term $b\mathbf{\hat{u}}_{B}$ is included to ensure that the control
strategy aims for a point in front of robot $C$ rather than behind it. As
such, the projected saddle straddle line segment $\hat{J}_{t}$ at each time
step is given by $\mathbf{p}_{c}=q_{C}+b\mathbf{u}_{C}$ with $\hat{J}_{t}$
orthogonal to $B_{S}$ at $q_{C}$ and $\|\hat{J}_{t}\|$ chosen to be in the
interval $[2d_{Min},2d_{Max}]$.

##  IV ANALYSIS

In this section, we discuss the theoretical feasibility of the proposed saddle
straddle control strategy. We begin with the following key assumption on the
robots’ initial positions.

###### Assumption 1

Given a team of three robots $\\{L,C,R\\}$, assume that
$d(\mathbf{x}_{C}(0),B_{S})<\epsilon$ for a small value of $\epsilon>0$,
$\|\mathbf{x}_{L}-\mathbf{x}_{C}\|=\|\mathbf{x}_{R}-\mathbf{x}\|=d_{Min}$ with
$d_{Min}>2r$, and robots $L$ and $R$ are on opposite sides of $B_{S}$.

In other words, assume that the robots initialize in a valid PIM triple
formation and their positions form a saddle straddle line segment orthogonal
to $B_{S}$. Our main result concerns the validity of the saddle straddle
control strategy.

###### Theorem 1

Given a team of $3$ robots with kinematics given by (1) and $\mathbf{u}_{i}$
given by (2), the feedback control strategy (3) and (4) maintains a valid
saddle straddle line segment in the time interval $[t,t+\Delta t]$ if the
initial positions of the robots, $\mathbf{x}(t)$, is a valid PIM triple.

###### Proof:

To show this, we must show that at time $t+\Delta t$, robots $L$ and $R$
remain on opposite sides of $B_{S}$. Consider the rate of change of the
following function

$\displaystyle H(\mathbf{x}_{L},\mathbf{x}_{R})$
$\displaystyle=\frac{1}{2}(\mathbf{x}_{L}-\mathbf{x}_{R})^{T}(\mathbf{x}_{L}-\mathbf{x}_{R}).$

The above expression is simply one half the square of the distance between
robots $L$ and $R$. Let $J_{t}$ denote the saddle straddle line segment
defined by $\mathbf{x}_{L}(t)$ and $\mathbf{x}_{R}(t)$ at $t$ and let
$\mathbf{p}_{B}$ be the intersection of $J_{t}$ and $B_{S}$. By construction,
if we linearize (2) about the point $\mathbf{p}_{B}$, then the Jacobian of (2)
at $\mathbf{p}_{B}$ will have one positive eigenvalue. Furthermore, the
linearized system can be diagonalized such that the direction of instability
lies along $J_{t}$ [20]. Thus, $\frac{d}{dt}H>0$ in the time interval
$[t,t+\Delta t]$ when $V_{i}=0$ in (3).

When $V_{i}\neq 0$ in (3) for $i=L,R$, $\frac{d}{dt}H<0$ if the robots $L$ and
$R$ are moving closer to robot $C$ after reaching the maximum allowable
separation distance. Recall $\rho_{min}(B_{S})>r$, the smallest radius of
curvature of $B_{S}$, and $d_{Min}>2r$. Furthermore, robot $C$ initializes
$\epsilon$-close to the boundary and (4) steers $C$ towards $\mathbf{p}_{C}$
on $\hat{J}_{t}$ where $\hat{J}_{t}$ is orthogonal to $B_{S}$ at
$\mathbf{x}_{C}$. This ensures that the rate of the change of the radius of
curvature of the manifold $B_{S}$ is small enough such that $\hat{J}_{t}$
intersects with $B_{S}$ only once. Since $d_{Min}>2r$, this ensures that even
if $\frac{d}{dt}H<0$, the straddling robots never cross the boundary as they
move closer to the tracking robot. ∎

The above theorem guarantees that for any given time interval $[t,t+\Delta t]$
the team maintains a valid PIM triple formation. As such, the iterative
application of the proposed control strategy leads to the following
proposition.

###### Proposition 1

Given a team of $3$ robots with kinematics given by (1) and $\mathbf{u}_{i}$
given by (2), the feedback control strategy (3) and (4) results in an estimate
of $B_{S}$, denoted as $\hat{B}_{S}$, such that
$\left<B_{S},\hat{B}_{S}\right>_{L_{2}}<W$ for some $W>0$ where
$\left<\cdot,\cdot\right>_{L_{2}}$ denotes the inner product (which provides
an $L_{2}$ measure between the $B_{S}$ and $\hat{B}_{S}$ curves).

From Thm. 1, since the team is able to maintain a valid PIM triple formation
across $B_{S}$ for any given time interval $[t,t+\Delta t]$, this ensures that
an estimate of $B_{S}$ in the given time interval also exists. Applying this
reasoning in a recursive fashion, one can show that an estimate of $B_{S}$ can
be obtained for any arbitrary time interval. The challenge, however, lies in
determining the bound on $W$ such that $\hat{B}_{S}$ results in a good enough
approximation since $W$ depends on the sensor and actuation noise, the vector
interpolation routine, the sampling frequency, and the time scales of the flow
dynamics. This is a direction for future work.

##  V RESULTS

###  V-A Simulations

We illustrate the proposed control strategy given by (3) and (4) with some
simulation results. Fig. 2 shows the trajectories of three robots tracking a
sinusoidal boundary while Fig. 2 shows the team tracking a 1D star-shaped
boundary. We note that throughout the entire length of the simulation, the
team maintains a saddle straddle formation across the boundary.

In both examples, $\mathbf{u}=-a\nabla\varphi-b\nabla\times\psi$ where $a,b>0$
and $\varphi$ is an artificial potential function such that
$\varphi(\mathbf{x})=0$ for all $\mathbf{x}\in B_{*}$ and
$\varphi(\mathbf{x})<0$ for any $\mathbf{x}\in\mathbb{R}^{2}/B_{*}$. The
vector $\psi$ is a $3\times 1$ vector whose entries are given by
$[0,\,0,\,\gamma(x,y)]^{T}$ where $\gamma(x,y)$ is the curve describing the
desired boundary [15]. Lastly, the estimated position of the boundary is given
by the position of the tracking robot, i.e., robot $C$. In these examples, we
filtered the boundary position using a simple first-order low pass filter.

Figure 2: Trajectories of $3$ robots tracking a (a) sinusoidal boundary and a
(b) star-shaped boundary. The red dashed line is the estimated position of the
desired boundary shown in solid black. The start positions are shown by
$\triangle$ and the end positions are shown by the circle-enclosed blue
triangles.

###  V-B Experiments

We also implemented the control strategy on our multi-robot testbed. The
testbed consisted of three mSRV-1 robots in a 4.8x5.4 meter workspace. The
mSRV-1 are differential-drive robots equipped with an embedded processor,
color camera, and 802.11 wireless capability. Localization for each robot was
provided via a network of overhead cameras. Fig. 3 shows the trajectories of
the robots tracking a star shaped boundary shown in black. Fig. 3 is a
snapshot of the experimental run. We refer the interested reader to the
attached multimedia file for a movie of the full simulation and experimental
runs.

Figure 3: Trajectories of the $3$ robot team tracking a (a) star shape. The
red dashed line is the estimated position of the desired boundary shown in
solid black. The start positions are shown by $\triangle$ and the end
positions are shown by $\bigcirc$. (b) Snapshot of the multi-robot experiment.

##  VI EXTENSION TO PERIODIC BOUNDARIES

In this section, we consider the system of $3$ robots with kinematics given by
(1) where $\mathbf{u}_{i}$ is determined by the wind-driven double-gyre flow
model with noise

$\displaystyle\dot{x}$ $\displaystyle=-\pi
A\sin(\pi\frac{f(x,t)}{s})\cos(\pi\frac{y}{s})-\mu x+\eta_{1}(t),$ (5a)
$\displaystyle\dot{y}$ $\displaystyle=\pi
A\cos(\pi\frac{f(x,t)}{s})\sin(\pi\frac{y}{s})\frac{df}{dx}-\mu
y+\eta_{2}(t),$ (5b) $\displaystyle f(x,t)$
$\displaystyle=\varepsilon\sin(\omega t+\psi)x^{2}+(1-2\varepsilon\sin(\omega
t+\psi))x.$ (5c)

When $\varepsilon=0$, the double-gyre flow is time-independent, while for
$\varepsilon\neq 0$, the gyres undergo a periodic expansion and contraction in
the $x$ direction. In (5), $A$ approximately determines the amplitude of the
velocity vectors, $\omega/2\pi$ gives the oscillation frequency, $\varepsilon$
determines the amplitude of the left-right motion of the separatrix between
the gyres, $\psi$ is the phase, $\mu$ determines the dissipation, $s$ scales
the dimensions of the workspace, and $\eta_{i}(t)$ describes a stochastic
white noise with mean zero and standard deviation $\sigma=\sqrt{2I}$, for
noise intensity $I$. In this work, $\eta_{i}(t)$ can be viewed as either
measurement or environmental noise. Fig. 4 shows the phase portrait of the
time-independent double-gyre model.

Figure 4: Phase portrait of the model given by (5) with $A=10$, $\mu=0.005$,
$\varepsilon=0$, $\psi=0$, $I=0$, and $s=50$.

Fig. 5 shows the use of the control strategy (3) and (4) to track the
Lagrangian coherent structures of the periodic double-gyre model with noise.
As mentioned in Section I, LCS are extensions of stable and unstable manifolds
to non-autonomous dynamical systems [21]. We note that while the control
strategy was based on techniques developed for time-independent systems, the
method performs surprisingly well in tracking LCS for slow time-varying
systems in the presence of noise. Details regarding LCS computation can be
found in [5] and we refer the interested reader to the attached multimedia
file for a movie of the full simulation run. While the control strategy was
developed for static flows, the movie shows the robustness of the strategy for
tracking LCS in time-varying flows.

(a) t=0.8

(b) t=1.4

(c) t=1.8

(d) t=2.6

(e) t=3.0

(f) t=3.2

(g) t=3.6

(h) t=4.0

Figure 5: Trajectories of the team of $3$ robots tracking the Lagrangian
coherent structures of the system described by (5) with $A=10$, $\mu=0.005$,
$\varepsilon=0.1$, $\psi=0$, $I=0.01$, and $s=50$. The trajectories of the
straddling robots are shown in black and the estimated LCS is shown in white.

##  VII DISCUSSION

In this paper, we have designed a control strategy that allows collaborating
robots to track coherent structures and manifolds on general static
conservative flows. In addition, we showed how the strategy can be used to
track LCS in time-dependent conservative flows with measurement noise. The
saddle straddle control strategy is based on the communication of local
velocity field measurements obtained by each robot. Using the local velocity
field information provided by the two straddling robots (the herders), one
robot (the tracker), is able to detect the coherent structures, a global
structure that delineates the phase space into different dynamical regions.
Our work is novel in that the robots are determining the location of a global
structure based solely on local information, and as far as we know, the
sensing of LCS in the ocean has never been performed using autonomous
vehicles. Moreover, only initial state knowledge of the LCS is required
locally to get an accurate prediction of the global structure.

While the cooperative control strategy was inspired by the PIM Triple
Procedure, a procedure that relies on the computation of escape times which is
a global property of the system, the controller itself only relies on
information provided by each robot’s onboard sensors. We also note that the
realization of the control strategy by the team of robots can be achieved
without the need for global localization information. As such, the strategy is
a purely local strategy. Furthermore, the cooperative control strategy was
derived to track the manifolds on a static flow, but performs surprisingly
well at tracking the LCS in the time-dependent double-gyre model in the
presence of noise.

Since realistic quasi-geostrophic ocean models exhibit double-gyre flow
solutions, our first attempt seems to suggest that our methods may be and
general enough to be applied to more complicated models, including multi-layer
PDE ocean models. As mentioned in Section IV the robustness of the control
strategy is dependent on numerous parameters in the system which includes
robots’ sensing and communication ranges, the bounds on the sensor and
actuation noise, the vector interpolation technique, the sampling frequency,
and the relative time scales of the AUV dynamics in relation to the
surrounding flow dynamics. While our initial results suggest that our approach
may be robust enough to measurement noise, a more thorough understanding of
the sensitivity of the proposed strategy to these various parameters is
instrumental in extending our approach to more realistic ocean models and for
field deployment.

##  VIII FUTURE WORK

In recent years, there has been significant interest in the use of AUVs to
collect scientific data in the ocean to improve our ability to forecast
harmful algae blooms and weather and climate patterns. One drawback to
operating sensors in time-dependent and stochastic environments like the ocean
is that the sensors will tend to escape from their monitoring region of
interest. As such, the ability to identify and track Lagrangian coherent
structures (LCS) in these dynamic environments is paramount in maintaining
appropriate sensor coverage in regions of interest. Additionally, since LCS
have been shown to coincide with optimal trajectories in the ocean which
minimize the energy and navigation time [3, 4], real-time knowledge of these
“super-highways” is key in planning efficient AUVs paths.

Of particular interest is the extension of our method to more realistic ocean
models. Specifically, can we extend our current cooperative tracking strategy
to a swarm of heterogenous mobile and stationary sensors? By increasing the
team size and incorporating both stationary and mobile sensing devices, it is
possible to refine our tracking technique to reveal the coherent structures at
various spatial and time scales. One immediate direction for future work is to
investigate how the proposed strategy scales to larger team sizes. Second,
underwater environments pose unique challenges in terms of wireless
communications. In general, acoustic transmissions generally have low data
rates and acoustic wave propagation can be further affected by the surrounding
fluid dynamics [22]. As such, a second direction for future work is to
investigate how communication delays and missed transmissions impact the
overall accuracy of the tracking methodology. In this work, we assume an
initial state knowledge of the LCS is required. This initial formation may be
difficult to achieve without any prior global knowledge of the flows. By
considering a team of both stationary and mobile sensors, one can potentially
obtain an initial estimate of a local LCS through the stationary sensing
network which can then be tracked and further refined by the mobile nodes. A
third direction for future work is to determine how one can strategically
place a combination of mobile and stationary sensors to provide real-time
updates on the locations of LCS.

## References

  * [1] G. Haller and G. Yuan, “Lagrangian coherent structures and mixing in two-dimensional turbulence,” _Phys. D_ , vol. 147, pp. 352–370, Dec 2000\.
  * [2] S. C. Shadden, F. Lekien, and J. E. Marsden, “Definition and properties of lagrangian coherent structures from finite-time lyapunov exponents in two-dimensional aperiodic flows,” _Physica D: Nonlinear Phenomena_ , vol. 212, no. 3-4, pp. 271 – 304, 2005.
  * [3] T. Inanc, S. Shadden, and J. Marsden, “Optimal trajectory generation in ocean flows,” in _American Control Conference, 2005. Proceedings of the 2005_ , 8-10, 2005, pp. 674 – 679.
  * [4] C. Senatore and S. Ross, “Fuel-efficient navigation in complex flows,” in _American Control Conference, 2008_ , june 2008, pp. 1244 –1248.
  * [5] E. Forgoston, L. Billings, P. Yecko, and I. B. Schwartz, “Set-based corral control in stochastic dynamical systems: Making almost invariant sets more invariant,” _Chaos_ , vol. 21, no. 013116, 2011.
  * [6] C. Hsieh, Z. Jin, D. Marthaler, B. Nguyen, D. Tung, A. Bertozzi, and R. Murray, “Experimental validation of an algorithm for cooperative boundary tracking,” in _Proc. of the 2005 American Control Conference._ , 2005, pp. 1078–1083.
  * [7] S. Susca, S. Martinez, and F. Bullo, “Monitoring environmental boundaries with a robotic sensor network,” _IEEE Trans. on Control Systems Technology_ , vol. 16, no. 2, pp. 288–296, 2008.
  * [8] I. Triandaf and I. B. Schwartz, “A collective motion algorithm for tracking time-dependent boundaries,” _Mathematics and Computers in Simulation_ , vol. 70, pp. 187–202, 2005.
  * [9] V. M. Goncalves, L. C. A. Pimenta, C. A. Maia, B. Dutra, and G. A. S. Pereira, “Vector fields for robot navigation along time-varying curves in n-dimensions,” _IEEE Trans. on Robotics_ , vol. 26, no. 4, pp. 647–659, 2010\.
  * [10] F. Zhang, D. M. Fratantoni, D. Paley, J. Lund, and N. E. Leonard, “Control of coordinated patterns for ocean sampling,” _Int. Journal of Control_ , vol. 80, no. 7, pp. 1186–1199, 2007.
  * [11] K. M. Lynch, P. Schwartz, I. B. Yang, and R. A. Freeman, “Decentralized environmental modeling by mobile sensor networks,” _IEEE Trans. Robotics_ , vol. 24, no. 3, pp. 710–724, 2008.
  * [12] W. Wu and F. Zhang, “Cooperative exploration of level surfaces of three dimensional scalar fields,” _Automatica, the IFAC Journall_ , vol. 47, no. 9, pp. 2044–2051, 2011.
  * [13] J. Spletzer and R. Fierro, “Optimal positioning strategies for shape changes in robot teams,” in _Proc. of the IEEE Int. Conf. on Robotics & Automation_, Barcelona, Spain, April 2005, pp. 754–759.
  * [14] S. Kalantar and U. R. Zimmer, “Distributed shape control of homogeneous swarms of autonomous underwater vehicles,” _to appear in Autonomous Robots (intl. journal)_ , 2006.
  * [15] M. A. Hsieh, S. Loizou, and V. Kumar, “Stabilization of multiple robots on stable orbits via local sensing,” in _Proc. of the Int. Conf. on Robotics & Automation (ICRA) 2007_, Rome, Italy, April 2007.
  * [16] H. E. Nusse and J. A. Yorke, “A procedure for finding numerical trajectories on chaotic saddles,” _Physica D Nonlinear Phenomena_ , vol. 36, pp. 137–156, jun 1989.
  * [17] J. C. Agui and J. Jimenez, “On the performance of particle tracking,” _Journal of Fluid Mechanics_ , vol. 185, pp. 447–468, 1987. [Online]. Available: http://dx.doi.org/10.1017/S0022112087003252
  * [18] C. Marchioli, V. Armenio, and A. Soldati, “Simple and accurate scheme for fluid velocity interpolation for eulerianlagrangian computation of dispersed flows in 3d curvilinear grids,” _Computers & Fluids_, vol. 36, pp. 1187–1198, 2007.
  * [19] E. J. Fuselier and G. B. Wright, “Stability and error estimates for vector field interpolation and decomposition on the sphere with rbfs,” _SIAM J. Numer. Anal._ , vol. 47, pp. 3213–3239, 2009.
  * [20] P. Hartman, _Ordinary Differential Equations_ , 2nd ed. Society for Industrial & Applied Math, 2002.
  * [21] S. Kent, “Lagrangian coherent structures: Generalizing stable and unstable manifolds to non-autonomous dynamical systems,” Tucson, AZ, 2008.
  * [22] D. Wang, P. F. Lermusiaux, P. J. Haley, D. Eickstedt, W. G. Leslie, and H. Schmidt, “Acoustically focused adaptive sampling and on-board routing for marine rapid environmental assessment,” _Journal of Marine Systems_ , vol. 78, pp. 393–407, 2009.