Quasi-Bi-Hamiltonian Structures of the 2-Dimensional Kepler Problem

Quasi-Bi-Hamiltonian Structures
of the 2-Dimensional Kepler Problem\star\star$\star$This paper is a
contribution to the Special Issue on Analytical Mechanics and Differential
Geometry in honour of Sergio Benenti. The full collection is available at
http://www.emis.de/journals/SIGMA/Benenti.html

Jose F. CARIÑENA and Manuel F. RAÑADA J.F. Cariñena and M.F. Rañada
Departamento de Física Teórica and IUMA, Universidad de Zaragoza, 50009
Zaragoza, Spain<EMAIL_ADDRESS><EMAIL_ADDRESS>
Received September 29, 2015, in final form January 25, 2016; Published online
January 27, 2016

The existence of quasi-bi-Hamiltonian structures for the Kepler problem is
studied. We first relate the superintegrability of the system with the
existence of two complex functions endowed with very interesting Poisson
bracket properties and then we prove the existence of a quasi-bi-Hamiltonian
structure by making use of these two functions. The paper can be considered as
divided in two parts. In the first part a quasi-bi-Hamiltonian structure is
obtained by making use of polar coordinates and in the second part a new
quasi-bi-Hamiltonian structure is obtained by making use of the separability
of the system in parabolic coordinates.

Kepler problem; superintegrability; complex structures; bi-Hamiltonian
structures; quasi-bi-Hamiltonian structures

37J15; 37J35; 70H06; 70H33

##  1 Introduction

###  1.1 A geometric introduction

In differential geometric terms, the phase space $M$ of a Hamiltonian system
is the $2n$-dimensional cotangent bundle $M=T^{*}Q$ of the $n$-dimensional
configuration space $Q$. Cotangent bundles are manifolds endowed, in a natural
or canonical way, with a symplectic structure $\omega_{0}$; if
$\\{(q_{i})\,|\,i=2,\ldots,n\\}$ are local coordinates in $Q$ and
$\\{(q_{j},p_{j});\,j=1,2,\dots,n\\}$ the induced coordinates in $T^{*}Q$,
then $\omega_{0}$ is given by

$\displaystyle\omega_{0}=dq_{j}\wedge
dp_{j},\qquad\omega_{0}=-d{\theta_{0}},\qquad\theta_{0}=p_{j}dq_{j},$

(we write all the indices as subscripts and summation convention on the
repeated index is used). Given a differentiable function $F$ in $T^{*}Q$,
$F=F(q,p)$, the vector field $X_{F}$ uniquely defined as the solution of the
equation

$\displaystyle i(X_{F})\omega_{0}=dF$

is called the Hamiltonian vector field of the function $F$. In particular,
given a Hamiltonian $H=H(q,p)$, the dynamics is given by the corresponding
Hamiltonian vector field $X_{H}$, that is, $i(X_{H})\omega_{0}=dH$.

A vector field $\Gamma\in\mathfrak{X}(T^{*}Q)$ is Hamiltonian if there is a
function $H$ such that $\Gamma=X_{H}$, i.e., $i(\Gamma)\omega_{0}=dH$, and
locally-Hamiltonian when $i(\Gamma)\omega_{0}$ is a closed 1-form. This is
equivalent to $\Gamma$-invariance of $\omega_{0}$, i.e.,
$\mathcal{L}_{\Gamma}\omega_{0}=0$.

A system of differential equations is called bi-Hamiltonian if it can be
written in two different ways in Hamiltonian form. Suppose a manifold $M$
equipped with two different symplectic structures $\omega_{0}$ and
$\omega_{1}$. A vector field $\Gamma$ on $T^{*}Q$ is said to be a bi-
Hamiltonian vector field if it is Hamiltonian with respect to both symplectic
structures, i.e.,

$\displaystyle i(\Gamma)\omega_{0}=dH_{0}\qquad{\rm and}\qquad
i(\Gamma)\omega_{1}=dH_{1}.$

The two functions, $H_{0}$ and $H_{1}$, are integrals of motion for $\Gamma$.
A weaker form of bi-Hamiltonian system is when the only symplectic form is the
first one ($\omega_{1}$ is a closed but nonsymplectic 2-form).

We point out that an important example of bi-Hamiltonian system is the
rational harmonic oscillator (non-central harmonic oscillator with rational
ratio of frequencies) [8]

$\displaystyle
H=\frac{1}{2}\bigl{(}p_{x}^{2}+p_{y}^{2}\bigr{)}+\frac{1}{2}\alpha_{0}^{2}\bigl{(}m^{2}x^{2}+n^{2}y^{2}\bigr{)}.$

A symplectic form determines a Poisson bivector $\Lambda$ that satisfies the
vanishing of the Schouten bracket $[\Lambda,\Lambda]=0$ (this property is
equivalent to the Jacobi identity); we note that there exist non-constant rank
Poisson structure not related with symplectic structures. The compatibility
condition between two different Poisson structures, $\Lambda_{0}$ and
$\Lambda_{1}$, means that the linear combination
$\Lambda_{\lambda}=\Lambda_{0}-{\lambda}\Lambda_{1}$ is a Poisson pencil, that
is, it is a Poisson bivector for every value of ${\lambda}$; therefore the
corresponding bracket
$\\{\cdot,\cdot\\}_{\lambda}=\\{\cdot,\cdot\\}_{0}-{\lambda}\\{\cdot,\cdot\\}_{1}$
is a pencil of Poisson brackets.

Bi-Hamiltonian systems are systems endowed with very interesting properties
but, in general, it is quite difficult to find a bi-Hamiltonian formulation
for a given Hamiltonian vector field, and for this reason it is useful to
introduce the concept of quasi-bi-Hamiltonian system.

######  Definition 1.1.

A vector field $X$ on a symplectic manifold $(M,\omega)$ is called quasi-
Hamiltonian if there exists a (nowhere-vanishing) function $\mu$ such that
$\mu X$ is a Hamiltonian vector field

$\displaystyle\mu X\in{\mathfrak{X}}_{H}(M).$

Thus $i(\mu X)\omega=dh$ for some function $h$.

We call $\mu$ an integrating factor of the quasi-Hamiltonian system, because
it is an integrating factor for the 1-form $i(X)\omega$, and we note that in
this case the function $h$ is a first integral of $X$. Note that this
condition can alternatively be written as as $i(X)(\mu\omega)=dh$, but the
point is that the 2-form $\mu\omega$ is not closed in the general case.

The scarcity of bi-Hamiltonian systems leads to relax the condition for the
vector field being Hamiltonian to a simpler situation of quasi-Hamiltonian
with respect to the second symplectic structure.

######  Definition 1.2.

A Hamiltonian vector field $X$ on a symplectic manifold $(M,\omega)$ is called
quasi-bi-Hamiltonian if there exist another symplectic structure $\omega_{1}$,
and a nowhere-vanishing function $\mu$, such that $\mu X$ is a Hamiltonian
vector field with respect to $\omega_{1}$.

This concept was first introduced in [4] in the particular case of systems
with two degrees of freedom and it was quickly extended in [23, 24] for a
higher-dimensional systems. Some recent papers considering properties of this
particular class of systems are [1, 2, 3, 4, 5, 6, 8, 14, 15, 23, 24, 25, 36].

The nondegeneracy of the canonical form $\omega_{0}$ provides a vector bundle
isomorphism $\widehat{\omega_{0}}$ of $T(T^{*}Q)$ on $T^{*}(T^{*}Q)$, inducing
an identification of vector fields and 1-forms on the phase space. A
consequence is that the pair $(\omega_{0},\omega_{1})$ determines a $(1,1)$
tensor field $R$ defined as
$R=\widehat{\omega_{0}}^{-1}\circ\widehat{\omega_{1}}$, that is,

$\displaystyle\omega_{1}(X,Y)=\omega_{0}(RX,Y),\qquad\forall\,X,Y\in\mathfrak{X}(T^{*}Q).$

Note that in the definition of $R$ the 2-form $\omega_{1}$ is necessarily
closed but it can be nonsymplectic. If $\Gamma$ is bi-Hamiltonian with respect
to $(\omega_{0},\omega_{1})$ then $R$ is $\Gamma$-invariant, that is,
${\mathcal{L}}_{\Gamma}R=0$, where ${\mathcal{L}}$ denotes the Lie derivative
(this means that the characteristic polynomial of $R$ is an invariant for
$\Gamma$, and consequently the coefficients of the polynomial are constants of
motion). The Nijenhuis tensor $N_{R}$ of the tensor field $R$ is defined by

$\displaystyle N_{R}(X,Y)=R^{2}[X,Y]+[RX,RY]-R[RX,Y]-R[X,RY].$

It has been proved that if $\Gamma$ is a Hamiltonian dynamical system of the
type described above and such that (i) The tensor $N_{R}$ of $R$ vanishes,
(ii) The tensor field $R$ has $n$ distinct eigenfunctions (that is, they are
maximally distinct), then the eigenfunctions of $R$ are in involution and the
system is therefore completely integrable [12, 16, 17, 18, 21]. It is
important to note that the eigenvalues of $R$ are constants of motion for
$\Gamma$ even in the case that the two properties (i) and (ii) are not
satisfied (but then the eigenfunctions are not in involution).

It is known that the Liouville formalism characterize the Hamiltonians that
are integrable but it does not provide a method for obtaining the integrals of
motion; therefore it has been necessary to elaborate different methods for
obtaining constants of motion (Hamilton–Jacobi separability, Lax pairs
formalism, Noether symmetries, Hidden symmetries, etc); the existence of a bi-
Hamiltonian structure with the above two mentioned properties (Nijenhuis
torsion condition and maximally distinct eigenvalues) can be considered as
method to establish the Liouville integrability of a system; because of this,
these two properties are frequently included in the definition of a bi-
Hamiltonian structure.

Most of systems admitting a bi-Hamiltonian structure are separable systems; so
these two properties (separability and double Hamiltonian structure) are
properties very close related (see [18] for a detailed discussion of this
question and [35] for the case of multiple separability). Quasi-bi-Hamiltonian
systems are very less known than the bi-Hamiltonian ones but it seems that
they are also related with separability. Nevertheless, in this case the tensor
field $R$ is not $\Gamma$-invariant and the eigenvalues of the tensor field
$R$ are not constants of motion. An interesting property is that it was shown
in [4] that in the particular case of two degrees of freedom the function
$\mu^{2}/\det R$ is a constant of the motion.

The potential of the Kepler problem is spherically symmetric and therefore it
admits alternative Lagrangians (the existence of alternative Lagrangians for
central potentials is studied in [13, 20, 26]); recall that if there exist
alternative Lagrangian descriptions then one can find non-Noether constants of
motion [7]. This system has been studied as a bi-Hamilltonian system by making
use of different approaches; Rauch-Wojciechowski proved the existence of a bi-
Hamiltonian formulation but introducing an extra variable so that the phase
space is odd-dimensional and the Poisson brackets are degenerate [34] and more
recently [19] a bi-Hamiltonian formulation for the perturbed Kepler problem
has also been studied by making use of Delaunay-type variables.

###  1.2 Structure and purpose of the paper

Now in this paper we will analyze certain properties of the Kepler problem
related with the existence of quasi-bi-Hamiltonian structures. The following
two points summarize the contents of the paper.

  * •

First, we will study the existence of certain complex functions with
interesting Poisson properties and then we will prove that the
superintegrability of the system is directly related with the properties of
these complex functions. As is well-known this system is multi-separable (it
separates in both polar and parabolic coordinates); so we first present the
study making use of polar coordinates and then we undertake a similar study by
employing parabolic coordinates (the parabolic complex functions are different
from the polar ones).

  * •

Second, we prove that the above mentioned complex functions determine the
existence of several quasi-bi-Hamiltonian structures. This is done in two
steps: first with complex 2-forms and then with several real 2-forms. The
properties of these geometric structures and of the associated recursion
operators are analyzed.

It is important to note that this study is concerned with the existence of
quasi-bi-Hamiltonian structures (instead of bi-Hamiltonian). So we recall that
if $i(\Gamma)\omega_{1}={\lambda}dH_{1}$ then
${\mathcal{X}}_{\Gamma}\omega_{1}=d{\lambda}\wedge dH_{1}\neq 0$. Consequently
the tensor field $R$ is not $\Gamma$-invariant and the eigenvalues of $R$ are
not constants of motion.

We must clearly say that the structures obtained by this method (wedge product
of the differentials of complex functions) do not satisfy the above mentioned
Nijenhuis torsion condition. So perhaps it is convenient to name them as weak
quasi-bi-Hamiltonian structures (in opposition to strong structures satisfying
the Nijenhuis condition). Nevertheless, the purpose in this paper is not to
prove the integrability of a system as consequence of a bi-Hamiltonian
structure, since it is perfectly known that the Kepler problem is not only
integrable but also superintegrable. The purpose is to study new and
interesting properties of the Kepler problem. In fact, it has been proved that
if a dynamical vector field satisfies certain properties (existence of
canonoid transformations [9, 10] or existence of non-symplectic symmetries
[27, 28]) then it is Hamiltonian with respect to two different structures
without satisfying necessarily the Nijenhuis condition.

The structure of the paper is as follows: In Section 2 we study the Kepler
problem by making use of polar coordinates $(r,\phi)$; we relate the
superintegrability of the system with the existence of two complex functions
$M_{r}$ and $N_{\phi}$ endowed with very interesting Poisson bracket
properties and then we prove the existence of a quasi-bi-Hamiltonian structure
making use of these two functions. Then in Section 3 we consider once more the
same system but in terms of parabolic coordinates $(a,b)$ and we obtain a new
quasi-bi-Hamiltonian structure (different to the previous one) making use of a
similar technique but with new complex functions $M_{a}$ and $M_{b}$. Finally
in Section 4 we make some final comments.

##  2 From superintegrability to quasi-bi-Hamiltonian structures

After these rather general comments we restrict our study to the Kepler
problem in the Euclidean plane that, as it is well known, is superintegrable
and multiseparable (polar and parabolic coordinates).

