# Entropic regularised optimal transport in a noncommutative setting

Lorenzo Portinale  Institut für Angewandte Mathematik, Universität Bonn,
Germany

###### Abstract

This survey has been written in occasion of the School and Workshop about
Optimal Transport on Quantum Structures at Erdös Center in September 2022. We
discuss some recent results on noncommutative entropic optimal transport
problems and their relation to the study of the ground-state energy of a
finite-dimensional composite quantum system at positive temperature, following
the work [FGP23]. In the first part, we review some of the classical primal-
dual formulations of optimal transport in the commutative setting, including
extensions to multimarginal problems and entropic regularisation. We discuss
the main features of the entropic problem and show how optimisers can be
efficiently computed via the so-called Sinkhorn algorithm. In the second part,
we discuss how to apply these ideas to a noncommutative setting, in particular
on the space of density matrices over finite dimensional Hilbert spaces. In
this framework, we present equivalences between primal and dual formulations,
and use them to characterise the optimisers. Despite the lack of explicit
formulas due to the noncommutative nature of the problem, one can also show
that a suitable quantum version of the Sinkhorn algorithm converges to the
minimiser of the entropic problem. In the final part of this work, we discuss
similar results for bosonic and fermionic systems.

##  1 Introduction

The goal of this paper is to discuss how the ideas from the theory of
classical optimal transport can be used and transferred to the study of
variational problems in noncommutative analysis. In particular, we focus on
the study of entropy-regularised transport problems (for a classical survey on
the topic, see [Léo14]) and their dual formulations, whose particular
structure provides the existence of a fast algorithm to compute the
corresponding optimisers. This is called Sinhorn’s algorithm and its origin
goes way back to [Sin64]. In the context of Computational Optimal Transport,
the entropic regularization and the Sinkhorn algorithm were introduced and
studied in [Cut13, GS10].

In this work, we are going to review the main features of the Sinkhorn
algorithm in the classical setting, and discuss how to extend its use and
application to the setting of quantum optimal transport, following the recent
work [FGP23]. Using a more physical point of view, a noncommutative version of
the entropy-regularised optimal transport problem naturally arises within the
study of the ground-state energy of a finite-dimensional composite quantum
system at positive temperature. In other words, the aim is to look for the
minimal energy levels of a composite system conditionally to the knowledge of
the states of all its subsystems. Or, in similar terms: which information can
we gather on the state of a composite system when we only have access (e.g.
experimentally) to measurements of the states of its subsystems?

In the final part of the paper, we also discuss the problem when enforcing
symmetry conditions, or more precisely, the fermionic and bosonic framework.
This creates a bridge to what in the physics literature is known as One-body
Reduced Density Matrix Functional Theory ($1$RDMFT), introduced in 1975 by
Gilbert [Gil75] as an extension of the Hohenberg-Kohn (Levy-Lieb) formulation
of Density Functional Theory (DFT) [HK64, Lev76, Lie02], and which can be
framed as a special case of [FGP23].

Some literature: The work [FGP23] is not the first attempt to apply the ideas
of optimal transport to the noncommutative setting. One of the first attempt
was carried out by E. Carlen and J. Maas [CM14], followed by many others. It
is important to observe a fundamental difference between two approaches:
commutative optimal transport can be cast equivalently as a static coupling
problem or as a dynamical optimisation problem. In the noncommutative setting,
the relation (if any) between the two interpretations is not clear. This draws
a line between the two groups of work, the ones that consider the dynamical
formulations (e.g. [CM14, CGT18, MM17, MV21]) and the ones which instead focus
on a static formulation (e.g. [CGP18, DPT21, GP15, MM17, PCVS19]). Concerning
duality formulations in the quantum setting, in [CGP18] the authors study the
case of $\varepsilon=0$ and prove a duality result for the noncommutative
problem in the very same spirit of the Kantorovich duality for the classical
Monge problem. In the dynamical framework, [Wir22] studies the entropic
quantum optimal transport problem as well, in constrast to the static approach
adopted in [FGP23]. For more details on the state of the art, both of quantum
optimal transport and 1RDMFT, see [FGP23] and references therein.

Organisation of the paper: In the first section we review part of the
classical theory of (regularised) optimal transport problems. We then move to
the noncommutative framework, focusing on similarities and differences with
the classical setting. Subsequently, we introduce the noncommutative Sinkhorn
algorithm and discuss the convergence toward optimisers (both in the primal
and dual problem). At last, we present the fermionic/bosonic setting, and
discuss a variational interpretation of the Pauli exclusion principle.

##  2 Classical Optimal Transport and Regularisation

The classical optimal transport problem consists in minimising the total cost
of transportation of a given probability measure $\mu\in\mathcal{P}(X)$ into a
target one $\nu\in\mathcal{P}(Y)$, with respect to some cost function
$c:X\times Y\to\mathbb{R}$, which for simplicity here we assume to be
continuous and bounded, i.e. $c\in C_{b}(X\times Y)$. In the most general
setting, $X$, $Y$ are Polish spaces, i.e. topological spaces that can be
endowed with a distance which makes them complete and separable. We refer to
the seminal book [Vil08] for a general treatment of this subject. The
Kantorovich formulation of the associated optimal transport problem then reads
as

