# A Non-Commutative Entropic Optimal Transport Approach to Quantum Composite
Systems at Positive Temperature

Dario Feliciangeli  IST Austria, Am Campus 1, 3400 Klosterneuburg, Austria.
Augusto Gerolin  Department of Theoretical Chemistry, Vrije Universiteit
Amsterdam, The Netherlands.  Lorenzo Portinale  IST Austria, Am Campus 1, 3400
Klosterneuburg, Austria.

(Dedicated to Tito.)

###### Abstract

This paper establishes new connections between many-body quantum systems, One-
body Reduced Density Matrices Functional Theory (1RDMFT) and Optimal Transport
(OT), by interpreting the problem of computing the ground-state energy of a
finite dimensional composite quantum system at positive temperature as a non-
commutative entropy regularized Optimal Transport problem. We develop a new
approach to fully characterize the dual-primal solutions in such non-
commutative setting. The mathematical formalism is particularly relevant in
quantum chemistry: numerical realizations of the many-electron ground state
energy can be computed via a non-commutative version of Sinkhorn algorithm.
Our approach allows to prove convergence and robustness of this algorithm,
which, to our best knowledge, were unknown even in the two marginal case. Our
methods are based on careful a priori estimates in the dual problem, which we
believe to be of independent interest. Finally, the above results are extended
in 1RDMFT setting, where bosonic or fermionic symmetry conditions are enforced
on the problem.

###### Contents

  1. 1 Introduction
  2. 2 Contributions and statements of the main results
    1. 2.1 Duality and minimization of $\mathfrak{F}^{\varepsilon}$
    2. 2.2 Non-commutative Sinkhorn algorithm
    3. 2.3 The symmetric case: one-body reduced density matrix functional theory
  3. 3 Preliminaries and a priori estimates
    1. 3.1 Non-commutative $(\operatorname{H},\varepsilon)$-transforms
    2. 3.2 Vectorial $(\operatorname{H},\varepsilon)$-transforms
    3. 3.3 One-step and Sinkhorn operators
  4. 4 Non-commutative multi-marginal optimal transport
    1. 4.1 Primal and dual functionals: lower bound and structure of the optimizers
    2. 4.2 Proof of Theorem 2.1
    3. 4.3 Stability and the functional derivative of $\mathfrak{F}^{\varepsilon}$
  5. 5 Non-commutative Sinkhorn algorithm
    1. 5.1 Definition of the algorithm
    2. 5.2 Convergence guarantees: proof of Theorem 2.2
  6. 6 One-body reduced density matrix functional theory
    1. 6.1 Fermionic dual problem and Pauli’s exclusion principle
    2. 6.2 Duality theorem for fermionic and bosonic systems

##  1 Introduction

In this work we are interested in studying the ground state energy of a finite
dimensional composite quantum system at positive temperature. In particular,
we focus on the problem of minimizing the energy of the composite system
conditionally to the knowledge of the states of all its subsystems.

The first motivation for this study is physical: it is useful to understand
how one could infer the state of a composite system when one only has
experimental access to the measurement of the states of its subsystems. The
second motivation is mathematical: indeed this problem can be cast as a non-
commutative optimal transport problem, therefore showcasing how several ideas
and concepts introduced in the commutative setting carry through to the non-
commutative framework. Finally, a third motivation comes from the fact that
one-body reduced density matrix functional theory, which is of interest on its
own, can be framed as a special case of our setting.

Let us consider a composite system with $N$ subsystems, each with state space
given by the complex Hilbert space $\mathfrak{h}_{j}$ of dimension
$d_{j}<\infty$, for $j=1,\dots,N$, and denote the state space of the composite
system
$\mathfrak{h}:=\mathfrak{h}_{1}\otimes\mathfrak{h}_{2}\otimes\dots\otimes\mathfrak{h}_{N}$
(with dimension $d=d_{1}\cdot d_{2}\cdot\dots d_{N}$). Further denote by
$\operatorname{H}$ the Hamiltonian to which the whole system is subject and
suppose that
$\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. Finally, suppose to have knowledge of the states
$\bm{\gamma}=(\gamma_{1},\dots,\gamma_{N})$ of the $N$ subsystems, where each
$\gamma_{j}$ is a density matrix over $\mathfrak{h}_{j}$.

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\\}$
$\displaystyle=\sum_{j=1}^{N}\operatorname{Tr}(\operatorname{H}_{j}\gamma_{j})+\mathfrak{F}^{\varepsilon}(\bm{\gamma})$
$\displaystyle:=\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\\},$ (1.1)

where the shorthand notation $\Gamma\mapsto\bm{\gamma}$ denotes the set of
density matrices over $\mathfrak{h}$ with $j$-th marginal equal to
$\gamma_{j}$, and
$S(\Gamma):=\operatorname{Tr}\left(\Gamma\log(\Gamma)\right)$ is the opposite
of the Von Neumann entropy of $\Gamma$ (note that we prefer to adopt the
mathematical sign convention).

Our approach for the study of $\mathfrak{F}^{\varepsilon}(\bm{\gamma})$
borrows ideas from optimal transport and convex analysis, and takes the
following observation as a starting point: the minimization appearing in
$\mathfrak{F}^{\varepsilon}$ can be cast as a non-commutative entropic optimal
transport problem. Indeed, one looks for an optimal non-commutative coupling
$\Gamma$, with fixed non-commutative marginals (i.e. partial traces)
$\bm{\gamma}$, which minimizes the sum of a transport cost (given by
$\operatorname{Tr}(\operatorname{H}_{\text{int}}\Gamma)$) and an entropic
term. In light of this interpretation, setting the quantum problem at positive
temperature $\varepsilon$ corresponds to consider an entropic optimal
transport problem with parameter $\varepsilon$.

Guided by this viewpoint, we first show that $\mathfrak{F}^{\varepsilon}$ has
a dual formulation (see Theorem 2.1 (i)), i.e. that the constrained
minimization appearing in its definition is in duality with an unconstrained
maximization problem (defined in (2.2)). We can then consider any vector
$(U^{\varepsilon}_{1},\dots,U^{\varepsilon}_{N})$ of self-adjoint matrices
which is a maximizer in the dual functional of $\mathfrak{F}^{\varepsilon}$,
whose existence and uniqueness up to trivial transformations we prove in
Theorem 2.1(ii). We refer to such $U^{\varepsilon}_{i}$-s as Kantorovich
potentials and show in Theorem 2.1(iii) that the unique minimizer
$\Gamma^{\varepsilon}$ realizing $\mathfrak{F}^{\varepsilon}(\bm{\gamma})$ can
be written in terms of them as

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

in the case of all the $\gamma_{j}$-s having trivial kernels (in the general
case a very similar formula holds). In this setting,
$\mathfrak{F}^{\varepsilon}$ is continuous and its functional derivative can
be computed in terms of the Kantorovich potentials as

$\displaystyle\frac{\,\mathrm{d}\mathfrak{F}^{\varepsilon}}{\,\mathrm{d}\gamma_{i}}(\bm{\gamma})=U_{i}^{\varepsilon},\quad\text{for
all}\quad i=1,\dots,N,$ (1.3)

as we show in Proposition 2.1.

Furthermore, we introduce the Non-Commutative Sinkhorn algorithm to compute
the optimizer realizing $\mathfrak{F}^{\varepsilon}(\bm{\gamma})$. This
algorithm exploits the shape of the minimizer obtained in (1.2), in order to
construct a sequence $\Gamma^{(k)}$ of density matrices converging to
$\Gamma^{\varepsilon}$ of the form

$\displaystyle\Gamma^{(k)}=\exp\left(\frac{\bigoplus_{i=1}^{N}U_{i}^{(k)}-\operatorname{H}_{\text{int}}}{\varepsilon}\right),$
(1.4)

where the vector $(U_{1}^{(k)},\dots,U_{N}^{(k)})$ is iteratively updated by
progressively imposing that $\Gamma^{(k)}$ has at least one correct marginal.
We prove the convergence and the robustness of this algorithm in Section 5.

It is important to note that studying
$\mathfrak{F}^{\varepsilon}(\bm{\gamma})$, i.e. the constrained minimization
at fixed marginals, can also help solving the unconstrained minimization of
the Hamiltonian $\operatorname{H}$ at positive temperature $\varepsilon$.
Indeed, denoting by $\mathfrak{P}(\mathfrak{h})$ the set of density matrices
over $\mathfrak{h}$, then

$\displaystyle
E^{\varepsilon}(\operatorname{H}):=\inf_{\Gamma\in\mathfrak{P}(\mathfrak{h})}\left\\{\operatorname{Tr}(\operatorname{H}\Gamma)+\varepsilon
S(\Gamma)\right\\}=\inf_{\bm{\gamma}}\left\\{\sum_{j=1}^{N}\operatorname{Tr}(\operatorname{H}_{j}\gamma_{j})+\mathfrak{F}^{\varepsilon}(\bm{\gamma})\right\\}.$
(1.5)

Combining (1.3) and (1.5) allows to write down the Euler–Lagrange equation of
(1.5) recovering its optimizer, i.e. the Gibbs state constructed with
$\operatorname{H}$ at temperature $\varepsilon$.

Our work is not the first to try to extend the theory of optimal transport to
the non-commutative setting. One of the first attempts was carried out by E.
Carlen and J. Maas [11], followed by many others (e.g. [6, 8, 9, 12, 14, 18,
19, 26, 38, 39, 43]). There is an important distinction to be made here.
Commutative optimal transport can be cast equivalently as a static coupling
problem or as a dynamical optimization problem. On the other hand, in the non-
commutative setting it is not clear what is the relation (if any) between the
two interpretations. This singles out a big difference between works that
consider the dynamical formulation of commutative optimal transport as a
starting point (e.g. [6, 11, 12, 14, 38, 39, 43]) and the ones which instead
focus on its static formulation (e.g. [17, 25, 32, 48, 54]).

This paper adopts an even different approach. We consider as a starting point
the Entropic regularization of optimal transport (which is to be considered as
an extension of static optimal transport, see e.g. the survey [32] and
references therein) and introduce its non-commutative counterpart. We carry
out this program by extending the method developed in [20, 21, 27]. See also
Section 5 for a detailed explanation of the multimarginal Sinkhorn algorithm
in the commutative setting, as studied in [20].

In the work [8], the authors study the case of $\varepsilon=0$ temperature and
prove a duality result for the non-commutative problem in the very same spirit
of the Kantorovich duality for the classical Monge problem. The recent work
[53] studies the entropic quantum optimal transport problem as well, adopting,
in constrast to our static approach, a dynamical formulation. Therein, the
author proves a dynamical duality result at positive and zero temperature. To
the best of our knowledge, the present work is the first complete analysis of
the quantum entropic transport problem in the static framework.

As for the Sinkhorn Algorithm, another concept which we borrow from the
commutative setting and extend to the quantum one, its convergence in the
commutative setting was first established in the $N=2$ marginal case [24, 49]
for discrete measures and in [45] for continuous measures (see also [13]). In
the multi-marginal setting, convergence guarantees were obtained for the
discrete case in [15, 31] and for continuous measures in [20, 21]. Other
variants of the Sinkhorn algorithm for (unbalanced) tensor-valued measures or
matrix optimal mass transport have been studied in [43, 46] and do not apply
to our setting. In the context of Computational Optimal Transport, the
entropic regularization and the Sinkhorn algorithm was introduced in [17, 25].

### Enforcing symmetry constraints: One-body Reduced Density Matrix Functional
Theory

We conclude this introduction by briefly discussing the case in which symmetry
conditions are enforced on the problem, either bosonic or fermionic, which we
can also treat (see Section 2.3). In this case, (1) makes sense only for
$\mathfrak{h}_{j}=\mathfrak{h}_{0}$ for all $j=1,\dots,N$ and
$\bm{\gamma}=(\gamma,\dots,\gamma)$ (i.e. the underlying Hilbert spaces and
the marginals must all be the same) and its study can be framed in the context
of One-body Reduced Density Matrix Functional Theory ($1$RDMFT), introduced in
1975 by Gilbert [29] as an extension of the Hohenberg-Kohn (Levy-Lieb)
formulation of Density Functional Theory (DFT) [30, 33, 34]. In the last
decades, DFT and $1$RDMFT have been standard methods for numerical electronic
structure calculations and are to be considered a major breakthrough in fields
ranging from materials science to chemistry and biochemistry.
In both these theories one tries to approximate a complicated N-particle
quantum system by studying one-particle objects, namely one-body densities in
the case of DFT and one-body reduced density matrices in the case of 1RDMFT,
by using a two-steps minimization analogous to the one introduced in (1.5).

It is interesting to see that the well-known Pauli principle (see e.g. [35,
Theorem 3.2]) , which provides necessary and sufficient conditions for
$\gamma$ to be the one-body reduced density matrix of an $N$-body
antisymmetric density matrix, finds a variational interpretation in our
discussion. Indeed, in the antisymmetric case we show (see Proposition 2.1)
that $\gamma$ satisfies the Pauli principle (resp. satisfies the Pauli
principle strictly) if and only if the supremum of the dual functional of
$\mathfrak{F}^{\varepsilon}$ is finite (resp. is attained), as it is to be
expected.

Other extensions of DFT have been considered, including Mermin’s Thermal
Density Functional Theory [37], Spin DFT [52], and Current DFT [51]. Physical
and computational aspects of 1RDMFT have been investigated in [1, 2, 3, 4, 7,
36, 40, 41, 42, 44, 47, 50]. A framework for 1RDMFT for Bosons at zero
temperature was recently introduced in [5] (see also [28] and references
therein for a recent review). In particular, the first exchange-correlation
energy in density-matrix functional theory was introduced by Müller [40],
leading to mathematical results [22, 23].

### Organisation of the paper

The paper is divided as follows: in Section 2 we introduce the framework, the
main definitions, and present our main results Theorem 2.1, Theorem 2.2, and
Theorem 2.7. In Section 3 we introduce and develop the technical tools needed
to prove our main results, in particular we define the notion of non-
commutative $(\operatorname{H},\varepsilon)$-transform (see Section 3.1) and
prove a stability and differentiability result for the primal problem in
Proposition 2.1. In Section 4, Section 5, and Section 6 we build upon Section
3 and prove our main results, respectively, Theorem 2.1, Theorem 2.2, and
Theorem 2.7.

##  2 Contributions and statements of the main results

The main contributions of this work consist in

  * •

Theorem 2.1, which represents a duality result for the functional
$\mathfrak{F}^{\varepsilon}$ (whose definition is recalled below in equation
(2.1)). Theorem 2.1 also includes the characterization of the optimizers of
$\mathfrak{F}^{\varepsilon}$ (and of its dual functional).

  * •

The introduction of a non-commutative Sinkhorn algorithm, which can be used to
compute the aforementioned optimizers. We also prove convergence and
robustness of this algorithm in Theorem 2.2.

  * •

The generalization of Theorem 2.1 to the case of bosonic or fermionic systems,
stated in Theorem 2.7. This also allows to give an interesting variational
characterization of the Pauli exclusion principle (see Proposition 2.1).

In our analysis, the main tool is the non-commutative
$(\operatorname{H},\varepsilon)$-transform (introduced in Section 3.1), which
allows to obtain a priori estimates on maximizing sequences of Kantorovich
potentials, yielding compactness. Although this approach is not strictly
necessary to prove duality (Theorem 2.1) in our finite dimensional setting, we
believe these estimates to be of independent interests. Moreover, they are
fundamental to prove the convergence of the non-commutative Sinkhorn algorithm
(Theorem 2.2).

We now proceed to introduce our setting and state our main contributions.

###  2.1 Duality and minimization of $\mathfrak{F}^{\varepsilon}$

We recall that in this case we simply work with a general composite system,
with no symmetry constraints enforced. For $d\in\mathbb{N}$, we shall denote
by $\mathcal{M}^{d}=\mathcal{M}^{d}(\mathbb{C})$ the set of all $d\times d$
complex matrices, by $\mathcal{S}^{d}$ the hermitian elements of
$\mathcal{M}^{d}$, and by $\mathcal{S}_{\geq}^{d}$ (respectively
$\mathcal{S}_{>}^{d}$) the set of all the positive semidefinite (resp.
positive definite) elements of $\mathcal{S}^{d}$. With a slight abuse of
notation, we denote by $\operatorname{Tr}$ the trace operator on
$\mathcal{M}_{d}$ for any dimension $d$. Furthemore, for any Hilbert space
$\mathfrak{h}$, we denote by $\mathfrak{P}(\mathfrak{h})$ the set of density
matrices over $\mathfrak{h}$, namely the positive self-adjoint operators with
trace one. For simplicity, we shall also use the notation
$\mathfrak{P}^{d}=\mathfrak{P}(\mathbb{C}^{d})$. For every $N\in\mathbb{N}$ we
adopt the notation $[N]:=\\{1,...,N\\}$.