Let us first notice that in some cases the two-dimensional Euclidean systems
possess certain interesting properties. For example, if the potential $V(x,y)$
takes the form $V=A(u)+B(v)$, $u=x+y$, $v=x-y$, then it admits a new
Hamiltonian structure (and also a Lax pair) [22]; unfortunately the new
structure is in most of cases constant (we mean that the new Poisson bracket
is determined by a symplectic form as $\omega_{1}=dx\wedge dp_{y}+dy\wedge
dp_{x}$).

In the next paragraphs we will relate the existence of bi-Hamiltonian
structures with the properties of the two complex functions with interesting
Poisson bracket properties.

###  2.1 Complex functions and superintegrabilty

It is well known that the Hamiltonian of the two-dimensional Kepler problem

$\displaystyle
H_{K}=\frac{1}{2}\left(p_{r}^{2}+\frac{p_{\phi}^{2}}{r^{2}}\right)+V_{K},\qquad
V_{K}=-\frac{g}{r},\qquad 0<g\in\mathbb{R},$

is Hamilton–Jacobi (H-J) separable in polar coordinates $(r,\phi)$ and it is,
therefore, Liouville integrable with the angular momentum $J_{2}=p_{\phi}$ as
the associate constant of motion. Moreover, it is also known that it is a
super-integrable system with the two components of the Laplace–Runge–Lenz
vector as additional integrals of motion. Now we will prove that this property
of superintegrability can be related to the existence of certain complex
functions with interesting Poisson bracket properties.

Let us denote by $M_{rj}$ and $N_{\phi j}$, $j=1,2$, the following real
functions

$\displaystyle M_{r1}=p_{r}p_{\phi},\qquad
M_{r2}=g-\frac{p_{\phi}^{2}}{r},\qquad\text{and}\qquad N_{\phi
1}=\cos\phi,\qquad N_{\phi 2}=\sin\phi.$

Then we have the following properties

(i) $\displaystyle\frac{d}{dt}M_{r1}=\\{M_{r1},H_{K}\\}=-\lambda
M_{r2},\qquad$ $\displaystyle\frac{d}{dt}M_{r2}=\\{M_{r2},H_{K}\\}=\lambda
M_{r1},$ (ii) $\displaystyle\frac{d}{dt}N_{\phi 1}=\\{N_{\phi
1},H_{K}\\}=-\lambda N_{\phi 2},\qquad$ $\displaystyle\frac{d}{dt}N_{\phi
2}=\\{N_{\phi 2},H_{K}\\}=\lambda N_{\phi 1},$

where $\lambda$ denote the following function

$\displaystyle\lambda=\frac{p_{\phi}}{r^{2}}.$

The property (ii), representing the behaviour of the angular functions
$N_{\phi j}$, is true for all the central potentials $V(r)$; but the property
(i), behaviour of the functions $M_{rj}$, is specific of the potential of the
Kepler problem.

Consider next the complexification of the linear space of functions on the
manifold and extend by bilinearity the Poisson bracket. If we denote $M_{r}$
and $N_{\phi}$ the complex functions

$\displaystyle M_{r}=M_{r1}+iM_{r2},\qquad N_{\phi}=N_{\phi 1}+iN_{\phi 2},$

then they have the following properties

$\displaystyle\\{M_{r},H_{K}\\}=i\lambda
M_{r},\qquad\\{N_{\phi},H_{K}\\}=i\lambda N_{\phi},$

and consequently the Poisson bracket of $M_{r}N_{\phi}^{*}$ with the Kepler
Hamiltonian vanishes

$\displaystyle\\{M_{r}N_{\phi}^{*},H_{K}\\}=\\{M_{r},H_{K}\\}N_{\phi}^{*}+M_{r}\\{N_{\phi}^{*},H_{K}\\}=(i\lambda
M_{r})N_{\phi}^{*}+M_{r}(-i\lambda N_{\phi}^{*})=0.$

We can summarize this result in the following proposition.

######  Proposition 2.1.

Let us consider the Hamiltonian of the Kepler problem

$\displaystyle
H_{K}=\frac{1}{2}\left(p_{r}^{2}+\frac{p_{\phi}^{2}}{r^{2}}\right)+V_{K},\qquad
V_{K}=-\frac{g}{r},$

Then, the complex function $J_{34}$ defined as

$\displaystyle J_{34}=M_{r}N_{\phi}^{*}$

is a $($complex$)$ constant of the motion.

Of course $J_{34}$ determines two real first-integrals

$\displaystyle
J_{34}=J_{3}+iJ_{4},\qquad\\{J_{3},H_{K}\\}=0,\qquad\\{J_{4},H_{K}\\}=0,$

whose coordinate expressions turn out to be

$\displaystyle J_{3}={\rm
Re}(J_{34})=p_{r}p_{\phi}\cos\phi-\frac{p_{\phi}^{2}}{r}\sin\phi+g\sin\phi,$
$\displaystyle J_{4}={\rm
Im}(J_{34})=p_{r}p_{\phi}\sin\phi+\frac{p_{\phi}^{2}}{r}\cos\phi-g\cos\phi.$

That is, the two functions $J_{3}$ and $J_{4}$ are just the two components of
the two-dimensional Laplace–Runge–Lenz vector.

Summarizing, we have got two interesting properties. First, the
superintegrability of the Kepler problem is directly related with the
existence of two complex functions whose Poisson brackets with the Hamiltonian
are proportional with a common complex factor to themselves, and second, the
two components of the Laplace–Runge–Lenz vector appear as the real and
imaginary parts of the complex first-integral of motion. Remark that
$N_{\phi}$ is a complex function of constant modulus one, while the modulus of
$M_{r}$ is a polynomial of degree four in the momenta given by

$\displaystyle
M_{r}M_{r}^{*}=(p_{r}p_{\phi})^{2}+\left(g-\frac{p_{\phi}^{2}}{r}\right)^{2}=2p_{\phi}^{2}H_{K}+g^{2}.$

###  2.2 Complex functions and quasi-bi-Hamiltonian structures

Let us denote by $Y_{34}$ the (complex) Hamiltonian vector field of $J_{34}$

$\displaystyle i(Y_{34})\omega_{0}=dJ_{34},$

that obviously satisfies $Y_{34}(H_{K})=\\{H_{K},J_{34}\\}=0$, and by $Y_{r}$
and $Y_{\phi}$ the Hamiltonian vector fields of $M_{r}$ and $N_{\phi}$:

$\displaystyle i(Y_{r})\omega_{0}=dM_{r},\qquad
i(Y_{\phi})\omega_{0}=dN_{\phi}.$

Their local coordinate expressions are, respectively, given by

$\displaystyle Y_{r}=p_{\phi}\frac{\partial}{\partial
r}+\left(p_{r}-2i\frac{p_{\phi}}{r}\right)\frac{\partial}{\partial\phi}-i\left(\frac{p_{\phi}^{2}}{r^{2}}\right)\frac{\partial}{\partial
p_{r}},$

and

$\displaystyle Y_{\phi}=(\sin\phi-i\cos\phi)\frac{\partial}{\partial
p_{\phi}}.$

Then, the vector field $Y_{34}$ appears as a linear combination of $Y_{r}$ and
$Y_{\phi}^{*}$; more specifically we have

$\displaystyle Y_{34}=N_{\phi}^{*}Y_{r}+M_{r}Y_{\phi}^{*}=Y+Y^{\prime},\qquad
Y=N_{\phi}^{*}Y_{r},\qquad Y^{\prime}=M_{r}Y_{\phi}^{*}.$

The vector field $Y_{34}$ is certainly a symmetry of the Hamiltonian system
$(T^{*}Q,\omega_{0},H_{K})$, but the two vector fields, $Y$ and $Y^{\prime}$,
are neither symmetries of the symplectic form $\omega_{0}$ (that is,
${\mathcal{X}}_{Y}\omega_{0}\neq 0$ and
${\mathcal{X}}_{Y^{\prime}}\omega_{0}\neq 0$) nor symmetries of the
Hamiltonian (that is, ${\mathcal{X}}_{Y}H_{K}\neq 0$ and
${\mathcal{X}}_{Y^{\prime}}H_{K}\neq 0$). Moreover, remark that they are not
symmetries of the dynamics, because

$\displaystyle[Y,\Gamma_{K}]\neq 0,\qquad[Y^{\prime},\Gamma_{K}]\neq 0,\qquad
i(\Gamma_{K})\omega_{0}=dH_{K}.$

More specifically:

######  Proposition 2.2.

The Lie bracket of the dynamical vector field $\Gamma_{K}$ with $Y$ is given
by

$\displaystyle[\Gamma_{K},Y]=iJ_{34}X_{\lambda},$

where $X_{\lambda}$ is the Hamiltonian vector field of the function $\lambda$.

######  Proof 2.3.

A direct computation leads to

$\displaystyle[\Gamma_{K},Y]=\Gamma_{K}(N_{\phi}^{*})Y_{r}+N_{\phi}^{*}[\Gamma_{K},Y_{r}]=-i{\lambda}N_{\phi}^{*}Y_{r}+N_{\phi}^{*}\bigl{(}{-}X_{\\{H_{K},M_{r}\\}}\bigr{)}.$

where we have used that the Lie bracket of two Hamiltonian vector fields
satisfies $[X_{f},X_{g}]=-X_{\\{f,g\\}}$. Note also that the Hamiltonian
vector field of a product $fg$ is given by $X_{fg}=fX_{g}+gX_{f}$, and then
the above Lie bracket becomes

$\displaystyle[\Gamma_{K},Y]=-i{\lambda}N_{\phi}^{*}Y_{r}+iN_{\phi}^{*}(X_{{\lambda}M_{r}})=-i{\lambda}N_{\phi}^{*}Y_{r}+iN_{\phi}^{*}({\lambda}Y_{r}+M_{r}X_{\lambda})=i(M_{r}N_{\phi}^{*})X_{\lambda}.$
∎

The vector field $X_{\lambda}$ on the right hand side represents an
obstruction for $Y$ to be a dynamical symmetry. Only when ${\lambda}$ be a
numerical constant the vector field $Y$ (and also $Y^{\prime}$) is a dynamical
symmetry of $\Gamma_{K}$.

In the following $\Omega$ will denote the complex 2-form defined as

$\displaystyle\Omega=dM_{r}\wedge dN_{\phi}^{*}.$

The two complex 2-forms $\omega_{Y}$ and $\omega_{Y}^{\prime}$ obtained by Lie
derivative of $\omega_{0}$, i.e.,

$\displaystyle{\mathcal{L}}_{Y}\omega_{0}=\omega_{Y},\qquad{\mathcal{L}}_{Y^{\prime}}\omega_{0}=\omega_{Y}^{\prime},$

are such

$\displaystyle{\mathcal{L}}_{Y}\omega_{0}=i_{Y}(d\omega_{0})+d(i_{Y}\omega_{0})=d(i_{Y}\omega_{0})=d(N_{\phi}^{*}dM_{r})=-\Omega,$
$\displaystyle{\mathcal{L}}_{Y^{\prime}}\omega_{0}=i_{Y^{\prime}}(d\omega_{0})+d(i_{Y^{\prime}}\omega_{0})=d(i_{Y^{\prime}}\omega_{0})=d(M_{r}dN_{\phi}^{*})=\Omega.$

Using the preceding results we can prove:

######  Proposition 2.4.

The Hamiltonian vector field $\Gamma_{K}$ of the Kepler problem is a quasi-
Hamiltonian system with respect to the complex $2$-form $\Omega$.

######  Proof 2.5.

The contraction of the vector field $\Gamma_{K}$ with the complex 2-form
$\Omega$ gives

$\displaystyle
i(\Gamma_{K})\Omega=\Gamma_{K}(M_{r})dN_{\phi}^{*}-\Gamma_{K}(N_{\phi}^{*})dM_{r},$

and recalling that

$\displaystyle\Gamma_{K}(M_{r})=\\{M_{r},H_{K}\\}=i\lambda
M_{r},\qquad\Gamma_{K}(N_{\phi}^{*})=\\{N_{\phi}^{*},H_{K}\\}=-i\lambda
N_{\phi}^{*},$

we arrive to

$\displaystyle i(\Gamma_{K})\Omega=(i\lambda M_{r})dN_{\phi}^{*}+(i\lambda
N_{\phi}^{*})dM_{r}=i\lambda d(M_{r}N_{\phi}^{*}).$ ∎

The complex 2-form $\Omega$ can be written as

$\displaystyle\Omega=\Omega_{1}+i\Omega_{2}$

where the two real 2-forms, $\Omega_{1}={\rm Re}(\Omega)$ and $\Omega_{2}={\rm
Im}(\Omega)$, take the form

$\displaystyle\Omega_{1}=dM_{r1}\wedge dN_{\phi 1}+dM_{r2}\wedge dN_{\phi
2}=d(p_{r}p_{\phi})\wedge
d(\cos\phi)+d\left(g-\frac{p_{\phi}^{2}}{r}\right)\wedge d(\sin\phi)$
$\displaystyle\hphantom{\Omega_{1}}{}=\alpha_{12}dr\wedge
d\phi+\alpha_{23}d\phi\wedge dp_{r}+\alpha_{24}d\phi\wedge dp_{\phi},$
$\displaystyle\Omega_{2}=-dM_{r1}\wedge dN_{\phi 2}+dM_{r2}\wedge dN_{\phi
1}=-d(p_{r}p_{\phi})\wedge
d(\sin\phi)+d\left(g-\frac{p_{\phi}^{2}}{r}\right)\wedge d(\cos\phi)$
$\displaystyle\hphantom{\Omega_{2}}{}=\beta_{12}dr\wedge
d\phi+\beta_{23}d\phi\wedge dp_{r}+\beta_{24}d\phi\wedge dp_{\phi}$

with $\alpha_{ij}$ and $\beta_{ij}$ being given by

$\displaystyle\alpha_{12}=\left(\frac{p_{\phi}^{2}}{r^{2}}\right)\cos\phi,\qquad\alpha_{23}=p_{\phi}\sin\phi,\qquad\alpha_{24}=p_{r}\sin\phi+2\left(\frac{p_{\phi}}{r}\right)\cos\phi,$

and

