text:
\section{Introduction} \label{sec:intro}
Let us recall the classical discrete optimal transport problem as stated in Hitchcock \cite{Hit41} and Kantorovich \cite{Kan60} (prepared in 1939), which is a variation of the classical transport problem initiated by Monge \cite{Mon81}.  Suppose we have $m$ factories producing $G$ amount of the same product that has to be distributed to $n$ consumers.  Assume that $x_{ij}^{AB}$ is the proportion of the goods sent from the factory $i$ to consumer $j$.
Then $x_i^A$ and $x_j^B$ are the proportions of the goods produced by factory $i$ and received by consumer $j$ respectively:
\begin{equation}\label{ccoupl}
x_i^A=\sum_{j=1}^n x_{ij}^{AB},\; i\in[m],\quad x_j^B=\sum_{i=1}^m x_{ij}^{AB}, \; j\in[n],
\end{equation}
where $[m]=\{1,2,\ldots,m\}$.
It is convenient to introduce the random variables $X^A,X^B$  such that
\begin{equation*}
x_i^A=\mathbb{P}(X^A=i),\;i\in[m], \quad x_j^B=\mathbb{P}(X^B=j),\;j\in[n].
\end{equation*}
Then the nonnegative matrix $X^{AB}=[x_{ij}^{AB}]\in\mathbb{R}_+^{m\times n}$ satisfying the above equalities is a joint distribution of the random variable $X^{AB}$: $x_{ij}^A=\mathbb{P}\big(X^{AB}=(i,j)\big)$.  The random variable $X^{AB}$, or the matrix $X^{AB}$, is called a coupling of $X^A$ and $X^B$.  Let $\mathbf{x}^A=(x_1^A,\ldots,x_m^A)^\top, \mathbf{x}^B=(x_1^B,\ldots,x_n^B)^\top$ be the probability vectors corresponding to $X^A$ and $X^B$ respectively.
The set of all coupling matrices $X^{AB}$  corresponding to $\mathbf{x}^A, \mathbf{x}^B$ is denoted by $\Gamma^{cl}(\mathbf{x}^A, \mathbf{x}^B)$.  Note that $X=\mathbf{x}^A(\mathbf{x}^B)^\top$, corresponding to independent coupling of $X^A$ and $X^B$, is in $\Gamma^{cl}(\mathbf{x}^A, \mathbf{x}^B)$.  Let $C=[c_{ij}]\in \mathbb{R}^{m\times n}_+$ be a nonnegative matrix where $c_{ij}$ is the transport cost of a unit of goods from the factory $i$ to the consumer $j$.  The classical optimal transport problem, abbreviated as OT, is
\begin{equation}\label{COTP}
\mathrm{T}_{C}^{cl}(\mathbf{x}^A,\mathbf{x}^B)=\min_{ X\in \Gamma^{cl}(\mathbf{x}^A,\mathbf{x}^B)}\tr C X ^\top.
\end{equation}
(Here $\tr $ denotes the trace of a square matrix, and $X^\top$ the transpose of $X$.)
The optimal transport problem is a linear programming problem (LP) which can be solved in polynomial time in the size of the inputs $\mathbf{x}^A,\mathbf{x}^B, C$ \cite{CCPS}.

Assume now that $m=n$.  Let $C=[c_{ij}]\in \mathbb{R}^{n\times n}_+$ be a symmetric nonnegative matrix with zero diagonal an positive off-diagonal entries such that $c_{ij}$ induces a distance on $[n]$: dist$(i,j)=c_{ij}$.  That is, in addition to the above conditions one has the triangle inequality $c_{ij}\le c_{ik}+c_{kj}$ for $i,j,k\in[n]$.
For $p>0 $ denote $C^{\circ p}=[c_{ij}^p]\in\mathbb{R}_+^{n\times n}$.
Then the quantity
\begin{equation}\label{Wasserp}
W_{C,p}^{cl}(\mathbf{x}^A,\mathbf{x}^B)=\big(\mathrm{T}_{C^{\circ p}}^{cl}(\mathbf{x}^A,\mathbf{x}^b)\big)^{1/p}, \quad p\ge 1
\end{equation}
is the Wasserstein-$p$ metric on the simplex of probability vectors, $\Pi_n\subset \mathbb{R}_+^n$. This follows from the continuous version of the Wasserstein-$p$ metric, as in \cite{Vas69}.  See \cite{Cut13} for $p=1$. It turns out that  $\mathrm{T}_{C}^{cl}(\mathbf{x}^A,\mathbf{x}^B)$ has many recent applications in machine learning \cite{AWR17,ACB17,LG15,MJ15,SL11}, statistics  \cite{BGKL17,FCCR18,PZ16,SR04} and computer vision \cite{BPPH11,RTG00,SGPCBNDG}.

Several attempts to generalize the notion of the Monge--Kantorovich
distance in quantum information theory (QIT) are known.
An early contribution defines the distance between any
 two quantum states  by the Monge distance between the corresponding Husimi
 functions \cite{ZS98,ZS01}.
As this approach depends on the choice of the set of coherent states,
other efforts were undertaken  \cite{AF17,GP18}
to introduce  the transport distance between quantum states
by applying the Kantorovich--Wasserstein optimization
over the set of bipartite quantum states with fixed marginals.
 Even though the matrix transport problem
 was often investigated in the recent literature
\cite{BV01,BGJ19,FGZ19,CGGT17,Fri20,Duv20,FECZ},
related to potential applications in quantum physics
\cite{CM20,DR20,CGP20,KdPMLL21},
 this aim has not been fully achieved until now \cite{Rie18,YZYY19,Ikeda20}.

The aim of this work is to present a constructive solution
of the  optimal transport problem  in the quantum finite-dimensional setting.
Furthermore, we show that the square root of the optimal transport cost
satisfies the triangle inequality and
 construct  a transport distance between arbitrary quantum states.

  Denote by $\Omega_m$ the convex set of density matrices, i.e., the set of $m\times m$ Hermitian positive semidefinite matrices of trace one.
Let $\rho^{A}\in\Omega_m, \rho^B\in \Omega_n$.  A quantum coupling of $\rho^A,\rho^B$ is a density matrix $\rho^{AB}\in\Omega_{mn}$, whose partial traces give $\rho^A, \rho^B$ respectively: $\tr_B\rho^{AB}=\rho^A$ and $\tr_A\rho^{AB}=\rho^B$.  The set of all quantum couplings of $\rho^{AB}$ is denoted by $\Gamma^{Q}(\rho^A, \rho^B)$.  Observe that $\rho^A\otimes\rho^B\in \Gamma^{Q}(\rho^A, \rho^B)$.
Let $C$ be a Hermitian matrix of order $mn$.  The {\em quantum optimal transport} problem,
 abbreviated as QOT, is defined as follows:
\begin{equation}\label{defkapCAB}
\mathrm{T}^Q_{C}(\rho^A,\rho^B)=\min_{\rho^{AB}\in \Gamma^{Q}(\rho^A,\rho^B)} \tr  C\rho^{AB}.
\end{equation}
The matrix $C$ can be viewed as a ``cost'' matrix in certain instances that will be explained later.
The quantum optimal transport has a simple operational interpretation. Suppose that Alice and Bob are two parties, who share a bipartite state $\rho^{AB}$. Their local detection statistics are fixed by the marginals $\rho^A = \tr_B\rho^{AB}$ and $\rho^B = \tr_A\rho^{AB}$. If $C$ is an effect, i.e. $0 \leq C \leq 1$, then $\mathrm{T}^Q_{C}(\rho^A,\rho^B)$ is the minimum probability of observing $C$ with fixed local states $\rho^A$, $\rho^B$. If $C$ is just positive semidefinite, then $\mathrm{T}^Q_{C}(\rho^A,\rho^B)$ is the minimum expected value of the observable $C$. For more details on the physical interpretation and applications we refer the Reader to the companion paper \cite{FECZ} and references therein.

Observe that finding the value of $\mathrm{T}^Q_{C}(\rho^A,\rho^B)$ is a semidefinite programming problem (SDP).  Using standard complexity results for SDP, as in  \cite[Theorem 5.1]{VB96},
we show that the complexity of finding
the value of $\mathrm{T}^Q_{C}(\rho^A,\rho^B)$ within a given precision $\varepsilon>0$ is polynomial in the size of the given data and $\log\frac{1}{\varepsilon}$.
There are quantum algorithms that offer a speedup for SDP \cite{brandao2017quantum}.

One of the aims of this paper is to study the properties of $\mathrm{T}_{C}^Q(\rho^A,\rho^B)$.  It is useful to compare $\mathrm{T}_{C}^Q$ with $\mathrm{T}_{C^{cl}}^{cl}$ defined as follows.  Observe that the diagonal entries of $\rho^A,\rho^B$ form probability vectors $\mathbf{p}^A,\mathbf{p}^B$.  (This corresponds to quantum decoherence, where the off-diagonal entries of $\rho^A$ and $\rho^B$ converge to zero.)  For $\mathbf{x}\in\mathbb{R}^n,X\in\mathbb{R}^{n\times n}$ denote by $\diag(\mathbf{x}),\diag(X)\in\mathbb{R}^{n\times n}$ the diagonal matrices induced by the entries  of $\mathbf{x}$ and the diagonal entries of $X$ and  respectively.
For $\mathbf{p}^A\in\Pi_m,\mathbf{p}^B\in\Pi_n$ denote by
$\Gamma_{de}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$
the convex subset of diagonal matrices in $\Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$.
We show that $\Gamma_{de}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$ is isomorphic to $\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^A)$.  Let $C^{cl}\in\mathbb{R}^{m\times n}$ the be the matrix induced by the diagonal entries of $C$ (see Section \ref{sec:diagdm}).  Then
\begin{equation}\label{qtrancheap}
\mathrm{T}_{C}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))\le \mathrm{T}_{C^{cl}}^{cl}(\mathbf{p}^A,\mathbf{p}^B) \; \textrm{ for } \;
\mathbf{p}^A\in\Pi_m, \mathbf{p}^B\in\Pi_n.
\end{equation}
We give examples where strict inequality holds.  Specific cases of this inequality were studied in \cite{CGP20}.

We now concentrate on the most important case $m=n$.  In this case we would like to find an analog of the Wasserstein-$p$ metric on $\Omega_n$.
A symmetric function sdist$:\Omega_n\times \Omega_n\to [0,\infty)$ is called a semi-metric when sdist$(\rho^A,\rho^B)=0$ if and only if $\rho^A=\rho^B$.
We show that $T^Q_C$ is a semi-distance if and only if $C$ is zero on $\mathcal{H}_S$ and $C$ is positive definite on $\mathcal{H}_A$, where $\mathcal{H}_S$ and $\mathcal{H}_A$ are the subspaces of symmetric and skew symmetric $n\times n$ matrices viewed as subspaces of $\mathbb{C}^n\otimes\mathbb{C}^n=\mathbb{C}^{n\times n}=\mathcal{H}_S\oplus\mathcal{H}_A$.  If $C$ is zero on $\mathcal{H}_S$ and positive definite on $\mathcal{H}_A$ then $\sqrt{\mathrm{T}_C^Q}$ is a weak distance: there is a metric $D'$ on $\Omega_n$ such that
$\sqrt{\mathrm{T}_C^Q(\rho^A,\rho^B)}\ge D'(\rho^A,\rho^B)$ for all $\rho^A,\rho^B\in\Omega_n$.
(One can choose $D'$ as the scaled Bures distance \cite{GLN05}.)
We show that in this case there exists a unique maximum metric $D'$ on $\Omega_n$, which can be called the quantum Wasserstein-2 metric
  and is given by the formula:
\begin{equation}\label{defQOTmet}
\mathrm{W}_C^Q(\rho^A,\rho^B)=\lim_{N\to\infty} \min_{\substack{\rho^{A_1},\ldots,\rho^{A_N}\in\Omega_n,\\ \rho^{A_0}=\rho^A, \, \rho^{A_{N+1}}=\rho^B}} \; \sum_{i=1}^{N+1}\sqrt{\mathrm{T}^Q_{C}(\rho^{A_{i-1}},\rho^{A_i})}.
\end{equation}
This metric does not seem to be easily computable for a general $C$.

The simplest example of such $C$ is $C^Q$---the orthogonal projection of $\mathbb{C}^{n\times n}$ on $\mathcal{H}_A$, as advocated in \cite{YZYY19,Duv20} and \cite{Rie18}.  It is straightforward   to show that $C^Q=\frac{1}{2}(\mathbb{I} -S)$, where $S$ is the SWAP operator $\mathbf{x}\otimes\mathbf{y}\mapsto \mathbf{y}\otimes\mathbf{x}$ and $\mathbb{I}$ is the identity operator on $\mathbb{C}^n\otimes\mathbb{C}^n$.
We show that $\big(\mathrm{T}_{C^Q}^Q\big)^{1/p}$ does not satisfy the triangle inequality for $p\in[1,2)$, and for the qubit case $n=2$, $\sqrt{\mathrm{T}_{C^Q}^Q}$ is a metric.
Hence $\mathrm{W}_{C^Q}^Q=\sqrt{\mathrm{T}_{C^Q}^Q}$ for qubits.
Furthermore  $\sqrt{\mathrm{T}_{C^Q}^Q}$ is a distance on pure states.
Numerical simulations point out that $\sqrt{\mathrm{T}_{C^Q}^Q}$ satisfies the triangle inequality for $n=3,4$ within numerical precision.  This was also noted in \cite{Rie18}.

A simple generalization of $C^Q $  is the following operator that vanishes on $\mathcal{H}_S$ and is positive definite on $\mathcal{H}_A$:
\begin{equation}\label{defCQE}
\begin{aligned}
C_E^Q=\sum_{1\le i<j\le n} e_{ij}\frac{1}{\sqrt{2}}\big(|i\rangle|j\rangle - |j\rangle|i\rangle\big)\big(\langle i|\langle j| -\langle j|\langle i|\big),\\
\text{with } e_{ij}>0 \textrm{ for } 1\le i<j \le n.
\end{aligned}
\end{equation}
Here $|1\rangle,\ldots,|n\rangle$ is any orthonormal basis in $\mathcal{H}_n$.
We show that decoherenc
of the marginal states, $\rho \to \diag (\rho)$, decreases the cost of QOT for $C^Q_E$:
\begin{equation}\label{decineq}
\mathrm{T}_{C^Q_E}^Q(\diag(\rho^A),\diag(\rho^B))\le \mathrm{T}_{C^Q_E}^Q(\rho^A,\rho^B) \textrm{ for } \rho^A, \rho^B\in\Omega_n.
\end{equation}

As in \cite{Fri20} we show that quantum transport can be defined on $d$-partite states.  In particular one can define an analog of $C^Q$ for multi-partite systems.  More precisely, $C^Q$ is the projection on the orthogonal complement of the boson subspace --- the subspace of symmetric tensors in $\otimes^d\mathbb{C}^n$.

\subsection{A brief survey of the main results}\label{subsec:survey}
Section \ref{sec:not} outlines our notation, which is a fusion of mathematical notation with Dirac's notation.  We do this to facilitate the reading of the paper by mathematicians.

In Section \ref{sec:QOT}
	we give some basic properties of the function $\mathrm{T}_C^Q$.   Proposition \ref{honconv} shows that this function is continuous and convex on $\Omega_n\times \Omega_n$.  Theorem \ref{QOTSDP}  states formally that the computation of $\mathrm{T}_C^Q$ is an SDP problem.  In particular, the computation of $\mathrm{T}_C^Q(\rho^A,\rho^B)$  within precession of $\varepsilon\in (0,1)$ is polynomial in the size of the data
	The complexity,  i.e.,  the computation time, depends on the value of $\varepsilon$:

	the smaller the $\varepsilon$ the more complex the computation, and in terms of time, the dependence is  polynomial in $\log 1/\varepsilon$.\

In Section \ref{sec:swap} we discuss QOT with respect to the SWAP operator
$S \in \mathrm{B}(\mathcal{H}_n\otimes \mathcal{H}_n)$ that swaps the two factors of
$\mathcal{H}_n\otimes \mathcal{H}_n$.  The operator $S$ has two invariant subspaces of $\mathcal{H}_n\otimes \mathcal{H}_n$, which is viewed as the set of $n\times n$ complex valued matrices $\mathbb{C}^{n\times n}$: the subspaces of symmetric and skew symmetric matrices, denoted as $\mathcal{H}_S$ and $\mathcal{H}_A$ respectively.  The subspaces $\mathcal{H}_S$ and $\mathcal{H}_A$ correspond to the eigenvalues $1$ and $-1$ of $S$ respectively.

In Section \ref{sec:metrics} we discuss metrics induced by QOT.
Theorem \ref{kapTABprop}  shows that $\mathrm{T}_C^Q$ is a semi-metric on $\Omega_n$ if and only if
$C$ is positive semidefinite and vanishes exactly on $\mathcal{H}_S$.  Furthermore,
for such $C$,
$\sqrt{\mathrm{T}_C^Q}$ is a weak metric, which induces the quantum Wasserstein-2 metric \eqref{defQOTmet}.

In Section \ref{sec:diagdm} we mainly compare the  classical and quantum optimal transports for diagonal density matrices.  For a given density matrix $\rho$ the diagonal density matrix $\diag(\rho)$ can be viewed as the decoherence of $\rho$.  Lemma \ref{diagdecr} shows that decoherence decreases the QOT for $C=C_E^Q$, cf. Formula \eqref{decineq}.  Lemma \ref{diaglemobs} gives a map of $\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$ to $\Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$.  Lemma \ref{clqotdiagdm} proves two fundamental results:  first, that the classical optimal transport is more expensive than the quantum optimal transport \eqref{qtrancheap}, and second, that $\mathrm{T}^Q_{C^Q}(\diag(\mathbf{p}^A),\diag( \mathbf{p}^B))$ can be stated as the minimum of a certain convex function on $\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$.  This shows that the computation of $\mathrm{T}^Q_{C^Q}(\diag(\mathbf{p}^A), \diag(\mathbf{p}^B))$ is simpler than the computation of $\mathrm{T}^Q_{C^Q}(\rho^A,\rho^B)$ for general $\rho^A,\rho^B$.  Theorem
\ref{QOTrhodiag} gives a closed formula for $\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)$ for two commuting qubits $\rho^A$ and $\rho^B$.

 In Section \ref{sec:decoh} we discuss the decoherence of
 the quantum cost matrix,
$C^{Q}_{\alpha}=\alpha C^Q+(1-\alpha)\diag(C^Q)$, where $\alpha\in[0,1]$.  Thus $\alpha=1$ and $\alpha=0$ correspond to QOT and OT respectively.

 Lemma \ref{Pstalphaform} gives an exact formula of the decoherence of two diagonal qubit density matrices.
 It yields that $\mathrm{T}_{C^Q_\alpha}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$ strictly decreases on the interval $[0,1]$, unless either of the states is pure or $\mathbf{p}^A = \mathbf{p}^B$.  In particular the cost of the classical optimal transport is bigger than the cost of the quantum optimal transport.

 In Section \ref{sec:dualprob} we discuss the dual problem of the SDP problem \eqref{defkapCAB}.
 Theorem~\ref{dualQOT} establishes the dual problem and shows that its resolution yields the value of $T^Q_C$.  (This was also shown in \cite{CHW19}.)
Furthermore, Theorem \ref{dualQOT} states  the complementary conditions in the case the supremum in the dual problem are achieved.  (This condition holds if $\rho^A$ and $\rho^B$ are positive definite.)  We found these complementary conditions to be very useful.  In Subsection \ref{subsec:T=T0}  we use these conditions to find a nice characterization for the cost of the quantum optimal transport for general qubit density matrices: Theorem \ref{mtqub}.  Corollary \ref{trinqub} to this theorem shows that $\sqrt{\mathrm{T}_{C^Q}^Q}$ is a metric on the qubit density matrices. Subsection~\ref{subsec:Blochb}
provides (Theorem \ref{thmC3}) a closed formula for $\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)$ in terms of solutions of the trigonometric equation \eqref{Phieq}.  Lemma \ref{6sollem} shows that this trigonometric equation is equivalent to a polynomial equation of degree
at most $6$.  Subsection \ref{subsec:isospectral} gives a nice closed formula for the value of QOT for two isospectral qubit density matrices. In Subsection \ref{subsec:nonexisF} we present a simple example where the supremum of the dual SDP problem to QOT is not achieved.  Subsection  \ref{subsec:lowbdT} gives a lower bound on $\sqrt{\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)}$ which is a metric on $\Omega_n$.  Furthermore, for $n=2$ the lower bound is equal  to $\sqrt{\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)}$.

Section \ref{subsec:diagqut} gives a closed formula for the QOT for almost all diagonal qutrits.

Section \ref{sec:QOTdpar} discusses the the quantum optimum transport for $d$-partite systems for $d\ge 3$, denoted as $\mathrm{T}_C^Q(\rho^{A_1},\ldots,\rho^{A_d})$.  The classical optimal transport of $d$-partite systems is discussed in \cite{Fri20}.  The most interesting case is where the density matrix is in $\otimes^d\mathcal{H}_n$.  Then the analog of $C^Q$ is  $C^B$---the projection on the complement of symmetric tensors.
The computation of $\mathrm{T}_{C^B}^Q(\rho^{A_1},\ldots,\rho^{A_d})$ is related to the permanent function on positive semidefinite matrices.
Assume that $d=2\ell$, where $\ell>1$. As in \cite{Fri20} one can define a Wasserstein-2 metric on the space of $\ell$-tuples of density matrices $\Omega_n^\ell$ and on the space of unordered $\ell$-tuples $\{\rho^{A_1},\ldots,\rho^{A_\ell}\}$.

We now summarize briefly the content of the Appendices.  In Appendix \ref{sec:partr} we recall briefly the basic properties of partial traces.  In Appendix \ref{sec:maxrankep}
we give an upper bound on the rank of the extreme points of the convex sets $\Gamma^Q(\rho^A,\rho^B)$, where $\rho^A\in\Omega_m,\rho^B\in\Omega_n$.  For $m=n$ our upper bound is equal to the upper bound of Parthasarathy \cite{Par05}.
Appendix \ref{sec:metrdm} discusses various metrics on density matrices.
Appendix \ref{sec:Lip} shows that $\mathrm{T}_ {C}^Q( \rho^A, \rho^B)$ is Lipschitz on the set of density matrices $\Omega_{n,a}=\{\rho\in\Omega_n, \rho\ge a\mathbb{I}_n\}$ for a fixed $a\in(0,1/n]$.   In Appendix
\ref{sec:YZYYres} we discuss the upper and lower bounds on QOT given in \cite{YZYY19}.  We reprove the lower bound for QOT since we use it in our paper.

\section{Notation}\label{sec:not}
In what follows we fuse mathematical and Dirac notations.
We view $\mathbb{C}^n$, the vector space of column vectors over the complex field $\mathbb{C}$, as a Hilbert space $\mathcal{H}_n$ with the inner product
\begin{equation*}
\langle \mathbf{y},\mathbf{x}\rangle=\mathbf{y}^\dagger\mathbf{x}=\langle \mathbf{y} |\mathbf{x}\rangle.
\end{equation*}
Then $|i\rangle\in\mathcal{H}_n$ is identified with the unit vector $\mathbf{e}_i=(\delta_{1i},\ldots,\delta_{ni})^\top$ for $i\in[n]$.
Let $\mathrm{B}(\mathcal{H}_n)\supset \mathrm{S}(\mathcal{H}_n)\supset \mathrm{S}_{+}\mathcal{H}_n\supset \Omega_n$ be  the space of linear operators, the real subspace of selfadjoint operators, the cone of positive semidefinite operators, and the convex set of density operators, respectively.   For $\rho\in \mathrm{B}(\mathcal{H}_n)$ we denote $|\rho|=\sqrt{\rho\rho^\dagger}\in \mathrm{S}_+(\mathcal{H}_n)$.  Then $\|\rho\|_1=\tr |\rho|$.
For $\rho,\sigma\in \mathrm{S}(\mathcal{H}_n)$ we write $\rho\geq \sigma$ and $\rho> \sigma$ if if the eigenvalues of $\rho-\sigma$ are all nonnegative or positive respectively.

The space of $n\times n$ complex valued matrices, denoted as  $\mathbb{C}^{n\times n}$, is a representation of $\mathrm{B}(\mathcal{H}_n)$, where the matrix $\rho=[\rho_{ij}]\in\mathbb{C}^{n\times n}$ represents the operator $\rho\in\mathrm{B}(\mathcal{H}_n)$.  The set  of density operators in $\mathrm{B}(\mathcal{H}_n)$ are viewed as  $\Omega_n$: the convex set of $n\times n$ Hermitian positive semidefinite trace-one matrices.
The tensor product $\mathcal{H}_m\otimes\mathcal{H}_n$ is represented by $\mathbb{C}^{m\times n}$.
An element is denoted by a matrix $X=[x_{ip}]=\sum_{i=p=1}^{m,n} x_{ip} |i\rangle |p\rangle$, which correspond to a bipartite state.  Observe that $\mathbf{x}\otimes\mathbf{y}=|\mathbf{x}\rangle|\mathbf{y}\rangle$ is represented by the rank-one matrix $\mathbf{x}\mathbf{y}^\top$.
We denote by $X^\dagger=\langle X|$ the complex conjugate of the transpose of  $X\in\mathbb{C}^{m\times n}$.
The inner product of bipartite states $X,Y\in\mathbb{C}^{m\times n}$ is $\langle X,Y\rangle=\langle X|Y\rangle=\tr X^\dagger Y$.
 We identify  $\mathrm{B}(\mathcal{H}_m\otimes \mathcal{H}_n)$ with $\mathbb{C}^{(mn)\times (mn)}$ as follows. An operator $\rho^{AB}\in \mathrm{B}(\mathcal{H}_m\otimes\mathcal{H}_n)$ is represented by a matrix $R\in\mathbb{C}^{(mn)\times (mn)}$, whose entries are indexed with two pairs of indices $r_{(i,p)(j,q)}$ where $i,j\in [m], p,q\in[n]$.  Then the partial traces of $R$ are defined as follows:

\begin{equation}\label{ptfrom}
\tr_A R=[\sum_{i=1}^m r_{(i,p)(i,q)}]=\rho^B\in\mathbb{C}^{n\times n}, \quad \tr_B R=[\sum_{p=1}^n r_{(i,p)(j,p)}]=\rho^A\in\mathbb{C}^{m\times m}.
\end{equation}
Recall that $\tr R=\tr (\tr_A R)=\tr (\tr_B R)$.
Some more known facts about partial traces that we use in this paper are discussed in the Appendix \ref{sec:partr}.

Let $\mathrm{M}:\mathrm{B}(\mathcal{H}_m\otimes\mathcal{H}_n)\to \mathrm{B}(\mathcal{H}_m)\oplus\mathrm{B}(\mathcal{H}_n)$ be the partial trace map: $\rho^{AB}\mapsto (\rho^A, \rho^B)$.  We identify $\mathrm{M}$ with the map $\mathrm{M}: \mathbb{C}^{(mn)\times (mn)}\to\mathbb{C}^{m\times m}\oplus\mathbb{C}^{n\times n}$.
For $ \rho^A\in\Omega_m, \rho^B\in\Omega_n$ we denote by $\Gamma^Q(\rho^A,\rho^B)$ the set of all bipartite density matrices whose partial traces are $\rho^A$ and $\rho^B$ respectively:
$$\Gamma^Q(\rho^A,\rho^B)=\{ \rho^{AB}\in\Omega_{mn}, \tr_B \rho^{AB}= \rho^A, \tr_A \rho^{AB}= \rho^B\}.$$
Then $\Omega_{mn}$ fibers over $\Omega_m\times \Omega_n$, that is, $\Omega_{mn}=\bigcup_{( \rho^A, \rho^B)\in \Omega_m\times\Omega_n} \Gamma^Q( \rho^A, \rho^B)$.  The Hausdorff distance between $\Gamma^Q( \rho^A, \rho^B)$ and $\Gamma^Q(\rho^C,\rho^D)$ is a complete metric on the fibers \cite{FGZ19}.

On the side, we note that bipartite density operators $\rho^{AB}$ play an important role in uniform continuity bounds for quantum entropies \cite{Win16}.

\section{Quantum Optimal Transport is a  Semidefinite programming  problem}\label{sec:QOT}
\begin{proposition}\label{honconv} For $C\in \mathrm{S}(\mathcal{H}_m\otimes \mathcal{H}_n)$  the function
$\mathrm{T}_{C}^Q(\cdot,\cdot)$ is a continuous convex function on $\Omega_m\times\Omega_n$: for any $0<a<1$,
\begin{align*}
\mathrm{T}_{C}^Q(a\rho^A+(1-a)\sigma^A,a\rho^B+(1-a)\sigma^B))\le
a\mathrm{T}_{C}^Q(\rho^A,\rho^B)+(1-a)\mathrm{T}_{C}^Q(\sigma^A,\sigma^B).
\end{align*}
Furthermore, if $C\ge  0$ then  $\mathrm{T}_{C}^Q(\cdot,\cdot)$ is nonnegative.
\end{proposition}
\begin{proof}
Assume that
\begin{align*}
&\mathrm{T}_{C}^Q(\rho^A,\rho^B)=\tr C \rho^{AB}, &&\hspace*{-2cm} \rho^{AB}\in \Gamma^Q(\rho^A,\rho^B),\\
&\mathrm{T}_{C}^Q(\sigma^A,\sigma^B)=\tr C \sigma^{AB}, &&\hspace*{-2cm} \sigma^{AB}\in \Gamma^Q(\sigma^A,\sigma^B).
\end{align*}
Let $\tau^{AB}=a\rho^{AB}+(1-a)\sigma^{AB}$. Then $\tau^{AB}\in\Gamma^Q(a\rho^A+(1-a)\sigma^A,a\rho^B+(1-a)\sigma^B)$.
Clearly $\tr C \tau^{AB}=a\mathrm{T}_{C}^Q(\rho^A,\rho^B)+(1-a)\mathrm{T}_{C}^Q(\sigma^A,\sigma^B)$.  The minimal characterization~\eqref{defkapCAB} of $\mathrm{T}$ yields the first inequality of the lemma.  Clearly if $C\ge 0$ then $\mathrm{T}_{C}^Q(\cdot,\cdot)$ is nonnegative.
This yields the second inequality of the lemma.

The continuity of $\mathrm{T}_{C}^Q(\cdot,\cdot)$ follows from the following argument.  First observe that for each $\rho^A\in\Omega_m,\rho^B\in\Omega_n$, the set $\Gamma^Q(\rho^A,\rho^B)$, viewed as a fiber over $(\rho_A,\rho_B)$, is a compact convex set.  Hence one can define the Hausdorff metric (distance) on the fibers.
It is shown in \cite[Theorem 5.2]{FGZ19} that the Hausdorff metric is a complete metric.  Furthermore the sequence $\Gamma^Q(\rho^{A,k},\rho^{B,k}),k\in\mathbb{N}$ converges to
$\Gamma^Q(\rho^A,\rho^B)$ in the Hausdorff distance if and only if $\lim_{k\to\infty}(\rho^{A,k},\rho^{A,k})=( \rho^A, \rho^B)$.  This proves the continuity of $\mathrm{T}_{C}^Q(\cdot,\cdot)$.
\end{proof}