Our main object of study is the minimisation problem for $N\in\mathbb{N}$,
$i\in[N]$, $\gamma_{i}\in\mathfrak{P}^{d_{i}}$,
$\operatorname{H}\in\mathcal{S}^{\bm{d}}$

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

where $d_{i}\in\mathbb{N}$, $\bm{d}:=\prod_{i=1}^{N}d_{i}$,
$\bm{\gamma}:=(\gamma_{i})_{i\in[N]}$, and $\Gamma\mapsto\bm{\gamma}$ means
that the $i$-th marginal (3.4) of $\Gamma$ is equal to $\gamma_{i}$. This
coincides with the Definition of $\mathfrak{F}^{\varepsilon}$ given in (1).

The natural space to work with is given by
$\mathcal{O}:=\bigotimes_{i=1}^{N}\big{(}\ker\gamma_{i}\big{)}^{\perp}$ where
for simplicity we set $\hat{d}_{i}:=(d_{i}-\dim\ker\gamma_{i})$ and
$\mathbf{\hat{d}}:=\prod_{i=1}^{N}\hat{d}_{i}$. We also denote by
$\operatorname{H}_{\mathcal{O}}$ the restriction of $\operatorname{H}$ to the
subspace $\mathcal{O}$. The corresponding dual problem is defined as

$\displaystyle\mathfrak{D}^{\varepsilon}(\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}_{\mathcal{O}}}{\varepsilon}\bigg{]}\bigg{)}\
:\ U_{i}\in\mathcal{S}^{\hat{d}_{i}}\Bigg{\\}}+\varepsilon\,,$ (2.2)

where $\bigoplus$ denotes the Kronecker sum (3.5).

Our first result is a duality result and serves also as a characterization of
the minimizers in (2.1). Note that, throughout the whole paper, when no
confusion can arise, we shall use the slightly imprecise notation $\alpha
1\\!\\!1=\alpha$ for $\alpha\in\mathbb{C}$.

######  Theorem 2.1 (Duality).

Let $\varepsilon>0$, $N\in\mathbb{N}$, and
$\operatorname{H}\in\mathcal{S}^{\bm{d}}$. For fixed
$\bm{\gamma}=(\gamma_{i}\in\mathfrak{P}^{d_{i}})_{i\in[N]}$, consider the
primal and dual problems $\mathfrak{F}^{\varepsilon}(\bm{\gamma})$,
$\mathfrak{D}^{\varepsilon}(\bm{\gamma})$ as in (2.1), (2.2) respectively. We
then have that

  * (i)

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

  * (ii)

$\mathfrak{D}^{\varepsilon}(\bm{\gamma})$ admits a maximizer
$\\{U_{i}^{\varepsilon}\in\mathcal{S}^{\hat{d}_{i}}\\}_{i=1}^{N}$, which is
unique up to trival translation. Precisely, if
$\\{\tilde{U}_{i}^{\varepsilon}\in\mathcal{S}^{\hat{d}_{i}}\\}_{i=1}^{N}$ is
another maximizer, then
$\tilde{U}_{i}^{\varepsilon}-U_{i}^{\varepsilon}=\alpha_{i}\in\mathbb{R}$ with
$\sum_{i}\alpha_{i}=0$.

  * (iii)

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

$\displaystyle{\Gamma}^{\varepsilon}=\exp\left(\frac{\bigoplus_{i=1}^{N}U_{i}^{\varepsilon}-\operatorname{H}_{\mathcal{O}}}{\varepsilon}\right)\quad\text{on}\;\mathcal{O}$
(2.3)

and $\Gamma^{\varepsilon}=0$ on $\mathcal{O}^{\perp}$.

The proof of the existence of maximizers for the dual problem follows the
direct method of Calculus of Variations. In analogy with [20, 21], where the
notion of commutative $(c,\varepsilon)$-trasform is introduced, we define the
non-commutative $(\operatorname{H},\varepsilon)$-transform (see Section 3.1).
We use this tool to obtain a priori estimates on $U$ and infer compactness of
the maximizing sequences of Kantorovich potentials.

As a byproduct of the a priori estimates obtained in Section 3.1, it is
possible to prove a stability result (with respect to the marginals) for the
Kantorovich potentials and compute the Frechét derivative of
$\mathfrak{F}^{\varepsilon}(\cdot)$. This is the content of the following
proposition, which is proved in Section 4. For simplicity, we here assume that
the marginals have trivial kernel. With a bit more effort, and arguing as in
Theorem 2.1 (see also Remark 3.6), one can obtain a similar result in the
general setting as well.

######  Proposition 2.1 (Stability and differentiability of
$\mathfrak{F}^{\varepsilon}(\cdot)$).

Fix $\varepsilon>0$ and assume $\ker(\gamma_{i})=\\{0\\}$.

  1. $(i)$

_Stability_ : if $\bm{\gamma}^{n}=(\gamma^{n}_{i})_{n\in\mathbb{N}}$,
$\gamma^{n}_{i}\subset\mathfrak{P}^{d_{i}}$ is a sequence of density matrices
converging to $\bm{\gamma}=(\gamma_{i})_{n\in\mathbb{N}}$ as $n\to\infty$,
then any sequence of Kantorovich potentials $\bm{U}^{\varepsilon,n}$
converges, up to subsequences and renormalisation, to a Kantorovich potential
$\bm{U}^{\varepsilon}$ for $\mathfrak{F}^{\varepsilon}(\bm{\gamma})$.
Therefore, the functional $\mathfrak{F}^{\varepsilon}(\cdot)$ is continuous.

  2. $(ii)$

_Frechét differential_ : $\mathfrak{F}^{\varepsilon}(\cdot)$ is Fréchet
differentiable and for every $i\in[N]$ it holds

$\displaystyle\Big{(}\dfrac{d\mathfrak{F}^{\varepsilon}}{d\gamma_{i}}\Big{)}_{\bm{\gamma}}(\sigma)=\operatorname{Tr}\big{(}\bm{U}^{\varepsilon}_{i}\sigma\big{)},\quad\forall\sigma\in\mathcal{S}^{d_{i}},\;\operatorname{Tr}(\sigma)=0\,,$
(2.4)

where $\bm{U}^{\varepsilon}$ is a Kantorovich potential for
$\mathfrak{F}^{\varepsilon}(\bm{\gamma})$.

As derived in [29] and explained, for instance, in [41], the relevance of the
functional derivative in the 1RDMFT case is to find an eigenvalue equation to
find an efficient optimization for the one-particle eigenvalue equations.

###  2.2 Non-commutative Sinkhorn algorithm

The second contribution of this work is to introduce and prove the convergence
of a non-commutative Sinkhorn algorithm (see Section 5), aimed at computing
numerically the optimal density matrix $\Gamma^{\varepsilon}$ and the
corresponding Kantorovich potentials $\\{U_{i}^{\varepsilon}\\}_{i}$.

For this purpose, we define non-commutative
$(\operatorname{H},\varepsilon)$-transform operators, which extend the notion
of $(c,\varepsilon)-$transforms as introduced in [20] (see also Section 5 for
a detailed explanation). Note that the
$(\operatorname{H},\varepsilon)$-transform also depends on $\bm{\gamma}$, but
we omit this dependence as $\bm{\gamma}$ is a fixed parameter of the problem.

For $i\in[N]$ and $\varepsilon>0$, we consider the operators
$\mathcal{T}^{\varepsilon}_{i}:{}_{j=1}^{N}\mathcal{S}^{\hat{d}_{j}}\to{}_{j=1}^{N}\mathcal{S}^{\hat{d}_{j}}$
of the form

$\displaystyle\bm{U}:=(U_{1},\dots,U_{N}),\quad\big{(}\mathcal{T}^{\varepsilon}_{i}(\bm{U})\big{)}_{j}=\begin{cases}U_{j}&\text{if
}j\neq i,\\\ \mathfrak{T}_{i}^{\varepsilon}(U_{1},\dots,U_{i-1},U_{i+1},\dots
U_{N})&\text{if }j=i\end{cases}$

where $\mathfrak{T}_{i}^{\varepsilon}$ is defined implicitly via

$\displaystyle\operatorname{P}_{i}\left[\exp\left(\frac{\bigoplus_{j=1}^{N}\big{(}\mathcal{T}^{\varepsilon}_{i}(\bm{U})\big{)}_{j}-\operatorname{H}_{\mathcal{O}}}{\varepsilon}\right)\right]=\gamma_{i}$
(2.5)

and $\operatorname{P}_{i}$ denotes the $i$-th marginal operator, obtained by
tracing out all but the $i$-th coordinate, see (3.4). In Section 5, we show
that the maps $\mathcal{T}^{\varepsilon}_{i}$ are well-defined, i.e. the
equation (2.5) admits a unique solution
$\mathcal{T}^{\varepsilon}_{i}(\bm{U})$.

Note that, by construction, the matrix
$\exp\left(\bigoplus_{i=1}^{N}\left(\left(\mathcal{T}^{\varepsilon}_{i}(\bm{U})\right)_{j}-\operatorname{H}_{\mathcal{O}}\right)/\varepsilon\right)\in\mathfrak{P}^{\mathbf{\hat{d}}}$
and it has the $i$-th marginal equal to $\gamma_{i}$. The non-commutative
Sinkhorn algorithm is then defined by iterating this procedure for every
$i\in[N]$. We define the one-step Sinkhorn map as

$\displaystyle\begin{gathered}\tau:{}_{j=1}^{N}\mathcal{S}^{\hat{d}_{j}}\to{}_{j=1}^{N}\mathcal{S}^{\hat{d}_{j}},\\\
\tau(\bm{U}):=(\mathcal{T}^{\varepsilon}_{N}\circ\dots\circ\mathcal{T}^{\varepsilon}_{1})(\bm{U}).\end{gathered}$

The Sinkhorn algorithm is obtained by iteration of the map $\tau$ and this is
sufficient to guarantee that the limit point of the resulting sequence is an
optimizer for the dual problem (2.2), as stated in the following Theorem.

######  Theorem 2.2 (Convergence of the non-commutative Sinkhorn algorithm).

Fix $\varepsilon>0$. The definition (2.5) of the operators
$\mathcal{T}^{\varepsilon}_{i}$ is well-posed. Additionally, for any initial
matrix
$\bm{U}^{(0)}=(U_{1},\dots,U_{N})\in{}_{j=1}^{N}\mathcal{S}^{\hat{d}_{j}}$,
there exist $\bm{\alpha}^{k}\in\mathbb{R}^{N}$ with
$\sum_{i=1}^{N}\bm{\alpha}^{k}_{i}=0$ such that

$\displaystyle\bm{U}^{(k)}:=\tau^{k}(\bm{U}^{(0)})+\bm{\alpha}^{k}\rightarrow{\bm{U}}^{\varepsilon}\quad\text{as
}k\to+\infty,$ (2.6)

where ${\bm{U}}^{\varepsilon}=(U_{1}^{\varepsilon},\dots,U_{N}^{\varepsilon})$
is optimal for the dual problem and
$\tau^{k}:=\underbrace{\tau\circ\dots\circ\tau}_{k\text{-times}}$.

Consequently, if one defines for $k\in\mathbb{N}$

$\displaystyle\Gamma^{(k)}:=\exp\left(\frac{\bigoplus_{i=1}^{N}(\bm{U}^{(k)})_{i}-\operatorname{H}_{\mathcal{O}}}{\varepsilon}\right)\quad\text{on
}\mathcal{O},$ (2.7)

and $0$ on $\mathcal{O}^{\perp}$, then $\Gamma^{(k)}\to{\Gamma^{\varepsilon}}$
as $k\to+\infty$ where ${\Gamma^{\varepsilon}}$, ${\bm{U}}^{\varepsilon}$
satisfy (2.3). In particular, $\Gamma^{\varepsilon}$ is optimal for
$\mathfrak{F}^{\varepsilon}(\bm{\gamma})$.

######  Remark 2.3 (Renormalisation).

In the previous theorem, a renormalisation procedure is needed in order to
obtain compactness for the dual potentials $\bm{U}^{k}$. Nonetheless, due to
the fact that $\sum_{i=1}^{N}(\bm{\alpha}^{k})_{i}=0$ and by the properties of
the operator $\bigoplus$, we observe that for $k\in\mathbb{N}$, the equality

$\displaystyle\Gamma^{(k)}=\exp\left(\frac{\bigoplus_{i=1}^{N}\tau^{k}(\bm{U})_{i}-\operatorname{H}_{\mathcal{O}}}{\varepsilon}\right)\quad\text{on
}\mathcal{O}$

is also satisfied. In fact, this shows that no renormalisation procedure is
needed at the level of the primal problem, i.e. for the density matrices
$\Gamma^{(k)}$.

######  Remark 2.4.

(Umegaki Relative entropies) Similar results can be obtained if instead of the
Von Neumann entropy one uses the quantum Umegaki relative entropy with respect
to a reference density matrix with trivial kernel. Specifically, suppose that
$m_{i}\in\mathcal{S}^{d_{i}}$ with $\ker m_{i}=\\{0\\}$. Then one can consider
the minimisation problem

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

where we set $\bm{m}:=\bigotimes_{i=1}^{N}m_{i}$ and
$S(\Gamma|\bm{m}):=\operatorname{Tr}(\Gamma(\log\Gamma-\log\bm{m})$ denotes
the relative entropy of $\Gamma$ with respect to $\bm{m}$. The functional
$\mathfrak{F}^{\varepsilon}$ defined in (2.1) corresponds to the case $\bm{m}$
equals the identity matrix. The corresponding dual functional
$\mathfrak{D}^{\varepsilon}$ as defined in (2.2) is replaced by

$\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}}{\varepsilon}\bigg{]}\bigg{)}\
:\ U_{i}\in\mathcal{S}^{d_{i}}\Bigg{\\}}+\varepsilon,$

for a modified matrix
$\operatorname{H}_{\bm{m}}^{\varepsilon}:=\operatorname{H}-\varepsilon\log\bm{m}$
(restricted to $\mathcal{O}$ in the case of non-trival kernels). It is easy to
see that our approach can also be used in this case. In particular, performing
a change of variables in the dual potentials of the form
$\tilde{\bm{U}}=\bm{U}+\varepsilon\log\bm{m}$ and using that
$S(\Gamma|\text{Id})=S(\Gamma|\bm{m})+\sum^{N}_{i=1}\left[S(\gamma_{i})-S(\gamma_{i}|m_{i})\right]$,
one readily derives the validity of the same results obtained in Theorem 2.1
and Theorem 2.2, with the substitution of $\operatorname{H}$ with
$\operatorname{H}_{\bm{m}}^{\varepsilon}$.

###  2.3 The symmetric case: one-body reduced density matrix functional theory

We are able to obtain the duality results stated above also in the symmetric
cases (either bosonic or fermionic). For given $d,N\in\mathbb{N}$, we set
$\bm{d}=d^{N}$. We consider the bosonic (resp. fermionic) projection operator
$\Pi_{+}$ (resp. $\Pi_{-}$)

$\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}\,,$
(2.8)

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}$, therefore
$\bigwedge_{i=1}^{N}\mathbb{C}^{d}\neq\\{0\\}$ if and only if $N\leq d$. We
denote by

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

the set of bosonic and fermionic density matrices. We fix
$\operatorname{H}\in\mathcal{S}^{\bm{d}}$ such that

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

where the $\operatorname{S}_{i}$ are the permutation operators in Definition
3.3. It is well-known that there exists $\Gamma\in\mathfrak{P}_{-}^{d}$ such
that $\Gamma\mapsto\gamma$ (where $\Gamma\mapsto\gamma$ means that $\Gamma$
has all marginals equal to $\gamma$) if and only if $\gamma$ satisfies the
Pauli exclusion principle, i.e. if and only if $\gamma\in\mathfrak{P}^{d}$ and
$\gamma\leq 1/N$ (see for example [35, Theorem 3.2]).

######  Definition 2.5 (Bosonic and fermionic primal problems).

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\\}\,.$ (2.11)

For any $\gamma\in\mathfrak{P}^{d}$ such that $\gamma\leq 1/N$, we define the
fermionic 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\\}\,.$ (2.12)