$\displaystyle\beta_{12}=-\left(\frac{p_{\phi}^{2}}{r^{2}}\right)\sin\phi,\qquad\beta_{23}=p_{\phi}\cos\phi,\qquad\beta_{24}=p_{r}\cos\phi-2\left(\frac{p_{\phi}}{r}\right)\sin\phi.$

Then we have

$\displaystyle i(\Gamma_{K})\Omega_{1}=-\lambda dJ_{4},\qquad
i(\Gamma_{K})\Omega_{2}=\lambda dJ_{3},$

what means that $\Gamma_{K}$ is also quasi-bi-Hamiltonian with respect to the
two real 2-forms $(\omega_{0},\Omega_{1})$ or $(\omega_{0},\Omega_{2})$.

Remark that the complex 2-form $\Omega$ is well defined but it is not
symplectic. In fact, from the above expressions in coordinates we have
$\Omega_{1}\wedge\Omega_{1}=0$, $\Omega_{2}\wedge\Omega_{2}=0$, and
$\Omega_{1}\wedge\Omega_{2}=0$, and therefore we obtain

$\displaystyle\Omega\wedge\Omega=(\Omega_{1}\wedge\Omega_{1}-\Omega_{2}\wedge\Omega_{2})+2i\Omega_{1}\wedge\Omega_{2}=0.$

The distribution defined by the kernel of $\Omega$, that is two-dimensional,
is given by

$\displaystyle\operatorname{Ker}\Omega=\big{\\{}f_{1}Z_{1}+f_{2}Z_{2}\,|\,f_{1},f_{2}\colon\mathbb{R}^{2}\times\mathbb{R}^{2}\to\mathbb{C}\big{\\}},$

where the vector fields $Z_{1}$ and $Z_{2}$ are

$\displaystyle Z_{1}=(\alpha_{23}+i\beta_{23})\frac{\partial}{\partial
r}+(\alpha_{12}+i\beta_{12})\frac{\partial}{\partial p_{r}},\qquad
Z_{2}=(\alpha_{24}+i\beta_{24})\frac{\partial}{\partial
r}+(\alpha_{12}+i\beta_{12})\frac{\partial}{\partial p_{\phi}}.$

Therefore it satisfies

$\displaystyle[\operatorname{Ker}\Omega,\Gamma_{K}]\subset\operatorname{Ker}\Omega.$

That is, $\Gamma_{K}$ preserves the distribution $\operatorname{Ker}\Omega$.

If $Y_{3}$ and $Y_{4}$ are the Hamiltonian vector fields (with respect to the
canonical symplectic form $\omega_{0}$) of the first integrals $J_{3}$ and
$J_{4}$, then the dynamical vector field $\Gamma_{K}$ is orthogonal to $Y_{4}$
with respect to the structure $\Omega_{1}$ and it is also orthogonal to
$Y_{3}$ with respect to the structure $\Omega_{2}$, that is,

$\displaystyle i(\Gamma_{K})i(Y_{4})\Omega_{1}=0,\qquad
i(\Gamma_{K})i(Y_{3})\Omega_{2}=0.$

Just to close the section we remark that had we applied this technique to the
isotropic two-dimensional harmonic oscillator with frequency $\alpha$ we had
obtained the function $M_{r}$ as

$\displaystyle
M_{r}=\left(\frac{2}{r}p_{r}p_{\phi}\right)+i\left(p_{r}^{2}-\frac{p_{\phi}^{2}}{r^{2}}+\alpha^{2}r^{2}\right),$

(the angular function $N_{\phi}$ would be the same) and the constants so
obtained are but the components of the Fradkin tensor. This shows that the
harmonic oscillator is an example of dynamical system both bi-Hamiltonian and
quasi-bi-Hamiltonian.

###  2.3 Recursion operators and some comments

The bi-Hamiltonian structure $(\omega_{0},\Omega)$ determines a complex
recursion operator $R$ defined as

$\displaystyle\Omega(X,Y)=\omega_{0}(RX,Y),\qquad\forall\,X,Y\in\mathfrak{X}(T^{*}Q).$

But as $\Omega$ and $R$ are complex, we can introduce two real recursion
operator $R_{1}$ and $R_{2}$ defined as

$\displaystyle\Omega_{1}(X,Y)=\omega_{0}(R_{1}X,Y),\qquad\Omega_{2}(X,Y)=\omega_{0}(R_{2}X,Y).$

We recall that $\widehat{\omega_{0}}$ is the map
$\widehat{\omega_{0}}\colon\mathfrak{X}(T^{*}Q)\to\wedge^{1}(T^{*}Q)$ given by
contraction, that is $\widehat{\omega_{0}}(X)=i(X)\omega_{0}$, and then the
nondegenerate character of $\omega_{0}$ means that the map
$\widehat{\omega_{0}}$ is a bijection. Using this notation we can write the
two operators $R_{1}$ and $R_{2}$ as follows

$\displaystyle R_{1}=\widehat{\omega_{0}}^{-1}\circ\widehat{\Omega_{1}},\qquad
R_{2}=\widehat{\omega_{0}}^{-1}\circ\widehat{\Omega_{2}}.$

Then we have the following properties

  * (i)

The coordinates expressions of $R_{1}$ and $R_{2}$ are

$\displaystyle R_{1}=-\alpha_{12}\frac{\partial}{\partial p_{\phi}}\otimes
dr+\left[\alpha_{23}\frac{\partial}{\partial
r}+\alpha_{24}\frac{\partial}{\partial\phi}+\alpha_{12}\frac{\partial}{\partial
p_{r}}\right]\otimes d\phi$
$\displaystyle\hphantom{R_{1}=}{}+\alpha_{23}\frac{\partial}{\partial
p_{\phi}}\otimes dp_{r}+\alpha_{24}\frac{\partial}{\partial p_{\phi}}\otimes
dp_{\phi}$

and

$\displaystyle R_{2}=-\beta_{12}\frac{\partial}{\partial p_{\phi}}\otimes
dr+\left[\beta_{23}\frac{\partial}{\partial
r}+\beta_{24}\frac{\partial}{\partial\phi}+\beta_{12}\frac{\partial}{\partial
p_{r}}\right]\otimes d\phi$
$\displaystyle\hphantom{R_{2}=}{}+\beta_{23}\frac{\partial}{\partial
p_{\phi}}\otimes dp_{r}+\beta_{24}\frac{\partial}{\partial p_{\phi}}\otimes
dp_{\phi}.$

  * (ii)

$R_{1}$ and $R_{2}$ have two different eigenvalues doubly degenerate and one
of them is null (that is, $\lambda_{1}=\lambda_{2}=0$,
$\lambda_{3}=\lambda_{4}\neq 0$). Therefore we have

$\displaystyle\det[R_{1}]=\det[R_{2}]=0,$

what is a consequence of the singular character of $\Omega_{1}$ and
$\Omega_{2}$.

We close this section summarizing the situation we have arrived. We have first
introduced two complex functions, $M_{r}$ and $N_{\phi}$, mainly because of
the behaviour of their Poisson brackets. Then we have proved that they are
interesting for two reasons: first because they determine the existence of
superintegrability (existence of additional constants of motion) and second
because they determine quasi-bi-Hamiltonian structures (first complex
$(\omega_{0},\Omega)$ and then real $(\omega_{0},\Omega_{1},\Omega_{2})$).

Concerning the first point, in this case the additional constants of motion
are just the components of the Runge–Lenz vector (that have been highly
studied making use of different approaches). Now we have arrived to a new
property: they can also be obtained as a consequence of this complex
formalism.

Concerning the second point, the two complex functions $M_{r}$ and $N_{\phi}$
determine the above mentioned geometric structures (first complex and then
real) but unfortunately they are degenerated (we recall that
$\Omega_{1}\wedge\Omega_{1}=0$, $\Omega_{2}\wedge\Omega_{2}=0$). This can be
considered as a limitation of these geometric structures. If a bi-Hamiltonians
structure satisfies all the appropriate properties (that is, symplectic forms,
vanishing of the Nijenhuis torsion of the recursion operator $R$,
diagonalizable recursion operator $R$ with functionally independent real
eigenvalues) then it determines the Liouville integrability of the system. In
fact, the aim of the approach presented in this paper is not to prove the
integrability of a system as consequence of a bi-Hamiltonian structure as we
start with a system known to be not only integrable but also superintegrable.
The existence of $(\Omega_{1},\Omega_{2})$ must be considered, not as a method
for arriving to the integrability of the system, but as a new and interesting
property of the Kepler problem (for the moment only of the two-dimensional
system, the generalization to the three-dimensional case must be considered as
an open question).

##  3 New complex functions and new quasi-bi-Hamiltonian
structures

In this section we will study the existence of new bi-Hamiltonian structures
for the Kepler dynamics by making use of parabolic coordinates $(a,b)$ defined
as

$\displaystyle x=\tfrac{1}{2}\big{(}a^{2}-b^{2}\big{)},\qquad y=ab.$

Of course all previous results can be translated to this new language in such
a way that the functions $M_{r1}$ and $M_{r2}$ are now given by

$\displaystyle
M_{r1}=\frac{(ap_{b}-bp_{a})(ap_{a}+bp_{b})}{\sqrt{a^{2}+b^{2}}},\qquad
M_{r2}=\frac{(ap_{b}-bp_{a})^{2}}{\sqrt{a^{2}+b^{2}}}-g,$

while functions $N_{\phi 1}$ and $N_{\phi 2}$ become

$\displaystyle N_{\phi 1}=\frac{a^{2}-b^{2}}{\sqrt{a^{2}+b^{2}}},\qquad
N_{\phi 2}=\frac{2ab}{\sqrt{a^{2}+b^{2}}}.$

But the important point is that the behaviour of the Kepler Hamiltonian in
these coordinates will permit us to obtain new results different from the
previous ones.

###  3.1 Complex functions and superintegrability

The general form of a natural Euclidean Hamiltonian is

$\displaystyle
H=\frac{1}{2m}\left(\frac{p_{a}^{2}+p_{b}^{2}}{a^{2}+b^{2}}\right)+V(a,b).$

in such a way that if the potential $V$ is of the form

$\displaystyle V(a,b)=\frac{A(a)+B(b)}{a^{2}+b^{2}},$

then the Hamiltonian is Hamilton–Jacobi separable and it is, therefore,
Liouville integrable with the following quadratic function

$\displaystyle
J_{2}=\frac{1}{(a^{2}+b^{2})}(ap_{b}-bp_{a})(ap_{b}+bp_{a})+2\left(\frac{a^{2}B-b^{2}A}{a^{2}+b^{2}}\right)$

as the second constant of motion (the first one is the Hamiltonian itself).

For simplifying the following expressions we introduce the following notation:

$\displaystyle J=ap_{b}-bp_{a},\qquad
P_{x}=\frac{ap_{a}-bp_{b}}{a^{2}+b^{2}},\qquad
P_{y}=\frac{ap_{b}+bp_{a}}{a^{2}+b^{2}}.$

The Hamiltonian of the Kepler problem when written in parabolic coordinates is

$\displaystyle
H_{K}=\frac{1}{2}\left(\frac{p_{a}^{2}+p_{b}^{2}}{a^{2}+b^{2}}\right)+V_{K},\qquad
V_{K}=-\frac{g}{a^{2}+b^{2}},$ (1)

and the Kepler dynamics is given by the following vector field

$\displaystyle\Gamma_{K}=\left(\frac{p_{a}}{a^{2}+b^{2}}\right)\frac{\partial}{\partial
a}+\left(\frac{p_{b}}{a^{2}+b^{2}}\right)\frac{\partial}{\partial
b}+\left(\frac{p_{a}^{2}+p_{b}^{2}-2g}{(a^{2}+b^{2})^{2}}\right)a\frac{\partial}{\partial
p_{a}}+\left(\frac{p_{a}^{2}+p_{b}^{2}-2g}{(a^{2}+b^{2})^{2}}\right)b\frac{\partial}{\partial
p_{b}},$

in such a way that, as $\omega_{0}=da\wedge dp_{a}+db\wedge dp_{b}$, we have

$\displaystyle i(\Gamma_{K})(da\wedge dp_{a}+db\wedge dp_{b})=dH_{K}.$

The Hamiltonian $H_{K}$ is Hamilton–Jacobi separable in coordinates $(a,b)$
and the associated quadratic constant of motion is the component $R_{x}$ of
the Laplace–Runge–Lenz vector

$\displaystyle R_{x}=JP_{y}-g\left(\frac{a^{2}-b^{2}}{a^{2}+b^{2}}\right).$

Let us now introduce the functions $M_{aj}$ and $M_{bj}$, $j=1,2$, defined by

$\displaystyle M_{a1}=\frac{Jp_{a}}{\sqrt{a^{2}+b^{2}}},\qquad
M_{a2}=\frac{2ga-Jp_{b}}{\sqrt{a^{2}+b^{2}}},$

and

$\displaystyle M_{b1}=\frac{Jp_{b}}{\sqrt{a^{2}+b^{2}}},\qquad
M_{b2}=\frac{2gb+Jp_{a}}{\sqrt{a^{2}+b^{2}}}.$

Then, the important property is that the Poisson bracket of the function
$M_{a1}$ with $H_{K}$ is proportional to $M_{a2}$ while the Poisson bracket of
$M_{a2}$ with $H_{K}$ is proportional to $M_{a1}$, but with the opposite sign:

$\displaystyle\\{M_{a1},H_{K}\\}=-\lambda
M_{a2},\qquad\\{M_{a2},H_{K}\\}=\lambda M_{a1},$

and the same is true for the functions $M_{b1}$ and $M_{b2}$,

$\displaystyle\\{M_{b1},H_{K}\\}=-\lambda
M_{b2},\qquad\\{M_{b2},H_{K}\\}=\lambda M_{b1},$

where now $\lambda$ denotes the following function

$\displaystyle\lambda=\frac{ap_{b}-bp_{a}}{(a^{2}+b^{2})^{2}}.$

Therefore the two complex functions $M_{a}$ and $M_{b}$ defined as

$\displaystyle M_{a}=M_{a1}+iM_{a2},\qquad M_{b}=M_{b1}+iM_{b2},$

satisfy