For a selfadjoint operator $\rho\in\mathrm{S}(\mathcal{H}_n)$ we denote by $\lambda_{\max}(\rho)=\lambda_1(\rho)\ge \cdots\ge\lambda_n(\rho)=\lambda_{\min}(\rho)$ the $n$ eigenvalues of $\rho$.  For $a\in[0,1/n]$ we denote
by $\Omega_{n,a}$ all density matrices that satisfy the inequality $\lambda_{\min}\ge a$.
Note that $\Omega_{n,0}=\Omega_n$.  In Appendix \ref{sec:Lip} we show that $\mathrm{T}_{C}^Q(\cdot,\cdot)$ is Lipschitz on $\Omega_{n,a}\times\Omega_{n,a}$ for $a\in(0,1/n)$.

The following Proposition shows that to compute $\mathrm{T}_{C}^Q( \rho^A, \rho^B)$ one can assume that the eigenvalues of $C$ are in the interval $[0,1]$:
 %
\begin{proposition}\label{propkapCC'}  Assume that $C\in \mathrm{S}(\mathcal{H}_m\otimes \mathcal{H}_n)$ is not a scalar operator ($C\ne c \mathbb{I}$).  Let
\begin{equation*}
\tilde C=\frac{1}{\lambda_{\max}(C)-\lambda_{\min}(C)} \big(C-\lambda_{\min}(C) \mathbb{I} \big).
\end{equation*}
Then $0\le \tilde C\le \mathbb{I}$.  Furthermore for $\rho^A\in\Omega_m,\rho^B\in\Omega_n$  the following equality holds:
\begin{equation}\label{kapCC'rel}
\mathrm{T}_{C}^Q(\rho^A,\rho^B)=(\lambda_{\max}(C)-\lambda_{\min}(C))\mathrm{T}_{\tilde C}^Q(\rho^A,\rho^B)+\lambda_{\min}(C).
\end{equation}
\end{proposition}
\begin{proof} Clearly $C=(\lambda_{\max}(C)-\lambda_{\min}(C))\tilde C+\lambda_{\min}(C)\mathbb{I}$.  Furthermore
\begin{align*}
\tr C \rho^{AB}=(\lambda_{\max}(C)-\lambda_{\min}(C))\tr \tilde C\rho^{AB}+\lambda_{\min}(\tilde C), \quad \rho^{AB}\in\Gamma^Q(\rho^A,\rho^B).
\end{align*}
As $\lambda_{\max}(C)-\lambda_{\min}(C)> 0$ we deduce \eqref{kapCC'rel}.
\end{proof}

We next observe that one can reduce the computation of $\mathrm{T}_{C}^Q(\rho^A,\rho^B)$ to a smaller
dimension problem if either $\rho^A$ or $\rho^B$ are not positive definite:
\begin{proposition}\label{redprop}    Assume that  $\rho^A\in\Omega_m, \rho^B\in\Omega_n$. Let $m'$ and $n'$ be the dimensions of\, $\range\rho^A=\mathcal{H}_{m'}$ and\, $\range\rho^B=\mathcal{H}_{n'}$ respectively.  Denote by $\,\rho^{A'}\in\Omega_{m'}$, and $\rho^{B'}\in\Omega_{n'}$ the restrictions of $\rho^A$ and $\rho^B$ to $\mathcal{H}_{m'}$ and $\mathcal{H}_{n'}$ respectively.  Assume that  $C\in\mathrm{S}(\mathcal{H}_m\otimes \mathcal{H}_n)$, and denote by $C'\in\mathrm{S}(\mathcal{H}_{m'}\otimes \mathcal{H}_{n'})$ the restriction of $C$ to $\mathcal{H}_{m'}\otimes \mathcal{H}_{n'}$.  Then
\begin{equation*}
\mathrm{T}^Q_{C}(\rho^A,\rho^B)=\mathrm{T}^Q_{C'}(\rho^{A'},\rho^{B'}).
\end{equation*}
\end{proposition}
\begin{proof}
Without loss of generality we can assume that we chose orthonormal bases in $\mathcal{H}_m$ and $\mathcal{H}_n$ to be the eigenvectors of $\rho^A$ and $\rho^B$ respectively.
Thus to prove the lemma it is enough to consider the following case: $\rho^A=\rho^C\oplus 0_{m-l}$ where $\rho^C\in \Omega_l, l<m$ and $0_l$ is an $l\times l$ zero matrix.  Let $\tilde C\in \mathrm{S}(\mathcal{H}_l\otimes\mathcal{H}_n)$ be the restriction of $C$ to $\mathcal{H}_l\otimes\mathcal{H}_n$.  We claim that
\begin{equation}\label{redprop1}
\mathrm{T}^Q_{C}(\rho^A,\rho^B)=\mathrm{T}^Q_{\tilde C}(\rho^{C},\rho^{B}).
\end{equation}
Let $R=[R_{(i,p)(j,q)}]\in\Gamma^Q(\rho^A,\rho^B)$.  As $R\ge 0$ it follows that the submatrix $R_{ii}=[R_{(i,p)(i,q)}], p,q\in[n]$ is positive semidefinite for each $i\in[m]$.
Since $\tr _B R=\rho^A$ we deduce that $\rho^A_{ii}=\sum_{p\in[n]} R_{(i,p)(i,p)}=\tr R_{ii}=0$ for $i>l$.   Therefore $R_{ii}=0$, that is, $R_{(i,p)(i,q)}=0$ for $p,q\in[n]$ and $i>l$.  Let $R'$ be the following submatrix of $R$: $[R_{(i,p)(j,q)}], i,j\in [l], p,q\in[n]$.  Then $R'\in\Gamma^Q(\rho^C,\rho^B)$.  Vice versa, given $R'\in\Gamma^Q(\rho^C,\rho^B)$, one can enlarge trivially $R'$ to $R$ in $\Gamma^Q(\rho^C, \rho^B) $.  Clearly $\tr CR=\tr \tilde C R'$.  Repeating the same process with $\rho^B$ establishes \eqref{redprop1}.
\end{proof}

As we point out in the next section it is natural to consider the case $m=n$.  However, if either $\rho^A$ or $\rho^B$ are singular density matrices then we can reduce the computation of $T_{C}^Q(\rho^A,\rho^B)$ to a lower-dimensional problem, and after this reduction it may happen that the dimensions are no longer equal.

One of the main results of this paper is the observation that the computation of the quantum transport is carried out efficiently using semidefinite programming \cite{VB96}.  We will sometimes use the abbreviation SDP for semidefinite programming.
\begin{theorem}\label{QOTSDP} Assume that $C\in \mathrm{S}(\mathcal{H}_m\otimes \mathcal{H}_n)$, $\rho^A\in\Omega_m, \rho^B\in\Omega_n$.
Then the computation of $\mathrm{T}_{C}^Q(\rho^A,\rho^B)$ is a semidefinite programming problem.  The value  of $\mathrm{T}_{C}^Q(\rho^A,\rho^B)$ can be approximated within precision  $\varepsilon>0$ in polynomial time in the size of the data and $\log1/\varepsilon$.
\end{theorem}
\begin{proof}
Assume that $\rho^A=[a_{ij}]\in\Omega_m, \rho^B=[b_{pq}]\in\Omega_n$.
Denote the entries of the Hermitian matrix $C$ by $c_{(i,p)(j,q)}$, i.e.,  $c_{(i,p)(j,q)}=\overline{c_{(j,q)(i,p)}}$. Let $\mathbf{i}=\sqrt{-1}$, and
\begin{align*}
& E_{ij}^A=|i\rangle\langle j|, && G_{ij}^A=\frac{1}{2}(E_{ij}^A+E_{ji}^A), && H_{ij}^A=\frac{1}{2}\mathbf{i} (E_{ij}^A-E_{ji}^A),&& i,j\in[m],\\
& E_{pq}^B=|p\rangle\langle q|, && G_{pq}^B=\frac{1}{2}(E_{pq}^B+E_{qp}^B), && H_{pq}^B=\frac{1}{2}\mathbf{i} (E_{pq}^B-E_{qp}^A), && p,q\in[n].
\end{align*}
Thus $|i\rangle,i\in[m]$, $E_{ij}^A,i,j\in[m]$,  $G_{ij}^A, 1\le i\le j\le m, H_{ij}^A , 1\le i< j\le m$ are the standard bases in $\mathbb{C}^m$,  $\mathbb{C}^{m\times m}$, and in the subspace of $m\times m$ Hermitian matrices respectively.  A similar observation applies when we replace $A$ and $m$ by $B$ and $n$.
The conditions $\tr_B \rho^{AB}=\rho^A, \tr_A \rho_{AB}=\rho^B$ are stated as the following linear conditions:
\begin{equation}\label{margcond}
\begin{aligned}
&\tr \rho^{AB}(G_{ij}\otimes \mathbb{I}_n)=\Re a_{ij}, && i\le j, && \tr \rho^{AB} (H_{ij}\otimes \mathbb{I}_n)=\Im a_{ij}, && i<j,\\
&\tr \rho^{AB}(\mathbb{I}_m\otimes G_{pq})=\Re b_{pq}, && p\le q, && \tr \rho^{AB} (\mathbb{I}_m\otimes H_{pq})=\Im b_{pq}, && p<q.
\end{aligned}
\end{equation}
Here $\Re z,\Im z$ are the real and the imaginary part of the complex number $z\in\mathbb{C}$.
We assume that $ \rho^{AB}\geq 0$.
Hence $\mathrm{T}_{C}^Q(\rho^A,\rho^B)$ is a semidefinite problem for $\rho^{AB}$.

Assume first that $\rho^A, \rho^B$ are positive definite.  Then $ \rho^A\otimes  \rho^B$, viewed as a Kronecker tensor product, is positive definite. Thus $\Gamma^Q(\rho^A,\rho^B)$ contains a positive definite operator $ \rho^A\otimes  \rho^B$.  The standard SDP theory \cite[Theorem 5.1]{VB96} yields that $\mathrm{T}_{C}^Q(\rho^A,\rho^B)$ can be computed in polynomial time with precision $\varepsilon>0$.

(Note that the standard SDP is stated for real symmetric positive semidefinite matrices.  It is well known that Hermitian positive semidefinite matrices can be encoded as special real symmetric matrices of double dimension.  See the proof of Theorem \ref{dualQOT} for details.)

Assume that $\rho^A,\rho^B\ge 0$.  Then the restrictions
$\rho^{A'}=\rho^A|_{\range\rho^A}$ and $\rho^{B'}=\rho^B|_{\range\rho^B}$ are positive definite.  Use
Proposition \ref{redprop}
to deduce that $\mathrm{T}_{C}^Q(\rho^A,\rho^B)$ can be computed in polynomial time in precision $\varepsilon>0$.
\end{proof}

We remark that one can try to generalize $\mathrm{T}_{C}^Q(\rho^A,\rho^B)$ to non-Hermitian matrices $C\in\mathrm{B}(\mathcal{H}_m\otimes\mathcal{H}_n)$ by defining $\mathrm{T}_C^Q(\rho^A, \rho^B)$ as the minimum of the real functional $\Re \tr C\rho^{AB}$ over all $\rho^{AB} \in \Gamma^Q(\rho^A, \rho^B)$.  Clearly
\begin{equation*}
\Re \tr C\rho^{AB}=\tr \hat C \rho^{AB}, \quad \hat C=\frac{1}{2}(C+C^\dagger), \quad \rho^{AB}\in\mathrm{S}(\mathcal{H}_m\otimes\mathcal{H}_n).
\end{equation*}
Hence
$\mathrm{T}_{C}^Q(\rho^A,\rho^B)= \mathrm{T}_{\hat C}^Q(\rho^A,\rho^B)$.

\section{Quantum transport problem induced by SWAP}\label{sec:swap}
When describing any two distinguishable physical objects one can
introduce an operation which exchanges them.
On the composite space $\mathcal{H}_n\otimes \mathcal{H}_n$
it corresponds to
 a natural isometry induced by swapping the two factors $\mathbf{x}\otimes \mathbf{y}\mapsto \mathbf{y}\otimes \mathbf{x}$. On the space of square matrices  the SWAP operator is the map $X\mapsto X^\top$.  This map is of fundamental importance in quantum information theory.  It allows to observe some interesting properties of bipartite system
 and is useful
 in the criterion for separability by Peres and Horodecki \cite{Per96, Hor96}.

 We will see below that if we let $S$ denote the SWAP operator,
 then it induces a cost matrix \[C^Q = \frac12(\mathbb I - S)\]
  for the quantum transport problem
   which enjoys several nice properties.

We identify $\mathcal{H}_n\otimes \mathcal{H}_n$ as the space of $n\times n$ complex valued matrices
$\mathbb{C}^{n\times n}$
as follows:   Let $\mathbf{e}_i=(\delta_{i1},\ldots,\delta_{in})^\top\equiv|i\rangle, i\in[n]$ be the standard basis in $\mathbb{C}^n\equiv \mathcal{H}_n$.  Then a state $|\psi\rangle\in \mathcal{H}_n\otimes \mathcal{H}_n$ is given by $|\psi\rangle=\sum_{i=j=1}^n x_{ij}|i\rangle|j\rangle$.  Thus we associate with $|\psi\rangle$ the matrix $X=[x_{ij}]\in \mathbb{C}^{n\times n}$.  Then $|\psi\rangle$ is a normalized state if and only if $\|X\|^2=\tr X  X^\dagger=1$.
Suppose we change the orthonormal basis $\mathbf{e}_1,\ldots,\mathbf{e}_n$ to an orthonormal basis $\mathbf{f}_1,\ldots,\mathbf{f}_n$, where $\mathbf{e}_i=\sum_{p=1}^n u_{pi}\mathbf{f}_p$.  Here $U=[u_{ip}]\in\mathbb{C}^{n\times n}$ is a unitary matrix.  Then $|\psi\rangle=\sum_{p=q=1}^n y_{pq}|\mathbf{f}_p\rangle|\mathbf{f}_q\rangle$, where $Y=UXU^\top$.

We now consider a pure state density operator
\begin{equation*}
|\psi\rangle \langle \psi|=\Big(\sum_{i=j=1}^n x_{ij}|i\rangle|j\rangle\Big) \Big(\sum_{p=q=1}^n \bar x_{pq}\langle p|\langle q|\Big)=\sum_{i=j=p=q=1}^n x_{ij}\bar x_{pq} |i\rangle|j\rangle \langle p|\langle q|.
\end{equation*}
We identify the coefficient matrix with the Kronecker product $X\otimes \bar X$.
Then
\begin{align*}
\rho^A & =\tr_B |\psi\rangle \langle \psi|=\sum_{i=p=1}^n (X X^\dagger)_{ip} |i\rangle \langle p|,\\
\rho^B & =\tr_A |\psi\rangle \langle \psi|=\sum_{j=q=1}^n (X^\top  \bar X)_{jq} |j\rangle \langle q|.
\end{align*}
Thus in the standard basis of $\mathcal{H}_n$ we can identify $\rho^A$ and $\rho^B$ with the density matrices
\begin{equation}\label{Xrhosigrel}
\rho^A = X X^\dagger, \quad \rho^B=X^\top  \bar X.
\end{equation}

Suppose we change from the standard basis $\mathbf{e}_1,\ldots,\mathbf{e}_n$ to the basis $\mathbf{f}_1,\ldots,\mathbf{f}_n$ using the unitary matrix $U$.  Then $\rho^A$ and $\rho^B$ are represented as the following density matrices
\begin{equation}\label{tilderhosigfor}
\begin{aligned}
\tilde \rho^A &= \tilde X \tilde X^\dagger \, = \, U(X X^\dagger) U^\dagger \, = \, U\rho^A U^\dagger, \\
\tilde \rho^B &= \tilde X^\top  \overline{\tilde X}=U(X^\top  \bar X) U^\dagger=U\rho^B U^\dagger.
\end{aligned}
\end{equation}

Note that if $\nu_1\ge \cdots\ge\nu_n\ge 0$ are the singular values of the matrix $X$ then $\lambda_1=\nu_1^2\ge \cdots \ge \lambda_n=\nu_n^2\ge 0$ are the eigenvalues of $\rho^A$ and $\rho^B$.  That is $\rho^A$ and $\rho^B$ are isospectral.  Vice versa:
\begin{proposition}\label{exrank1}  Let $\rho^A,\rho^B\in\Omega_n$.  Then $\Gamma^Q(\rho^A,\rho^B)$ contains a matrix $R$ of rank one if and only if $\rho^A$ and $\rho^B$ are isospectral.
\end{proposition}
\begin{proof}  Suppose first that $\rho^A$ and $\rho^B$ are isospectral, i.e., have the same eigenvalues $\lambda_1\ge\cdots\ge\lambda_n\ge 0$.  Assume that $\rho^A$ and $\rho^B$ have the following spectral decompositions:
\begin{align}\label{specdecrho}
\begin{split}
\rho^A & =\sum_{i=1}^n \lambda_i |\mathbf{x}_i\rangle\langle \mathbf{x}_i|, \quad \langle \mathbf{x}_i,\mathbf{x}_j\rangle =\delta_{ij},\\
\rho^B & =\sum_{j=1}^n \lambda_i |\mathbf{y}_j\rangle\langle \mathbf{y}_j|, \quad\! \langle \mathbf{y}_i,\mathbf{y}_j\rangle =\delta_{ij}.
\end{split}
\end{align}
Then $\Gamma^Q(\rho^A, \rho^B)$ contains the rank-one matrix

\begin{equation}\label{simpur}
R= \Big(\sum_{i=1}^n \sqrt{\lambda_i} |\mathbf{x}_i\rangle |\mathbf{y}_i\rangle \Big) \Big(\sum_{j=1}^n \sqrt{\lambda_j}\langle \mathbf{x}_j|\langle \mathbf{y}_j| \Big).
\end{equation}
Vice versa, if $R$ is a pure bipartite state in $\mathrm{S}_+(\mathcal{H}_n\otimes \mathcal{H}_n)$ then it has the above decomposition, when using the Schmidt, also known as Singular Value Decomposition (SVD) \cite{Frb16}.  Hence $\tr_A R$ and $\tr_B R$ are isospectral density matrices.
\end{proof}

For $\mathcal{H}_n\otimes\mathcal{H}_n$ the SWAP operation $S\in\mathrm{B}(\mathcal{H}_n\otimes\mathcal{H}_n)$ acts on the product states as follows: $S(|\mathbf{x}\rangle|\mathbf{u}\rangle)=|\mathbf{u}\rangle| \mathbf{x}\rangle$.  So $S$ is both unitary and an involution operator:
$S^\dagger S=I,S^2=I$. Hence the eigenvalues of $S$ are $\pm 1$ and $S$ is selfadjoint, $S^\dagger=S$.  The invariant subspaces of $S$ corresponding to the eigenvalues $1$ and $-1$ are the symmetric and skew-symmetric tensors respectively,
which can be identified with the symmetric  $\mathcal{H}_A=\mathrm{S}^2\mathbb{C}^n$ and skew symmetric $\mathcal{H}_A=\mathrm{A}^2\mathbb{C}^n$  matrices in $\mathbb{C}^{n\times n}$ respectively.
Note that the decomposition of a  matrix $X$ into a sum of symmetric and skew symmetric matrices $X=(1/2)(X+X^\top)+(1/2)(X-X^\top)$ is an orthogonal decomposition.  That is
\begin{equation*}
\mathcal{H}_n\otimes\mathcal{H}_n=\mathcal{H}_S\oplus \mathcal{H}_A=\mathbb{C}^{n\times n}=\mathrm{S}^2\mathbb{C}^n\oplus \mathrm{A}^2\mathbb{C}^n
\end{equation*}
is an orthogonal decomposition.  Observe that $S(X)=X^\top$.   Hence the action of $S$ on a rank-one operator $|X\rangle\langle Y|$ in $\mathrm{B}(\mathcal{H}_n\otimes \mathcal{H}_n)$ is $S(|X\rangle \langle Y|)=|X^\top\rangle \langle Y|$.  Therefore the action of $S$ on rank one product operator in $\mathrm{B}(\mathcal{H}_n\otimes\mathcal{H}_n)$ is given by
\begin{eqnarray*}
S(|\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{y}|\langle\mathbf{v}|)=S(|\mathbf{x}\rangle|\mathbf{u}\rangle)\langle \mathbf{y}|\langle \mathbf{v}|=|\mathbf{u}\rangle |\mathbf{x}\rangle \langle \mathbf{y}|\langle\mathbf{v}|.
\end{eqnarray*}
Hence
\begin{eqnarray*}
\tr S(|\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{y}|\langle\mathbf{v}|)=
(\langle \mathbf{y}|\langle\mathbf{v}|)(|\mathbf{u}\rangle|\mathbf{x}\rangle)=\langle \mathbf{y}|\mathbf{u}\rangle\langle \mathbf{v}|\mathbf{x}\rangle.
\end{eqnarray*}
Similarly
\begin{eqnarray*}
S(|\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{y}|\langle\mathbf{v}|)S^\dagger=
|\mathbf{u}\rangle |\mathbf{x}\rangle\langle\mathbf{v}|\langle\mathbf{y}|.
\end{eqnarray*}
Use the identity \eqref{tenprodiden} and the above results to deduce that
\begin{eqnarray*}
&&\tr S(|\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{y}|\langle\mathbf{v}|)=
\langle \mathbf{y}|\mathbf{u}\rangle \langle \mathbf{v}|\mathbf{x}\rangle=\tr ((|\mathbf{x}\rangle \langle\mathbf{y}|)\otimes(|\mathbf{u}\rangle \langle\mathbf{v}|)),\\
&&\tr_A S(|\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{y}|\langle\mathbf{v}|)S^\dagger=
\langle \mathbf{v}|\mathbf{u}\rangle|\mathbf{x}\rangle\langle \mathbf{y}|=\tr_B |\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{y}|\langle\mathbf{v}|,\\
&&\tr_B S(|\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{y}|\langle\mathbf{v}|)S^\dagger=
\langle \mathbf{y}|\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{v}|=\tr_A |\mathbf{x}\rangle|\mathbf{u}\rangle\langle \mathbf{y}|\langle\mathbf{v}|.
\end{eqnarray*}
Use \eqref{tenprodiden} to deduce
\begin{equation*}
S((|\mathbf{x}\rangle \langle \mathbf{y}|)\otimes (|\mathbf{u}\rangle\langle \mathbf{v}|))=|\mathbf{u}\rangle |\mathbf{x}\rangle \langle\mathbf{y}|\langle\mathbf{v}|=(|\mathbf{u}\rangle\langle\mathbf{y}|)\otimes (|\mathbf{x}\rangle\langle \mathbf{v}|).
\end{equation*}
Combine the above equalities to obtain the following identities:
\begin{align}\label{Werid}
\begin{split}
& \tr S(\rho^A\otimes\rho^B)=\tr \rho^A\rho^B, \quad
\rho^A,\rho^B\in\mathrm{B}(\mathcal{H}_n),\\
& \tr_A S\rho^{AB}S^\dagger=\tr_A \rho^{AB}, \;\; \tr_B S\rho^{AB}S^\dagger=\tr_B \rho^{AB}, \quad \rho^{AB}\in\mathrm{B}(\mathcal{H}_n\otimes\mathcal{H}_n).
\end{split}
\end{align}
The first identity is due to Werner \cite{Wer89}, see also \cite{MPHUZ}.

Denote by $\ker C$ the kernel of a linear operator $C:\mathcal{H}_n\otimes\mathcal{H}_n\to :\mathcal{H}_n\otimes\mathcal{H}_n$.  An operator $C$  is said to vanish exactly on symmetric matrices if $\ker C=\mathcal{H}_S$. Thus a positive semidefinite $C$  vanishes exactly on $\mathcal{H}_S$ if and only if it has $n(n-1)/2$ positive eigenvalues (counting with multiplicities) with the corresponding skew symmetric eigenvectors.

Let $|1\rangle,\ldots,|n\rangle$ be an orthonormal basis in $\mathcal{H}_n$.  Define (as in \cite{FECZ}) the maximally entangled singlet states spanned on two dimensional subspaces:
\begin{equation}\label{SOBHA}
|\psi_{ij}^-\rangle=\frac{1}{\sqrt{2}}\big(|i\rangle|j\rangle -|j\rangle|i\rangle\big) \textrm{ for } 1\le i<j\le n.
\end{equation}
Given a matrix $E = [e_{ij}]_{i, j = 1}^n$ with $e_{ij} > 0$ for all $1 \le i < j \le n$, the following operator is positive semidefinite and vanishes exactly on
the symmetric subspace, $\mathrm{S}^2\mathbb{C}^n$
\cite[(11)]{FECZ}:
\begin{equation}\label{defCEop}
C^Q_E=\sum_{1\le i<j\le n} e_{ij} |\psi_{ij}^-\rangle \langle \psi_{ij}^-|,
\end{equation}
Consider the operator
\begin{equation}\label{defT}
C^Q=\frac{1}{2}(\mathbb{I} -S).
\end{equation}
Then $C^Q$ is an orthogonal projection of $\mathbb{C}^{n\times n}$ onto
antisymmetric subspace, $\mathrm{A}^2\mathbb{C}^n$.
 Hence $C^Q$ is of the form \eqref{defCEop} , where $e_{ij}=1$ for all $i < j$.
conjunit
Denote by $\mathrm{U}(n)\subset \mathbb{C}^{n\times n}$ the group of unitary matrices. The following lemma shows that $\mathrm{T}_{C^Q}^Q$ is invariant under conjugation by a unitary matrix:
\begin{proposition}\label{conjunit}  Assume that $\rho^A,\rho^B\in\Omega_n$ and $\rho^{AB}\in\Gamma^Q(\rho^A,\rho^B)$.  Then for $U\in\mathrm{U}(n)$ the following equalities hold:
\begin{equation}\label{conjunit1}
\begin{aligned}
& \tr_B((U\otimes U) \rho^{AB} (U^\dagger\otimes U^\dagger)) = U\rho^A U^\dagger, \\
& \tr_A ((U\otimes U) \rho^{AB} (U^\dagger\otimes U^\dagger)) = U\rho^B U^\dagger,\\
& (U\otimes U) \Gamma^Q(\rho^A,\rho^B) (U^\dagger\otimes U^\dagger) =\Gamma^Q(U\rho^A U^\dagger,U \rho^B U^\dagger),\\
& \mathrm{T}_{C}^Q(\rho^A,\rho^B) =\mathrm{T}_{(U\otimes U)C(U^\dagger\otimes U^\dagger)}^Q(U\rho^A U^\dagger,  U\rho^B U^\dagger).
\end{aligned}
\end{equation}
In particular
\begin{equation}\label{CQinvQOT}
\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)=\mathrm{T}_{C^Q}^Q(U\rho^A U^\dagger,  U\rho^B U^\dagger).
\end{equation}
\end{proposition}
\begin{proof}
Assume that $R$ is a pure state $R=|\psi\rangle \langle \psi|$.  The state $|\psi\rangle$ corresponds to a matrix $X\in\mathbb{C}^{n\times n}$ with $\tr X X^\dagger=1$.  Then $\tr_B R=X X^\dagger$ and $\tr_A R=X^\top \bar X$.
 Recall that $(U\otimes U)|\psi\rangle $ is represented by $\tilde X=UX U^\top$

   Now use \eqref{tilderhosigfor} to deduce the first two equalities in \eqref{conjunit1} if $R\in\Gamma^Q(\rho^A,\rho^B)$.
Recall that any $\rho^{AB}\in \Gamma^Q(\rho^A,\rho^B)$ is a convex combination of pure states $R_i=|\psi_i\rangle \langle \psi_i|, i \in[k]$.  That is $R=\sum_{i=1}^k a_i R_i$, where $a_i>0$ and $\sum_{i=1}^k a_i=1$.  Then $\tr_B R_i=\rho_i^A, \tr_A R_i=\rho^B_i$.
Now use the above results for $R_i$ to deduce the first two equalities in \eqref{conjunit1}.
The other equalities of \eqref{conjunit1} are deduced easily from the first two equalities in \eqref{conjunit1}.  Equality \eqref{CQinvQOT} is deduced from the equality
\begin{equation}\label{invCqt}
(U\otimes U) C^{Q}(U^\dagger\otimes U^\dagger)=C^{Q}.\qedhere
\end{equation}
\end{proof}
\section{Metrics induced by the quantum optimal transport}\label{sec:metrics}
Let $X$ be a set of points.
Assume that $D:X\times X\to\mathbb{R}_+(=[0,\infty))$.  Then $D(\cdot,\cdot)$ is called a metric on $X$ if it satisfies the following three properties:
 \begin{enumerate}[(a)]
 \item Symmetry:  $D(x,y)=D(y,x)$;
 \item Positivity: $D(x,y)\ge 0$, and equality holds if and only if $x=y$.
 \item Triangle inequality: $D(x,y)+D(y,z)\ge D(x,z)$.
 \end{enumerate}

 We call $D(\cdot,\cdot)$ a semi-metric if it satisfies the above first two conditions.  A semi-metric is called a weak metric if there exists a metric $D'(\cdot,\cdot)$ such that
\begin{equation}\label{DmajD'}
D'(x,y)\le D(x,y) \textrm{ for all }x,y\in X.
\end{equation}
\begin{proposition}\label{inddist}
Assume that $D$ is a weak metric on the space $X$ satisfying \eqref{DmajD'}, where $D'$ is a metric on $X$.  For each positive integer $N$ define the following function:
\begin{equation*}
D_N(x,y)=\inf_{\substack{z_1,\ldots,z_N\in X,\\mathbf{z}_0=x, \; z_{N+1}=y}} \; \sum_{i=0}^N D(z_i,z_{i+1}) \textrm{ for } x,y\in X.
\end{equation*}
Then
\begin{enumerate}[(a)]
\item For each $N$ the function $D_N(\cdot,\cdot)$ is a weak metric
that satisfies the inequality \eqref{DmajD'}.
\item For each $x,y\in X$ and $N$ we have the inequalities $0\le D_{N+1}(x,y)\le D_N(x,y)\le D(x,y)$.
\item For each $M,N\ge 1$ we have the inequality
\begin{equation*}
D_M(x,u)+D_N(u,y)\ge D_{M+N+1}(x,y) \textrm{ for } x,y,u\in X.
\end{equation*}
\item Denote by $D_\infty(x,y)=\lim_{N\to\infty}D_N(x,y)$.  Then $D_\infty(x,y)$ is a metric, called the induced metric of $D$.
Furthermore, $D_\infty$ is the maximum metric $D'$ that satisfies \eqref{DmajD'}.
\end{enumerate}
\end{proposition}
\begin{proof}
(a) Clearly $D_N(x,y)\ge 0$.  As $D(x,y)=D(y,x)$ it follows that
\begin{equation*}
\ D(z_0,z_1)+\cdots+D(z_N,z_{N+1})=D(z_{N+1},z_N)+\cdots+D(z_1,z_0).
\end{equation*}
Hence $D_N(x,y)=D_N(y,x)$.  Assume that $y=x$.  Choose $z_1=\cdots=z_N=x$.
As $D(x,x)=0$ we deduce that $\sum_{i=0}^N D(z_i,z_{i+1})=0$.  Hence $D_N(x,x)$=0.
As $D'$ is a metric we deduce
\begin{equation*}
\sum_{i=0}^N D'(z_i,z_{i+1})\ge D'(z_0, z_{N+1})=D'(x,y).
\end{equation*}
Use \eqref{DmajD'} to deduce that
\begin{eqnarray*}
\sum_{i=0}^N D(z_i,z_{i+1})\ge \sum_{i=0}^N D'(z_i,z_{i+1})\ge D'(x,y).
\end{eqnarray*}
Hence $D_N$ satisfies the inequality \eqref{DmajD'}.  In particular, if $x\ne y$ then
$D_N(x,y)\ge D'(x,y)>0$.  Therefore $D_N$ is a weak metric.

\noindent
(b)  Assume that $z_1=\ldots = z_{N}=x, z_{N+1}=y$.  Then
$\sum_{i=0}^N D(z_i,z_{i+1})=D(x,y)$.  Hence $D_N(x,y)\le D(x,y)$. Now let $z_{N+1}=z_{N+2}=y$.
Then
\begin{equation*}
\sum_{i=0}^N D(z_i,z_{i+1})=\sum_{i=0}^{N+1} D(z_i,z_{i+1}).
\end{equation*}
Hence $D_{N+1}(x,y)\le D_N(x,y)$.

\noindent
(c)  Choose $z_0 = x$, $z_{M + 1} = u$, $z_{M + N + 2} = y$, and $z_1,\ldots,z_{M+N+1}$ arbitrarily.  Then $\sum_{i=0}^{M+N+1} D(z_i,z_{i+1})\ge D_{M+N+1}(x,y)$.
Compare that with the definitions of $D_{M}(x,u)$ and $D_{N}(u,y)$ to deduce the inequality $D_M(x,u)+D_N(u,y)\ge D_{M+N+1}(x,y)$.

\noindent
(d) As $\{D_N(x,y)\}$ is a nonincreasing sequence such that $D_N(x,y)\ge D'(x,y)$ we deduce that the limit $D_\infty(x,y)$ exists and $D(x,y)\ge D_\infty(x,y)\ge D'(x,y)$. Since $D_N(x,y)=D_N(y,x)$ it follows that $D_\infty(x,y)=D_\infty(y,x)$. Hence $D_\infty(x,y)\ge 0$ and equality holds if and only if $x=y$.  In the inequality $D_M(x,u)+D_N(u,x)\ge D_{M+N+1}(x,y)$ let $M=N\to\infty$ to deduce that $D_\infty$ satisfies the triangle inequality.  Hence $D_\infty$ is a metric.  The inequality  $D(x,y)\ge D_\infty(x,y)\ge D'(x,y)$ yields that $D_\infty$ is a maximum metric $D'$ that satisfies \eqref{DmajD'}.
\end{proof}
\begin{theorem}\label{kapTABprop}  Let $C\in\mathrm{S}(\mathcal{H}_n\otimes\mathcal{H}_n)$.  Then $\mathrm{T}^Q_C$ is a semi-distance on $\Omega_n\times \Omega_n$ if and only if $C$ is positive semidefinite and $\ker(C)=\mathcal{H}_S$.
Assume that $C$ is positive semidefinite and $\ker(C)=\mathcal{H}_S$.
Then $\sqrt{\mathrm{T}_C^Q}$ is a weak metric.  Furthermore,
for $\rho^A,\rho^B\in\Omega_n$ the following statements hold:
\begin{enumerate}[(a)]
\item $\mathrm{T}_{C}^Q(\rho^A,\rho^B)=\mathrm{T}_{C}^Q(\rho^B,\rho^A)$.
\item $\mathrm{T}_{C}^Q(\rho^A,\rho^B)\ge 0$.
\item $\mathrm{T}_{C}^Q(\rho^A,\rho^B)=0$ if and only if $\rho^A=\rho^B$.
\item$\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)\le \frac{1}{2}(1-\tr \rho^A\rho^B$).
Furthermore
\begin{equation}\label{QOTrankone}
\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)= \frac{1}{2}(1-\tr \rho^A\rho^B) \textrm{ if either }
\rho^A \textrm{ or } \rho^B \textrm{ is a pure state}.
\end{equation}
\item $\sqrt{\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)}$ is a distance on pure states.
\end{enumerate}
\end{theorem}
\begin{proof}  We first show the second part of the theorem.   Assume that $C$ is positive semidefinite and  vanishes exactly on symmetric matrices.

\noindent
(a)   As $S$ is an involution with the eigenspaces $\mathrm{S}^2\mathbb{C}^n$ and $\mathrm{A}^2\mathbb{C}^n$ corresponding to the eigenvalues $1$ and $-1$ respectively, and $C\mathrm{S}^2\mathbb{C}^n=0$, it follows that $SC=CS=-C$.  Hence $SC S^\dagger=C$.
The second equality in \eqref{Werid} yields that $S\Gamma^Q(\rho^A,\rho^B)S^\dagger=\Gamma^Q(\rho^B,\rho^A)$.  As $\tr C \rho^{AB}=\tr C S\rho^{AB}S^\dagger$ we deduce (a).

\noindent
(b)  Since $C\ge 0$, for any $\rho^{AB}\in \Omega_{n^2}$ we get that $\tr C\rho^{AB}\ge 0$.  This proves (b).

\noindent
(c) Suppose that $\rho^A=\rho^B=\rho$.   Consider the spectral  decomposition of $\rho$ given by \eqref{specdecrho}.
Then a purification of $\rho$ is
\begin{equation}\label{purifA}
R= \Big( \sum_{i=1}^n \sqrt{\lambda_i} |\mathbf{x}_i\rangle|\mathbf{x}_i\rangle \Big) \Big( \sum_{j=1}^n \sqrt{\lambda_i} \langle\mathbf{x}_j|\langle\mathbf{x}_j| \Big)\in\Omega_{n^2}.\end{equation}
Clearly $R\in \Gamma^Q(\rho,\rho)$.  As $X=\sum_{i=1}^n \sqrt{\lambda_i} |\mathbf{x}_i\rangle|\mathbf{x}_i\rangle$ is a symmetric matrix it follows that $C X=0$.  Hence $\tr CR=0$ and  $\mathrm{T}_{C}^Q(\rho,\rho)=0$.

Assume now that $\mathrm{T}_{C}(\rho^A,\rho^B)=0$.  Hence $\tr C\rho^{AB}=0$ for some $\rho^{AB}\in\Gamma^Q(\rho^A,\rho^B)$.  That is, the eigenvectors of $\rho^{AB}$  are symmetric matrices.  Therefore $\rho^{AB}=\sum_{j=1}^k p_j|\psi_j\rangle \langle \psi_j|$
where each $ |\psi_j\rangle $ is a symmetric matrix and $p_j>0$.   We claim that each $|\psi_j\rangle \langle \psi_j|$ is of the form \eqref{purifA}.  This is equivalent to the Autonne--Takagi  factorization theorem \cite[Corollary 4.4.4, part (c)]{HJ13}
that any symmetric $X\in\mathbb{C}^{n\times n}$ is of the form
\begin{eqnarray*}
X=\sum_{i=1}^n d_i|\mathbf{x}_i\rangle|\mathbf{x}_i\rangle=UDU^\top, \quad D=\diag(\mathbf{d}), \quad U\in \mathrm{U}(n),
\end{eqnarray*}
where the columns of $U$ represent vectors,  $\mathbf{x}_1,\ldots,\mathbf{x}_n$.
Clearly $\tr_A |\psi_j\rangle \langle \psi_j|=\tr_B |\psi_j\rangle \langle \psi_j|$.
Hence $\rho^B=\tr_A \rho^{AB}=\tr_B \rho^{AB}=\rho^A$.

\noindent
(d)  As $\rho^A\otimes \rho^B\in \Gamma^Q(\rho^A,\rho^B)$ it follows the $\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)\le \tr C^{Q} (\rho^A\otimes\rho^B)$.  Clearly $\tr \mathbb{I}(\rho^A\otimes\rho^B)=1$.  The first part of \eqref{Werid} yields that $\tr S(\rho^A\otimes\rho^B)=\tr(\rho^A\rho^B)$.
Hence
$\tr C^{Q}(\rho^A\otimes\rho^B)=\frac{1}{2}\big(1-\tr \rho^A\rho^B\big)$, and
$\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\le \frac{1}{2}\big(1-\tr \rho^A\rho^B\big)$. Assume that either $\rho^A$ or $\rho^B$ is a pure state.  Lemma \ref{rangecont}  yields that $\Gamma^Q(\rho^A \rho^B)=\{\rho^A\otimes\rho^B\}$.  Hence \eqref{QOTrankone} holds.

\noindent
(e) It is known that if $\rho^A,\rho^B$ are pure state then \cite{Ren}
\begin{equation}\label{psiden}
\begin{aligned}
& \sqrt{1-\tr \rho^A\rho^B}=\frac{1}{2}\|\rho^A-\rho^B\|_1, \\
& \rho^A=|\mathbf{x}\rangle\langle \mathbf{x}|,\; \rho^B=|\mathbf{y}\rangle\langle \mathbf{y}|,\quad \langle \mathbf{x}|\mathbf{x}\rangle=\langle \mathbf{y}|\mathbf{y}\rangle=1.
\end{aligned}
\end{equation}
(Observe that $\sqrt{1-\tr \rho^A\rho^B}$ is the root infidelity if one of the states is pure.)
We give a short proof for completeness.
By changing the orthonormal basis in $\mathcal{H}_n$ we can assume that $n=2$ and
\begin{equation*}
\rho^A=\begin{bmatrix}1&0\\mathbf{0}&0\end{bmatrix}, \quad \rho^B=\begin{bmatrix}b&c\\c&1-b\end{bmatrix},\quad 0\le b\le 1,\,0\le  c, \,c^2=b(1-b).
\end{equation*}
As $\tr(\rho^A-\rho^B)=0$ it follows that the two eigenvalues of $\rho^A-\rho^B$
are
\begin{equation*}
\pm \sqrt{-\det(\rho^A-\rho^B)}=\pm \sqrt{(1-b)^2 +c^2}=\pm\sqrt{1-b}=\pm\sqrt{1-\tr\rho^A\rho^B}.
\end{equation*}
This proves \eqref{psiden}.  Hence
$\frac{1}{2}\|\rho^A-\rho^B\|_1+\frac{1}{2}\|\rho^B-\rho^C\|_1\ge\frac{1}{2}\|\rho^A-\rho^C\|_1$.  Combine that with (d) to deduce (e).

We now show the first part of the theorem.  Suppose that $C$ is positive semidefinite and  vanishes exactly on symmetric matrices. Then parts (a)-(c) of the theorem
show that  $\mathrm{T}_{C}^Q$ is a semi-distance.  Next observe that $C\ge aC^Q$ for some $a>0$.
  Hence $\mathrm{T}_C^Q(\rho^A,\rho^B)\ge a\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)$.  The inequality \eqref{YZYYin} proven in \cite{YZYY19} yields that
\begin{equation*}
\sqrt{\mathrm{T}_C^Q(\rho^A,\rho^B)}\ge \sqrt{a\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)}\ge D'(\rho^A,\rho^B)= \sqrt{a\frac{1-\sqrt{F(\rho^A,\rho^B)}}{2}},
\end{equation*}
where $F$ is the quantum fidelity \eqref{fidelity}.
As $D'$ is a
scaled Bures distance \cite{GLN05}, we deduce that $\sqrt{\mathrm{T}_C^Q}$ is a weak metric.

Assume now that $C\in\mathrm{S}(\mathcal{H}_n\otimes\mathcal{H}_n)$ and $\mathrm{T}_C^Q$ is a semi-distance.
For $n=1$ it is straightforward to see that $C=0$.  Assume that $n>1$.
As $\mathrm{T}_C^Q(\rho^A,\rho^B)>0$ for $\rho^A\ne \rho^B\in\Omega_n$ it follows that $C\ne 0$.   Let $R\in\mathrm{S}(\mathcal{H}_n\otimes \mathcal{H}_n)$ be nonzero and positive semidefinite. We claim that $\tr CR\ge 0$.  It is enough to assume that $\tr R=1$.  Set $\rho^A=\tr_B R, \rho^B=\tr_A R$.  Then $R\in \Gamma^Q(\rho^A,\rho^B)$.  Thus $0\le \mathrm{T}_C^Q(\rho^A,\rho^B)\le \tr CR$.  Suppose that $C=\sum_{k=1}^{n^2}\mu_k|\psi_k\rangle\langle \psi_k|$, where $|\psi_1\rangle,\ldots,|\psi_{n^2}\rangle$ is an orthonormal basis for $\mathcal{H}_n\otimes\mathcal{H}_n$.
Choose rank-one $R_k=|\psi_k\rangle\langle \psi_k|\ge 0$.  Thus $\mu_k=\tr CR_k\ge 0$ for $k\in[n^2]$.  Hence $C\ge 0$.
Let $\rho=|\mathbf{x}\rangle \langle\mathbf{x}|$ be a pure state.  Lemma \ref{rangecont} yields that $\Gamma^Q(\rho,\rho)=\{\rho\otimes\rho\}$.  Hence
$0=\mathrm{T}_C^Q(\rho,\rho)=\tr C(\rho\otimes\rho)$.
Noting that $\rho\otimes \rho=(|\mathbf{x}\rangle|\mathbf{x}\rangle)(\langle \mathbf{x}|\langle \mathbf{x}|)$,  as $C$ is positive semidefinite we deduce that $C(|\mathbf{x}\rangle|\mathbf{x}\rangle)=0$.  So $C$ vanishes on all rank one symmetric matrices, hence $C\mathcal{H}_S=0$.

It is left to show that $C|Y\rangle \ne 0$ if $Y$ is a nonzero skew-symmetric matrix.
Assume to the contrary that $C|Y\rangle =0$ for some nonzero skew-symmetric matrix $Y$.
Let $Z\in \mathrm{S}^2\mathbb{C}^n$ be the unique symmetric matrix with zero diagonal  such that $X=Z+Y$ is a nonzero lower triangular matrix with zero diagonal.  Note that $C|X\rangle =0$.  Normalize $X$ such that $\tr X X^\dagger=1$.
Let $R=|X\rangle\langle X|$, $\rho^A=\tr_B R, \rho^B =\tr_A R\in\Omega_n$.
Clearly $\tr CR=0$.  Hence $0\le \mathrm{T}_C^Q(\rho^A,\rho^B)\le \tr CR=0$.  As $\mathrm{T}_C^Q$ is a semi-distance we deduce that $\rho^A=\rho^B$.  We now contradict this equality.
Indeed, consider the equality \eqref{Xrhosigrel}.  As $X$ is lower triangular with zero diagonal its first row is zero.  Hence $\rho^A_{11}=0$.  Hence $\rho^B_{11}=0$.  Note that $\rho_{11}^B$ is the norm squared of the first column of $X$.  Hence the first column of $X$ is zero.  Therefore the second row of $X$ is zero.
Thus $\rho^A_{22}=0$, which yields that $\rho^B_{22}=0$.  Therefore the second column of $X$ is zero.  Repeat this argument to deduce that $X=0$ which contradicts our assumption that $\tr X X^\dagger=1$.
\end{proof}

\begin{definition}\label{defWasQ}  For a positive semidefinite $C$ with $\ker C=\mathcal{H}_S$
we define the metric \eqref{defQOTmet} induced by $\sqrt{\mathrm{T}_C^Q}$ as the
quantum Wasserstein-2 metric,
and denote it by  $W_C^Q(\rho^A,\rho^B)$.
\end{definition}
The key problem concerning the quantum Wasserstein-2 metric
is how to compute it.
If $\sqrt{T_{C^Q}^Q}$ is a metric then $W_{C^Q}=\sqrt{T_{C^Q}^Q}$, and in this case $W_{C^Q}$ can be computed within $\varepsilon$ precision in polynomial time.

We now give a variation of the inequality stated in part (d) of Theorem  \ref{kapTABprop}.
We start with the following (whose first part is well known \cite{Frb16}):
\begin{proposition}\label{purlem}
Assume that a normalized  $|\psi\rangle\in \mathcal{H}_n\otimes \mathcal{H}_n$ has Schmidt decomposition
\begin{align*}
|\psi\rangle=\sum_{i=1}^n \sqrt{\lambda_i}|\mathbf{x}_i\rangle |\mathbf{y}_i \rangle, && \lambda_1\ge \cdots\ge \lambda_n\ge 0, && \sum_{i=1}^n\lambda_i=1, && \langle \mathbf{x}_i,\mathbf{x}_j\rangle= \langle \mathbf{y}_i,\mathbf{y}_j\rangle=\delta_{ij}.
\end{align*}
Then $\tr_B |\psi\rangle\langle\psi|=\rho^A,\tr_A= |\psi\rangle\langle\psi|=\rho^B$, where $\rho^A$ and $\rho^B$ are two isospectral density matrices that are given by \eqref{specdecrho}.
Furthermore,
\begin{equation}\label{trCqupurrs}
\begin{aligned}
&\tr S |\psi\rangle\langle\psi| =\frac{1}{4}\Big(\|\sum_{i=1}^n\sqrt{\lambda_i}( |\mathbf{x}_i\rangle|\mathbf{y}_i\rangle+|\mathbf{y}_i\rangle|\mathbf{x}_i\rangle)\|^2-\|\sum_{i=1}^n \sqrt{\lambda_i}(|\mathbf{x}_i\rangle|\mathbf{y}_i\rangle-|\mathbf{y}_i\rangle|\mathbf{x}_i\rangle)\|^2\Big),\\
&\tr C^Q|\psi\rangle\langle\psi| = \frac{1}{4}\|\sum_{i=1}^n\sqrt{\lambda_i}( |\mathbf{x}_i\rangle|\mathbf{y}_i\rangle-|\mathbf{y}_i\rangle|\mathbf{x}_i\rangle)\|^2.
\end{aligned}
\end{equation}
\end{proposition}
\begin{proof}
Let us view $|\psi\rangle$ as a matrix $X\in\mathbb{C}^{n\times n}$.  Recall that $S$ is a selfadjoint involution with eigenvalue $1$ on the subspace of symmetric matrices and with eigenvalue $-1$ on the subspace of skew-symmetric matrices.  Moreover the orthogonal decomposition of $X$ is $(1/2)(X+X^\top)+(1/2)(X-X^\top)$, which corresponds to
\begin{equation*}
\sum_{i=1}^n \sqrt{\lambda_i}|\mathbf{x}_i\rangle|\mathbf{y}_i\rangle=\frac{1}{2}\Big(\sum_{i=1}^n \sqrt{\lambda_i}(|\mathbf{x}_i\rangle|\mathbf{y}_i\rangle+|\mathbf{y}_i\rangle|\mathbf{x}_i\rangle) +\sum_{i=1}^n \sqrt{\lambda_i}(|\mathbf{x}_i\rangle|\mathbf{y}_i\rangle-|\mathbf{y}_i\rangle|\mathbf{x}_i\rangle)\Big).
\end{equation*}
This gives the first part of \eqref{trCqupurrs}. The second part of  \eqref{trCqupurrs}
follows from the first part.
\end{proof}
Observe that the second part of of \eqref{trCqupurrs} gives an upper bound on
$\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)$ for isospectral $\rho^A,\rho^B$:
\begin{equation*}
\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\le \frac{1}{4}\|\sum_{i=1}^n\sqrt{\lambda_i}( |\mathbf{x}_i\rangle|\mathbf{y}_i\rangle-|\mathbf{y}_i\rangle|\mathbf{x}_i\rangle)\|^2.
\end{equation*}
However, this upper bound is not tight.  Indeed, if $|\psi\rangle $ corresponds to a skew symmetric matrix then this upper bound is $2$, while part (d) of Theorem \ref{kapTABprop} yields that $\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\le 1$.

The following lemma seems to be an improvement of part (d) of Theorem \ref{kapTABprop} for the case where $\rho^A$ and $\rho^B$ are isospectral:
\begin{lemma}\label{ubdkapAB}
Let $\rho^A,\rho^B\in\Omega_n$ be isospectral, with the spectral decompositions \eqref{specdecrho}.
Then
\begin{equation*}
\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\le \frac{1}{2}\Big(1-\sum_{i=1}^n  \lambda_i|\langle \mathbf{x}_i|\mathbf{y}_i\rangle|^2\Big).
\end{equation*}
Equality holds if $\rho^A$ and $\rho^B$ are pure states.
\end{lemma}
\begin{proof}
Set $\rho^{i,A}=|\mathbf{x}_i\rangle\langle\mathbf{x}_i|, \rho^{i,B}=|\mathbf{y}_i\rangle\langle\mathbf{y}_i|$.  Then
part  (d) of Theorem \ref{kapTABprop} yields that $\mathrm{T}_{C^{Q}}^Q(\rho^{i,A},\rho^{i,B})=\frac{1}{2}\big(1-|\langle\mathbf{x}_i|\mathbf{y}_i\rangle|^2\big)$.  The convexity of $\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)$ yields
\begin{equation*}
\sum_{i=1}^n \lambda_i\mathrm{T}_{C^{Q}}^Q(\rho^{i,A},\rho^{i,B})= \sum_{i=1}^n \frac{1}{2}\big(\lambda_i(1-|\langle\mathbf{x}_i|\mathbf{y}_i\rangle|^2)\big)\ge \mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B).  \qedhere
\end{equation*}
\end{proof}
Note that if $\rho^A=\rho^B$, we can take $\mathbf{y}_i=\mathbf{x}_i$.  Then the upper estimate in Lemma \ref{ubdkapAB} is $0$.
Thus if $\rho^A$ and $\rho^B$ are close one can choose the spectral decompositions of $\rho^A$ and $\rho^B$ such that the the upper estimate in Lemma \ref{ubdkapAB} is close $0$.