An analysis of the extremal points and the existence of the minimizer in
(2.11) and (2.12) have been carried out in [16] for the zero temperature case,
and in [28] in the positive temperature case. As in the non-symmetric case, we
consider the associated bosonic and fermionic dual problems. 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}$ .

######  Definition 2.6 (Bosonic and fermionic dual problems).

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\,.$
(2.13)

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\\}\,.$ (2.14)

We note that a priori $\mathfrak{D}_{\varepsilon}^{-}(\gamma)$ can be defined
for any $\gamma\in\mathfrak{P}^{d}$, whereas
$\mathfrak{F}^{\varepsilon}_{-}(\gamma)$ is only well defined for
$\gamma\in\mathfrak{P}^{d}$ such that $\gamma\leq 1/N$. This constraint on the
primal problem naturally translates to an admissibility condition in order to
have $\mathfrak{D}_{\varepsilon}^{-}(\gamma)<\infty$. To ensure the existence
of a maximizer for $\text{D}_{\gamma}^{-,\varepsilon}$ we further need to
impose $\gamma<1/N$. The following proposition gives an interesting and
variational point of view of the Pauli principle, and it is proved in Section
6.1.

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

We have the following equivalences:

  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 maximiser $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}$.

Finally we state the duality result in the fermionic and bosonic setting.

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

Let $\operatorname{H}\in\mathcal{S}^{\bm{d}}$ satisfying (2.10).

  * (i)

For any given $\gamma\in\mathfrak{P}^{d}$, such that $\gamma\leq\frac{1}{N}$,
the fermionic primal and dual problems coincide, thus
$\mathfrak{F}^{\varepsilon}_{-}(\gamma)=\mathfrak{D}_{\varepsilon}^{-}(\gamma)$.
Moreover, if $0<\gamma<\frac{1}{N}$ then $\emph{D}_{\gamma}^{-,\varepsilon}$
admits a unique maximizer $U_{-}^{\varepsilon}$ such that

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

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

  * (ii)

For any given $\gamma\in\mathfrak{P}^{d}$, the bosonic primal and dual
problems coincide, thus
$\mathfrak{F}^{\varepsilon}_{+}(\gamma)=\mathfrak{D}_{\varepsilon}^{+}(\gamma)$.
Moreover, if $\gamma>0$, $\emph{D}_{\gamma}^{+,\varepsilon}$ admits a unique
maximizer $U_{+}^{\varepsilon}$ such that

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

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

##  3 Preliminaries and a priori estimates

We start this section by recalling the setting and the notation. For
$d\in\mathbb{N}$, we denote by $\mathcal{M}^{d}=\mathcal{M}^{d}(\mathbb{C})$
the set of all $d\times d$ complex matrices, by $\mathcal{S}^{d}$ the
hermitian elements of $\mathcal{M}^{d}$, and by $\mathcal{S}_{\geq}^{d}$
(respectively $\mathcal{S}_{>}^{d}$) the set of all the positive semidefinite
(positive definite) elements of $\mathcal{S}^{d}$. With a slight abuse of
notation, we denote by $\operatorname{Tr}$ the trace operator on
$\mathcal{M}_{d}$ for any dimension $d$. Furthemore, we denote by
$\mathfrak{P}^{d}$ the set of $d\times d$ density matrices, namely the
matrices in $\mathcal{S}_{\geq}^{d}$ with trace one. For the sake of notation,
for every $N\in N$ we denote by $[N]:=\\{1,...,N\\}$.

For a given $N\in\mathbb{N}$ and $(d_{i})_{i=1}^{N}\subset\mathbb{N}$, we
consider for any $i\in[N]$ the injective maps

$\displaystyle\begin{gathered}Q_{i}:\mathcal{M}^{d_{i}}\to\mathcal{M}^{\bm{d}}=\bigotimes_{j=1}^{N}\mathcal{M}^{d_{j}},\quad\bm{d}:=\prod_{j=1}^{N}d_{j},\\\
\forall
A\in\mathcal{M}^{d_{i}},\quad\operatorname{Q}_{i}(A):=\bigotimes_{j=1}^{N}A_{j},\quad
A_{j}=\begin{cases}A&\text{if }j=i,\\\ 1\\!\\!1&\text{if }j\neq
i.\end{cases}\end{gathered}$ (3.3)

We shall use the same notation also for subsets of $\mathbb{C}^{d}$. I.e., we
also denote by $\operatorname{Q}_{i}$ the map
$\operatorname{Q}_{i}:\mathbb{C}^{d_{i}}\to\mathbb{C}^{\bm{d}}$ defined as

$\displaystyle\forall
K\subset\mathbb{C}^{d_{i}},\quad\operatorname{Q}_{i}(K):=\bigotimes_{j=1}^{N}K_{j}\subset\mathbb{C}^{\bm{d}},\quad
K_{j}=\begin{cases}K&\text{if }j=i,\\\ \mathbb{C}^{d_{j}}&\text{if }j\neq
i.\end{cases}$

The marginal operators are the left-inverse of the $\operatorname{Q}_{i}$,
namely $\operatorname{P}_{i}:\mathcal{M}^{\bm{d}}\to\mathcal{M}^{d_{i}}$,
where for every $\Gamma\in\mathcal{M}^{\bm{d}}$,
$\operatorname{P}_{i}(\Gamma)\in\mathcal{M}^{d_{i}}$ is defined by duality as

$\displaystyle\operatorname{Tr}(\operatorname{P}_{i}(\Gamma)A)=\operatorname{Tr}\Big{(}\Gamma\operatorname{Q}_{i}(A)\Big{)},\quad\forall
A\in\mathcal{M}^{d_{i}}.$ (3.4)

######  Remark 3.1.

Observe that $\operatorname{Tr}(\operatorname{P}_{i}(A))=\operatorname{Tr}(A)$
for every $i=1,\dots,N$ and $A\in\mathcal{M}^{\bm{d}}$. Furthermore, if
$A=\bigotimes_{i=1}^{N}A_{i}$ with $\operatorname{Tr}(A_{i})=1$, then
$\operatorname{P}_{i}(A)=A_{i}$.

For a given family of density matrices $\gamma_{i}\in\mathfrak{P}^{d_{i}}$ ,
we use the notation $\bm{\gamma}:=(\gamma_{i})_{i\in[N]}$ and we write
$\Gamma\mapsto\bm{\gamma}=(\gamma_{1},\dots,\gamma_{N})$ whenever
$\Gamma\in\mathfrak{P}^{\bm{d}}$ and $\operatorname{P}_{i}(\Gamma)=\gamma_{i}$
for every $i=[N]$. With the next definitions, we introduce the Kronecker sum
and Permutation operators.

######  Definition 3.2 (Kronecker sum).

For $A_{i}\in\mathcal{M}^{d_{i}}$, we call their Kronecker sum the matrix

$\displaystyle\bigoplus_{i=1}^{N}A_{i}:=\sum_{i=1}^{N}\operatorname{Q}_{i}(A_{i})\in\mathcal{M}^{\bm{d}}$
(3.5)

where $\operatorname{Q}_{i}$ is defined in (3.3).

######  Definition 3.3 (Permutation operators).

For any $i\in[N]$, we introduce the permutation operator
$\operatorname{S}_{i}:\mathcal{M}^{\bm{d}}\approx\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_{i}\in\mathcal{M}^{d_{i}}$ and extended to the whole
$\mathcal{M}^{\bm{d}}$ by linearity.

######  Remark 3.4.

The permutation operators preserve the spectral properties of any operator.
Precisely, $\sigma(\operatorname{S}_{i}(A))=\sigma(A)$ for every $i\in[N]$,
$A\in\mathcal{S}^{\bm{d}}$, where $\sigma(A)$ denotes the spectrum of $A$. In
particular, for every continuous function $f:\mathbb{R}\to\mathbb{R}$, we have
that $\operatorname{Tr}(f(\operatorname{S}_{i}(A)))=\operatorname{Tr}(f(A))$,
for every $A\in\mathcal{S}^{\bm{d}}$.

###  3.1 Non-commutative $(\operatorname{H},\varepsilon)$-transforms

For this section, we specify to the simply case of a two-fold tensor product
and introduce the notion of non-commutative
$(\operatorname{H},\varepsilon)$-transform, which is a central object in our
discussion. We shall see in Section 3.2 how it is then easy to extend this
notion to a general $N$-fold tensor product. We fix
$d,d^{\prime}\in\mathbb{N}$, $0<\alpha\in\mathfrak{P}^{d^{\prime}}$,
$\operatorname{H}\in\mathcal{S}^{dd^{\prime}}$ and $\varepsilon>0$ and define
the map
$\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}:\mathcal{S}^{d}\times\mathcal{S}^{d^{\prime}}\to\mathbb{R}$
as

$\displaystyle\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,V):=\operatorname{Tr}(V\alpha)-\varepsilon\operatorname{Tr}\left(\exp\left[\frac{U\oplus
V-\operatorname{H}}{\varepsilon}\right]\right).$ (3.6)

The $(\operatorname{H},\varepsilon)$-transform of any $U\in\mathcal{S}^{d}$ is
obtained as the maximiser of the map
$\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,\cdot)$.

######  Definition 3.5 ($(\operatorname{H},\varepsilon)$-transform).

We call the unique maximizer of
$\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,\cdot)$ the
$(\operatorname{H},\varepsilon)$-transform of $U\in\mathcal{S}^{d}$. We use
the notation

$\displaystyle\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}:\mathcal{S}^{d}\to\mathcal{S}^{d^{\prime}},\quad\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)=\text{arg}~{}\text{max}\\{\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,V)\
:\ V\in\mathcal{S}^{d^{\prime}}\\}.$ (3.7)

The following lemma shows that the definition of
$(\operatorname{H},\varepsilon)$-transform is indeed well-posed.

######  Lemma 3.1.

Let $U\in\mathcal{S}^{d}$. Then there exists a unique maximizer
$\bar{V}\in\mathcal{S}^{d^{\prime}}$ of
$\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,\cdot)$.

###### Proof.

Fix $U\in\mathcal{S}^{d}$. For every $V\in\mathcal{S}^{d^{\prime}}$, we write
$V=V_{+}-V_{-}$ where $V_{+}$,$V_{-}\in\mathcal{S}^{d^{\prime}}$ denote
respectively the positive and the negative part of $V$ (with respect to its
spectrum). We begin by observing that

$\displaystyle\begin{gathered}\operatorname{Tr}\left(\exp\left[\frac{U\oplus
V-\operatorname{H}}{\varepsilon}\right]\right)\geq\operatorname{Tr}\left(\exp\left[\frac{U\oplus
V-\|\operatorname{H}\|_{\infty}}{\varepsilon}\right]\right)\\\
=\operatorname{Tr}\left(\exp\left(\frac{V}{\varepsilon}\right)\right)\operatorname{Tr}\left(\exp\left(\frac{U}{\varepsilon}\right)\right)\exp\left[\frac{-\|\operatorname{H}\|_{\infty}}{\varepsilon}\right]=:\kappa\operatorname{Tr}\left(\exp\left(\frac{V}{\varepsilon}\right)\right)\geq\kappa
e^{\varepsilon^{-1}\|V_{+}\|_{\infty}},\end{gathered}$ (3.10)

where in the second step we used that $\exp(U\oplus V)=\exp(U)\otimes\exp(V)$
and $\kappa=\kappa(U,\varepsilon,\operatorname{H})$ is a _finite_ constant
depending on $U$, $\varepsilon$, and $\operatorname{H}$. On the other hand, it
clearly holds $\operatorname{Tr}(V\alpha)\leq\|V_{+}\|_{\infty}$ which
combined with (3.10) yields for every $V\in\mathcal{S}^{d^{\prime}}$

$\displaystyle\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,V)\leq\|V_{+}\|_{\infty}-\kappa
e^{\varepsilon^{-1}\|V_{+}\|_{\infty}}.$ (3.11)

Moreover, it is immediate to obtain that

$\displaystyle\begin{aligned}
\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,V)\leq\operatorname{Tr}(V\alpha)=\operatorname{Tr}(V_{+}\alpha)-\operatorname{Tr}(V_{-}\alpha)\leq\|V_{+}\|_{\infty}-\sigma_{\min}(\alpha)\|V_{-}\|_{\infty},\end{aligned}$
(3.12)

where $\sigma_{\min}(\alpha)$ is the spectral gap of $\alpha$, which is
strictly positive by assumption. Let $V_{n}$ be a maximizing sequence for
$\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,\cdot)$, then the
bounds (3.11) and (3.12) imply that $(V_{n})_{+}$, $(V_{n})_{-}$ (and hence
$V_{n}$) are uniformly bounded. Therefore, we can obtain a subsequence (which
we do not relabel) such that
$V_{n}\rightarrow\bar{V}\in\mathcal{S}^{d^{\prime}}.$ The optimality of
$\bar{V}$ follows from the fact that
$\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,\cdot)$ is
continuous and strictly concave (see for example [10]), which also implies
uniqueness. ∎

In the following lemma we use the fact that the
$(\operatorname{H},\varepsilon)$-transform is obtained through a maximization
to show that it can be characterized as the solution of the associated
Euler–Lagrange equation. This property is crucial for the proof of our main
results.

######  Lemma 3.2 (Optimality conditions for the
$(\operatorname{H},\varepsilon)$-transforms).

Given $d,d^{\prime}\in\mathbb{N}$, $0<\alpha\in\mathfrak{P}^{d^{\prime}}$,
$\operatorname{H}\in\mathcal{S}^{dd^{\prime}}$, $\varepsilon>0$, the operator
$\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}$ can be characterized
implicitly by the fact that, for any $U\in\mathcal{S}^{d}$,
$\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)$ is the unique
solution of

$\displaystyle\alpha=\operatorname{P}_{2}\left(\exp\left[\frac{U\oplus\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)-\operatorname{H}}{\varepsilon}\right]\right).$
(3.13)

###### Proof.

Let us pick any $\Lambda\in\mathcal{S}^{d}$ and define
$V_{s}:=\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)+s\Lambda$. By
construction, due to the optimality of
$\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)$, the map

$\displaystyle s\mapsto
g(s):=\operatorname{Tr}(V_{s}\alpha)-\varepsilon\operatorname{Tr}\left(\exp\left[\frac{U\oplus
V_{s}-\operatorname{H}}{\varepsilon}\right]\right)$

must have vanishing derivative at $s=0$. This can be computed [10, Section
2.2] as

$\displaystyle
g^{\prime}(0)=\operatorname{Tr}(\Lambda\alpha)-\varepsilon\operatorname{Tr}\left((I\otimes\Lambda)\exp\left[\frac{U\oplus\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)-\operatorname{H}}{\varepsilon}\right]\right).$
(3.14)

Using the definition of partial trace and the previous formula, we infer

$\displaystyle\operatorname{Tr}\left(\Lambda\Big{(}\alpha-\operatorname{P}_{1}\Big{(}\exp\left[\frac{U\oplus\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)-\operatorname{H}}{\varepsilon}\right]\Big{)}\Big{)}\right)=0$
(3.15)

for every $\Lambda\in\mathcal{S}^{d}$. Note that
$\alpha,U,V_{s},\operatorname{H}$ being self-adjoint, it follows that the
operator

$\displaystyle\alpha-\operatorname{P}_{1}\left(\exp\left[\frac{U\oplus\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)-\operatorname{H}}{\varepsilon}\right]\right)$

is self-adjoint as well. Together with (3.15), this shows (3.13). On the other
hand, since (3.15) is the Euler Lagrange equation associated to the
maximization of the strictly concave functional
$\operatorname{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U,\cdot)$, any
solution of (3.15) is necessarily a maximizer and hence coincides with
$\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)$, by uniqueness (see
Lemma 3.1).

∎

The next step is to obtain some regularity estimates on
$\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)$. To do so, we
extrapolate information from the optimality conditions proved in Lemma 3.2.

######  Proposition 3.1 (Regularity of the
$(\operatorname{H},\varepsilon)$-transform).

Given $d,d^{\prime}\in\mathbb{N}$, $0<\alpha\in\mathfrak{P}^{d^{\prime}}$,
$\operatorname{H}\in\mathcal{S}^{dd^{\prime}}$, $\varepsilon>0$, we define for
all $A\in\mathcal{S}^{d}$ (or $A\in\mathcal{S}^{d^{\prime}}$)

$\displaystyle\lambda_{\varepsilon}(A):=\varepsilon\log\Big{(}\operatorname{Tr}\Big{[}\exp\left(\frac{A}{\varepsilon}\right)\Big{]}\Big{)}.$
(3.16)

Then for every $U\in\mathcal{S}^{d}$ it holds