$\displaystyle\\{M_{a},H_{K}\\}=i\lambda
M_{a},\qquad\\{M_{b},H_{K}\\}=i\lambda M_{b}.$

######  Proposition 3.1.

The complex function $K_{34}$ defined as

$\displaystyle K_{34}=M_{a}M_{b}^{*}$

is a $($complex$)$ constant of the motion for the dynamics of the Kepler
problem described by the Hamiltonian (1).

We omit the proof because it is quite similar to the proof of the previous
Proposition 2.1.

Note that the modulus of the complex functions $M_{a}$ and $M_{b}$ are given
by

$\displaystyle M_{a}M_{a}^{*}=2\big{(}J^{2}H_{K}-gR_{x}+g^{2}\big{)},\qquad
M_{b}M_{b}^{*}=2\big{(}J^{2}H_{K}+gR_{x}+g^{2}\big{)},$

and then

$\displaystyle M_{a}M_{a}^{*}+M_{b}M_{b}^{*}=4\big{(}J^{2}H_{K}+g^{2}\big{)}.$

Of course the complex function $K_{34}$ determines two real functions that are
first integrals for the Kepler problem

$\displaystyle
K_{34}=K_{3}+iK_{4},\qquad\\{K_{3},H_{K}\\}=0,\qquad\\{K_{4},H_{K}\\}=0,$

with $K_{3}$ and $K_{4}$ given by

$\displaystyle K_{3}={\rm
Re}(K_{34})=M_{a1}M_{b1}+M_{a2}M_{b2}=JP_{x}+g\left(\frac{2ab}{a^{2}+b^{2}}\right),$
$\displaystyle K_{4}={\rm Im}(K_{34})=M_{a2}M_{b1}-M_{a1}M_{b2}=-2J^{2}H_{K}.$

Remark that the function $K_{3}$ is the other Laplace–Runge–Lenz constant,
while $K_{4}$, that is a fourth order polynomial in the momenta, determines
the angular momentum $J$ as a factor.

###  3.2 Complex functions and quasi-bi-Hamiltonian structures

First we recall that the complex functions $M_{a}$ y $M_{b}$ are given by

$\displaystyle
M_{a}=\left(\frac{Jp_{a}}{\sqrt{a^{2}+b^{2}}}\right)+i\left(\frac{2ga-
Jp_{b}}{\sqrt{a^{2}+b^{2}}}\right),\qquad
M_{b}=\left(\frac{Jp_{b}}{\sqrt{a^{2}+b^{2}}}\right)+i\left(\frac{2gb+Jp_{a}}{\sqrt{a^{2}+b^{2}}}\right).$

Let us now denote by $Z_{34}$ the Hamiltonian vector field of the function
$K_{34}$, i.e., $i(Z_{34})\omega_{0}=dZ_{34}$, such that $Z_{34}(H_{K})=0$,
and by $Z_{a}$ and $Z_{b}$ the Hamiltonian vector fields of the complex
functions $M_{a}$ and $M_{b}$, that is,

$\displaystyle i(Z_{a})\omega_{0}=dM_{a},\qquad i(Z_{b})\omega_{0}=dM_{b}.$

Their coordinate expressions are given by

$\displaystyle Z_{a}=\left(\frac{\partial M_{a}}{\partial
p_{a}}\right)\frac{\partial}{\partial a}+\left(\frac{\partial M_{a}}{\partial
p_{b}}\right)\frac{\partial}{\partial b}-\left(\frac{\partial M_{a}}{\partial
a}\right)\frac{\partial}{\partial p_{a}}-\left(\frac{\partial M_{a}}{\partial
b}\right)\frac{\partial}{\partial p_{b}}=Z_{a1}+iZ_{a2},$

with $Z_{a1}$ and $Z_{a2}$ given by

$\displaystyle
Z_{a1}=\frac{1}{\sqrt{a^{2}+b^{2}}}\left((ap_{b}-2bp_{a})\frac{\partial}{\partial
a}+ap_{a}\frac{\partial}{\partial
b}-\frac{(ap_{a}+bp_{b})}{a^{2}+b^{2}}\left(bp_{a}\frac{\partial}{\partial
p_{a}}-ap_{a}\frac{\partial}{\partial p_{b}}\right)\right),$ $\displaystyle
Z_{a2}=\frac{1}{\sqrt{a^{2}+b^{2}}}\left(bp_{b}\frac{\partial}{\partial
a}+(bp_{a}-2ap_{b})\frac{\partial}{\partial
b}-\frac{(ap_{a}+bp_{b})p_{b}-2gb}{a^{2}+b^{2}}\left(-b\frac{\partial}{\partial
p_{a}}+a\frac{\partial}{\partial p_{b}}\right)\right),$

and

$\displaystyle Z_{b}=\left(\frac{\partial M_{b}}{\partial
p_{a}}\right)\frac{\partial}{\partial a}+\left(\frac{\partial M_{b}}{\partial
p_{b}}\right)\frac{\partial}{\partial b}-\left(\frac{\partial M_{b}}{\partial
a}\right)\frac{\partial}{\partial p_{a}}-\left(\frac{\partial M_{b}}{\partial
b}\right)\frac{\partial}{\partial p_{b}}=Z_{b1}+iZ_{b2},$

with $Z_{b1}$ and $Z_{b2}$ given by

$\displaystyle
Z_{b1}=\frac{1}{\sqrt{a^{2}+b^{2}}}\left(-bp_{b}\frac{\partial}{\partial
a}+(2ap_{b}-bp_{a})\frac{\partial}{\partial
b}-\frac{(ap_{a}+bp_{b})}{a^{2}+b^{2}}\left(bp_{b}\frac{\partial}{\partial
p_{a}}-ap_{b}\frac{\partial}{\partial p_{b}}\right)\right),$ $\displaystyle
Z_{b2}=\frac{1}{\sqrt{a^{2}+b^{2}}}\left((ap_{b}-2bp_{a})\frac{\partial}{\partial
a}+ap_{a}\frac{\partial}{\partial
b}-\frac{(ap_{a}+bp_{b})p_{a}-2ga}{a^{2}+b^{2}}\left(b\frac{\partial}{\partial
p_{a}}-a\frac{\partial}{\partial p_{b}}\right)\right).$

Now recalling that

$\displaystyle dZ_{34}=d(M_{a}M_{b}^{*})=M_{b}^{*}d(M_{a})+M_{a}d(M_{b}^{*}),$

we obtain

$\displaystyle
Z_{34}=M_{b}^{*}Z_{a}+M_{a}Z_{b}^{*}=Z+Z^{\prime},\qquad\textrm{where}\qquad
Z=M_{b}^{*}Z_{a},\qquad Z^{\prime}=M_{a}Z_{b}^{*}.$

The following proposition is similar to that of Proposition 2.2 and we omit
the proof.

######  Proposition 3.2.

The Lie bracket of the Kepler dynamical vector field $\Gamma_{K}$ with the
vector field $Z$ is given by

$\displaystyle[\Gamma_{K},Z]=iK_{34}X_{\lambda},$

where $X_{\lambda}$ is the Hamiltonian vector field of $\lambda$ solution of
the equation $i(X_{\lambda})\omega_{0}=d\lambda$.

In the following we will denote by $\Omega_{ab}$ the complex 2-form defined as

$\displaystyle\Omega_{ab}=dM_{a}\wedge dM_{b}^{*}$
$\displaystyle\hphantom{\Omega_{ab}}{}=d\left[\left(\frac{Jp_{a}}{\sqrt{a^{2}+b^{2}}}\right)+i\left(\frac{2ga-
Jp_{b}}{\sqrt{a^{2}+b^{2}}}\right)\right]\wedge
d\left[\left(\frac{Jp_{b}}{\sqrt{a^{2}+b^{2}}}\right)-i\left(\frac{2gb+Jp_{a}}{\sqrt{a^{2}+b^{2}}}\right)\right].$

Then the two 2-forms $\omega_{Z}$ and $\omega_{Z}^{\prime}$ obtained by Lie
derivation of $\omega_{0}$ with respect to $Z$ and $Z^{\prime}$ are given by

$\displaystyle{\mathcal{L}}_{Z}\omega_{0}=\omega_{Z}=-\Omega_{ab},\qquad{\mathcal{L}}_{Z^{\prime}}\omega_{0}=\omega_{Z}^{\prime}=\Omega_{ab}.$

######  Proposition 3.3.

The Hamiltonian vector field $\Gamma_{K}$ of the Kepler problem is quasi-
Hamiltonian with respect to the complex $2$-form $\Omega_{ab}$.

######  Proof 3.4.

This can be proved by a direct computation

$\displaystyle
i(\Gamma_{K})\Omega_{ab}=\Gamma_{K}(M_{a})dM_{b}^{*}-\Gamma_{K}(M_{b}^{*})dM_{a}=(i\lambda
M_{a})dM_{b}^{*}+(i\lambda M_{b}^{*})dM_{a}=i\lambda
d(M_{a}M_{b}^{*}).\\!\\!\\!\\!\\!\\!$ ∎

The complex 2-form $\Omega_{ab}$ can be decomposed as

$\displaystyle\Omega_{ab}=\Omega_{ab1}+i\Omega_{ab2},$

where the two real 2-forms, $\Omega_{ab1}={\rm Re}(\Omega_{ab})$ and
$\Omega_{ab2}={\rm Im}(\Omega_{ab})$, take the form

$\displaystyle\Omega_{ab1}=dM_{a1}\wedge dM_{b1}+dM_{a2}\wedge
dM_{b2},\qquad\Omega_{ab2}=-dM_{a1}\wedge dM_{b2}+dM_{a2}\wedge dM_{b1},$

Considering the real and imaginary parts we obtain

$\displaystyle i(\Gamma_{K})\Omega_{ab1}=-\lambda dK_{4},\qquad
i(\Gamma_{K})\Omega_{ab2}=\lambda dK_{3},$

what means that $\Gamma_{K}$ is also quasi-bi-Hamiltonian with respect to the
two real 2-forms $(\omega_{0},\Omega_{ab1})$ and $(\omega_{0},\Omega_{ab2})$.

Once more we obtain that the factor $\lambda$ determines that the system is
quasi-bi-Hamiltonian instead of just bi-Hamiltonian.

The complex 2-form $\Omega_{ab}$ is well defined but it is not symplectic. The
kernel is two-dimensional and it is invariant under the action of $\Gamma_{K}$

$\displaystyle[\operatorname{Ker}\Omega_{ab},\Gamma_{K}]\subset\operatorname{Ker}\Omega_{ab}.$

The coordinate expressions of $\Omega_{ab1}$ and $\Omega_{ab2}$ are

$\displaystyle\Omega_{ab1}=\frac{2J}{(a^{2}\\!+b^{2})^{2}}(\alpha_{13}da\wedge
dp_{a}+\alpha_{14}da\wedge dp_{b}+\alpha_{23}db\wedge
dp_{a}+\alpha_{24}db\wedge dp_{b}+\alpha_{34}dp_{a}\wedge dp_{b}),$
$\displaystyle\Omega_{ab2}=\frac{2g}{(a^{2}\\!+b^{2})^{2}}(\beta_{13}da\wedge
dp_{a}+\beta_{14}da\wedge dp_{b}+\beta_{23}db\wedge dp_{a}+\beta_{24}db\wedge
dp_{b}),$

with $\alpha_{ij}$ and $\beta_{ij}$ being given by

$\displaystyle\alpha_{13}=\bigl{(}2gb-ap_{a}p_{b}-bp_{b}^{2}\bigr{)}b,\qquad$
$\displaystyle\alpha_{14}=-\bigl{(}2ga-ap_{a}^{2}-bp_{a}p_{b}\bigr{)}b,$
$\displaystyle\alpha_{23}=\bigl{(}-2gb+ap_{a}p_{b}+bp_{b}^{2}\bigr{)}a,\qquad$
$\displaystyle\alpha_{24}=\bigl{(}2ga-ap_{a}^{2}-bp_{a}p_{b}\bigr{)}a,$
$\displaystyle\alpha_{34}=2J\big{(}a^{2}+b^{2}\big{)},$

and

$\displaystyle\beta_{13}=\bigl{(}2abp_{a}-a^{2}p_{b}-b^{2}p_{b}\bigr{)}b,\qquad$
$\displaystyle\beta_{14}=\bigl{(}2abp_{b}-a^{2}p_{a}-b^{2}p_{a}\bigr{)}b,$
$\displaystyle\beta_{23}=\bigl{(}-2abp_{a}+a^{2}p_{b}+b^{2}p_{b}\bigr{)}a,\qquad$
$\displaystyle\beta_{24}=\bigl{(}-2abp_{b}-a^{2}p_{a}-b^{2}p_{a}\bigr{)}a.$

We close this section with the following properties:

  * (i)

The two real 2-forms are not symplectic. In fact we have verified that
$\Omega_{1}\wedge\Omega_{1}=0$, $\Omega_{2}\wedge\Omega_{2}=0$, and also
$\Omega_{1}\wedge\Omega_{2}=0$.

  * (ii)

These two 2-forms, $\Omega_{ab1}$ and $\Omega_{ab2}$, determine two recursion
operators ($(1,1)$ tensor fields) $R_{ab1}$ and $R_{ab2}$ defined as

$\displaystyle\Omega_{ab1}(X,Y)=\omega_{0}(R_{ab1}X,Y),\qquad\Omega_{ab2}(X,Y)=\omega_{0}(R_{ab2}X,Y),$

or in an equivalent way

$\displaystyle
R_{ab1}=\widehat{\omega_{0}}^{-1}\circ\widehat{\Omega_{ab1}},\qquad
R_{ab2}=\widehat{\omega_{0}}^{-1}\circ\widehat{\Omega_{ab2}}.$

As in Section 2.3, a consequence of the singular character of $\Omega_{ab1}$
and $\Omega_{ab2}$ is that

$\displaystyle\det[R_{ab1}]=\det[R_{ab2}]=0.$

  * (iii)

If we denote by $Z_{3}$ and $Z_{4}$ the Hamiltonian vector fields (with
respect to the canonical symplectic form $\omega_{0}$) of the integrals
$K_{3}$ and $K_{4}$, then the dynamical vector field $\Gamma_{K}$ is
orthogonal to $Z_{4}$ with respect to the structure $\Omega_{1}$ and it is
also orthogonal to $Z_{3}$ with respect to the structure $\Omega_{2}$, that
is,