We now give a very general metric on positive semidefinite matrices, inspired by our lower bound on  $\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)$, which is exact on qubit density matrices.
\begin{proposition}\label{metricdenmat}  Let $\nu:\mathbb{R}^n \to [0,\infty)$ be a norm.  Assume that $f:[0,\infty)\to [0,\infty)$ is a continuous,  strictly increasing function.  For $\rho^A,\rho^B$ positive semidefinite define
\begin{equation}\label{defDrhosigma}
\begin{aligned}
D(\rho^A,\rho^B)=
\max_{U\in\mathrm{U}(n)}
\nu\Big( & \big( f((U^\dagger\rho^A U)_{11}),\ldots,f((U^\dagger\rho^A U)_{nn}) \big)^\top
\\
& - \big(f((U^\dagger\rho^B U)_{11}),\ldots,f((U^\dagger\rho^B U)_{nn}) \big)^\top\Big).
\end{aligned}
\end{equation}
Then $D(\rho^A,\rho^B)$ is a metric on positive semidefinite matrices.  In particular,
\begin{equation}\label{defD0rhosigma}
D_0(\rho^A,\rho^B)=
\max_{U\in\mathrm{U}(n)} \big|f((U^\dagger\rho^A U)_{11})-f((U^\dagger\rho^B U)_{11}) \big|
\end{equation}
is a metric on positive semidefinite matrices.
\end{proposition}
\begin{proof}  By definition $D(\rho^A,\rho^B)=D(\rho^B,\rho^A)\ge 0$.
Assume that $D(\rho^A,\rho^B)=0$.  Then $f((U^\dagger\rho^A U)_{ii})=f((U^\dagger\rho^B U)_{ii})$ for each $i\in[n]$ and $U\in U(n)$.  As $f$ is strictly increasing we deduce that $(U^\dagger\rho^A U)_{ii}=(U^\dagger\rho^B U)_{ii}$ for $i\in[n]$.  That is for each $U\in \mathrm{U}(n)$ the diagonal entries of $U^\dagger(\rho^A-\rho^B) U$ are 0.  Choose a unitary $V$ so that $V^\dagger(\rho^A-\rho^B)V$ is diagonal.  Then $V^\dagger(\rho^A-\rho^B)V=0$.  Hence $\rho^A=\rho^B$.  It is left to show the triangle inequality.

Denote by $\mathbf{f}(\rho)$ the vector $(f(\rho_{11}),\ldots,f(\rho_{nn}))^\top$.
Since $f$ is continuous there exists $V\in U(n)$ such that
$D(\rho^A,\rho^B)=\nu\big(\mathbf{f}(V^\dagger\rho^A V)-\mathbf{f}(V^\dagger\rho^B V)\big)$.  Hence
\begin{align*}
D(\rho^A,\rho^B) & = \nu\big(\mathbf{f}(V^\dagger\rho^A V)-\mathbf{f}(V^\dagger\rho^B V)\big)\\
 & \le \nu\big(\mathbf{f}(V^\dagger\rho^A V)-\mathbf{f}(V^\dagger\rho^C V)\big)+\nu\big(\mathbf{f}(V^\dagger\rho^C V)-\mathbf{f}(V^\dagger\rho^B V)\big) \\
 & \le D(\rho^A,\rho^C)+D(\rho^C,\rho^B).
\end{align*}

To show that $D_0(\cdot,\cdot)$ is a metric we observe that $D_0(\rho^A,\rho^B)=D(\rho^A,\rho^B)$ where $\nu\big((x_1,\ldots,x_n)^\top\big)=\max_{i\in[n]}|x_i|$.
Indeed, let $\mathrm{P}_n\subset \mathrm{U}(n)$ denote the group of permutation matrices.  Then
\begin{align*}
& \max_{i\in[n]} \big|f((U^\dagger\rho^A U)_{ii}) - f((U^\dagger\rho^B U)_{ii})\big| \\
& \qquad\qquad\qquad = \max_{P\in \mathrm{P}_n} \big| f \big(((UP)^\dagger\rho^A (UP))_{11} \big)-f \big(((UP)^\dagger\rho^B (UP))_{11} \big)\big|. \qquad \qedhere
\end{align*}
\end{proof}
\section{Comparison of classical and quantum optimal transports for diagonal density matrices}\label{sec:diagdm}
\begin{lemma}\label{diagdecr} Assume that $\rho^A,\rho^B\in\Omega_n$ and $C_E^Q$ is defined by \eqref{defCEop}.
Then
\begin{equation}\label{diagdecr1}
 \mathrm{T}_{C^Q_E}^Q(\diag(\rho^A),\diag(\rho^B))\le\mathrm{T}_{C^Q_E}(\rho^A,\rho^B).
\end{equation}
\end{lemma}
\begin{proof}  Without loss of generality we can assume that the basis
$|1\rangle,\ldots,|n\rangle$ used in \eqref{defCEop} is the standard orthonormal basis in $\mathcal{H}_n=\mathbb{C}^n$.
Denote  by $\mathcal{D}_n\subset \mathbb{C}^{n\times n}$ the subgroup of diagonal matrices whose diagonal entries are $\pm 1$.  Note that $|\mathcal{D}_n|=2^n$ and $\mathcal{D}_n$ is a subgroup of unitary matrices.   Observe next that, for $D\in\mathcal{D}_n$,
\begin{align*}
(D\otimes D) |\psi_{ij}^-\rangle \langle\psi_{ij}^-| (D\otimes D)= |\psi_{ij}^-\rangle \langle\psi_{ij}| && \Rightarrow && (D\otimes D)C^Q_E(D\otimes D)=C^Q_E.
\end{align*}
Hence $\mathrm{T}_{C^{Q}_E}^Q(\rho^A,\rho^B)=\mathrm{T}_{C^{Q}_E}^Q(D\rho^A D,D\rho^B D)$ for each $D\in\mathcal{D}_n$.  Clearly,
\begin{equation*}
\diag(\rho^A)=2^{-n}\sum_{D\in\mathcal{D}_n} D\rho^A D, \qquad \diag(\rho^B)=2^{-n}\sum_{D\in\mathcal{D}_n} D\rho^B D.
\end{equation*}
Use the convexity of $\mathrm{T}_{C^{Q}_E}^Q(\rho^A,\rho^B)$ to obtain
\begin{equation*}
\mathrm{T}_{C^{Q}_E}(\diag(\rho^A),\diag(\rho^B))\le 2^{-n}\sum_{D\in\mathcal{D}_n} \mathrm{T}_{C^{Q}_E}^Q(D\rho^A D,D\rho^B D)=\mathrm{T}_{C^{Q}_E}^Q(\rho^A ,\rho^B ).\qedhere
\end{equation*}
\end{proof}

Assume that $\mathbf{p}^A\in\Pi_m,\mathbf{p}^B\in \Pi_n$.
The following lemma gives the  isomorphism of $\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$ to $\Gamma^Q_{de}(\diag(\mathbf{p}^A),\diag(\mathbf{p}^A))$  described in the Introduction.
Furthermore it describes special $\rho^{AB}\in\Gamma^Q(\diag(\mathbf{p}^A), \diag(\mathbf{p}^B))$ induced by $\mathbf{p}^{AB}\in \Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$.
\begin{lemma}\label{diaglemobs}  Let $\rho^A \in\Omega_m, \rho^B\in\Omega_n$ and assume that  $\mathbf{p}^A\in\Pi^m,\mathbf{p}^B\in\Pi_n$ are induced by the diagonal entries of $\rho^A,\rho^B$ respectively.  Then
\begin{enumerate}[(a)]
\item Each matrix $X=[x_{ip}]_{i\in[m],p\in[n]}\in \Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$ induces the following two matrices
\begin{flalign*}
&& R=[r_{(i,p)(j,q)}],\tilde R=[\tilde r_{(i,p)(j,q)}] \in\Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B)), \;\;\, i,j\in[m], p,q\in[n].
\end{flalign*}
The matrix $R$ is diagonal with $r_{(i,p)(i,p)}=x_{ip}$ for $i\in[m], \, p\in[n]$, and $ \tilde R- R$ is a matrix whose only possible nonzero entries are the entries $((i,p)(p,i))$ for $i,p\in[\min(m,n)]$ and  $i\ne p$ which are equal to $\sqrt{x_{ip}x_{pi}}$. Furthermore, $\rank \tilde R\le mn-\min(m,n)(\min(m,n)-1)/2$.
\item Each matrix $R=[r_{(i,p)(j,q)}]\in\Gamma^Q(\rho^A,\rho^B)$ induces the following two matrices: First,
$X=[x_{ip}]\in \Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$, where $x_{ip}=r_{(i,p)(i,p)}$ for $i\in [m],p\in[n]$.
Second, $\hat R\in \Gamma^Q(\diag(\rho^A),\diag(\rho^B))$, which is obtained by replacing the entries of $R$ at places $((i,p)(j,q))$ by zero unless either $((i,p)(j,q))= ((i,p)(i,p))$ for $i\in[m], p\in[n]$ or $((i,p)(j,q))= ((i,p)(p,i))$ for $i,p\in[\min(m,n)], i\ne p$.
\end{enumerate}
\end{lemma}
\begin{proof} (a) As $X\in \Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$ we deduce
\begin{equation*}
\sum_{j=1}^n x_{ij}=p^{A}_i,\,i\in[m], \qquad \sum_{i=1}^m x_{ij}=p^B_j, j\in[n].
\end{equation*}
Assume that $R$ is a diagonal matrix with $r_{(i,p)(i,p)}=x_{ip}$.  Use
\eqref{ptfrom} to deduce that  $R\in\Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$.

Consider now the matrix $\tilde R$.  In view of \eqref{ptfrom} we deduce that $\tr_B \tilde R=\diag(\mathbf{p}^A)$ and  $\tr_A \tilde R=\diag(\mathbf{p}^B)$.  It is left to show that $\tilde R$ is positive semidefinite.
Observe that $\tilde R$ is a direct sum of $\big(mn-\min(m,n)(\min(m,n)-1)\big)$ $1\times 1$ matrices and
$\big(\min(m,n)(\min(m,n)-1)/2\big)$  $2\times 2$ matrices: $[x_{ii}]$ for $i\in[\min(m,n)]$, $[x_{ip}]$ for $i\in[m],p\in[n],\max(i,p)>\min(m,n)$, and
\begin{equation}\label{defXij}
 X_{ip}= \begin{bmatrix}x_{ip}&\sqrt{x_{ip}x_{pi}}\\\sqrt{x_{ip}x_{pi}}&x_{pi}\end{bmatrix}, \quad \textrm{ for } 1\le i<p\le \min(m,n).
 \end{equation}
   As $X\ge 0$ each block is positive semidefinite

    and has rank at most $1$.  Hence $\rank \tilde R\le mn- \min(n,n)(\min(m,n)-1)/2$.