$\displaystyle\Big{|}\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)-\varepsilon\log\alpha+\lambda_{\varepsilon}(U)1\\!\\!1\Big{|}\leq\|{\rm
H}\|_{\infty}1\\!\\!1,$ (3.17)
$\displaystyle\Big{|}\lambda_{\varepsilon}(U)+\lambda_{\varepsilon}\big{(}\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)\big{)}\Big{|}\leq\|{\rm
H}\|_{\infty},$ (3.18)
$\displaystyle\Big{|}\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)-\varepsilon\log\alpha-\lambda_{\varepsilon}(\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U))1\\!\\!1\Big{|}\leq
2\|{\rm H}\|_{\infty}1\\!\\!1.$ (3.19)

where the inequalities are understood as two-sided quadratic forms bounds.

###### Proof.

Note that (3.19) is an immediate consequence of (3.17) and (3.18) and we shall
therefore only prove the latter two. Let us start with the proof of (3.17). We
know from Lemma 3.2 that for every $U\in\mathcal{S}^{d}$,
$\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)$ satisfies equation
(3.13). By the properties of the partial trace (Remark 3.1) and
$\operatorname{H}\leq\|{\rm H}\|_{\infty}1\\!\\!1$, it follows that

$\displaystyle\begin{aligned} \alpha&\leq e^{\frac{\|{\rm
H}\|_{\infty}}{\varepsilon}}\operatorname{P}_{1}\left(\exp\left[\frac{U\oplus\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)}{\varepsilon}\right]\right)\\\
&=e^{\frac{\|{\rm
H}\|_{\infty}}{\varepsilon}}\operatorname{P}_{1}\left(\exp\left(\frac{U}{\varepsilon}\right)\otimes\exp\left(\frac{\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)}{\varepsilon}\right)\right)\\\
&=e^{\frac{\|{\rm
H}\|_{\infty}}{\varepsilon}}\operatorname{Tr}\left(\exp\left(\frac{U}{\varepsilon}\right)\right)\exp\left(\frac{\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)}{\varepsilon}\right)\,,\end{aligned}$
(3.20)

where in the first inequality we used that $\exp(A\oplus B)=\exp A\otimes\exp
B$. Similarly, using instead the lower bound $\operatorname{H}\geq-\|{\rm
H}\|_{\infty}1\\!\\!1$, from (3.13) we can also obtain

$\displaystyle\alpha\geq e^{\frac{-\|{\rm
H}\|_{\infty}}{\varepsilon}}\operatorname{Tr}\left(\exp\left(\frac{U}{\varepsilon}\right)\right)\exp\left(\frac{\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)}{\varepsilon}\right).$
(3.21)

We can put together the two bounds in (3.20), (3.21) to obtain

$\displaystyle\alpha e^{\frac{-\|{\rm
H}\|_{\infty}}{\varepsilon}}\leq\operatorname{Tr}\left(\exp\left(\frac{U}{\varepsilon}\right)\right)\exp\left(\frac{\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)}{\varepsilon}\right)\leq\alpha
e^{\frac{\|{\rm H}\|_{\infty}}{\varepsilon}}.$ (3.22)

Taking the log in the latter inequalities we conclude the proof of (3.17). If
we instead take the trace of both sides in (3.22), we obtain

$\displaystyle e^{\frac{-\|{\rm
H}\|_{\infty}}{\varepsilon}}\leq\operatorname{Tr}\left(\exp\left(\frac{U}{\varepsilon}\right)\right)\operatorname{Tr}\left(\exp\left(\frac{\mathfrak{T}_{\alpha,\operatorname{H}}^{\varepsilon}(U)}{\varepsilon}\right)\right)\leq
e^{\frac{\|{\rm H}\|_{\infty}}{\varepsilon}},$

and then applying the log, we conclude the proof $\eqref{eq:control_lambda}$.
∎

###  3.2 Vectorial $(\operatorname{H},\varepsilon)$-transforms

In this section, we consider a vectorial version of the
$(\operatorname{H},\varepsilon)$-transforms introduced in the previous
section. This turns out to be a key object in the proof of Theorem 2.1 and
Theorem 2.2, necessary to deal with the multi-marginal setting.

Let us first introduce the general framework, which remains in force
throughout the section. Let $N\in\mathbb{N}$ and $[N]$ be a index set of $N$
elements. For all $i\in[N]$, let $d_{i}\in\mathbb{N}$ and
$\gamma_{i}\in\mathfrak{P}^{d_{i}}$ be density matrices. Set
$\bm{\gamma}:=(\gamma_{i})_{i\in[N]}$, $\bm{d}=\prod_{j=1}^{N}d_{i}$. Finally,
consider a Hamiltonian $\operatorname{H}\in\mathcal{S}^{\bm{d}}$.

######  Remark 3.6.

(Kernels) Without loss of generality, we can assume $\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}$ and
consider the matrix
$\operatorname{H}_{\mathcal{O}}=\Pi_{\mathcal{O}}\operatorname{H}\Pi_{\mathcal{O}}$,
where $\Pi_{\mathcal{O}}$ is the projector onto $\mathcal{O}$.

We therefore assume that $\ker\gamma_{i}=\\{0\\}$ for all $i\in[N]$. In this
section we extend the notion of $(\operatorname{H},\varepsilon)$-transform as
introduced in previous section 3.1 to the multi-marginal setting, and we apply
it to our specific setting. We are interested in the maximization (2.2) of the
dual functional, that we introduce below.

######  Definition 3.7 (Dual Functional).

For any $\bm{U}=(U_{1},\dots,U_{N})\in{}_{j=1}^{N}\mathcal{S}^{d_{j}}$, we
define

$\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.$

######  Remark 3.8.

Note that $\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ is invariant by
translation for any vector $\bm{a}=(a_{1},\dots,a_{N})\in\bm{\mathbb{R}}^{N}$
such that $\sum_{k=1}^{N}a_{k}=0$, i.e.

$\displaystyle\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}+\bm{a})=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}).$

As a consequence of this property, we see in Section 5 that the set of
maximizers is invariant by such transformations (Lemma 4.1).

With the following definition, we introduce the vectorial
$(\operatorname{H},\varepsilon)$-transforms.

######  Definition 3.9 (Vectorial $(\operatorname{H},\varepsilon)$-transform).

For any $i\in[N]$, we define the $i$-th vectorial
$(\operatorname{H},\varepsilon)$-transform $\mathfrak{T}_{i}^{\varepsilon}$ as
the map

$\displaystyle\begin{gathered}\mathfrak{T}_{i}^{\varepsilon}:{}_{j=1,\,j\neq
i}^{N}\mathcal{S}^{d_{j}}\to\mathcal{S}^{d_{i}},\\\
\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})=\mathop{\mathrm{argmax}}_{V\in\mathcal{S}^{d_{i}}}\left\\{\operatorname{Tr}(V\gamma_{i})-\varepsilon\operatorname{Tr}\left(\exp\left[\frac{1}{\varepsilon}\left(U_{1}\oplus\dots\oplus
U_{i-1}\oplus V\oplus U_{i+1}\oplus\dots\oplus
U_{N}-\operatorname{H}\right)\right]\right)\right\\},\end{gathered}$

where for $\bm{U}\in{}_{j=1}^{N}\mathcal{S}^{d_{j}}$, we set
$\hat{\bm{U}}_{i}$ to be the product of all the $U_{j}$ but the $i$-th one,
namely

$\displaystyle\hat{\bm{U}}_{i}:=\big{(}U_{1},\dots,U_{j-1},U_{j+1},\dots,U_{N}\big{)}\in{}_{j=1,\,j\neq
i}^{N}\mathcal{S}^{d_{j}}.$ (3.23)

######  Remark 3.10.

Observe that we can identify the $i$-th vectorial
$(\operatorname{H},\varepsilon)$-transforms with a particular case of the
operators $\mathfrak{T}_{\varepsilon,\operatorname{H},\alpha}$ as introduced
in Section 3.1. Indeed, as a consequence of Remark 3.4 it is straightforward
to see that for $i\in[N]$

$\displaystyle\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})=\mathfrak{T}_{\gamma_{i},\operatorname{S}_{i}(\operatorname{H})}^{\varepsilon}\left(\bigoplus_{j=1,\,j\neq
i}^{N}U_{i}\right),\quad\mathfrak{T}_{\gamma_{i},\operatorname{S}_{i}(\operatorname{H})}^{\varepsilon}:\bigotimes_{j=1,\,j\neq
i}^{N}\mathcal{S}^{d_{j}}\approx\mathcal{S}^{\tilde{d}_{i}}\to\mathcal{S}^{d_{i}},$
(3.24)

where we set $\tilde{d}_{i}:=\prod_{j\neq i}d_{j}$ and the
$\operatorname{S}_{i}$ are the permutation operators in Definition 3.3. This
shows that the definition is well posed (i.e. that the
$\mathop{\mathrm{argmax}}$ appearing in the definition exists and is unique).
Moreover it allows us to extend the validity of the properties of the
$(\operatorname{H},\varepsilon)$-transform shown in Section 3.1 to the
operators $\mathfrak{T}_{i}^{\varepsilon}$, as we shall see in Lemma 3.3 and
Proposition 3.2 below. Note that the dependence on the specific entry $i$ is
reflected in both the use of $\gamma_{i}$ and in the fact that the transform
is performed w.r.t. $\operatorname{S}_{i}(\operatorname{H})$.

###  3.3 One-step and Sinkhorn operators

We use the vectorial $(\operatorname{H},\varepsilon)$-transforms to define
what we call one-step operators and Sinkhorn operators. The first ones map a
vector of $N$ potentials into a vector of $N$ potentials, exchanging its
$i$-th entry with the $i$-th vectorial
$(\operatorname{H},\varepsilon)$-transform applied to the other $N-1$. The
second is simply obtained by composing all the different $N$ one-step
operators.

######  Definition 3.11 (One-step operators).

For $i\in[N]$, we introduce the one-step operators
$\mathcal{T}^{\varepsilon}_{i}$, which are defined by

$\displaystyle\mathcal{T}^{\varepsilon}_{i}:{}_{j=1}^{N}\mathcal{S}^{d_{j}}\to{}_{j=1}^{N}\mathcal{S}^{d_{j}}$
$\displaystyle\bm{U}:=(U_{1},\dots,U_{N})\longmapsto(U_{1},\dots,U_{i-1},\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i}),U_{i+1},\dots,U_{N})=:\mathcal{T}^{\varepsilon}_{i}(\bm{U}).$

The Sinkhorn operator is simply the composition of the $N$ one-step operators
$\mathcal{T}^{\varepsilon}_{i},i\in[N]$.

######  Definition 3.12 (Sinkhorn Operator).

We introduce the Sinkhorn operator $\tau$, defined by

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

######  Remark 3.13.

Note that, by definition of $\tau$, it follows immediately that, for any
$\bm{U}\in{}_{j=1}^{N}\mathcal{S}^{d_{j}}$

$\displaystyle\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\tau(\bm{U}))\geq\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}),$

i.e. applying $\tau$ to any vector increases its energy. Moreover, any
maximizer of $\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ is a fixed point
of $\tau$ (as a consequence of the uniqueness proved in Lemma 3.1). The
converse is also true and implies that the set of maximizers of
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ coincides with the set of fixed
points of $\tau$, see Remark 4.2.

######  Remark 3.14.

Note that for any vector $\bm{a}\in\bm{\mathbb{R}}^{N}$ such that
$\sum_{k=1}^{N}\bm{a}_{k}=0$, one has

$\displaystyle\mathcal{T}^{\varepsilon}_{i}(\bm{U}+\bm{a})=\mathcal{T}^{\varepsilon}_{i}(\bm{U})+\bm{a},$

i.e. $\mathcal{T}^{\varepsilon}_{i}$ commutes with translations by vectors
whose coordinates sum up to zero (notice that this fact is particularly
interesting in light of Remark 3.8). This is a straightforward consequence of
the fact that

$\displaystyle\mathfrak{T}_{i}^{\varepsilon}\left((\widehat{\bm{U}+\bm{a}})_{i}\right)=\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})+a_{i},$

which can be readily verified from the definitions. Trivially, this also
implies

$\displaystyle\tau(\bm{U}+\bm{a})=\tau(\bm{U})+\bm{a}.$

We now take advantage of the observations in Remark 3.10 to deduce properties
for the vectorial $(\operatorname{H},\varepsilon)$-transforms, the one-step
operators, and the Sinkhorn operator. First of all, as a corollary of Lemma
3.2, we characterize the vectorial $(\operatorname{H},\varepsilon)$-transforms
as solutions of implicit equations.

######  Lemma 3.3 (Optimality conditions for vectorial
$(\operatorname{H},\varepsilon)$-transforms).

Let $i\in[N]$, $\varepsilon>0$, $\gamma_{i}\in\mathfrak{P}^{d_{i}}$,
$\operatorname{H}\in\mathcal{S}^{d}$, with $\ker\gamma_{i}=\\{0\\}$. For any
$\bm{U}\in{}_{j=1}^{N}\mathcal{S}^{d_{j}}$, the one step-operator
$\mathcal{T}^{\varepsilon}_{i}(\bm{U})$ (or equivalently the $i$-th vectorial
$(\operatorname{H},\varepsilon)$-transform
$\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})$) is implicitly
characterized as the unique solution of the equation

$\displaystyle\gamma_{i}=\operatorname{P}_{i}\left(\exp\left[\frac{1}{\varepsilon}\left(\bigoplus_{j=1}^{N}(\mathcal{T}^{\varepsilon}_{i}(\bm{U}))_{j}-\operatorname{H}\right)\right]\right)\,.$
(3.25)

###### Proof.

As a consequence of (3.24), we can apply Lemma 3.2 and deduce

$\displaystyle\gamma_{i}$
$\displaystyle=\operatorname{P}_{i}\left(\exp\left[\frac{1}{\varepsilon}\left(\left(\bigoplus_{j=1,\,j\neq
i}^{N}U_{i}\right)\oplus\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})-\operatorname{S}_{i}(\operatorname{H})\right)\right]\right)$
$\displaystyle=\operatorname{P}_{i}\left(\exp\left[\frac{1}{\varepsilon}\left(\bigoplus_{j=1}^{N}(\mathcal{T}^{\varepsilon}_{i}(\bm{U}))_{j}-\operatorname{H}\right)\right]\right)\,,$

where $\operatorname{S}_{i}$ is the $i$-th permutation operator, as defined in
3.3, and in the last equality we used Remark 3.4 and that

$\displaystyle\left(\bigoplus_{j=1,\,j\neq
i}^{N}U_{i}\right)\oplus\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})=\operatorname{S}_{i}\left(\bigoplus_{j=1}^{N}(\mathcal{T}^{\varepsilon}_{i}(\bm{U}))_{j}\right)$

for every $i\in[N]$ and $\bm{U}\in{}_{j=1}^{N}\mathcal{S}^{d_{j}}$. ∎

The next proposition collects the regularity properties of the
$(\operatorname{H},\varepsilon)$-transforms. Once again, they are direct
consequence of the properties proved in the two marginals case, in particular
in Proposition 3.1.

######  Proposition 3.2 (Regularity of the
$(\operatorname{H},\varepsilon)$-transforms).

Let $i\in[N]$, $\varepsilon>0$, $\gamma_{i}\in\mathfrak{P}^{d_{i}}$,
$\operatorname{H}\in\mathcal{S}^{d}$, with $\ker\gamma_{i}=\\{0\\}$. Then for
every $\bm{U}\in{}_{j=1}^{N}\mathcal{S}^{d_{j}}$, for every $i\in[N]$ it holds

$\displaystyle\left|\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})-\varepsilon\log\gamma_{i}+\sum_{j=1,j\neq
i}^{N}\lambda_{\varepsilon}(U_{j})1\\!\\!1\right|\leq\|{\rm
H}\|_{\infty}1\\!\\!1\,,$ (3.26) $\displaystyle\left|\sum_{j=1,j\neq
i}^{N}\lambda_{\varepsilon}(U_{j})+\lambda_{\varepsilon}\left(\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})\right)\right|\leq\|{\rm
H}\|_{\infty}\,,$ (3.27)
$\displaystyle\left|\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})-\varepsilon\log\gamma_{i}-\lambda_{\varepsilon}\left(\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})\right)1\\!\\!1\right|\leq(2\|{\rm
H}\|_{\infty})1\\!\\!1,$ (3.28)

where $\lambda_{\varepsilon}$ is defined in (3.16).

###### Proof.