$\displaystyle i(\Gamma_{K})i(Y_{4})\Omega_{ab1}=0,\qquad
i(\Gamma_{K})i(Y_{3})\Omega_{ab2}=0.$

##  4 Final comments

The Kepler problem is separable in two different coordinate systems, and
because of this, it is superintegrable with quadratic integrals of motion. Now
we have proved that this superintegrability is directly related with the
existence of certain complex functions possessing very interesting Poisson
bracket properties and also that these functions are also related with the
existence of quasi-bi-Hamiltonian structures.

We finalize with some questions for future work. First, as stated in the
Introduction, quasi-bi-Hamiltonian structures is a matter that still remain as
slightly studied (in contrast to the bi-Hamiltonian systems); so the
particular case of the Kepler problem can be a good motivation to undertake a
better study of these systems. Second, the existence of these complex
functions is not a specific characteristic of the Kepler problem; in fact, it
has been proved that the superintegrability of certain systems recently
studied (as the isotonic oscillator or the TTW or PW systems) [11, 29, 30, 31,
32, 33] is also related with such a class of complex functions. Therefore, an
interesting open question is whether these other superintegrable systems are
also endowed with bi-Hamiltonian structures or with quasi-bi-Hamiltonian
structures; the starting point for this study must be a deeper analysis of the
properties of these complex functions making use of the geometric (symplectic)
formalism as an approach.

### Acknowledgments

This work has been supported by the research projects MTM–2012–33575 (MICINN,
Madrid) and DGA-E24/1 (DGA, Zaragoza).

## References

  * [1] Błaszak M., Bi-Hamiltonian representation of Stäckel systems, Phys. Rev. E 79 (2009), 056607, 9 pages, arXiv:0904.2070.
  * [2] Boualem H., Brouzet R., Rakotondralambo J., Quasi-bi-Hamiltonian systems: why the Pfaffian case?, Phys. Lett. A 359 (2006), 559–563.
  * [3] Boualem H., Brouzet R., Rakotondralambo J., About the separability of completely integrable quasi-bi-Hamiltonian systems with compact levels, Differential Geom. Appl. 26 (2008), 583–591.
  * [4] Brouzet R., Caboz R., Rabenivo J., Ravoson V., Two degrees of freedom quasi-bi-Hamiltonian systems, J. Phys. A: Math. Gen. 29 (1996), 2069–2076.
  * [5] Cariñena J.F., Guha P., Rañada M.F., Quasi-Hamiltonian structure and Hojman construction, J. Math. Anal. Appl. 332 (2007), 975–988.
  * [6] Cariñena J.F., Guha P., Rañada M.F., Hamiltonian and quasi-Hamiltonian systems, Nambu–Poisson structures and symmetries, J. Phys. A: Math. Theor. 41 (2008), 335209, 11 pages.
  * [7] Cariñena J.F., Ibort L.A., Non-Noether constants of motion, J. Phys. A: Math. Gen. 16 (1983), 1–7.
  * [8] Cariñena J.F., Marmo G., Rañada M.F., Non-symplectic symmetries and bi-Hamiltonian structures of the rational harmonic oscillator, J. Phys. A: Math. Gen. 35 (2002), L679–L686, hep-th/0210260.
  * [9] Cariñena J.F., Rañada M.F., Canonoid transformations from a geometric perspective, J. Math. Phys. 29 (1988), 2181–2186.
  * [10] Cariñena J.F., Rañada M.F., Generating functions, bi-Hamiltonian systems, and the quadratic-Hamiltonian theorem, J. Math. Phys. 31 (1990), 801–807.
  * [11] Cariñena J.F., Rañada M.F., Santander M., Two important examples of nonlinear oscillators, in Proceedings of 10th International Conference in Modern Group Analysis (Larnaca, Cyprus, 2004), Editors N.H. Ibragimov, C. Sophocleous, P.A. Damianou, 2005, 39–46, math-ph/0505028.
  * [12] Crampin M., Marmo G., Rubano C., Conditions for the complete integrability of a dynamical system admitting alternative Lagrangians, Phys. Lett. A 97 (1983), 88–90.
  * [13] Crampin M., Prince G.E., Alternative Lagrangians for spherically symmetric potentials, J. Math. Phys. 29 (1988), 1551–1554.
  * [14] Crampin M., Sarlet W., A class of nonconservative Lagrangian systems on Riemannian manifolds, J. Math. Phys. 42 (2001), 4313–4326.
  * [15] Crampin M., Sarlet W., Bi-quasi-Hamiltonian systems, J. Math. Phys. 43 (2002), 2505–2517.
  * [16] Crampin M., Sarlet W., Thompson G., Bi-differential calculi and bi-Hamiltonian systems, J. Phys. A: Math. Gen. 33 (2000), L177–L180.
  * [17] De Filippo S., Vilasi G., Marmo G., Salerno M., A new characterization of completely integrable systems, Nuovo Cimento B 83 (1984), 97–112.
  * [18] Fernandes R.L., Completely integrable bi-Hamiltonian systems, J. Dynam. Differential Equations 6 (1994), 53–69.
  * [19] Grigoryev Y.A., Tsiganov A.V., On bi-Hamiltonian formulation of the perturbed Kepler problem, J. Phys. A: Math. Theor. 48 (2015), 175206, 7 pages, arXiv:1504.02211.
  * [20] Henneaux M., Shepley L.C., Lagrangians for spherically symmetric potentials, J. Math. Phys. 23 (1982), 2101–2107.
  * [21] Marmo G., A geometrical characterization of completely integrable systems, in Proceedings of the International Meeting on Geometry and Physics (Florence, 1982), Pitagora, Bologna, 1983, 257–262.
  * [22] McLenaghan R.G., Smirnov R.G., A class of Liouville-integrable Hamiltonian systems with two degrees of freedom, J. Math. Phys. 41 (2000), 6879–6889.
  * [23] Morosi C., Tondo G., Quasi-bi-Hamiltonian systems and separability, J. Phys. A: Math. Gen. 30 (1997), 2799–2806, solv-int/9702006.
  * [24] Morosi C., Tondo G., On a class of dynamical systems both quasi-bi-Hamiltonian and bi-Hamiltonian, Phys. Lett. A 247 (1998), 59–64, solv-int/9811007.
  * [25] Rakotondralambo J., A cohomological obstruction for global quasi-bi-Hamiltonian fields, Phys. Lett. A 375 (2011), 1064–1068.
  * [26] Rañada M.F., Alternative Lagrangians for central potentials, J. Math. Phys. 32 (1991), 2764–2769.
  * [27] Rañada M.F., Dynamical symmetries, bi-Hamiltonian structures, and superintegrable $n=2$ systems, J. Math. Phys. 41 (2000), 2121–2134.
  * [28] Rañada M.F., A system of $n=3$ coupled oscillators with magnetic terms: symmetries and integrals of motion, SIGMA 1 (2005), 004, 7 pages, nlin.SI/0511007.
  * [29] Rañada M.F., A new approach to the higher order superintegrability of the Tremblay–Turbiner–Winternitz system, J. Phys. A: Math. Theor. 45 (2012), 465203, 9 pages, arXiv:1211.2919.
  * [30] Rañada M.F., Higher order superintegrability of separable potentials with a new approach to the Post–Winternitz system, J. Phys. A: Math. Theor. 46 (2013), 125206, 9 pages, arXiv:1303.4877.
  * [31] Rañada M.F., The Tremblay–Turbiner–Winternitz system on spherical and hyperbolic spaces: superintegrability, curvature-dependent formalism and complex factorization, J. Phys. A: Math. Theor. 47 (2014), 165203, 9 pages, arXiv:1403.6266.
  * [32] Rañada M.F., The Post–Winternitz system on spherical and hyperbolic spaces: a proof of the superintegrability making use of complex functions and a curvature-dependent formalism, Phys. Lett. A 379 (2015), 2267–2271.
  * [33] Rañada M.F., Rodríguez M.A., Santander M., A new proof of the higher-order superintegrability of a noncentral oscillator with inversely quadratic nonlinearities, J. Math. Phys. 51 (2010), 042901, 11 pages, arXiv:1002.3870.
  * [34] Rauch-Wojciechowski S., A bi-Hamiltonian formulation for separable potentials and its application to the Kepler problem and the Euler problem of two centers of gravitation, Phys. Lett. A 160 (1991), 149–154.
  * [35] Tempesta P., Tondo G., Generalized Lenard chains, separation of variables, and superintegrability, Phys. Rev. E 85 (2012), 046602, 11 pages, arXiv:1205.6937.
  * [36] Zeng Y.B., Ma W.X., Families of quasi-bi-Hamiltonian systems and separability, J. Math. Phys. 40 (1999), 4452–4473, solv-int/9906005.


# Optimal transport maps as flows of control-affine systems ††thanks: This
work is supported by the MIUR Excellence Department Project MatMod@TOV awarded
to the Department of Mathematics, University of Rome Tor Vergata (CUP
E83C23000330006), by the projects PRIN 2022-E53D23005910006 and PRIN 2022
PNRR- E53D23017910001 founded by the EU - Next Generation EU, and by Gnampa -
INdAM.

Marco Caponigro Dipartimento di Matematica, Università di Roma “Tor Vergata”,
Via della Ricerca Scientifica 1, 00133 Roma, Italy<EMAIL_ADDRESS>Arianna Vicari Dipartimento di Matematica, Università di Roma “Tor Vergata”,
Via della Ricerca Scientifica 1, 00133 Roma, Italy<EMAIL_ADDRESS>
###### Abstract

We consider the controllability problem for the continuity equation where the
dynamics is governed by a driftless control-affine system. Under suitable
regularity conditions, the controllability of the system is a sufficient
condition for the controllability of the continuity equation by means of time-
varying feedback controls. Moreover, we show that there exist controls such
that the flow of the control system is the optimal transport map, for the
$2$-Wasserstein distance, between two given probability measures.

##  1 INTRODUCTION

The analysis of the continuity equation represents the fundamental step in the
analysis of many physical equations as e.g. the Boltzmann equation, the
Fokker-Planck equation, the Vlasov equation, arising in several areas of
sciences (see, for instance, [C+69, BCGM15, ABBG+12]). The continuity
equation, or Liouville model, boasts many possible applications ranging from
crowd dynamics [CHM11, CPT14] to opinion formation [BHT13], from herd analysis
[TB04] to autonomous multi-vehicle navigation [RB08] and, in general, it can
be used to model uncertainties on the initial state of a deterministic
controlled system [MQ18, CQ08, JJZ22a, JJZ22b]. The controlled continuity
equation can be interpreted as the application to multiple trials involving a
single system or, again, as applying to the control of many identical copies
of a single system; e.g., control of a flock.

The problem is to study the evolution of a large number of agents whose
interactions depend on the distribution of the whole group of agents. Adopting
a _macroscopic_ viewpoint, we consider the collection of all agents in a state
space $X$, to be a mass evolving through time. If the total number of the
agents is constant, then, up to a renormalization, the total mass of the
system is constantly equal to $1$. Hence, the evolution of the system is
described by a curve $t\mapsto\mu_{t}$, where $\mu_{t}$ is the distribution of
probability of the system at time $t$. The conservation of the mass of the
system for any time $t$ is described by the Continuity, or Liouville, equation

$\partial_{t}\mu_{t}+\mathrm{div}(F\cdot\mu_{t})=0,$ (1)

where $F$ is a smooth vector field. To model the action of a controller on the
system we consider a vector field $F$ which is linear-in-control, namely

$F(x)=\sum_{i=1}^{m}u_{i}f_{i}(x),$ (2)

with $m<n$, where $u_{i}$ are the scalar controls, and $f_{1},\ldots,f_{m}$
are smooth vector fields on $X$.

The role of the controller is to steer the solution of (1) from a given
initial measure $\mu$ to a final one $\nu$ in finite time.

This controllability problem can be seen as a continuous analogous of
“ensemble controllability”: i.e. a problem involving a large number of more or
less identical subsystems which are being manipulated by a single source of
command signals [BK00, ABS16].

In order to study these kind of problems it is natural to work with time-
varying feedback controls. Indeed, if the $u_{i}$’s are continuous feedback
controls not depending on time even the realization of locally asymptotically
stable equilibria for a control system $\dot{x}=f(x,u)$ is, in general, not
possible [Bro83] For this kind of problem the use of periodic time-varying
feedback controls has been introduced by J.-M. Coron (see Coron [Cor92, Cor95]
or [Cor07, Section 11.2]).

In this paper we prove that, under suitable regularity conditions, the
controllability of

$\dot{x}=\sum_{i=1}^{m}u_{i}f_{i}(x),$ (3)

by means of open-loop controls is a sufficient condition for the
controllability of (1) by means of time-varying feedback controls (Theorem 3
below). Moreover, in Theorem 8, we show that it is possible to write the
optimal transport map, for the $2$-Wasserstein distance (see Definition 6
below), between two sufficiently regular measures as flow at time $1$ of
system (3).

The rest of this section contains a literature review on the subject and the
main notations useful in the statement of the main results, presented in
Section 2. The proofs of the main results are in Section 3

###  1.1 State of the art

The problem of the controllability of the Liouville Equation (1) has been
proposed by R.W. Brockett in [Bro07] as a problem of ensemble controllability
applied to robotics (see also [ABBG+12]). Brockett studied the density
equation associated with linear systems. A result, phrased in terms of
controlling diffeomorphisms, is published in [AC09]. The result is valid for
driftless control-affine systems on compact manifolds under the sole
assumption of controllability of the system. In [AS22], the authors prove
approximate controllability for the density equation using open-loop controls,
under the additional hypothesis of _strong Lie Algebra rank condition_. The
result has then been used in [SF23] to approximate the optimal transport map
between two smooth densities with the flow of system 3.

The first result linking the controllability of the continuity equation and
the optimal transport is due to Khesin and Lee [KL09] (see also [AL09]). They
prove that, under the Lie Algebra rank condition the optimal transport map is
theflow of a driftless control-affine system associated with time-varying
feedback controls. In this paper, Theorem 3, in fact extends this result to
controllable systems without any assumption on the Lie brackets.