\noindent
(b) Assume that $R\in\Gamma^Q(\rho^A,\rho^B)$.
As $R$ is positive semidefinite we deduce that $r_{(i,p)(i,p)}\ge 0$.  The above arguments yield that
the matrix $X=[r_{(i,p)(i,p)}]\in\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$.  Observe next that $\hat R$ is a direct sum of $1\times 1$ and $2 \times 2$ matrices: $[r_{(i,i)(i,i)}]$ for $i\in[\min(m,n)]$,
$[r_{(i,p)(i,p)}]$ for $\max(i,p) >[\min(m,n)]$, and
\begin{equation}\label{defRij}
R_{ip}=\begin{bmatrix}r_{(i,p)(i,p)}&r_{(i,p)(p,i)}\\r_{(p,i)(i,p)}&r_{(p,i)(p,i)}\end{bmatrix}, \quad  \textrm{ for } 1\le i<p\le \min(m,n).
\end{equation}
Clearly all these $1\times 1$ and $2\times 2$ submatrices are principal submatrices of $R$.  As $R$ is positive semidefinite, each such submatrix is positive semidefinite.  Hence $\hat R$ is positive semidefinite.
Use \eqref{ptfrom} to deduce that $\tr_B \hat R
 \diag(\mathbf{p}^A)$ and
  $\tr_A \hat R
 \diag(\mathbf{p}^B)$.
\end{proof}
\begin{lemma}\label{clqotdiagdm}  Assume that $\mathbf{p}^A\in\Pi_m, \mathbf{p}^B\in\Pi_n$ are induced by the diagonal entries of  $\rho^A\in\Omega_m,\rho^B\in\Omega_n$ respectively.   Let $C=[C_{(i,p)(j,q)}]$ for $i,j\in[m],p,q\in[n]$ be a Hermitian matrix.  Define $C^{cl}=[C^{cl}_{ip}]$ by $C^{cl}_{ip}=C_{(i,p)(i,p)}$ for $i\in[m],p\in[n]$.   Let $\Gamma_{de}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B)) \subset \Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$ be the subset of diagonal matrices.  Define
\begin{equation*}
\mathrm{T}^Q_{C,de}(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))=\min_{R\in\Gamma^Q_{de}(\mathbf{p}^A,\mathbf{p}^B)}\tr CR.
\end{equation*}
Then
\begin{equation}\label{OTgeQOT}
\begin{aligned}
\mathrm{T}_{C^{cl}}^{cl}(\mathbf{p}^A,\mathbf{p}^B) & =\mathrm{T}_{C,de}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))=
\mathrm{T}_{\diag(C)}^Q(\rho^A,\rho^B) \\
& \ge \mathrm{T}_{C}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B)).
\end{aligned}
\end{equation}
Assume that $m\le n$, and $C^Q=[C_{(i,p)(j,q)}^Q]\in \mathrm{S}_+(\mathcal{H}_n\otimes \mathcal{H}_n)$.
Denote by $C^Q_{m,n}\in\mathrm{S}_+(\mathcal{H}_m\otimes\mathcal{H}_n)$ the submatrix of $C^Q$ whose entries are $C_{(i,p)(j,q)}^Q$ for $i,j\in[m], p,q\in[n]$.  Let $C^{cl}_{m,n}$ be the $m\times n$ nonnegative matrix induced by the diagonal entries of $C^Q_{m,n}$.
Then
\begin{align}\label{clqotdiagdm0}
& \mathrm{T}_{C^{cl}_{m,n}}^{cl}(\mathbf{p}^A,\mathbf{p}^B) =
\frac{1}{2}\min_{X\in\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)}\Big(\sum_{1\le i<p\le m} ( x_{ip}+x_{pi} ) \;\, +\sum_{\substack{1\le i \le m,\\ m+1\le p\le n}} x_{ip}\Big) , \\
& \mathrm{T}_{C^Q_{m,n}}^Q \big(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B) \big)  \notag\\
& \qquad\;\; =\frac{1}{2}\min_{X\in\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)}\Big(\sum_{1\le i<p\le m} \big(x_{ip}+x_{pi}-2\sqrt{x_{ip}x_{pi}} \big) \; +\sum_{\substack{1\le i \le m, \\m+1\le p\le n}} x_{ip}\Big). \notag
\end{align}
\end{lemma}
\begin{proof}  Let $X=[x_{ij}]\in\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B))$ correspond to a diagonal matrix

\noindent
$R\in \Gamma_{de}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$ as in Lemma \ref{diaglemobs}.  Then $\tr C^{cl}X^\top=\tr C R$. This shows the first equality in \eqref{OTgeQOT}.
To show the second equality in \eqref{OTgeQOT} observe that for $R\in \Gamma^Q(\rho^A,\rho^B)$ we have  $\tr \diag(C) R=\tr \diag(C)\diag(R)$.
Next observe that $\diag(R)\in\Gamma_{de}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$.
As $$\Gamma_{de}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))\subset \Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$$ we deduce the inequality
\begin{equation*}
\mathrm{T}_{C,de}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))\ge \mathrm{T}_{C}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B)).
\end{equation*}
The proof of \eqref{OTgeQOT} is complete.

We now show  \eqref{clqotdiagdm0}.  Let  $R\in \Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$.  Define $X\in \Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$ and $\hat R\in\Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$  as in part (b) of Lemma \ref{diaglemobs}.  Furthermore, let $\tilde R\in\Gamma^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$  be defined as in part (a) of Lemma \ref{diaglemobs}.
It is straightforward to show that
\begin{eqnarray*}
\tr \diag(C^Q_{m,n})R=\tr C^{cl}_{m,n} X^\top,\quad
\tr C^Q_{m,n} R=\tr C^Q_{m,n}\hat R.
\end{eqnarray*}
Use the equalities in \eqref{OTgeQOT} to deduce the first equality in  \eqref{clqotdiagdm0}.

We now show the second equality in  \eqref{clqotdiagdm0}.
 As each $R_{ip}$ in \eqref{defRij} is positive semidefinite we deduce
\begin{align*}
& \tr C^Q_{m,n}\hat R \\
& \quad = \frac{1}{2}\Big( \sum_{1\le i<p\le m} \big( r_{(i,p)(i,p)} +r_{(p,i)(p,i)}-2\Re r_{(i,p)(p,i)}\big) \: +\!\!\!\!
\sum_{\substack{1\le i\le m,\\ m+1\le p\le n}} r_{(i,p)(i,p)}\Big) \\
& \quad \ge \frac{1}{2}\Big( \sum_{1\le i<j\le m} \big( r_{(i,p)(i,p)}+ r_{(p,i)(p,i)}-2\sqrt{r_{(i,p)(i,p)}r_{(p,i)(p,i)}}\big) \: + \!\!\!\!
\sum_{\substack{1\le i\le m,\\ m+1\le p\le n}} r_{(i,p)(i,p)}\Big) \\
& \quad \ge \frac{1}{2}\Big(\sum_{1\le i<p\le m} \big( x_{ij}+x_{ji}-2\sqrt{x_{ij}x_{ji}}\big) \: + \!\!\!\!
\sum_{\substack{1\le i\le m,\\ m+1\le p\le n}} x_{ip}\Big)=\tr C^Q_{m,n}\tilde R.
\end{align*}
This establishes the second equality in \eqref{clqotdiagdm0}.
\end{proof}
Observe that \eqref{OTgeQOT} generalizes the result in \cite{CGP20}, which claims that the cost of quantum optimal transport is cheaper than the cost of the classical optimal transport.

On the set of rectangular matrices matrices $\mathbb{R}^{m\times n}$, where $m\le n$, define
\begin{equation}\label{deff(X)}
f(X)=\frac{1}{2}\Big( \sum_{1\le i<p\le m} \big( x_{ip}+x_{pi}-2\sqrt{x_{ip}x_{pi}}\big)+\sum_{\substack{1\le i\le m,\\ m+1\le p\le n}} x_{ip} \Big), \quad X=[x_{ip}]\in \mathbb{R}_+^{m\times n}.
\end{equation}
As the function $\sqrt{xy}$ is a concave function on $\mathbb{R}_+^2$ it follows that $f(X)$ is a convex function on $\mathbb{R}^{m\times n}_+$.  Hence $\mathrm{T}_{C^{Q}_{m,n}}^Q(\diag(\mathbf{p}^A),\diag(\mathbf{p}^B))$ is the minimum of the convex function $f(X)$ on $\Gamma^{cl}(\mathbf{p}^A,\mathbf{p}^B)$.  Therefore this minimum can be computed in polynomial time within precision $\varepsilon>0$.
\begin{remark}\label{remextCQE}  We remark that we can extend the second equality in \eqref{clqotdiagdm0} to $C^Q_E$, which is given by \eqref{defCEop}.
\end{remark}

Lemma 11 in \cite{YZYY19} shows that
\begin{equation}\label{YZYYdub}
\mathrm{T}_{C^Q}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))\le\frac{1}{2}\Big( \sum_{i=1}^n (\sqrt{s_i}-\sqrt{t_i})^2-\min_{j\in[n]} (\sqrt{s_j}-\sqrt{t_j})^2\Big), \quad \mathbf{s},\mathbf{t}\in\Pi_n.
\end{equation}
Moreover, Algorithm 1 in \cite{YZYY19} gives $X\in \Gamma^{cl}(\mathbf{s},\mathbf{t})$ such that $f(X)$ is bounded from above by the right hand side of \eqref{YZYYdub}.

We now show that for $n=2$ the inequality \eqref{YZYYdub} is sharp.
\begin{lemma}\label{QOTdiagform2} Assume that $\mathbf{s}=(s_1,s_2)^\top, \mathbf{t}=(t_1,t_2)^\top$ are two probability vectors.  Then
\begin{equation}\label{QOTfdiagorm2a}
\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=\frac{1}{2}
\begin{cases}
(\sqrt{s_1}-\sqrt{t_1})^2, \; \textrm{ if } s_2\ge t_1,\\
(\sqrt{s_2}-\sqrt{t_2})^2, \; \textrm{ if } s_2< t_1.
\end{cases}
\end{equation}
Furthermore
\begin{equation}\label{QOTfdiagorm2a}
\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=\frac{1}{2}\max\big((\sqrt{s_1}-\sqrt{t_1})^2,(\sqrt{s_2}-\sqrt{t_2})^2\big).
\end{equation}
\end{lemma}
\begin{proof}
Assume that $s_2\ge t_1$.  Then $A=\begin{bmatrix}0&s_1\\t_1&s_2-t_1\end{bmatrix}\in\Gamma^{cl}(\mathbf{s},\mathbf{t})$.  Therefore
\begin{equation*}
\mathrm{T}_{C^{Q}}^Q(\mathbf{s},\mathbf{t})\le \frac{1}{2}(t_1+s_1-2\sqrt{t_1s_1})=\frac{1}{2}(\sqrt{s_1}-\sqrt{t_1})^2.
\end{equation*}
If $s_1t_1=0$ then $\Gamma^{cl}(\mathbf{s},\mathbf{t})=\{A\}$, and $\mathrm{T}_{C^{Q}}^Q(\mathbf{s},\mathbf{t})= \frac{1}{2}(\sqrt{s_1}-\sqrt{t_1})^2$.

Assume that $s_1t_1>0$.  Then $\Gamma^{cl}(\mathbf{s},\mathbf{t})$ is an interval $[A,B]$.
Indeed, let $C=\begin{bmatrix}1&-1\\-1&1\end{bmatrix}$.  So $A+tC\in \Gamma_{cl}(\mathbf{s},\mathbf{t})$ for $t$ small and positive, and $B=A+t_0C$ for some $t_0>0$.
Let $g(t)=f(A+tC)$ for $t\in[0,t_0]$.  Recall that $g(t)$ is a convex function on $[0,t_0]$.
Observe next that
$$g'(0+)=\frac{1}{2}\big(-2+s_1^{-1/2}t_1^{1/2}+s_1^{1/2}t_1^{-1/2}\big)=\frac{1}{2}s_1^{-1/2}t_1^{-1/2}\big(\sqrt{s_1}-
\sqrt{t_1})^2\ge 0.$$
Hence $g(t)\ge g(0)$ for $t\in[0,t_0]$.

It is left to show that  $\big(\sqrt{s_1}-\sqrt{t_1})^2\ge \big(\sqrt{s_2}-\sqrt{t_2})^2$.  Let $x\in[0,1/2]$.  Observe that the function $\sqrt{1/2+x}+\sqrt{1/2-x}$ is strictly  decreasing on $[0,1/2]$.   Hence
\begin{align*}
&\sqrt{s_1}+\sqrt{s_2}\le \sqrt{t_1}+\sqrt{t_2} \iff \max(s_1,s_2)\ge \max(t_1,t_2),\\
&\sqrt{s_1}+\sqrt{s_2}\ge \sqrt{t_1}+\sqrt{t_2} \iff \max(s_1,s_2)\le \max(t_1,t_2).
\end{align*}

Suppose first that $s_2\ge t_2$.  Hence $s_2\ge \max(t_1,t_2)$, and $s_1=1-s_2\le1-t_2=t_1$.  Thus
\begin{equation*}
|\sqrt{s_1}-\sqrt{t_1}|=\sqrt{t_1}-\sqrt{s_1}\ge \sqrt{s_2}-\sqrt{t_2}=|\sqrt{s_2}-\sqrt{t_2}|.
\end{equation*}

Suppose second that $s_2<t_2$.  Hence $t_2\ge s_1>t_1$.  Thus $\max(t_1,t_2)\ge \max(s_1,t_1)$.  Hence
 \begin{equation*}
|\sqrt{s_1}-\sqrt{t_1}|=\sqrt{s_1}-\sqrt{t_1}\ge \sqrt{t_2}-\sqrt{s_2}=|\sqrt{s_2}-\sqrt{t_2}|.
\end{equation*}
This proves the lemma in the case $s_2\geq t_1$.  Similar arguments prove the lemma in the case $s_2<t_1$.
\end{proof}
\begin{theorem}\label{QOTrhodiag} Let $\rho^A,\rho^B\in\Omega_2$ be two commuting density matrices of the form
\begin{equation*}
\rho^A=U\diag(s,1-s)U^\dagger \quad \rho^B=U\diag(t,1-t)U^\dagger,  \quad s,t\in [0,1],
\end{equation*}
for some unitary $U$.
Then
\begin{equation}\label{lowbdQOT1}
\begin{aligned}
\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B) & =
\mathrm{T}_{C^{Q}}^Q(\diag(s,1-s), \diag(t,1-t)) \\
& = \frac{1}{2}\max( (\sqrt{s}-\sqrt{t})^2, (\sqrt{1-s}-\sqrt{1-t})^2).
\end{aligned}
\end{equation}
Furthermore, the quantity $\sqrt{\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)}$ is a distance on the set of commuting density matrices in $\Omega_2$.
\end{theorem}
\begin{proof}
The first equality in \eqref{lowbdQOT1} follows from Corollary \ref{corA2}.  The second equality in \eqref{lowbdQOT1} follows from
\eqref{QOTfdiagorm2a}.

Let $\mathcal{C}\subset \Omega_2$ be a variety of commuting matrices.  Then there exists a unitary $U\in\mathbb{C}^{2\times 2}$ such that $\mathcal{C}=U\mathcal{D} U^\dagger$, where $\mathcal{D}$ is the variety of diagonal density matrices in $\Omega_2$. In view of \eqref{invCqt} it is enough to show that $\sqrt{\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)}$ is a distance on $\mathcal{D}$.
As $\sqrt{\mathrm{T}_{C^{Q}}(\rho^A,\rho^B)}$ is a semi-distance we need to show the triangle inequality on $\mathcal{D}$.
Assume that $\mathbf{r}=(r_1,r_2)^\top, \mathbf{s}=(s_1,s_2)^\top, \mathbf{t}=(t_1,t_2)^\top$ are probability vectors.  Then
\begin{align*}
& \sqrt{\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{r}),\diag(\mathbf{t}))} =\frac{1}{\sqrt{2}}\max\big(|\sqrt{r_1}-\sqrt{t_1}|,|\sqrt{r_2}-\sqrt{t_2}|\big) \\
&  \qquad \le
\frac{1}{\sqrt{2}}\max\big(|\sqrt{r_1}-\sqrt{s_1}|+|\sqrt{s_1}-\sqrt{t_1}|,|\sqrt{r_2}-\sqrt{s_2}|+|\sqrt{s_2}-\sqrt{t_1}|\big) \\
&  \qquad \leq
\frac{1}{\sqrt{2}}\big[\max\big(|\sqrt{r_1}-\sqrt{s_1}|,|\sqrt{r_2}-\sqrt{s_2}|\big)+
\max\big(|\sqrt{s_1}-\sqrt{t_1}|,|\sqrt{s_2}-\sqrt{t_2}|\big)\big] \\
& \qquad =
\sqrt{\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{r}),\diag(\mathbf{s}))} +\sqrt{\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))}. \qedhere
\end{align*}
\end{proof}

We now give a lower bound for $\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)$ on $\Omega_2$ and we will show later that this lower bound is sharp.
\begin{lemma}\label{lowbdQgam}
Assume that $\rho^A,\rho^B\in\Omega_2$.
Then
\begin{equation}\label{lowbdQgam1}
\begin{aligned}
\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\ge
\mathrm{T}_0(\rho^A,\rho^B) = \frac12\max_{U\in\mathrm{U}(2)} \bigg(\sqrt{(U^\dagger\rho^AU)_{11}}-\sqrt{(U^\dagger\rho^B U)_{11}}\bigg)^2.
\end{aligned}
\end{equation}
If $\rho^A$ and $\rho^B$ commute
then the inequality is sharp.
Furthermore the quantity $\sqrt{\mathrm{T}_0(\rho^A,\rho^B)}$ is a distance on $\Omega_2$.
\end{lemma}
\begin{proof}  First recall the equality \eqref{CQinvQOT}.
 Combine that with
 Lemma \ref {diagdecr} and \eqref{QOTfdiagorm2a} to deduce:
\begin{align*}\label{lowbdQgam0}
\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B) & \ge\mathrm{T}_{C^{Q}}^Q(\diag(U^\dagger\rho^A U),\diag(U^\dagger\rho^B U)) \\
& \ge
\frac12\bigg(\sqrt{(U^\dagger\rho^AU)_{11}}-\sqrt{(U^\dagger\rho^B U)_{11}}\bigg)^2.
\end{align*}
This proves \eqref{lowbdQgam1}.

Assume that $\rho^A$ and $\rho^B$ commute.  Without loss of generality we can assume that $\rho^A$ and $\rho^B$ are diagonal.  Choose $U=I$ to deduce that
$\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\ge \big(\sqrt{\rho^A_{11}}-\sqrt{\rho^B_{11}}\big)^2$.  Now choose $U$ to be the permutation matrix $A=\begin{bmatrix}0&1\\mathbf{1}&0\end{bmatrix}$.  Then \eqref{lowbdQgam1} yields that $\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\ge\frac{1}{2} \big(\sqrt{\rho^A_{22}}-\sqrt{\rho^B_{22}}\,\big)^2$.
Hence
\begin{equation*}
\textstyle
\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\ge\frac12\max\Big[ \big(\sqrt{\rho^A_{11}}-\sqrt{\rho^B_{11}}\big)^2,\big(\sqrt{\rho^A_{22}}-\sqrt{\rho^B_{22}}\big)^2\Big].
\end{equation*}
Theorem \ref{QOTrhodiag} yields that we have equality in the above inequality.
Hence we have equality in \eqref{lowbdQgam1}.

Finally observe that $\sqrt{T_0(\rho^A,\rho^B)}$ is the quantity $D(\rho^A,\rho^B) $ given in \eqref{defD0rhosigma} on $\Omega_2$, where  $f(x)=\sqrt{x}$.  Proposition \ref{metricdenmat} yields that  $\sqrt{T_0(\rho^A,\rho^B)}$ is a distance.
\end{proof}

\section{Decoherence of the quantum cost matrix}
\label{sec:decoh}
Let us denote
\begin{equation}\label{defCqalpha}
C^{Q}_{\alpha}=\frac{1}{2}\begin{bmatrix}0&0&0&0\\mathbf{0}&1&-\alpha&0\\mathbf{0}&-\alpha&1&0\\mathbf{0}&0&0&0
\end{bmatrix}=\alpha C^Q+(1-\alpha)\diag(C^Q), \quad \alpha\in[0,1].
\end{equation}
Assume that $\mathbf{s}=(s_1,s_2)^\top, \mathbf{t}=(t_1,t_2)^\top$ are probability vectors.
Then the quantity $\mathrm{T}_{C^Q_{\alpha}}(\diag(\mathbf{s}),\diag(\mathbf{t}))$ describes a continuous decoherence from $\alpha=1$ to $\alpha=0$.  We will show that, as expected, this function of $\alpha$ is decreasing on $[0,1]$ and give an exact formula for it.
\begin{lemma}\label{Pstalphaform} Let $\mathbf{s},\mathbf{t}$ be two probability vectors in $\mathbb{R}^2$.
Assume that $0\le \alpha\le 1$ and denote
\begin{equation*}
f_\alpha(X)=\frac{1}{2} \big( x_{12}+x_{21}-2\alpha\sqrt{x_{12}x_{21}} \, \big), \quad X=[x_{ij}]\in  \Gamma_{cl}(\mathbf{s},\mathbf{t}).
\end{equation*}
Then
\begin{equation}\label{defPalphast}
\mathrm{T}^Q(\mathbf{s},\mathbf{t},\alpha)=\mathrm{T}^Q_{C^{Q}_{\alpha}} \big(\diag(\mathbf{s}),\diag(\mathbf{t}) \big)=\min_{X\in \Gamma^{cl}(\mathbf{s},\mathbf{t})} f_\alpha(X).
\end{equation}
Let $\mathrm{T}^Q(\mathbf{s},\mathbf{t},1)=\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s},\mathbf{t}))$ be given by \eqref{QOTfdiagorm2a}.
Assume that $\mathrm{T}^Q(\mathbf{s},\mathbf{t},1)=(\sqrt{s_i}-\sqrt{t_i})^2$.
If either $\min(s_i,t_i)=0$ or $\mathbf{s}=\mathbf{t}$ then
\begin{eqnarray*}
\mathrm{T}^Q(\mathbf{s},\mathbf{t},\alpha)= \mathrm{T}^Q(\mathbf{s},\mathbf{t},1) \textrm{ for all }\alpha\in[0,1].
\end{eqnarray*}
Otherwise $\mathrm{T}^Q(\mathbf{s},\mathbf{t},\alpha)$ is a strictly decreasing function for $\alpha\in[0,1]$ given by the formula
\begin{equation}\label{Pstalphaform1}
\mathrm{T}^Q(\mathbf{s},\mathbf{t},\alpha)=\frac{1}{2}\begin{cases}
\sqrt{1-\alpha^2}|s_i-t_i|, &  \textrm{ for } 0\le \alpha<\frac{2\sqrt{s_it_i}}{s_i+t_i},\\
\mathrm{T}^Q(\mathbf{s},\mathbf{t},1)+2(1-\alpha)\sqrt{s_it_i}, & \textrm{ for } \frac{2\sqrt{s_it_i}}{s_i+t_i}\le \alpha\le 1.
\end{cases}
\end{equation}
\end{lemma}
\begin{proof}  The equality \eqref{defPalphast} is deduced as the second equality in \eqref{clqotdiagdm0}.
Observe next that $C^{Q}_{\alpha}$ is positive semidefinite.  Hence $\mathrm{T}^Q(\mathbf{s},\mathbf{t},\alpha)\ge 0$.  Therefore for $\mathbf{s}=\mathbf{t}$ we choose $X=\mathbb{I}\in\Gamma^{cl}(\mathbf{s},\mathbf{t})$ to deduce from  \eqref{defPalphast} that $\mathrm{T}^Q(\mathbf{s},\mathbf{t},\alpha)=0$.  Assume that $\min(s_i,t_i)=0$.  Then $\Gamma^{cl}(\mathbf{s},\mathbf{t})=\{B\}$, where $B$ has one zero off-diagonal element, and $\mathrm{T}^Q(\mathbf{s},\mathbf{t},\alpha)= \mathrm{T}^Q(\mathbf{s},\mathbf{t},1)$.

Assume that $\min(s_i,t_i)>0$ and $\mathbf{s}\ne \mathbf{t}$.  Suppose first that
that $s_2\ge t_1$.    Then for $\alpha=1$ Eq. \eqref{QOTfdiagorm2a} yields that $\mathrm{T}^Q(\mathbf{s},\mathbf{t},1)=\frac{1}{2}(\sqrt{s_1}-\sqrt{t_1})^2$, i.e., $i=1$.  Thus $\min(s_1,t_1)>0$.  The proof of Lemma \ref{QOTdiagform2}  yields that the minimum of $f_1(X)$ is achieved at the matrix $A=\begin{bmatrix}0&s_1\\t_1&s_2-t_1\end{bmatrix}$, which is an extreme point of $\Gamma^{cl}(\mathbf{s},\mathbf{t})$.  As $s_1,t_1>0$ it follows that $\Gamma^{cl}(\mathbf{p},\mathbf{t})$ is an interval, where the the second extreme matrix is $C=\begin{bmatrix}\min(s_1,t_1)&s_1-\min(s_1,t_1)\\t_1-\min(s_1,t_1)&s_2-t_1+\min(s_1,t_1)\end{bmatrix}$.  Thus we can move from $A$ to the relative interior of $\Gamma_{cl}(\mathbf{s},\mathbf{t})$ by considering $A(x)=A+xB$,
 where where $B=\begin{bmatrix}1&-1\\-1&1\end{bmatrix}$ and $x>0$.
Denoting
\begin{equation*}
g_{\alpha}(x)=f_{\alpha}(A(x))=\frac{1}{2} \big(s_1+t_1-2x -2\alpha\sqrt{s_1-x}\sqrt{t_1-x} \big),
\end{equation*}
one obtains
\begin{equation*}
g_\alpha'(0+)=\frac{1}{2}\Big[-2+\alpha\Big(\frac{\sqrt{t_1}}{\sqrt{s_1}} +\frac{\sqrt{s_1}}{\sqrt{t_1}}\Big)\Big].
\end{equation*}
Hence this derivative is nonnegative for
$\alpha\ge \frac{2\sqrt{s_1t_1}}{s_1+t_1}$
and negative for $0\le \alpha<\frac{2\sqrt{s_1t_1}}{s_1+t_1}$.
As $g_{\alpha}(x)$ is convex on the interval $[0, \min(s_1,t_1)]$ we obtain that for $ \frac{2\sqrt{s_it_i}}{s_i+t_i}\le \alpha\le 1$ the minimum of $g_{\alpha}$ for $\frac{2\sqrt{s_1t_1}}{s_1+t_1}$
is achieved at $x=0$.  This proves the second part of \eqref{Pstalphaform1}.
So assume that $0\le \alpha<
\frac{2\sqrt{s_1t_1}}{s_1+t_1}$.  Clearly the minimum of $f_0(X)$ on $\Gamma^{cl}(\mathbf{s},\mathbf{t})$ is achieved at $A(\min(s_1,t_1))$. For $\alpha>0$ we have $g'_{\alpha}(\min(s_1,t_1)-)=\infty$.  Hence for $0< \alpha<\frac{2\sqrt{s_1t_1}}{s_1+t_1}$ the minimum $g_{\alpha}(x)$ is achieved at a critical point $x\in (0, \min(s_1,t_1))$.  This critical point is unique, as $g_{\alpha}(x)$ is strictly convex on $(0,\min(s_1,t_1))$ and satisfies the quadratic equation
\begin{equation}\label{defxalpha}
4(s_1-x)(t_1-x)-\alpha^2(s_1+t_1-2x)^2=0, \quad  0\le \alpha<\frac{2\sqrt{s_1t_1}}{s_1+t_1}.
\end{equation}
We claim that the critical point is given by
\begin{equation*}
x(\alpha)=\frac{1}{2}\Big(s_1+t_1 -\frac{|s_1-t_1|}{\sqrt{1-\alpha^2}}\Big), \quad
0\le \alpha<\frac{2\sqrt{s_1t_1}}{s_1+t_1}.
\end{equation*}
A direct computation shows that $x(\alpha)$ satisfies \eqref{defxalpha}.
Next observe that as $s_1\ne t_1$ the function $x(\alpha)$ is a strictly decreasing function on $[0,1)$.  Clearly
\begin{equation*}
x(0)=\min(s_1,t_1), \quad
x \Big(\frac{2\sqrt{s_it_i}}{s_i+t_i} \Big)=0.
\end{equation*}
Hence $x(\alpha)\in (0,\min(s_1,t_1)]$.
Note that for $x(\alpha)$ we have equality
\begin{equation*}
2\sqrt{s_1-x(\alpha)}\sqrt{t_1-x(\alpha)}=\alpha\big(s_1+t_1 -2x(\alpha)\big).
\end{equation*}
This proves the first part of \eqref{Pstalphaform1} in the case for $i=1$.  Similar arguments show the first part of \eqref{Pstalphaform1} in the case for $i=2$.
Clearly for $s_i\ne t_i$ and $\min(s_i,t_i)>0$ the function $\mathrm{T}^Q(\mathbf{s},\mathbf{t},\alpha)$ is strictly decreasing on the interval $[0,1]$.
\end{proof}
\section{The dual problem}\label{sec:dualprob}
\begin{theorem}\label{dualQOT} Assume that $\rho^A\in\Omega_m,\rho^B\in\Omega_n$ and $C\in \mathrm{S}(\mathcal{H}_m\otimes \mathcal{H}_n)$.  Then the dual problem to \eqref{defkapCAB} is
\begin{equation}\label{dualQOT1}
\sup\{\tr \sigma^A  \rho^A +\tr \sigma^B  \rho^B, \; \sigma^A\in\mathrm{S}(\mathcal{H}_m), \: \sigma^B\in \mathrm{S}(\mathcal{H}_n), \:  C-\sigma^A\otimes \mathbb{I}_n-\mathbb{I}_m\otimes \sigma^B\ge 0\}.
\end{equation}
Furthermore, the above supremum is equal to $\mathrm{T}_{C}^Q( \rho^A, \rho^B)$.
Moreover, for $\rho^{AB}\in \Gamma^Q(\rho^A,\rho^B)$ and $F=C-\sigma^A\otimes \mathbb{I}_n-\mathbb{I}_m\otimes \sigma^B\ge 0$ the following complementary implication holds:
\begin{equation}\label{orthcond}
\tr F\rho^{AB}=0 \iff \tr C \rho^{AB}=\tr \sigma^A  \rho^A +\tr \sigma^B  \rho^B=\mathrm{T}_{C}^Q( \rho^A, \rho^B).
\end{equation}
In particular, if $\tr F\rho^{AB}=0$ then $\rank F\le mn -\rank \rho^{AB}$.