The proof is a direct application of Proposition 3.1 and the considerations in
Remark 3.10. Precisely, the estimate (3.26) follows from (3.17), (3.27)
follows from (3.18) and (3.28) follows from (3.19), together with the fact
that

$\displaystyle\lambda_{\varepsilon}\left(\bigoplus_{j=1,\,j\neq
i}^{N}U_{j}\right)=\sum_{j=1,j\neq
i}^{N}\lambda_{\varepsilon}\left(U_{j}\right),\quad\forall\bm{U}\in{}_{j=1}^{N}\mathcal{S}^{d_{j}}.$

∎

In light of Remark 3.13, it is reasonable to check whether sequences of the
form $\tau^{k}(\bm{U}_{0})$ are maximizing for
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ and compact. On the other hand,
a priori it is not clear how to obtain compactness for such sequences and
Remark 3.8 shows that there could even exist sequences ‘converging’ to the set
of maximizers which are not compact. It is therefore natural to introduce a
suitable renormalization operator, aimed at retrieving compactness. Note that
any such operator should increase or leave invariant the value of
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ and therefore, by Remark 3.8,
any translation by vectors whose coordinates sum up to zero is a good
candidate.

######  Definition 3.15 (Renormalisation).

Let $\lambda_{\varepsilon}$ be defined as in (3.16). We define the
renormalisation map
$\operatorname{Ren}:{}_{i=1}^{N}\mathcal{S}^{d_{i}}\to{}_{i=1}^{N}\mathcal{S}^{d_{i}}$
as the function

$\displaystyle\operatorname{Ren}(\bm{U})_{i}=\begin{cases}U_{i}-\lambda_{\varepsilon}(U_{i}),&\text{
if }i\in\\{1,\dots,N-1\\}\\\ \displaystyle
U_{N}+\sum_{j=1}^{N-1}\lambda_{\varepsilon}(U_{j}),&\text{ if
}i=N.\end{cases}$

In the following proposition we show that
$\operatorname{Ren}\left(\tau\left({}_{i=1}^{N}\mathcal{S}^{d_{i}}\right)\right)$
is bounded and therefore compact. This shows that the map $\operatorname{Ren}$
is indeed a reasonable renormalization operator for our purposes.

######  Proposition 3.3 (Renormalisation of
$(\operatorname{H},\varepsilon)$-transforms and uniform bounds).

Let $i\in[N]$, $\varepsilon>0$, $\gamma_{i}\in\mathfrak{P}^{d_{i}}$,
$\operatorname{H}\in\mathcal{S}^{d}$, with $\ker\gamma_{i}=\\{0\\}$. Then, for
any $\bm{U}\in\mathcal{S}^{\bm{d}}$, one has that
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\operatorname{Ren}\tau(\bm{U}))\geq\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}),$
and the following bounds hold true:

$\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].$ (3.29)

###### Proof.

First of all, Remark 3.8 and Remark 3.13 trivially yield
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\operatorname{Ren}\tau(\bm{U}))\geq\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}).$
To show (3.29), note that for any $i\in[N]$, $(\tau(\bm{U}))_{i}$ is obtained
applying $\mathfrak{T}_{i}^{\varepsilon}$ to some element of ${}_{j=1,\,j\neq
i}^{N}\mathcal{S}^{d_{j}}$. Therefore, applying (3.28) from Proposition 3.2,
we obtain

$\displaystyle\left\|(\operatorname{Ren}\tau(\bm{U}))_{i}-\varepsilon\log\gamma_{i}\right\|_{\infty}=\left\|(\tau(\bm{U}))_{i}-\varepsilon\log\gamma_{i}-\lambda_{\varepsilon}((\tau(\bm{U}))_{i})\right\|_{\infty}\leq
2\|\operatorname{H}\|_{\infty}$

for every $i\in[N-1]$. Moreover,
$(\tau(\bm{U}))_{N}=\mathfrak{T}_{\varepsilon}^{N}(\widehat{\tau(\bm{U})}_{N})$
and hence, applying (3.26) from Proposition 3.2, we arrive at

$\displaystyle\left\|(\operatorname{Ren}\tau(\bm{U}))_{N}-\varepsilon\log\gamma_{N}\right\|_{\infty}=\left\|(\tau(\bm{U}))_{N}-\varepsilon\log\gamma_{N}+\sum_{j=1}^{N-1}\lambda_{\varepsilon}((\tau(\bm{U}))_{j})\right\|_{\infty}\leq\|\operatorname{H}\|_{\infty},$

which completes the proof. ∎

##  4 Non-commutative multi-marginal optimal transport

In this section we prove Theorem 2.1, our first main result stated in Section
2, exploiting the tools developed in Section 3. Again, we fix the setup, which
remains in force throughout the whole Section 4 and Section 5. Let
$N\in\mathbb{N}$, and for $i\in[N]$ we consider density matrices
$\gamma_{i}\in\mathfrak{P}^{d_{i}}$. Set
$\bm{\gamma}:=(\gamma_{i})_{i\in[N]}$, $\bm{d}=\prod_{j=1}^{N}d_{i}$, and
assume that $\ker\gamma_{i}=\\{0\\}$ (see Remark 3.6). We also fix
$\operatorname{H}\in\mathcal{S}^{\bm{d}}$. In this section, we prove the
Theorem 2.1.

We begin by introducing the primal functional, which appears in the
minimisation (2.1).

######  Definition 4.1 (Primal Functional).

Let $\Gamma\in\mathfrak{P}^{{\bf d}}$ the primal functional is defined by

$\operatorname{F}^{\varepsilon}(\Gamma)=\operatorname{Tr}(\operatorname{H}\Gamma)+\varepsilon
S(\Gamma)=\operatorname{Tr}(\operatorname{H}\Gamma)+\varepsilon\operatorname{Tr}(\Gamma\log\Gamma).$
(4.1)

We also recall the definitions of the primal and the dual problem

$\mathfrak{F}^{\varepsilon}(\bm{\gamma})=\inf\left\\{\operatorname{F}^{\varepsilon}(\Gamma)\
:\ \Gamma\in\mathfrak{P}^{{\bf d}}\text{ and
}\Gamma\mapsto(\gamma_{1},\dots,\gamma_{N})\right\\},$ (4.2)
$\mathfrak{D}^{\varepsilon}(\bm{\gamma})=\sup\left\\{\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U})\
:\ \bm{U}\in{}_{i=1}^{N}\mathcal{S}^{d_{i}}\right\\},$ (4.3)

where the dual functional $\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ is
given in Definition 3.7.

###  4.1 Primal and dual functionals: lower bound and structure of the
optimizers

We begin with the proof of the lower bound for the primal functional (4.1), in
terms of the dual functional (3.7).

######  Proposition 4.1 (Lower bound).

Fix $N\in\mathbb{N}$ and $\varepsilon>0$. For all $i\in[N]$, let
$\gamma_{i}\in\mathfrak{P}^{d_{i}}$ be density matrices,
$\operatorname{H}\in\mathcal{S}^{\bm{d}}$. Then, for all
$\bm{U}\in{}_{i=1}^{N}\mathcal{S}^{d_{i}}$ and every
$\Gamma\in\mathfrak{P}^{\bm{d}}$, $\Gamma\mapsto\bm{\gamma}$ we have that

$\displaystyle\operatorname{F}^{\varepsilon}(\Gamma)\geq\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U})\,.$

###### Proof.

For any $\bm{U}\in{}_{i=1}^{N}\mathcal{S}^{d_{i}}$ and any admissible
$\Gamma\in\mathfrak{P}^{\bm{d}}$, $\Gamma\mapsto\bm{\gamma}$, we can write

$\displaystyle\operatorname{F}^{\varepsilon}(\Gamma)$
$\displaystyle=\operatorname{F}^{\varepsilon}(\Gamma)+\sum_{j=1}^{N}\operatorname{Tr}(U_{j}\gamma_{j})-\operatorname{Tr}\left(\left(\bigoplus_{j=1}^{N}U_{j}\right)\Gamma\right)$
$\displaystyle=\sum_{j=1}^{N}\operatorname{Tr}(U_{j}\gamma_{j})+\varepsilon
S(\Gamma)-\operatorname{Tr}\left(\Gamma\left(\bigoplus_{j=1}^{N}U_{j}-\operatorname{H}\right)\right).$

Let us denote the Hilbert-Schmidt scalar product (on $\mathcal{M}^{\bm{d}}$)
by $\langle\cdot,\cdot\rangle_{HS}$. It follows that

$\displaystyle\operatorname{F}^{\varepsilon}(\Gamma)=\sum_{j=1}^{N}\operatorname{Tr}(U_{j}\gamma_{j})+\varepsilon\left[S(\Gamma)-\langle\Gamma,\overline{Y}\rangle_{HS}\right]\geq\sum_{j=1}^{N}\operatorname{Tr}(U_{j}\gamma_{j})-\varepsilon
S^{*}(\overline{Y}),$ (4.4)

where
$\overline{Y}=\varepsilon^{-1}\left(\bigoplus_{j=1}^{N}U_{j}-\operatorname{H}\right)\in\mathcal{S}^{\bm{d}}$
and, for any $Y\in\mathcal{S}^{\bm{d}}$

$\displaystyle S^{*}(Y):=\sup_{\Gamma\in\mathcal{S}_{\geq}^{\bm{d}}}\\{\langle
Y,\Gamma\rangle_{HS}-S(\Gamma)\\}$

denotes the Legendre transform of $S$ on the subspace
$\mathcal{S}_{\geq}^{\bm{d}}$. This can be explicitly computed as

$\displaystyle S^{*}(Y)=\operatorname{Tr}\left[\exp(Y-1)\right],\quad\forall
Y\in\mathcal{S}^{\bm{d}}.$ (4.5)

For the sake of completeness, let us explain how to prove (4.5). First of all
we show that for any $Y\in\mathcal{S}^{\bm{d}}$ the supremum appearing in the
definition of $S^{*}(Y)$ is attained at some
$\bar{\Gamma}\in\mathcal{S}_{>}^{\bm{d}}$. Indeed, for any $\Gamma\geq 0$
define $\sigma_{+}$ to be the maximum of its spectrum, then it holds

$\displaystyle\langle
Y,\Gamma\rangle_{HS}-S(\Gamma)\leq\bm{d}^{2}\|Y\|_{\infty}\sigma_{+}-\sigma_{+}\log\sigma_{+}-\min_{\mathbb{R}_{+}}\\{x\log
x\\}(\bm{d}^{2}-1)\xrightarrow{\sigma_{+}\to\infty}-\infty.$

This implies that the super-levels of $\langle y,\Gamma\rangle_{HS}-f(\Gamma)$
are bounded and hence pre-compact and allows us to conclude the existence of a
maximizer $\bar{\Gamma}$. Moreover, it is straightforward to show that
$\bar{\Gamma}>0$, otherwise one would have a contradiction by perturbing
$\bar{\Gamma}$ with $\Pi_{\ker\bar{\Gamma}}$ (the projector onto
$\ker\bar{\Gamma}$).

Let us derive the optimality conditions for $\bar{\Gamma}$. Define
$\Gamma_{s}:=\bar{\Gamma}+s\Gamma^{\prime}$ with
$\Gamma^{\prime}\in\mathcal{S}^{\bm{d}}$ (note that for any
$\Gamma^{\prime}\in\mathcal{S}^{\bm{d}}$ for $s$ sufficiently small
$\Gamma_{s}$ is positive since $\bar{\Gamma}>0$), then the Euler-Lagrange
equation for the maximization problem reads

$\displaystyle 0=\frac{d}{ds}\Big{|}_{s=0}\left(\langle
Y,\Gamma_{s}\rangle_{HS}-S(\Gamma_{s})\right)=\langle
Y,\Gamma^{\prime}\rangle_{HS}-\operatorname{Tr}\left[\Gamma^{\prime}(\log\bar{\Gamma}+1)\right].$

This yields $\bar{\Gamma}=\exp(Y-1)$. Substituting in the expression for
$S^{*}$, we arrive at (4.5).

Plugging this into (4.4) with $Y=\overline{Y}$ and recalling the definition of
$Y$, we obtain

$\displaystyle\operatorname{F}^{\varepsilon}(\Gamma)\geq\sum_{j=1}^{N}\operatorname{Tr}(U_{j}\gamma_{j})-\varepsilon\operatorname{Tr}\left(\exp\left(\frac{\bigoplus_{j=1}^{N}U_{j}-\operatorname{H}-\varepsilon}{\varepsilon}\right)\right).$

Changing the variable $U_{1}$ to $\tilde{U}_{1}:=U_{1}+\varepsilon$, we
conclude the proof. ∎

######  Remark 4.2 (The non-commutative Schrödinger problem).

Suppose that $\bm{U}\in{}_{i=1}^{N}\mathcal{S}^{d_{i}}$ is a fixed point for
$\tau$, namely $\tau(\bm{U})=\bm{U}$. This can be equivalently recast as
$\mathfrak{T}_{i}^{\varepsilon}(\hat{\bm{U}}_{i})=U_{i}$, $\forall i\in[N].$
Then Lemma 3.3, (3.25) imply that the density matrix defined by

$\displaystyle\Gamma:=\exp\left(\frac{\bigoplus_{i=1}^{N}U_{i}-\operatorname{H}}{\varepsilon}\right)$
(4.6)

has the correct marginals $\Gamma\mapsto(\gamma_{1},\dots,\gamma_{N})$ and
thus it is admissible for the primal problem. In particular, it has trace $1$
and we have

$\displaystyle\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(U_{1},\dots,U_{N})=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U})=\sum_{i=1}^{N}\operatorname{Tr}(U_{i}\gamma_{i})=\operatorname{Tr}\left(\left(\bigoplus_{i=1}^{N}U_{i}\right)\Gamma\right).$

On the other hand, directly from formula (4.6), we compute

$\displaystyle\Gamma\operatorname{H}+\varepsilon\,\Gamma\log\Gamma=\Gamma\operatorname{H}+\Gamma\left(\bigoplus_{i=1}^{N}U_{i}-\operatorname{H}\right)=\Gamma\left(\bigoplus_{i=1}^{N}U_{i}\right)$

and thus

$\displaystyle\operatorname{F}^{\varepsilon}(\Gamma)=\operatorname{Tr}\left(\left(\bigoplus_{i=1}^{N}U_{i}\right)\Gamma\right)=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(U_{1},\dots,U_{N}).$
(4.7)

In light of Proposition 4.1, this shows that if we are able to find a fixed
point of $\tau$, then this must be optimal for the dual problem (note that any
maximizer is also a fixed point for $\tau$ as discussed in Remark 3.13) and
the corresponding $\Gamma$ as obtained in (4.6) must be optimal for the primal
problem.

Another consequence of the above observations is that the set of maximizers
for the dual problem is invariant under translations.

######  Lemma 4.1 (Structure of the maximizers).

Let $\bm{U}$ and $\bm{V}$ be two maximizers of
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$, then there exists
$\bm{\alpha}\in\mathbb{R}^{N}$ such that $\sum_{i=1}^{N}\bm{\alpha}_{i}=0$ and
$\bm{U}=\bm{V}+\bm{\alpha}.$

###### Proof.

Thanks to Remark 4.2 and using that the primal functional admits an unique
minimizer by strict convexity, we find

$\displaystyle\exp\left(\frac{\bigoplus_{i=1}^{N}(\bm{U})_{i}-\operatorname{H}}{\varepsilon}\right)=\exp\left(\frac{\bigoplus_{i=1}^{N}(\bm{V})_{i}-\operatorname{H}}{\varepsilon}\right)\;\Longrightarrow\;\bigoplus_{i=1}^{N}(\bm{U})_{i}=\bigoplus_{i=1}^{N}(\bm{V})_{i}\
.$ (4.8)

Applying the partial traces to the latter equality, we obtain

$\displaystyle(\bm{U})_{i}=(\bm{V})_{i}+\sum_{j=1,j\neq
i}^{N}\operatorname{Tr}(\bm{V})_{j}-\operatorname{Tr}(\bm{U})_{j}=:(\bm{V})_{i}+\bm{\alpha}_{i}.$

Using (4.8) once again, one sees that

$\displaystyle\sum_{i=1}^{N}\bm{\alpha}_{i}=(N-1)\left(\operatorname{Tr}\left(\bigoplus_{i=1}^{N}(\bm{U})_{i}\right)-\operatorname{Tr}\left(\bigoplus_{i=1}^{N}(\bm{V})_{i}\right)\right)=0,$

which concludes the proof. ∎

###  4.2 Proof of Theorem 2.1