Beside controllability, we mention also optimal controllability results for
(1) with open loop controls: In [BBFR19] for a particular class of vector
fields of the form $F=a(t,x)+u_{1}(t)+xu_{2}(t),$ and in [Pog16] for optimal
control without terminal constraints for nonlinear systems.

We also refer to the recent [CGP21] for a review of optimal mass transport
results casted as stochastic control problems.

###  1.2 Notations and definitions

We consider an open subset $X$ of $\mathbb{R}^{n}$ and a control affine system

$\dot{x}(t)=\displaystyle\sum_{i=1}^{m}u_{i}f_{i}(x(t)),$ (4)

where $\\{f_{1},\dots,f_{m}\\}\subset\mathrm{Vec}(X)$ is a family of smooth,
i.e. $C^{\infty}$ vector fields on $X$ and $u_{1},\dots,u_{m}$ are scalar
controls with values in $\mathbb{R}$.

Admissible controls $u_{1},\ldots,u_{m}$ are taken in a class $\mathscr{U}$
that is either the one of open-loop controls or the class of time-varying
feedback controls. More precisely, we denote with the letter $\mathbf{u}$ the
$m$-tuple of open-loop controls $(u_{1}(t),\ldots,u_{m}(t))$ with
$u_{j}:[0,\infty)\to\mathbb{R}$ measurable bounded. We use the letter
$\mathbf{v}=(v_{1}(t,x),\ldots,v_{m}(t,x))$ for controls that are measurable
bounded with respect to $t$ and smooth with respect to $x\in X$.

Moreover we assume that the vector fields $f_{1},\dots,f_{m}$ satisfy the
growth condition $f_{j}(x)\leq\alpha(1+|x|)$, for $\alpha$ constant.

These conditions guarantee global existence and uniqueness of the solution of
the Cauchy problem

$\begin{cases}\dot{x}(t)&=\displaystyle\sum_{i=1}^{m}u_{i}f_{i}(x(t)),\quad\mbox{
a.e. }t\geq 0\\\ x(0)&=x_{0},\end{cases}$ (5)

for any $x_{0}\in X$, in the Carathéodory sense (see e.g. [Hal09, Theorem
5.3]), associated with controls in $\mathscr{U}$.

Throughout the paper, by a diffeomorphism, we mean a $C^{\infty}$ invertible
map with $C^{\infty}$ inverse. We denote by $\mathrm{Diff}_{0}(X)$ the
(infinite dimensional) group of orientation preserving diffeomorphisms, and by
$C^{\infty}(X)$ the space of smooth functions from $X$ to $\mathbb{R}$. Both
$\mathrm{Diff}_{0}(X)$ and $C^{\infty}(X)$ are endowed with the standard
topology of the uniform convergence of the partial derivatives of any order on
any compact subset of $X$, (see e.g. [Rud91]).

######  Definition 1.

We say that (4) is _controllable_ in $X$ if for any pair of points
$x_{0},x_{1}$ in $X$, there exists an open-loop control $\mathbf{u}$, such
that the associated solution $x(t;\mathbf{u})$ of (5) satisfies
$x(1;\mathbf{u})=x_{1}$.

######  Remark 2.

It is not restricting to consider the controllability in time $1$. Indeed the
controls take values in $\mathbb{R}$. Since the system is linear in the
control, up to a time-rescaling, this is equivalent to the controllability in
large time with bounded controls.

Given a Borel probability measure $\mu\in\mathscr{P}(X)$ we consider the
Cauchy problem

$\begin{cases}\partial_{t}\mu_{t}&+\mathrm{div}\left(\sum_{i=1}^{m}u_{i}f_{i}\cdot\mu_{t}\right)=0,\\\
\mu(0)&=\mu.\end{cases}$ (6)

A sufficient condition for the existence and uniqueness of the solutions for
the Cauchy problem (6) is the existence and uniqueness of solutions for (5)
(see [AGS05]). Moreover, if $\Phi^{t}$ is the flow at time $t$ of (5) then the
solution of (6) is