Assume that $ \rho^A, \rho^B> 0$. Then the above supremum is achieved:  There exist $\sigma^A\in\mathrm{S}(\mathcal{H}_m), \sigma^B\in \mathrm{S}(\mathcal{H}_n)$ such that
\begin{equation}\label{maxsig12}
\mathrm{T}_{C}^Q( \rho^A, \rho^B)=\tr (\sigma^A\rho^A+\sigma^B\rho^B),\quad C-\sigma^A\otimes \mathbb{I}_n-\mathbb{I}_m\otimes\sigma^B\ge  0.
\end{equation}
\end{theorem}
\begin{proof}
Let us first consider the simplified case where $ \rho^A, \rho^B,  C$ are real symmetric.  Let $\mathrm{S}_k\supset \mathrm{S}_{k,+}\supset \mathrm{S}_{k,+,1}$ be the space of $k\times k$ real symmetric matrices, the cone of positive semidefinite matrices and the convex set of real density matrices.
Define
\begin{align*}
& \Gamma^Q( \rho^A, \rho^B,\mathbb{R})=\mathrm{S}_{mn,+,1}\cap \Gamma^Q( \rho^A, \rho^B),\\
& \mathrm{T}_{C}^Q( \rho^A, \rho^B,\mathbb{R})=\min_{  \rho^{AB}\in \Gamma^Q( \rho^A, \rho^B,\mathbb{R})}\tr  C  \rho^{AB}.
\end{align*}
We claim that the dual problem to $\mathrm{T}_{C}^Q( \rho^A, \rho^B,\mathbb{R})$ is given by
\begin{equation}\label{dualQOT1R}
\sup\{\tr \sigma^A \rho^A +\tr \sigma^B \rho^B, \sigma^A\in\mathrm{S}_m, \sigma^B\in \mathrm{S}_n, C-\sigma^A\otimes \mathbb{I}_n-\mathbb{I}_m\otimes \sigma^B\ge 0\}.
\end{equation}
Indeed, the conditions  $\tr_B  \rho^{AB}= \rho^A, \tr_A  \rho^{AB}= \rho^B$ for $ \rho^{AB}\in \mathrm{S}_{mn,+}$ are stated as the linear conditions given by the first part of \eqref{margcond}.  Assume that $ \rho^A=[a_{ij}]\in \Omega_m,  \rho^B=[b_{ij}]\in\Omega_n$.
Then the standard dual characterization of the above semidefinite problem over   $\Gamma^Q( \rho^A, \rho^B,\mathbb{R})$ has the following form (see~\cite[Theorem 3.1]{VB96} or~\cite[(2.4)]{FriSDP}):
\begin{align*}
 \max \Big\{ & \sum_{1\le i \le j\le m} a_{ij} \tilde u_{ij} + \sum_{1\le p \le q\le n} b_{pq}\tilde v_{pq}, \quad \tilde u_{ij},\tilde v_{pq}\in\mathbb{R}, \\
& \Big(\sum_{1\le i \le j\le m} \tilde u_{ij}(G_{ij,m}\otimes \mathbb{I}_n) + \sum_{1\le p \le q\le n}\tilde v_{pq}(\mathbb{I}_m\otimes G_{pq,n})\Big) \le   C \Big\}.
\end{align*}
Let
\begin{eqnarray*}
\sigma^A=\sum_{1\le i \le j\le m} \tilde u_{ij}G_{ij,m}, \quad \sigma^B= \sum_{1\le p \le q\le n}\tilde v_{pq} G_{pq,n}.
\end{eqnarray*}
Then the last condition of the above maximum is $\sigma^A\otimes \mathbb{I}_n+ \mathbb{I}_m\otimes \sigma^B\le  C$.  Next observe that
\begin{eqnarray*}
\tr \sigma^A \rho^A+\tr \sigma^B \rho^B=
\bigl(\sum_{1\le i \le j\le m} a_{ij} \tilde u_{ij}\bigr) + \bigl(\sum_{1\le p \le q\le n} b_{pq}\tilde v_{pq}\bigr)
\end{eqnarray*}
Hence the dual to $\mathrm{T}_{C}^Q( \rho^A, \rho^B,\mathbb{R})$ is given by \eqref{dualQOT1R}.
Observe that we can choose $\sigma^A=-a\mathbb{I}_m, \sigma^B=0$, where $a$ is a positive big
number such that $$C -\sigma^A \otimes \mathbb{I}_n-\mathbb{I}_m\otimes\sigma^B=C+a\mathbb{I}_{mn}>0 .$$ Hence the duality theorem \cite[Theorem 3.1]{VB96} yields that the supremum \eqref{dualQOT1R} is equal  to $\mathrm{T}^Q_{C}( \rho^A, \rho^B,\mathbb{R})$.  Assume that $ \rho^A, \rho^B> 0$.  Then $0< \rho^A\otimes  \rho^B\in\Gamma^Q( \rho^A, \rho^B,\mathbb{R})$.  Theorem 3.1 in \cite{VB96} yields that the supremum \eqref{dualQOT1R} is achieved.

We now discuss the Hermitian case.  Let $\mathbf{i}=\sqrt{-1}$. There is a standard injective map $L:\mathrm{S}(\mathcal{H}_m) \to \mathrm{S}_{2m}$:
\begin{eqnarray*}
L(X+\mathbf{i} Y)=\begin{bmatrix} X&Y\\-Y&X\end{bmatrix}, \quad X,Y\in \mathbb{R}^{m\times m}, X^\top=X, Y^\top=-Y.
\end{eqnarray*}
Note that $L(X+\mathbf{i} Y) \ge 0\iff X+\mathbf{i} Y \ge 0$ and $L(X+\mathbf{i} Y)> 0\iff X+\mathbf{i} Y> 0$.
Hence it is possible to translate an SDP problem over Hermitian matrices to an SDP problem over reals.  This yields the proof that the supremum in \eqref{dualQOT1} is equal to $\mathrm{T}_{C}^Q( \rho^A, \rho^B)$.

Assume that  $\rho^{AB}\in \Gamma^Q(\rho^A,\rho^B)$ and $F=C-\sigma^A\otimes \mathbb{I}_n-\mathbb{I}_m\otimes \sigma^B\ge 0$.  As $\rho^{AB}$ and $F$ are positive semidefinite we obtain
\begin{eqnarray*}
0\le \tr F\rho^{AB}=\tr C \rho^{AB} -\tr \sigma^A\rho^A-\tr\sigma^B\rho^B
\end{eqnarray*}
The characterization \eqref{dualQOT1} yields the implication \eqref{orthcond}.
As $F$ and $\rho^{AB}$ are positive semidefinite the condition $\tr F\rho^{AB}=0$  yields that
$\rank F+\rank \rho^{AB}\le mn$.

Assume that $ \rho^A, \rho^B> 0$.  Then the above arguments show that the supremum in \eqref{dualQOT1} is achieved.
\end{proof}

We remark that the equality \eqref{dualQOT1} is stated in \cite[(4.2)]{CHW19}.

 In Subsection \ref{subsec:nonexisF} we give an example of $\rho^A,\rho^B\in\Omega_2$, where $\rho^A$ is a pure state,  for which the supremum  \eqref{dualQOT1} is not achieved.
Note that the dual problem has an advantage over the original problem, as we are not constrained by linear conditions \eqref{margcond}.   Also the number of variables is smaller, as the supremum is restricted to $\mathrm{S}(\mathcal{H}_m)\times \mathrm{S}(\mathcal{H}_n)$.  However we have to deal with the condition $  \sigma^A\otimes \mathbb{I}_n+\mathbb{I}_m\otimes  \sigma^B\le C$.

We now give a few applications of Theorem \ref{dualQOT}.
\subsection{The equality $\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)=\mathrm{T}_0(\rho^A,\rho^B)$ for qubits}\label{subsec:T=T0}
First we show that the inequality \eqref{lowbdQgam1} is sharp.
\begin{theorem}\label{mtqub}  Assume that $\rho^A,\rho^B\in\Omega_2$.  Then
\begin{equation}\label{mtqub1}
\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)=\frac12\max_{U\in\mathrm{U}(2)}\bigg(\sqrt{(U^\dagger\rho^A U)_{11}}-\sqrt{(U^\dagger\rho^B U)_{11}}\,\bigg)^2.
\end{equation}
\end{theorem}
\begin{proof}
First observe that $F$ that is given in \eqref{dualQOT}  is of the form:
\begin{equation}\label{sigmaFform}
\begin{aligned}
& \sigma^A=-\begin{bmatrix}a&b\\\bar b&c\end{bmatrix}, \quad \sigma^B=-\begin{bmatrix}e&f\\\bar f&g\end{bmatrix}, & a,c,e,g\in\mathbb{R}, \; b,f\in\mathbb{C},\\
& F \;\, = \; \begin{bmatrix}a+e&f&b&0\\\bar f&a+g+1/2&-1/2&b\\\bar b&-1/2&c+e+1/2&f\\mathbf{0}&\bar b&\bar f&c+g
\end{bmatrix}. &
\end{aligned}
\end{equation}

We now assume that $\rho^A, \rho^B$ are positive definite and non-isospectral.  Proposition \ref{exrank1} yields that $\Gamma^Q(\rho^A,\rho^B)$ does not contain a matrix of rank one.    Let $\rho^{AB}$ and $F$ be the matrices for which  \eqref{orthcond} holds.  Our assumptions yield that $\rank \rho^{AB}\ge 2$.  Proposition \ref{dualQOT} yields that
$\tr F\rho^{AB}=0$.  Hence $\rank F\le 4-2=2$.  Note that
the second and the third columns of $F$ are nonzero.  Hence $\rank F\ge 1$.

For $U\in\mathrm{U}(2)$ we have the equalities
\begin{align*}
&& \mathrm{T}_{C^Q}^Q(\rho^A,\rho^B) & = \mathrm{T}_{C^Q}^Q(U^\dagger\rho^A U,U^\dagger\rho^B U)=\tr \big(\sigma^A \rho^B+ \sigma^B\rho^B\big) \\
&&& = \tr \big((U^\dagger\sigma^A U)( U^\dagger\rho^B U)+(U^\dagger\sigma^B U)(U^\dagger\rho^A U)\big)
\end{align*}
\begin{align*}
\underline{ F} = (U^\dagger\otimes U^\dagger)F(U\otimes U^\dagger)=C^{Q}-(U^\dagger\sigma^A U)\otimes \mathbb{I}_2 -\mathbb{I}_2\otimes (U^\dagger\sigma^B U)\ge 0.
\end{align*}
We now choose $V\in\mathrm{U}(2)$ so that $V^\dagger\sigma^A V$ is a diagonal matrix. Let
\begin{equation*}\label{sigmaFform1}
\begin{aligned}
& \underline{\rho}^A=V^\dagger\rho^A V, \quad \underline{\rho}^B=V^\dagger\rho^B V, \\
& \underline{\sigma}^A=V^\dagger\sigma^A V=-\begin{bmatrix} \underline{a}&0\\mathbf{0}& \underline{c}\end{bmatrix}, \quad \underline{\sigma}^B=V^\dagger\sigma^B V=-\begin{bmatrix} \underline{e}& \underline {f}\\ \bar{\underline{f}} &\underline{g}\end{bmatrix}, \quad \underline{a},\underline {c},\underline {e},\underline {g}\in\mathbb{R}, \underline {f}\in\mathbb{C},\\
& \underline{F}=\begin{bmatrix}\underline{a}+\underline{e}&\underline{f}&0&0\\\bar {\underline{f}}&\underline{a}+\underline{g}+1/2&-1/2&0\\mathbf{0}&-1/2&\underline{c}+\underline{e}
+1/2&\underline{f}\\mathbf{0}&\bar 0&\bar {\underline{f}}&\underline{c}+\underline{g}
\end{bmatrix}.
\end{aligned}
\end{equation*}
Clearly $\rank \underline{F}=\rank F\le 2$.
We claim that $\rank \underline{F}=2$.  Assume to the contrary that $\rank \underline{F}=1$.  As the third column is nonzero we deduce that the fourth column is a multiple of the third column.  Hence the fourth column is zero.  That is, $\underline{f}=\underline{c}+\underline{g}=0$.  Similarly $\underline{a}+\underline{e}=0$.   Next observe that we can replace $\sigma^A,\sigma^B$ by $\sigma^A-\underline{a} \mathbb{I}_2, \sigma^B+\underline{a} \mathbb{I}_2$ without affecting the supremum in~\eqref{dualQOT1}.  This is equivalent to the assumption that $\underline{a}=0$.  Hence $\underline{e}=0$ and
$\underline{g}=-\underline{c}$.
As $\underline{F}$ is Hermitian and $\rank\underline F = 1$ we have the condition
\begin{eqnarray*}
0=(-\underline{c}+1/2)(\underline{c}+1/2)-1/4=-\underline{c}^2.
\end{eqnarray*}
Hence $\underline{c}=-\underline{g}=0$.  Thus
we can assume that $\sigma^A=\sigma^B=0$.  Equality \eqref{maxsig12} yields that $\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)=0$, which implies that $\rho^A=\rho^B$.  This contradicts our assumption that $\rho^A$ and $\rho^B$ are  not similar.  Hence $\rank \underline{F}=\rank F=2$.

We claim that either $\underline{x}=\underline{a}+\underline{e}$ or $\underline{z}=\underline{c}+\underline{g}$ are zero.  Assume to the contrary that $\underline{x},\underline{z}>0$.  (Recall that $\underline{F}> 0$.)  Let $\mathbf{c}_1,\mathbf{c}_2,\mathbf{c}_3,\mathbf{c}_4$ be the four columns of $\underline{F}$.  Clearly $\mathbf{c}_1,\mathbf{c}_4$ are linearly independent.  Hence $\mathbf{c}_2=u\mathbf{c}_1+v\mathbf{c}_4$.  As the fourth coordinate of $\mathbf{c}_2$ is zero we deduce that $v=0$.  Hence $\mathbf{c}_2=u\mathbf{c}_1$. This is impossible since  the third coordinate of $\mathbf{c}_1$ is $0$ and the third coordinate of $\mathbf{c}_2$ is $-1/2$.    Hence  either $\underline{x}=\underline{a}+\underline{e}$ or $\underline{z}=\underline{c}+\underline{g}$ are zero.
Suppose that $\underline{x}=0$.  As $\underline{F}$ is positive semidefinite  we deduce that the first row and column of  $\underline{F}$ is zero.  Hence $f=0$.
Similarly, if $\underline{z}=0$ we deduce that $f=0$.  Thus $\underline{\sigma}^A$ and $\underline{\sigma}^B$ are diagonal matrices.
Therefore
\begin{eqnarray*}
\mathrm{T}_{C^Q}^Q(\underline{\rho}^A,\underline{\rho}^B)=\tr\big(\underline{\sigma}^A\underline{\rho}^A+
\underline{\sigma}^B\underline{\rho}^B\big)=\tr\big(\underline{\sigma}^A
\diag({\underline{\rho}}^A)+
\underline{\sigma}^B\diag({\underline{\rho}}^B)\big).
\end{eqnarray*}
As $\underline{F}\ge 0$, the maximum dual characterization yields
\begin{align*}
\tr\big(\underline{\sigma}^A
\diag({\underline{\rho}}^A)+
\underline{\sigma}^B\diag({\underline{\rho}}^B)\big) \le \mathrm{T}_{C^{Q}}^Q(\diag({\underline{\rho}}^A),\diag({\underline{\rho}}^B)).
\end{align*}
Hence $\mathrm{T}_{C^Q}^Q(\underline{\rho}^A,\underline{\rho}^B)\le \mathrm{T}_{C^Q}^Q(\diag({\underline{\rho}}^A),\diag({\underline{\rho}}^B))$.  Compare that with \eqref{diagdecr1} to deduce the equalities
\begin{equation*}\label{4QOTeq}
\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)=\mathrm{T}_{C^Q}^Q(\underline{\rho}^A,\underline{\rho}^B)= \mathrm{T}_{C^Q}^Q(\diag({\underline{\rho}}^A),\diag({\underline{\rho}}^B)).
\end{equation*}
Use \eqref{QOTfdiagorm2a} to deduce
\begin{align*}
\mathrm{T}_{C^Q}^Q \underline{\rho}^A,\underline{\rho}^B) &= \mathrm{T}_{C^Q}^Q \big(\diag(\underline{\rho})^A,\diag({\underline{\rho}}^B) \big) \\
& =\frac{1}{2}\max\Big[ \big(\sqrt{\underline{\rho}^A_{11}}-
\sqrt{\underline{\rho}^B_{11}}\, \big)^2, \, \big(\sqrt{\underline{\rho}^A_{22}}-
\sqrt{\underline{\rho}^B_{22}}\, \big)^2\Big].
\end{align*}
The inequality \eqref{lowbdQgam1} yields the theorem
for $\rho^A$ and $\rho^B$ positive definite and non-isospectral.  Clearly every pair $\rho^A,\rho^B\in\Omega_2$ can be approximated by $\hat\rho^A,\hat\rho^B\in\Omega_2$ which are positive definite and non-isospectral.  Use the continuity of $T_{C^Q}^Q(\rho^A,\rho^B)$ on $\Omega_2\times\Omega_2$ (Proposition \ref{honconv})  to deduce the theorem in the general case.
\end{proof}

Combine this theorem with the last part of Lemma \ref{lowbdQgam} to deduce:
\begin{corollary}\label{trinqub}  The quantity $\sqrt{\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)}$ is a distance on $\Omega_2$.
\end{corollary}

\subsection{A semi-analytic formula for the single-qubit optimal transport}
\label{subsec:Blochb}
We now introduce a convenient notation for qubits in the $y=0$ section of the Bloch ball \cite[Section 5.2]{BZ17}. Let $O$ denote the rotation matrix
\begin{align*}
O(\theta)=\begin{bmatrix}\cos(\theta/2)&-\sin(\theta/2)\\\sin(\theta/2)&\cos(\theta/2)
\end{bmatrix}, \quad \text{for } \theta \in [0,2\pi),
\end{align*}
and define, for $r\in [0,1]$,
\begin{align*}
\rho(r,\theta) & = O(\theta) \begin{bmatrix}r&0\\mathbf{0}&1-r\end{bmatrix}O(\theta)^\top
\end{align*}
Because of unitary invariance \eqref{CQinvQOT}, the quantum transport problem between two arbitrary qubits $\rho^A, \rho^B \in \Omega_2$ can be reduced to the case $\rho^A = \rho(s,0)$ and $\rho^B = \rho(r,\theta)$, with three parameters, $s,r \in [0,1]$ and $\theta \in [0,2\pi)$. The parameter $\theta$ is the angle between the Bloch vectors associated with $\rho^A$ and $\rho^B$. With such a parametrization we can further simplify the single-qubit transport problem.

Observe first that if $s \in \{0,1\}$ then $\rho^A$ is pure, and if $r \in \{0,1\}$ then $\rho^B$ is pure. In any such case an explicit solution of the qubit transport problem is given \eqref{QOTrankone}.
\begin{theorem}\label{thmC3}
Let $\rho^A = \rho(s,0), \rho^B = \rho(r,\theta)$ and assume that $0<r,s<1$. Then
\begin{equation}\label{defT0}
\mathrm{T}^Q_{C^Q}(\rho^A,\rho^B) =
\max_{\phi\in \Phi(s,r,\theta)} \frac{1}{4}\Big(\sqrt{1+ (2s-1)\cos\phi} - \sqrt{1+(2r-1)\cos(\theta+\phi)}\Big)^2,\notag
\end{equation}
where $\Phi(s,r,\theta)$ is the set of all $\phi\in[0,2\pi)$ satisfying the equation
\begin{equation}\label{Phieq}
\frac{(2s-1)^2\sin^2\phi}{1+(2s-1)\cos\phi}=\frac{(2r-1)^2\sin^2(\theta+\phi)}{1+(2r-1)\cos(\theta+\phi)}.
\end{equation}
\end{theorem}
\begin{proof}
A unitary $2 \times 2$ matrix $U$ can be parametrized, up to a global phase, with three angles $\alpha, \beta, \phi \in [0,2\pi)$,
\begin{align*}
U = \begin{bmatrix} e^{\mathbf{i} \alpha} & 0\\ 0 & e^{-\mathbf{i} \alpha} \end{bmatrix} O(\phi) \begin{bmatrix} e^{\mathbf{i} \beta} & 0\\ 0 & e^{-\mathbf{i} \beta} \end{bmatrix}.
\end{align*}
Thus, setting $f(r,\theta;\alpha,\phi) = (U^{\dagger}\rho(r,\theta) U)_{11}$, we have
\begin{align*}
f(r,\theta;\alpha,\phi) =  \frac{1}{2} \Big( 1+ (2 r-1) \big( \cos (\theta ) \cos (\phi ) + \cos (2 \alpha ) \sin (\theta ) \sin (\phi ) \big) \Big).
\end{align*}
This quantity does not depend on the parameter $\beta$, so we can set $\beta = 0$. Note also that $f(s,0;\alpha,\phi)$ does not depend on $\alpha$. With $\rho^A = \rho(s,0), \rho^B = \rho(r,\theta)$, Theorem \ref{mtqub}  yields
\begin{align*}
\mathrm{T}^Q_{C^Q}(\rho^A,\rho^B) = \frac12\max_{\alpha,\phi\in [0,2\pi)} \Big( \sqrt{f(s,0;0,\phi)} - \sqrt{f(r,\theta;\alpha,\phi)} \Big)^2.
\end{align*}
Now, note that the equation $\partial_\alpha f(r,\theta;\alpha,\phi) = 0$ yields the extreme points $\alpha_0 = k \pi /2$, with $k \in \mathbb{Z}$. Since $f(r,\theta;\alpha + \pi,\phi) = f(r,\theta;\alpha,\phi)$ we can take just $\alpha_0 \in \{0,\pi/2\}$. Consequently,
\begin{align*}
\mathrm{T}^Q_{C^Q}(\rho^A,\rho^B) =
\max_{\phi \in [0,2\pi)} \{ g_-(s,r,\theta;\phi), g_+(s,r,\theta;\phi) \},
\end{align*}
where we introduce the auxilliary functions
\begin{align}\label{g}
 g_\pm(s,r,\theta;\phi) = \frac{1}{4}\Big(\sqrt{1+ (2s-1)\cos\phi} - \sqrt{1+(2r-1)\cos(\theta\pm\phi)}\Big)^2.
\end{align}
But since $g_-(s,r,\theta;2\pi - \phi) = g_+(s,r,\theta;\phi)$ we can actually drop the $\pm$ index in the above formula. In conclusion, we have shown that it is sufficient to take $U = O(\phi)$ for $\phi \in [0,2\pi)$ in formula \eqref{mtqub1}.

Finally, it is straightforward to show that the equation $\partial_\phi g(s,r,\theta;\phi) = 0$ is equivalent to \eqref{Phieq}. Hence, $\Phi(s,r,\theta)$ is the set of extreme points, and \eqref{defT0} follows.
\end{proof}

\begin{lemma}\label{6sollem}  The equation \eqref{Phieq} has at most six solutions $\phi\in[0,2\pi)$ for given $r,s\in(0,1), \theta\in[0,2\pi)$.  Moreover there is an open set of $s,r\in (0,1),\theta\in[0,2\pi)$ where there are exactly $6$ distinct solutions.
\end{lemma}
\begin{proof}
Write $z=e^{\mathbf{i}\phi}, \zeta=e^{\mathbf{i}\theta}$.  Then
\begin{align*}
& 2\cos\phi =z+\frac{1}{z}, &&  2\mathbf{i}\sin \phi=z-\frac{1}{z}, \\
& 2\cos(\theta+\phi)=\zeta z+\frac{1}{\zeta z}, && 2\mathbf{i}\sin(\theta+\phi)=\zeta z-\frac{1}{\zeta z}.
\end{align*}
Thus \eqref{Phieq} is equivalent to
\begin{multline}
(1-2 r)^2 \left[ (2 s-1) \left(z^2+1\right)+2 z\right] \left(\zeta ^2 z^2-1\right)^2 \label{z6}\\
 -\zeta  (1-2 s)^2 \left(z^2-1\right)^2 \left[ (2 r-1) \left(\zeta ^2 z^2+1\right)+2 \zeta  z\right] = 0.
\end{multline}
This a 6th order polynomial equation in the variable $z$, so it has at most 6 real solutions.  Since we must have $\vert z \vert = 1$, not every complex root of \eqref{z6} will yield a real solution to the original  \eqref{Phieq}. Nevertheless, it can be shown  that there exist open sets in the parameter space $s,r \in (0,1)$, $\theta \in [0,2\pi)$ on which \eqref{Phieq} does have 6 distinct solutions.

Observe that if $\theta=0$ and $s,r\in(0,1)$ and $s\ne r$  then two  solutions to the equality \eqref{Phieq} are $\phi\in\{0,\pi\}$, which means that $z=\pm 1$.
In this case the equality \eqref{Phieq} is
\begin{equation*}
\sin^2\phi \: \bigg(\frac{(2s-1)^2}{1+(2s-1)\cos\phi}-\frac{(2r-1)^2}{1+(2r-1)\cos(\phi)}\bigg)=0.
\end{equation*}
As $\sin^2\phi=-(1/4)z^{-2}(z^2-1)^2$ we see that $z=\pm 1$ is a double root.

Another solution $\phi\notin\{0,\pi\}$ is given by
\begin{equation*}
\cos\phi=\frac{(2s-1)^2-(2r-1)^2}{(2r-1)^2(2s-1)-(2r-1)(2s-1)^2}=\frac{2(1-r-s)}{(2r-1)(2s-1)}.
\end{equation*}
Assume that $r+s=1$.  Then $\cos\phi=0$, so $\phi\in\{\pi/2, 3\pi/2\}$.  Thus if $r+s$ is close to $1$ we have that $\phi$ has two values close to $\pi/2$ and $3\pi/2$ respectively.  Hence in this case we have $6$ solutions counting with multiplicities.

We now take a small $|\theta|>0$.  The two simple solutions $\phi$ are close to $\pi/2$ and $3\pi/2$.  We now need to show that the double roots $\pm 1$ split to two pairs of solutions on the unit disc: one pair close to $1$ and the other pair close to $-1$.
Let us consider the pair close to $1$, i.e., $\phi$ close to zero.  Then the equation \eqref{Phieq} can be written in the form
\begin{multline*}
(2s-1)^2\big(1+(2r-1)\cos(\theta+\phi)\big)\sin^2\phi\\
- (2r-1)^2\big(1+(2s-1)\cos\phi\big)\sin^2(\theta+\phi)=0.
\end{multline*}
Replacing $\sin\phi, \sin(\theta+\phi)$ by $\phi, \theta+\phi$ respectively we see that the first term gives the equation:
$(2s-1)^2(2r)\phi^2-(2r-1)^2 2s(\theta+\phi)^2=0$. Then we obtain two possible Taylor series of $\phi$ in terms of $\theta$:
\begin{align*}
\phi_1(\theta) & =\frac{(2r-1)\sqrt{s}\theta}{(2s-1)\sqrt{r}-(2r-1)\sqrt{s}} + \theta^2 E_1(\theta), \\
\phi_2(\theta) & =-\frac{(2r-1)\sqrt{s}\theta}{(2s-1)\sqrt{r}+(2r-1)\sqrt{s}}+\theta^2E_2(\theta).
\end{align*}
Use the implicit function theorem to show that $E_1(\theta)$ and $E_2(\theta)$ are analytic in $\theta$ in the neighborhood of $0$.
Hence in this case we have $6$ different solutions.
\end{proof}

We have thus shown that the general solution of the quantum transport problem of a single qubit with cost matrix $C^Q = \tfrac{1}{2} \big(\mathbb{I}_{4} - S\big)$ is equivalent to solving a 6th degree polynomial equation with certain parameters. For some specific values of these parameters an explicit analytic solution can be given. This is discussed in the next subsection.

\subsection{Two isospectral qubit density matrices}
\label{subsec:isospectral}
In view of unitary invariance \eqref{CQinvQOT} and the results of the previous section we can assume that  two isospectral qubits have the following form:
 $\rho^A = \rho(s,0)$ and $\rho^B = \rho(s,\theta)$ for some $s \in [0,1]$ and $\theta \in [0,2\pi)$.
\begin{theorem}\label{thmE1}
 For any $s \in [0,1]$ and $\theta \in [0,2\pi)$ we have
\begin{equation}\label{Tiso}
\mathrm{T}^Q_{C^Q} \big(\rho(s,0),\rho(s,\theta) \big) = \Big( \tfrac{1}{2} -\sqrt{s(1-s)} \Big) \sin^2 (\theta/2).
\end{equation}
\end{theorem}
\begin{proof}
Note first that if the states $\rho^A,\rho^B$ are pure, i.e. $s = 0$ or $s=1$, formula \eqref{Tiso} gives $\mathrm{T}^Q_{C^Q} \big(\rho(s,0),\rho(s,\theta) \big) = \tfrac{1}{2} \sin^2 (\theta / 2)$, which agrees with \eqref{QOTrankone}.

From now on we assume that that $\rho^A, \rho^B$ are not pure.  When $r = s$, \eqref{z6} simplifies to the following:
\begin{multline}
 (\zeta -1) (1-2 s)^2 \left(\zeta  z^2-1\right) \times  \\
\quad \times  \left[4 s (\zeta +1) \left(\zeta  z^2+1\right) z +(2 s-1) (z-1)^2 (\zeta  z-1)^2 \right]= 0. \label{zeta_iso}
\end{multline}

Eq.\ \eqref{zeta_iso} is satisfied when $z = \pm \zeta^{-1/2}$. This corresponds to $\phi_0 = -\theta/2$ or $\phi_0' = \pi - \theta/2$. Observe, however, that we have $g(s,s,\theta;\phi_0) = g(s,s,\theta;\phi_0') = 0$, so we can safely ignore $\phi_0, \phi_0' \in \Phi(s,s,\theta)$ in the maximum in  \eqref{defT0}.

Hence, we are left with a 4th order equation
\begin{align}\label{zeta4}
4 s (\zeta +1) \left(\zeta  z^2+1\right) z +(2 s-1) (z-1)^2 (\zeta  z-1)^2 = 0,
\end{align}
which
reads
\begin{equation}\label{phi4}
(2 s-1) \big[ 2 + \cos (\theta +2 \phi )+ \cos (\theta ) \big]
+2 \big[ \cos (\theta +\phi )+ \cos (\phi ) \big] = 0.
\end{equation}
Now, observe that if $\phi$ satisfies  \eqref{phi4}, then so does $\phi' = -\phi - \theta$. This translates to the fact that if $z$ satisfies  \eqref{zeta4}, then so does $(z \zeta)^{-1}$. Furthermore, $g(s,s,\theta;\phi) = g(s,s,\theta;\phi')$. Hence, in the isospectral case we are effectively taking the maximum over just two values of $\phi$.

Let us now seek an angle $\phi_1 \in [0,2\pi)$ such that $g(s,s,\theta;\phi_1)$ equals the righthand side of  \eqref{Tiso}. The latter equation reads
\begin{align*}
& \Big\{ (2 s-1) \big[\cos \left(\theta +\phi _1\right)+\cos \left(\phi _1\right)\big]
 -\big(2 \sqrt{s(1-s)}-1\big) \big(\cos (\theta )-1\big)+2\Big\}^2 \\
& \qquad\quad = 4 \big[(2 s-1) \cos \left(\phi _1\right)+1\big] \big[(2 s-1) \cos
   \left(\theta +\phi _1\right)+1\big].
\end{align*}
In terms of $z$ and $\zeta$, the above is equivalent to a 4th order polynomial equation in $z$, which can be recast in the following form:
\begin{align}\label{zeq}
\Big[ \zeta  (1-2 s) z^2+(\zeta +1) \big(2 \sqrt{s(1-s)}-1\big) z-2 s+1 \big]^2 = 0.
\end{align}
Hence,  \eqref{zeq} has two double roots:
\begin{multline*}
 z_1^{\pm} = \big[ 2 \zeta  (1-2 s) \big]^{-1} \bigg\{ (\zeta +1) \big( 1-2 \sqrt{s(1-s)} \, \big) \\
 \pm \sqrt{(\zeta +1)^2 \big(1-2 \sqrt{s(1-s)} \,\big)^2-4 \zeta  (1-2s)^2} \bigg\}.