We are finally ready to prove the equivalence between dual and primal problem,
and to characterise the optimisers of the two problems. For the sake of
clarity, recall that

$\displaystyle\operatorname{Ren}(\bm{U})_{i}=\begin{cases}U_{i}-\lambda_{\varepsilon}(U_{i}),&\text{
if }i\in\\{1,\dots,N-1\\}\\\ \displaystyle
U_{N}+\sum_{j=1}^{N-1}\lambda_{\varepsilon}(U_{i}),&\text{ if
}i=N,\end{cases}$

as in Definition 3.15 and $\lambda_{\varepsilon}$ is defined in (3.16) as
$\lambda_{\varepsilon}(A):=\varepsilon\log\Big{(}\operatorname{Tr}\Big{[}\exp\left(\frac{A}{\varepsilon}\right)\Big{]}\Big{)}$,
for every $A\in\mathcal{S}^{d}$, $d\in\mathbb{N}$.

###### Proof of Theorem 2.1.

$(ii)$. Take a maximizing sequence $\bm{U}_{n}$ for the dual problem and
consider $\tilde{\bm{U}}_{n}:=\operatorname{Ren}\tau(\bm{U}_{n})$, where
$\tau=\mathcal{T}^{\varepsilon}_{N}\circ\dots\circ\mathcal{T}^{\varepsilon}_{1}$
is the Sinkhorn operator as introduced in Definition 3.12. Thanks to
Proposition 3.3, $\tilde{\bm{U}}_{n}$ is again a maximizing sequence that
satisfies

$\displaystyle\big{\|}\tilde{\bm{U}}_{n}\big{\|}_{\infty}\leq
2\|\operatorname{H}\|_{\infty}+\varepsilon\sup_{i\in[N]}\|\log\gamma_{i}\|_{\infty}<\infty,\quad\forall
n\in\mathbb{N},$

and it is therefore compact. Pick any
$\bm{U}^{\varepsilon}\in{}_{i=1}^{N}\mathcal{S}^{d_{i}}$ limit point of
$\tilde{\bm{U}}_{n}$. By continuity of the dual functional we infer

$\displaystyle\mathfrak{D}^{\varepsilon}(\bm{\gamma})=\lim_{N\to\infty}\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\tilde{\bm{U}}_{n})=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}^{\varepsilon})$

which shows that $\bm{U}^{\varepsilon}$ is a maximizer for
$\mathfrak{D}^{\varepsilon}(\bm{\gamma})$. The fact that any other maximizer
must coincide with $\bm{U}^{\varepsilon}$ follows from Lemma 4.1.

$(i)\&(iii)$ Proposition 4.1 proves one of the inequalities. To show the other
inequality, we take any maximizer $\bm{U}^{\varepsilon}$ (which exists by the
previous proof of $(ii)$). By construction of the Sinkhorn map,
$\bm{U}^{\varepsilon}$ must be a fixed point of $\tau$. Thanks to Remark 4.2,
we conclude that

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

satisfies
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}^{\varepsilon})=\operatorname{F}^{\varepsilon}(\Gamma^{\varepsilon})\geq\mathfrak{F}^{\varepsilon}(\bm{\gamma})$.
Hence $\Gamma^{\varepsilon}$ is optimal for $\operatorname{F}^{\varepsilon}$
and
$\mathfrak{F}^{\varepsilon}(\bm{\gamma})=\mathfrak{D}^{\varepsilon}(\bm{\gamma})$.
∎

###  4.3 Stability and the functional derivative of
$\mathfrak{F}^{\varepsilon}$

In this last section, we show stability of the Kantorovich potentials with
respect to the marginals $\bm{\gamma}$ and compute the Fréchet differential of
$\mathfrak{F}^{\varepsilon}(\bm{\gamma})$ (or simply the differential in our
finite dimensional setting). A similar result was first obtained by Pernal in
[41] at zero temperature and in [28] in the positive temperature 1RDMFT case,
i.e. considering also the fermionic and bosonic symmetry constraints. In [41],
the result follows by a direct computation via chain rule, by taking the
partial derivatives with respect to the eigenvalues and eigenvectors of a
density matrix $\Gamma$. On the other hand, [28] uses tools from convex
analysis and exploits the regularity of $\mathfrak{F}^{\varepsilon}$.

Our strategy is based on the Kantorovich formulation of (2.1) and follows
ideas contained in [21].

###### Proof of Proposition 2.1.

Consider $\bm{\gamma}^{n}\xrightarrow[]{n\to\infty}\bm{\gamma}$ and pick any
sequence of Kantorovich potentials $\bm{U}^{\varepsilon,n}$ for
$\mathfrak{F}^{\varepsilon}(\bm{\gamma}^{n})$. By optimality, they must be a
fixed point for $\tau$ and hence, thanks to Proposition 3.3,
$\operatorname{Ren}(\bm{U}^{\varepsilon,n})$ is uniformly bounded. Note that
$\operatorname{Ren}(\bm{U}^{\varepsilon,n})$ are also maximizers for
$\mathfrak{D}^{\varepsilon}(\bm{\gamma}^{n})$. This implies that any limit
point of $\operatorname{Ren}(\bm{U}^{\varepsilon,n})$ must be a maximizer for
$\mathfrak{D}^{\varepsilon}(\bm{\gamma})$. The continuity of
$\mathfrak{F}^{\varepsilon}(\cdot)$ directly follows from this stability
property.

Let us prove the differentiability. Fix $\sigma\in\mathcal{S}^{d_{i}}$, with
$\operatorname{Tr}(\sigma)=0$, and denote by $\bm{\gamma}^{h}$ the pertubation
of $\bm{\gamma}$ with $+h\sigma$ in the $i$th entry. Denote by
$\bm{U}^{\varepsilon}$ any Kantorovich potential for
$\mathfrak{F}^{\varepsilon}(\bm{\gamma})$. From duality (Theorem 2.1) we can
estimate

$\displaystyle\frac{1}{h}\left(\mathfrak{F}^{\varepsilon}(\bm{\gamma}^{h})-\mathfrak{F}^{\varepsilon}(\bm{\gamma})\right)\geq\frac{1}{h}\left(\sum_{i=1}^{N}\operatorname{Tr}\big{(}\bm{U}_{i}^{\varepsilon}\gamma_{i}^{h}-\bm{U}_{i}^{\varepsilon}\gamma_{i}\big{)}\right)=\operatorname{Tr}(\bm{U}_{i}^{\varepsilon}\sigma)$
(4.9)

for every $h\in\mathbb{R}$. Reversely, denote by $\bm{U}^{\varepsilon,h}$ any
sequence of Kantorovich potentials for
$\mathfrak{F}^{\varepsilon}(\bm{\gamma}^{h})$. Then for every $h>0$ we obtain

$\displaystyle\frac{1}{h}\left(\mathfrak{F}^{\varepsilon}(\bm{\gamma}^{h})-\mathfrak{F}^{\varepsilon}(\bm{\gamma})\right)\leq\frac{1}{h}\left(\sum_{i=1}^{N}\operatorname{Tr}\big{(}\bm{U}_{i}^{\varepsilon,h}\gamma_{i}^{h}-\bm{U}_{i}^{\varepsilon,h}\gamma_{i}\big{)}\right)=\operatorname{Tr}(\bm{U}_{i}^{\varepsilon,h}\sigma).$
(4.10)

From the first part of the proof, we know that any limit point of
$\operatorname{Ren}(\bm{U}_{i}^{\varepsilon,h})$ is a Kantorovich potential,
which up to translation (Lemma 4.1) must coincide with
$\bm{U}_{i}^{\varepsilon}$. Therefore, passing to the limit in (4.9) and
(4.10), we obtain (2.4). ∎

##  5 Non-commutative Sinkhorn algorithm

In this section we introduce and prove convergences guarantees (Theorem 2.2)
of the non-commutative version of the Sinkhorn algorithm, allowing us to
compute numerically the minimiser (2.3) of the non-commutative multi-marginal
optimal transport problem (4.1).

The idea of the Sinkhorn algorithm is to fix the shape of an ansatz

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

since it is the actual shape of the minimizer in (2.3), and alternately
project the Kantorovich potentials $U_{i}^{(k)}$ via the
$(\operatorname{H},\varepsilon)$-transforms (Definition 3.9) to approximately
reach the constraints $\Gamma^{(k)}\mapsto(\gamma_{1},\dots,\gamma_{N})$.
Recall that for $i\in[N]$, the one-step operators
$\mathcal{T}^{\varepsilon}_{i}:{}_{i=1}^{N}\mathcal{S}^{d_{j}}\to{}_{i=1}^{N}\mathcal{S}^{d_{j}}$
are given by

$\displaystyle\bm{U}:=(U_{1},\dots,U_{N}),\quad\big{(}\mathcal{T}^{\varepsilon}_{i}(\bm{U})\big{)}_{j}=\begin{cases}U_{j}&\text{if
}j\neq i,\\\ \mathfrak{T}_{i}^{\varepsilon}(U_{1},\dots,U_{i-1},U_{i+1},\dots
U_{N})&\text{if }j=i\end{cases}$

where $\mathfrak{T}_{i}^{\varepsilon}$ can be implicitly defined (Lemma 3.3)
solving the equation

$\displaystyle\operatorname{P}_{i}\left[\exp\left(\frac{\bigoplus_{i=1}^{N}\big{(}\mathcal{T}^{\varepsilon}_{i}(\bm{U})\big{)}_{j}-\operatorname{H}}{\varepsilon}\right)\right]=\gamma.$
(5.1)

_Connection with the multi-marginal Sinkhorn algorithm:_ let us shortly
describe what is the corresponding picture in the commutative setting [20,
21]. For every $i\in[N]$, let $X_{i}$ be Polish Spaces,
$\rho_{i}\mathfrak{m}_{i}\in\mathcal{P}(X_{i})$ be probability measures with
reference measures $\mathfrak{m}_{i}$. The Hamiltonian $\operatorname{H}$
corresponds to a bounded cost function $c:X_{1}\times\dots\times
X_{N}\to\mathbb{R}$.

The Sinkhorn iterates define recursively the sequences
$(a^{n}_{j})_{n\in\mathbb{N}},j\in[N]$ by

$\begin{array}[]{lcl}\displaystyle a^{0}_{j}(x_{j})&=&\rho_{j}(x_{j}),\quad
j\in\\{2,\dots,N\\},\\\ \displaystyle
a^{n}_{j}(x_{j})&=&\dfrac{\rho_{j}(x_{j})}{\int\otimes^{N}_{i<j}a_{i}^{n}(x_{i})\otimes^{N}_{i>j}a_{i}^{n-1}(x_{i})e^{-c(x_{1},\dots,x_{N})/\varepsilon}{\rm
d}(\otimes^{N}_{i\neq j}\mathfrak{m}_{i})},\,\forall n\in\mathbb{N}\text{ and
}j\in[N].\end{array}$ (5.2)

Via the new variables $u^{n}_{j}=\varepsilon\ln(a^{n}_{j}),\,j\in[N]$, one can
rewrite the Sinkhorn sequences (5.2) as

$\displaystyle u^{n}_{j}(x_{j})$
$\displaystyle=-\varepsilon\log\left(\int_{\Pi_{i\neq
j}X_{i}}\exp\left(\frac{\sum_{i\neq
j}u^{n}_{i}(x_{i})-c(x_{1},\dots,x_{N})}{\varepsilon}\right){\rm
d}\left(\otimes^{N}_{i\neq
j}\mathfrak{m}_{i}\right)\right)+\varepsilon\log(\rho_{j})$
$\displaystyle=(\hat{u^{n}_{j}})^{(N,c,\varepsilon)}(x_{j}).$

Or, more generally, for every $j\in[N]$, $u^{n}_{j}(x_{j})$ corresponds to the
solution of the maximisation

$\mathop{\mathrm{argmax}}_{u_{i}\in
L^{\infty}(X_{i})}\left\\{\sum^{N}_{i=1}\int_{X_{j}}u_{i}\rho_{j}{\rm
d}\mathfrak{m}_{j}-\varepsilon\int_{\Pi^{N}_{i\neq
j}X_{i}}\exp\left(\frac{\sum_{i\neq j}u^{n}_{i}+u-c}{\varepsilon}\right){\rm
d}\left(\otimes^{N}_{i\neq
j}\mathfrak{m}_{i}\right)\right\\}+\varepsilon\log(\rho_{j})$

which corresponds to the commutative counterpart of the $i$-th vectorial
$(\operatorname{H},\varepsilon)$-transform in Definition 3.9.

###  5.1 Definition of the algorithm

The non-commutative Sinkhorn algorithm is then defined iterating the
$(\operatorname{H},\varepsilon)$-transforms as in (5.1) for every $i\in[N]$.
Note that, by construction, the matrix
$\exp\left(\bigoplus_{i=1}^{N}\left(\left(\mathcal{T}^{\varepsilon}_{i}(\bm{U})\right)_{j}-\operatorname{H}\right)/\varepsilon\right)\in\mathfrak{P}^{\bm{d}}$
and its $i$-th marginal coincide with $\gamma_{i}$. We define the one-step
Sinkhorn map as

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

Note that this is the non-commutative counterpart of the iteration defined in
(5.2). The Sinkhorn algorithm is obtained iterating the map $\tau$ in the
following way.

Step $0$. We fix $\bm{U}^{(0)}\in\times_{i=1}^{N}\mathcal{S}^{d_{i}}$ an
initial vector of potentials and define the density matrix

$\displaystyle\Gamma^{(0)}:=\exp\left(\frac{\bigoplus_{i=1}^{N}\bm{U}^{(0)}_{i}-\operatorname{H}}{\varepsilon}\right)\in\mathfrak{P}^{\bm{d}}.$

Step $k$. For every $k\in\mathbb{N}$, we define the $k$-th density matrix via
the formula

$\displaystyle\Gamma^{(k)}:=\exp\left(\frac{\bigoplus_{i=1}^{N}\tau^{k}(\bm{U}^{(0)})_{i}-\operatorname{H}}{\varepsilon}\right)\in\mathfrak{P}^{\bm{d}}\,,$
(5.3)

where we write $\tau^{k}:=\tau\circ\dots\circ\tau$ the composition of $\tau$
for $k$-times.

Our goal is to prove the convergence $\Gamma^{(k)}\to\Gamma^{\varepsilon}$
where $\Gamma^{\varepsilon}$ is optimal for
$\mathfrak{F}^{\varepsilon}(\bm{\gamma})$. To do so, our plan is to obtain
compactness at the level of the corresponding dual potentials. Nonetheless,
the vectors $\tau^{k}(\bm{U}^{(0)})$ do not enjoy good a priori estimates and
a renormalisation procedure is needed. For any given sequence
$(\bm{\alpha}^{k})_{k\in\mathbb{N}}\subset\mathbb{R}^{N}$ such that
$\sum_{i=1}^{N}\bm{\alpha}^{k}_{i}=0$, we define

$\displaystyle\bm{U}^{(k)}:=\tau^{k}(\bm{U}^{(0)})+\bm{\alpha}^{k},\quad
k\in\mathbb{N},$ (5.4)

and observe that, by the properties of $\bigoplus$, the correspond density
matrix does not change, thus

$\displaystyle\Gamma^{(k)}=\exp\left(\frac{\bigoplus_{i=1}^{N}\bm{U}^{(k)}_{i}-\operatorname{H}}{\varepsilon}\right)\in\mathfrak{P}^{\bm{d}},\quad\forall
k\in\mathbb{N}.$ (5.5)

Thanks to the good property of the renormalisation map and the Sinkhorn
operator, we claim we can find a sequence $\bm{\alpha}^{k}$ such that the
corresponding potentials $\bm{U}^{(k)}$ as defined in (5.4) do enjoy good a
priori estimates and they can be used to prove the convergence of the
algorithm, as we see in the next section.

###  5.2 Convergence guarantees: proof of Theorem 2.2

We are ready to prove our main result Theorem 2.2, which follows from the next
Proposition.

######  Proposition 5.1 (Convergence of non-commutative Sinkhorn algorithm).

Fix $N\in\mathbb{N}$ and $\varepsilon>0$. For all $i\in[N]$, let
$\gamma_{i}\in\mathfrak{P}^{d_{i}}$ be density matrices,
$\operatorname{H}\in\mathcal{S}^{\bm{d}}$, with $\ker\gamma_{i}=\\{0\\}$. For
any initial potential $\bm{U}^{(0)}\in{}_{i=1}^{N}\mathcal{S}^{d}$, we
consider the sequence $\Gamma^{(k)}\in\mathfrak{P}^{\bm{d}}$ as defined in
(5.3).

  1. 1.