$\mu_{t}=\Phi^{t}_{\\#}\mu$ (7)

that is, for any Borel set $A\subset X$

$\Phi^{t}_{\\#}\mu(A):=\mu(\Phi^{-t}(A)).$

We denote by $\mathscr{L}^{n}$ the Lebesgue measure of $\mathbb{R}^{n}$. If
$\mu$ is absolutely continuous with respect to $\mathscr{L}^{n}$ then so is
the solution $\mu_{t}$ of (6) for every $t\geq 0$. Namely if
$\mu=\rho_{0}\mathscr{L}^{n}$ for some Lipschitz function
$\rho_{0}:X\to\mathbb{R}$ then the solution $\mu_{t}=\rho_{t}\mathscr{L}^{n}$
where $\rho_{t}:X\to\mathbb{R}$ is Lipschitz for every $t$ and satisfies the
_density equation_

$\rho_{t}(x)=\frac{1}{\mathrm{det}\nabla\Phi^{t}(x)}\rho_{0}\left(\Phi^{-t}(x)\right),$
(8)

for almost every $t\geq 0$ and for every $x\in X$.

##  2 MAIN RESULTS

###  2.1 Controllability of the continuity equation

Our main result is a controllability result for the continuity equation, in
the sense that every pair of measures with smooth positive densities can be
joined by an admissible solution of (6) under the sole assumption of
controllability of the underlying dynamics on $\mathbb{R}^{n}$.

######  Theorem 3.

Let $\Omega$ be a smooth, bounded, connected open subset of $\mathbb{R}^{n}$
compactly contained in $X$. If (4) is controllable in $X$ then for
$\mu=\rho_{\mu}\mathscr{L}^{n}$ and $\nu=\rho_{\nu}\mathscr{L}^{n}$, two
probability measures on $\Omega$ with smooth, strictly positive densities
$\rho_{\mu}$ and $\rho_{\nu}$, there exists a time-varying feedback control
such that the associated solution $\mu_{t}$ of (6) satisfies

$\mu_{0}=\mu,\qquad\mbox{ and }\qquad\mu_{1}=\nu.$

The proof of Theorem 3 relies mainly on the following result of
controllability on the group of diffeomorphisms, in the spirit of [AC09]. The
main difference with the results in literature is the lack of compactness of
the set $\Omega$ of definition of the measures. To overcome this difficulty we
must assume that $X$ contains the closure of $\Omega$ and study the problem on
$\overline{\Omega}$. The proof is rather technical and, for the sake of
redability is postponed in Section 3.

######  Proposition 4.

Let $\Omega$ be a smooth, bounded, connected open subset of $\mathbb{R}^{n}$
compactly contained in $X$. If (4) is controllable in $X$ then for any
$P\in\mathrm{Diff}_{0}(\Omega)$ there exists a time-varying feedback control
$\mathbf{v}:(t,x)\mapsto\Big{(}v_{1}(t,x),...,v_{m}(t,x)\Big{)}$ such that

$P=\Phi^{1}_{\mathbf{v}}\Big{|}_{\Omega},$

where $\Phi^{1}_{\mathbf{v}}$ is the flow at time $1$ of (4) associated with
the control $\mathbf{v}$.

###  2.2 Optimal transport via time-varying feedback

We introduce some basic notion of optimal transportation useful to state the
main result of the Section. For further details we refer to the book [V+09].

######  Definition 5.

Given $\mu,\nu\in\mathscr{P}(X)$, a measurable map $T:X\to X$ is called a
_transport map_ from $\mu$ to $\nu$ if $T_{\\#}\mu=\nu$, that is,

$\nu(A)=\mu(T^{-1}(A))$

for any $A\subset X$ Borel. The measure $T_{\\#}\mu$ is called _pushforward
measure_ of $\mu$ through the map $T$.

The pushforward measure satisfies

$\int_{\mathbb{R}^{m}}\varphi(x)\hskip 2.84526ptdT_{\\#}\mu(x)\hskip
2.84526pt\hskip 2.84526pt=\hskip 2.84526pt\hskip
2.84526pt\int_{\mathbb{R}^{n}}\varphi\circ T(x)\hskip 2.84526ptd\mu(x),$

for every test function $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ smooth with
compact support. We denote by $\mathscr{P}_{2}(\mathbb{R}^{n})$ the set of
Borel probability measures having finite second moment, namely

$\mathscr{P}_{2}(\mathbb{R}^{n}):=\Big{\\{}\mu\in\mathscr{P}(\mathbb{R}^{n})\hskip
2.84526pt:\hskip 2.84526pt\int_{\mathbb{R}^{n}}|x|^{2}\hskip
2.84526ptd\mu(x)\hskip 2.84526pt<\hskip 2.84526pt+\infty\Big{\\}}.$

For every pair of probability measures $\mu,\nu\in\mathscr{P}(\mathbb{R}^{n})$
we define the set of _admissible transport plans_ between $\mu$ and $\nu$ as

$\displaystyle\textrm{Adm}(\mu,\nu)$
$\displaystyle:=\\{\gamma\in\mathscr{P}(\mathbb{R}^{n}\times\mathbb{R}^{n}):$
$\displaystyle(P_{1})_{\\#}\gamma=\mu,(P_{2})_{\\#}\gamma=\nu\\},$

where
$P_{1},P_{2}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ are
the canonical projections on the first and the second component, respectively.

######  Definition 6.

For every pair of probability measures
$\mu,\nu\in\mathscr{P}_{2}(\mathbb{R}^{n})$, the $2$-_Wasserstein distance_
$W_{2}(\mu,\nu)$ between $\mu$ and $\nu$ is defined as

$\displaystyle W^{2}_{2}(\mu,\nu)$
$\displaystyle:=\mathrm{inf}\Big{\\{}\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}|x-y|^{2}d\gamma(x,y):$
$\displaystyle\qquad\gamma\in\mathrm{Adm}(\mu,\nu)\Big{\\}}.$ (9)

We denote by $\mathrm{Opt}(\mu,\nu)$ the set of admissible transport plans
realizing the infimum in (9).

$\displaystyle\mathrm{Opt}(\mu,\nu)$
$\displaystyle:=\Big{\\{}\gamma\in\textrm{Adm}(\mu,\nu)\hskip 2.84526pt:$
$\displaystyle\qquad\int_{\mathbb{R}^{n}\times\mathbb{R}^{n}}|x-y|^{2}\hskip
2.84526ptd\gamma(x,y)=W_{2}^{2}(\mu,\nu)\Big{\\}}.$

######  Definition 7.

Given a Borel map $T:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, we say that $T$
is an optimal transport map between $\mu,\nu\hskip
2.84526pt\in\mathscr{P}_{2}(\mathbb{R}^{n})$ if
$\gamma_{T}:=(\mathrm{Id},T)_{\\#}\mu\in\mathrm{Opt}(\mu,\nu)$.

######  Theorem 8.

Let $\Omega$ be a smooth, bounded, connected, uniformly convex, open subset of
$\mathbb{R}^{n}$ compactly contained in $X$. Let (4) be controllable in $X$.
For every $\mu,\nu\in\mathscr{P}_{2}(\Omega)$ such that

  * $(i)$

$\mu=\rho_{\mu}\mathscr{L}_{\mathbb{R}^{n}}$ and
$\nu=\rho_{\nu}\mathscr{L}_{\mathbb{R}^{n}}$;

  * $(ii)$

$\rho_{\mu},\rho_{\nu}\in C^{\infty}(\overline{\Omega},\mathbb{R})$;

  * $(iii)$

$\rho_{\mu},\rho_{\nu}$ are strictly positive on $\overline{\Omega}$;

there exists a time-varying feedback control such that associated flow at time
$1$ is the optimal transport map between $\mu$ and $\nu$.

##  3 PROOFS

###  3.1 Proof of Proposition 4

First, note that $\mathrm{Diff}_{0}(\mathbb{R}^{n})$ coincides with the set of
orientation preserving diffeomorphisms of $\mathbb{R}^{n}$. Indeed a
diffeomorphism isotopic to the identity clearly preserves the orientation.
Conversely, let $P$ be an orientation preserving diffeomorphism of
$\mathbb{R}^{n}$. Without loss of generality, we can suppose that $P$ fixes
the origin possibly taking the time dependent translation $P-tP(0),t\in[0,1]$.
For simplicity rename $P-P(0)$ by $P$ and consider the isotopy
$H(t,x)=P(tx)/t,$ for $t\in(0,1]$ and $H(0,x)=\lim_{t\rightarrow 0}P(tx)/t.$
Since $P$ is orientation preserving then $H(0,\cdot)$ belongs to the connected
component of the identity of the group of linear invertible operators on
$\mathbb{R}^{n}$.

For every vector field $f\in\mathrm{Vec}(X)$, the map which associates with
any $x_{0}\hskip 2.84526pt\in X$ the evaluation at time $t$ of the solution of
$\dot{x}(t)\hskip 2.84526pt=\hskip 2.84526ptf(x(t))$, with initial condition
$x(0)=x_{0}$, is a diffeomorphism from $X$ into itself, denoted by
$e^{tf}:x_{0}\rightarrow e^{tf}(x_{0})$ and called the flow of $f$ at time
$t$. Given the family of vector fields
$\mathcal{F}:=\\{f_{1},\dots,f_{m}\\}\hskip 2.84526pt\subset\mathrm{Vec}(X)$
appearing in (4), we associate to $\mathcal{F}$ the subgroup of
$\mathrm{Diff}_{0}(X)$

$\mathrm{Gr}\mathcal{F}:=\big{\\{}e^{t_{k}f_{k}}\circ\dots\circ
e^{t_{1}f_{1}}:k\in\mathbb{N},\hskip 2.84526ptt_{i}\in\mathbb{R},\hskip
2.84526ptf_{i}\in\mathcal{F}\big{\\}}.$

We also define $\mathcal{P}\subseteq\mathrm{Diff}_{0}(X)$ to be the set

$\displaystyle\mathcal{P}$
$\displaystyle:=\big{\\{}e^{a_{l}f_{j_{l}}}\circ\dots\circ
e^{a_{1}f_{j_{1}}}:$ $\displaystyle l\in\mathbb{N},\hskip
2.84526pta_{1},...,a_{l}\in C^{\infty}(X),\hskip
2.84526ptf_{j_{1}},...,f_{j_{l}}\in\mathcal{F}\big{\\}}.$ (10)

More intuitively, $\mathcal{P}$ is the set of all diffeomorphisms from $X$
into itself that can be written as composition of finitely many flows of
vector fields in $\\{f_{1},\dots,f_{m}\\}$ rescaled by smooth real functions
on $X$.

If $U$ is an open subset of $X$, we define $\mathcal{P}|_{U}:=\\{P|_{U}\hskip
2.84526pt:\hskip 2.84526ptP\in\mathcal{P}\\}$. Morover, for any $x\in X$, we
define $\mathcal{P}_{x}:=\\{P\in\mathcal{P}:P(x)=x\\}$. Finally, for any
$U\subset X$ and $x\in U$, we denote by $C^{\infty}_{x}(U,X)$ the set of all
smooth functions from $U$ to $X$, fixing $x$.
Since system (4) is controllable in the manifold $X$, by [AC09, Lemma 5.1] we
have that any $x\in X$ possesses an open neighborhood $U_{x}\subset X$ such
that the set

$\\{P|_{U_{x}}:P\in\mathcal{P}_{x}\\}$

has nonempty interior in $C^{\infty}_{x}(U_{x},X)$.

This result relies on a geometric idea that goes back to the Orbit Theorem of
Sussmann [Sus73] (see also [AS13]) which allows us to state that for any $x\in
X$ there exist $P_{1},\dots,P_{n}$ diffeomorphisms in
$\mathrm{Gr}(\mathcal{F})$, and $f_{1},\dots,f_{n}$ smooth vector fields in
$\mathcal{F}$ such that, setting $X_{i}(\cdot):=(P_{i})_{*}f_{i}\hskip
2.84526pt\in\mathrm{Vec}(X)$ for $i=1,\dots,n$, we have

$\mathrm{span}\big{\\{}X_{1}(x),\dots,X_{n}(q)\big{\\}}=\mathbb{R}^{n}.$

Moreover, from the smoothness of the vector fields, there exists $V$
neighborhood of $x$ such that
$\mathrm{span}\big{\\{}X_{1}(\xi),\dots,X_{n}(\xi)\big{\\}}=\mathbb{R}^{n}$,
for every $\xi\in\overline{V}$. More precisely $U_{x}$ is the intersection
between $V$ and the set

$\displaystyle\Big{\\{}e^{t_{n}(x)X_{n}}\circ\dots\circ
e^{t_{1}(x)X_{1}}(x)\hskip 2.84526pt:$ $\displaystyle\hskip 2.84526ptt_{i}\in
C^{\infty}(V,\mathbb{R}),\hskip 2.84526pt||t_{i}||_{C^{1}(V)}<\eta_{x}\hskip
2.84526pt\hskip 2.84526pt\forall i=1,\dots,n\Big{\\}},$ (11)

where $\eta_{x}$ has to be chosen small enough so that if
$||t_{i}||_{C^{1}(V)}<\eta_{x}\hskip 2.84526pt\hskip 2.84526pt\forall
i=1,\dots,n$, then
$||\mathrm{D}_{x}\Phi_{h}-\mathrm{Id}||_{C^{0}(\mathbb{R}^{n})}<1$ for any
$h=1,\dots,n$, being $\Phi_{h}$ the smooth diffeomorphism on $V$ defined by
$\Phi_{h}(p):=e^{t_{h}(p)X_{h}}\circ\dots\circ e^{t_{1}(p)X_{1}}(p)$.
Since $\overline{\Omega}$ is compact in $X$, there exist
$x_{1},\dots,x_{k}\hskip 2.84526pt\in\hskip 2.84526pt\overline{\Omega}$ such
that

$\overline{\Omega}\subseteq\bigcup\limits_{i=1}^{k}U_{x_{i}}.$ (12)

We claim that we can assume $x_{1},\dots,x_{k}\hskip 2.84526pt\in\hskip
2.84526pt\Omega$.

Indeed from (11), it follows that for any $i=1,\dots,k$ we can imagine the
neighborhood $U_{x_{i}}$ of $x_{i}$ as the squared ball centered in $x_{i}$,
consisting in all the possible points that we can reach traveling along the
trajectories given by $X_{1},\dots,X_{n}$ for a time shorter than $\eta_{x}$.
Now if $x\in\partial\Omega$, we consider $\overline{x}\in U_{x}\cap\Omega$
close enough to $x$ in a way that will be specified later. The corresponding
neighborhood

$\displaystyle U_{\overline{x}}$
$\displaystyle=\Big{\\{}e^{t_{n}(\overline{x})X_{n}}\circ\dots\circ
e^{t_{1}(\overline{x})X_{1}}(\overline{x})\hskip 2.84526pt:\hskip
2.84526ptt_{i}\in C^{\infty}(V,\mathbb{R}),$
$\displaystyle||t_{i}||_{C^{1}(V)}<\eta_{\overline{x}}\hskip 2.84526pt\hskip
2.84526pt\forall i=1,\dots,n\Big{\\}}\cap V$ (13)

contains a small ball in $\mathbb{R}^{n}$ centered at $\overline{x}$. To be
more precise, we can define

$m_{i}:=\min_{\xi\in\overline{V}}||X_{i}(\xi)||_{2}.$

We observe that since $\mathrm{span}\\{X_{1},\dots,X_{n}\\}=\mathbb{R}^{n}$ in
$\overline{V}$, holds $m_{i}>0$. Now let’s take the minimum over all the
$i=1,\dots,n$

$m:=\min_{i=1,\dots,n}m_{i}\hskip 2.84526pt>\hskip 2.84526pt0.$

Hence the ball $B(\overline{x};\delta)$ is a subset of $U_{\overline{x}}$,
with $\delta=\frac{\eta_{\overline{x}}}{2}m$. Thus, if we take $\overline{x}$
such that the distance between $\overline{x}$ and $x$ is smaller than
$\frac{\delta}{2}$, there exist $N\in\mathbb{N}$ such that

$\displaystyle\frac{U_{x}}{N}$
$\displaystyle:=\Big{\\{}e^{t_{n}(x)X_{n}}\circ\dots\circ
e^{t_{1}(x)X_{1}}(x)\ :\hskip 2.84526ptt_{i}\in C^{\infty}(V,\mathbb{R}),$
$\displaystyle||t_{i}||_{C^{1}(V)}<\frac{\eta_{x}}{N}\forall
i=1,\dots,n\Big{\\}}\cap V$

is contained in $U_{\overline{x}}$.

This allows us to replace in (12) any open set $U_{x_{i}}$ having
$x_{i}\in\partial{\Omega}$ with another open set $U_{\overline{x_{i}}}$ having
$\overline{x_{i}}\in\Omega$. We will omitt the overline sign for the sake of
readability.
Now, setting $V_{x_{i}}:=\Omega\cap U_{x_{i}}$ for any $i\in\\{1,\dots,k\\}$,
we have $\Omega=\bigcup\limits_{i=1}^{k}V_{x_{i}}$.
Notice that $V_{x_{i}}$ is an open subset of $U_{x_{i}}$ and that $x_{i}\in
U_{x_{i}}\cap V_{x_{i}}$ for any $i\in\\{1,\dots,k\\}$.
Let’s fix $i\in\\{1,\dots,k\\}$. Since
$\\{P|_{U_{x_{i}}}:P\in\mathcal{P}_{x_{i}}\\}$ has nonempty interior in
$C^{\infty}_{x_{i}}(U_{x_{i}},X)$, also

$\\{P|_{V_{x_{i}}}:P\in\mathcal{P}_{x_{i}}\\}$ (14)

has nonempty interior in $C^{\infty}_{x_{i}}(V_{x_{i}},X)$.
In addition, the interior of set (14) contains the identical map
$\mathrm{Id}|_{V_{x_{i}}}$.
This follows from the fact that, given an open subset $\mathcal{A}$ of
$C^{\infty}_{x_{i}}(V_{x_{i}},X)$ contained in the interior of (14) and
considered $P_{0}:X\rightarrow X$ such that
$P_{0}|_{V_{x_{i}}}\in\mathcal{A}$, then the set

$P_{0}^{-1}\circ\mathcal{A}=\\{P_{0}^{-1}\circ Q\hskip 2.84526pt:\hskip
2.84526ptQ\in\mathcal{A}\\}$ (15)

is a neighborhood of the identity $\mathrm{Id}|_{V_{x_{i}}}$ contained in
(14). Indeed $\mathrm{Id}|_{V_{x_{i}}}=P_{0}^{-1}\circ P_{0}|_{V_{x_{i}}}$,
hence (15) is a neighborhood of the identity. Moreover for any
$Q\in\mathcal{A}$, $Q\in\mathcal{P}_{x_{i}}|_{V_{x_{i}}}$, thus from (10)
there exist $l\in\mathbb{N},\hskip 2.84526pta_{1},...,a_{l}\in
C^{\infty}(X),\hskip 2.84526ptf_{j_{1}},...,f_{j_{l}}\in\mathcal{F}$ such that
$Q=e^{a_{l}f_{j_{l}}}\circ\dots\circ e^{a_{1}f_{j_{1}}}|_{V_{q_{i}}}$, but
also $P_{0}=e^{b_{w}f_{r_{w}}}\circ\dots\circ e^{b_{1}f_{r_{1}}}$ for some
$w\in\mathbb{N},\hskip 2.84526ptb_{1},...,b_{w}\in C^{\infty}(X),\hskip
2.84526ptf_{r_{1}},...,f_{r_{w}}\in\mathcal{F}$. Thus

$P_{0}^{-1}\circ Q=e^{-b_{1}f_{r_{1}}}\circ\dots\circ e^{-b_{w}f_{r_{w}}}\circ
e^{a_{l}f_{j_{l}}}\circ\dots\circ e^{a_{1}f_{j_{1}}}|_{V_{q_{i}}}$

belongs to (14), so $P_{0}^{-1}\circ\mathcal{A}$ is a subset of (14).

Now let’s define the support of a diffeomorphism $P\in\mathrm{Diff}(\Omega)$
as the set

$\mathrm{supp}P:=\overline{\\{x\in\Omega:P(x)\neq x\\}}.$

According to [AC09, Lemma 5.3.] applied to $\Omega$, it holds that for any
$\mathcal{O}$ neighborhood of the identity $\mathrm{Id}|_{\Omega}$ in
$C^{\infty}(\Omega,X)$, for any $x\in\Omega$ and any neighborhood
$V_{x}\subset\Omega$, we have

$x\in\mathrm{int}\\{P(x):P\in\mathcal{O}\cap\mathcal{P}|_{\Omega},\hskip
2.84526pt\mathrm{supp}P\subset V_{x}\\}.$

Thus, putting together what we obtained so far, we have that there exists
$\mathcal{O}\subset C^{\infty}(\Omega,X)$ neighborhood of the identity
$\mathrm{Id}|_{\Omega}$ such that if $P\in\mathcal{O}$ and
$\mathrm{supp}P\subset V_{x_{i}}$ for some $i\in\\{1,\dots,k\\}$, then
$P\in\mathcal{P}|_{\Omega}$.
At this stage let’s consider $\mathcal{U}$, a neighborhood of the identity in
$\mathrm{Diff}_{0}(\Omega)$. We observe that to any diffeomorphism
$P\in\mathcal{U}$ it is possible to associate a smooth map
$\tilde{P}:\Omega\rightarrow X$ in $C^{\infty}(\Omega,X)$ just replacing the
original codomain with $X$. Furthermore, we have that the norm is preserved,
i.e. $||P||_{\mathrm{Diff}(\Omega)}=||\tilde{P}||_{C^{\infty}(\Omega,X)}$.
Hence we can consider a neighborhood $\mathcal{U}$ of the identity in
$\mathrm{Diff}_{0}(\Omega)$ such that $\tilde{\mathcal{U}}:=\\{\tilde{P}\hskip
2.84526pt:\hskip 2.84526ptP\in\mathcal{U}\\}$ is contained in $\mathcal{O}$.

Now let $P\in\mathcal{U}$ and define the isotopy
$H:\Omega\times[0,1]\rightarrow\Omega$ such that $H(\cdot,0)=\mathrm{Id}$ and
$H(\cdot,1)=P$. Without loss of generality we can assume that
$H(\cdot,t)\in\mathcal{U}$ for any $t\in[0,1]$. Consider a partition of unity

$\\{\lambda_{i}:\Omega\rightarrow\mathbb{R}\hskip
2.84526pt:\mathrm{supp}\lambda_{i}\subset V_{x_{i}}\\}_{i=1}^{k}$

subordinated to the open cover $\\{V_{x_{i}}\\}_{i=1}^{k}$ of $\Omega$.
For any $i\in\\{0,1,\dots,k\\}$, let
$\mu_{i}:\Omega\rightarrow\Omega\times[0,1]$ the map such that for any
$x\in\Omega$

$\mu_{i}(x):=(x,\lambda_{1}(x)+\dots+\lambda_{i}(x)),\hskip 2.84526pt\hskip
2.84526pti\in\\{1,\dots,k\\},$ $\mu_{0}(x):=(x,0).$

Consider $Q_{i}=H\circ\mu_{i}$, then $Q_{i}\in\mathcal{U}$ for any
$i\in\\{0,1,\dots,k\\}$. By construction, $Q_{k}=P$ and $Q_{i}=Q_{i-1}$ on
$\Omega\setminus\mathring{\mathrm{supp}\lambda_{i}}$ for any $i\geq 1$.
Finally, setting $P_{i}:=Q_{i}\circ Q_{i-1}^{-1}$, we have
$P=P_{k}\circ\dots\circ P_{1}$ and $\mathrm{supp}P_{i}\subset V_{x_{i}}$ for
any $i\in\\{1,\dots,k\\}$.
Hence for any $P\in\mathcal{U}$, we proved that there are
$P_{1},\dots,P_{k}\hskip 2.84526pt\in\mathcal{U}$ such that
$P=P_{k}\circ\dots\circ P_{1}$ and $\mathrm{supp}P_{i}\subset V_{x_{i}}$ for
any $i\in\\{1,\dots,k\\}$.
Now let’s consider the corresponding smooth functions
$\tilde{P_{k}},\dots,\tilde{P_{1}}$ in $C^{\infty}(\Omega,X)$. We know by
assumption on $\mathcal{U}$ that
$\tilde{P_{i}}\in\tilde{\mathcal{U}}\subset\mathcal{O}$ for any $i=1,\dots,k$.
Thus $P_{i}\in\mathcal{P}|_{\Omega}$ for any $i=1,\dots,k$ and
$P\in\mathcal{P}|_{\Omega}$, that is to say that the neighborhood of the
identity $\mathcal{U}$ in $\mathrm{Diff}_{0}(\Omega)$ is contained in
$\mathcal{P}|_{\Omega}$.
This concludes the proof, indeed for any $P\in\mathrm{Diff}_{0}(\Omega)$ there
exists a path $\\{P^{t}\hskip 2.84526pt:\hskip 2.84526ptt\in[0,1]\\}\hskip
2.84526pt\subset\hskip 2.84526pt\mathrm{Diff}_{0}(\Omega)$ such that $P^{1}=P$
and $P^{0}=\mathrm{Id}$. For every $N\in\mathbb{N}$, consider the
diffeomorphism $P^{k/N}\circ\big{(}P^{(k-1)/N}\big{)}^{-1}$ for $k=1,...,N$.
For $N\in\mathbb{N}$ sufficiently large
$P^{k/N}\circ\big{(}P^{(k-1)/N}\big{)}^{-1}\hskip 2.84526pt\in\mathcal{U}$ for
every $k=1,...,N$ and, therefore, there exist $a_{1},\dots,a_{\mu N}\in
C^{\infty}(X)$ such that

$\displaystyle P$ $\displaystyle=P\circ\Big{(}P^{(N-1)/N}\Big{)}^{-1}\circ
P^{(N-1)/N}\circ\dots\circ P^{1/N}$ $\displaystyle=e^{a_{\mu N}f_{j_{\mu
N}}}\circ\dots\circ e^{a_{1}f_{j_{1}}}|_{\Omega}$

with $j_{1},\dots,j_{\mu N}\in\\{1,...,m\\}$. Now consider, for $i=1,...,m$
the time-varying feedback controls piecewise constant with respect to $t$
defined on every interval $\Big{[}\frac{k-1}{\mu N},\frac{k}{\mu
N}\Big{)},\hskip 2.84526ptk=1,...,\mu N$ by

$v_{i}(t,x)=\begin{cases}a_{k}(x)\hskip 2.84526pt\hskip 2.84526pt\hskip
2.84526pt\hskip 2.84526pt\hskip 2.84526pt\mathrm{if}\hskip
2.84526ptf_{i}=f_{j_{k}},\\\ 0\hskip 2.84526pt\hskip 2.84526pt\hskip
2.84526pt\hskip 2.84526pt\hskip 2.84526pt\hskip 2.84526pt\hskip
2.84526pt\hskip 2.84526pt\hskip 2.84526pt\hskip 2.84526pt\hskip
2.84526pt\hskip 2.84526pt\mathrm{if}\hskip 2.84526ptf_{i}\neq
f_{j_{k}}.\end{cases}$

Then the flow of (4) at time $1$ corresponding to the controls
$v_{1},\dots,v_{m}$ is $e^{a_{\mu N}f_{j_{\mu N}}}\circ\dots\circ
e^{a_{1}f_{j_{1}}}$.
This completes the proof.

###  3.2 Proof of Theorem 3

We are now ready to prove Theorem 3. Given any pair of probability measures
there exist a transport map between $\mu$ and $\nu$ which is an orientation
preserving diffeomorphism thanks to classical results on the density equation
(8) (see, for instance, [DM90] and [Mos65]). The result then follows from
Proposition 4.

###  3.3 Proof of Theorem 8

Theorem 8 follows from Proposition 4 provided that the optimal transport map
is an orientation preserving diffeomorphism. The result follows from Brenier’s
Theorem [Bre91].

The statement as reported below, for the sake of completeness, together with a
self-contained proof is [SF23, Proposition 3.2].

######  Proposition 9.

Let $\mu=\rho_{\mu}\mathscr{L}_{\mathbb{R}^{n}}$ and
$\nu=\rho_{\nu}\mathscr{L}_{\mathbb{R}^{n}}$ be two probability measures, with
$\rho_{\mu}:\Omega_{1}\rightarrow\mathbb{R}$ and
$\rho_{\nu}:\Omega_{2}\rightarrow\mathbb{R}$, where $\Omega_{1}$ and
$\Omega_{2}$ are open and bounded subsets of $\mathbb{R}^{n}$. Let us assume
that there exists a constant $C>1$ such that $C\geq\rho_{\mu}\geq 1/C$ on
$\Omega_{1}$ and $C\geq\rho_{\nu}\geq 1/C$ on $\Omega_{2}$, and in addition
that

  * •

$\rho_{\mu}\in C^{\infty}(\overline{\Omega_{1}},\mathbb{R})$ and
$\rho_{\nu}\in C^{\infty}(\overline{\Omega_{2}},\mathbb{R})$;

  * •

$\Omega_{1},\Omega_{2}$ are smooth and uniformly convex.

Let $T:\overline{\Omega}_{1}\rightarrow\overline{\Omega}_{2}$ be the optimal
transport map between $\mu$ and $\nu$. Then $T$ is isotopic to the identity.

Let’s call $T:\Omega\rightarrow\Omega$ the optimal transport map between $\mu$
and $\nu$. From classical regularity results in optimal transportation (see
e.g. [AG13, Theorem 1.6], [DPF14, Theorem 3.3], [Caf96] [Caf91]) it follows
that $T$ is a smooth diffeomorphism of $\Omega$. According to Proposition 9
applied to the case $\Omega_{1}=\Omega_{2}$ , $T$ is isotopic to the identity,
hence $T\in\mathrm{Diff}_{0}(\Omega)$. Since (4) is controllable in $X$, we
can apply Proposition 4 and conclude that there exist $m$ time-varying
feedback controls $(t,x)\mapsto\Big{(}v_{1}(t,x),...,v_{m}(t,x)\Big{)}$,
piecewise constant with respect to $t$ and smooth with respect to $x$, such
that the restriction to $\Omega$ of the $1$-time flow of (4) corresponding to
$v_{1},\dots,v_{m}$ is exactly $T$.

## ACKNOWLEDGEMENTS

The authors are grateful to Alfio Borzì and Alessandro Scagliotti for the
precious discussions that inspired this work.

## References

  * [ABBG+12] Fatiha Alabau-Boussouira, Roger Brockett, Olivier Glass, Jérôme Le Rousseau, Enrique Zuazua, and Roger Brockett.  Notes on the control of the liouville equation.  Control of Partial Differential Equations: Cetraro, Italy 2010, Editors: Piermarco Cannarsa, Jean-Michel Coron, pages 101–129, 2012.
  * [ABS16] Andrei Agrachev, Yuliy Baryshnikov, and Andrey Sarychev.  Ensemble controllability by lie algebraic methods.  ESAIM: Control, Optimisation and Calculus of Variations, 22(4):921–938, 2016.
  * [AC09] Andrei A Agrachev and Marco Caponigro.  Controllability on the group of diffeomorphisms.  In Annales de l’Institut Henri Poincaré C, Analyse non linéaire, volume 26, pages 2503–2509. Elsevier, 2009.
  * [AG13] Luigi Ambrosio and Nicola Gigli.  A user’s guide to optimal transport.  Modelling and Optimisation of Flows on Networks Lecture Notes in Mathematics, pages 1–155, 2013.
  * [AGS05] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré.  Gradient flows: in metric spaces and in the space of probability measures.  Springer Science & Business Media, 2005.
  * [AL09] Andrei Agrachev and Paul Lee.  Optimal transportation under nonholonomic constraints.  Transactions of the American Mathematical Society, 361(11):6019–6047, 2009.
  * [AS13] Andrei A Agrachev and Yuri Sachkov.  Control theory from the geometric viewpoint, volume 87.  Springer Science & Business Media, 2013.
  * [AS22] Andrei Agrachev and Andrey Sarychev.  Control on the manifolds of mappings with a view to the deep learning.  Journal of Dynamical and Control Systems, 28(4):989–1008, 2022\.
  * [BBFR19] Jan Bartsch, Alfio Borzì, Francesco Fanelli, and Souvik Roy.  A theoretical investigation of brockett’s ensemble optimal control problems.  Calculus of Variations and Partial Differential Equations, 58:1–34, 2019.
  * [BCGM15] Raul Borsche, Rinaldo M Colombo, Mauro Garavello, and Anne Meurer.  Differential equations modeling crowd interactions.  Journal of Nonlinear Science, 25:827–859, 2015.
  * [BHT13] Nicola Bellomo, Miguel A Herrero, and Andrea Tosin.  On the dynamics of social conflicts: Looking for the black swan.  Kinetic and Related Models, 6(3):459–479, 2013.
  * [BK00] Roger Brockett and Navin Khaneja.  On the stochastic control of quantum ensembles.  In System theory: modeling, analysis and control, pages 75–96. Springer, 2000.
  * [Bre91] Yann Brenier.  Polar factorization and monotone rearrangement of vector-valued functions.  Communications on pure and applied mathematics, 44(4):375–417, 1991\.
  * [Bro83] R. W. Brockett.  Asymptotic stability and feedback stabilization.  In Differential geometric control theory (Houghton, Mich., 1982), volume 27 of Progr. Math., pages 181–191. Birkhäuser Boston, Boston, MA, 1983.
  * [Bro07] Roger W Brockett.  Optimal control of the liouville equation.  AMS IP Studies in Advanced Mathematics, 39:23, 2007.
  * [C+69] Carlo Cercignani et al.  Mathematical methods in kinetic theory, volume 1.  Springer, 1969.
  * [Caf91] Luis A Caffarelli.  Some regularity properties of solutions of monge amp re equation.  Communications on pure and applied mathematics, 44(8-9):965–969, 1991.
  * [Caf96] Luis A Caffarelli.  Boundary regularity of maps with convex potentials–ii.  Annals of mathematics, 144(3):453–496, 1996.
  * [CGP21] Yongxin Chen, Tryphon T Georgiou, and Michele Pavon.  Optimal transport in systems and control.  Annual Review of Control, Robotics, and Autonomous Systems, 4:89–113, 2021.
  * [CHM11] Rinaldo M Colombo, Michael Herty, and Magali Mercier.  Control of the continuity equation with a non local flow.  ESAIM: Control, Optimisation and Calculus of Variations, 17(2):353–379, 2011.
  * [Cor92] Jean-Michel Coron.  Global asymptotic stabilization for controllable systems without drift.  Math. Control Signals Systems, 5(3):295–312, 1992.
  * [Cor95] Jean-Michel Coron.  On the stabilization in finite time of locally controllable systems by means of continuous time-varying feedback law.  SIAM J. Control Optim., 33(3):804–833, 1995.
  * [Cor07] Jean-Michel Coron.  Control and nonlinearity, volume 136 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 2007.
  * [CPT14] Emiliano Cristiani, Benedetto Piccoli, and Andrea Tosin.  Multiscale modeling of pedestrian dynamics, volume 12.  Springer, 2014.
  * [CQ08] Pierre Cardaliaguet and Marc Quincampoix.  Deterministic differential games under probability knowledge of initial condition.  International Game Theory Review, 10(01):1–16, 2008.
  * [DM90] Bernard Dacorogna and Jürgen Moser.  On a partial differential equation involving the jacobian determinant.  7(1):1–26, 1990.
  * [DPF14] Guido De Philippis and Alessio Figalli.  The monge–ampère equation and its link to optimal transportation.  Bulletin of the American Mathematical Society, 51(4):527–580, 2014\.
  * [Hal09] Jack K Hale.  Ordinary differential equations.  Courier Corporation, 2009.
  * [JJZ22a] F. Jean, O. Jerhaoui, and H. Zidani.  A mayer optimal control problem on wasserstein spaces over riemannian manifolds.  IFAC-PapersOnLine, 55(16):44–49, 2022.  18th IFAC Workshop on Control Applications of Optimization CAO 2022.
  * [JJZ22b] Frédéric Jean, Othmane Jerhaoui, and Hasnaa Zidani.  Deterministic optimal control on riemannian manifolds under probability knowledge of the initial condition.  2022\.
  * [KL09] Boris Khesin and Paul Lee.  A nonholonomic moser theorem and optimal transport.  Journal of Symplectic Geometry, 7(4):381–414, 2009.
  * [Mos65] Jürgen Moser.  On the volume elements on a manifold.  Transactions of the American Mathematical Society, 120(2):286–294, 1965.
  * [MQ18] Antonio Marigonda and Marc Quincampoix.  Mayer control problem with probabilistic uncertainty on initial positions.  Journal of Differential Equations, 264(5):3212–3252, 2018.
  * [Pog16] Nikolay Pogodaev.  Optimal control of continuity equations.  Nonlinear Differential Equations and Applications NoDEA, 23:1–24, 2016.
  * [RB08] Wei Ren and Randal W Beard.  Distributed consensus in multi-vehicle cooperative control, volume 27.  Springer, 2008.
  * [Rud91] Walter Rudin.  Functional analysis 2nd ed.  International Series in Pure and Applied Mathematics. McGraw-Hill, Inc., New York, 1991.
  * [SF23] Alessandro Scagliotti and Sara Farinelli.  Normalizing flows as approximations of optimal transport maps via linear-control neural odes.  arXiv preprint arXiv:2311.01404, 2023.
  * [Sus73] Héctor J Sussmann.  Orbits of families of vector fields and integrability of distributions.  Transactions of the American Mathematical Society, 180:171–188, 1973.
  * [TB04] Chad M Topaz and Andrea L Bertozzi.  Swarming patterns in a two-dimensional kinematic model for biological groups.  SIAM Journal on Applied Mathematics, 65(1):152–174, 2004.
  * [V+09] Cédric Villani et al.  Optimal transport: old and new, volume 338.  Springer, 2009.