\end{multline*}
Furthermore, one can check that $z_1^{-} = (\zeta z_1^{+})^{-1}$.

Now, it turns out that $z_1^{\pm}$ are also solutions to \eqref{zeta4}, as one can quickly verify using \textsc{Mathematica}~\cite{Mathematica}. We thus conclude that $\phi_1, \phi_1' \in \Phi(s,s,\theta)$.

We now divide the polynomial in \eqref{zeta4} by $(z-z_1^{+})(z-z_1^{-})$. We are left with the following quadratic equation
\begin{align*}
\zeta  \Big[ (2 s-1) \left(\zeta  z^2+1\right)+(\zeta +1) \big(2 \sqrt{(1-s) s}+1\big) z\Big] = 0.
\end{align*}
Its solutions are
\begin{multline*}
z_2^{\pm} = \big[ 2 \zeta  (1-2 s) \big]^{-1} \bigg\{ (\zeta +1) \big( 1+2 \sqrt{s(1-s)} \, \big) \\
\pm \sqrt{(\zeta +1)^2 \big(1+2 \sqrt{s(1-s)} \,\big)^2-4 \zeta  (1-2s)^2} \bigg\}.
\end{multline*}
Again, we have $z_2^{-} = (\zeta z_2^{+})^{-1}$, in agreement with the symmetry argument. Setting $z_2^+ = \vcentcolon e^{\mathbf{i} \phi_2}$ and $z_2^- = \vcentcolon e^{\mathbf{i} \phi_2'}$ we have $\phi_2, \phi_2' \in \Phi(s,s,\theta)$. Then we deduce, with the help of \textsc{Mathematica}, that
\begin{align*}
g(s,s,\theta;\phi_2) & = g(s,s,\theta;\phi_2')  = \tfrac{1}{4} \Big[ (1-6 \sqrt{(1-s) s} - \big(1+2 \sqrt{(1-s) s} \, \big) \cos (\theta ) \Big].
\end{align*}

Finally, we observe that
\begin{align*}
g(s,s,\theta;\phi_1) - g(s,s,\theta;\phi_2) = \sqrt{(1-s) s} \, \big(1+\cos (\theta ) \big) \geq 0.
\end{align*}
This shows that, for any $s \in (0,1)$, $\theta \in [0,2\pi)$,
\begin{align*}
\mathrm{T}^Q_{C^Q} \big(\rho(s,0),\rho(s,\theta) \big) = g(s,s,\theta;\phi_1),
\end{align*}
and \eqref{Tiso} follows.
\end{proof}

Note that $g(s,s,\theta;\phi_2)$ can become negative for certain values of $s$ and $\theta$. This means that for such values $\Phi(s,s,\theta) = \{\phi_0,\phi_0',\phi_1,\phi_1'\}$.

\subsection{An example where the supremum  \eqref{dualQOT1} is not achieved}\label{subsec:nonexisF}
Assume that $m=n=2$, $ C=C^Q$, $ \rho^A=\diag((1,0)^\top)$ and $ \rho^B=(1/2)\mathbb{I}_2$.
Recall that in such a case, $\Gamma^Q(\rho^A,\rho^B)=\{\rho^A\otimes\rho^B\}$ and
\begin{eqnarray*}
 \rho^A\otimes \rho^B=\left[\begin{array}{rrrr}\frac{1}{2}&0&0&0\\mathbf{0}&\frac{1}{2}&0&0\\mathbf{0}&0&0&0\\mathbf{0}&0&0&0\end{array}\right].
\end{eqnarray*}
We can easily see that the supremum in~\eqref{dualQOT1} is not attained in this case.  Let $F$ be of the form \eqref{sigmaFform}.  Suppose that there exists $\sigma^A,\sigma^B\in\mathrm{S}(\mathcal{H}_2)$ such that $F\ge 0$ and $\mathrm{T}_{C^Q}^Q(\rho^A,\rho^B)=\tr (\sigma^A\rho^A+\sigma^B\rho^B)$. As in the proof of Proposition \ref{dualQOT} we deduce that
$\tr F( \rho^A\otimes \rho^B)=0$.  Hence the $(1,1)$ and $(2,2)$ entries of $F$ are zero.  Since $F\ge 0$ it follows that the first and the second row and column of $F$ are zero.
Observe next that the $(2,3)$ and $(3,2)$ entries of $F$ are $-1/2$. Hence such $\sigma^A,\sigma^B$ do not exist.
\subsection{A lower bound on $\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)$}\label{subsec:lowbdT}
We first give some complementary optimality conditions for the minimum QOT problem and the maximum dual problem for positive definite diagonal density matrices.  Let $f(X)$ be defined as in \eqref{deff(X)}.  The following lemma will be extremely useful for proving closed forms for QOT for diagonal qubits and qutrits.
\begin{lemma}\label{compcondlem}  Assume that $\mathbf{s},\mathbf{t}\in\mathbb{R}^n$ are nonnegative probability vectors and $\rho^A = \diag(\mathbf{s})$, $\rho^B = \diag(\mathbf{t})$.  Then the dual supremum problem \eqref{dualQOT1} can be restricted to diagonal matrices $\sigma^A=-\diag(\mathbf{a}),\sigma^B=-\diag(\mathbf{b})$ for $\mathbf{a}, \mathbf{b}\in\mathbb{R}^n$ which satisfy the condition that $F=C^{Q}+\diag(\mathbf{a})\otimes \mathbb{I}_n +\mathbb{I}_n\otimes \diag(\mathbf{b})$ is positive semidefinite.

Let $X^\star=[x_{ij}^\star]\in \Gamma^{cl}(\mathbf{s},\mathbf{t})$ be a solution to the second minimum problem in \eqref{clqotdiagdm0}, where $\mathbf{p}^A=\mathbf{s},\mathbf{p}^B=\mathbf{t}$ and $m=n$.
Assume that the maximum in the dual supremum problem \eqref{dualQOT1} is achieved by a matrix
of the form $F^\star=C^{Q}+\diag(\mathbf{a}^\star)\otimes \mathbb{I}_n +\mathbb{I}_n\otimes \diag(\mathbf{b}^\star)$, where $\rho^A=\diag(\mathbf{s}),\rho^B=\diag(\mathbf{t}),\sigma^A=-\diag(\mathbf{a}),\sigma^B=-\diag(\mathbf{b})$.
Then the following equalities hold:
\begin{equation}\label{compcond}
\begin{aligned}
& x^\star_{ii} (a_i^\star+b_i^\star)=0, \, \textrm{ for }i\in[n],\\
& x_{ij}^\star(a_{i}^\star +b_{j}^\star +1/2) + x_{ji}^\star(a_{j}^\star +b_{i}^\star+1/2) -\sqrt{x_{ij}^\star x_{ji}^\star}=0, \, \textrm{ for }1\le i< j\le n.
\end{aligned}
\end{equation}
Furthermore the following conditions are satisfied
\begin{enumerate}[(a)]
\item For $i\ne j$ either $x^\star_{ij} x^\star_{ji}>0$ or $x^\star_{ij} =x^\star_{ji}=0$.
\item  Assume that $x^\star_{ii} x^\star_{jj}>0$.  Then $x^\star_{ij} =x^\star_{ji}$.
Let $X(t)$ be obtained from $X^\star$ by replacing the entries $x^\star_{ii}, x^\star_{ij},x^\star_{ji}, x^\star_{jj}$ with $x^\star_{ii}-t, x^\star_{ij}+t,x^\star_{ji}+t, x^\star_{jj}-t$.  Then $X(t)$ is also a solution to the second minimum problem in \eqref{clqotdiagdm0} for $t\in[-x_{ij}^\star, \min(x^\star_{ii},x^\star_{jj})]$.  Furthermore, $a_i^\star=a_j^\star=-b_i^\star=-b_j^\star$.
\item Suppose that $x_{ip}^\star,x_{iq}^\star,x_{jp}^\star,x_{jq}^\star$ are positive for $i\ne j, p\ne q$, where $i,j,p,q\in[n]$.  Then
\begin{equation}\label{compcond1}
\begin{aligned}
& \frac{\sqrt{x_{pi}^\star}}{\sqrt{x_{ip}^\star}}+\frac{\sqrt{x_{qj}^\star}}{\sqrt{x_{jq}^\star}}-\frac{\sqrt{x_{qi}^\star}}{\sqrt{x_{iq}^\star}}-\frac{\sqrt{x_{jp}^\star}}{\sqrt{x_{pj}^\star}}=0, && \textrm{ if } i\ne p, i\ne q, j\ne i, j\ne q,\\
& 1+\frac{\sqrt{x_{qj}^\star}}{\sqrt{x_{jq}^\star}}-\frac{\sqrt{x_{qi}^\star}}{\sqrt{x_{iq}^\star}}-\frac{\sqrt{x_{jp}^\star}}{\sqrt{x_{pj}^\star}}=0, && \textrm{ if } i=p,i\ne q, i\ne j, j\ne q.
\end{aligned}
\end{equation}
\end{enumerate}
Furthermore, there exists a minimizing matrix $X^\star$, satisfying the above conditions, such that it has at most one nonzero diagonal entry even if a maximizing $F^\star$ does not exist.
\end{lemma}
\begin{proof} Let $\mathbf{a}=(a_1,\ldots,a_n)^\top,\mathbf{b}=(b_1,\ldots,b_n)^\top\in\mathbb{R}^n$, and consider the matrix $F=C^{Q}+\diag(\mathbf{a})\otimes\mathbb{I}_n +\mathbb{I}_n\otimes \diag(\mathbf{b})$ .  Then
$F$ is a direct sum of $n$ blocks of size one of the form $a_i+b_i$ corresponding to the diagonal entries $((i,i),(i,i))$ and $n(n-1)/2$ blocks of size two corresponding to
the entries $((i,j)(i,j)), ((i,j)(j,i))$, $((j,i)(i,j)), ((j,i),(j,i))$:
\begin{equation}\label{defMij}
M_{ij}=\begin{bmatrix}a_i+b_j+1/2& -1/2\\-1/2& a_j+b_i+1/2\end{bmatrix},  \quad i\in[n]
\end{equation}
Hence $F\ge 0$ if and only if the following inequalities hold:
\begin{align}\label{posdeFcond}
& a_i+b_i\ge 0, \quad \textrm{ for } i\in[n],\\
& a_i+b_j+1/2\ge 0, \, a_j+b_i+1/2\ge 0,\, (a_i+b_j+1/2)(a_j+b_i+1/2)\ge 1/4, \quad i\ne j. \notag
\end{align}
Assume that $G=C^{Q}-\sigma^A\otimes \mathbb{I}_n - \mathbb{I}_n\otimes \sigma^B\ge 0$.  Let $\mathbf{a},\mathbf{b}\in\mathbb{R}^n$ be the vectors obtained from the diagonal entries of $-\sigma^A,-\sigma^B$ respectively.  Observe that the $n$ $1\times 1$  and $n(n-1)/2$   diagonal blocks of $F$ and $G$ discussed above are identical.  As $G$ is positive semidefinite
then $F$ is positive semidefinite.  Clearly
\begin{equation*}
\tr \sigma^A\diag(\mathbf{s})=-\tr\diag(\mathbf{a})\diag(\mathbf{s}), \quad \tr \sigma^B\diag(\mathbf{t})=-\tr\diag(\mathbf{b})\diag(\mathbf{t}).
\end{equation*}
Hence the dual supremum problem \eqref{dualQOT1} can be restricted to diagonal matrices $\sigma^A=-\diag(\mathbf{a}),\sigma^B=-\diag(\mathbf{b})$ for $\mathbf{a}, \mathbf{b}\in\mathbb{R}^n$ that satisfy the condition that $F$ is positive semidefinite.

Recall that $X^\star$ induces a solution to the original  SDP
$R^\star\in \Gamma^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))$ of the form
described in part (a) of Lemma \ref{diaglemobs}.   That is, the diagonal entries of $R^\star$
are $R^\star_{(i,j)(i,j)}=x^\star_{ij}$ with additional nonnegative entries:
$R^\star_{(i,j)(j,i)}=\sqrt{x^\star_{ij}x^\star_{ji}}$ for $i\ne j$.
Clearly, $R^\star $ is a direct sum of $n$ submatrices of order $1$ and $n(n-1)/2$ of order $2$ as above.
The implication \eqref{orthcond} yields that $\tr F^\star R^\star=0$.

 As $F^\star$ is positive semidefinite we deduce
the conditions \eqref{posdeFcond} for $\mathbf{a}^\star, \mathbf{b}^\star$.
The blocks  $[x_{ii}^\star]$ and $[a_i^\star+b_i^\star]$  contribute $1$ to the ranks of $R^\star$ and $F^\star$ if and only if $x_{ii}^\star>0$ and $a_i^\star+b_i^\star>0$.
Each $2\times 2$ block of $R^\star$ is of the form $\begin{bmatrix} x_{ij}^\star&\sqrt{x^\star_{ij}x^\star_{ji}}\\\sqrt{x^\star_{ij}x^\star_{ji}}& x_{ji}^\star\end{bmatrix}$ for $1\le i <j\le n$.   Note that the rank of this block is either zero or one.
Each corresponding $2\times 2$ submatrix  of $F^\star$ is of the form
$M_{ij}^\star$ given by \eqref{defMij}.
Thus $M_{ij}^\star$ is positive semidefinite with rank at least one.
This matrix has rank one if and only if  the following quadratic condition holds:
\begin{equation}\label{3quadcond}
(a_i^\star+ b_j^\star+1/2)(a_j^\star+b_i^\star+1/2)-1/4=0, \textrm{ for } 1\le i < j\le n.
\end{equation}

Recall the complementary condition
\begin{align*}
0 & =\tr R^\star F^\star \\
& =\sum_{i=1}^n x_{ii}^\star(a_{i}^\star+b_{i}^\star)+\!\!\sum_{1\le i <j\le n}\!\big(x_{ij}^\star(a_{i}^\star +b_{j}^\star +1/2) + x_{ji}^\star(a_{j}^\star +b_{i}^\star+1/2) -\sqrt{x_{ij}^\star x_{ji}^\star}\big).
\end{align*}
As all three $1\times 1$ and  $2\times 2 $ corresponding blocks of $R^\star$ and $F^\star$ are positive semidefinite, it follows that we have the complementary conditions \eqref{compcond}.

We now show the second part of the lemma.

\noindent
(a) Assume that $x_{ij}^\star=0$ for $i\ne j$.  Then the second part of \eqref{compcond} yields $x_{ji}^\star(a_{j}^\star +b_{i}^\star+1/2)=0$.  The second condition in \eqref{defMij} yield that $x^\star_{ji}=0$.

\noindent
(b) Observe that $X(t)\in \Gamma^{cl}(\mathbf{s},\mathbf{t})$ for $t\in[-\min(x_{ij}^\star,x_{ji}^\star), \min(x^\star_{ii},x^\star_{jj})]$.  Assume first that $x_{ij}^\star x_{ji}^\star>0$.
As $t=0$ is an interior point of this interval, and $X(0)=X^\star$ we have the critical condition $\left.\frac d{dt}f(X(t))\right|_{t = 0}$, with $f$ given by \eqref{deff(X)}.  This yields the equality
$2 -\frac{\sqrt{x_{ij}^\star}}{\sqrt{x_{ji}^\star}} -\frac{\sqrt{x_{ji}^\star}}{\sqrt{x_{ij}^\star}}=0$.  Hence $x_{ij}^\star=x_{ji}^\star$ and thus $f(X(t))=f(X(0))$ for
$t\in[-x_{ij}^\star, \min(x^\star_{ii},x^\star_{jj})]$.

Assume now that $x_{ij}^\star=x_{ji}^\star=0$.  Then $f(X(t))=f(X(0))$ for $t\in [0,\min(x^\star_{ii},x^\star_{jj})]$.

It is left to show that $a_i^{\star}=a_j^\star=-b_i^\star=-b_j^\star$.  First observe that
the first set of conditions of \eqref{compcond} yield that $a_i^\star+b_i^\star=a_j^\star+b_j^\star=0$.  By replacing $\mathbf{a}^\star, \mathbf{b}^\star$ by
$\mathbf{a}^\star-c\mathbf{1}, \mathbf{b}^\star+c\mathbf{1}$ we do not change $F^\star$.  Hence we can assume that $a_j^\star=b_j^\star=0$.  Set $b_i^\star=-a_i^\star$.  Then the assumption that the diagonal entries of $M_{ij}^\star$ are nonnegative yields that $|a_i^\star|\le 1/2$.  Use the assumption that $\det M_{ij}^\star\ge 0$ to deduce that $0=a_i^\star=-b_i^\star$.

\noindent
(c)  Let $X(t)$ be the matrix obtained from $X^\star$ by replacing  $x_{ip}^\star,x_{iq}^\star,x_{jp}^\star,x_{jq}^\star$ with  $x_{ip}^\star-t,x_{iq}^\star+t,x_{jp}^\star +t,x_{jq}^\star-t$. Then for $t\in [-\min(x_{iq}^\star,x_{jp}^\star), \min(x^\star_{ip},x^\star_{jq})]$ we have $X(t)\in\Gamma^{cl}(\mathbf{s},\mathbf{t})$.  As $t=0$ is an interior point of this interval we deduce that $\left.\frac d{dt}f(X(t))\right|_{t = 0}$.

Suppose first that $i\ne p, i\ne q, j\ne i, j\ne q$. Then Eq. \eqref{deff(X)} yields
\begin{multline*}
\!\!\!\!\! f(X(t))=-\Big(\sqrt{(x^\star_{ip}-t)x^\star_{pi}}+\sqrt{(x^\star_{iq}+t)x^\star_{qi}}+\sqrt{(x^\star_{jp}+t)x^\star_{pj}}+\sqrt{(x^\star_{jq}-t)x^\star_{qj}} \,\Big) \\
\qquad + C,
\end{multline*}
where $C$ is a term that does not depend on $t$.
The condition $\left.\frac d{dt}f(X(t))\right|_{t = 0}$ yields the first condition \eqref{compcond1}.

Assume now that $i=p$ and $i\ne q, j\ne i, j\ne q$.
Then we have
\begin{eqnarray*}
f(X(t))=t/2-\big(\sqrt{(x_{iq}^\star+t)x^\star_{qi}}+\sqrt{(x^\star_{jp}+t)x^\star_{pj}}+\sqrt{(x^\star_{jq}-t)x^\star_{qj}} \,\,\big)+ C,
\end{eqnarray*}
where $C$ does not depend on $t$.
Now, the condition $\left.\frac d{dt}f(X(t))\right|_{t = 0}$ yields the second condition in \eqref{compcond1}.

Finally, we need to prove the existence of an $X^\star$ with at most one nonzero entry that satisfies the conditions of the lemma.
Assume first that $\mathbf{s},\mathbf{t}>\mathbf{0}$.  Then Theorem~\ref{dualQOT} yields that there exists a maximizing matrix $F^\star$ to the dual supremum problem.  As we showed above we can assume that $F^\star=C^{Q}+\diag(\mathbf{a}^\star)\otimes \mathbb{I}_n +\mathbb{I}_n\otimes \diag(\mathbf{b}^\star)$.  Let $X^\star$ be a minimizing matrix with
at most $k$ zeros on the diagonal.
Assume to the contrary that  $x^\star_{ii} x^\star_{jj}>0$ for $1\le i < j\le n$.
Part (b) yields that for $t\in[-\min(x_{ij}^\star,x_{ji}^\star), \min(x^\star_{ii},x^\star_{jj})]$ the matrix  $X(t)$ minimizes $f$.  Choose $t^\star=\min(x^\star_{ii},x^\star_{jj})$.  Then $X(t^\star)$ is a minimizing matrix with at least $k+1$
zeros on the diagonal,
which contradicts our choice of $X^\star$.

Assume now that $\mathbf{s},\mathbf{t}$ are nonnegative. Let $\mathbf{s}_k,\mathbf{t}_k>\mathbf{0},k\in\mathbb{N}$ be two sequences that converge to $\mathbf{s},\mathbf{t}$ respectively.  Let $X_k^\star$ be a minimizing matrix of $f(X)$ corresponding to $\mathbf{s}_k,\mathbf{t}_k$ that has at most one nonzero diagonal element.  Clearly, there exists a subsequence $X^\star_{k_l}$ which has either all zero diagonal elements or exactly one positive diagonal element in a fixed diagonal entry.  Choose a subsequence $[\tilde x^\star_{ij,l}], l\in\mathbb{N}$ of this subsequence which converges to $X^\star$.  Clearly $X^\star$ is a minimizing matrix of $f(X)$ corresponding to $\mathbf{s},\mathbf{t}$.  If $x^\star_{ij}>0$ then $\tilde x^\star_{ij,l}>0$ for $l\gg 1$.  Hence $X^\star$ satisfies the conditions of the lemma.
\end{proof}
\begin{theorem}\label{lowbdTdiagn} Assume that $\mathbf{s}=(s_1,\ldots,s_n)^\top$, $\mathbf{t}=(t_1,\ldots,t_n)^\top\in\mathbb{R}^n_+$ are probability vectors and $U\in\mathrm{U}(n)$.  Then
\begin{equation}\label{eqlowbdiagn}
\mathrm{T}_{C^{Q}}^Q \big(U^\dagger \diag(\mathbf{s}) U, U^\dagger\diag(\mathbf{t}) U \big)\ge\frac{1}{2} \max_{i\in[n]} \, \big( \sqrt{s_i}-\sqrt{t_i} \big)^2
\end{equation}
Equality holds
if and only the exists $i\in[n]$ such that
\begin{equation}\label{eqlowbdiagn1}
\begin{aligned}
\textrm{ either } s_j\ge t_j \textrm{ and } t_it_j\ge s_i s_j \textrm{ for all } j\ne i, \\
\textrm{ or } t_j\ge s_j \textrm{ and } s_is_j\ge t_i t_j \textrm{ for all } j\ne i.
\end{aligned}
\end{equation}
\end{theorem}
\begin{proof}  Without loss of generality we can assume that  $U=\mathbb{I}_n$.  Suppose first that $\mathbf{s},\mathbf{t}>\mathbf{0}$.  Lemma \ref{compcondlem}  yields that $\mathrm{T}_{C^{Q}}^Q$ is the maximum of the dual problem where $F=C^{Q}+\diag(\mathbf{a})\otimes \mathbb{I}_n + \mathbb{I}_n\otimes \diag(\mathbf{b})$ is positive semidefinite.
Choose $i\in[n]$.  Assume that the coordinates of $\mathbf{a},\mathbf{b}$ are given as follows:
\begin{equation}\label{abichoice}
a_{i}=\frac{1}{2}\Big(\frac{\sqrt{t_i}}{\sqrt{s_i}}-1\Big), \, b_{i}=\frac{1}{2}\Big(\frac{\sqrt{s_i}}{\sqrt{t_i}}-1\Big), \quad a_{j}=b_{j}=0 \textrm{ for } j\ne i.
\end{equation}
Clearly
\begin{align*}
& a_{i}+b_{i}=\frac{(\sqrt{s_i}-\sqrt{t_i})^2}{2\sqrt{s_it_i}}\ge 0,\quad a_{j}+b_{j}=0, && \textrm{ for } j\ne i,\\
& 1/2+a_{i}>0, \,1/2+b_{i}>0, \quad 1/2+a_{j}=1/2+b_{j}=1/2, && \textrm{ for } j\ne i,\\
& (a_{i}+b_{j}+1/2)(a_{j}+b_{i}+1/2)=(a_{i}+1/2)(b_{i}+1/2)=1/4, && \textrm{ for } j\ne i,\\
& (a_{j}+b_{p}+1/2)(a_{p}+b_{j}+1/2)=1/2\times 1/2=1/4, && \textrm{ for }p\ne j\in [n]\setminus\{i\}.
\end{align*}
Thus $F\ge 0$.  Therefore
\begin{align*}
\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t})) &\ge  -\tr\big(\diag(\mathbf{a})\diag(\mathbf{s})+\diag(\mathbf{b})\diag(\mathbf{t})\big) \\
& = \frac{1}{2}\Big[ \Big(1-\frac{\sqrt{t_i}}{\sqrt{s_i}} \Big)s_i + \Big(1-\frac{\sqrt{s_i}}{\sqrt{t_i}}\Big) t_i\Big]=\frac{1}{2} \big(\sqrt{s_i}-\sqrt{t_i} \big)^2.
\end{align*}
As we let $i\in[n]$ we deduce the inequality \eqref{eqlowbdiagn}.
Since $\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))$ is continuous on $\Pi_n\times \Pi_n$ we deduce the inequality  \eqref{eqlowbdiagn} for all $(\mathbf{s},\mathbf{t})\in \Pi_n\times \Pi_n$.

 We now discuss the equality case in \eqref{eqlowbdiagn}.
Clearly $\max_{i\in[n]}(\sqrt{s_i}-\sqrt{t_i})^2=0$ if and only if $\mathbf{s}=\mathbf{t}$, in which case $\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=0$.  Assume that  $\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))> 0$.  Suppose first that equality holds in \eqref{eqlowbdiagn}. Then there exists an index $i\in[n]$ such that \mbox{$\mathrm{T}_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))=\frac{1}{2}(\sqrt{s_i}-\sqrt{t_i})^2>0$}.  By renaming indices and interchanging $\mathbf{s}$ and $\mathbf{t}$ if needed we can assume that $t_1>s_1$ and  $\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=\frac{1}{2}(\sqrt{s_1}-\sqrt{t_1})^2$.
Let $X=X^\star$ be a solution to the second minimum problem in \eqref{clqotdiagdm0}.
Recall that $f(X^\star)=\frac{1}{2}(\sqrt{t_1}-\sqrt{s_1})^2$.
Suppose first that $s_1=0$.  Then the first row of each $X\in \Gamma^{cl}(\mathbf{s},\mathbf{t})$ is zero.  Hence
\begin{eqnarray*}
2f(X)=\sum_{j=2}^n x_{j1}+\sum_{2\le j<k\le n}(\sqrt{x_{jk}}-\sqrt{x_{kj}})^2=t_1 +
\sum_{2\le j <k\le n}(\sqrt{x_{jk}}-\sqrt{x_{jk}})^2,
\end{eqnarray*}
for $ X\in\Gamma^{cl}(\mathbf{s},\mathbf{t})$.  As $f(X^\star)=t_1$ we deduce that the submatrix $Y=[x_{jk}^\star]_{j,k\ge 2}$ is a nonnegative symmetric matrix.  Thus for $j\ge 2$
\begin{eqnarray*}
s_j=\sum_{k=1}^n x_{jk}^\star=x_{j1}^\star +\sum_{k=2}^n x^\star_{jk}=x_{j1}^\star +\sum_{k=2}^n x^\star_{kj}=x_{j1}^\star+t_j.
\end{eqnarray*}
Therefore $s_j\ge t_j$ and $t_1t_j\ge 0=s_1s_j$ for $j\ge 2$.  Hence the conditions \eqref{eqlowbdiagn1} hold.

Assume now that $s_1>0$.
Let $F$ be defined as above for $i=1$.  Our assumption is that $F=F^\star$ is a solution to the maximum dual problem.   Lemma \ref{compcondlem} yields the equalities
\eqref{compcond}.  Hence $x_{11}^\star=0$.  Next consider the second part of the equalities \eqref{compcond} for $i=1$ and $j\ge 2$:
\begin{equation*}
\frac{\sqrt{t_1}}{\sqrt{s_1}}x_{1j}^\star=\frac{\sqrt{s_1}}{\sqrt{t_1}}x_{j1}^\star=c_j\ge 0 \; \textrm{ for } j\ge 2.
\end{equation*}
Observe next that
\begin{eqnarray*}
s_1=\sum_{j=2}^n x_{1j}^\star =\frac{\sqrt{s_1}}{\sqrt{t_1}}\sum_{j=2}^n c_j \quad \Rightarrow \quad \sum_{j=2}^n c_j =\sqrt{s_1t_1}.
\end{eqnarray*}
Therefore
\begin{eqnarray*}
\sum_{j=2}^n \big( x_{1j}^\star+x_{j1}^\star -2\sqrt{x_{ij}^\star x_{ji}^\star} \big) =s_1+t_1-2\sum_{j=2}^n c_j=s_1+t_1-2\sqrt{s_1 t_1}=(\sqrt{s_1}-\sqrt{t_1})^2.
\end{eqnarray*}
Hence
\begin{equation*}
2f(X^\star)=(\sqrt{s_1}-\sqrt{t_1})^2+ \sum_{2\le j<k\le n}(\sqrt{x_{jk}}-\sqrt{x_{kj}})^2=(\sqrt{s_1}-\sqrt{t_1})^2.
\end{equation*}
Therefore the submatrix $Y=[x_{jk}^\star]_{j,k\ge 2}$ is a nonnegative symmetric matrix.  Observe next that
\begin{align*}
& s_j=x^\star_{j1}+\sum_{k=2}^n x^\star_{jk}=\frac{\sqrt{t_1}}{\sqrt{s_1}} c_j+\sum_{k=2}^n x^\star_{jk},\\
& t_j\,=x^\star_{1j}+\sum_{k=2}^n x^\star_{kj}=\frac{\sqrt{s_1}}{\sqrt{t_1}} c_j+\sum_{k=2}^n x^\star_{kj}, \; \textrm{ for }j\ge 2.
\end{align*}
As $Y$ is symmetric we obtain that
\begin{eqnarray*}
s_j-t_j=\frac{(t_1-s_1)c_j}{\sqrt{s_1 t_1}}\ge 0 \quad \Rightarrow \quad  c_j=\frac{(s_j-t_j)\sqrt{s_1t_1}}{t_1-s_1}.
\end{eqnarray*}
As
\begin{eqnarray*}
s_j\ge x^\star_{j1}=\frac{\sqrt{t_1}}{\sqrt{s_1}}c_j=\frac{(s_j-t_j)t_1}{t_1-s_1}
\end{eqnarray*}
we deduce that $t_1t_j\ge s_1 s_j$.  Hence the conditions \eqref{eqlowbdiagn1} hold.