There exist $\bm{\alpha}^{k}\in\mathbb{R}^{N}$ with
$\sum_{i=1}^{N}\alpha_{i}^{k}=0$ such that

$\displaystyle\bm{U}^{(k)}=\tau^{k}(\bm{U})+\bm{\alpha}^{k}\rightarrow\bm{U}^{\varepsilon}\quad\text{as
}k\to+\infty.$ (5.6)

  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 (4.3).

  3. 3.

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

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

###### Proof.

For any $\bm{U}^{(0)}\in{}_{i=1}^{N}\mathcal{S}^{d_{i}}$, we define the
sequence $\bm{U}_{k}:=\operatorname{Ren}\tau^{k}(\bm{U}^{(0)}).$ Note that
$\bm{U}_{k}$ is of the form (5.4), for some $\bm{\alpha}^{k}$. Thanks to
Proposition 3.3, we infer that $\bm{U}_{k}$ is uniformly bounded and hence
compact. Therefore, there exists a subsequence
$\bm{U}_{k_{j}}\to\bm{U}^{\varepsilon}$. We first show that
$\bm{U}^{\varepsilon}$ is a maximizer for the dual problem. Indeed, using the
properties of $\operatorname{Ren}$ and $\tau$, it holds

$\displaystyle\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\tau(\bm{U}_{k_{j}}))=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\tau^{k_{j}+1}(\bm{U}^{(0)}))\leq\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\tau^{k_{j+1}}(\bm{U}^{(0)}))=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}_{k_{j+1}}).$

Passing to the limit the previous inequality, using the continuity of
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ and $\tau$ and recalling that
for any $\bm{U}$ we have
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\tau(\bm{U}))\geq\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U})$,
we obtain

$\displaystyle\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\tau(\bm{U}^{\varepsilon}))=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U}^{\varepsilon}).$

By definition, this means that $\bm{U}^{\varepsilon}$ is a fixed point for
$\tau$ and therefore a maximizer (Remark 4.2).

In order to prove $(1)$, we show there exists a choice $\bm{\alpha}^{k}$ such
that $\bm{U}_{k}+\bm{\alpha}^{k}\to\bm{U}^{\varepsilon}$. For $k=k_{j}$ for
some $j$, we pick $\bm{\alpha}^{k}=0$, for all the others $k$, we instead pick
$\bm{\alpha}^{k}$ defined by

$\displaystyle\bm{\alpha}^{k}=\operatorname{argmin}_{\bm{\alpha}}\left\\{\|\bm{U}_{k}+\bm{\alpha}-\bm{U}^{\varepsilon}\|_{\infty}\
:\ \sum_{i=1}^{N}\bm{\alpha}_{i}=0\right\\}.$

Note that, by Lemma 4.1, this is equivalent to picking $\bm{\alpha}^{k}$ such
that $\bm{U}^{\varepsilon}$ is the closest maximizer to
$\bm{U}_{k}+\bm{\alpha}^{k}$. We claim this is the right choice. Suppose
indeed by contradiction that there exists a subsequence
$\bm{U}_{k^{\prime}_{j}}$ such that
$\|\bm{U}_{k^{\prime}_{j}}+\bm{\alpha}_{k_{j}}-\bm{U}^{\varepsilon}\|_{\infty}\geq\delta>0$,
then by construction
$\|\bm{U}_{k^{\prime}_{j}}+\bm{\alpha}_{k_{j}}-\bm{U}^{\prime}\|_{\infty}\geq\delta$
for any other maximizer $\bm{U}^{\prime}$. By compactness, this is a
contradiction, since there exists a further subsequence
$\bm{U}_{k^{\prime\prime}_{j}}$ of $\bm{U}_{k^{\prime}_{j}}$ converging to a
maximizer $\bm{U}^{\prime}$ (by the same reasoning carried out above). This
proves $(1)$ and by optimality of $\bm{U}^{\varepsilon}$, $(2)$ as well. The
convergence of $\Gamma^{(k)}$ follows from the compactness of $\bm{U}^{(k)}$
and (5.5), whereas the optimality of the limit point $\Gamma^{\varepsilon}$
and (5.7) are consequence of the optimality of $\bm{U}^{\varepsilon}$ and
Remark 4.2. ∎

##  6 One-body reduced density matrix functional theory

In this last section, we prove Proposition 2.1 and consequently Theorem 2.7.

For given $d,N\in\mathbb{N}$, we set $\bm{d}=d^{N}$ and consider the space of
bosonic (resp. fermionic) density matrices $\mathfrak{P}_{+}^{d}$ (resp.
$\mathfrak{P}_{-}^{d}$) as introduced in (2.9). Recall as well that 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}$, where $\Pi_{\pm}$ are defined in
(2.8).

The universal functional in the bosonic and in the fermionic case is then
given as in Definition 2.5, which we recall here for simplicity is given by

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

whereas the corresponding dual functional and problem (see Definition 2.6) are
given by

$\displaystyle\text{D}_{\gamma}^{\pm,\varepsilon}(U):=\operatorname{Tr}(U\gamma)-\varepsilon\operatorname{Tr}\left(\exp\left[\frac{1}{\varepsilon}\bigg{(}\frac{1}{N}\bigotimes_{i=1}^{N}U-\operatorname{H}\bigg{)}_{\pm}\right]\right)+\varepsilon\,,$
$\displaystyle\mathfrak{D}_{\pm}^{\varepsilon}(\gamma):=\sup\left\\{\text{D}_{\gamma}^{\pm,\varepsilon}(U)\
:\ U\in\mathcal{S}^{d}\right\\}\,.$

We are interested in fully characterizing the existence of the optimizers in
the primal and the dual problems, for both bosonic and fermionic cases.
Proceeding in a similar way as in the proof of Lemma 3.2, one can prove that
every maximizer $U_{\pm}^{\varepsilon}$ of the dual functional
$\text{D}_{\gamma}^{\pm,\varepsilon}(\cdot)$ must satisfy the corresponding
Euler-Lagrange equation given by

$\displaystyle\gamma=\operatorname{P}_{1}\left(\exp\left[\frac{1}{\varepsilon}\bigg{(}\frac{1}{N}\bigoplus_{i=1}^{N}U_{\pm}^{\varepsilon}-\operatorname{H}\bigg{)}_{\pm}\right]\right)\,.$
(6.1)

###  6.1 Fermionic dual problem and Pauli’s exclusion principle

The aim of this section is to prove Proposition 2.1. For simplicity we assume,
with no loss of generality, that $\varepsilon=1$ and set
$\text{D}_{\gamma}^{-}:=\operatorname{D}_{\gamma}^{-,1}$.

For any $U\in\mathcal{S}^{d}$, we fix a basis of normalized eigenvectors of
$U$, denoted by $\\{\psi_{j}\\}_{j}$, and consider the decomposition

$\displaystyle U=\sum_{j=1}^{d}u_{j}|\psi_{j}\rangle\langle\psi_{j}|,\quad
u_{j}\in\mathbb{R}\quad(\text{eigenvalues})\,.$ (6.2)

We also denote by $\gamma_{j}:=\langle\psi_{j}|\gamma|\psi_{j}\rangle$. In
particular, the linear terms read

$\displaystyle\operatorname{Tr}(U\gamma)=\sum_{j=1}^{d}\gamma_{j}u_{j}\,.$

For any such basis $\\{\psi_{i}\\}_{i}$, we obtain a basis of the fermionic
tensor product

$\displaystyle\psi^{\emph{as}}_{\bm{j}}:=\bigwedge_{i=1}^{N}\psi_{j_{i}}\,,\quad\bm{j}=(j_{i})_{i=1}^{N}\in\Theta_{-}\,,$
$\displaystyle\Theta_{-}:=\big{\\{}(j_{1},\dots,j_{N})\ :\
j_{i}\in\\{1,\dots,d\\},\;j_{i}\neq j_{k},\,\text{if }i\neq
k\big{\\}}/\mathfrak{S}_{N}\,,$

where $\mathfrak{S}_{N}$ denotes the set of permutations of $N$ elements. With
respect to this basis, we can write

$\displaystyle\frac{1}{N}\left(\bigoplus_{i=1}^{N}U\right)_{-}=\sum_{\bm{j}\in\Theta_{-}}\bigg{(}\frac{1}{N}\sum_{i=1}^{N}u_{j_{i}}\bigg{)}|\psi_{\bm{j}}^{\emph{as}}\rangle\langle\psi_{\bm{j}}^{\emph{as}}|\,.$
(6.3)

Using the monotonicity of the exponential and the trace, we obtain the
following result.

######  Lemma 6.1 (Bounds for $\text{D}_{\gamma}^{-}(U)$).

Fix $U\in\mathcal{S}^{d}$ with eigenvalues $u_{j}$ and eigenvectors
$\\{\psi_{j}\\}_{j}$. For $\gamma\in\mathfrak{P}(d)$, set
$\gamma_{j}:=\langle\psi_{j}|\gamma|\psi_{j}\rangle$. Then one has

$\displaystyle\begin{aligned}
\sum_{j=1}^{d}\gamma_{j}u_{j}-C\sum_{\bm{j}\in\Theta_{-}}\exp\left(\frac{1}{N}\sum_{i=1}^{N}u_{j_{i}}\right)&\leq\emph{D}_{\gamma}^{-}(U)-1\\\
&\leq\sum_{j=1}^{d}\gamma_{j}u_{j}-\frac{1}{C}\sum_{\bm{j}\in\Theta_{-}}\exp\left(\frac{1}{N}\sum_{i=1}^{N}u_{j_{i}}\right),\end{aligned}$
(6.4)

where $C=\exp\big{(}\|\operatorname{H}\|_{\infty}\big{)}\in(0,+\infty)$.

Before moving to the proof of Proposition 2.1, we need the following technical
lemma.

######  Lemma 6.2 (Linear term estimates).

Consider $\\{u_{j}\\}_{j=1}^{d}\subset\mathbb{R}$ and
$\\{\gamma_{j}\\}_{j=1}^{d}$ such that

$\displaystyle\gamma_{j}\in\Big{(}\delta,\frac{1}{N}-\delta\Big{)},\quad\sum_{j=1}^{d}\gamma_{j}=1\,,$
(6.5)

for some $\displaystyle\delta\in\Big{[}0,\frac{1}{2N}\Big{)}$. Suppose that
$u_{j}\leq u_{k}$ if $j\leq k$. Then we have

$\displaystyle\sum_{j=1}^{d}\gamma_{j}u_{j}\leq\frac{1}{N}\sum_{i=1}^{N}u_{j}-\delta(u_{1}-u_{d})\,.$
(6.6)

###### Proof.

Thanks to the fact the $u_{j}$ are ordered, we have the inequality

$\displaystyle\sum_{j=1}^{d}\bar{\gamma}_{j}u_{j}\leq\frac{1}{N}\sum_{i=1}^{N}u_{j},\quad\forall\;0\leq\bar{\gamma}_{j}\leq\frac{1}{N},\quad\sum_{j=1}^{d}\bar{\gamma}_{j}=1\,.$

Then (6.6) follows applying the above inequality to

$\displaystyle\bar{\gamma}_{1}:=\gamma_{1}+\delta\in\Big{(}0,\frac{1}{N}\Big{)},\quad\bar{\gamma}_{j}:=\gamma_{j},\quad\bar{\gamma}_{d}:=\gamma_{d}-\delta\in\Big{(}0,\frac{1}{N}\Big{)}\,,$

for every $j\in\\{2,\dots,d-1\\}$. ∎

We are ready to prove Proposition 2.1.

###### Proof.

($\gamma\leq 1/N\Rightarrow$ $\sup\text{D}_{\gamma}^{-}<\infty$). This is
consequence of Proposition 6.2 with $\delta=0$. More precisely, pick
$U\in\mathcal{S}^{d}$ and consider a decomposition in eigenfunctions as in
(6.2). Assume that $\\{u_{j}\\}_{j}$ are non increasing in $j$ (with no loss
of generality). We can then apply Proposition 6.2 with $\delta=0$ and from
(6.6) and (6.4) we deduce

$\displaystyle\text{D}_{\gamma}^{-}(U)-1\leq\frac{1}{N}\sum_{i=1}^{N}u_{j}-\frac{1}{C}\sum_{\bm{j}\in\Theta_{-}}\exp\left(\frac{1}{N}\sum_{i=1}^{N}u_{j_{i}}\right)\leq\frac{1}{N}\sum_{i=1}^{N}u_{j}-\frac{1}{C}\exp\left(\frac{1}{N}\sum_{i=1}^{N}u_{j}\right),$

where in the last inequality we used the positivity of the exponential.
Therefore

$\displaystyle\sup_{U\in\mathcal{S}^{d}}\text{D}_{\gamma}^{-}(U)\leq\sup_{x\in\mathbb{R}}\,(x-\frac{1}{C}e^{x})+1=\log
C<\infty\,.$

($\sup\text{D}_{\gamma}^{-}<\infty$ $\Rightarrow\gamma\leq 1/N$). Suppose by
contradiction that the Pauli’s principle is not satisfied. With no loss of
generality, we can assume that

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

We build the sequence of bounded operators $U^{n}\in\mathcal{S}^{d}$ given by

$\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\,.$ (6.7)

Observe that by construction, we can estimate the non-linear part of
$\text{D}_{\gamma}^{-}(U)$ as

$\displaystyle\forall\bm{j}\in\Theta_{-},\quad\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}$

It follows that we can bound from below $\text{D}_{\gamma}^{-}(U^{n})$ as

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

We claim that the linear contribution goes to $+\infty$ as $n\to+\infty$. To
see that, note that

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

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

(Equation for the maximizer and uniqueness). If a maximizer exists, then it
solves the equation (6.1). Thanks to the Peierls inequality, we also know that
$\text{D}_{\gamma}^{-}$ is strictly concave (because the exponential is
strictly convex), hence the uniqueness of the maximizer.

(Existence of $\mathop{\mathrm{argmax}}\text{D}_{\gamma}^{-}$ $\Rightarrow
0<\gamma<1/N$ ). We proceed as in the latter proof. By contradiction, assume
that

$\displaystyle\gamma=\sum_{j=1}^{d}\gamma_{j}|\psi_{j}\rangle\langle\psi_{j}|,\quad\gamma_{1}=\frac{1}{N}\,,\quad\gamma_{j}\in\Big{(}0,\frac{1}{N}\Big{)}\,,\quad\forall
j\geq 2\,.$

The case $\gamma_{j}=0$ can be directly ruled out from the Euler-Lagrange
equation for a maximizer (6.1). We can then consider the very same sequence
$U^{n}$ as defined in (6.7). From (6.8), (6.9), and the first part of Theorem
2.1, on one hand we deduce

$\displaystyle-C\binom{d}{N}\leq\text{D}_{\gamma}^{-}(U^{n})\leq\log
C\,,\quad\forall n\in\mathbb{N}\,.$

On the other hand, $\|U^{n}\|_{\infty}\to+\infty$ as $n\to\infty$, which means
that $\text{D}_{\gamma}^{-}$ is not coercive. Thanks to Peierls inequality, we
also know that $\text{D}_{\gamma}^{-}$ is strictly concave, which implies that
$\text{D}_{\gamma}^{-}$ can not attain its maximum.

($0<\gamma<1/N\Rightarrow$ existence of
$\mathop{\mathrm{argmax}}\text{D}_{\gamma}^{-}$). Let $U\in\mathcal{S}^{d}$
and consider a decomposition in eigenfunctions as in (6.2). Assume that
$\\{u_{j}\\}_{j}$ are non increasing in $j$ (with no loss of generality) and
denote by $\gamma_{j}:=\langle\psi_{j}|\gamma|\psi_{j}\rangle$. By assumption,
there exists $\delta\in\big{(}0,\frac{1}{N}\big{)}$ such that

$\displaystyle\sum_{j=1}^{d}\gamma_{j}=1,\quad\gamma_{j}\in\Big{(}\delta,\frac{1}{N}-\delta\Big{)}\,,\quad\forall
j\in\\{1,\dots,d\\}\,.$ (6.10)

We can then apply Proposition 6.2 and (6.4) to obtain

$\displaystyle\text{D}_{\gamma}^{-}(U)-1$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}u_{j}-\frac{1}{C}\sum_{\bm{j}\in\Theta_{-}}\exp\left(\frac{1}{N}\sum_{i=1}^{N}u_{j_{i}}\right)-\delta(u_{1}-u_{d})$
(6.11)
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}u_{j}-\frac{1}{C}\exp\left(\frac{1}{N}\sum_{i=1}^{N}u_{j}\right)-\delta(u_{1}-u_{d}),$
(6.12)