$\displaystyle\inf_{\pi\in\mathcal{P}(X\times Y)}\left\\{\int
c(x,y)\,\mathrm{d}\pi(x,y)\ :\
\big{(}\text{Proj}_{X}\big{)}_{\\#}\pi=\mu\quad\text{and}\quad\big{(}\text{Proj}_{Y}\big{)}_{\\#}\pi=\nu\right\\}\,,$
(1)

where $\text{Proj}_{X}:(X,Y)\to X$ (resp. $\text{Proj}_{Y}:(X,Y)\to Y$)
denotes the projection operator onto $X$ (resp. Y). A probability measure
$\pi\in\mathcal{P}(X\times Y)$ which satisfies the marginal constraint as in
(1) is called an admissible transport plan between $\mu$ and $\nu$. When the
involved marginal measures are atomic (or discrete), i.e. of the form

$\displaystyle\mu=\sum_{i=1}^{L}\mu_{i}\delta_{x_{i}}\quad\text{and}\quad\nu=\sum_{j=1}^{M}\nu_{j}\delta_{y_{j}}\,,$
(2)

then the variational problem (1) is a fully discrete constrained linear
program which can be equivalently recast as

$\displaystyle\inf_{\pi_{ij}\geq
0}\left\\{\sum_{i=1}^{L}\sum_{j=1}^{M}c(x_{i},y_{j})\pi_{ij}\ :\
\sum_{j=1}^{M}\pi_{ij}=\mu_{i}\quad\text{and}\quad\sum_{i=1}^{L}\pi_{ij}=\nu_{j}\right\\}\,.$
(3)

An important feature of the variational problem (1) is the fact it admits a
dual formulation, using the natural duality between the space of measures and
the space continuous and bounded functions. The equivalent, dual formulation
of (1) is given by

$\displaystyle\sup_{\varphi,\psi}\left\\{\int\varphi(x)\,\mathrm{d}\mu(x)+\int\psi(y)\,\mathrm{d}\nu(y)\
:\ \varphi(x)+\psi(y)\leq c(x,y)\right\\}\,,$ (4)

where the supremum is taken over $\varphi\in C_{b}(X)$ and $\psi\in C_{b}(Y)$.

As we are going to see several times in this survey, a natural technique in
order to look for maximising couples $(\varphi,\psi)$ in (4) consists in
fixing one of two functions and maximising the target functional with respect
to the second one. In more precise terms, for a given $\varphi\in C_{b}(X)$,
we define $\varphi^{c}\in C_{b}(Y)$ as the solution to

$\displaystyle\varphi^{c}(y):=\inf_{x\in
X}\big{\\{}c(x,y)-\varphi(x)\big{\\}}=\mathop{\mathrm{argmax}}_{\psi}\left\\{\int\psi(y)\,\mathrm{d}\nu(y)\
:\ \varphi(x)+\psi(y)\leq c(x,y)\right\\}\,.$ (5)

We call $\varphi^{c}$ the $c$-transform of $\varphi$. Similarly, for any given
$\psi\in C_{b}(Y)$, one defines $\psi^{c}\in C_{b}(X)$ by maximising (4) over
$\varphi$, or equivalently

$\displaystyle\psi^{c}(x):=\inf_{y\in Y}\big{\\{}c(x,y)-\psi(y)\big{\\}}\in
C_{b}(X)\,.$ (6)

We say that $\psi$ is c-concave if $\psi=\varphi^{c}$ for some $\varphi$. One
can prove that $\varphi^{ccc}=\varphi^{c}$ and that $\psi$ is c-concave if and
only if $\psi^{cc}=\psi$ (see e.g. [Vil08, Prop. 5.8]). This means that the
optimisers in (4) are necessarily a c-concave function and its c-trasform, or
in other words we can write (4) equivalently as

$\displaystyle\sup_{\varphi\in
C_{b}(X)}\left\\{\int\varphi(x)\,\mathrm{d}\mu(x)+\int\varphi^{c}(y)\,\mathrm{d}\nu(y)\
:\ \varphi\text{ is c-concave}\right\\}\,.$ (7)

The optimal $\varphi$ in the dual problem is usually called Kantorovich
potential.

######  Remark 2.1 (Euclidean spaces).

A widely studied example of optimal transport is the one corresponding to the
euclidean space $X=Y=\mathbb{R}^{d}$ and the quadratic cost
$c(x,y)=|x-y|^{2}$. In this case, the square root of (1) defines a distance on
the space of probability measures with finite second moment known in
literature as Kantorovich–Rubenstein–Wasserstein distance $W_{2}$. One of the
consequence of the duality formula (4) (and in particular the characterisation
of $c$-concave maps) in this framework is the existence of the so-called
optimal transport map $T:\mathbb{R}^{d}\to\mathbb{R}^{d}$, whenever $\mu$ is
absolutely continuous with respect to the Lebesgue measure in
$\mathbb{R}^{d}$. This means that the optimal plan $\pi$ in (1) is of the form
$\pi=(\text{id},T)_{\\#}\mu$, which corresponds to the original formulation of
the problem, as introduced by Monge [Mon81] back in the 18th century. In
particular, $\pi$ is supported on the graph of the map $T$. Note that this
phenomenon is quite in contrast with what happens in the fully discrete case
(3), where in most cases no admissible plan is induced by a map $T$.

#### Multimarginal Optimal transport

Up to this point, we have discussed two marginals problems, namely where the
optimal transport problem happens between two spaces and two probability
measures $\mu$ and $\nu$. A natural, interesting generalisation are the
multimarginal optimal transport problems, where instead of an initial and a
target space, one considers $N$ spaces $X_{1},\dots,X_{N}$ and $N$ given
marginals $\\{\mu_{i}\in\mathcal{P}(X_{i})\,:\,i=1,\dots,N\\}$. This time a
cost function is given by a map from the product of all the spaces, namely

$\displaystyle\mathbf{X}:=\bigtimes_{i=1}^{N}X_{i}\quad\text{and}\quad
c:\mathbf{X}\to\mathbb{R}\quad\text{continuous and bounded}\,.$ (8)

Using the notation $\mathbf{x}:=(x_{1},\dots,x_{N})$, the variational problem
associated to $c$ and $\\{\mu_{i}\\}_{i}$ is then given by

$\displaystyle\inf_{\pi\in\mathcal{P}(\mathbf{X})}\left\\{\int
c(\mathbf{x})\,\mathrm{d}\pi(\mathbf{x})\ :\
\big{(}\text{Proj}_{X_{i}}\big{)}_{\\#}\pi=\mu_{i}\quad\text{for}\
i=1,\dots,N\right\\}\,.$ (9)

In this picture, we lose the original interpretation of transport from a given
origin to a target location. Nonetheless, in some special cases, one can
reduce this formulation to the ”Monge” one, by considering admissible
multimarginal plans associated with $N-1$ maps, of the form
$\pi=(\text{id},T_{1},\dots,T_{N-1})_{\\#}\mu_{1}$, see e.g. [GS98]. In a
similar spirit as for the two-marginals case, the primal problem (9) admits a
dual formulation, which is given by

$\displaystyle\sup\left\\{\sum_{i=1}^{N}\int_{X_{i}}\varphi_{i}(x_{i})\,\mathrm{d}\mu_{i}(x_{i})\
:\ \varphi_{j}\in
C_{b}(X_{j})\quad\text{and}\quad\sum_{i=1}^{N}\varphi_{i}(x_{i})\leq
c(\mathbf{x})\right\\}\,.$ (10)

Moreover, we can define the associated $c$-transform operators: this time,
they take as input $N-1$ functions, and find the optimal $N$-th function that
satisfies the cost constraint. More precisely, for a fixed index
$j\in\\{1,\dots,n\\}$, the $c$-transform of the continuous functions
$(\varphi_{1},\dots,\varphi_{j-1},\varphi_{j+1},\dots,\varphi_{N})$ is the
continuous function given by

$\displaystyle(\varphi_{1},\dots,\varphi_{j-1},\varphi_{j+1},\dots,\varphi_{N})^{c}(x_{j}):=\inf_{\mathbf{x}^{(j)}}\left\\{c(\mathbf{x})-\sum_{i\neq
j}\varphi_{i}(x_{i})\right\\}\,,$ (11)

where we used the notation
$\mathbf{x}^{(j)}:=(x_{1},\dots,x_{j-1},x_{j+1},\dots,x_{N})\in\bigtimes_{i\neq
j}X_{i}$.

#### Regularisation of optimal transport

Especially in the two marginal case, the optimal transport plans can be very
singular with respect to diffused measures on the product, for example the
Lebesgue measure on $\mathbb{R}^{d}\times\mathbb{R}^{d}$. They are instead
often concentrated on the graph of functions, and therefore less ”diffused”
(cfr. Remark 2.1). In order to obtain less singular optimisers, one
possibility is to introduce a penalisation term in the minimisation that
forces the plan to be absolutely continuous with respect to a given reference
measure. A typical energy that penalises singular measures with respect to a
given reference is the relative entropy functional. For
$\sigma\in\mathcal{P}(\mathbf{X})$ a given reference measure, the relative
entropy w.r.t. $\sigma$ is defined as

$\displaystyle\text{Ent}(\cdot|\sigma):\mathcal{P}(\mathbf{X})\to[0,+\infty]\,,\qquad\text{Ent}(\pi|\sigma):=\begin{cases}\int\rho\log\rho\,\mathrm{d}\sigma&\text{if
}\,\mathrm{d}\pi=\rho\,\mathrm{d}\sigma\,,\\\
+\infty&\text{otherwise}\,.\end{cases}$ (12)

In particular, if ${\rm Ent}(\pi|\sigma)<\infty$, then $\pi$ is necessarily
absolutely continuous with respect to $\sigma$. Note that we are adopting the
mathematical convention, which means that the entropy of a probability measure
is a positive number (and not negative, as more commonly assumed in physics).
For a given parameter $\varepsilon\in(0,+\infty)$, we can then define a
$\varepsilon$-regularised multimarginal optimal transport problem by adding an
entropy penalisation, that is

$\displaystyle\inf_{\pi\in\mathcal{P}(\mathbf{X})}\left\\{\int
c(\mathbf{x})\,\mathrm{d}\pi(\mathbf{x})+\varepsilon\text{Ent}(\pi|\sigma)\ :\
\big{(}\text{Proj}_{X_{i}}\big{)}_{\\#}\pi=\mu_{i}\quad\text{for}\
i=1,\dots,N\right\\}\,.$ (13)

In order to have a well-posed problem, we assume
$\sigma=\bigotimes_{i}\mathfrak{m}_{i}$ and
$\text{Ent}(\mu_{j}|\mathfrak{m}_{j})<\infty$, for every $j$. Note that the
minimiser of the regularised problem (13) with $\varepsilon>0$ must have
finite entropy, therefore it is necessarily absolutely continuous with respect
to $\sigma$. Furthermore, by means of a simple computation, we see that

$\displaystyle\int
c(\mathbf{x})\,\mathrm{d}\pi(\mathbf{x})+\varepsilon\text{Ent}(\pi|\sigma)$
$\displaystyle=\varepsilon\int\rho\big{(}\log
e^{\frac{c}{\varepsilon}}+\log\rho\big{)}\,\mathrm{d}\sigma$ (14)
$\displaystyle=\varepsilon\int e^{\frac{c}{\varepsilon}}\rho\log(\rho
e^{\frac{c}{\varepsilon}})e^{-\frac{c}{\varepsilon}}\,\mathrm{d}\sigma=\varepsilon\text{Ent}(\pi|e^{-\frac{c}{\varepsilon}}\sigma)\,,$
(15)

which in particular shows that we can equivalently recast (13) as

$\displaystyle\inf_{\pi\in\mathcal{P}(\mathbf{X})}\left\\{\text{Ent}(\pi|e^{-\frac{c}{\varepsilon}}\sigma)\
:\ \big{(}\text{Proj}_{X_{i}}\big{)}_{\\#}\pi=\mu_{i}\quad\text{for}\
i=1,\dots,N\right\\}\,.$ (16)

In other words, the $\varepsilon$-regularised optimal transport problem is
equivalent to finding the minimiser of the relative entropy with respect to
the weighted measure
$\mathcal{K}_{\varepsilon}:=e^{-\frac{c}{\varepsilon}}\sigma$. This is known
in the mathematical literature as Schrödinger problem associated with $c$ and
$\sigma$, and $\mathcal{K}_{\varepsilon}$ is often referred to as Gibbs
measure.

Even with a positive regularisation parameter $\varepsilon>0$, the problem
(13) admits a dual formulation. This time, the presence of a superlinear term
(the entropy) causes a particular phenomenon: the dual problem is not a
linearly constrained problem (as for $\varepsilon=0$) where the cost appears
in the constraint. On the contrary, it is obtained as a free maximisation
problem where the marginals appear linearly in the dual functional, whereas
the cost appears as part of a (concave) nonlinear term. The dual formulation
of (13) indeed reads as

$\displaystyle\sup\left\\{\sum_{i=1}^{N}\int_{X_{i}}\varphi_{i}(x_{i})\,\mathrm{d}\mu_{i}(x_{i})-\varepsilon\int_{\mathbf{X}}\exp\left\\{\frac{1}{\varepsilon}\bigg{(}\sum_{i=1}^{N}\varphi_{i}(x_{i})-c(\mathbf{x})\bigg{)}\right\\}\,\mathrm{d}\sigma(\mathbf{x})\right\\}+\varepsilon\,,\
$ (17)

where the supremum runs over all $\varphi_{i}\in C_{b}(X_{i})$. At least
formally, by sending $\varepsilon\to 0$, we recover the dual formula (10).
This is related to the so-called Laplace principle, which informally states
that

$\displaystyle\int_{\mathbf{X}}\exp\left\\{\frac{1}{\varepsilon}\bigg{(}\sum_{i=1}^{N}\varphi_{i}(x_{i})-c(\mathbf{x})\bigg{)}\right\\}\,\mathrm{d}\sigma(\mathbf{x})\underset{\varepsilon\to
0}{\approx}\exp\left\\{\frac{1}{\varepsilon}\sigma\text{-esssup}\bigg{(}\sum_{i=1}^{N}\varphi_{i}(x_{i})-c(\mathbf{x})\bigg{)}\right\\}\,,$
(18)

showing that, in the limit as $\varepsilon\to 0$, in order to have a
competitor which has energy value which is not $-\infty$, one must necessarily
enforce the constraint as in (10).

######  Remark 2.2 (Schrödinger problem in $\mathbb{R}^{d}$).

Let us consider the Euclidean and two-marginals setup, i.e. $N=2$,
$X=Y=\mathbb{R}^{d}$, and consider $\sigma:=\,\mathrm{d}x\,\mathrm{d}y$ the
Lebesgue measure on the product, together with the quadratic cost (up to
translation by a constant)
$c(x,y)=\frac{1}{2}|x-y|^{2}-\log\sqrt{2\pi\varepsilon^{2}}$. In this case,
the Gibbs measure coincides with the Gaussian kernel
$\mathcal{K}_{\varepsilon}=\frac{1}{\sqrt{2\pi\varepsilon^{2}}}e^{-\frac{|x-y|^{2}}{2\varepsilon}}\,\mathrm{d}x\,\mathrm{d}y$
and

$\displaystyle\inf_{\pi}\big{\\{}\text{Ent}\big{(}\pi|\mathcal{K}_{\varepsilon}\big{)}\
:\
\big{(}\text{Proj}_{X_{i}}\big{)}_{\\#}\pi=\mu_{i}\big{\\}}=\inf_{P}\big{\\{}\text{Ent}\big{(}P|W_{\varepsilon}\big{)}\
:\ P(\cdot|t=i)=\mu_{i}\big{\\}}\,,$

where $P\in\mathcal{P}\big{(}C([0,1];\mathbb{R}^{d})\big{)}$ and
$W_{\varepsilon}\in\mathcal{P}\big{(}C([0,1];\mathbb{R}^{d})\big{)}$ is the
Wiener measure. This represents a bridge between the theory of (regularised)
optimal transport and Large Deviations. Indeed, as described by Sanov’s
Theorem, the relative entropy functional $\text{Ent}(\cdot|W_{\varepsilon})$
naturally appears as rate functional for the large deviation principle
satisfied by i.i.d. random variables equally distributed as $W_{\varepsilon}$.

As a typical feature of primal-dual problems, we have that the minimisers and
maximisers of the respective variational problems are related to each other.
In this case, we can write an explicit formula for the minimal coupling
$\pi_{\varepsilon}\in\mathcal{P}(\mathbf{X})$ in (13) in terms of the optimal
Kantorovich potentials
$(\varphi_{\varepsilon,1},\dots,\varphi_{\varepsilon,N})$,
$\varphi_{\varepsilon,i}\in C_{b}(X_{i})$ in (17). It turns out that
$\pi_{\varepsilon}$ is given the absolutely continuous measure with respect to
$\mathcal{K}_{\varepsilon}$ with density given by

$\displaystyle\frac{\,\mathrm{d}\pi_{\varepsilon}}{\,\mathrm{d}\mathcal{K}_{\varepsilon}}(\mathbf{x})=\exp\left(\frac{1}{\varepsilon}\left(\sum_{i=1}^{N}\varphi_{\varepsilon,i}(x_{i})-c(\mathbf{x})\right)\right)=\prod_{i=1}^{N}\exp\left(\frac{\varphi_{\varepsilon,i}(x_{i})-N^{-1}c(\mathbf{x})}{\varepsilon}\right)\,.$
(19)

In particular, once evaluated on $\pi_{\varepsilon}$, the nonlinear
(exponential) integral term in (17) simply gives $1$ (cause
$\pi_{\varepsilon}$ is a probability measure) and value of the primal (13) and
dual (17) problem coincides, and it is given by

$\displaystyle\sum_{i=1}^{N}\int_{X_{i}}\varphi_{\varepsilon,i}(x_{i})\,\mathrm{d}\mu_{i}(x_{i})=\int
c(\mathbf{x})\,\mathrm{d}\pi_{\varepsilon}(\mathbf{x})+\varepsilon\text{Ent}(\pi|\sigma)\,.$
(20)

The entropic regularised problem also admits a natural notion of
$c$-transform, that the authors in [DMG20] call $(c,\varepsilon)$-trasform. In
a similar way as for the case $\varepsilon=0$, it is a trasformation which
corresponds to optimising over a single potential, say $\varphi_{j}$, while
fixing the other $N-1$, that is $\\{\varphi_{i}\,:\,i\neq j\\}$. In the
regularised setting ($\varepsilon>0$), one can write explicitly the solution
to this maximisation problem. In particular, for a fixed index
$j\in\\{1,\dots,n\\}$, the $(c,\varepsilon)$-transform of the continuous
functions $(\varphi_{1},\dots,\varphi_{j-1},\varphi_{j+1},\dots,\varphi_{N})$
is the continuous function denoted by
$(\varphi_{1},\dots,\varphi_{j-1},\varphi_{j+1},\dots,\varphi_{N})^{(c,\varepsilon)}$
and given by

$\displaystyle x_{j}\mapsto-\varepsilon\log\left(\int_{\bigtimes_{i\neq
j}X_{i}}\exp\left(\frac{1}{\varepsilon}\left(\sum_{i\neq
j}\varphi_{i}(x_{i})-c(\mathbf{x})\right)\right){\rm d}\bigotimes^{N}_{i\neq
j}\mathfrak{m}_{i}\right)+\varepsilon\log\Big{(}\frac{\,\mathrm{d}\mu_{j}}{\,\mathrm{d}\mathfrak{m}_{j}}\Big{)}\,.$
(21)

Once again, by the Laplace principle, in the limit as $\varepsilon\to 0$, one
recovers (at least formally) the definition of $c$-transform as given in (11).

#### Sinkhorn algorithm

In order to describe the Sinkhorn algorithm, let us start the discussion in
the case of the two-marginals problems. Namely, for $\varepsilon>0$, we seek
solutions to

$\displaystyle\inf_{\pi\in\mathcal{P}(X\times
Y)}\left\\{\text{Ent}\big{(}\pi|\mathcal{K}_{\varepsilon}\big{)}\ :\
\big{(}\text{Proj}_{X}\big{)}_{\\#}\pi=\mu\quad\text{and}\quad\big{(}\text{Proj}_{Y}\big{)}_{\\#}\pi=\nu\right\\}\,,\quad\mathcal{K}_{\varepsilon}:=e^{\frac{-c}{\varepsilon}}\,\mathrm{d}\sigma\,.$
(22)

The key idea behind Sinkhorn algorithm is to solve alternating optimisations
problems in order to find a sequence of plans $\\{\pi^{(n)}\\}_{n}$ having one
correct marginal at a time, in such a way that $\pi^{(n)}\to\pi_{\varepsilon}$
(the optimal solution) as $n\to\infty$. Moreover, from the optimality
conditions described in (19), we know that the optimiser must be in product
form. In fact, one can prove that it completely characterises the minimiser of
the primal problem: if $\pi_{\varepsilon}$ has the right marginals, and it is
in product form, then it must be optimal. Driven by these considerations, one
constructs the following algorithm.

Initialisation. We initialise $\pi^{(0)}=(f^{(0)}\otimes
g^{(0)})\mathcal{K}_{\varepsilon}\in\mathcal{P}(X\times Y)$ in product form.

Iteration. We iteratively define $\pi^{(n)}=(f^{(n)}\otimes
g^{(n)})\mathcal{K}_{\varepsilon}\in\mathcal{P}(X\times Y)$ by imposing the
correct marginals $\mu$, $\nu$ one at a time, i.e. given $f^{(n-1)}$,
$g^{(n-1)}$, we solve

$\displaystyle(\text{Proj}_{X})_{\\#}(f^{(n)}\otimes
g^{(n-1)})\mathcal{K}_{\varepsilon}=\mu\,,\quad(\text{Proj}_{Y})_{\\#}(f^{(n)}\otimes
g^{(n)})\mathcal{K}_{\varepsilon}=\nu\,.$ (23)

Explicitly, for $\,\mathrm{d}\mu=\rho\,\mathrm{d}\mathfrak{m}_{1}$,
$\,\mathrm{d}\nu=\eta\,\mathrm{d}\mathfrak{m}_{2}$, the iterations are given
by

$\displaystyle
f^{(n)}(x)=\frac{\rho(x)}{\int_{Y}g^{(n-1)}(y)e^{-c(x,y)/\varepsilon}\,\mathrm{d}\mathfrak{m}_{2}(y)}\,,\quad
g^{(n)}(y)=\frac{\eta(y)}{\int_{X}f^{(n)}(x)e^{-c(x,y)/\varepsilon}\,\mathrm{d}\mathfrak{m}_{1}(x)}\,.\quad$
(24)

In this setting, one can prove the following result: assume that $c\in
L^{\infty}(X\times Y)$. Up to subsequence, $f^{(n)}/\lambda^{(n)}\to
f_{\varepsilon}$ in $L^{1}(\mathfrak{m}_{0})$ and $g^{(n)}/\lambda^{(n)}\to
g_{\varepsilon}$ in $L^{1}(\mathfrak{m}_{1})$, for some
$\\{\lambda^{(n)}>0\\}_{n}$. In particular, $\pi^{(n)}\to\pi_{\varepsilon}$ in
$L^{1}(\sigma)$, where $\pi_{\varepsilon}$ is optimal for (13). The
convergence of this procedure was first provided by J. Franklin and J. Lorenz
[FL89] in the discrete case, via the introduction of a Hilbert metric on the
projection cone of the Sinkhorn iterations, and by L. Ruschendorf [Rus95] in
the continuous one. An alternative proof based on the Franklin-Lorenz approach
was also provided by Y. Chen, T. Georgiou and M. Pavon [CGP16], which is
particular yields a linear convergence rate of the algorithm (in the Hilbert
metric).

It turns out that solving (23) is equivalent to compute the
$(c,\varepsilon)$-transform of the previous step. More precisely, one can
write

$\displaystyle f^{(n)}:=\exp(\varphi^{(n)})\,,\
g^{(n)}:=\exp(\psi^{(n)})\,,\quad\text{and}\ \ \
\varphi^{(n)}:=\big{(}\psi^{(n-1)}\big{)}^{(c,\varepsilon)}\,,\
\psi^{(n)}:=\big{(}\varphi^{(n)}\big{)}^{(c,\varepsilon)}.$ (25)

In the work [DMG20], the authors proved the convergence the Sinkhorn algorithm
in the setting of multimarginal optimal transport. In this case, one
initialise as
$\pi^{(0)}=\big{(}\bigotimes_{j=1}^{N}f_{j}\big{)}\mathcal{K}_{\varepsilon}$
and consider the iterations
$\pi^{(n)}=\big{(}\bigotimes_{j=1}^{N}f_{j}^{(n)}\big{)}\mathcal{K}_{\varepsilon}$
where the sequence
$\\{f_{j}^{(n)}:=\exp(\varphi_{j}^{(n)})\,:\,j=1,\dots,N\\}_{n}$ is defined by

$\displaystyle(\text{Proj}_{X_{i}})_{\\#}\Big{(}\bigotimes_{i=1}^{j}f_{j}^{(n)}\bigotimes_{i=j+1}^{N}f_{j}^{(n-1)}\Big{)}\mathcal{K}_{\varepsilon}=\mu_{i}\,\Leftrightarrow\,\varphi_{j}^{(n)}=\big{(}\varphi_{1}^{(n)},...\,,\varphi_{j-1}^{(n)},\varphi_{j+1}^{(n-1)},...\big{)}^{(c,\varepsilon)}\,,$
(26)

using the notion of $c$-transform in the sense of (21). In particular, it has
been showed in [DMG20] that for bounded cost functions, there exist
renormalisation parameters $\\{\lambda_{j}^{(n)}>0\,:\,j=1,\dots,N\\}$ so that
$f_{j}^{(n)}/\lambda_{j}^{(n)}\to\exp(\varphi_{\varepsilon,j})$ in
$L^{1}(\mathfrak{m}_{j})$ as well as $\pi^{(n)}\to\pi_{\varepsilon}$ in
$L^{1}(\sigma)$. In particular, $\pi_{\varepsilon}$ is optimal for (13),
$(\varphi_{\varepsilon,1},\dots,\varphi_{\varepsilon,N})$ is optimal for the
dual problem (17), and they are linked by the formula

$\displaystyle\pi_{\varepsilon}=\Big{(}\bigotimes_{i=1}^{N}\exp(\varphi_{\varepsilon,i})\Big{)}\mathcal{K}_{\varepsilon}\,.$
(27)

Note one can write an explicit formula for each iteration $\varphi_{j}^{(n)}$
in terms of the previous ones, in the same spirit of (24). As we are going to
see, this is in constrast with the noncommutative setting, where no explicit
formulas are available (which makes the study of the convergence of the
algorithm more involved).

##  3 Noncommutative optimal transport

In this section, we are going to discuss a noncommutative analog of the
multimarginal optimal transport problem (9). As in [FGP23], we only focus on
the finite dimensional case, i.e. a noncommutative version of the discrete
optimal transport problem as in (3).

A short discussion about the notation is due: for a given Hilbert space
$\mathfrak{g}$, we denote by $\mathcal{M}(\mathfrak{g})$ the set of all
bounded operators over $\mathfrak{g}$, by $\mathcal{S}(\mathfrak{g})$ the
hermitian elements of $\mathcal{M}(\mathfrak{g})$, and by
$\mathcal{S}_{\geq}(\mathfrak{g})$ (resp. $\mathcal{S}_{>}(\mathfrak{g})$) the
set of all the positive semidefinite (resp. positive definite) elements of
$\mathcal{S}(\mathfrak{g})$. With a slight abuse of notation, we denote by
$\operatorname{Tr}$ the trace operator on $\mathcal{M}(\mathfrak{g})$ for any
Hilbert space $\mathfrak{g}$. Furthemore, we denote by
$\mathfrak{P}(\mathfrak{g})$ the set of density matrices over $\mathfrak{g}$,
namely the elements of $\mathcal{S}_{\geq}(\mathfrak{g})$ with trace one. We
also set $[N]:=\\{1,...,N\\}$.

Instead of considering transport problems in metric spaces, the noncommutative
framework consists of $N$ complex Hilbert spaces $\mathfrak{h}_{j}$ of
dimension $d_{j}<\infty$, for $j\in[N]$, which physically represent the state
spaces of $N$ subsystems. The state space of the composite system is given by
$\mathfrak{h}:=\mathfrak{h}_{1}\otimes\mathfrak{h}_{2}\otimes\dots\otimes\mathfrak{h}_{N}$,
which is also a finite dimensional Hilbert space of dimension $d=d_{1}\cdot
d_{2}\cdot\dots\cdot d_{N}$, and plays the role of $\mathbf{X}$ in the
classical picture. Concerning the other quantities involved in (9), we have
that:

  1. 1.

The role of the cost function $c:\mathbf{X}\to\mathbb{R}$ is played by a
Hamiltonian $\operatorname{H}\in\mathcal{S}(\mathfrak{h})$.

  2. 2.

The set of probability measures over $\mathbf{X}$ are replaced by density
matrices $\mathfrak{P}(\mathfrak{h})$. In particular, having total mass equals
one becomes having trace one.

  3. 3.

The notion of marginal via the projection operators $\text{Proj}_{X_{i}}$ is
replaced by the one of reduced density matrix. We say that
$\Gamma\in\mathfrak{P}(\mathfrak{h})$ has marginals
$\bm{\gamma}:=(\gamma_{1},..,\gamma_{N})$ if

$\displaystyle\operatorname{Tr}(\gamma_{i}A)=\operatorname{Tr}\Big{(}\Gamma\big{(}1\\!\\!1\otimes...\otimes
1\\!\\!1\otimes\underbrace{A}_{\text{i-th}}\otimes 1\\!\\!1..\otimes
1\\!\\!1\big{)}\Big{)},\quad\forall A\in\mathcal{M}(\mathfrak{h}_{i}).$ (28)

We use the short-hand notation $\Gamma\mapsto\bm{\gamma}$. In general, if
$\Gamma$ satisfies (28) for a given $i$, we write
$\operatorname{P}_{i}(\Gamma)=\gamma_{i}$. Here
$\operatorname{P}_{i}:\mathfrak{P}(\mathfrak{h})\to\mathfrak{P}(\mathfrak{h}_{i})$
is the $i$-th partial trace operator, which plays the same role of the
projection $(\text{Proj}_{X_{i}})_{\\#}(\cdot)$ in the commutative setting.

The noncommutative analog of the multimarginal transport problem (9) is given
by

$\displaystyle\inf\left\\{\operatorname{Tr}(\operatorname{H}\Gamma)\ :\
\Gamma\in\mathfrak{P}(\mathfrak{h})\text{ and
}\Gamma\mapsto\bm{\gamma}\right\\}\,,$ (29)

for some given marginal operators (or reduced density matrices)
$\\{\gamma_{i}\in\mathfrak{P}(\mathfrak{h}_{i})\\}_{i\in[N]}$.

Note that without any constraint on the marginals, the optimal density matrix
is a pure state (i.e. rank one projection) onto the eigenspace of
$\operatorname{H}$ corresponding to its minimal eigenvalue. More precisely,
using the spectral theorem on $\operatorname{H}$, we can always write

$\displaystyle\operatorname{H}=\sum_{i=1}^{d}\lambda_{i}|\psi_{i}\rangle\langle\psi_{i}|\,,\quad\text{where}\
\sigma(\operatorname{H})=\\{\lambda_{i}\\}_{i=1}^{d}\ \text{is the spectrum of
}\operatorname{H}\,,$ (30)

where $\psi_{i}\in\mathfrak{h}$ are the eigenfunctions of $\operatorname{H}$
(normalised to have norm one, i.e. $\langle\psi_{i}|\psi_{i}\rangle=1$, for
every $i$). We are using the bracket notation
$|\psi_{i}\rangle\langle\psi_{i}|$ to denote the rank-one projection onto the
linear subspace spanned by $\psi_{i}$. The solution to the unconstrained
problem of minimising $\operatorname{Tr}(\Gamma\operatorname{H})$ is then
given by the pure state $\Gamma=|\psi_{i_{0}}\rangle\langle\psi_{i_{0}}|$,
where $i_{0}$ corresponds to the minimal eigenvalue
$\lambda_{i_{0}}=\min_{i}\lambda_{i}$ of $\operatorname{H}$. This
concentration phenomenon is close in spirit to what discussed in Remark 2.1,
which is typical of the unregularised setting. In general, for given marginals
$\bm{\gamma}$, the rank-one projection onto the span of $\psi_{i_{0}}$ might
not be admissible.

#### Noncommutative entropic optimal transport

From a physical perspective, one can study the system at zero or at positive
temperature. It turns out that the quantum problem at positive temperature
$\varepsilon>0$ corresponds to a regularised version of the noncommutative
optimal transport problem as introduced in (29). The regularisation is
performed once again by adding an entropy term, which penalises concentration,
namely for which pure states are not optimal anymore. For a given reference
positive operator $\bm{m}\in\mathcal{S}_{>}(\mathfrak{h})$, the relative
entropy as introduced in (12) is here replaced by (minus) the Umegaki relative
entropy of $\Gamma$, which is given by

$\displaystyle
S(\Gamma|\bm{m}):=\operatorname{Tr}(\Gamma(\log\Gamma-\log\bm{m}))\,,\quad\Gamma\in\mathfrak{P}(\mathfrak{h})\,.$
(31)

In the simpler case when $\bm{m}=1\\!\\!1$, we obtain (minus) the Von Neumann
entropy functional. Using the spectral decomposition of $\Gamma$, we can
compute it as

$\displaystyle
S(\Gamma):=S(\Gamma|1\\!\\!1)=\sum_{i=1}^{d}\lambda_{i}\log\lambda_{i}\,,\quad\text{if}\quad\Gamma=\sum_{i=1}^{d}\lambda_{i}|\psi_{i}\rangle\langle\psi_{i}|\,,$
(32)

which shows that $S(\Gamma)$ can be interpreted as the entropy functional (12)
evaluated in the empirical measure with support on the spectrum of $\Gamma$
relative to the counting measure.

In general, we assume $\bm{m}$ to be of product form, that is
$\bm{m}=\otimes_{i=1}^{N}m_{i}$, for some
$\\{m_{i}\in\mathcal{S}_{>}(\mathfrak{h}_{i})\\}$. The entropic noncommutative
optimal transport problem with parameter $\varepsilon>0$ is then given by

$\displaystyle\mathfrak{F}^{\varepsilon}_{\bm{m}}(\bm{\gamma}):=\inf\left\\{\operatorname{Tr}(\operatorname{H}\Gamma)+\varepsilon
S(\Gamma|\bm{m})\ :\ \Gamma\in\mathfrak{P}(\mathfrak{h})\text{ and
}\Gamma\mapsto\bm{\gamma}\right\\}\,,$ (33)

for a given set of marginals
$\\{\gamma_{i}\in\mathfrak{P}(\mathfrak{h}_{i})\\}_{i\in[N]}$.

######  Remark 3.1 (Many-body problems).

Typically in quantum many-body problems, the Hamiltonian $\operatorname{H}$ to
which the whole system is subject is of the form
$\operatorname{H}=\operatorname{H}_{0}+\operatorname{H}_{\text{int}}$, where
$\operatorname{H}_{0}$ is the non-interacting part of the Hamiltonian, i.e.
$\operatorname{H}_{0}=\bigoplus_{j=1}^{N}\operatorname{H}_{j}:=\operatorname{H}_{1}\otimes
1\\!\\!1\dots\otimes 1\\!\\!1+1\\!\\!1\otimes\operatorname{H}_{2}\otimes
1\\!\\!1\dots\otimes 1\\!\\!1+\dots+1\\!\\!1\otimes\dots\otimes
1\\!\\!1\otimes\operatorname{H}_{N}$, with $\operatorname{H}_{j}$ acting on
$\mathfrak{h}_{j}$, and $\operatorname{H}_{\text{int}}$ is its interacting
part. Now, suppose we have knowledge of the states
$\bm{\gamma}=(\gamma_{1},\dots,\gamma_{N})$ of the $N$ subsystems. Then the
energy of the composite system at temperature $\varepsilon>0$ is given by

$\displaystyle\inf_{\Gamma\mapsto\bm{\gamma}}\left\\{\operatorname{Tr}(\operatorname{H}\Gamma)+\varepsilon
S(\Gamma)\right\\}=\sum_{j=1}^{N}\operatorname{Tr}(\operatorname{H}_{j}\gamma_{j})+\inf_{\Gamma\mapsto\bm{\gamma}}\left\\{\operatorname{Tr}(\operatorname{H}_{\text{int}}\Gamma)+\varepsilon
S(\Gamma)\right\\}\,,$ (34)

where we clearly see in the second term of the right-hand side an entropic
problem (33).

######  Remark 3.2 (Lack of concentration).

Denote by $\bar{\Gamma}_{\varepsilon}\in\mathfrak{P}(\mathfrak{h})$ the
density matrix which minimises the energy in (33) with $\bm{m}=1\\!\\!1$
without the marginals constraints. We claim that $\bar{\Gamma}_{\varepsilon}$
is never a pure state; in fact, we shall show that
$\bar{\Gamma}_{\varepsilon}$ has trivial kernel. Indeed, suppose that
$\ker\bar{\Gamma}_{\varepsilon}\neq\\{0\\}$ and denote by
$\Pi_{\ker\bar{\Gamma}_{\varepsilon}}$ the projection onto it. Define
$\tilde{\Gamma}_{s}:=\bar{\Gamma}_{\varepsilon}+s\Pi_{\ker\bar{\Gamma}_{\varepsilon}}$
for $s>0$, and observe that $\tilde{\Gamma}_{s}$ has the same eigenspaces of
$\bar{\Gamma}_{\varepsilon}$ but eigenvalue $s>0$ for every element in
$\ker\bar{\Gamma}_{\varepsilon}$. Computing the entropy, we see that

$\displaystyle S(\tilde{\Gamma}_{s})=S(\bar{\Gamma}_{\varepsilon})+s\log
s\,d_{\ker}\,,\quad\text{where}\quad
d_{\ker}:=\dim{\ker\bar{\Gamma}_{\varepsilon}}\,.$

Therefore, by construction we obtain that

$\displaystyle\operatorname{Tr}(\operatorname{H}\tilde{\Gamma}_{s})-\varepsilon
S(\tilde{\Gamma}_{s})$
$\displaystyle=\operatorname{Tr}(\operatorname{H}\bar{\Gamma}_{\varepsilon})-\varepsilon
S(\bar{\Gamma}_{\varepsilon})+s\left(\operatorname{Tr}(\operatorname{H}\Pi_{\ker\bar{\Gamma}_{\varepsilon}})-\varepsilon\log
s\,d_{\ker}\right)$
$\displaystyle>\operatorname{Tr}(\operatorname{H}\bar{\Gamma}_{\varepsilon})-\varepsilon
S(\bar{\Gamma}_{\varepsilon})$

as soon as $s<\exp\left((\varepsilon
d_{\ker})^{-1}\operatorname{Tr}(\operatorname{H}\Pi_{\ker\bar{\Gamma}_{\varepsilon}})\right)\in\mathbb{R}_{+}$,
which contradicts the optimality of $\bar{\Gamma}_{\varepsilon}$.

#### Dual formulation in the noncommutative setting

Another common feature between the classical and quantum setting is the
existence of an equivalent, dual formulation of the problem (33). To simplify
the discussion, we assume from now on that $\ker\gamma_{i}=\\{0\\}$, for every
$i\in[N]$. In the general case, it suffices to consider the restriction to the
set $\mathcal{O}:=\bigotimes_{i=1}^{N}\big{(}\ker\gamma_{i}\big{)}^{\perp}$,
see [FGP23, Remark 3.9]. In this case, the corresponding dual problem is
defined as

$\displaystyle\mathfrak{D}^{\varepsilon}_{\bm{m}}(\bm{\gamma})=\sup\Bigg{\\{}\sum_{i=1}^{N}\operatorname{Tr}(U_{i}\gamma_{i})-\varepsilon\operatorname{Tr}\bigg{(}\exp\bigg{[}\frac{\bigoplus_{i=1}^{N}U_{i}-\operatorname{H}_{\bm{m}}}{\varepsilon}\bigg{]}\bigg{)}\
:\ U_{i}\in\mathcal{S}(\mathfrak{h}_{i})\Bigg{\\}}+\varepsilon\,,\quad$ (35)

where $\operatorname{H}_{\bm{m}}:=\operatorname{H}-\varepsilon\log\bm{m}$ and
$\bigoplus$ denotes the Kronecker sum of operators, defined as

$\displaystyle\bigoplus_{i=1}^{N}U_{i}:=\sum_{i=1}^{N}\big{(}1\\!\\!1\otimes...\otimes
1\\!\\!1\otimes\underbrace{U_{i}}_{\text{i-th}}\otimes 1\\!\\!1..\otimes
1\\!\\!1\big{)}\,.$ (36)

The first result proved in [FGP23] is the following duality result.

######  Theorem 3.3 (Duality).

For any given
$\bm{\gamma}=(\gamma_{i}\in\mathfrak{P}(\mathfrak{h}_{i}))_{i\in[N]}$,
$\operatorname{H}\in\mathcal{S}(\mathfrak{h})$, we have that

  * (i)

the primal and dual problems coincide, i.e.
$\mathfrak{F}^{\varepsilon}_{\bm{m}}(\bm{\gamma})=\mathfrak{D}^{\varepsilon}_{\bm{m}}(\bm{\gamma})$.

  * (ii)

there exists a maximiser
$(U_{i}^{\varepsilon}\in\mathcal{S}^{\hat{d}_{i}})_{i=1}^{N}$ in (35), which
is unique up to trival translation. More precisely, if
$(\tilde{U}_{i}^{\varepsilon}\in\mathcal{S}^{\hat{d}_{i}})_{i=1}^{N}$ is
another maximiser, then
$\tilde{U}_{i}^{\varepsilon}-U_{i}^{\varepsilon}=\alpha_{i}1\\!\\!1$ with
$\alpha_{i}\in\mathbb{R}$ and $\sum_{i}\alpha_{i}=0$.

  * (iii)

There exists a unique ${\Gamma}^{\varepsilon}\in\mathfrak{P}^{\bm{d}}$ with
${\Gamma}^{\varepsilon}\mapsto\bm{\gamma}$ minimiser in (33). Moreover,
${\Gamma}^{\varepsilon}$ and $(U_{i}^{\varepsilon})_{i}$ are related via the
formula

$\displaystyle{\Gamma}^{\varepsilon}=\exp\left(\frac{\bigoplus_{i=1}^{N}U_{i}^{\varepsilon}-\operatorname{H}_{\bm{m}}}{\varepsilon}\right)\,.$
(37)

One should compare the shape of the optimisers given in (37) with the explicit
formula given in (19). The main difference in the noncommutative setting is
that, a priori, we cannot write the optimiser in product form. Indeed, we have
that

$\displaystyle\exp(A+B)=\exp(A)\exp(B)\qquad\text{if}\qquad[A,B]:=AB-BA=0\,,$
(38)

but it is not clear whether the optimal potentials $U_{i}^{\varepsilon}$
commute (i.e. $[U_{i}^{\varepsilon},\operatorname{H}]=0$) with the Hamiltonian
$\operatorname{H}$. This observation is crucial when one wants to construct a
noncommutative version of the Sinkhorn algorithm: we do not expect anymore to
be able to write explicit iterations as in (24), hence also direct
manipulations are not available anymore, already at the level of the
definition of a noncommutative notion of the $(c,\varepsilon)$-transform.

#### A noncommutative $(\operatorname{H},\varepsilon)$-transform and the
Sinkhorn algorithm

In this section, we discuss the noncommutative version of the
$(c,\varepsilon)$-transform as introduced in [FGP23]. For simplicity, we
discuss the case $\bm{m}=1\\!\\!1$, the general problem following by simply
considering the modified Hamiltonian
$\operatorname{H}_{\bm{m}}^{\varepsilon}:=\operatorname{H}-\varepsilon\log\bm{m}$,
see also [FGP23, Remark 2.5]. We denote by
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ the maximising functional of
the dual problem, namely for every
$\bm{U}=(U_{1},\dots,U_{N})\in\bigtimes_{i=1}^{N}\mathcal{S}(\mathfrak{h}_{i})$,
we set

$\displaystyle\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}):=\sum_{i=1}^{N}\operatorname{Tr}(U_{i}\gamma_{i})-\varepsilon\operatorname{Tr}\bigg{(}\exp\bigg{[}\frac{\bigoplus_{i=1}^{N}U_{i}-\operatorname{H}}{\varepsilon}\bigg{]}\bigg{)}+\varepsilon\,.$
(39)

The $j$-th $(\operatorname{H},\varepsilon)$-transform of $\bm{U}$ is then the
solution to the variational problem

$\displaystyle\mathop{\mathrm{argmax}}\Big{\\{}\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(U_{1},\dots,U_{j-1},V,U_{j+1},\dots,U_{N})\
:\ V\in\mathcal{S}(\mathfrak{h}_{j})\Big{\\}}\,.$ (40)

One can show that the above maximisation problem admits a unique solution,
that we denote by
$\mathfrak{T}_{j}^{\varepsilon}\big{(}\bm{U}\big{)}\in\mathcal{S}(\mathfrak{h}_{j})$.
A direct analysis of the optimality conditions associated with the
maximisation (40) shows that
$\mathfrak{T}_{j}^{\varepsilon}\big{(}\bm{U}\big{)}$ can be characterised as
the unique solution to

$\displaystyle\operatorname{P}_{j}\left[\mathcal{T}^{\varepsilon}(\bm{U})\right]=\gamma_{j}\quad\text{for}\quad\mathcal{T}^{\varepsilon}(\bm{U}):=\exp\left(\frac{\bigoplus_{i=1}^{j-1}U_{i}\oplus\mathfrak{T}_{j}^{\varepsilon}(\bm{U})\bigoplus_{i=j+1}^{N}U_{i}-\operatorname{H}}{\varepsilon}\right)\,,$
(41)

where recall that $\operatorname{P}_{j}$ denotes the $j$-th marginal operator
(28). In other words, the density matrix
$\mathcal{T}^{\varepsilon}(\bm{U})\in\mathfrak{P}(\mathfrak{h})$ has the right
$j$-th marginal operator $\gamma_{j}$, similarly to the commutative setting.
On the other hand, a crucial difference is that, due to the possible lack of
commutativity between the involved operators, one has no way to provide an
explicit formula for $\mathfrak{T}_{i}^{\varepsilon}\big{(}\bm{U}\big{)}$.
Nonetheless, a careful analysis of the maximisation problem still provides
uniform bounds in terms of spectral bounds for the Hamiltonian
$\operatorname{H}$. More precisely, in [FGP23] the authors show that there
exists a renormalisation operator
$\operatorname{Ren}:\bigtimes_{i=1}^{N}\mathcal{S}(\mathfrak{h}_{i})\to\bigtimes_{i=1}^{N}\mathcal{S}(\mathfrak{h}_{i})$
with the following properties:

  1. 1.

$\operatorname{Ren}$ is an idempotent operator, i.e.
$\operatorname{Ren}\circ\operatorname{Ren}=\operatorname{Ren}$.

  2. 2.

$\operatorname{Ren}(\bm{U})=\bm{U}+\big{(}\alpha_{1}(\bm{U})1\\!\\!1,\dots,\alpha_{N}(\bm{U})1\\!\\!1\big{)}$,
where $\alpha_{i}(\bm{U})\in\mathbb{R}$ is so that

$\displaystyle\sum_{i=1}^{N}\alpha_{i}(\bm{U})=0\,,\quad\text{for
every}\quad\bm{U}\in\bigtimes_{i=1}^{N}\mathcal{S}(\mathfrak{h}_{i}).$ (42)

Therefore, we have that
$\bigoplus_{i}(\operatorname{Ren}(\bm{U}))_{i}=\bigoplus_{i}U_{i}$. In
particular, it does not change the value of the dual functional, namely
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\operatorname{Ren}(\bm{U}))=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U})$.

  3. 3.

Once we compose $\operatorname{Ren}$ with the
$(\operatorname{H},\varepsilon)$-transform operators, we obtain a bounded
nonlinear map from $\bigtimes_{i=1}^{N}\mathcal{S}(\mathfrak{h}_{i})$ onto
itself. More precisely, we have that

$\displaystyle\big{|}(\operatorname{Ren}\tau(\bm{U}))_{i}-\varepsilon\log\gamma_{i}\big{|}\leq
2\|{\rm H}\|_{\infty}1\\!\\!1,\quad\forall i\in[N]\,,$ (43)

where $\tau$ denotes the Sinkhorn operator, given by

$\displaystyle\begin{gathered}\tau:\bigtimes_{i=1}^{N}\mathcal{S}(\mathfrak{h}_{i})\to\bigtimes_{i=1}^{N}\mathcal{S}(\mathfrak{h}_{i}),\\\
\tau(\bm{U}):=(\mathcal{T}^{\varepsilon}_{N}\circ\dots\circ\mathcal{T}^{\varepsilon}_{1})(\bm{U})\,.\end{gathered}$
(46)

There are two important observations worth discussion at this point: the first
is that a renormalisation procedure is necessary in order to obtain uniform
bounds of the type (43). Indeed, the dual functional being invariant by
translation of operators as in (42) implies that no uniform bound is in
general expected on maximising sequences for
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$. Secondly, an interesting
feature of the bound in (43) is that it is dimension free, in the sense that
the right-hand side only depends on $\operatorname{H}$ but not on $\bm{d}$.
This possibly represents a hint to the fact that some of the estimates and
properties of the $(\operatorname{H},\varepsilon)$-transform obtained in
[FGP23] in the finite dimensional setting ($\bm{d}<\infty$) might be
extendable to an infinite dimensional setting. Of course, in this case the
estimate would be useful only for bounded operators $\operatorname{H}$. In
case of unbounded operators, a different type of analysis would be certainly
required.

We are finally ready to introduce the noncommutative version of the Sinkhorn
algorithm. Its main goal is to provide a sequence of density matrices
$\\{\Gamma^{(n)}\\}_{n}$ having one correct marginal at a time, in such a way
that $\Gamma^{(n)}\to\Gamma^{\varepsilon}$ (the optimal solution) as
$n\to\infty$. In contrast with the commutative framework, no explicit formulas
are available, hence no product form for $\Gamma^{(n)}$is ensured. On the
other hand, we use the $(\operatorname{H},\varepsilon)$-transforms in order to
ensure the correct marginals (one at a time) and the renormalisation operator
$\operatorname{Ren}$ to obtain the sought compactness to pass to the limit in
the number of iterations going to $\infty$.

Initialisation. We initialise
$\Gamma^{(0)}=\exp\left(\frac{\bigoplus_{i=1}^{N}U_{i}^{(0)}-\operatorname{H}}{\varepsilon}\right)\in\mathfrak{P}(\mathfrak{h})$
in exponential form, for a given, initial vector of self-adjoint operators
$\bm{U}^{(0)}=(U_{1}^{(0)},\dots,U_{N}^{(0)})$, with
$U_{i}^{(0)}\in\mathcal{S}(\mathfrak{h}_{i})$.

Iteration. We iteratively define
$\Gamma^{(n)}=\exp\left(\frac{\bigoplus_{i=1}^{N}U_{i}^{(n)}-H}{\varepsilon}\right)\in\mathfrak{P}(\mathfrak{h})$
by applying the $(\operatorname{H},\varepsilon)$-transforms for each component
$i\in[N]$. Precisely, we iteratively define

$\displaystyle\bm{U}^{(n)}:=\tau\big{(}\bm{U}^{(n-1)}\big{)}\,,\quad
n\in\mathbb{N}\,,$

where $\tau$ is the Sinkhorn operator introduced in (46). Taking advantage of
the uniform estimates for the renormalisation operator (43), one gets the
following convergence result.

######  Theorem 3.4.

_([FGP23, Prop. 5.1])_ The previous described algorithm converges as
$n\to\infty$, in the sense that:

  1. 1.

There exist
$\bm{U}^{\varepsilon}\in\bigtimes_{i=1}^{N}\mathcal{S}(\mathfrak{h}_{i})$ and
$\bm{\alpha}^{n}\in\mathbb{R}^{N}$ with $\sum_{i=1}^{N}\alpha_{i}^{n}=0$ such
that

$\displaystyle\bm{U}^{(n)}+\bm{\alpha}^{n}\rightarrow\bm{U}^{\varepsilon}\quad\text{as
}n\to\infty\,.$ (47)

  2. 2.

$\bm{U}^{\varepsilon}=(\bm{U}_{1}^{\varepsilon},\dots,\bm{U}_{N}^{\varepsilon})$
is optimal for the dual problem $\mathfrak{D}^{\varepsilon}(\gamma)$ as
defined in (35).

  3. 3.

$\Gamma^{(n)}$ converges as $n\to\infty$ to some
$\Gamma^{\varepsilon}\in\mathfrak{P}(\mathfrak{h})$ which is optimal for the
primal problem $\mathfrak{F}^{\varepsilon}(\bm{\gamma})$ as defined in (33).
In particular, it holds

$\displaystyle\Gamma^{\varepsilon}=\exp\left(\frac{\bigoplus_{i=1}^{N}U^{\varepsilon}_{i}-\operatorname{H}}{\varepsilon}\right)\,.$
(48)

Note that the renormalisation with $\bm{\alpha}^{n}$ leaves $\Gamma^{(n)}$
invariant, i.e.

$\displaystyle\Gamma^{(n)}=\exp\left(\frac{\bigoplus_{i=1}^{N}U^{(n)}_{i}-\operatorname{H}}{\varepsilon}\right)=\exp\left(\frac{\bigoplus_{i=1}^{N}\big{(}U^{(n)}_{i}+\alpha_{i}^{n}\big{)}-\operatorname{H}}{\varepsilon}\right)\,,$

therefore the convergence of $\bm{U}^{(n)}+\bm{\alpha}^{n}$ ensures the
convergence of $\Gamma^{(n)}$.

#### The symmetry-constrained problems: fermionic and bosonic setting

It is interesting to see that the theory presented up to this point can be
also extended in the physically relevant case where additional symmetry
constraints on the space of density matrices are added. Two typical types of
particles in physics are bosons and fermions, which from a mathematical
perspective are described by density operators that are respectively symmetric
and antisymmetric with respect to the underlined product space. Particles are
hereby considered indistinguishable, in particular the Hibert spaces
$\\{\mathfrak{h}_{i}\\}_{i}$ are to be considered the same. In the finite
dimensional setting, this can be identified with the complex space
$\mathbb{C}^{d}$ for a given dimension $d\in\mathbb{N}$.

From a mathematical perspective, we shall consider the symmetric and
antisymmetric tensor product of the one-particle Hilbert space (instead of the
free tensor product in the setting without symmetries), namely we work with

$\displaystyle\left(\bigodot_{i=1}^{N}\mathbb{C}^{d}\right)\,,\,\left(\bigwedge_{i=1}^{N}\mathbb{C}^{d}\right)\subset\bigotimes_{i=1}^{N}\mathbb{C}^{d}\,,$
(49)

where $\odot$ (resp. $\wedge$) denotes the symmetric (resp. antisymmetric)
tensor product. Note that the cardinality of
$\bigwedge_{i=1}^{N}\mathbb{C}^{d}$ is $\binom{d}{N}$ if $d\geq N$, and
$\bigwedge_{i=1}^{N}\mathbb{C}^{d}=\\{0\\}$ if $d<N$, due to the fact that
$e\wedge e=0$. In particular, in the fermionic setting, we assume to work with
a number of particles $N$ which is smaller than the considered dimension $d$.
We denote by

$\displaystyle\mathfrak{P}_{+}^{\bm{d}}:=\mathfrak{P}\left(\bigodot_{i=1}^{N}\mathbb{C}^{d}\right),\,\
\mathfrak{P}_{-}^{\bm{d}}:=\mathfrak{P}\left(\bigwedge_{i=1}^{N}\mathbb{C}^{d}\right)\,\
\
\subset\,\mathfrak{P}^{\bm{d}}:=\mathfrak{P}\left(\bigotimes_{i=1}^{N}\mathbb{C}^{d}\right)\,,$
(50)

the set of bosonic and fermionic density matrices. To simplify the notation,
we also write $\mathcal{M}^{d}:=\mathcal{M}(\mathbb{C}^{d})$,
$\mathcal{S}^{d}:=\mathcal{S}(\mathbb{C}^{d})$. It is key to observe that, in
both cases, the marginals operators of a given density matrix
$\Gamma\in\mathfrak{P}_{\pm}^{d}$ all coincide.

The fermionic symmetry constraint has an important consequence on the type of
operators that one can obtain as marginals, namely not every operator
$\gamma\in\mathfrak{P}^{d}$ can be obtained as marginal of a fermionic
operator. It is indeed well-established that the existence of a
$\Gamma\in\mathfrak{P}_{-}^{\bm{d}}$ so that $\Gamma\mapsto\gamma$ (where
$\Gamma\mapsto\gamma$ means that $\Gamma$ has all marginals equal to $\gamma$)
is equivalent to $\gamma$ satisfying the so-called Pauli exclusion principle,
i.e. if and only if $\gamma\in\mathfrak{P}^{d}$ and $\gamma\leq 1/N$ (see for
example [LS10, Theorem 3.2]). As we are going to see in the final part of
these notes, the dual formulation of the entropic transport problem provides,
in a natural way, a more variational interpretation of the Pauli principle.

In order to talk about optimal transport, we must fix a reference Hamiltonian.
In the symmetric setting, it is natural to consider an operator which also
satisfies some additional symmetry conditions. To this purpose, for $i\in[N]$,
we introduce the permutation operators
$\operatorname{S}_{i}:\mathcal{M}^{\bm{d}}\simeq\bigotimes_{j=1}^{N}\mathcal{M}^{d_{j}}\to\mathcal{M}^{\bm{d}}$
as the map defined by

$\displaystyle\operatorname{S}_{i}\bigg{(}\bigotimes_{j=1}^{N}A_{j}\bigg{)}=A_{1}\otimes\dots\otimes
A_{i-1}\otimes A_{i+1}\otimes\dots\otimes A_{N}\otimes A_{i},$

for any $A_{j}\in\mathcal{M}^{d_{j}}$ and extended to the whole
$\mathcal{M}^{\bm{d}}$ by linearity. The natural assumption on
$\operatorname{H}\in\mathcal{S}^{\bm{d}}$ is then to assume that

$\displaystyle\operatorname{S}_{i}\circ\operatorname{H}\circ\operatorname{S}_{i}=\operatorname{H}\,,\quad\forall
i=1,\dots,N\,.$ (51)

An interpretation of the latter condition is that the quantum systems is
driven by an Hamiltonian which does not distinguish between the different $N$
particles.

We are now in position to define the bosonic and the fermionic entropy-
regularised optimal transport problem. For any $\gamma\in\mathfrak{P}^{d}$, we
define the bosonic primal problem as

$\displaystyle\mathfrak{F}^{\varepsilon}_{+}(\gamma):=\inf\left\\{\operatorname{Tr}(\operatorname{H}\Gamma)+\varepsilon\operatorname{Tr}(\Gamma\log\Gamma)\
:\ \Gamma\in\mathfrak{P}_{+}^{\bm{d}}\text{ and
}\Gamma\mapsto\gamma\right\\}\,.$ (52)

In light of the Pauli principle for fermionic density operators, we define the
fermionic primal problem for a given $\gamma\in\mathfrak{P}^{d}$ satisfying
$\gamma\leq 1/N$ as

$\displaystyle\mathfrak{F}^{\varepsilon}_{-}(\gamma):=\inf\left\\{\operatorname{Tr}(\operatorname{H}\Gamma)+\varepsilon\operatorname{Tr}(\Gamma\log\Gamma)\
:\ \Gamma\in\mathfrak{P}_{-}^{\bm{d}}\text{ and
}\Gamma\mapsto\gamma\right\\}\,.$ (53)

In order to introduce and discuss the dual formulation of the
bosonic/fermionic primal problems, we first introduce two natural projection
operators $\Pi_{+}$ (resp. $\Pi_{-}$), given by

$\displaystyle\Pi_{+}:\bigotimes_{i=1}^{N}\mathbb{C}^{d}\to\bigodot_{i=1}^{N}\mathbb{C}^{d}\,,\qquad\Pi_{-}:\bigotimes_{i=1}^{N}\mathbb{C}^{d}\to\bigwedge_{i=1}^{N}\mathbb{C}^{d}\,.$
(54)

For any given operator $A\in\mathcal{S}^{\bm{d}}$, we denote by $A_{\pm}$ the
corresponding projection onto the symmetric space, obtained as
$A_{\pm}:=\Pi_{\pm}\circ A\circ\Pi_{\pm}$ . Consequently, for any
$\gamma\in\mathfrak{P}^{d}$, we define the bosonic dual functional
$\text{D}_{\gamma}^{+,\varepsilon}$ and the fermionic dual functional
$\text{D}_{\gamma}^{-,\varepsilon}$ as

$\displaystyle\text{D}_{\gamma}^{\pm,\varepsilon}:\mathcal{S}^{d}\to\mathbb{R}\,,\quad\text{D}_{\gamma}^{\pm,\varepsilon}(U):=\operatorname{Tr}(U\gamma)-\varepsilon\operatorname{Tr}\left(\exp\bigg{[}\frac{1}{\varepsilon}\bigg{(}\frac{1}{N}\bigoplus_{i=1}^{N}U-H\bigg{)}_{\pm}\bigg{]}\right)+\varepsilon\,.$
(55)

The corresponding dual problems are given by

$\displaystyle\mathfrak{D}_{\pm}^{\varepsilon}(\gamma):=\sup\left\\{\text{D}_{\gamma}^{\pm,\varepsilon}(U)\
:\ U\in\mathcal{S}^{d}\right\\}\,.$ (56)

Note that a priori $\mathfrak{D}^{\varepsilon}_{-}(\gamma)$ is well-defined
for $\gamma\in\mathfrak{P}^{d}$ without any additional assumption, whereas
$\mathfrak{F}^{\varepsilon}_{-}(\gamma)$ is only well defined for
$\gamma\in\mathfrak{P}^{d}$ such that $\gamma\leq 1/N$, due to the presence of
the Pauli principle. Nonetheless, once we try to maximise the dual functional,
the constraint provided by the Pauli principle naturally arises as an
admissibility condition in order to have
$\mathfrak{D}^{\varepsilon}_{-}(\gamma)<\infty$. Another interesting feature
concerns the existence of maximisers: to ensure the existence of a maximiser
for $\text{D}_{\gamma}^{-,\varepsilon}$, one has to further impose that
$\gamma<1/N$, a condition that the authors in [FGP23] call strict Pauli
principle. This variational interpretation of the principle based on the dual
formulation is the content of the next theorem [FGP23, Prop. 2.8].

######  Proposition 3.5 (Pauli’s principle and duality).

Let $N\geq d$. The following are equivalent:

  1. 1.

$\mathfrak{D}^{\varepsilon}_{-}(\gamma)<+\infty$ if and only if
$\gamma\in\mathfrak{P}^{d}$ and $\gamma\leq\frac{1}{N}$.

  2. 2.

There exists a maximizer $U_{0}\in\mathcal{S}^{d}$ of
$\emph{D}_{\gamma}^{-,\varepsilon}$ if and only if $\gamma\in\mathfrak{P}^{d}$
and $0<\gamma<\frac{1}{N}$.

The proof of the necessity of the Pauli principle to have the dual functional
bounded from above provided in [FGP23] is constructive. In particular,
reasoning by contradiction, one can explicitly construct a sequence of
operators $U^{n}\in\mathcal{S}^{d}$ along which the functional
$\text{D}_{\gamma}^{-,\varepsilon}$ diverges to $+\infty$. To see this,
suppose by contradiction that the Pauli principle is not satisfied. By
spectral calculus, and up to reordering the terms in the sum, we can without
loss of generality assume that

$\displaystyle\gamma=\sum_{i=1}^{d}\gamma_{i}|\psi_{i}\rangle\langle\psi_{i}|,\quad\gamma_{1}>\frac{1}{N}\,.$

The sought sequence of operators $U^{n}\in\mathcal{S}^{d}$ is then obtained as
matrices having the same eigenspaces of $\gamma$, choosing a very positive
eigenvalue in correspondence to first eigenspace (where the condition
described by the Pauli principle does not hold), and negative eigenvalues
otherwise. Precisely, we define

$\displaystyle
U^{n}:=\sum_{i=1}^{d}u_{i}^{n}|\psi_{i}\rangle\langle\psi_{i}|,\quad
u_{1}^{n}:=n\,,\quad u_{j}^{n}:=-\frac{n}{N-1},\quad\forall j\geq 2\,.$ (57)

This sequence is constructed in such a way that, for any choice of indexes
$\bm{j}=(j_{1},\dots,j_{N})$ with $j_{i}\in\\{1,\dots,d\\}$ without repetition
(namely $j_{i}\neq j_{k}$ if $i\neq k$), we have that

$\displaystyle\exp\left(\frac{1}{N}\sum_{i=1}^{N}u_{j_{i}}\right)\begin{cases}=1&\text{if
}j_{i}=1\text{ for some }i\,,\\\ \leq 1&\text{otherwise}\,.\end{cases}$

This ensures that the nonlinear part of $\text{D}_{\gamma}^{-,\varepsilon}(U)$
can be estimated from above by the dimension of the antisymmetric product
space. In particular, we have that

$\displaystyle\text{D}_{\gamma}^{-,\varepsilon}(U^{n})\geq\sum_{j=1}^{d}\gamma_{j}u_{j}^{n}-C\binom{d}{N}\,.$
(58)

The claim is that the linear part in the bound above goes to $+\infty$ as
$n\to+\infty$. Indeed

$\displaystyle\sum_{j=1}^{d}\gamma_{j}u_{j}^{n}=n\bigg{(}\gamma_{1}-\frac{1}{N-1}\sum_{i=2}^{d}\gamma_{i}\bigg{)}=\frac{n}{N-1}\Big{(}N\gamma_{1}-1\Big{)}\,,$
(59)

where we used that $\sum_{i}\gamma_{i}=1$. From this, using
$\gamma_{1}>\frac{1}{N}$ and (58), we deduce
$\text{D}_{\gamma}^{-,\varepsilon}(U^{n})\to+\infty$ as $n\to+\infty$, which
provides the sought contradiction.

Finally, putting together the duality result in the non-symmetric case and the
relation to the Pauli principle, one can prove a duality result also in the
bosonic and fermionic setting, see [FGP23, Thm 2.9].

######  Theorem 3.6 (Fermionic and bosonic duality).

Let $\operatorname{H}\in\mathcal{S}^{\bm{d}}$ satisfy the condition (51).

  * (i)

Assume that $N\geq d$. For any $\gamma\in\mathfrak{P}^{d}$ such that
$\gamma\leq\frac{1}{N}$, the fermionic primal and dual problems coincide, i.e.
$\mathfrak{F}^{\varepsilon}_{-}(\gamma)=\mathfrak{D}^{\varepsilon}_{-}(\gamma)$.
If additionally $\gamma$ satisfies the strict Pauli principle
$0<\gamma<\frac{1}{N}$, then $\emph{D}_{\gamma}^{-,\varepsilon}$ admits a
unique maximiser $U_{-}^{\varepsilon}\in\mathcal{S}^{d}$ and

$\displaystyle\Gamma_{-}^{\varepsilon}=\exp\left(\frac{1}{\varepsilon}\bigg{[}\frac{1}{N}\bigoplus_{i=1}^{N}U_{-}^{\varepsilon}-\operatorname{H}\bigg{]}_{-}\right)$
(60)

is the unique optimal fermionic solution to the primal problem
$\mathfrak{F}^{\varepsilon}_{-}(\gamma)$.

  * (ii)

Fix $N\in\mathbb{N}$. For any $\gamma\in\mathfrak{P}^{d}$, the bosonic primal
and dual problems coincide, i.e.
$\mathfrak{F}^{\varepsilon}_{+}(\gamma)=\mathfrak{D}^{\varepsilon}_{+}(\gamma)$.
Moreover, if $\gamma>0$ (or equivalently $\ker\gamma=\\{0\\}$), then
$\emph{D}_{\gamma}^{+,\varepsilon}$ admits a unique maximiser
$U_{+}^{\varepsilon}$ and

$\displaystyle\Gamma_{+}^{\varepsilon}=\exp\left(\frac{1}{\varepsilon}\bigg{[}\frac{1}{N}\bigoplus_{i=1}^{N}U_{+}^{\varepsilon}-\operatorname{H}\bigg{]}_{+}\right)$
(61)

is the unique minimiser in the primal problem
$\mathfrak{F}^{\varepsilon}_{+}(\gamma)$.

### Acknowledgments

We are very thankful to the organisers of the School and Workshop on Optimal
Transport on Quantum Structures at Erdös Center: J. Maas, S. Rademacher, T.
Titkos, D. Virosztek, for the nice scientific experience. We are grateful to
the anonymous reviewer for their careful reading and useful suggestions. The
author also acknowledges fundings from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy \- GZ
2047/1, Projekt-ID 390685813.

## References

  * [CGP16] Yongxin Chen, Tryphon Georgiou, and Michele Pavon.  Entropic and displacement interpolation: a computational approach using the Hilbert metric.  SIAM Journal on Applied Mathematics, 76(6):2375–2396, 2016.
  * [CGP18] Emanuele Caglioti, François Golse, and Thierry Paul.  Towards optimal transport for quantum densities.  Available on https://hal.archives-ouvertes.fr/hal-01963667/, 2018\.
  * [CGT18] Y. Chen, T. T. Georgiou, and A. Tannenbaum.  Matrix Optimal Mass Transport: A Quantum Mechanical Approach.  IEEE Transactions on Automatic Control, 63(8):2612–2619, 2018.
  * [CM14] Eric A Carlen and Jan Maas.  An analog of the 2-Wasserstein metric in non-commutative probability under which the Fermionic Fokker–Planck equation is gradient flow for the entropy.  Communications in mathematical physics, 331(3):887–926, 2014.
  * [Cut13] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In Advances in Neural Information Processing Systems, pages 2292–2300, 2013.
  * [DMG20] Simone Di Marino and Augusto Gerolin.  An Optimal Transport approach for the Schrödinger bridge problem and convergence of Sinkhorn algorithm.  Journal of Scientific Computing, 85(27), 2020.
  * [DPT21] Giacomo De Palma and Dario Trevisan.  Quantum optimal transport with quantum channels.  Ann. Henri Poincaré, 22(10):3199–3234, 2021.
  * [FGP23] Dario Feliciangeli, Augusto Gerolin, and Lorenzo Portinale.  A non-commutative entropic optimal transport approach to quantum composite systems at positive temperature.  J. Funct. Anal., 285(4):Paper No. 109963, 39, 2023.
  * [FL89] Joel Franklin and Jens Lorenz.  On the scaling of multidimensional matrices.  Linear Algebra and its applications, 114:717–735, 1989.
  * [Gil75] Thomas L Gilbert.  Hohenberg-Kohn theorem for nonlocal external potentials.  Physical Review B, 12(6):2111, 1975.
  * [GP15] Tryphon T Georgiou and Michele Pavon.  Positive contraction mappings for classical and quantum Schrödinger systems.  Journal of Mathematical Physics, 56(3):033301, 2015.
  * [GS98] Wilfrid Gangbo and Andrzej Świȩch.  Optimal maps for the multidimensional Monge-Kantorovich problem.  Comm. Pure Appl. Math., 51(1):23–45, 1998.
  * [GS10] Alfred Galichon and Bernard Salanié.  Matching with trade-offs: Revealed preferences over competing characteristics.  2010\.
  * [HK64] Pierre Hohenberg and Walter Kohn.  Inhomogeneous electron gas.  Phys. Rev., 136(B864), 1964.
  * [Léo14] Christian Léonard.  A survey of the Schrödinger problem and some of its connections with optimal transport.  Discrete & Continuous Dynamical Systems-A, 34(4):1533–1574, 2014\.
  * [Lev76] Mel Levy.  Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem.  Proceedings of the National Academy of Sciences, 12:6062–6065, 76\.
  * [Lie02] Elliott H Lieb.  Density functionals for Coulomb systems.  In Inequalities, pages 269–303. Springer, 2002.
  * [LS10] Elliott H. Lieb and Robert Seiringer.  The stability of matter in quantum mechanics.  Cambridge University Press, Cambridge, 2010.
  * [MM17] Markus Mittnenzweig and Alexander Mielke.  An entropic gradient structure for Lindblad equations and couplings of quantum systems to macroscopic models.  Journal of Statistical Physics, 167(2):205–233, 2017.
  * [Mon81] Gaspard Monge.  Mémoire sur la théorie des déblais et des remblais.  De l’Imprimerie Royale, 1781.
  * [MV21] Léonard Monsaingeon and Dmitry Vorotnikov.  The Schrödinger problem on the non-commutative Fisher-Rao space.  Calculus of Variations and Partial Differential Equations, 60(1):1–48, 2021.
  * [PCVS19] Gabriel Peyré, Lénaïc Chizat, François-Xavier Vialard, and Justin Solomon.  Quantum entropic regularization of matrix-valued optimal transport.  European Journal of Applied Mathematics, 30(6):1079–1102, 2019\.
  * [Rus95] Ludger Ruschendorf.  Convergence of the iterative proportional fitting procedure.  The Annals of Statistics, pages 1160–1174, 1995.
  * [Sin64] Richard Sinkhorn.  A relationship between arbitrary positive matrices and doubly stochastic matrices.  The annals of mathematical statistics, 35(2):876–879, 1964.
  * [Vil08] Cédric Villani.  Optimal transport: old and new, volume 338.  Springer Science & Business Media, 2008.
  * [Wir22] Melchior Wirth.  A dual formula for the noncommutative transport distance.  J. Stat. Phys., 187(2):Paper No. 19, 18, 2022.


# Relativistic motion under constant force: velocity and acceleration behavior

Jhonnatan G. Pereira, Victor E. Mouchrek-Santos, Manoel M. Ferreira Jr
Departamento de Física, Universidade Federal do Maranhão (UFMA), Campus
Universitário do Bacanga, São Luís-MA, 65085-580 - Brasil
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In relativistic dynamics, force and acceleration are no longer parallel. In
this article, we revisit the relativistic motion of a particle under the
action of a constant force, $\boldsymbol{f}$. For a two-dimensional motion,
the final velocity in each axis is $(f_{i}/f)c,$ independently of the initial
velocities, yielding an asymptotic velocity always parallel to the force. The
particular case in which the force is applied in a single axis is analyzed in
detail, with the behavior of velocity and acceleration being exhibited for
several configurations. Some previous results of the literature concerning
velocity and acceleration behavior are improved and better explored.
Differently from which was previously claimed, it is shown that a negative
acceleration component can exist in the direction of the biggest force
component and that acceleration does not decrease monotonically to zero.

##  I Introduction

In relativistic dynamics Rindler ,Aharoni ,French ,Hagedorn , the usual
concept of mass needs to be changed in order to preserve conservation of
linear momentum, in such a way it becomes a function of velocity,
$m(\mathrm{v})=m_{0}\gamma(\mathrm{v}),$ with
$\gamma(\mathrm{v})=1/\sqrt{1-\mathbf{v}^{2}/c^{2}}$. The 3-momentum and
energy are accordingly redefined, $\boldsymbol{p}=m(\mathrm{v})\mathbf{v}$,
$E=m_{0}\gamma(\mathrm{v})c^{2},$ with the 3-force given as
$\boldsymbol{f}=d\boldsymbol{p}/dt$, or $\
\boldsymbol{f}=m\boldsymbol{a}+\dot{m}\mathbf{v}.$ Relativistic dynamics has
been studied in several interesting respects, involving methods for deriving
force and momentum transformations Martins ,Stewart ,Momentum , relativistic
mass Mass , and the description of relativistic collisions Collision1 ,
Collision2 .

The mass dependence on the velocity yields
$\dot{m}=m_{0}\dot{\gamma}=m\gamma^{2}(\mathbf{v}\cdot\boldsymbol{a})\boldsymbol{/}\boldsymbol{c}^{2}$
($\dot{m}=dm/dt)$, so that the relativistic 3-force presents two components,

$\boldsymbol{f}=m\boldsymbol{a}+m\gamma^{2}(\mathbf{v}\cdot\boldsymbol{a})\mathbf{v}\boldsymbol{/}c^{2},$
(1)

one in the direction of acceleration, other parallel to velocity, becoming the
force, in general, no longer parallel to the acceleration engendered by
itself. This is a relativistic feature with no analogue in Newtonian dynamics,
opening a diversity of dynamic scenarios very different from the usual ones.
Taking the internal product of the force $\boldsymbol{f}$ with $\mathbf{v},$
in eq. (1), it results
$(\boldsymbol{f}\cdot\mathbf{v)}=m\gamma^{2}(\mathbf{v}\cdot\boldsymbol{a})$.
Thus, eq. (1) can be rewritten as:

$\boldsymbol{f}=m\boldsymbol{a}+(\boldsymbol{f}\cdot\mathbf{v})\mathbf{v}\boldsymbol{/}c^{2}.$
(2)

Taking the internal product of eq. (2) with force, it implies
$\boldsymbol{f}^{2}-m(\boldsymbol{a\cdot
f)}=(\boldsymbol{f}\cdot\mathbf{v})^{2}\boldsymbol{/}c^{2}.$ Since
$(\boldsymbol{f}\cdot\mathbf{v})^{2}\boldsymbol{/}c^{2}\leq\boldsymbol{f}^{2},$
one has $(\boldsymbol{a\cdot f)}>0,$ revealing that $\boldsymbol{a}$ and
$\boldsymbol{f}$ define an acute angle, as demonstrated in Ref. Redding .
Another contribution investigated the orientation of the acceleration
components relative to the force and its dependence on the system of reference
Lock .

One of the first studies examining the intricate relation between force and
acceleration in relativistic dynamics is due to Tolman Tolman , which
investigated the condition to have acceleration only in one axis. Tolman
noticed that applying force in one direction causes acceleration in this axis
and mass increases. Since the 3-momentum is conserved in the direction where
there is no applied force, the velocity in that direction should decrease in
order to assure a constant momentum. The x and y components of eq. (1) are

$\displaystyle f_{x}$
$\displaystyle=ma_{x}+m\gamma^{2}(\mathrm{v}_{x}a_{x}+\mathrm{v}_{y}a_{y})\mathrm{v}_{x}\boldsymbol{/}c^{2},$
(3) $\displaystyle f_{y}$
$\displaystyle=ma_{y}+m\gamma^{2}(\mathrm{v}_{x}a_{x}+\mathrm{v}_{y}a_{y})\mathrm{v}_{y}\boldsymbol{/}c^{2}.$
(4)

Supposing an applied force in two dimensions,
$\boldsymbol{f}=(f\boldsymbol{{}_{x}},f\boldsymbol{{}_{y}}),$ for acceleration
only in the y-axis, $a_{x}=0,$ Tolman wrote the ratio,

$\frac{f_{x}}{f_{y}}=\frac{\gamma^{2}\mathrm{v}_{y}\mathrm{v}_{x}\boldsymbol{/}c^{2}}{1+\gamma^{2}\mathrm{v}_{y}^{2}\boldsymbol{/}c^{2}}=\frac{\mathrm{v}_{y}\mathrm{v}_{x}}{c^{2}-\mathrm{v}_{x}^{2}},$
(5)

which states the relation to be fulfilled to guarantee $a_{x}=0,$ involving
force applied in both axis. Tolman also observed Tolman2  the existence of two
situations where force and acceleration remain parallel: (i) when force is
applied longitudinally to the direction of the motion
($\boldsymbol{f}\parallel\mathbf{v}$), (ii) when force is applied in a
direction orthogonal to the velocity of the particle, making it describe a
circular motion ($\boldsymbol{f}\bot\mathbf{v}$). In both situations, the
initial geometrical configuration is maintained in time: the force is kept
parallel to acceleration as the motion evolves. For force parallel to
velocity, eq. (2) provides $\boldsymbol{f}=m_{\parallel}\boldsymbol{a}$, with
$m_{\parallel}=m_{0}\gamma^{3}$ being the longitudinal mass; for force
orthogonal to velocity, it holds $\boldsymbol{f}=m_{\bot}a$, where
$m_{\bot}=m_{0}\gamma$ is orthogonal mass, with ($m_{\bot}<m_{\parallel}$).

In most relativistic scenarios, force and velocity are not parallel. The
richness of these situations was already addressed in some respects Redding .
One of the remarkable consequences is the existence of a negative acceleration
component Lock ,Ficken , in two-dimensional systems with force applied in both
orthogonal axis, $f_{x}>0,$ $f_{y}>0.$ This motion was examined by plotting
the value of the ratio $a_{y}/a_{x\text{ }}$for several values of velocity
($\mathrm{v}_{x\text{ }}$and $\mathrm{v}_{y})$ in the interval
$\left[0,\sqrt{c^{2}/2}\right],$ for different values of the ratio
$R=f_{x}/f_{y}$ Ficken . A similar analysis was developed in three dimensions
Gonzalez , so that the ratios $a_{x}/a_{y}$ and $a_{z}/a_{y}$ were plotted for
values of velocities $\mathrm{v}_{x}$, $\mathrm{v}_{y}$ and $\mathrm{v}_{z}$
in the interval $[0,\sqrt{c^{2}/3}]$, for different values of the ratios
$R_{1}=f_{x}/f_{y}$, $R_{2}=f_{y}/f_{z}$. These authors asserted that the
negative acceleration appears in the direction where the applied force is
smaller, which is not always true, however. This kind of analysis, based on
the ratio $R=f_{x}/f_{y},$ is restricted in some respects, namely: (i) the
components of velocity are made equal ($\mathrm{v}_{x}=\mathrm{v}_{y})$ all
the time in order to simplify the analysis, which exclude interesting
situations where the velocities differ from each other; (ii) the behavior of
the velocity components is not properly discussed; (iii) by studying the ratio
between accelerations, it is not clear what happens with each acceleration
component.

In this work, we reassess the description of a two dimensional relativistic
motion under the action of a constant force, combining two different
approaches, an analytical and numerical procedure, in order to shed light on
the three main aspects highlighted above, which remain barely discussed until
this moment.

##  II Two-dimensional Motion of a particle under the action of a relativistic
3-force

In this section, we develop an analysis that goes beyond the aspects
scrutinized in the previous works about force and acceleration in special
relativity Ficken , Gonzalez . Remembering that $\
\boldsymbol{f}=(\boldsymbol{f}_{x},\boldsymbol{f}_{y})=(d\boldsymbol{p}_{x}/dt,d\boldsymbol{p}_{y}/dt),$
and regarding the case the force is constant, one achieves

$\displaystyle p_{x}(t)$ $\displaystyle=p_{0x}+f_{x}t,$ (6) $\displaystyle
p_{y}(t)$ $\displaystyle=p_{0y}+f_{y}t,$ (7)

in which $p_{0x}$ and $p_{0y}$ are the initial x and y components of momentum,
respectively. Knowing that the particle velocity is
$\mathbf{v}=c^{2}\mathbf{p/E,}$ where $\mathbf{p}$ is the 3-momentum and
$E=\sqrt{m_{0}^{2}c^{4}+\mathbf{p}^{2}c^{2}}$ is the relativistic energy, with
$\mathbf{p}^{2}=p_{x}^{2}+p_{y}^{2}$, the velocity components are given as

$\mathrm{v}_{i}(t)=\frac{(p_{0i}+f_{i}t)c^{2}}{\sqrt{m_{0}^{2}c^{4}+\mathbf{p}^{2}c^{2}}},$
(8)

where $i=1,2$ for $x$ and $y$ components. The direction of the force is given
by the angle $\theta$ it forms with the $x$-axis, where
$\tan\theta=f_{y}/f_{x}$. In the same way, the velocity of the particle points
in the direction of the $\alpha$ angle, written as

$\tan\alpha=\frac{\mathrm{v}_{y}}{\mathrm{v}_{x}}=\frac{p_{0y}+f_{y}t}{p_{0x}+f_{x}t}.$
(9)

In the asymptotic limit ($t\rightarrow\infty)$, eq. (9) tends to
$f_{y}/f_{x}$, therefore $\alpha\rightarrow\theta$: the direction of the
velocity approaches asymptotically the direction of the applied force. It is
thus demonstrated that the velocity will always tend to be parallel to force
after a long enough time. If the same limit $t\rightarrow\infty$ is applied to
the velocity components (8), we attain the following asymptotic values:

$\mathrm{v}_{x}(\infty)=\frac{f_{x}c}{f},\text{
}\mathrm{v}_{y}(\infty)=\frac{f_{y}c}{f},$ (10)

where $f=\sqrt{f_{x}^{2}+f_{y}^{2}}$ is the 3-force magnitude. Eq. (10) shows
that the final value of the velocity is always $c$, that is,
$\mathrm{v}(\infty)=\sqrt{\mathrm{v}_{x}^{2}(\infty)+\mathrm{v}_{y}^{2}(\infty)}=c.$
Further, the final velocity in each axis, $(f_{i}/f)c,$ depends on the value
of the corresponding component of force, independently of the initial
velocities. If there is force applied only in one axis, the final velocity in
this axis will be equal to $c$, being null in the other axis. The cases where
there is force in only one axis are well appropriate to demonstrate the
manifestation of the negative acceleration component, and will be discussed in
the next sections.

###  II.1 Velocity behavior

A particle of rest mass $m_{0},$ with initial velocity
$\mathbf{v}_{0}=(\mathrm{v}_{0x},\mathrm{v}_{0y})$, when accelerated under the
action of a constant 3-force $f=(f_{x},f_{y}),$ reaches asymptotic velocity
components equal to $(f_{i}/f)c$ with $f_{i}$ being the force component in the
$i$-axis. We know the initial and final velocities, but some aspects of the
motion evolution need to be better elucidated and discussed. As observed by
Tolman Tolman ,Tolman2 , if a force is applied in a single direction, there
will appear a (negative) acceleration in the orthogonal axis. The situation
($f_{x}>0,f_{y}=0)$, corresponding to $R=\infty,$ was in principle depicted in
the graph of Ref. Ficken , showing the existence of negative acceleration in
the y-axis. No additional details were explored, though, and are the goal of
the present work.

Consider Eq. (2), written for the acceleration components:

$\displaystyle a_{x}$
$\displaystyle=\frac{1}{m}\left[f_{x}(1-\mathrm{v}_{x}^{2}/c^{2})-f_{y}\mathrm{v}_{x}\mathrm{v}_{y}/c^{2}\right],$
(11) $\displaystyle a_{y}$
$\displaystyle=\frac{1}{m}\left[f_{y}(1-\mathrm{v}_{y}^{2}/c^{2})-f_{x}\mathrm{v}_{x}\mathrm{v}_{y}/c^{2}\right].$
(12)

Now, we specify the case where force is applied only along the y-axis,
$\boldsymbol{f}=(0,f\boldsymbol{{}_{y}}),$ with $\boldsymbol{f_{y}}>0$,
because this choice allows us to notice clearly the effects of negative
acceleration. Eq. (10) shows that the particle will be progressively
accelerated in the y-axis, having asymptotic velocity
$\mathrm{v}_{y}(\infty)=c$, while its component $\mathrm{v}_{x\text{
}}$decreases, tending to zero. For this situation Eqs. (11) and (12) simplify
to:

$\displaystyle\frac{d\mathrm{v}_{x}}{dt}$
$\displaystyle=-\frac{f_{y}\mathrm{v}_{x}\mathrm{v}_{y}}{mc^{2}},$ (13)
$\displaystyle\frac{d\mathrm{v}_{y}}{dt}$
$\displaystyle=\frac{f_{y}(1-\mathrm{v}_{y}^{2}/c^{2})}{m},$ (14)

wherein it was replaced $a_{x}=d\mathrm{v}_{x}/dt,$
$a_{y}=d\mathrm{v}_{y}/dt$. These equations will be used to reveal some
aspects not yet explored of this two-dimensional motion.

From the analysis of eqs. (13) and (14), we see that the y-acceleration is
always positive regardless of the velocity. The x-acceleration, however,
behaves differently, depending on the product $\mathrm{v}_{x}\mathrm{v}_{y}$.
If $\mathrm{v}_{x}(0)>0$ and $\mathrm{v}_{y}(0)=0$ the component $a_{x}$ is
initially null and negative afterwards. This makes the velocity
$\mathrm{v}_{x}$ to be initially constant and decreasing later. This initial
behavior of velocity is demonstrated by the little plateau at $t=0$, see Fig.
1.

Figure 1: Graph of velocities $\mathrm{v}_{x}$, $\mathrm{v}_{y}$ $\times$
time. The red and blue line represent $\mathrm{v}_{y}$ and $\mathrm{v}_{x}$,
respectively. The initial velocities are $\mathrm{v}_{x}(0)=0.8c$ and
$\mathrm{v}_{y}(0)=0$.

If $\mathrm{v}_{x}(0)>0$ and $\mathrm{v}_{y}(0)>0$, the only difference in
relation to the first case is that $\mathrm{v}_{x}$ decreases immediately from
the start, as displayed in Fig. 2. When $\mathrm{v}_{x}(0)<0$,
$\mathrm{v}_{y}(0)>0$, Eq. (13) shows that $a_{x}$ is positive. Since
$\mathrm{v}_{x}$ is negative, the particle will decelerate in the x-axis,
tending continuously to zero, as shown in Fig. 3.

Figure 2: Graph of velocities $\mathrm{v}_{x}$, $\mathrm{v}_{y}\times$ time.
The red and blue line represent $\mathrm{v}_{y}$ and $\mathrm{v}_{x}$,
respectively. The initial velocities are $\mathrm{v}_{x}(0)=0.8c$ and
$\mathrm{v}_{y}(0)=0.5c$. Figure 3: Graph of velocities $\mathrm{v}_{x}$,
$\mathrm{v}_{y}$ $\times$ time. The red and blue line represent
$\mathrm{v}_{y}$ and $\mathrm{v}_{x}$, respectively. The initial velocities
are $\mathrm{v}_{x}(0)=-0.8c$ and $\mathrm{v}_{y}(0)=0.5c$. Figure 4: Graph of
velocities $\mathrm{v}_{x}$, $\mathrm{v}_{y}$ x time. The red line represents
$\mathrm{v}_{y}$, while the blue one represents $\mathrm{v}_{x}$. The initial
velocities are $\mathrm{v}_{x}(0)=0.8c$ and $\mathrm{v}_{y}(0)=-0.5c$. Figure
5: Graph of velocities $\mathrm{v}_{x}$, $\mathrm{v}_{y}$ x time. The red and
blue line represent $\mathrm{v}_{y}$ and $\mathrm{v}_{x}$, respectively. The
initial velocities are $\mathrm{v}_{x}(0)=-0.8c$ and
$\mathrm{v}_{y}(0)=-0.5c$. Figure 6: Graph of $|\mathbf{v}|$ against time.
This graph corresponds to the situation depicted in Fig. 1. Figure 7: Graph of
$(\boldsymbol{a}\cdot\mathbf{v})$ against time. This graph corresponds to the
situation depicted in Fig. 1.

Another interesting situation is the one where $\mathrm{v}_{y}(0)<0$ and
$\mathrm{v}_{x}(0)>0$. In that case, we see from eq. (13) that $a_{x}$ is
initially positive and therefore the particle initially accelerates in the
x-axis, as shown in Fig. 4. With the force acting on the positive sense of
y-axis, it will reverse the direction of motion in the y-axis:
$\mathrm{v}_{y}$ becomes positive and the particle starts to decelerate in the
x-axis, similarly to the other cases previously discussed. For the case where
both initial velocities are negative, $\mathrm{v}_{y}(0)<0$ and
$\mathrm{v}_{x}(0)<0$, $a_{x}$ is initially negative. The particle will first
accelerate in x and then decelerate, but with $\mathrm{v}_{x}$ being negative,
as illustrated in Fig. 5. These last two cases exhibit an initial acceleration
regime followed by a deceleration one, behavior that was not reported in the
literature yet.

Figure 8: Graph of $|\mathbf{v}|$ against time. This graph corresponds to the
situation depicted in Figs. 4 and 5. Figure 9: Graph of
$\boldsymbol{a}\cdot\mathbf{v}$ against time. This graph corresponds to the
situation depicted in Figs. 4 and 5.

As shown by Redding Redding , the expression
$m_{0}\gamma^{3}(\mathbf{v}\cdot\mathbf{a})=(\boldsymbol{f}\cdot\mathbf{v})$
states the angle between force and velocity is of the same type as the one
between acceleration and velocity. By analyzing the situation depicted in Fig.
1, we see that acceleration and velocity are initially perpendicular. After
the initial moment, the particle presents a non null $\mathrm{v}_{y}$ and the
angle between force and velocity becomes acute. The same holds for the angle
between acceleration and velocity, in such a way that
$(\mathbf{v}\cdot\mathbf{a})>0,$ which implies increasing $|\mathbf{v}|,$
since $d(\mathbf{v}^{2})/dt=2(\mathbf{v}\cdot\mathbf{a}).$ This is the
scenario displayed in Figs. 6 and 7, which shows the behavior of
$(\mathbf{v}\cdot\mathbf{a})$ and $|\mathbf{v}|$ for the situation of Fig. 1.
When $\mathrm{v}_{y}(0)<0,$ force and velocity make initially an obtuse angle,
implying $(\mathbf{v}\cdot\mathbf{a})<0$, which yields a decreasing
$|\mathbf{v}|$. Such scenario can be seen in the initial part of the graphs of
Figs. 8 and 9, corresponding to the situation represented in Figs. 4 and 5.
The absolute value $|\mathbf{v}|$ decreases while
$(\mathbf{v}\cdot\mathbf{a})<0$, behavior that changes when
$\mathrm{v}_{y}=0$. After $\mathrm{v}_{y}$ changes sign, becoming positive,
one has $(\boldsymbol{f}\cdot\mathbf{v})>0$ and
$(\mathbf{v}\cdot\mathbf{a})>0,$ so that $|\mathbf{v}|$ increases continuously
to its asymptotic value.

###  II.2 Acceleration behavior

In Ref. Redding , it is claimed that the acceleration $\boldsymbol{a}$ is a
decreasing function of velocity, tending monotonically to zero when
$\beta=\mathrm{v}/c$ tends to unity. More specifically, we can say that
$\boldsymbol{a}$ magnitude goes to zero when the velocity components tend to
their respective asymptotic values, but it does not go in a monotonic way.
This fact, already highlighted in the last section, was not noticed in Ref.
Redding  and will be explored in this section. Taking the scalar product of
Eq. 1 with the acceleration,

$\boldsymbol{f}\cdot\boldsymbol{a}=m\boldsymbol{a}^{2}+m\gamma^{2}(\mathbf{v}\cdot\boldsymbol{a})^{2}\boldsymbol{/}c^{2},$
(15)

and
$(\boldsymbol{f}\cdot\mathbf{a}\boldsymbol{)}=\left|\boldsymbol{f}\right|\left|\mathbf{a}\right|\cos\phi,$
$(\mathbf{v}\cdot\mathbf{a})=\left|\mathbf{v}\right|\left|\mathbf{a}\right|\cos\sigma,$
the resulting magnitude of the acceleration is

$\left|\boldsymbol{a}\right|=\frac{\left|\boldsymbol{f}\right|\cos\phi}{m[1+\gamma^{2}\left|\mathbf{v}\right|^{2}\cos\sigma\boldsymbol{/}c^{2}]},\text{
}$ (16)

which depends on $\left|\mathbf{v}\right|^{2}$ and on the factors $\cos\phi$
and $\cos\sigma,$ that also change with time. As $\left|\mathbf{v}\right|^{2}$
does not have a monotonic behavior, the same can be expected for
$\left|\boldsymbol{a}\right|.$ Further, Ficken Ficken  and Gonzalez Gonzalez
claimed that under the action of a relativistic 3-force, there will exist a
negative acceleration in the direction of the smallest force component, which
is not exactly correct. To investigate the conditions in which the
acceleration becomes negative, we differentiate the velocity components (8),
obtaining:

$\displaystyle a_{x}$
$\displaystyle=\frac{c}{(m_{0}^{2}c^{2}+\mathbf{p}^{2})^{3/2}}[\alpha
f_{y}t+f_{x}m_{0}^{2}c^{2}+\alpha p_{0y}],$ (17) $\displaystyle a_{y}$
$\displaystyle=\frac{c}{(m_{0}^{2}c^{2}+\mathbf{p}^{2})^{3/2}}[-\alpha
f_{x}t+f_{y}m_{0}^{2}c^{2}-\alpha p_{0x}],$ (18)

with

$\alpha=\left(f_{x}p_{0y}-f_{y}p_{0x}\right).$ (19)

Eqs. (17) and (18) reveal that the appearance of negative acceleration is not
determined only by the force magnitude or the ratio $R=f_{x}/f_{y},$ as
claimed in refs. Gonzalez  and Ficken , but by a combination of factors
involving the initial momentum components as well. In the general two-
dimensional case, $\boldsymbol{f}=(f_{x},f_{y}),$ with $f_{x}\geq 0,f_{y}\geq
0$, the factor $\alpha$ can take on both signs.

####  II.2.1 Force in two directions

Case 1: we begin analyzing the situation $\alpha>0$ or $R>p_{0x}/p_{0y}$, that
is,

$\text{\ }\frac{f_{x}}{f_{y}}>\frac{p_{0x}}{p_{0y}},$ (20)

for which there are two main possibilities.

Situation (a): for $p_{0y}>0,$ we have $a_{x}>0$ for any time ($t\geq 0$). In
the y-axis, $a_{y}$ can begin positive or negative, depending on the sign of
the factor $(f_{y}m_{0}^{2}c^{2}-\alpha p_{0x}).$ If $a_{y}(0)>0$, this
acceleration component will become negative after
$t=(f_{y}m_{0}^{2}c^{2}-\alpha p_{0x})/\alpha f_{x}.$ Note that $p_{0x}<0$ and
$p_{0y}>0$ also fulfills condition (20). One notices that one can have
$a_{y}(t)<0$ even for $f_{y}>f_{x}$, thus constituting a scenario of negative
acceleration in the direction of the major force.

Situation (b): For $p_{0y}<0,$ we can have $a_{x}(0)<0,$ if
$f_{x}m_{0}^{2}c^{2}<\alpha\left|p_{0y}\right|.$ In such a situation, the
acceleration $a_{x}$ will be negative until $t=-(f_{x}m_{0}^{2}c^{2}+\alpha
p_{0y})/\alpha f_{y},$becoming positive afterwards. Again, one can have
$p_{0x}<0$ or $p_{0x}>0$ in accordance with Eq. (20). So $a_{y}$ can begin
positive or negative, respectively. For $a_{y}(0)>0$, this component will
become negative after the moment $t=(f_{y}m_{0}^{2}c^{2}-\alpha p_{0x})/\alpha
f_{x}.$ This takes place even for $f_{y}>f_{x}.\vskip 12.0pt plus 4.0pt minus
4.0pt$

Case 2: The situation $\alpha<0$ or $R<p_{0x}/p_{0y}$, that is,

$\text{\ }\frac{f_{x}}{f_{y}}<\frac{p_{0x}}{p_{0y}},$ (21)

is similar to case 1, with $a_{x}$ and $a_{y}$ displaying the behavior of one
another. The component $p_{0x}$ determines whether $a_{y}$ is always positive
or starts negative and become positive in the sequel. Similarly, $p_{0y}$ sign
dictates if $a_{x}$ starts positive $\left(p_{0y}<0\text{ or }p_{0y}>0\text{
with }f_{x}m_{0}^{2}c^{2}>-\alpha p_{0y}\right)$ and then becomes negative or
if $a_{x}$ is always negative$\left(p_{0y}>0\text{ and
}f_{x}m_{0}^{2}c^{2}<\left|\alpha p_{0y}\right|\right)$. The time in which the
acceleration components switch signs can be obtained by making the
substitutions $f_{x}\rightarrow f_{y}$ and $p_{0x}\rightarrow p_{0y}$ in the
previous equations of case 1.

Case 3: An interesting situation happens for $\alpha=0,$ that is,
$R=p_{0x}/p_{0y}$, so that

$\displaystyle a_{x}$
$\displaystyle=\frac{c}{(m_{0}^{2}c^{2}+\mathbf{p}^{2})^{3/2}}[f_{x}m_{0}^{2}c^{2}],$
(22) $\displaystyle a_{y}$
$\displaystyle=\frac{c}{(m_{0}^{2}c^{2}+\mathbf{p}^{2})^{3/2}}[f_{y}m_{0}^{2}c^{2}],$
(23)

and no negative acceleration appears at all for any positive force components
( $f_{x}>0,f_{y}>0)$. This case is not included in previous analysis, as well.

####  II.2.2 Force in only one direction

In the circumstance where there is force applied in only one axis, the
acceleration becomes negative (not always) in the axis where there is no
force. For $\boldsymbol{f}=(0,f_{y})$, one has $\alpha=-f_{y}p_{0x},$ and

$\displaystyle a_{x}$
$\displaystyle=\frac{c\alpha}{(m_{0}^{2}c^{2}+\mathbf{p}^{2})^{3/2}}[f_{y}t+p_{0y}],$
(24) $\displaystyle a_{y}$
$\displaystyle=\frac{c}{(m_{0}^{2}c^{2}+\mathbf{p}^{2})^{3/2}}[f_{y}m_{0}^{2}c^{2}-\alpha
p_{0x}].$ (25)

Concerning such a scenario, there are 4 configurations for the acceleration,
to be analysed as follows.

Situation 1: For $p_{0y}\geq 0$ and $p_{0x}\geq 0$ $\left(\alpha\leq
0\right),$ we notice that $a_{x}\leq 0$ and $a_{y}(t)>0$ for any time. This
case is depicted in Fig. [10], corresponding to the situation of Fig. 1, with
$\mathrm{v}_{x}(0)=0.8c,$ $\mathrm{v}_{y}(0)=0.$ Fig. 10 shows clearly that
the component $a_{y}$ decreases monotonically, but the same does not happen
for $a_{x}.$ Since $\mathrm{v}_{y}(0)=0$, from eq. (13) we have $a_{x}(0)=0$.
After the movement begins, we have $\mathrm{v}$${}_{y}>0$, $\mathrm{v}_{x}>0$,
and Eq. (13) provides $a_{x}$ negative, that is, it starts from zero,
decreases, reaches a minimum value and then rises (absolute value diminishes),
tending continuously to zero. The component $a_{y}$, on the other hand, is
always positive, in accordance with Eq. (14).

Figure 10: Graphs of $a_{x}$ $\times$ time (blue), $a_{y}$ $\times$ time
(red), corresponding to the scenario represented in Fig. 1;
$\mathrm{v}_{x}(0)=0.8c,$ $\mathrm{v}_{y}(0)=0.$

Situation 2: For $p_{0y}>0$ and $p_{0x}<0$ $\left(\alpha>0\right),$ we notice
that $a_{x}(t)>0$ and $a_{y}(t)>0$ for any time. As the product
$\mathrm{v}_{x}\mathrm{v}_{y}$ starts negative, we see from eq. (13) that
$a_{x}(0)>0$. The signs of both acceleration components are kept unaltered
during the whole movement in this situation. In this circumstance, despite the
positive $a_{x}(t),$ the absolute value of velocity $\mathrm{v}_{x}$ decreases
steadily, due to the fact that $\mathrm{v}_{x}$ is negative, and therefore the
positive x-acceleration slows down the particle in that direction. This
scenario is depicted in the graph of Fig. 3 and the corresponding acceleration
graph is shown in Fig. 11.

Figure 11: Graphs of $a_{x}$ $\times$ time (blue), $a_{y}$ $\times$ time
(red), corresponding to the scenario represented in Fig. 3;
$\mathrm{v}_{x}(0)=-0.8c,$ $\mathrm{v}_{y}(0)=0.5c.$

Situation 3: For $p_{0y}<0$ and $p_{0x}>0$ $\left(\alpha<0\right),$ we notice
that $a_{x}(0)>0$, but $a_{x}$ becomes negative from $t=-p_{0y}/f_{y}$ on,
while $a_{y}(t)>0$ for any time. This case is depicted in Fig. [12],
corresponding to the scenario of Fig. 4, that is, $\mathrm{v}_{x}(0)=0.8c,$
$\mathrm{v}_{y}(0)=-0.5c.$ As the product $\mathrm{v}_{x}\mathrm{v}_{y}$
begins negative, eq. (13) shows that x-acceleration is initially positive,
$a_{x}(0)>0$, being positive while $\mathrm{v}_{x}\mathrm{v}_{y}<0.$ As
$\mathrm{v}_{y}$ becomes less negative (increases) under the action of the
force, $a_{x}$ decreases, becoming null when $\mathrm{v}_{y}=0.$ From Eq.
(14), one notices that the smaller $\mathrm{v}_{y}^{2}$ is, the bigger is the
factor $(1-\mathrm{v}_{y}^{2}/c^{2})$ and the component $a_{y}$. For this
reason, $a_{y}$ reaches a maximum when $\mathrm{v}_{y}=0$, at the same point
in which $a_{x}=0$. From this point on the product
$\mathrm{v}_{x}\mathrm{v}_{y}$ is positive, so that the acceleration profiles
recover the same behavior of Fig. 10: the component $a_{y}$ diminishes
monotonically; $a_{x}$ decreases, reaches a minimum value and then rises,
tending continuously to zero.

Figure 12: Graphs of $a_{x}$ $\times$ time (blue), $a_{y}$ $\times$ time
(red), corresponding to the scenario represented in Fig. 4;
$\mathrm{v}_{x}(0)=0.8c,$ $\mathrm{v}_{y}(0)=-0.5c.$

Situation 4: For $p_{0y}<0$ and $p_{0x}<0$ $\left(\alpha>0\right),$ we notice
that $a_{x}(0)<0$, but $a_{x}$ becomes positive from $t=-p_{0y}/f_{y}$ on; and
$a_{y}(t)>0$ for any time. This case is depicted in Fig. [13], corresponding
to the scenario of Fig. 5, where $\mathrm{v}_{x}(0)=-0.8c,$
$\mathrm{v}_{y}(0)=-0.5c.$ The product $\mathrm{v}_{x}\mathrm{v}_{y}$ is
initially positive and $a_{x}$ starts negative. Due to the action of the force
$f$, $\mathrm{v}_{x}\mathrm{v}_{y}$ starts to decrease, until $\mathrm{v}_{y}$
and $a_{x}$ become null. From this point on, $a_{x}>0$, reaches a maximum and
then tends continuously to zero. The graph of $a_{x}$ in Fig. 13 is the
negative of the one in Fig. 12, while the acceleration $a_{y}$ displays the
same behavior of the previous case - because it depends only on the absolute
value of velocities. The overall behavior of movement in the x-direction is
the same, acceleration followed by deceleration, with the difference being the
sign of $\mathrm{v}_{x}$ and $a_{x}$. That explains why the graph of $a_{x}$
in Fig. 13 is the negative of the one in Fig. 12.

In all of the previously discussed scenarios, the time at which $a_{x}$
achieves its minimum or maximum can be found by differentiating twice eq. (8)
and setting it equal to zero, which yields

$t=\frac{\sqrt{2m_{0}^{2}c^{2}+2p_{0x}^{2}}-2p_{0y}}{2f_{y}}$ (26)

Figure 13: Graphs of $a_{x}$ $\times$ time (blue), $a_{y}$ $\times$ time
(red), corresponding to the scenario represented in Fig. 5.

The acceleration behavior for the case $\boldsymbol{f}=(0,f_{y})$, shown in
Figs. [10,11,12,13], corresponding to the situations of Figs. [1,3,4,5],
reveal a non-monotonic behavior, in general. By analyzing such graphs, it is
correct to claim that the acceleration absolute value becomes monotonically
decreasing after it achieves a maximum or minimum point, but not since the
beginning of movement.

##  III Final Remarks

In this article, we have revisited some aspects of the accelerated motion of a
particle under the action of a constant relativistic force. The asymptotic
velocity in each axis, $(f_{i}/f)c,$ is proportional to the force component,
$f_{i},$ so that after a very long time, force and velocity always tend to be
parallel. In order to fulfill momentum conservation, it is possible to have
acceleration in one axis where there is no force. Previous works in literature
claimed that a particle will exhibit a negative acceleration component in the
direction of the smallest force component, which is not always true. It was
argued that a negative acceleration can exist in the direction of the biggest
force. We have also shown that the acceleration magnitude does not decreases
monotonically, as claimed before. The acceleration can reach a maximum or a
minimum peak before decaying to zero. Several graphs were plotted illustrating
the behavior of velocity and acceleration components for the case of a single
force in the y-axis, $\boldsymbol{f}=(0,f_{y}).$

## References

  * (1) W. Rindler, “Introduction to Special Relativity”, Claredon Press, Oxford , 1991.
  * (2) J. Aharoni, “The special theory of relativity”, Oxford University Press, 1959.
  * (3) A. P. French, “Special Relativity”, W.W. Norton & Company, New York, 1966.
  * (4) R. Hagedorn, “Relativistic Kinematics”, W.A. Benjamin, INC. New York, 1963.
  * (5) R. A. Martins, ”Force measurement and force transformation in special relativity”, Am. J. Phys. 50, 1008 (1982).
  * (6) B. U. Stewart, ”Simplified relativistic force transformation”, Am. J. Phys. 57, 50 (1979); O. D. Jefimenko, ”Derivation of relativistic force transformation equations from Lorentz force law”, Am. J. Phys. 64, 618 (1996).
  * (7) P. Finkler, ”Relativistic momentum”, Am. J. Phys. 64, 655 (1996); L. Sartori, ”On the derivation of the formula for relativistic momentum”, Am. J. Phys. 62, 280 (1994); P. C. Peters, ”An alternate derivation of relativistic momentum”, Am. J. Phys. 54, 804 (1986).
  * (8) Leslie H. Bryant, ”Simplification of Relativistic Mass Equations”, Am. J. Phys. 25, 484 (1957); L. C. Baird, ”Relativistic mass”, Am. J. Phys. 48, 779 (1980); P. D. Gupta, ”Relativistic mass”, Am. J. Phys. 49, 890 (1981).
  * (9) Y. Chow,” The Kinematics of Two-Particle Collisions in Special Relativity”, Am. J. Phys. 37, 397 (1969); R. H. Bacon, ”Quanta and Material Particles”, Am. J. Phys. 8, 354 (1940).
  * (10) I. R. Lapidus, ”Relativistic collision problem”, Am. J. Phys. 53, 786 (1985); V. J. Menon, R. K. Dubey, and M. Mishra, B. K. Patra, ”Momentum-energy change in elastic relativistic binary collisions”, Am. J. Phys. 77, 447 (2009).
  * (11) J. L. Redding, ”Spatial relations between force and acceleration in relativistic mechanics”, Am. J. Phys. 50, 163 (1982).
  * (12) J. A. Lock, ”Negative acceleration components for a relativistic particle”, Am. J. Phys. 49, 693 (1981).
  * (13) R. C. Tolman, Relativity Thermodynamics and Cosmology, Claredon, Oxford, 1934.
  * (14) R. C. Tolman, Philos. Mag. 22, 458 (1911).
  * (15) P. F. González-Díaz, ”Relativistic negative acceleration components”, American Journal of Physics 46, 932 (1978); https://doi.org/10.1119/1.11520
  * (16) G. W. Ficken Jr., ”A relativity paradox: The negative acceleration component”, American Journal of Physics 44, 1136 (1976); https://doi.org/10.1119/1.10210.