Assume now that the conditions \eqref{eqlowbdiagn1} hold.  To be specific we assume that $t_1\ge s_1$ and $s_j\ge t_j$ for $j\ge 2$.  If $s_j=t_j$ for $j\ge 2$ then $\mathbf{s}=\mathbf{t}$ and equality holds in \eqref{eqlowbdiagn}.  Hence we assume that $t_1>s_1$.
Define $X=[x_{ij}]$ as follows:
\begin{eqnarray*}
x_{11}=0, \, x_{1j}=\frac{s_1(s_j-t_j)}{t_1-s_1},\, x_{j1}=\frac{t_1(s_j-t_j)}{t_1-s_1},\, x_{jk}=\frac{t_1t_j-s_1s_j}{t_1-s_1}\delta_{jk} \; \textrm{ for } j,k\ge 2.
\end{eqnarray*}
Then $X\in \Gamma^{cl}(\mathbf{s},\mathbf{t})$.  Furthermore $2f(X)=s_1+t_1-2\sqrt{s_1t_1}=(\sqrt{s_1}-\sqrt{t_1})^2$.  Therefore  $2\mathrm{T}_{C^{Q}}^Q(\mathbf{s},\mathbf{t})\le (\sqrt{s_1}-\sqrt{t_1})^2$. On the other hand, inequality \eqref{eqlowbdiagn} yields that $2\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))\ge  (\sqrt{s_1}-\sqrt{t_1})^2$.  Consequently, we conclude that \mbox{$\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))= \frac{1}{2}(\sqrt{s_1}-\sqrt{t_1})^2$.}
\end{proof}
\begin{corollary}\label{lowbdT}
For $\rho^A,\rho^B\in\Omega_n$ let $D_0(\rho^A,\rho^B)$ be defined as in \eqref{defD0rhosigma}, where $f(x)=\sqrt{x}$ for $x\ge 0$.
Then
\begin{equation}\label{lowbdT0}
\mathrm{T}_{C^{Q}}^Q(\rho^A,\rho^B)\ge D_0^2(\rho^A,\rho^B).
\end{equation}
Furthermore for $n=2$ equality holds in \eqref{lowbdT0}.
\end{corollary}
\begin{proof}
Let $\rho^A,\rho^B\in\Omega_n$.  Recall the equality $\mathrm{T}_{C^{Q}}^Q( \rho^A, \rho^B)=\mathrm{T}_{C^{Q}}^Q(U^\dagger \rho^A U, U^\dagger\rho^B U)$ for $U\in U(n)$, and the
inequality \eqref{diagdecr1}.  Use the inequality \eqref{eqlowbdiagn} for $U=\mathbb{I}_n$  to deduce
\begin{align*}
\mathrm{T}_{C^{Q}}^Q( \rho^A, \rho^B) &= \mathrm{T}_{C^{Q}}^Q(U^\dagger \rho^A U, U^\dagger\rho^B U)\ge \mathrm{T}_{C^{Q}}^Q(\diag(U^\dagger \rho^A U), \diag(U^\dagger\rho^B U))\\
& \ge\frac{1}{2}\max_{i\in[n]}\bigg(\sqrt{(U^\dagger\rho^A U)_{ii}}-\sqrt{(U^\dagger\rho^B U)_{ii}}\bigg)^2.
\end{align*}
Take the maximum over $U\in \mathrm{U}(n)$ and use the proof of Proposition \ref{metricdenmat} to deduce \eqref{lowbdT0}. Theorem \ref{mtqub}  yields the equality in \eqref{lowbdT0} for $n=2$.
\end{proof}
\section{Diagonal qutrits}\label{subsec:diagqut}
In this section we provide a closed formula for $\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))$ for diagonal qutrits, $n=3$.
\begin{theorem}\label{diaggen3}  Let $\mathbf{s}=(s_1,s_2,s_3)^\top,\mathbf{t}=(t_1,t_2,t_3)^\top\in\mathbb{R}^3$ be  probability vectors.
Then the quantum optimal transport problem
for diagonal qutrits
is determined by the given formulas in the following cases:
\begin{enumerate}[(a)]
\item
\begin{equation*}
\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=\frac{1}{2}\max_{p\in [3]}(\sqrt{s_p}-\sqrt{t_p})^2
\end{equation*}
if and only if the conditions \eqref{eqlowbdiagn1} hold for $n=3$.
\item Suppose that
there exists $\{p,q,r\}=\{1,2,3\}$ such that
\begin{equation}\label{condminb}
\begin{aligned}
& t_r\ge s_p+s_q \textrm{ and}\\
& \textrm{either } s_p\ge t_p>0, s_q\ge t_q>0 \; \textrm{ or } \; t_p\ge s_p>0, t_q\ge s_q >0.
\end{aligned}
\end{equation}
Then
\begin{equation}\label{QOTbchar}
\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=\frac{1}{2}\Big((\sqrt{s_p}-\sqrt{t_p})^2 +(\sqrt{s_q}-\sqrt{t_q})^2\Big).
\end{equation}
\item  Suppose that
there exists $\{p,q,r\}=\{1,2,3\}$ such that
\begin{equation}\label{condmin1}
s_p>t_q>0, \quad t_p>s_q>0, \quad s_q+s_r\ge t_p,
\end{equation}
and
\begin{equation}\label{posa+b2}
\begin{aligned}
& 1 +\frac{\sqrt{t_q}}{\sqrt{s_q}}-\sqrt{\frac{s_p-t_q}{t_p-s_q}}\ge  0, \qquad
1+\frac{\sqrt{s_q}}{\sqrt{t_q}}-\sqrt{\frac{t_p-s_q}{s_p-t_q}}\ge 0,\\
& \left(1+\frac{\sqrt{t_q}}{\sqrt{s_q}}-\sqrt{\frac{s_p-t_q}{t_p-s_q}} \,\right)
\left(1+\frac{\sqrt{s_q}}{\sqrt{t_q}}-\sqrt{\frac{t_p-s_q}{s_p-t_q}} \, \right)\ge 1,
\\
& \max \Big(\frac{s_q}{t_q},\frac{t_q}{s_q}\Big)\ge \max \Big(\frac{s_p-t_q}{t_p-s_q},\frac{t_p-s_q}{s_p-t_q} \Big).
\end{aligned}
\end{equation}
Then
\begin{equation}\label{QOPTd3char1}
\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))= \frac{1}{2}\Big((\sqrt{s_q}-\sqrt{t_q})^2 +(\sqrt{s_p-t_q}-\sqrt{t_p-s_q})^2\Big).
\end{equation}
\item Assume that $\mathbf{s}=(s_1,s_2,0)^\top,\mathbf{t}=(t_1,t_2,t_3)^\top$ are probability vectors.  Then
\begin{equation}\label{QOTspqutritf}
\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))=
\begin{cases}
\frac{1}{2}\big((\sqrt{t_1}-\sqrt{t_2})^2+t_3\big), & \textrm{ if } s_1\ge t_2 \textrm{ and } s_2\ge t_1,\\
\frac{1}{2}\big((\sqrt{t_1}-\sqrt{s_1})^2+t_3\big), & \textrm{ if } s_1< t_2 ,\\
\frac{1}{2}\big((\sqrt{t_2}-\sqrt{s_2})^2+t_3 \big), & \textrm{ if } s_2< t_1.
\end{cases}
\end{equation}
If $\mathbf{s}=(s_1,s_2,s_3)^\top,\mathbf{t}=(t_1,t_2,0)^\top$, then formula \eqref{QOTspqutritf} holds after the swapping $s_i \leftrightarrow t_i$.
\end{enumerate}
\end{theorem}
\begin{proof} (a) This follows from Theorem \ref{lowbdTdiagn}.

\noindent
(b)   Suppose that the condition \eqref{condminb} holds.  By relabeling the coordinates and interchanging $\mathbf{s}$ and $\mathbf{t}$ if needed we can assume the conditions \eqref{condminb} are satisfied with $p=1,q=2, r=3$:
\begin{equation*}
s_1\ge t_1>0,\quad s_2\ge t_2>0, \quad t_3\ge s_1+s_2.
\end{equation*}
Hence
\begin{equation}\label{Bopt}
X^\star=\begin{bmatrix}0&0&s_1\\mathbf{0}&0&s_2\\t_1&t_2&t_3-(s_1+s_2)\end{bmatrix}
\in\Gamma^{cl}(\mathbf{s},\mathbf{t}).
\end{equation}
We claim that the conditions \eqref{condminb} yield that $X^\star$ is a minimizing matrix for

\noindent
$\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))$ as given in \eqref{clqotdiagdm0}.  To show that we use the complementary conditions in Lemma \ref{compcondlem}.
Let $R^\star\in \Gamma^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))$ be the matrix induced by $X^\star$ of the form
described in part (a) of Lemma \ref{diaglemobs}.   That is, the diagonal entries of $R^\star$
are $R^\star_{(i,j)(i,j)}=x^\star_{ij}$ with additional nonnegative entries:
$R^\star_{(i,j)(j,i)}=\sqrt{x^\star_{ij}x^\star_{ji}}$ for $i\ne j$.
Clearly, $R^\star $ is a direct sums of $3$ submatrices of order $1$ and $3$ of order $2$ as above.    Let  $F^\star$ be defined as in Lemma~\ref{compcondlem} with the following parameters:
\begin{equation}\label{ab3id1}
\begin{aligned}
& a_{1}^\star=\frac{1}{2}\Big(\frac{\sqrt{t_1}}{\sqrt{s_1}}-1\Big), \qquad b_{1}^\star=\frac{1}{2}\Big(\frac{\sqrt{s_1}}{\sqrt{t_1}}-1\Big),\\
& a_{2}^\star=\frac{1}{2}\Big(\frac{\sqrt{t_2}}{\sqrt{s_2}}-1\Big), \qquad b_{2}^\star=\frac{1}{2}\Big(\frac{\sqrt{s_2}}{\sqrt{t_2}}-1\Big),\\
& a_3^\star=b_3^\star=0.
\end{aligned}
\end{equation}
We claim that  the conditions \eqref{condminb} yield that $F^\star$ is positive semidefinite.
We verify that the three blocks of size one and the three  blocks of size two
of $F^\star$ are positive semidefinite.  The condition $a_i^\star+b_i^\star\ge 0$ for $i\in[3]$ is straightforward.  The conditions for $M_{12}^\star$ and $M_{13}^\star$ are straightforward.  We now show that $M_{12}^\star$ is positive semidefinite.
First note that as $s_1\ge t_1$ and $s_2\ge t_2$ we get that $b_1^\star\ge 0$ and $b_2^\star\ge 0$.  Clearly $a_1^\star>-1/2$ and $a_2^\star>-1/2$.  Hence the diagonal entries of  $M_{12}^\star$ are positive.  It is left to show that $\det M_{12}^\star\ge 0$.  Set $u=\sqrt{t_1}/{\sqrt{s_1}}\le 1$ and $v=\sqrt{s_2}/{\sqrt{t_2}}\ge 1$.
Then
\begin{align*}
2(a_1^\star+b_2^\star+1/2)=u+v-1, && 2(a_2^\star+b_1^\star +1/2)=1/u+1/v-1,
\end{align*}
\begin{align*}
4\det M_{12}^\star & =(u+v-1)(1/u+1/v-1)-1 \\
& =\big(1/(uv\big) \big)\big(u+v-1)(u+v-uv)-uv\big) \\
& =\big(1/(uv)\big) (u+v)(1-u)(v-1)\ge 0.
\end{align*}
We next observe that equalities \eqref{compcond} hold.  The first three equalities hold
as $x_{11}^\star=x_{22}^\star=(a_3^\star +b_3^\star)=0$.  The equality of $i=1,j=2$ holds as $x_{12}^\star=x_{21}^\star=0$.  The equalities for $i=1, j=3$ and $i=2,j=3$ follow from the following equalities:
\begin{align*}
x_{13}^\star(a_1^\star+b_3^\star +1/2)+x_{31}^\star(a_3^\star+b_1^\star +1/2)=
\tfrac{1}{2}\big(s_1\tfrac{\sqrt{t_1}}{\sqrt{s_1}}+t_1\tfrac{\sqrt{s_1}}{\sqrt{t_1}}\big)=\sqrt{s_1t_1}=\sqrt{x_{13}^\star x_{31}^\star},\\
x_{23}^\star(a_2^\star+b_3^\star +1/2)+x_{32}^\star(a_3^\star+b_2^\star +1/2)=
\tfrac{1}{2}\big(s_2\tfrac{\sqrt{t_2}}{\sqrt{s_2}}+t_2\tfrac{\sqrt{s_2}}{\sqrt{t_2}}\big)=\sqrt{s_2t_2}=\sqrt{x_{23}^\star x_{32}^\star}.
\end{align*}
Hence $\tr R^\star F^\star=0$ and $X^\star$ is a minimizing matrix.  Therefore \eqref{QOTbchar} holds for $p=1$, $q=2$.

\bigskip

\noindent
(c)  Suppose that the condition \eqref{condmin1} holds.  By relabeling the coordinates we can assume the conditions \eqref{condmin1} are satisfied with $p=1,q=2, r=3$:
\begin{equation*}
s_1>t_2,\quad t_1> s_2, \quad s_2+s_3-t_1\ge 0.
\end{equation*}
Hence
\begin{equation}\label{A2opt}
X^\star=\begin{bmatrix}0&t_2&s_1-t_2\\mathbf{s}_2&0&0\\t_1-s_2&0&s_2+s_3-t_1\end{bmatrix}\in\Gamma^{cl}(\mathbf{s},\mathbf{t}).
\end{equation}
We claim that the conditions \eqref{posa+b2} yield that $X^\star$ is a minimizing matrix for

\noindent
$\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))$ as given in \eqref{clqotdiagdm0}.  To show this we use the complementary conditions in Lemma \ref{compcondlem}.
Let $R^\star\in \Gamma^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))$ be the matrix induced by $X^\star$ of the form
described in part (a) of Lemma \ref{diaglemobs}.
Recall that $R^\star $ is a direct sum of $3$ submatrices of order $1$ and $3$ of order $2$ as above.
 Let  $F^\star$ correspond to
\begin{equation}\label{abstarceq}
\begin{aligned}
a^\star_1=\frac{1}{2}\Big(\frac{\sqrt{t_1-s_2}}{\sqrt{s_1-t_2}}-1\Big), && a_2^\star=
\frac{1}{2}\Big(\frac{\sqrt{t_2}}{\sqrt{s_2}}-\sqrt{\frac{s_1-t_2}{t_2-s_1}} \, \Big),&& a_3^\star=0,\\
b^\star_1=\frac{1}{2}\Big(\frac{\sqrt{s_1-t_2}}{\sqrt{t_1-s_2}}-1\Big), &&
b_2^\star=\frac{1}{2}\Big(\frac{\sqrt{s_2}}{\sqrt{t_2}}-\sqrt{\frac{t_1-s_2}{s_1-t_2}} \,\Big), && b_3^\star=0.
\end{aligned}
\end{equation}
We claim that  \eqref{posa+b2} yield that $F^\star$ is positive semidefinite.
We verify that the three blocks of size one and the three blocks of size two matrices of $F^\star$ are positive semidefinite.  The condition $a_1^\star+b_1^\star\ge 0$ is straightforward.
To show the condition $a_2^\star+b_2^\star\ge 0$ we argue as follows.
Let
\begin{equation*}
u=\frac{\sqrt{t_1}}{\sqrt{s_1}}, \quad v=\sqrt{\frac{s_1-t_2}{t_2-s_1}}.
\end{equation*}
Then $2(a_2^\star+b_2^\star)=u+1/u-(v+1/v)$.  The fourth condition of \eqref{posa+b2}
is $\max(u,1/u)\ge \max(v,1/v)$.  As $w+1/w$ increases on $[1,\infty)$ we deduce that
$a_2^\star+b_2^\star\ge 0$.  Clearly $a_3^\star+b_3^\star=0$.
We now show that the matrices \eqref{defMij} are positive semidefinite, where the last three inequalities follow from the first three inequalities of \eqref{posa+b2}:
\begin{align*}
& 2(a_1^\star+b_2^\star +1/2)=\frac{\sqrt{s_2}}{\sqrt{t_2}}>0, \qquad\qquad 2(a_2^\star+b_1^\star +1/2)=\frac{\sqrt{t_2}}{\sqrt{s_2}}>0,\\
& (a_1^\star+b_2^\star +1/2)(a_2^\star+b_1^\star +1/2)-1/4=0,\\
& 2(a_1^\star+b_3^\star +1/2)=\frac{\sqrt{t_1-s_2}}{\sqrt{s_1-t_2}}>0, \qquad
 2(a_3^\star+b_1^\star+1/2)=\frac{\sqrt{s_1-t_2}}{\sqrt{t_1-s_2}}>0,\\
& (a_1^\star+b_3^\star +1/2)(a_3^\star+b_1^\star +1/2)-1/4=0,\\
& 2(a_2^\star+b_3^\star +1/2)=\frac{\sqrt{s_2}}{\sqrt{t_2}}-\sqrt{\frac{t_1-s_2}{s_1-t_2}}+1\ge 0,\\
& 2(a_3^\star+b_2^\star +1/2)=\frac{\sqrt{t_2}}{\sqrt{s_2}}-\sqrt{\frac{s_1-t_2}{t_1-s_2}}+1\ge 0,\\
& (a_2^\star+b_3^\star +1/2)(a_3^\star+b_2^\star +1/2)-1/4\ge 0.
\end{align*}

Moreover, the conditions \eqref{compcond} hold: As $x_{11}^\star=x_{22}^\star= a_3^\star+b^\star_3=0$ the first three conditions of \eqref{compcond} hold.
As $x_{23}^\star=x_{32}^\star=0$ the second conditions of  \eqref{compcond} for $p=2,q=3$ trivially hold.  The other two conditions follow from the following equalities:
\begin{align*}
& x_{12}^\star(a_{1}^\star +b_{2}^\star +1/2) + x_{21}^\star(a_{2}^\star +b_{1}^\star+1/2) -\sqrt{x_{12}^\star x_{21}^\star} \\
& \hspace*{3.5cm }=
t_2\frac{\sqrt{s_2}}{2\sqrt{t_2}}+s_2\frac{\sqrt{t_2}}{2\sqrt{s_2}}-\sqrt{t_2 s_2}=0,\\
& x_{13}^\star(a_{1}^\star +b_{3}^\star +1/2) + x_{31}^\star(a_{3}^\star +b_{1}^\star+1/2) -\sqrt{x_{13}^\star x_{31}^\star}\\
& \hspace*{3.5cm }=
(s_1-t_2)\frac{\sqrt{t_1-s_2}}{2\sqrt{s_1-t_2}}+s_2\frac{\sqrt{t_2}}{2\sqrt{s_2}}-\sqrt{(s_1-t_2)(t_1- s_2)}=0.
\end{align*}
$\tr F^\star R^\star=0$.   Therefore
\begin{multline*}
\mathrm{T}^Q_{C^Q}\big(\diag(\mathbf{s}),\diag(\mathbf{t})\big)=\tr C^Q R^\star\\
= \frac{1}{2}\big(t_2+s_2+(s_1-t_2)+(t_1-s_2)\big)-\sqrt{t_2 s_2}-\sqrt{(s_1-t_2)(t_1-s_2)}.
\end{multline*}
This proves \eqref{QOPTd3char1}.

\bigskip

\noindent
(d)  Observe that the third row of every matrix in $\Gamma^{cl}(\mathbf{s},\mathbf{t})$ is a zero row.  Let $\mathbf{s}'=(s_1,s_2)^\top$.  Thus $\Gamma^{cl}(\mathbf{s}',\mathbf{t})$ is obtained from $\Gamma^{cl}(\mathbf{s},\mathbf{t})$ by deleting the third row in each matrix in $\Gamma^{cl}(\mathbf{s},\mathbf{t})$.    Proposition \ref{redprop} yields that
\begin{equation*}
\mathrm{T}_{C^Q}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=\mathrm{T}_{C^Q_{2,3}}^Q(\diag(\mathbf{s}'),\diag(\mathbf{t})).
\end{equation*}
(See Lemma \ref{clqotdiagdm} for the definition of $C_{2, 3}^Q$.)   We use now the minimum characterization of $\mathrm{T}_{C^Q_{2,3}}^Q(\diag(\mathbf{s}'),\diag(\mathbf{t}))$ given in \eqref{clqotdiagdm0}.  Assume that the minimum is achieved for $X^\star=[x_{il}^\star]\in\Gamma^{cl}(\mathbf{s}',\mathbf{t}), i\in[2],l\in[3]$.  We claim that either $x_{11}^\star=0$ or $x_{22}^\star=0$.

Let $Y=[x_{il}^\star], i,l\in[2]$.  Suppose first that $Y=0$.  Then $t_1=t_2=0$ and $t_3=1$.  So $\diag(\mathbf{t})$ is a rank-one matrix and $\tr \big( \diag(\mathbf{s})\diag(\mathbf{t}) \big)=0$.  The equality  \eqref{QOTrankone} yields that $\mathrm{T}_{C^Q}^Q \big(\diag(\mathbf{s}),\diag(\mathbf{t}) \big)=1$.  Clearly, $s_1\ge t_2=0, s_2\ge t_1=0$.  Hence \eqref{QOTspqutritf} holds.

Suppose second that $Y\ne 0$. Then $t_1+t_2$, the sum of the entries of $Y$, is positive.  Using continuity arguments it is enough to consider the case $t_1,t_2,t_3>0$.
Denote by $\Gamma'$ the set of all matrices $X=[x_{il}] \in \Gamma^{cl}(\mathbf{s}',\mathbf{t})$ such that $x_{i3}=x_{i3}^\star$ for $i=1,2$.
Clearly $\min_{A\in\Gamma'} f(A)=f(Y)$.  We now translate this minimum to the minimum problem we studied above.

Let $Z=\frac{1}{t_1+t_2} Y$.  The vectors corresponding to the  row sums and the column sums $Z$ are the probabilty vectors $\hat \mathbf{s} = (\hat s_1,\hat s_2)^\top
 and
$\hat \mathbf{t}=\frac{1}{t_1+t_2}(t_1,t_2)^\top$ respectively.  Consider the minimum problem
$\min_{W\in \Gamma^{cl}(\hat\mathbf{s},\hat\mathbf{t})} f(W)$.  The proof of Lemma \ref{QOTdiagform2} yields that this minimum is achieved at $W^\star$ which has at least one zero diagonal element.  Hence $Y$ has at least one zero diagonal element.

Assume first that $Y$ has two zero diagonal elements.  Then $X^\star=\begin{bmatrix}0&t_2&s_1-t_2\\t_1&0&s_2-t_1\end{bmatrix}$.
This corresponds to the first case of  \eqref{QOTspqutritf}.  It is left to show that $X^\star$ is a minimizing matrix.  Using the continuity argument we may assume that  $s_1>t_2, s_2>t_1$.  Let $B\in \mathbb{R}^{2\times 3}$ be a nonzero matrix such that $X^\star+cB\in \Gamma^{cl}(\mathbf{s}',\mathbf{t})$ for $c\in[0,\varepsilon]$ for some small positive $\varepsilon$.  Then $B=\begin{bmatrix}a&-b&-a+b\\-a&b&a-b\end{bmatrix}$, where $a,b\ge 0$ and $a^2+b^2>0$. It is clear that
$f(X^\star)<f(X+cB)$ for each $c\in (0,\varepsilon]$.  This proves the first case of \eqref{QOTspqutritf}.

Assume second that $x_{11}^\star=0$ and $x_{22}^\star>0$.  Observe that $x_{21}^\star=t_1>0$.
We claim that $x_{13}^\star=0$.  Indeed, suppose that it is not the case.  Let $B=\begin{bmatrix}0&1&-1\\mathbf{0}&-1&1\end{bmatrix}$.  Then $X^\star + cB\in\Gamma^{cl}(\mathbf{s}',\mathbf{t})$ for $c\in[0,\varepsilon]$ for some positive $\varepsilon$.
Clearly $f(X^\star+cB)<f(X^\star)$ for $c\in(0,\varepsilon]$.  Thus contradicts the minimality of $X^\star$.  Hence $x_{13}^\star=0$.  Therefore $X^\star=\begin{bmatrix}0&s_1&0\\t_1&t_2-s_1&t_3\end{bmatrix}$.  This corresponds to the second case of \eqref{QOTspqutritf}.

The third case is when $x_{11}^\star >0$ and $x_{22}^\star=0$. We show, as in the second case, that $x_{23}^\star=0$.  Then $X^\star=\begin{bmatrix}t_1-s_2&t_2&t_3\\mathbf{s}_2&0&0\end{bmatrix}$.  This corresponds to the third case of \eqref{QOTspqutritf}.

The case $\mathbf{s}=(s_1,s_2,s_3)^\top,\mathbf{t}=(t_1,t_2,0)^\top$ is completely analogous, hence the proof is complete.
\end{proof}

Basing on the numerical studies we conjecture that the cases (a)--(d) exhaust the parameter space $\Pi_3 \times \Pi_3$. Nevertheless, we include for completeness an analysis of the quantum optimal transport $\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))$ under the assumption that this is not the case. The employed techniques might prove useful when studying more general qutrit states or diagonal ququarts.

\begin{proposition}
Let $O\subset\Pi_3\times \Pi_3$ be the set of pairs $\mathbf{s},\mathbf{t}$, which do not meet neither of conditions (a)--(d) from Theorem  \ref{diaggen3}. Suppose that $O$ is nonempty. Then
each minimizing $X^\star$ in the characterization \eqref{clqotdiagdm0} of $\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))$ has zero diagonal.
Let $O'\subset O$ be an open dense subset of $O$ such that for each $(\mathbf{s},\mathbf{t})\in O'$ and each triple $\{i,j,k\}=[3]$ the  inequalities $s_p\ne t_q$ and $s_p+s_q\ne t_r$ hold.  Assume that $(\mathbf{s},\mathbf{t})\in O'$.
The set of matrices in $\Gamma^{cl}(\mathbf{s},\mathbf{t})$ with zero diagonal is an interval spanned by two distinct extreme points $E_1,E_2$, which have exactly five positive off-diagonal elements.
Let $Z(u)=uE_1+(1-u)E_2$ for $u\in[0,1]$.  Then the minimum of the function $f(Z(u)),  u\in[0,1]$, where $f$ is defined by \eqref{deff(X)},  is attained at a unique point $u^\star\in(0,1)$.  The point $u^\star$ is the unique solution in the interval $(0,1)$ to a polynomial equation of degree at most $12$.  The matrix $X^\star=Z(u^\star)$ is the minimizing matrix for the second minimum problem in \eqref{clqotdiagdm0}, and $\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=f(X^\star)$.
\end{proposition}
\begin{proof}
Assume first that the set $O\subset \Pi_3\times \Pi_3$ is  nonempty and satisfies the conditions \emph{(i)-(iv)}.
Combine Theorem \ref{lowbdTdiagn} with part (a) of the theorem to deduce that if
the conditions \eqref{eqlowbdiagn1} do hold for $n=3$ then
\begin{equation}\label{negcond(a)}
\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))>\max_{p\in[3]} \frac{1}{2}(\sqrt{s_p}-\sqrt{t_p})^2.
\end{equation}
In view of our assumption the above inequality holds.
We first observe that $s_p\ne t_p$ for each $p\in[3]$.  Assume to the contrary that $s_p=t_p$.  Without loss of generality we can assume that $s_3=t_3$.  Assume that in addition $s_q=t_q$ for some $q\in[2]$.  Then $\mathbf{s}=\mathbf{t}$ and
\begin{equation*}
\mathrm{T}_{C^{Q}}^Q(\diag(\mathbf{s}),\diag(\mathbf{t}))=\frac{1}{2}\max_{p\in [3]}(\sqrt{s_p}-\sqrt{t_p})^2=0
\end{equation*}
This contradicts \eqref{negcond(a)}. Hence there exists $q\in[2]$ such that $s_q>t_q$ for $q\in[2]$.  Without loss of generality we can assume that $s_2>t_2$, therefore $s_1<t_1$, as $s_1+s_2=t_1+t_2=1-s_3=1-t_3$.   Hence for $Y=\begin{bmatrix}s_1&0\\t_1-s_1&t_2\end{bmatrix}$ we have
$X=Y\oplus [s_3]\in \Gamma^{cl}(\mathbf{s},\mathbf{t})$.  Recall that $\mathbf{s},\mathbf{t}>\mathbf{0}$.
We replace $Y$ by $Y^\star=Y+u^\star\begin{bmatrix}-1&1\\mathbf{1}&-1\end{bmatrix}$ such that $u^\star>0, Y^\star\ge 0$ and one of the diagonal elements of $Y^\star$ is zero.
By relabeling $\{1,2\}$ if necessary we can assume that $Y^\star=\begin{bmatrix}0& s_1\\t_1&t_2-s_1\end{bmatrix}$ So $t_2\ge s_1$ and $X^\star=Y^\star\oplus [s_3]\in\Gamma^{cl}(\mathbf{s},\mathbf{t})$.  The minimal characterization  \eqref {clqotdiagdm0} of $\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))$ yields
\begin{eqnarray*}
\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))\le f(X^\star)=\frac{1}{2}(\sqrt{s_1}-\sqrt{t_1})^2.
\end{eqnarray*}
This contradicts \eqref{negcond(a)}.

As $\mathbf{s},\mathbf{t}>\mathbf{0}$ there exists a maximizing matrix $F^\star$ to the dual problem of the form given by Lemma  \ref{compcondlem}.
Let $X^\star$ be the corresponding minimizing matrix.
We claim that $X^\star$ has zero diagonal.  Assume first that $X^\star$ has a positive diagonal.  Then the arguments in part (b) of Lemma  \ref{compcondlem} yield that $X^\star$ is a symmetric matrix. Thus $\mathbf{s}=\mathbf{t}$, and this contradicts~\eqref{negcond(a)}.

Assume second that $X^\star$ has two positive diagonal entries.  By renaming the indices we can assume that $x_{11}^\star =0$, $x_{22}^\star, x_{33}^\star>0$.
Part (b) of Lemma  \ref{compcondlem} and the arguments of its proof yield that we can assume that $a_2^\star=a_3^\star=b_2^\star=0$.  Let $u^\star=a_1^\star+1/2, v^\star=b_1^\star+1/2$.  As $M_{12}^\star$ is positive semidefinite we have the inequalities: $u^\star\ge 0, v^\star\ge 0, u^\star v^\star\ge 1/4$.  Hence $x^\star>0, y^\star>0$.   Recall that $F^\star$ is a maximizing matrix for the dual problem \eqref{dualQOT1}.  Hence
\begin{align*}
\mathrm{T}^Q_{C^Q}\big(\diag(\mathbf{s}),\diag(\mathbf{t})\big) & = -(u^\star -1/2)s_1-(v^\star-1/2)t_1 \\
& = -u^\star s_1-v^\star t_1+(s_1+t_1)/2 \\
& \le -u^\star s_1-t_1/(4u^\star) +(s_1+t_1)/2\\
& \leq -\sqrt{s_1t_1} +(s_1+t_1)/2=(\sqrt{s_1}-\sqrt{t_1})^2/2.
\end{align*}
This contradicts \eqref{negcond(a)}.

We now assume that $X^\star$ has one positive diagonal entry.  Be renaming the indices $1,2,3$ we can assume that $x_{11}^\star=x_{22}^\star=0, x_{33}^\star>0$.
The conditions \eqref{compcond} yield that $a_3^\star+b_3^\star=0$.  Since we can choose $b_3^\star=0$ we assume that $a_3^\star=b_3^\star=0$.