where we used the positivity of the exponential. Set
$S:=\sup_{x}(x-\frac{e^{x}}{C})+1<\infty$, and infer

$\displaystyle\text{D}_{\gamma}^{-}(U)\leq
S-\delta(u_{max}-u_{min}),\quad\forall U\in\mathcal{S}^{d},\quad
U=\sum_{j=1}^{d}u_{j}|\psi_{j}\rangle\langle\psi_{j}|\,,$ (6.13)

where $u_{max}$ and $u_{min}$ denotes respectively the maximum/minimum
eigenvalue of $U$. Let us use this estimate to prove to existence of a
maximizer for $\text{D}_{\gamma}^{-}$. Consider a maximizing sequence $U^{n}$
of bounded operators. In particular, we can assume that
$-I:=\inf_{n}\text{D}_{\gamma}^{-}(U^{n})\geq-\infty$. If the sequence
$\\{U^{n}\\}_{n}$ is bounded in $\mathcal{S}^{d}$, then any limit point is a
maximum for $\text{D}_{\gamma}^{-}$, by concavity and continuity of
$\text{D}_{\gamma}^{-}$, and the proof is complete. Suppose by contradiction
that $\|U^{n}\|_{\infty}\to+\infty$ as $n\to+\infty$. Note that from (6.13) we
deduce

$\displaystyle\sup_{n\in\mathbb{N}}\big{(}u_{max}^{n}-u_{min}^{n}\big{)}\leq\frac{S+I}{\delta}<\infty\,,$
(6.14)

therefore we deduce that either $u_{j}^{n}\to-\infty$ or $u_{j}^{n}\to+\infty$
for every $j\in\\{1,\dots,d\\}$. In the first case, we would have a
contradiction, because

$\displaystyle-I\leq\text{D}_{\gamma}^{-}(U^{n})\leq\sum_{j=1}^{d}\gamma_{j}u_{j}^{n}+1\to-\infty\quad\text{as
}n\to+\infty\,.$

In the second case, we can use (6.11) to find a contradiction, because

$\displaystyle-I\leq\text{D}_{\gamma}^{-}(U^{n})\leq\frac{1}{N}\sum_{i=1}^{N}u_{j}^{n}-\frac{1}{C}\exp\left(\frac{1}{N}\sum_{i=1}^{N}u_{j}^{n}\right)+1\to-\infty\,,$

where we used that $\displaystyle\lim_{x\to+\infty}(x-C^{-1}e^{x})=-\infty$.
The proof is complete. ∎

###  6.2 Duality theorem for fermionic and bosonic systems

In this section we prove Theorem 2.7. The proof relies on the use of Theorem
2.1 and the existence of maximizers for $\text{D}_{\gamma}^{-,\varepsilon}$,
proved in Proposition 2.1, and $\text{D}_{\gamma}^{+,\varepsilon}$. The latter
can be proven easily by noting that the spectrum of
$\left(\bigoplus_{j=1}^{N}U_{j}\right)_{+}$ contains the spectrum of $U$ and,
applying similar computations to the ones used in the case of
$\text{D}_{\gamma}^{-,\varepsilon}$, deducing the coercivity of
$\text{D}_{\gamma}^{+,\varepsilon}$. We also need the following observation.

######  Remark 6.1.

If $\operatorname{H}$ satisfies (2.10) and $\bm{\gamma}=(\gamma_{i})_{i}$ ,
$\gamma_{i}=\gamma$, then the minimizers of
$\operatorname{D}_{\bm{\gamma}}^{\varepsilon}$ (the dual functional without
symmetry constraints) can be taken to satisfy $U_{i}\equiv U$, for some
$U\in\mathcal{S}^{d}$. In particular

$\displaystyle\mathfrak{D}^{\varepsilon}(\bm{\gamma})=\sup_{\bm{U}\in(\mathcal{S}^{\bm{d}})^{N}}\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U})=\sup_{U\in\mathcal{S}^{d}}\left\\{\operatorname{Tr}(U\gamma)-\varepsilon\operatorname{Tr}\left(\exp\bigg{[}\frac{1}{\varepsilon}\bigg{(}\frac{1}{N}\bigotimes_{i=1}^{N}U-\operatorname{H}\bigg{)}\bigg{]}\right)\right\\}+\varepsilon\,.$

This follows from the observation that if
$\bm{U}\in(\mathcal{S}^{\bm{d}})^{N}$, then we obtain a symmetric competitor
$\tilde{\bm{U}}$

$\displaystyle(\tilde{\bm{U}})_{i}=\frac{1}{N}\sum_{j=1}^{N}U_{j}\,,\quad\text{such
that}\quad\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\tilde{\bm{U}})=\operatorname{D}_{\bm{\gamma}}^{\varepsilon}(\bm{U})\,.$

###### Proof of Theorem 2.7.

Let us assume that $\gamma>0$ in the bosonic case ($0<\gamma<\frac{1}{N}$ in
the fermionic case). The general duality result (including the case $\gamma$
in which does not satisfy the above strict inequalities) can be handled by
decomposition of the space, in the same way as in Remark 3.6.

Under these assumptions, thanks to Proposition 2.1, we know that a maximizer
$U_{\pm}^{\varepsilon}$ exists and satisfies (6.1).

We then define the $N$-particle density matrix

$\displaystyle\tilde{\Gamma}_{\pm}^{\varepsilon}:=\exp\bigg{[}\frac{1}{\varepsilon}\bigg{(}\frac{1}{N}\bigoplus_{i=1}^{N}U_{\pm}^{\varepsilon}-\operatorname{H}\bigg{)}_{\pm}\bigg{]}\in\mathcal{S}^{\bm{d}}\,,$

and thanks to Remark 4.2, we know that $\tilde{\Gamma}_{\pm}^{\varepsilon}$ is
optimal for the problem
$\mathfrak{F}^{\varepsilon}\big{(}\operatorname{P}_{1}(\tilde{\Gamma}_{\pm}^{\varepsilon})\big{)}$
without symmetry contraints. Observing that
$(\tilde{\Gamma}_{\pm}^{\varepsilon})_{\pm}=\Gamma_{\pm}^{\varepsilon}$
(defined in (2.15),(2.16)), we deduce that $\Gamma_{\pm}^{\varepsilon}$ must
be optimal for the primal problem $\mathfrak{F}^{\varepsilon}_{\pm}(\gamma)$
with symmetry constraints. This also proves the equality between primal and
dual problems and concludes the proof. ∎

## Acknowledgments

This work started when A.G. was visiting the Erwin Schrödinger Institute and
then continued when D.F. and L.P visited the Theoretical Chemistry Department
of the Vrije Universiteit Amsterdam. The authors thanks the hospitality of
both places and, especially, P. Gori-Giorgi and K. Giesbertz for fruitful
discussions and literature suggestions in the early state of the project.
Finally, the authors also thanks J. Maas and R. Seiringer for their feedback
and useful comments to a first draft of the article.

L.P. acknowledges support by the Austrian Science Fund (FWF), grants No W1245
and No F65. D.F acknowledges support by the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreements No 716117 and No 694227). A.G. acknowledges funding by the
European Research Council under H2020/MSCA-IF “OTmeetsDFT” [grant ID: 795942].

## References

  * [1] P. W. Ayers and M. Levy, Generalized density-functional theory: Conquering the n-representability problem with exact functionals for the electron pair density and the second-order reduced density matrix, Journal of Chemical Sciences, 117 (2005), pp. 507–514.
  * [2] T. Baldsiefen, A. Cangi, and E. K. Gross, Reduced-density-matrix-functional theory at finite temperature: Theoretical foundations, Physical Review A, 92 (2015), p. 052514.
  * [3] T. Baldsiefen, F. G. Eich, and E. K. Gross, Reduced density matrix functional theory at finite temperature. ii. application to the electron gas: Exchange only, arXiv preprint arXiv:1208.4705, (2012).
  * [4] T. Baldsiefen and E. K. Gross, Reduced density matrix functional theory at finite temperature. iii. application to the electron gas: Correlation effects, arXiv preprint arXiv:1208.4707, (2012).
  * [5] C. L. Benavides-Riveros, J. Wolff, M. A. L. Marques, and C. Schilling, Reduced density matrix functional theory for bosons, Physical Review Letters, 124 (2020), p. 180603.
  * [6] Y. Brenier and D. Vorotnikov, On optimal transport of matrix-valued measures, SIAM Journal on Mathematical Analysis, 52 (2020), pp. 2849–2873.
  * [7] M. Buijse and E. Baerends, An approximate exchange-correlation hole density as a functional of the natural orbitals, Molecular Physics, 100 (2002), pp. 401–421.
  * [8] E. Caglioti, F. Golse, and T. Paul, Towards optimal transport for quantum densities, Available on https://hal.archives-ouvertes.fr/hal-01963667/, (2018).
  * [9] , Quantum optimal transport is cheaper, Journal of Statistical Physics (2020), 181 (2020), pp. 149–162.
  * [10] E. Carlen, Trace inequalities and quantum entropy: an introductory course, Entropy and the quantum, 529 (2010), pp. 73–140.
  * [11] E. A. Carlen and J. 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 (2014), pp. 887–926.
  * [12] Y. Chen, W. Gangbo, T. Georgiou, and A. Tanenbaum, On the matrix Monge–Kantorovich problem, European Journal of Applied Mathematics, (2019), pp. 1–27.
  * [13] Y. Chen, T. Georgiou, and M. Pavon, Entropic and displacement interpolation: a computational approach using the Hilbert metric, SIAM Journal on Applied Mathematics, 76 (2016), pp. 2375–2396.
  * [14] Y. Chen, T. T. Georgiou, and A. Tannenbaum, Matrix optimal mass transport: a quantum mechanical approach, IEEE Transactions on Automatic Control, 63 (2018), pp. 2612–2619.
  * [15] L. Chizat, G. Peyré, B. Schmitzer, and F.-X. Vialard, Scaling algorithms for unbalanced optimal transport problems, Mathematics of Computation, 87 (2018), pp. 2563–2609.
  * [16] A. J. Coleman, Structure of fermion density matrices, Reviews of modern Physics, 35 (1963), p. 668.
  * [17] M. Cuturi, Sinkhorn distances: Lightspeed computation of optimal transport, in Advances in Neural Information Processing Systems, 2013, pp. 2292–2300.
  * [18] G. De Palma and D. Trevisan, Quantum optimal transport with quantum channels, arXiv preprint arXiv:1911.00803, (2019).
  * [19] G. De Palma, D. Trevisan, V. Giovannetti, and L. Ambrosio, Gaussian optimizers for entropic inequalities in quantum information, Journal of Mathematical Physics, 59 (2018), p. 081101.
  * [20] S. Di Marino and A. Gerolin, An Optimal Transport approach for the Schrödinger bridge problem and convergence of Sinkhorn algorithm, Journal of Scientific Computing, 85 (2020).
  * [21] , Optimal Transport losses and Sinkhorn algorithm with general convex regularization, arXiv preprint arXiv:2007.00976, (2020).
  * [22] R. L. Frank, E. H. Lieb, R. Seiringer, and H. Siedentop, Müller’s exchange-correlation energy in density-matrix-functional theory, Physical Review A, 76 (2007), p. 052517.
  * [23] R. L. Frank, P. T. Nam, and H. Van Den Bosch, The maximal excess charge in Müller density-matrix-functional theory, in Annales Henri Poincaré, vol. 19, Springer, 2018, pp. 2839–2867.
  * [24] J. Franklin and J. Lorenz, On the scaling of multidimensional matrices, Linear Algebra and its applications, 114 (1989), pp. 717–735.
  * [25] A. Galichon and B. Salanié, Matching with trade-offs: Revealed preferences over competing characteristics, (2010).
  * [26] T. T. Georgiou and M. Pavon, Positive contraction mappings for classical and quantum Schrödinger systems, Journal of Mathematical Physics, 56 (2015), p. 033301.
  * [27] A. Gerolin, A. Kausamo, and T. Rajala, Multi-marginal Entropy-Transport with repulsive cost, Calc. Var., 59 (2020).
  * [28] K. J. Giesbertz and M. Ruggenthaler, One-body reduced density-matrix functional theory in finite basis sets at elevated temperatures, Physics Reports, 806 (2019), pp. 1–47.
  * [29] T. L. Gilbert, Hohenberg-Kohn theorem for nonlocal external potentials, Physical Review B, 12 (1975), p. 2111.
  * [30] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev., 136 (1964).
  * [31] J. Karlsson and A. Ringh, Generalized Sinkhorn iterations for regularizing inverse problems using optimal mass transport, SIAM Journal on Imaging Sciences, 10 (2017), pp. 1935–1962.
  * [32] C. Léonard, A survey of the Schrödinger problem and some of its connections with optimal transport, Discrete & Continuous Dynamical Systems-A, 34 (2014), pp. 1533–1574.
  * [33] M. 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 (76), pp. 6062–6065.
  * [34] E. H. Lieb, Density functionals for Coulomb systems, in Inequalities, Springer, 2002, pp. 269–303.
  * [35] E. H. Lieb and R. Seiringer, The stability of matter in quantum mechanics, Cambridge University Press, Cambridge, 2010.
  * [36] L. Mentel, Reduced Density Matrix inspired approaches to electronic structure theory, PhD thesis, Vrije Universiteit Amsterdam:, 2015.
  * [37] N. D. Mermin, Thermal properties of the inhomogeneous electron gas, Physical Review, 137 (1965), p. A1441.
  * [38] M. Mittnenzweig and A. Mielke, An entropic gradient structure for Lindblad equations and couplings of quantum systems to macroscopic models, Journal of Statistical Physics, 167 (2017), pp. 205–233.
  * [39] L. Monsaingeon and D. Vorotnikov, The Schrödinger problem on the non-commutative Fisher-Rao space, arXiv preprint arXiv:2007.09042, (2020).
  * [40] A. Müller, Explicit approximate relation between reduced two-and one-particle density matrices, Physics Letters A, 105 (1984), pp. 446–452.
  * [41] K. Pernal, Effective potential for natural spin orbitals, Physical review letters, 94 (2005), p. 233002.
  * [42] K. Pernal and K. J. Giesbertz, Reduced density matrix functional theory (RDMFT) and linear response time-dependent rdmft (TD-RDMFT), in Density-Functional Methods for Excited States, Springer, 2015, pp. 125–183.
  * [43] G. Peyré, L. Chizat, F.-X. Vialard, and J. Solomon, Quantum entropic regularization of matrix-valued optimal transport, European Journal of Applied Mathematics, 30 (2019), pp. 1079–1102.
  * [44] R. Requist and O. Pankratov, Generalized Kohn-Sham system in one-matrix functional theory, Physical Review B, 77 (2008), p. 235121.
  * [45] L. Ruschendorf, Convergence of the iterative proportional fitting procedure, The Annals of Statistics, 23 (1995), pp. 1160–1174.
  * [46] E. K. Ryu, Y. Chen, W. Li, and S. Osher, Vector and matrix optimal mass transport: theory, algorithm, and applications, SIAM Journal on Scientific Computing, 40 (2018), pp. A3675–A3698.
  * [47] R. Schade, New methods for the ab-initio simulation of correlated systems, PhD thesis, Niedersächsische Staats-und Universitätsbibliothek Göttingen, 2019.
  * [48] E. Schrödinger, Über die umkehrung der naturgesetze, Verlag Akademie der wissenschaften in kommission bei Walter de Gruyter u. Company, 1931\.
  * [49] R. Sinkhorn, A relationship between arbitrary positive matrices and doubly stochastic matrices, The annals of mathematical statistics, 35 (1964), pp. 876–879.
  * [50] R. van Leeuwen, Functionals of the one-particle density matrix: Uniqueness theorems and v-representability in equilibrium and time-dependent systems, University of Jyväskylä, Jyväskylä, Finland, private communication (unpublished), (2007).
  * [51] G. Vignale and M. Rasolt, Density-functional theory in strong magnetic fields, Physical review letters, 59 (1987), p. 2360.
  * [52] U. von Barth and L. Hedin, A local exchange-correlation potential for the spin polarized case. i, Journal of Physics C: Solid State Physics, 5 (1972), p. 1629.
  * [53] M. Wirth, A dual formula for the noncommutative transport distance, arXiv preprint arXiv:2104.11923, (2021).
  * [54] J. Zambrini, Stochastic mechanics according to E. Schrödinger, Physical review A, 33 (1986), p. 1532.