Let us assume, case (A1), that $X^\star$ has six positive off-diagonal entries.  We first claim that either $x^\star_{13}=x^\star_{31}$ or $x^\star_{23}=x^\star_{32}$.  (Those are equivalent conditions if we interchange the indices $1$ and $2$.)  We deduce these conditions and an extra condition using the second conditions of \eqref{compcond1}.  First we consider $x^\star_{12}, x_{13}^\star, x_{32}^\star,x_{33}^\star$, that is $i=p=3$, $j=1, q=2$.  By replacing these entries by $x^\star_{12}-v, x_{13}^\star+v, x_{32}^\star+v,x_{33}^\star-v$
we obtain the equalities
\begin{align*}
1 +x=y+z, \qquad x=\frac{\sqrt{x_{21}^\star}}{\sqrt{x_{12}^\star}}, \quad y=\frac{\sqrt{x_{31}^\star}}{\sqrt{x_{13}^\star}}, \quad z=\frac{\sqrt{x_{23}^\star}}{\sqrt{x_{32}^\star}}.
\end{align*}
Second we consider $x^\star_{21}, x_{23}^\star, x_{31}^\star,x_{33}^\star$.
By replacing these entries by $x^\star_{21}-v, x_{23}^\star+v, x_{31}^\star+v,x_{33}^\star-v$ we obtain the equality:
\begin{equation*}
1+\frac{1}{x}=\frac{1}{z} +\frac{1}{y}.
\end{equation*}
Multiply the first and the second equality
to deduce
\begin{eqnarray*}
x+\frac{1}{x}=u+\frac{1}{u} , \quad u=\frac{y}{z}\Rightarrow \textrm{ either } x=u \textrm{ or } x=\frac{1}{u}.
\end{eqnarray*}
Assume first that $x=u=\frac{y}{z}$.  Substitute that into the first equality to deduce that $z=1$, which implies that $x_{23}^\star=x_{32}^\star$.  Similarly, if $x=1/u$ we deduce that $y=1$, which implies that $x_{13}^\star=x_{31}^\star$.  Let us assume for simplicity of exposition that $x_{23}^\star=x_{32}^\star$.  Let $X(w)$ be obtained from $X^\star$ by replacing $x_{22}^\star=0,x_{23}^\star, x_{32}^\star, x_{33}^\star$ with $x_{22}^\star+w,x_{23}^\star-w, x_{32}^\star-w, x_{33}^\star+w$ for $0<w<x_{23}^\star$.  Then $X(w)$ is a minimizing matrix and has two positive diagonal entries.  This contradicts our assumption

 that $X^\star$ has only one positive diagonal entry.

We now consider the case (A2) that $x^\star_{ij}=0$ for some $i\ne j$. Part (a) of Lemma \ref{compcondlem}  yields that $x^\star_{ji}=0$.
We claim that all four off-diagonal entries are positive.  Assume to the the contrary that $x^\star_{pq}=0$ for some $p\ne q$ and $\{p,q\}\ne \{i,j\}$.  Then $x_{qp}^\star=0$.
As $\mathbf{s},\mathbf{t}>\mathbf{0}$ we must have that $x^\star_{12}x^\star_{21}>0$ and all four other off-diagonal entries are zero.  But then $s_1=t_2, t_1=s_2, s_3=t_3$.   This is impossible since we showed that $s_3\ne t_3$.
Hence $X^\star$
has exactly four positive off-diagonal entries.

Let us assume first that  $x_{12}^\star=x_{21}^\star=0$.  Then $X^\star$ is of the form
given by \eqref{Bopt}, where $t_3>s_1+s_2$.
We now recall again the conditions \eqref{compcond}.  As we already showed, we can assume that $a_3^\star=b_3^\star=0$.  As $x_{11}^\star=x_{22}^\star=0$ all of the first three conditions of  \eqref{compcond} hold.  As $x_{12}^\star=x_{21}^\star=0$ the second condition of  \eqref{compcond} holds trivially for $i=1,j=2$. The conditions for $i=1,j=3$ and $i=2, j=3$ are
\begin{align*}
&s_1(a_1^\star+1/2)+t_1(b_1^\star+1/2)=\sqrt{s_1t_1},\\
&s_2(a_2^\star+1/2)+t_2(b_2^\star+1/2)=\sqrt{s_2t_2}.
\end{align*}
We claim that \eqref{ab3id1} holds.
Using the assumption that $\det M_{13}^\star\ge 1/4$ and the inequality of arithmetic and geometric means we deduce that $\det M_{13}^\star= 1/4$.  Hence
\begin{align*}
& a_1^\star+1/2=u, \quad b_1^\star+1/2=1/(4u), \qquad \textrm{ for some }u>0, \\
& s_1u+t_1/(4u)t_1\ge \sqrt{s_1t_1}.
\end{align*}
Equality holds if and only if $u=\sqrt{t_1}/(2\sqrt{s_1})$.  This shows the first equality in
\eqref{ab3id1}.   The second equality in \eqref{ab3id1} is deduced similarly.  We now show that the conditions \eqref{condminb} hold for $i=1,j=2,k=3$.  As $t_3>s_1+s_2$ the first condition of \eqref{ab3id1} holds.  We use the conditions that $M_{12}^\star $ is positive semidefinite.  Let $u=\sqrt{t_1}/{\sqrt{s_1}}, v=\sqrt{s_2}{\sqrt{t_2}}$.
Then the arguments of the proof of part (b) yield
\begin{align*}
& 2(a_1^\star +b_2^\star+1)=u+v-1 >0, \quad 2(a_2^\star +b_1^\star+1)=(1/u+1/v-1)>0,\\
& 4\det M_{12}^\star =\big(1/(uv)\big)(1-u)(v-1).
\end{align*}
So either $u\ge 1$ and $v\le 1$, or $u\le 1$ and $v\ge 1$.  Hence \eqref{condminb} holds for $i=1,j=2,k=3$.  This contradicts our assumption that \eqref{condminb} does not hold.

Let us assume second that  $x_{12}^\star>0, x_{21}^\star>0$.  Then either $x_{13}^\star=x_{31}^\star=0$ or $x_{23}^\star=x_{32}^\star=0$. By relabeling $1,2$ we can assume that $x_{23}^\star=x_{32}^\star=0$.  Hence $X^\star$ is of the form
\eqref{A2opt}, where $s_1>t_2>0, t_1>s_2>0, s_2+s_3>t_1$.  Hence the conditions
\eqref{condmin1} are satisfied with $i=1,j=2, k=3$.  We now obtain a contradiction by showing that the conditions \eqref{posa+b2} are satisfied.  This is done using the same arguments as in the previous case as follows.
First observe that the second  nontrivial conditions of  \eqref{compcond} are:
\begin{align*}
& t_2(a_1^\star+b_2^\star +1/2)+s_2(a_2^\star +b_1^\star+1/2)=\sqrt{s_2t_2},\\
& (s_1-t_2)(a_1^\star+1/2)+t_1(b_1^\star+1/2)=\sqrt{(s_1-t_2)(t_1-s_2)}.
\end{align*}
As in the previous case we deduce that
\begin{align*}
& a_1^\star+b_2^\star +1/2=\sqrt{s_2}/(2\sqrt{t_2}), && b_1^\star+a_2^\star +1/2=\sqrt{t_2}/(2\sqrt{s_2}),\\
& a_1^\star +1/2=\sqrt{t_1-s_2}/(2\sqrt{s_1-t_2}), && b_1^\star +1/2=\sqrt{s_1-t_2}/(2\sqrt{t_1-s_2}).
\end{align*}
 Hence \eqref{abstarceq} holds.  We now recall the proof of part (c) of the theorem.

 We have thus shown that the minimizing matrix $X^\star$ has zero diagonal

We now show that $O$ is an open set.  Clearly, the set of all pairs  of probability vectors  $O_1\subset \Pi_3\times \Pi_3$ such that at least one of them has a zero coordinate is a closed set. Let  $O_2, O_3, O_4\subset \Pi_3\times \Pi_3$ be the sets which satisfy the conditions (a), (b),(c) of the theorem respectively.   It it straightforward to show: $O_2$ is a closed set, and  Closure$(O_3)\subset (O_3\cup O_1)$.  We now show that Closure$(O_4)\subset O_4\cup O_1\cup O_2$.   Indeed, assume that we have a sequence $(\mathbf{s}_l,\mathbf{t}_l)\in O_4, l\in\mathbb{N}$ that converges to $(\mathbf{s},\mathbf{t})$.  It is enough to consider the case where $\mathbf{s},\mathbf{t}>\mathbf{0}$.  Again we can assume for simplicity that each  $(\mathbf{s}_l,\mathbf{t}_l)$ satisfies the conditions
\eqref{condmin1} and \eqref{posa+b2} for $i=1, j=2,k=3$.    Then we deduce that the limit of the minimizing matrices $X^\star_l$ is of the form \eqref{A2opt}.  Hence $\lim_{l\to\infty}X^\star_l=X^\star$, where $X^\star$ is of the form \eqref{A2opt}.
Also $X^\star$ is a minimizing matrix for $\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))$.  Recall that $s_2,t_2>0$.  If $s_1-t_2>0,t_1-s_2>0$ then $(\mathbf{s},\mathbf{t})\in O_4$.  So assume that $(s_1-t_2)(t_1-s_2)=0$.  As $X^\star$ is minimizes $\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))$ and $\mathbf{s},\mathbf{t}>\mathbf{0}$, part (a) of Lemma \ref{compcondlem} yields that $s_1=t_2, t_1=s_2$.  Hence $s_3=t_3$. As $X^\star$ is minimizes $\mathrm{T}^Q_{C^Q}(\diag(\mathbf{s}),\diag(\mathbf{t}))$ we get that $\mathrm{T}^Q_{C^Q}=\frac{1}{2}(\sqrt{s_2}-\sqrt{t_2})^2$.  Hence $(\mathbf{s},\mathbf{t})\in O_2$.  This shows that $O_1\cup O_2\cup O_3\cup O_4$ is a closed set.
Therefore $O=\Pi_3\times \Pi_3\setminus(O_1\cup O_2\cup O_3\cup O_4)$ is an  open set.  If $O$ is an empty set then proof of the theorem is concluded.

Assume that $O$ is a nonempty set.
Let $O'\subset O$ be an open dense subset of $O$ such that for each $(\mathbf{s},\mathbf{t})\in O'$ and each triple $\{p,q,r\}=[3]$ the  inequality $s_p\ne t_q$
and $s_p+s_q\ne t_r$ hold.

Assume that $(\mathbf{s},\mathbf{t})\in O'$.  Let $\Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ be the convex subset of $\Gamma^{cl}(\mathbf{s},\mathbf{t})$  of  matrices with zero diagonal.
We claim that any $X\in\Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ has at least $5$ nonzero entries.
Indeed, suppose that $X\in\Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ has two zero
off-diagonal entries.
As $\mathbf{s},\mathbf{t}>0$ they cannot be in the same row or column. By relabeling the rows we can assume that the two zero elements are in the first and the second row.
Suppose first that $x_{12}^\star=x_{23}^\star=0$.  Then $X=\begin{bmatrix}0&0&s_1\\
s_2&0&0\\t_1-s_2&t_2&0\end{bmatrix}$.  Thus $s_1=t_3$ which is impossible.  Assume now that $x_{12}^\star=x_{21}^\star=0$.  Then $s_1+s_2=t_3$ which is impossible.  All other choices also are impossible.

We claim that $ \Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ is spanned by two distinct extreme points $E_1,E_2$, which have exactly five positive off-diagonal elements.  Suppose first that there exists $X\in\Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ which has six postive off-diagonal elements.
Let
\begin{equation*}
B=\begin{bmatrix}0&1&-1\\-1&0&1\\mathbf{1}&-1&0\end{bmatrix}.
\end{equation*}
Then all matrices in $\Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ are of the form  $X^\star+uB, u\in [u_1,u_2]$ for some $u_1<u_2$.  Consider the matrix $E_1=X^\star +u_1B$.  It has at least one zero off-diagonal entry
 hence we conclude
that $E_1$ has exactly five off-diagonal positive elements.  Similarly $E_2=X+u_2B$ has five positive off-diagonal elements.  Assume now that $E\in\Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ has
five positive off-diagonal elements.
Hence there exits a small $u>0$ such that either $E+uB$ or $E-uB$ has six positive off-diagonal elements.  Hence $ \Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ contains a matrix with six positive diagonal elements.  Therefore $ \Gamma^{cl}_0(\mathbf{s},\mathbf{t})$ is an interval spanned by $E_1\ne E_2\in  \Gamma^{cl}_0(\mathbf{s},\mathbf{t})$, where $E_1$ and $E_2$ have five positive off-diagonal elements.  Part (a) of Lemma \ref{compcondlem} yields that $X^\star$ has six positive off-diagonal elements.  Consider $E_1$ and assume that the $(1,2)$ entry of $E_1$ is zero.  Then
\begin{equation*}
E_1=\begin{bmatrix}0&0&s_1\\mathbf{s}_1+s_2-t_3&0&t_3-s_1\\mathbf{s}_3-t_2&t_2&0 \end{bmatrix}.
\end{equation*}

As $f(E_1+uB)$ is strictly convex on $[0,u_3]$, there exists a unique $u^\star\in (0, u_3)$ which satisfies the equation
\begin{multline*}
-\frac{\sqrt{s_1+s_2-t_3-u}}{\sqrt{u}}+\frac{\sqrt{u}}{\sqrt{s_1+s_2-t_3-u}} -\frac{\sqrt{s_1-u}}{\sqrt{s_3-t_2+u}} \\
 +\frac{\sqrt{s_3-t_2+u}}{\sqrt{s_1-u}}-
\frac{\sqrt{t_2-u}}{\sqrt{t_3-s_1+u}}+\frac{\sqrt{t_3-s_1+u}}{\sqrt{t_2-u}}=0.
\end{multline*}
It is not difficult to show that the above equation is equivalent to a polynomial equation of degree at most $12$ in $u$.  Indeed, group the six terms into three groups, multiply by the common denominator, and pass the last group to the other side of the equality to obtain the equality:
\begin{align*}
&\sqrt{(s_1-u)(s_3-t_2+u)(t_3-s_1+u)(t_2-u)}(2u+t_3-s_1-s_2)\\
& \, +\sqrt{u(s_1+s_2-t_3-u)(t_3-s_1+u)(t_2-u)}(2u+t_3-s_1-s_2)(2u+s_3-s_1-t_2)\\
& \hspace*{1.85cm} =\sqrt{u(s_1+s_2-t_3-u)(s_3-t_2+u)(t_3-s_1+u)}(-2u+s_1+t_2-t_3).
\end{align*}
Raise this equality to the second power.  Put all polynomial terms of degree $6$ on the left hand side, and the one term with a square radical on the other side.
Raise to the second power to obtain a polynomial equation in $u$ of degree at most 12.   Hence $X^\star=E_1+u^\star B$.  This completes the proof of (e).
\end{proof}

\section{Quantum optimal transport for $d$-partite systems}\label{sec:QOTdpar}
We now explain briefly how to state the quantum optimal transport problem for $d$-partite system, where $d\ge 3$, similarly to what was done in \cite{Fri20}.  Let $\mathcal{H}_{n_j}$ be a Hilbert space of dimension $n_j$ for $j\in[d]$.
We consider the $d$-partite tensor product space $\otimes_{j=1}^d \mathcal{H}_{n_j}$.
A product state in Dirac's notation is $\otimes_{i=1}^d |\mathbf{x}_i\rangle$.  Then
\begin{eqnarray*}
\langle \otimes_{i=1}^d \mathbf{x}_i,\otimes_{j=1}^d \mathbf{y}_j\rangle=(\otimes_{i=1}^d \langle \mathbf{x}_i|)
(\otimes_{j=1}^d |\mathbf{y}_j\rangle)=
\prod_{j=1}^d \langle \mathbf{x}_j| \mathbf{y}_j
\rangle.
\end{eqnarray*}
Consider the space $\mathrm{B}(\otimes_{j=1}^d\mathcal{H}_{n_j})$ of linear operators from $\otimes_{j=1}^d\mathcal{H}_{n_j}$ to itself.  A rank-one product operator is of the form
$(\otimes_{i=1}^d |\mathbf{x}_i\rangle)(\otimes_{j=1}^d\langle\mathbf{y}_j|)$ and
acts on a product state as follows:
\begin{eqnarray*}
(\otimes_{i=1}^d |\mathbf{x}_i\rangle)(\otimes_{j=1}^d\langle\mathbf{y}_j|)(\otimes_{k=1}^d |\mathbf{z}_k\rangle)=(\prod_{j=1}^d\langle \mathbf{y}_j| \mathbf{z}_j\rangle) (\otimes_{i=1}^d|\mathbf{x}_i\rangle).
\end{eqnarray*}
Given $\rho^{A_1,\ldots ,A_d}\in \mathrm{B}(\otimes_{j=1}^d\mathcal{H}_{n_j})$ one can define a $k$-partial trace on $k\in[d]$:
\begin{align*}
&\tr_k: \mathrm{B} \big(\otimes_{j=1}^d\mathcal{H}_{n_j} \big)\to \mathrm{B} \big(\otimes_{j\in[d]\setminus\{k\}}\mathcal{H}_{n_j} \big),\\
&\tr_k (\otimes_{i=1}^d |\mathbf{x}_i\rangle)(\otimes_{j=1}^d\langle\mathbf{y}_j|)=\langle \mathbf{y}_k|\mathbf{x}_k\rangle
(\otimes_{i\in[d]\setminus\{k\}} |\mathbf{x}_i\rangle)(\otimes_{j\in[d]\setminus\{k\}}^d\langle \mathbf{y}_j|).
\end{align*}
We will denote $\tr_k\rho^{A_1,\ldots,A_d}$ by $\rho^{A_1,\ldots,A_{k-1},A_{k+1},\ldots,A_d}$.  Let $\rho^{A_k}\in\mathrm{B}(\mathcal{H}_{n_k})$ be the operator obtained from $\rho^{A_1,\ldots,A_d}$ by tracing out all but the $k$-th component.
Thus we have the map
\begin{align*}
& \widetilde{\mathrm{Tr}}:  \mathrm{B} \big(\otimes_{j=1}^d\mathcal{H}_{n_j} \big)\to \oplus_{j=1}^d \mathrm{B} \big(\mathcal{H}_{n_j} \big),\\
& \widetilde{\mathrm{Tr}}(\rho^{A_1,\ldots,A_d})=(\rho^{A_1},\ldots,\rho^{A_d}).
\end{align*}
Let $N=\prod_{j=1}^d n_j$ and view the set of density matrices $\Omega_N$ as a subset of selfadjoint operators on $\mathcal{H}_N=\otimes_{j=1}^d \mathcal{H}_{n_j}$.
For $\rho^{A_i}\in\Omega_{n_i}, i\in[d]$ denote
\begin{eqnarray*}
\Gamma^{Q}(\rho^{A_1},\ldots,\rho^{A_d})=\{\rho^{A_1,\ldots,A_d}\in \Omega_N,
\widetilde{\mathrm{Tr}} (\rho^{A_1,\ldots,A_d})=(\rho^{A_1},\ldots,\rho^{A_d})\}.
\end{eqnarray*}

Assume that  $C$ is a selfadjoint operator on $\mathcal{H}_N$.  We define the quantum optimal transport as
\begin{equation}\label{defdQT}
\mathrm{T}_{C}^Q(\rho^{A_1},\ldots,\rho^{A_d})=\min_{\rho^{A_1,\ldots,A_d}\in\Gamma^Q(\rho^{A_1},\ldots,\rho^{A_d})}\tr C\rho^{A_1,\ldots,A_d}.
\end{equation}

We now give an analog of
a result in \cite{Fri20}.  Assume that $d=2\ell\ge 4$, and
$n_1=\cdots=n_d=n$.  Then $\mathcal{H}_n^{\otimes d}=\otimes^d\mathcal{H}_n$.
We want to give a semidistance between two ordered $\ell$-tuples of density matrices $( \rho^{A_1},\cdots,\rho^{A_\ell}),(\rho^{A_{\ell+1}},\cdots,\rho^{A_{2\ell}})\in\Omega_n^\ell$.  We view $\mathcal{H}_n^{\otimes (2\ell)}$ as bipartite states $\mathcal{H}_n^{\otimes \ell}\otimes \mathcal{H}_n^{\otimes \ell}$.  Let $S\in\mathrm{B}(\mathcal{H}_n^{\otimes (2\ell)})$ be the SWAP operator:
\begin{eqnarray*}
S(\otimes_{j=1}^{2\ell} |\mathbf{x}_j)\rangle=(\otimes_{j=1}^{\ell} |\mathbf{x}_{j+\ell}\rangle)\otimes (\otimes_{j=1}^\ell |\mathbf{x}_j\rangle).
\end{eqnarray*}
Denote by
$C^{Q}=\frac{1}{2}(\mathbb{I}-S)$.
Then $\mathrm{T}_{C^{Q}}^Q(\rho^{A_1},\ldots,\rho^{A_{2\ell}})\ge 0$. Equality holds if and only if
$(\rho^{A_1},\ldots,\rho^{A_\ell})=(\rho^{A_{1+\ell}},\ldots,\rho^{A_{2\ell}})$.
Also
\begin{equation*}
\mathrm{T}_{C^{Q}}^Q(\rho^{A_1},\ldots,\rho^{A_{2\ell}})=\mathrm{T}_{C^{Q}}^Q(\rho^{A_{1+\ell}},\ldots,\rho^{A_{2\ell}}, \rho^{A_1},\ldots,\rho^{A_{\ell}}).
\end{equation*}
Hence $\mathrm{T}_{C^{Q}}^Q(\rho^{A_1},\ldots,\rho^{A_{2\ell}})$ is a semi-metric on $\Omega_n^\ell$.  As in the case of $\ell=1$ we can show that $\sqrt{\mathrm{T}_{C^{Q}}^Q(\rho^{A_1},\ldots,\rho^{A_{2\ell}})}$ is a weak metric.
Denote by

\noindent
$W^Q_{C^Q}((\rho^{A_1},\ldots,\rho^{A_\ell}),(\rho^{A_{\ell+1}},\ldots,\rho^{A_{2\ell}}))$ the Wasserstein-2 metric on $\Omega_n^\ell$ induced by the weak metric $\sqrt{\mathrm{T}_{C^{Q}}^Q(\rho^{A_1},\ldots,\rho^{A_{2\ell}})}$.

Let $\Sigma_\ell$ be the group of bijections $\pi:[\ell]\to[\ell]$.
Then
\begin{eqnarray*}
\min_{\pi\in\Sigma_{\ell}}W_{C^{Q}}^Q \big( (\rho^{A_{\pi(1)}},\ldots,\rho^{A_{\pi(\ell)}}), \rho^{1+\ell},\ldots,\rho^{2\ell}) \big)
\end{eqnarray*}
gives a metric on unordered $\ell$-tuples of density matrices.
We call this metric the quantum Wasserstein-2 metric on the set of unordered $\ell$-tuples $\{\rho^{A_1},\ldots,\rho^{A_\ell}\}$.

On $\mathcal{H}_n^{\otimes d}$ we define for two integers $1\le p<q\le d$ the SWAP operator $S_{pq}\in \mathrm{B}(\mathcal{H}_{n})^{\otimes d}$, which swaps $\mathbf{x}_p$ with $\mathbf{x}_q$ in the tensor product $|\mathbf{x}_1\rangle\otimes\cdots\otimes|\mathbf{x}_d\rangle$.  Note that $S_{pq}$ is unitary and involutive.  Hence $S_{pq}$ is selfadjoint with eigenvalues $\pm 1$.
The common invariant subspace of $\mathcal{H}_n^{\otimes d}$ for all $S_{pq}$ is the
the subspace of symmetric tensors\
---``bosons''---,  denoted as $\mathrm{S}^d\mathcal{H}_n$.
 Let $C^{B}\in \mathrm{S}_+(\mathcal{H}_n^{\otimes d})$ be the projection on the orthogonal complement of $\mathrm{S}^d\mathcal{H}_n$.  Note that  $C^{B}=C^Q$ for $d=2$.   We now have a partial analog of Theorem \ref{kapTABprop}:
\begin{theorem}\label{kapTd}
 Let $\rho^{A_1},\ldots,\rho^{A_d}\in\Omega_n$.  Then
 \begin{enumerate}[(a)]
 \item
 $\mathrm{T}_{C^{B}}^Q(\rho^{A_1},\ldots,\rho^{A_d})\ge 0$.
 \item
 $\mathrm{T}_{C^{B}}^Q(\rho^{A_1},\ldots,\rho^{A_d})=0$ if and only if $\rho^{A_1}=\cdots=\rho^{A_d}$.
 \item
 Assume that at least $d-1$ out of $\rho^{A_1},\ldots,\rho^{A_d}$ are pure states.
 Then
 \begin{equation*}
 \mathrm{T}_{C^{B}}^Q(\rho^{A_1},\ldots,\rho^{A_d})=\tr C^{B}(\otimes_{j=1}^d\rho^{A_j}).
 \end{equation*}
 \end{enumerate}
 \end{theorem}
 \begin{proof} (a) This follows from the fact that $\tr C^{B}\rho^{A_1,\ldots,A_d}\ge 0$

 \noindent
 (b) Assume that $\mathrm{T}_{C^{B}}^Q(\rho^{A_1},\ldots,\rho^{A_d})=\tr C^{B}\rho^{A_1,\ldots,A_d}=0$.  Hence all the eigenvectors of $\rho^{A_1,\ldots,A_d}$ corresponding to positive eigenvalues are symmetric tensors.  So $S_{pq} \rho^{A_1,\ldots,A_d} S_{pq}=\rho^{A_1,\ldots,A_d}$.  Therefore $\widetilde{\mathrm{Tr}}(\rho^{A_1,\ldots,A_d})=(\rho,\ldots,\rho)$.
 Thus $ \rho^{A_1}=\cdots=\rho^{A_d}=\rho$.  We now show that $\mathrm{T}_{C^{B}}^Q(\rho,\ldots,\rho)=0$.
 Suppose that $\rho$ has the spectral decomposition \eqref{specdecrho}.  Let us take
 a $d$-purification of $\rho$
 \begin{equation*}
 \rho^{\mathrm{pur},d}=\big(\sum_{i=1}^n \sqrt{\lambda_i}\otimes^d|\mathbf{x}_i\rangle \big) \big( \sum_{j=1}^n \sqrt{\lambda_j}\otimes^d\langle\mathbf{x}_j| \big).
 \end{equation*}
 Clearly we have $\rho^{\mathrm{pur},d}\in\Gamma^Q(\rho,\ldots,\rho)$.  As $\rho^{\mathrm{pur},d}$ is a pure state whose eigenvector corresponding to its positive eigenvalue is
 a symmetric tensor  we deduce that $\tr C^{B}\rho^{\mathrm{pur},d}=0$.

 \noindent (c) Assume for simplicity of the exposition that $\rho^{A_2},\ldots,\rho^{A_d}$ are pure states.  Then $\rho^B=\otimes_{j=2}^d \rho^{A_j}$ is a pure state.    Lemma \ref{rangecont} yields that $\Gamma^Q(\rho^{A_1},\rho^{B})=\{\rho^{A_1}\otimes \rho^{B}\}$.  Hence $\Gamma^{Q}(\rho^{A_1},\ldots,\rho^{A_d})=\{\otimes_{j=1}^d \rho^{A_j}\}$. This proves part (c) of the theorem.
 \end{proof}

The next question concerns the optimal technique
 to compute $\tr C^{B}(\otimes_{j=1}^d\rho^{A_j})$.
This problem is related to the permanent function. Assume first that each $\rho^{A_j}$ is a pure state $|\mathbf{x}_j\rangle\langle\mathbf{x}_j|$, where $\langle\mathbf{x}_j|\mathbf{x}_j\rangle=1$.  Then $\otimes_{j=1}^d \rho^{A_j}$ is a pure product state with the positive eigenvector $\otimes_{j=1}^d |\mathbf{x}_j\rangle$.  A symmetrization of  $\otimes_{j=1}^d |\mathbf{x}_j\rangle$ is the orthogonal projection on the subspace of symmetric tensors, given by
\begin{eqnarray*}
(\mathbb{I}-C^{B})(\otimes_{j=1}^d |\mathbf{x}_j\rangle)=\frac{1}{d!}\sum_{\pi\in\Sigma_d} \otimes_{j=1}^d |\mathbf{x}_{\pi(j)}\rangle.
\end{eqnarray*}
Hence
\begin{eqnarray*}
\big\|(\mathbb{I}-C^{B})(\otimes_{j=1}^d |\mathbf{x}_j\rangle)\big\|^2=\frac{1}{d!} \sum_{\pi\in\Pi_d}\prod_{j=1}^d \langle \mathbf{x}_j|\mathbf{x}_{\pi(j)}\rangle.
\end{eqnarray*}
Let $X=[\mathbf{x}_1\cdots \mathbf{x}_d]\in\mathbb{C}^{n\times d}$ be the matrix whose columns are the vectors $[\mathbf{x}_1,\ldots,\mathbf{x}_d]$.  The $G(\mathbf{x}_1,\ldots,\mathbf{x}_d)=X^\dagger X$ is the Gramian matrix $[\langle \mathbf{x}_i|\mathbf{x}_j\rangle]\in \mathrm{H}_{d,+}$.  Note that since $\|\mathbf{x}_1\|=\cdots=\|\mathbf{x}_d\|=1$ the diagonal entries of $G(\mathbf{x}_1,\ldots,\mathbf{x}_d)$ are all 1, and $G(\mathbf{x}_1, \ldots, \mathbf{x}_d)$  is called a complex covariance matrix.
It now follows that $\|(\mathbb{I}-C^{B})\otimes_{j=1}^d|\mathbf{x}_j\rangle\|^2$ is $\frac{1}{d!}$ times the permanent of $G(\mathbf{x}_1,\ldots,\mathbf{x}_d)$, denoted as per$\, G(\mathbf{x}_1,\ldots,\mathbf{x}_d)$.
Hence
\begin{align*}
\tr C^{B}(\otimes_{i=1}^d|\mathbf{x}_i\rangle)(\otimes_{j=1}^d\langle\mathbf{x}_i|)=1-\frac{1}{d!}\textrm{per}\,G(\mathbf{x}_1,\ldots,\mathbf{x}_d),&& \|\mathbf{x}_1\|=\cdots=\|\mathbf{x}_d\|=1.
\end{align*}
\begin{lemma}\label{trTprodps}
Assume that $\rho^{A_1},\ldots,\rho^{A_d}\in\Omega_n$ have the following spectral decomposition:
\begin{eqnarray*}
\rho^{A_j}=\sum_{i=1}^n \lambda_{i,j}|\mathbf{x}_{i,j}\rangle\langle\mathbf{x}_{i,j}|, \quad j\in[d].
\end{eqnarray*}
Then
\begin{equation}\label{trTprodps1}
\tr C^{B}(\otimes_{j=1}^d\rho^{A_j})=1-\frac{1}{d!}
\sum_{i_1\ldots,i_d\in[n]}\prod_{j=1}^d\lambda_{i_j,j}\; \mathrm{per}\,G(\mathbf{x}_{i_1,1},\ldots,\mathbf{x}_{i_d,d}).
\end{equation}
\end{lemma}
The proof of this lemma follows straightforwardly from the multilinearity of $\otimes_{j=1}^d\rho^{A_j}$.

We now state the analog to part (d) of Theorem \ref{kapTABprop} , which is a corollary to the above lemma:
\begin{corollary}\label{QTupbdcase}
Let $\rho^{A_1},\ldots,\rho^{A_d}$ be density matrices with the spectral decomposition given by Lemma \ref{trTprodps}.  Then
\begin{eqnarray*}
\mathrm{T}_{C^{B}}(\rho^{A_1},\ldots,\rho^{A_d})\le 1- \frac{1}{d!}
\sum_{i_1\ldots,i_d\in[n]}\prod_{j=1}^d\lambda_{i_j,j}\; \mathrm{per}\,G(\mathbf{x}_{i_1,1},\ldots,\mathbf{x}_{i_d,d}).
\end{eqnarray*}
If at least $d-1$ density matrices are pure states then equality holds.
\end{corollary}
\emph{Acknowledgements:}
It is a pleasure to thank W.~S{\l}omczy{\'n}ski
for inspiring discussions and helpful remarks.
Financial support by Simons collaboration grant for mathematicians, Narodowe Centrum Nauki
under the Maestro grant number DEC-2015/18/A/ST2/00274
and by Foundation for Polish Science
under the Team-Net project no. POIR.04.04.00-00-17C1/18-00
is gratefully acknowledged.
\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv