Thermodynamically optimal creation of correlations

Let us start our investigation with a detailed description of the problem: we consider the process of correlating two previously uncorrelated quantum systems. More precisely, we study the controlled interaction of two quantum systems A and B with Hamiltonians $\newcommand{\SA}{{_{A}}} H\SA$ and $\newcommand{\SB}{{_{B}}} H\SB$, respectively, which are initially in a tensor product state $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \varrho\SA\otimes\varrho\SB$, with the goal of increasing the correlations between them. As we want to obtain optimal bounds, we consider optimal coherently controlled systems, i.e. the ability to externally engineer and time any interaction Hamiltonian between the two quantum systems, resulting in unitary dynamics on the system and we thus study the global unitary orbit of product states. Furthermore, we are interested in the energy needed to establish correlations between these initially uncorrelated quantum systems.

One could generalise the question from perfect system control and thus unitaries on the system, to correlating maps generated by unitaries on the system and arbitrary auxiliary systems, inducing completely positive and trace preserving (CPTP) maps on the system. However, some transformations (e.g. cooling to the ground state) may only be achievable asymptotically (i.e. using infinite time or energy [21–24]). Moreover, the implementation of CPTP maps generally requires access to and control over auxiliary systems, resulting in an implicit energy cost for preparation and manipulation of the auxiliaries that is obscured by the use of the CPTP maps instead of the explicit description of the corresponding unitaries on the total Hilbert space. Consequently, it is of interest to understand the limits of unitary correlating processes, which allow for a complete account of all the energetic changes within the system. Moreover, to give a fair account of the energy that needs to be supplied to the system, it is important to consider initial states that do not implicitly supply free energy. That is, we assume that the initial states of the subsystems are thermal states at the same temperature T, i.e. $\newcommand{\SA}{{_{A}}} \varrho\SA=\tau\SA(\beta)$ and $\newcommand{\SB}{{_{B}}} \varrho\SB=\tau\SB(\beta)$ with $\beta=1/T$ (we use units where $k_{\mathrm{B}}=1$ throughout), where $\tau(\beta):={\rm e}^{-\beta H}/\mathcal{Z}$ and $\newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \mathcal{Z}=\tr({\rm e}^{-\beta H})$. In this way, the initial state $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \tau\SA(\beta)\otimes\tau\SB(\beta)$ minimises the local energies for given entropies and is also completely passive [25], i.e. a state from which no energy can be extracted unitarily even when taking multiple copies.

Correlations can be quantified in different ways. Qualitatively, one may describe correlations between two subsystems as the feature that more information is available about properties of the joint system than about properties of the subsystems. If there is no preference as to the specific properties that are to be correlated, i.e. if one is interested in quantifying any correlations, it is most useful to consider entropic measures. The most paradigmatic among them is the mutual information $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\SAB}{{_{AB}}} \mathcal{I}(\varrho\SAB)=S(\varrho\SA)+S(\varrho\SB)-S(\varrho\SAB)$, based on the von Neumann entropy $\newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\bi}{\boldsymbol} S(\varrho):=-\tr\bigl(\varrho\ln\varrho\bigr)$. Although we focus on the von Neumann entropy, in principle one can replace it with any other measure of local 'mixedness'. And indeed, our methods are framed in the context of majorisation relations, and thus in principle also open an avenue towards studying other entropic measures.

When focusing on the mutual information, the problem of unitarily creating correlations can be reduced to the task of maximising the sum of the marginal entropies under local energy constraints. That is, the invested energy $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\SAB}{{_{AB}}} \newcommand{\re}{{\rm Re}} \newcommand{\bi}{\boldsymbol} \renewcommand{\dag}{\dagger}\Delta E:={\rm Tr}\bigl[H\SAB\bigl(U\SAB\tau\SA(\beta)\,\otimes$ $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\SAB}{{_{AB}}} \newcommand{\re}{{\rm Re}} \newcommand{\bi}{\boldsymbol} \renewcommand{\dag}{\dagger} \tau\SB(\beta)U^{\dagger}\SAB-\tau\SA(\beta)\otimes\tau\SB(\beta)\bigr)\bigr]$, where $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\SAB}{{_{AB}}} H\SAB: = H\SA + H\SB,$ only depends on the local marginals, while the global entropy is invariant under unitary operations. The change in mutual information thus reduces to $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \Delta \mathcal{I}=\Delta S\SA+\Delta S\SB$. We thus arrive at the following question:

Question 1 (Maximal mutual information under energy constraint). For a given pair of local Hamiltonians $\newcommand{\SA}{{_{A}}} H\SA$ and $\newcommand{\SB}{{_{B}}} H\SB$ and initial temperature $T=1/\beta$, what is the maximum value of

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \Delta \mathcal{I}=\Delta S\SA+\Delta S\SB \label{Ques:DeltaMuUuuuuuh} \nonumber \end{align} \tag{ 1 }$$

for an invested energy of at most $\Delta E$?

The fact that thermal states maximise the local entropies under local energy constraints suggests that the maximal value of $\Delta \mathcal{I}$ is always achieved when the final state marginals are thermal at a higher temperature $T'\geqslant T$ (or, equivalently, for $\beta'\leqslant \beta$), i.e. $\newcommand{\SA}{{_{A}}} \tilde{\varrho}\SA=\tau\SA(\beta')$ and $\newcommand{\SB}{{_{B}}} \tilde{\varrho}\SB=\tau\SB(\beta')$, provided that the corresponding STUs exist. This is indeed the case as shown in appendix A.1 and is what leads us to the question of the existence of STUs. While when $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} H\SA=H\SB$, one might still hope for STUs to exist for all dimensions, Hamiltonians and temperatures, such a statement cannot be made when $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} H\SA\neq H\SB$. Before introducing our framework for constructing STUs, let us therefore briefly discuss general bounds on the achievable local entropies.

2.1.1. Asymmetric case.

Let us consider question 1 in the case $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} H\SA\neq H\SB$. Although the problem is in general complex, note that nontrivial constraints on marginal spectra for given global spectra exist in the form of entropy inequalities. In particular, both the von Neumann entropy $S(\varrho)$ and the Rényi 0-entropy $S_{0}(\varrho):=-\log_2({\rm rank}(\varrho))$ satisfy subadditivity [26, 27], while no other entropy satisfies (dimension independent) linear inequalities [28]. Since we are interested in thermodynamic questions, where the rank of the involved states is typically full, only subadditivity of the von Neumann entropy, which we can write as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SAB}{{_{AB}}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle S(\varrho\SAB) \geqslant |S(\varrho\SA)-S(\varrho\SB)|, \label{triangle ineq} \nonumber \end{align} \tag{ 2 }$$

provides a nontrivial constraint. Also observe that, given an initial state $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \tau\SA(\beta)\otimes \tau\SB (\beta)$ subject to a unitary evolution, the entropy of the total system is fixed to $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\SAB}{{_{AB}}} S(\varrho\SAB) = S(\tau\SA(\beta))+S(\tau\SB (\beta))$. If we assume that the final marginals are thermal states with equal higher temperature i.e. $\beta'\leqslant \beta$, inequality (2) can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle |S(\tau\SA(\beta^{\prime}))-S(\tau\SB (\beta^{\prime}))| \leqslant\, S(\tau\SA(\beta))+ S(\tau\SB (\beta)). \nonumber \end{align} \tag{ 3 }$$

This already implies a nontrivial constraint, and also provides a counterargument to the naive assumption that two thermal marginals at the same temperature are reachable in general.

While this makes the structure complicated to deal with in general, one can still solve the optimisation in some special cases. In particular, for a system initialised in the ground state, i.e. a product of two pure states $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\SAB}{{_{AB}}} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} \tau\SAB (\beta \to \infty)= \ket{0\SA\,0\SB}\!\bra{0\SA\,0\SB}$, the total state remains pure during the unitary evolution and hence has a Schmidt decomposition. The final state can therefore be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ket}[1]{| #1\rangle} \displaystyle \ket{\tilde{\psi}} = \sum_{i=0}^{d} \sqrt{p_i}\ket{\varphi^{A}_{i}, \varphi^{B}_{i}}, \label{eq:schmidt decomp} \nonumber \end{align} \tag{ 4 }$$

for some probabilities p _i, with $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} d= \mathrm{min}\{d\SA,d\SB\}$, and $\newcommand{\SA}{{_{A}}} d\SA$ and $\newcommand{\SB}{{_{B}}} d\SB$ the dimensions of the respective Hilbert spaces $\newcommand{\SA}{{_{A}}} \mathcal{H}\SA$ and $\newcommand{\SB}{{_{B}}} \mathcal{H}\SB$. And the final marginals are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\langle #1 |} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tilde{\varrho}\SA =\sum_{i=0}^{d-1} p_i \ket{\varphi^{A}_{i}}\!\bra{\varphi^{A}_{i}},\ \ \ \tilde{\varrho}\SB= \sum_{i=0}^{d-1} p_i \ket{\varphi^{B}_{i}}\!\bra{\varphi^{B}_{i}}. \label{substateA} \nonumber \end{align} \tag{ 5 }$$

This already means that the marginals (since they have equal rank) cannot both be thermal states if $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} d\SA\neq d\SB$.

However, in this case it is still possible to maximise the amount of correlation subject to a given maximal amount of energy consumption. In other words, one can find optimal points of the following optimisation problem:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \label{eq:puremaxS} &\max_{\tilde{\varrho}\SA, \tilde{\varrho}\SB} S(\tilde{\varrho}\SA)+S(\tilde{\varrho}\SB)\nonumber \\ &{\rm subject~ to}~ \tr(\tilde{\varrho}\SA H\SA)+ \tr(\tilde{\varrho}\SB H\SB) \leqslant c\nonumber \\ & \hphantom{{\rm subject~ to}}\ \tilde{\varrho}\SA\ ~{\rm and}~ \tilde{\varrho}\SB\ ~{\rm as~ in~ equation~(5).}\nonumber \end{align} \tag{ 6 }$$

To do so one first solves the (loosely speaking) inverse problem of minimising the energy consumption for a given amount of correlation. This problem can be turned into two instances of the well-known passivity problem [25], for which passive states are the solutions. Using this fact, the problem reduces to a simple problem that can be solved using Lagrange multipliers. One then shows that the solution found for that problem is strictly monotonically increasing in the constraint, allowing us to reverse the problem and show that it is a solution of the problem in (6). For more details see appendix A.2. See also [29] for an alternative approach to the energy minimisation problem. The solutions obtained for (6) are thermal states at inverse temperature $\beta(c)$ (uniquely defined by the constant c) of the Hamiltonians $\newcommand{\SA}{{_{A}}} \tilde{H}\SA$ and $\newcommand{\SB}{{_{B}}} \tilde{H}\SB$, respectively, defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\langle #1 |} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tilde{H}\SA&= \sum_{i=0}^{d-1} (E_{i}^{A}+E_{i}^{B}) \ket{i}\!\bra{i}\SA\nonumber \\ \tilde{H}\SB&= \sum_{i=0}^{d-1} (E_{i}^{A}+E_{i}^{B}) \ket{i}\!\bra{i}\SB, \nonumber \end{align} \tag{ 7 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\langle #1 |} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle H\SA&= \sum_{i=0}^{d\SA-1} E_{i}^{A} \ket{i}\!\bra{i}\SA, \quad E_{i}^{A} \leqslant E_{i+1}^{A},\nonumber \\ H\SB&= \sum_{i=0}^{d\SB-1} E_{i}^{B} \ket{i}\!\bra{i}\SB, \quad E_{i}^{B} \leqslant E_{i+1}^{B}. \nonumber \end{align} \tag{ 8 }$$

The solutions hence take the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tilde{\varrho}_{\text{opt},\tiny{A}} (\beta(c)) &= \frac{\tilde{\Pi}\SA e^{-\beta(c) \tilde{H}\SA}\tilde{\Pi}\SA}{\tr(\tilde{\Pi}\SA e^{-\beta(c) \tilde{H}\SA})},\nonumber \\ \tilde{\varrho}_{\text{opt},\tiny{B}} (\beta(c)) &= \frac{\tilde{\Pi}\SB e^{-\beta(c) \tilde{H}\SB}\tilde{\Pi}\SB}{\tr(\tilde{\Pi}\SB e^{-\beta(c) \tilde{H}\SB})}. \nonumber \end{align} \tag{ 9 }$$

where .

The maximum amount of correlation achievable subject to a given maximal amount of energy c is then

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \newcommand{\SA}{{_{A}}} \displaystyle 2 S(\tilde{\varrho}\SA) = 2 \beta(c) c+ 2 \ln \left[\tr({\rm e}^{-\beta(c) \tilde{H}\SA})\right]. \label{finmut1} \nonumber \end{align} \tag{ 10 }$$

2.1.2. Symmetric case.

Let us now consider the symmetric case $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} H\SA=H\SB$. Following our considerations in the previous section and appendix A.1, it is clear that an upper bound for $\Delta\mathcal{I}$ is in this case achieved when $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \Delta S\SA=\Delta S\SB$, and since the states with maximal entropy given a fixed energy are thermal states, question 1 can be substantially simplified to that of the existence of symmetrically thermalizing unitaries.

Question 2 (Existence of STUs). Does there exist a unitary $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} U\SAB$ on $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \mathcal{H}\SAB$ such that

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\SAB}{{_{AB}}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tilde{\varrho}\SA &=\tr\SB\bigl(U\SAB\tau\SAB(\beta) U\SAB^{\dagger}\bigr)=\tau\SA(\beta'),\nonumber \\ \tilde{\varrho}\SB &=\tr\SA\bigl(U\SAB\tau\SAB(\beta) U\SAB^{\dagger}\bigr)=\tau\SB(\beta'), \label{Ques:Uandmarginals} \nonumber \end{align} \tag{ 11 }$$

for every pair of local Hamiltonians $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} H\SA=H\SB$, for all final temperatures $T'=1/\beta'$ and all initial temperatures $T=1/\beta\leqslant T'$?

If such STUs exist, then they are automatically the optimally correlating unitaries for a given system. It is already known that question 1 can be answered affirmatively when the subsystem Hamiltonians are equally spaced, i.e. $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} H\SA=H\SB=\sum\nolimits_{n=0}^{d-1}E_{n}\ket{n}\!\bra{n}$ with $E_{n+1}-E_{n}=\omega ~\forall n\in\{0,1,\ldots,d-1\}$ for some constant $\omega$ (with appropriate units) [11] (see also the more recent formulation of the derivation in [30]). In particular, this implies that such optimally correlating unitaries always exist for two qubits, i.e. when d  =  2. However, there is (yet) no answer for the general situation of arbitrary Hamiltonians and dimensions. To fill this gap, we present a framework that generalises previous approaches [11] and allows investigating this general question in the symmetric cases for any dimension. This approach is based on unitary operations in LCSs, as we will explain in section 2.2, and enables us to show that STUs exist for the simplest nontrivial cases of two qutrits (d  =  3) and two ququarts (d  =  4).

Let us now discuss our general framework for marginal transformations. The central idea of this approach is to decompose the diagonal elements of the marginals into elements originating from different subspaces, with the property that any unitary that leaves the division into these subspaces invariant never creates local off-diagonal elements.

More precisely, we consider a pair of d-dimensional systems A and B with matching local Hamiltonians $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} H\SA=H\SB=\sum\nolimits_{i=0}^{d-1}E_{i}\ket{i}\!\bra{i}$, where we set E₀  =  0 without loss of generality. Moreover, let the energy eigenvalues, sorted in non-decreasing order (with $E_{i+1}\geqslant E_{i}\ \forall\,i$), be measured in units of E₁. The initial thermal states of A and B are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\langle #1 |} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tau\SA(\beta) = \tau\SB(\beta)\, = \sum_{i=0}^{d-1}\,p_{i}(\beta)\,\ket{i}\!\bra{i}, \label{eq:thermal state} \nonumber \end{align} \tag{ 12 }$$

where $p_{i}(\beta)={\rm e}^{-\beta E_i}/\mathcal{Z}(\beta)$ and $\mathcal{Z}(\beta)=\sum\nolimits_{i}{\rm e}^{-\beta E_{i}}$. For convenience, we use the shorthand $\newcommand{\e}{{\rm e}} p_{i}\equiv p_{i}(\beta)$. The joint initial state $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\SAB}{{_{AB}}} \tau\SAB=\tau\SA(\beta)\otimes\tau\SB (\beta)$ is also diagonal with respect to the energy eigenbasis and its diagonal entries are products of the diagonal entries of the reduced states, i.e. we write $p_{ij}:=p_{i}p_{j}$ such that $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\ketbra}[2]{{|{#1}\rangle\!\langle{#2}|}} \tau\SAB=\sum\nolimits_{i,j=0}^{d-1}p_{ij}\,\ketbra{ij}{ij}$. Since the unitaries that we use throughout the manuscript do not create coherence in the marginals, it is convenient to introduce a vectorised notation for the diagonal entries of the reduced states. That is, for arbitrary diagonal joint states $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \newcommand{\diag}{{\rm diag}} {\bf p}\SAB=\diag\{\,p_{ij}\}_{i,j=0}^{d-1}$ where the p _ij need not factorise with respect to i and j , the reduced states $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SAB}{{_{AB}}} \varrho\SA=\tr\SB(\varrho\SAB)$ and $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SAB}{{_{AB}}} \varrho\SB=\tr\SA(\varrho\SAB)$ are diagonal, and the diagonal entries can be collected into vectors $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} {\bf p}\SA,{\bf p}\SB\in\mathbb{R}^{d}$ with nonnegative components and unit one-norm, $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\dv}{{||}} \dv {\bf p}\SA\dv=\tr(\varrho\SA)=1=\dv {\bf p}\SB\dv$.

2.2.1. Locally classical subspaces.

Having expressed the diagonals of the marginals in this vectorised form, we further choose particular subspace vector decompositions, as we illustrate for d  =  3 in section 2.2.2. For general local dimension $d\geqslant3$, we write $\newcommand{\SA}{{_{A}}} {\bf p}\SA$ and $\newcommand{\SB}{{_{B}}} {\bf p}\SB$ as sums of d vectors according to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle {\bf p}\SA = \sum_{i=0}^{d-1}{\bf r}^{A}_{i} = \sum_{i=0}^{d-1}{\bf r}_{i}, \ \ \ {\bf p}\SB = \sum_{i=0}^{d-1}{\bf r}^{B}_{i}=\sum_{i=0}^{d-1}\Pi^i\,{\bf r}_{i}, \label{deco. marginal} \nonumber \end{align} \tag{ 13 }$$

where $\newcommand{\re}{{\rm Re}} {\bf r}^{A}_{i} ={\bf r}_{i}= \sum\nolimits_{j=0}^{d-1} p_{j\, j+i}\,\mathbf {e}_{j}$ and $\newcommand{\re}{{\rm Re}} {\bf r}^{B}_{i}=\sum\nolimits_{j=0}^{d-1} p_{j-i\, j}\,\mathbf {e}_{j}$. Here, $\newcommand{\re}{{\rm Re}} \{\mathbf {e}_{j}\}_{j=0}^{d-1}$ denotes the chosen orthonormal basis of $\mathbb{R}^{d}$ and all indices are to be taken modulo d. In the second equality, we have further used the fact that the ith vector ${\bf r}^{B}_{i}$ in the decomposition of $\newcommand{\SB}{{_{B}}} {\bf p}\SB$ can be related to the ith vector ${\bf r}^{A}_{i}$ in the decomposition of $\newcommand{\SA}{{_{A}}} {\bf p}\SA$ via the ith power of the cyclic permutation matrix $\Pi=(\Pi_{ij})$ with components $\Pi_{ij} = \delta_{i,j+1~{\rm mod}~d}$.

One further observes that the decomposition of $\newcommand{\SA}{{_{A}}} {\bf p}\SA$ and $\newcommand{\SB}{{_{B}}} {\bf p}\SB$ into these vectors corresponds to the selection of a total of d subspaces $\mathcal{H}_{q}$ and $\{\mathcal{H}_{r_i}\}_{i=1}^{d-1}$ of the joint Hilbert space $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \mathcal{H}\SAB$,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ket}[1]{| #1\rangle} \displaystyle \mathcal{H}_{q} := \mathcal{H}_{r_0} = {\rm span}\{\ket{j\,j}\}_{j=0}^{d-1},\ \ \ \mathcal{H}_{r_i} = {\rm span}\{\ket{j\,j+i}\}_{j=0}^{d-1}, \label{eq:subspace dec} \nonumber \end{align} \tag{ 14 }$$

with $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \mathcal{H}\SAB=\mathcal{H}_{q}\oplus_{i=1}^{d-1}\mathcal{H}_{r_i}$, such that arbitrary unitaries U_q and $U_{r_i}$ applied in either of the d subspaces cannot lead to nonzero off-diagonal elements in the reduced states $\newcommand{\SA}{{_{A}}} \varrho\SA$ or $\newcommand{\SB}{{_{B}}} \varrho\SB$, provided that none are present to begin with. In other words, any state on $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \mathcal{H}\SAB$ that has support on only one of these subspaces, or which has a block-diagonal structure with respect to this subspace decomposition, is locally classical, i.e. has diagonal marginals with respect to the local bases $\newcommand{\ket}[1]{| #1\rangle} \{\ket{j}\}_{j=0}^{d-1}$. We hence call these subspaces locally classical. However, note that the choice of LCSs we make here is not unique, as discussed below in section 2.2.3.

A unitary transformation $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\varrho\SAB\mapsto\tilde{\varrho}\SAB=U\varrho\SAB U^{\dagger}$ with $U=U_{q}\oplus_{i=1}^{d-1} U_{r_i}$ hence does not preclude off-diagonal elements from appearing in the joint state $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \tilde{\varrho}\SAB$ but leaves the reduced states diagonal, implying that we can directly read off the spectra of the marginals. In such a case, it is still useful to describe the marginals by d-component vectors $\newcommand{\re}{{\rm Re}} \newcommand{\SA}{{_{A}}} {\tilde{\mathbf {p}}}\SA$ and $\newcommand{\re}{{\rm Re}} \newcommand{\SB}{{_{B}}} {\tilde{\mathbf {p}}}\SB$ collecting their diagonal elements, and the transformation of these vectors can be represented as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SA}{{_{A}}} \displaystyle {\bf p}\SA \mapsto\,\tilde{{\bf p}}\SA = M_{q}{\bf q}+\sum_{i=1}^{d-1}M_{r_i}{\bf r}_{i},\label{eq: gen evo. marginal A} \nonumber \end{align} \tag{ 15 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \displaystyle {\bf p}\SB \mapsto\,\tilde{{\bf p}}\SB = M_{q}{\bf q}+\sum_{i=1}^{d-1}\Pi^i M_{r_i}{\bf r}_{i}, \label{gen. evolution marginal} \nonumber \end{align} \tag{ 16 }$$

where ${\bf q}={\bf r}^{A}_{0}={\bf r}^{B}_{0}$, and M_q and $M_{r_i}$ are unistochastic $d\times d$ matrices, i.e. doubly-stochastic (entries in rows and columns sum to 1) square matrices $M_{\alpha}$ (with $\alpha\in\{q,r_i\}$) whose components can be understood as the moduli squared of unitary matrices, $(M_{\alpha})_{kl}=|(U_{\alpha})_{kl}|^{2}$.

As a technical remark, note that the Schur–Horn theorem implies a weak convexity property for unistochastic matrices, namely that for any $\newcommand{\re}{{\rm Re}} \mathbf {v}\in\mathbb{R}^{d}$, the set of vectors obtained by applying the set of unistochastic $d\times d$ matrices to $\newcommand{\re}{{\rm Re}} \mathbf {v}$ is equivalent to the set of vectors obtained by applying the set of doubly-stochastic $d\times d$ matrices to $\newcommand{\re}{{\rm Re}} \mathbf {v}$, see [31]. This means that for any doubly-stochastic matrix M and vector $\newcommand{\re}{{\rm Re}} \mathbf {v}$, there exists a unistochastic matrix $\newcommand{\re}{{\rm Re}} M_{\mathbf {v}}$ such that $\newcommand{\re}{{\rm Re}} M\,\mathbf {v}=M_{\mathbf {v}}\mathbf {v}$. Here, this property permits us to conclude that the vectorised marginals of any state reachable by unitaries of the form $U=U_{q}\oplus_{i=1}^{d-1} U_{r_i}$ can be written in terms of the action of doubly-stochastic matrices M_q and $M_{r_{i}}$ on ${\bf q}$ and ${\bf r}_{i}$, respectively.

In general, equation (15) can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SA}{{_{A}}} \displaystyle \tilde{{\bf p}}\SA = \left\{\begin{array}{@{}ll@{}} M_{q}{\bf q} + \sum\limits_{i=1}^{k} (M_{r_{i}}{\bf r}_{i} + M_{r_{d-i}}{\bf r}_{d-i}) & d\ ~{\rm odd} \nonumber \\ M_q{\bf q} + \sum\limits_{i=1}^{k-1} (M_{r_{i}}{\bf r}_{i} + M_{r_{d-i}}{\bf r}_{d-i})+M_{r_{d/2}} {\bf r}_{d/2} & d\ ~{\rm even} \end{array}\right., \label{gen evol marginal A} \nonumber \end{align} \tag{ 17 }$$

where $k=\frac{d-1}{2}$ if d is odd and $k=\frac{d}{2}$ if d is even. Further taking into account the symmetry $p_{ij} = p_{ji}$, which implies $\newcommand{\re}{{\rm Re}} {\bf r}_{d-i}= \sum\nolimits_{j=0}^{d-1} p_{j\, j+d-i} \Pi^{i}\mathbf {e}_{j-i}=\Pi^{i} {\bf r}_{i}$, the vectorised form of the marginal A can be written in a compact way as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SA}{{_{A}}} \displaystyle \tilde{{\bf p}}\SA = M_{q}{\bf q} + \sum\limits_{i=1}^{k} \left(\lfloor\frac{2i}{d}\rfloor+1\right)^{-1}(M_{r_i} + M_{r_{d-i}}\Pi^i){\bf r}_{i}, \label{eq:gen evo. marg. A 1} \nonumber \end{align} \tag{ 18 }$$

where $\lfloor x \rfloor$ denotes the floor function of x, and the prefactor $(\lfloor\frac{2i}{d}\rfloor+1){}^{-1}$ is equal to 1 unless d is even and i  =  k, in which case it is $\frac{1}{2}$. Using equation (16), one may obtain the transformed marginal B as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \displaystyle \tilde{{\bf p}}\SB = M_q{\bf q} + \sum\limits_{i=1}^{k} \left(\lfloor\frac{2i}{d}\rfloor+1\right)^{-1} (\Pi^i M_{r_i} + \Pi^{-i}M_{r_{d-i}}\Pi^i)\, {\bf r}_{i}, \nonumber \end{align} \tag{ 19 }$$

where we have used the property $\Pi^{-i}=\Pi^{d-i}$ of the d-dimensional cyclic permutation. To satisfy the conditions of question 2, we need to further restrict ourselves to a subgroup of transformations which generate the same marginals, $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \tilde{{\bf p}}\SA=\tilde{{\bf p}}\SB$. This requirement results in the condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{r_{d-i}}=\Pi^{i} M_{r_i} \Pi^{-i}\ \ \ \ {\rm for}\ \ \ i=1,\dots,k, \nonumber \end{align} \tag{ 20 }$$

for the doubly-stochastic matrices mentioned in equation (18). With this, we arrive at the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}}={\tilde{\mathbf {p}}}\SA={\tilde{\mathbf {p}}}\SB = M_q {\bf q} + \sum_{i=1}^{k} \left(\lfloor\frac{2i}{d}\rfloor+1\right)^{-1} (\mathbb{1}+ \Pi^i)M_{r_i} {\bf r}_i \label{eq:EqualEvoMarginal C} \nonumber \end{align} \tag{ 21 }$$

for both vectorised marginals. Now that we have ensured that both marginals transform in the same way, we may concentrate on one of them, say $\newcommand{\SA}{{_{A}}} \tilde{{\bf p}}\SA$ (dropping the subscript A for ease of notation) and use the established framework to investigate the existence of STUs as specified in question 2.

2.2.2. Box 1 (Unitary operations on locally classical subspaces for two qutrits).

For any tensor product of two identical states $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\ketbra}[2]{{|{#1}\rangle\!\langle{#2}|}} \tau\SAB=\sum\nolimits_{i,j=0}^{2}p_{ij}\,\ketbra{ij}{ij}$, we can separate the local diagonal elements (i.e. the eigenvalues, since the matrix is diagonal) through an economic notation where we denote the diagonal entries of the marginals as entries in vectors

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \newcommand{\re}{{\rm Re}} \displaystyle {\bf p}\SA:= \mathbf {\tau}\SA= \left(\begin{array}{@{}c@{}} {p_{00}}+{p_{01}}+{p_{02}} \nonumber \\ {p_{11}}+{p_{12}}+{p_{10}} \nonumber \\ {p_{22}}+{p_{20}}+{p_{21}}\end{array}\right)={{\bf q}}+{{\bf r}_1}+{{\bf r}_2}, \quad {\bf p}\SB:=\mathbf {\tau}\SB= \left(\begin{array}{@{}c@{}} {p_{00}}+{p_{10}}+{p_{20}} \nonumber \\ {p_{11}}+{p_{21}}+{p_{01}} \nonumber \\ {p_{22}}+{p_{02}}+{p_{12}}\end{array}\right)={{\bf q}}+{\Pi^2\,{\bf r}_2}+{\Pi\,{\bf r}_1}, \nonumber \end{align} \tag{ 22 }$$

where $\Pi=(\Pi_{ij})$ is a cyclic permutation matrix with components $\Pi_{ij} = \delta_{i,j+1~{\rm mod}~3}$. Employing the unitary $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} U\SAB={U_q}\oplus{U_{r_1}}\oplus {U_{r_2}}$, where U_q, $U_{r_{1}}$, and $U_{r_{2}}$ act unitarily on the subspaces

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ket}[1]{| #1\rangle} \displaystyle {\mathcal{H}_{q}:=\mathrm{span}\{\ket{00},\ket{11},\ket{22}\}},\quad {\mathcal{H}_{r_1}:=\mathrm{span}\{\ket{01},\ket{12},\ket{20}\}},\quad {\mathcal{H}_{r_2}:=\mathrm{span}\{\ket{02}, \ket{10},\ket{21}\}}, \nonumber \end{align} \tag{ 23 }$$

respectively, the transformation of the marginals can be described by unistochastic matrices

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}}\SA={M_q\,{\bf q}}+{M_{r_1}{\bf r}_1}+{M_{r_2}{\bf r}_2}, \quad {\tilde{\mathbf {p}}}\SB={M_{q}{\bf q}}+{\Pi^2\,M_{r_2}{\bf r}_2}+{\Pi \,M_{r_1}{\bf r}_1}. \nonumber \end{align} \tag{ 24 }$$

The components of the unistochastic matrices are determined by the moduli squared of corresponding unitary matrices, $(M_{\alpha})_{kl}=|(U_{\alpha})_{kl}|^{2}$ for $\alpha=q,r_{1},r_{2}$. Due to the symmetry $p_{ij}=p_{ji}$, we can further identify ${{\bf r}_2}=\Pi {{\bf r}_1}$ to obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}}\SA={M_q\,{\bf q}}+{M_{r_1}{\bf r}_1}+{M_{r_2}\Pi {\bf r}_1}, \quad {\tilde{\mathbf {p}}}\SB={M_{q}{\bf q}}+{\Pi^2\,M_{r_2}\Pi\,{\bf r}_1}+{\Pi\,M_{r_1}{\bf r}_1}. \nonumber \end{align} \tag{ 25 }$$

2.2.3. General locally classical subspaces.

While we focus here on the decomposition into the specific LCSs from equation (14), we note that this is by far not the only option. Indeed, let $\{\mathcal{P}_{i}\}_{i=0,\dots,d-1}$ be a set of permutations of d elements $j=0,\dots,d-1$. Then a sufficient condition for obtaining an LCS decomposition $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \newcommand{\bi}{\boldsymbol} \mathcal{H}\SAB=\bigoplus\nolimits_{i=0}^{d-1}\tilde{\mathcal{H}}_{\tilde{r}_i}$ into LCSs of the form $\newcommand{\ket}[1]{| #1\rangle} \tilde{\mathcal{H}}_{\tilde{r}_i}={\rm span}\{\ket{j\,\mathcal{P}_{i}(\,j)}\}_{j=0}^{d-1}$ is that the $d\times d$ matrix $\Gamma$ with components $\Gamma_{ij}=\mathcal{P}_{i}(\,j)$ is a Latin square, i.e. every entry $j=0,\dots,d-1$ appears exactly once in each row and in each column. For every such choice of LCS decomposition, one may then apply unitaries $\newcommand{\bi}{\boldsymbol} \tilde{U}=\bigoplus\nolimits_{i=0}^{d-1}U_{\tilde{r}_i}$ that leave the LCSs invariant, i.e. where $U_{\tilde{r}_i}$ acts nontrivially only on the subspace $\tilde{\mathcal{H}}_{\tilde{r}_i}$. Denoting the corresponding vector decomposition of the first vectorised marginal as $\newcommand{\SA}{{_{A}}} {\bf p}\SA=\sum\nolimits_{i=0}^{d-1}\tilde{{\bf r}}_{i}$, with $(\tilde{{\bf r}}_{i})_{j}=p_{j\mathcal{P}_{i}(\,j)}$, we have $\newcommand{\SB}{{_{B}}} {\bf p}\SB=\sum\nolimits_{i=0}^{d-1}\mathcal{P}_{i}^{-1}\tilde{{\bf r}}_{i}$ since $({\bf r}_{i}^{B})_{j}=p_{\mathcal{P}^{-1}_{i}(\,j){j}}=(\mathcal{P}_{i}^{-1}\tilde{{\bf r}}_{i})_{j}$. In this case, unitaries that leave the LCSs invariant transform the vectorised marginals according to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SA}{{_{A}}} \newcommand{\re}{{\rm Re}} \displaystyle {\bf p}\SA \mapsto\,{\tilde{\mathbf {p}}}\SA = \sum_{i=0}^{d-1}M_{\tilde{r}_i}\tilde{{\bf r}}_{i}, \label{eq: gen evo. marginal A general} \nonumber \end{align} \tag{ 26 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\re}{{\rm Re}} \displaystyle {\bf p}\SB \mapsto\,{\tilde{\mathbf {p}}}\SB = \sum_{i=0}^{d-1}\mathcal{P}^{-1}_{i} M_{\tilde{r}_i}\tilde{{\bf r}}_{i}, \label{gen. evolution marginal general} \nonumber \end{align} \tag{ 27 }$$

where the matrices $M_{\tilde{r}_i}$ are the unistochastic matrices corresponding to the unitaries $U_{\tilde{r}_i}$. This implies the existence of symmetric marginal transformations, e.g. for unistochastic matrices M that commute with all the $\mathcal{P}_{i}$, such that the marginals take the form $\newcommand{\re}{{\rm Re}} \newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} M{\bf p}\SA =M{\bf p}\SB={\tilde{\mathbf {p}}}$. For the remainder of this work we focus again on the specific special case we consider in equation (14), where $\mathcal{P}_{i}=(\Pi^{-1}){}^{i}$.

2.2.4. Asymptotic case.

In the particular case of d  =  3 the final expression for the (equal) marginals from equation (21) is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}}={\tilde{\mathbf {p}}}\SB = {\tilde{\mathbf {p}}}\SA = M_{q}{\bf q}+(\mathbb {1}+\Pi)M_{r}{\bf r}, \label{eq:gen. trans. 3d} \nonumber \end{align} \tag{ 28 }$$

where, as in equation (16), we have introduced the vectors ${\bf q}$ and ${\bf r}$, with components

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf q} = \left(\begin{array}{@{}c@{}} p_{00},\ p_{11},\ p_{22} \end{array}\right)^T,\ \ {\bf r} = \left(\begin{array}{@{}c@{}} p_{01},\ p_{12},\ p_{20} \end{array}\right)^T. \label{vecqrs} \nonumber \end{align} \tag{ 29 }$$

Based on the HLP theorem and equation (28), we are going to prove the existence of STUs in d  =  3 via the following lemma:

Lemma 1. For any $3\times 3$ doubly-stochastic matrix M, there exists a $3\times 3$ doubly-stochastic matrix $\tilde{M}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:commutingM} M (\mathbb{1}+\Pi)=(\mathbb{1}+\Pi) \tilde{M}. \nonumber \end{align} \tag{ 30 }$$

The proof of lemma 1 is presented in appendix A.4. We are now ready to state the first theorem for d  =  3.

Theorem 1 (Majorised marginals in d  =  3). For every pair of states $\varrho$ and $\tilde{\varrho}$ in a 3-dimensional Hilbert space which satisfy the condition $\newcommand{\re}{{\rm Re}} \mathbf {\lambda}(\tilde{\varrho})\prec \mathbf {\lambda}({\varrho})$, where $\newcommand{\re}{{\rm Re}} mathbf \mathbf {\lambda}(\varrho)$ is the vector of eigenvalues of $\varrho$, there exists a unitary $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} U\SAB$ on $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \mathcal{H}\SAB$ such that

$$\begin{align*} \renewcommand{\dag}{\dagger} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\SAB}{{_{AB}}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SA}{{_{A}}} \displaystyle \tilde{\varrho}_{{{A}}} &=\tr_{{{B}}}\bigl(U\SAB\varrho \otimes \varrho U\SAB^{\dagger}\bigr)= \tilde{\varrho},\nonumber \\ \tilde{\varrho}_{{{B}}}&=\tr_{{{A}}}\bigl(U\SAB\varrho \otimes \varrho U\SAB^{\dagger}\bigr)=\tilde{\varrho}. \nonumber \end{align*}$$

Proof. The unitary is given by $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} U\SAB=U_{\rm loc}U_{\rm ent}$, i.e. by a product of an entangling unitary $U_{\rm ent}$ which acts on the global system preserving the equal marginals as discussed in section 2.2, and a local unitary of the form $U_{\rm loc}=U\otimes U$, such that the marginals are kept equal. From lemma 1 we then see that, given any doubly-stochastic matrix M_q, one can find M_r such that $M_q (\mathbb{1} + \Pi)= (\mathbb{1}+\Pi) M_r$ and consequently a mapping to each probability vector $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {p}}}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}}= M_q {\bf q}+ (\mathbb{1}+\Pi) M_r {\bf r} = M_q\, [{\bf q}+ (\mathbb{1}+\Pi){\bf r}]=M_q \, {\bf p}. \nonumber \end{align} \tag{ 31 }$$

From the HLP theorem, we thus know that it is possible to produce all the vectors $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {p}}}$ that are majorised by the vector ${\bf p}$ corresponding to the initial state. Thus, we proved that through the entangling unitary, all the marginals $\newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} {\varrho^\prime}=\sum\nolimits_{i=0}^{2}\tilde{p}_i \ket{i}\!\bra{i}$ such that $\newcommand{\re}{{\rm Re}} mathbf {\tilde{\vmathbf {p}}}\prec \mathbf {p}$ can be reached from the initial state. Then, if we can also arbitrarily change their eigenstates in a symmetric way, theorem 1 is proven. For that, we use a local unitary $U \otimes U$ with $\newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} U =\sum\nolimits_{i=0}^2 \ket{\phi_i}\!\bra{i}$ to change the eigenstates of the marginal to an arbitrary set of eigenstates $\newcommand{\ket}[1]{| #1\rangle} \{\ket{\phi_i}\}_{i=0}^{2}$, such that finally the marginals become $\newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} \tilde{\varrho}=\sum\nolimits_{i=0}^{2}\tilde{p}_i \ket{\phi_i}\!\bra{\phi_i}$ in which the set $\newcommand{\ket}[1]{| #1\rangle} \{\ket{\phi_i}\}$ can be any orthonormal basis in d  =  3. □

As a corollary of theorem 1, the existence of STUs is proven for the two-qutrit case for initial and final thermal states with inverse temperatures $\beta$ and $\beta'\leqslant \beta$, respectively, since $\newcommand{\re}{{\rm Re}} {\bf p}(\beta')\prec \mathbf {p}(\beta)$ holds whenever $\beta' \leqslant\beta$.

Unfortunately, generalisations of lemma 1 to higher dimensions fail, as we explain in more detail in appendix A.5. Despite the general statement of theorem 1, any attempts at generalisations must hence be based on a different approach, which does not require the existence of a doubly-stochastic matrix $\tilde M$ as in equation (30) for all doubly-stochastic matrices M.

Here, we explore another possible approach that makes use of the HLP theorem and tries to overcome the difficulties encountered for generalising the previous approach. Let us first discuss the method for two d-dimensional systems, before turning to the special cases d  =  3 and d  =  4.

To begin, we once again employ the transformation presented in equation (21) to generate equal marginals. Recalling that the marginals of our final target state are thermal states at inverse temperature $\beta'<\beta$, we decompose the final state vector into the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}} = \mathbf {a} + \sum_{i=1}^{k} (\lfloor\frac{2i}{d}\rfloor+1)^{-1} (\mathbb{1}+ \Pi^i) \mathbf {b}_{i}, \label{eq:EqualEvoMarginal Cprime} \nonumber \end{align} \tag{ 32 }$$

with $\newcommand{\re}{{\rm Re}} \mathbf {a}:= {\bf q}(\beta^\prime)= \sum\nolimits_{j=0}^{d-1} p_{j{j}}(\beta^\prime)\,\mathbf {e}_{j}$ and $\newcommand{\re}{{\rm Re}} \mathbf {b}_i: = {\bf r}_i(\beta^\prime)= \sum\nolimits_{j=0}^{d-1} p_{j {j+i}}(\beta^\prime)\,\mathbf {e}_{j}$. Further, recall that the initial states have the same decomposition with inverse temperature $\beta$. One way of transforming ${\bf p}$ into $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {p}}}$ is to transform ${\bf q}$ into $\newcommand{\re}{{\rm Re}} \mathbf {a}$ and ${\bf r}_i$ into $\newcommand{\re}{{\rm Re}} \mathbf {b}_i$ for all $i=1,\dots,k$. Looking at equation (21), one realises that this is possible if

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle M_q {\bf q} =\mathbf {a},\label{eq:intuitivePassNorm} \nonumber \end{align} \tag{ 33 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle M_{r_i} {\bf r}_i =\mathbf {b}_i, \quad i=1,\dots,k. \label{eq:intuitivePassNorm2} \nonumber \end{align} \tag{ 34 }$$

Unfortunately, for any $\beta'<\beta$ we have $\newcommand{\re}{{\rm Re}} \Vert {\bf q} \Vert > \Vert \mathbf {a} \Vert$, where $\Vert \cdot \Vert$ denotes the one-norm, see appendix A.7, and since doubly-stochastic transformations conserve the norm of vectors, equation (33) cannot hold true. The norms of the ${\bf r}_i$ also generically differ from those of the $\newcommand{\re}{{\rm Re}} \mathbf {b}_i$, meaning equation (34) does not hold either. However, it is not clear whether $\Vert {\bf r}_i(\beta) \Vert$ is a monotonic function of $\beta$, and so $\Vert {\bf r}_i \Vert$ may in principle be smaller or larger than $\newcommand{\re}{{\rm Re}} \Vert \mathbf {b}_i \Vert$, depending on $\beta$, $\beta'$ and i. Loosely speaking, the vector ${\bf q}$ has 'too much norm'. Intuitively, one may thus make use of the excessive norm of ${\bf q}$ by transforming the required 'amount' of ${\bf q}$ into $\newcommand{\re}{{\rm Re}} \mathbf {a}$, and redistributing its excessive norm to the rest of $\tilde{{\bf p}}$, or in other words by passing the excessive norm of ${\bf q}$ on to the rest of $\tilde{{\bf p}}$, giving rise to the name of the approach. Formally, one achieves this by splitting M_q into a convex combination of two doubly-stochastic matrices, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_q= \alpha_0 M_{q \rightarrow a}+ \alpha_1 M_{q \rightarrow b}, \nonumber \end{align} \tag{ 35 }$$

where $\alpha_{0}, \alpha_{1} \geqslant 0$, $\alpha_0 +\alpha_1 =1$, and $M_{q \rightarrow a}$ and $M_{q \rightarrow b}$ are doubly-stochastic matrices such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle M_{q \rightarrow a} \frac{{\bf q}}{\Vert {\bf q} \Vert}= \frac{\mathbf {a}}{\Vert \mathbf {a} \Vert}, \nonumber \end{align} \tag{ 36 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle M_{q \rightarrow b} \frac{{\bf q}}{\Vert \mathbf {q}\Vert}= \mathbf {f}(\mathbf {b}_{1}, \dots, \mathbf {b}_{k}), \nonumber \end{align} \tag{ 37 }$$

where $\newcommand{\re}{{\rm Re}} \mathbf {f}$ depends on the vectors $\newcommand{\re}{{\rm Re}} \mathbf {b}_i$. The main idea of this approach is that if one chooses $\newcommand{\re}{{\rm Re}} \mathbf {f}$ and $M_{q \rightarrow b}$ in an appropriate way, the $M_{r_i}$ may potentially be chosen such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle M_{r_i} \frac{{\bf r}_i}{\Vert {\bf r}_i \Vert} = \frac{\mathbf {b}_i}{\Vert \mathbf {b}_i \Vert}, \quad i=1,\dots,k, \nonumber \end{align} \tag{ 38 }$$

in order for ${\bf p}$ to transform into $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {p}}}$ as desired. One way of 'passing on the norm' of ${\bf q}$ in this way, i.e. choosing $\newcommand{\re}{{\rm Re}} \mathbf {f}$ and $M_{q \rightarrow b}$, is to pass it to each $\newcommand{\re}{{\rm Re}} \mathbf {b}_i$ individually

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{q \rightarrow b}= \sum_{i=1}^{k} \tilde{\alpha}_i (\mathbb{1}+\Pi^i) M_{q \rightarrow b_i}, \nonumber \end{align} \tag{ 39 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\rm with}\ \ M_{q \rightarrow b_i} \frac{{\bf q}}{\Vert {\bf q} \Vert}= \frac{\mathbf {b}_i}{\Vert \mathbf {b}_i \Vert}, \quad i=1,\dots,k, \nonumber \end{align} \tag{ 40 }$$

which implies $\newcommand{\re}{{\rm Re}} \mathbf {f}(\mathbf {b}_{1}, \dots, \mathbf {b}_{k})=\sum\nolimits_{i=1}^{k} \tilde{\alpha}_i \frac{\mathbf {b}_i}{\Vert \mathbf {b}_i \Vert}$, where the $\{\tilde{\alpha}_i\}_{i}$ constitutes a probability distribution. However, this is not always possible as the condition $\newcommand{\re}{{\rm Re}} {\bf q} / \Vert {\bf q} \Vert \succ \mathbf {b}_i / \Vert \mathbf {b}_i \Vert$, necessary for the existence of $M_{q \rightarrow b_{i}}$, cannot be satisfied in general. See appendix A.6 for details. One is therefore forced to pass the norm of ${\bf q}$ in a more clever way. A good candidate is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{q \rightarrow b}= \sum_{i=1}^{k} \tilde{\alpha}_i M_{q \rightarrow 2b_i}, \nonumber \end{align} \tag{ 41 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\rm with}\ \ M_{q \rightarrow 2b_i} \frac{{\bf q}}{\Vert {\bf q} \Vert}= \frac{(\mathbb{1}+\Pi^i)\mathbf {b}_i}{2 \Vert \mathbf {b}_i \Vert}, \quad i=1,\dots,k. \nonumber \end{align} \tag{ 42 }$$

However, this separation of M_q into one part 'passing on the norm' from ${\bf q}$ to $\newcommand{\re}{{\rm Re}} \mathbf {a}$ and another part from ${\bf q}$ to the $\newcommand{\re}{{\rm Re}} \mathbf {b}_i$ is not strictly needed. One may instead only require the existence of a doubly-stochastic matrix M_q such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:generalMq} M_q {\bf q}= \mathbf {a} + \tilde{\mathbf {f}}({\bf r}_1,\mathbf {b}_1,\dots,{\bf r}_k,\mathbf {b}_k), \nonumber \end{align} \tag{ 43 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \tilde{\mathbf {f}}({\bf r}_1,\mathbf {b}_1,\dots,{\bf r}_k,\mathbf {b}_k)=\sum_{i=1}^{k} (\lfloor\frac{2i}{d}\rfloor+1)^{-1} (1-\frac{\Vert {\bf r}_i \Vert}{\Vert \mathbf {b}_i \Vert}) (\mathbb{1}+\Pi^i) \mathbf {b}_i. \nonumber \end{align} \tag{ 44 }$$

Unfortunately, equation (43) is challenging to check, e.g. even confirming whether $\newcommand{\re}{{\rm Re}} \mathbf {a} + \tilde{\mathbf {f}}({\bf r}_1,\mathbf {b}_1,\dots,{\bf r}_k,\mathbf {b}_k)$ is a vector of nonnegative components is already complicated. In summary, the approach requires checking the two following conditions:

• (i)  
  For any pair of vectors $\newcommand{\re}{{\rm Re}} ({\bf r}_i,\mathbf {b}_i)$, there exists a doubly-stochastic matrix $M_{r_i}$ such that

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{\mathbf {b}_i}{\Vert \mathbf {b}_i\Vert} =M_{r_i} \frac{{\bf r}_i}{\Vert {\bf r}_i \Vert},\; i=0,\dots,k, \nonumber \end{align} \tag{ 45 }$$

  while at the same time $\newcommand{\re}{{\rm Re}} \Vert {\bf q} \Vert \geqslant \Vert \mathbf {a} \Vert$, where $\newcommand{\re}{{\rm Re}} \mathbf {a}:={\bf q}(\beta')$.
• (ii)  
  There exists a doubly-stochastic matrix M_q such that

  $$\begin{align} \newcommand{\e}{{\rm e}} mathbf \newcommand{\re}{{\rm Re}} \displaystyle M_q {\bf q}= \mathbf {a} + \sum_{i=1}^{k} (\lfloor\frac{2i}{d}\rfloor+1)^{-1} (1-\frac{\Vert {\bf r}_i \Vert}{\Vert \mathbf {b}_i \Vert}) (\mathbb{1}+\Pi^i) \mathbf {b}_i, \nonumber \end{align} \tag{ 46 }$$

  or the stronger version: $\newcommand{\re}{{\rm Re}} \Vert {\bf r}_i\Vert \leqslant \Vert \mathbf {b}_i\Vert$ for all $i=1,\dots,k$ and there exist doubly-stochastic matrices $M_{q \rightarrow 2 b_i},\; i=1,\dots,k$ such that

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:passnormcondII2} M_{q \rightarrow 2b_i} \frac{{\bf q}}{\Vert {\bf q} \Vert}= \frac{(\mathbb{1}+\Pi^i)\mathbf {b_i}}{2 \Vert \mathbf {b_i} \Vert} ~~\, i=1,\dots,k. \nonumber \end{align} \tag{ 47 }$$

Condition (i) indeed holds, and we thus state it below as the following lemma, the proof of which is presented in appendix A.7.

Lemma 2. For any d-dimensional system with arbitrary Hamiltonian, for any $\beta'< \beta$, and for all $i=0,\dots, d-1$ we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \Vert {\bf q} \Vert \geqslant\, \Vert \mathbf {a} \Vert , \label{lemma:rismallerbi} \nonumber \end{align} \tag{ 48 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{{\bf r}_i}{\Vert {\bf r}_i \Vert} \succ\, \frac{\mathbf {b}_i}{\Vert \mathbf {b}_i \Vert}. \label{lemma:rimajbi} \nonumber \end{align} \tag{ 49 }$$

Given the above lemma and assuming that Condition (ii) also holds, one can prove the existence of STUs for all dimensions by using the following majorisation relations. The majorisation relation (49) from lemma 2 ensures the existence of doubly-stochastic matrices $M_{r_i}$ which satisfy $\newcommand{\re}{{\rm Re}} M_{r_i}\, {\bf r}_i / \Vert {\bf r}_i\Vert =\mathbf {b}_i / \Vert \mathbf {b}_i\Vert$. To reach any thermal state with higher temperature, i.e. $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {p}}} = {\bf p}(\beta')$, we then need the existence of M_q such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\bf p}(\beta^{\prime}) = M_{q}\,{\bf q}+\sum_{i=1}^{k} (\lfloor\frac{2i}{d}\rfloor+1)^{-1}\frac{\Vert {\bf r}_i \Vert}{\Vert \mathbf {b}_i \Vert} (\mathbb{1}+\Pi^i)\, \mathbf {b}_i, \nonumber \end{align} \tag{ 50 }$$

therefore compensating for the required norm of the vectors associated to the subspace $\{\mathcal{H}_{r_i}\}_{i=1}^{d-1}$. This is precisely what is ensured by Condition (ii).

However, since Condition (ii) is in general cumbersome to prove, one can check whether equation (47) holds instead. See appendix A.8 for more details about how to prove the existence of STUs in general in this way. In the particular case of d  =  3, the strong version of Condition (ii) can be reduced to a ('norm passing') requirement that we formulate in the following lemma:

Lemma 3. In d  =  3, for every choice of E₁ and E₂ with $E_{2}\geqslant E_{1}$ and for any $\beta'\leqslant\beta$, the following majorisation relation holds:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{(\mathbb{1} + \Pi)\mathbf {b}}{2 \Vert \mathbf {b} \Vert}:=\frac{\mathbf {c}}{\Vert \mathbf {c} \Vert}. \nonumber \end{align} \tag{ 51 }$$

The proof of lemma 3 is presented in appendix A.9. Together, lemmas 2 and 3 then confirm the existence of STUs for the d  =  3 case.

Now, let us examine the case d  =  4. Consider a bipartite system with equal local Hamiltonians $\newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} H=\sum\nolimits_{i=0}^{3}E_i\ket{i}\!\bra{i}$ in which the energy eigenvalues are ordered in increasing order and E₀  =  0. Following equation (13), the vectorised form of the marginals for the initial uncorrelated thermal state is given by ${\bf p}= {\bf q} +{\bf r}_1+ {\bf r}_2+{\bf r}_3$ with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf q}=\left(\begin{array}{@{}c@{}}p_{00} \nonumber \\ p_{11} \nonumber \\ p_{22} \nonumber \\ p_{33} \end{array}\right), \; {\bf r}_1=\left(\begin{array}{@{}c@{}}p_{01} \nonumber \\ p_{12} \nonumber \\ p_{23} \nonumber \\ p_{30} \end{array}\right),\;{\bf r}_2=\left(\begin{array}{@{}c@{}}p_{02} \nonumber \\ p_{13} \nonumber \\ p_{20} \nonumber \\ p_{31} \end{array}\right),\; {\bf r}_3=\left(\begin{array}{@{}c@{}}p_{03} \nonumber \\ p_{10} \nonumber \\ p_{21} \nonumber \\ p_{32} \end{array}\right)=\Pi {\bf r}_1, \label{eq:d4 vectors} \nonumber \end{align} \tag{ 52 }$$

where ${\bf r}_i$ is a shorthand for ${\bf r}_i(\beta)$. Furthermore, we again denote the vector decomposition of any thermal state with higher temperature $\beta' \leqslant \beta$ as $\newcommand{\re}{{\rm Re}} {\bf p}(\beta')=\mathbf {a}+\sum\nolimits_{i}\mathbf {b}_{i}$ with $\newcommand{\re}{{\rm Re}} \mathbf {a}:={\bf q}(\beta')$ and $\newcommand{\re}{{\rm Re}} \mathbf {b}_i:={\bf r}_i(\beta')$. Unitaries on the LCSs that generate the same marginals, according to equation (21), lead to the transformation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\bf p}\mapsto{\tilde{\mathbf {p}}} = M_{q}\,{\bf q}+(\mathbb {1}+\Pi)\,M_{r_1}{\bf r}_1+\frac{1}{2}(\mathbb {1}+\Pi^2)M_{r_2}\,{\bf r}_2, \label{eq:gen. trans. 4d} \nonumber \end{align} \tag{ 53 }$$

where the $M_{r_i}$ are arbitrary doubly-stochastic matrices. To achieve a thermal marginal with higher temperature, one can transform each ${\bf r}_i$ to $\newcommand{\re}{{\rm Re}} \mathbf {b}_i$ as prescribed by Condition (i) and pass the extra norm of the vector ${\bf q}$ as dictated by the stronger version of Condition (ii). This is possible at least under some restrictions on the energy level spacings $\delta_i:=E_{i+1}-E_i$. Let us phrase this statement more precisely in the following lemma 4, a detailed proof of which is presented in appendix A.9.

Lemma 4. In d  =  4, for every choice of the set $\{E_{i}\}_{i=0}^{3}$ with $E_{i+1}\geqslant E_i$ and $\delta_{i+1} \leqslant \delta_{i}$ and for any $\beta'\leqslant\beta$, the following relations hold for $i=1,2$:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{{\bf q}}{\Vert {\bf q} \Vert} \succ\, \frac{\mathbb{1}+\Pi^i}{2} \frac{\mathbf {b}_i}{\Vert \mathbf {b}_i \Vert}, \label{lemma:qmaj2ri4d} \nonumber \end{align} \tag{ 54 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \Vert {\bf r}_i \Vert \leqslant\, \Vert \mathbf {b}_i \Vert . \label{lemma:rismallerbi4d} \nonumber \end{align} \tag{ 55 }$$

With lemma 4 at hand, we are ready to state the result achieved for the case d  =  4 with this method.

Theorem 2 (Existence of STUs in d  =  4). In any four-dimensional system, for every set of energy eigenvalues $\{E_{i}\}_{i=0}^3$ with $E_{i+1}\geqslant E_i$ and $\delta_{i+1} \leqslant \delta_{i}$, and for any $\beta'\leqslant\beta$, there exists a set of doubly-stochastic matrices $\{M_{r_i}\}_{i=0}^2$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf p}(\beta') = M_{q}\,{\bf q}+(\mathbb {1}+\Pi)\,M_{r_1}{\bf r}_1+\frac{1}{2}(\mathbb {1}+\Pi^2)M_{r_2}\,{\bf r}_2, \nonumber \end{align} \tag{ 56 }$$

which implies the existence of STUs in d  =  4 for symmetric Hamiltonians with decreasing energy gaps.

Proof. The statement follows from lemmas 2 and 4. □

We have thus seen that the approach discussed in this section does, at least partially, generalise to local dimension 4. However, the proof of theorem 2, which exploits the stronger version of Condition (ii) fails for $E_3 \gg 1$. Note that this failure is not an artifact of even dimensions, but rather continues to persist in subsequent higher dimensions. Still, the possibility remains that above results can be generalised by proving the weaker version of Condition (ii).

Our attempts to generalise the approaches of sections 3.1 and 3.2 to higher dimensions have shown that the problem can be recast in terms of different sets of conditions. However, checking these conditions has proven to be increasingly complex with growing dimension, and has thus only provided partial results even for the case d  =  4. We therefore now turn to a third approach which at least provides a complete proof for d  =  4. This approach is centred around the geometric structure generated by doubly-stochastic matrices. More specifically, recall that the Schur–Horn theorem implies that for any $v\in\mathbb{R}^{n}$, the set of vectors obtained by applying the set of unistochastic $n\times n$ matrices to $v$ is a convex polytope given by the convex hull of all permutations of the entries of $v$, see, e.g. [31]. Here, we can apply this idea to the vectors ${\bf q}$ and ${\bf r}_{i}$ and the matrices M_q and $M_{r_{i}}$, respectively, to generate matching diagonal marginals according to equation (21). In other words, the set of all possible symmetric marginal vectors reachable by unitaries that are block-diagonal with respect to the chosen LCS decomposition is the polytope with vertices given by the set of points

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}}^{(i_{0},i_{1},\ldots,i_{k})} = \Pi^{(i_{0})} \mathbf{q} + \sum_{n=1}^{k} (\lfloor\frac{2n}{d}\rfloor+1)^{-1} (\mathbb{1}+ \Pi^{n})\,\Pi^{(i_{n})} \mathbf{r}_{n}, \label{eq:EqualEvoMarginal} \nonumber \end{align} \tag{ 57 }$$

with $i_{j}\in\{1,2,\ldots,d!\}$ for all $j=0,\ldots,k$ with k as in equation (17). Here, $\Pi^{(i)}$ for $i=1,\ldots,d!$ are the possible permutations of d elements. The question about the existence of STUs is then equivalent to asking if the curve defined by the set of points $\{{\bf p}(\beta')|\beta\geqslant\beta'\geqslant0\}$ is enclosed within the polytope corresponding to the convex hull of the (d!)^k+1 points $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {p}}}^{(i_{0},i_{1},\ldots,i_{k})}$.

To answer this question, we then proceed in the following way. First, we note that the d-component vectors $\mathbf{p}=(\,p_{i})$ only have d  −  1 independent entries due to the normalisation condition $\sum\nolimits_{i=0}^{d-1}p_{i}=1$. The problem can thus be reduced to d  −  1 dimensions by choosing coordinates $\{x_{i}\}_{i=0,1,\ldots,d-2}$ with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x_{n} = (n+1)\,p_{n+1}-\sum_{i=0}^{n}p_{i}, \label{eq:new coords} \nonumber \end{align} \tag{ 58 }$$

for $n=0,1,\ldots,d-2$, while the additional last coordinate $x_{d-1}=-\sum\nolimits_{i=0}^{d-1}p_{i}=-1$ is fixed by normalisation and can thus be disregarded. In these coordinates, the point obtained for infinite temperature ($\beta\rightarrow0$) is the origin, $p(\beta\rightarrow0)=(0,0,\ldots,0,-1)$, the point for vanishing temperature is $p(\beta\rightarrow\infty)=(-1,-1,\ldots,-1)$ and all thermal states lie on a continuous curve connecting these points that is strictly confined to negative coordinate regions, $x_i \,\leq\, 0~\forall i$. However, note that there are points corresponding to reachable marginals that have positive values for some of the new coordinates.

For any given initial temperature $1/\beta$ and dimension d, a sufficient set of conditions for the inclusion of the curve defined by the points $\{{\bf p}(\beta')|\,\beta\geqslant\beta'\geqslant0\}$ in the reachable polytope is as follows:

• (I)  
  Inclusion of the vertices

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \mathbf {v}_0&:={\bf p}(\beta)=(x_{0}(\beta),x_{1}(\beta),\ldots,x_{d-3}(\beta),x_{d-2}(\beta),-1) , \nonumber \\ \mathbf {v}_1&:=(0,x_{1}(\beta),\ldots,x_{d-3}(\beta),x_{d-2}(\beta),-1) , \nonumber \\ \mathbf {v}_2&:=(0,0,x_2(\beta)\ldots,x_{d-3}(\beta),x_{d-2}(\beta),-1) , \nonumber \\ &\ \ \vdots \nonumber \\ \mathbf {v}_{d-2}&:=(0,0,\ldots,0,x_{d-2}(\beta),-1) , \nonumber \\ \mathbf {v}_{d-1}&:={\bf p}(\beta\rightarrow 0) = (0,\dots,0,-1) , \label{eq:polyvertgeom} \nonumber \end{align} \tag{ 59 }$$

  in the polytope of achievable marginals.
• (II)  
  Inclusion of all points ${\bf p}(\beta')$ with $\beta'\leqslant \beta$ in the simplex corresponding to convex hull of the set of vertices $\newcommand{\re}{{\rm Re}} \{\mathbf {v}_{i}|\,i=0,\ldots,d-1\}$.

In general, confirming Condition (II) is relatively straightforward, either by proving the positivity of the partial derivatives $\newcommand{\bi}{\boldsymbol} \frac{\partial x_{m}}{\partial x_{n}} =\bigl(\frac{\partial x_{n}}{\partial \beta}\bigr){}^{-1}\frac{\partial x_{m}}{\partial \beta}\geqslant0$ and $\newcommand{\bi}{\boldsymbol} \frac{\partial^{2} x_{m}}{\partial x_{n}^{2}} =\frac{\partial^{2} x_{m}}{\partial \beta^{2}}- \bigl(\frac{\partial x_{m}}{\partial \beta}\bigr)\bigl(\frac{\partial x_{n}}{\partial \beta}\bigr){}^{-1} \frac{\partial^{2} x_{n}}{\partial \beta^{2}}\geqslant0$ (which we show explicitly for d  =  4 in appendix A.10), or by showing that all points ${\bf p}(\beta')$ can be written as a convex combination of the vertices (59), which we prove in general for dimension d in appendix A.10. Let us therefore state Condition (II) as the following lemma:

Lemma 5. In the coordinates defined by equation (58), the curve of thermal states at inverse temperature $\beta'\leqslant\beta$ satisfies Condition (II) for all $\beta$ in all dimensions d.

Besides satisfying Condition (II), it is also easy to see that the vertices $\newcommand{\re}{{\rm Re}} \mathbf {v}_{0}$ and $\newcommand{\re}{{\rm Re}} \mathbf {v}_{d-1}$ from Condition (I), which represent the initial thermal state and maximally mixed state in the new coordinates, are reachable for all dimensions. However, showing the inclusion of the rest of the points in (I) is increasingly difficult, due to the rapidly growing number of possible polytope vertices and the difficulty of visualizing the $(d-1)$-dimensional polytope beyond d  =  4. In the following, we prove that Condition (I) holds (at least) in the particular cases of d  =  3 and d  =  4. See also figure 1 for an illustration in dimension 3 and appendices A.11 and A.12 for further details.

Zoom In Zoom Out Reset image size

Theorem 3 (Geometric approach in d  =  3 and d  =  4). In d  =  3 and d  =  4 systems, for every choice of Hamiltonians and initial inverse temperature $\beta$, the set of thermal states with $\beta' \leqslant \beta$ is contained within the polytope with vertices defined in equation (59), which proves the existence of STUs in the symmetric two-qutrit and two-quqart cases.

Proof. For d  =  3 we have to prove that the point $\newcommand{\re}{{\rm Re}} \mathbf {v}^{(3)}_{1}=(0,x_1(\beta),-1)$ can be reached with transformations of the type ${\bf p}\mapsto\tilde{{\bf p}}=M_{q}{\bf q}+(\mathbb {1} +\Pi)M_r{\bf r}$ for some doubly-stochastic $3\times3$ matrices M_q and M_r, where ${\bf q}$ and ${\bf r}$ are as in equation (29). This is indeed the case, e.g. for $M_r=\mathbb {1}$ being the identity matrix and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{q} = \left(\begin{array}{@{}ccc@{}} m & 1-m & 0 \nonumber \\ 1-m & m & 0 \nonumber \\ 0 & 0 & 1 \end{array}\right), \nonumber \end{align} \tag{ 60 }$$

with $m=1-1/[2(\,p_{0}+p_{1})]$, which is a doubly-stochastic matrix since $p_{0}+p_{1}\geqslant 1/2$ for d  =  3. Note that, geometrically, this choice of M_q and M_r represents a convex combination of the points ${\bf p}(\beta)$ and $\tilde{{\bf p}}^{(6,1)}$ shown in figure 1. For this transformation, we have $\newcommand{\re}{{\rm Re}} \tilde{{\bf p}}=\sum\nolimits_{i=0}^{2}\tilde{p}_{i}\mathbf {e}_{i}$, where the components with respect to the original basis $\newcommand{\re}{{\rm Re}} \{\mathbf {e}_{j}\}_{j=0,1,2}$ are $\tilde{p}_{0}=\tilde{p}_{1}=(\,p_{0}+p_{1})/2$ and $\tilde{p}_{2}=p_{2}$, which means that the components of $\tilde{{\bf p}}$ with respect to to the new coordinates of equation (58) are $\newcommand{\re}{{\rm Re}} \tilde{{\bf p}}=(0,2-3(\,p_{0}+p_{1}),-1){}^{T}=(0,x_{1}(\beta),-1){}^{T}=\mathbf {v}^{(3)}_{1}$. This concludes the proof for d  =  3.

For d  =  4, we have to show that the two points $\newcommand{\re}{{\rm Re}} \mathbf {v}^{(4)}_{1}=(0,x_1(\beta),x_2(\beta),-1)$ and $\newcommand{\re}{{\rm Re}} \mathbf {v}^{(4)}_{2}=(0,0,x_2(\beta),-1)$ can be reached with transformations of the type ${\bf p}\mapsto\tilde{{\bf p}}=M_q{\bf q}+(\mathbb {1}\,+$ $\Pi)M_{r_1}{\bf r}_1+\frac{1}{2} (\mathbb {1}+\Pi^2)M_{r_2}{\bf r}_2$ for some doubly-stochastic $4\times4$ matrices M_q, $M_{r_{1}}$, and $M_{r_{2}}$, where ${\bf q}$ and ${\bf r}_{1}$ are as in equation (52).

The point $\newcommand{\re}{{\rm Re}} \mathbf {v}^{(4)}_{1}$ can be reached with the equivalent of the transformation used to reach $\newcommand{\re}{{\rm Re}} \mathbf {v}^{(3)}_{1}$ above, that is, using $M_{r_{1}}=M_{r_{2}}=\mathbb {1}$ and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{q} = \left(\begin{array}{@{}cccc@{}} m & 1-m & 0 & 0 \nonumber \\ 1-m & m & 0 & 0 \nonumber \\ 0 & 0 & 1 & 0 \nonumber \\ 0 & 0 & 0 & 1 \end{array}\right) , \nonumber \end{align} \tag{ 61 }$$

with $m=1-1/[2(\,p_{0}+p_{1})]$, which is again doubly-stochastic since $p_{0}+p_{1}\geqslant 1/2$ also holds for d  =  4.

To prove that $\newcommand{\re}{{\rm Re}} \mathbf {v}^{(4)}_{2}$ can be reached, one can see that it can be obtained as a convex combination of (at most) 5 of the vertices in equation (57), as we show in detail in appendix A.12. □

For higher dimensions, as we already mentioned, the problem becomes more and more complex, but one can try to build a recursive approach based on the above lower dimensional proofs and borrowing some ideas from the 'passing on the norm' approach. In particular, we outline a possible route for such an approach for the case d  =  5 in appendix A.13, where we can show the existence of STUs for d  =  5 for a subset of all possible Hamiltonians.

In this appendix, we show that if STUs exist in general, i.e. in particular for the desired temperature, they provide an upper bound for the amount of correlation that can be achieved unitarily. Using the same notation as in the main text, we are interested in solving the problem

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \max_U S(\tilde{\varrho}\SA) + S(\tilde{\varrho}\SB) \quad ~{\rm s.t.}\ \ \tr(\tilde{\varrho}\SA H\SA) + \tr(\tilde{\varrho}\SB H\SB) \leqslant c. \nonumber \end{align} \tag{ A.1 }$$

This can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle & \max_U S(\tilde{\varrho}\SA \otimes \tilde{\varrho}\SB)\nonumber \\ &{\rm s.t.}\ \ \tr(\tilde{\varrho}\SA \otimes \tilde{\varrho}\SB (H\SA \otimes \mathbb{1} + \mathbb{1} \otimes H\SB)) \leqslant c. \nonumber \end{align} \tag{ A.2 }$$

According to Jaynes' principle, the maximum is obtained for

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tau\SA(\bar{\beta}) \otimes \tau\SB(\bar{\beta}), \nonumber \end{align} \tag{ A.3 }$$

where $\bar{\beta}$ is chosen such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tr (\tau\SA(\bar{\beta}) \otimes \tau\SB(\bar{\beta}) (H\SA \otimes \mathbb{1} + \mathbb{1} \otimes H\SB)) = c. \nonumber \end{align} \tag{ A.4 }$$

This solution can be found using Lagrange multipliers by considering the $n\times n$ matrix $\newcommand{\SA}{{_{A}}} \newcommand{\SAB}{{_{AB}}} \varrho\SAB$ as an vector with n² components. We therefore have as desired

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle &\max_U S(\tilde{\varrho}\SA) + S(\tilde{\varrho}\SB)\ {\rm s.t.}\ \tr (\tilde{\varrho}\SA H\SA) + \tr (\tilde{\varrho}\SB H\SB) \leqslant c. \nonumber \\ &\leqslant \, S(\tau\SA(\bar{\beta}))+S(\tau\SB(\bar{\beta})), \nonumber \end{align} \tag{ A.5 }$$

with $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \tr (\tau\SA(\bar{\beta}) H\SA)+ \tr (\tau\SB(\bar{\beta}) H\SB) = c$.

Here, we discuss and solve the problem of maximising the correlations under an energy constraint in the asymmetric case with an initial pure state. That is, we want to solve

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \max_{\tilde{\varrho}\SA, \tilde{\varrho}\SB} S(\tilde{\varrho}\SA)+ S(\tilde{\varrho}\SB), \nonumber \end{align} \tag{ A.6 }$$

subject to the constraint $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \tr(\tilde{\varrho}\SA H\SA) + \tr(\tilde{\varrho}\SB H\SB) \leqslant c$, where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\langle #1 |} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tilde{\varrho}\SA=\sum_{i=0}^{d-1} p_{i} \ket{\varphi^{A}_{i}}\!\bra{\varphi^{A}_{i}},\ \ \ \tilde{\varrho}\SB=\sum_{i=0}^{d-1} p_{i} \ket{\varphi^{B}_{i}}\!\bra{\varphi^{B}_{i}}, \label{sub state A} \nonumber \end{align} \tag{ A.7 }$$

with $p_i \geqslant 0$ for $i=0, \dots, d-1$, $\sum\nolimits_{i=0}^{d-1} p_{i}=1$, $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} d=\min\{d\SA,d\SB\}$, and where $\newcommand{\SA}{{_{A}}} \newcommand{\ket}[1]{| #1\rangle} \{\ket{\varphi^{A}_{i}}\}_{i=0}^{d\SA-1}$ and $\newcommand{\SB}{{_{B}}} \newcommand{\ket}[1]{| #1\rangle} \{\ket{\varphi^{B}_{i}}\}_{i=0}^{d\SB-1}$ are orthonormal bases of $\newcommand{\SA}{{_{A}}} \mathcal{H}\SA$ and $\newcommand{\SB}{{_{B}}} \mathcal{H}\SB$, respectively. Without loss of generality we then assume $\newcommand{\SA}{{_{A}}} d=d\SA$. Note that $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} S(\tilde{\varrho}\SA)=S(\tilde{\varrho}\SB)$ and $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\tilde{\varrho}\SB = U \tilde{\varrho}\SA U^{\dagger}$, where we write $\newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} U= \sum\nolimits_{i=0}^{d-1} \ket{\varphi^{B}_{i}}\!\bra{\varphi^{A}_{i}}$ (in a slight abuse of notation) such that the problem may be rewritten as

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \label{equ:maxS} \max_{\varrho, U} S(\varrho) \quad {\rm s.t.}\ \ \ \tr[\varrho (H\SA+U^{\dagger} H\SB U)]\,\leqslant\, c, \nonumber \end{align} \tag{ A.8 }$$

where we have dropped the tilde and subscript A on $\varrho$ for ease of notation. To solve this problem we consider the converse problem

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \label{equ:minEn} \min_{\varrho, U} \tr [\varrho (H\SA+U^{\dagger} H\SB U)] \quad {\rm s.t.}\ \ \ S(\varrho)= \kappa, \nonumber \end{align} \tag{ A.9 }$$

and show that (at least a family of) optimal points of (A.9) are optimal points of (A.8). To simplify the notation further let us write

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\langle #1 |} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\SA}{{_{A}}} \displaystyle H\SA = \sum_{i=0}^{d-1} E_{i}^{A} \ket{i}\!\bra{i}\SA \nonumber \end{align} \tag{ A.10 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\langle #1 |} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\SB}{{_{B}}} \displaystyle H\SB = \sum_{i=0}^{d_B-1} E_{i}^{B} \ket{i}\!\bra{i}\SB, \nonumber \end{align} \tag{ A.11 }$$

with $E_{i}^{A} \leqslant E_{i+1}^{A}$ and $E_{i}^{B} \leqslant E_{i+1}^{B}$.

Proposition A.1. The pair $(\varrho_{\mathrm{opt}}(\kappa), U_{\mathrm{opt}})$, given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\langle #1 |} \newcommand{\ket}[1]{| #1\rangle} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle U_{\mathrm{opt}} &:= \sum_{i=0}^{d-1} \ket{i}\SB \bra{i}\SA, \nonumber \end{align} \tag{ A.12 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varrho_{\mathrm{opt}}(\beta(\kappa))&:= \frac{\tilde{\Pi} e^{-\beta(\kappa) \tilde{H}}\tilde{\Pi}}{\tilde{Z}},\nonumber \end{align} \tag{ A.13 }$$

where $\newcommand{\suptiny}[3]{{^{\hspace{#1 pt}\protect\raisebox{#2 pt}{\tiny{$ #3$}}}}} \tilde{Z}=\mathrm{Tr}(\tilde{\Pi} e^{-\beta(\kappa) \tilde{H}})$ and , is a solution of the minimisation in (A.9).

Proof. Denoting the spectrum of $\varrho$ as $\lambda_{\varrho}=(\lambda_0,\dots,\lambda_{d-1})$, we first show that $(\varrho_{\mathrm{opt}}(\beta(\kappa)), U_{\mathrm{opt}})$ is a solution of the following minimisation problem

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \min_{\lambda} \left(\min_{\varrho\ {\rm s.t.}\ \lambda_{\varrho}=\lambda} \tr (\varrho H\SA) + \min_{U, \varrho\ {\rm s.t.}\ \lambda_{\varrho}=\lambda} \tr (U \varrho U^{\dagger} H\SB) \right), \label{equ:minEdecomp} \nonumber \end{align} \tag{ A.14 }$$

subject to the constraint $H(\lambda)=\kappa$, where H($\lambda$) denotes the Shannon entropy of probability distribution $\lambda$ with length d. Since all density matrices of a given spectrum are unitarily related we have

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SA}{{_{A}}} \displaystyle \min_{\varrho\ {\rm s.t.}\ \lambda_{\varrho}=\lambda} \tr (\varrho H\SA)= \min_{V} \tr(V {\rm diag}\{\lambda_\varrho\} V^{\dagger} H\SA), \nonumber \end{align} \tag{ A.15 }$$

where the minimisation on the right-hand side is over all unitaries $V$. The passive state with spectrum $\lambda$ is well-known to solve this minimisation. Adopting the notation $v^{\downarrow}=(v^{\downarrow}_{i})$ to denote the vector obtained by arranging the components of the vector $v=(v_1,\dots,v_n)$ in decreasing order, i.e. such that $v_1^{\downarrow} \geqslant v_2^{\downarrow} \geqslant \dots \geqslant v_n^{\downarrow}$, we thus have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \newcommand{\SA}{{_{A}}} \displaystyle \min_{\varrho\ {\rm s.t.}\ \lambda_{\varrho}=\lambda} \tr (\varrho H\SA)= \sum_{i=0}^{d-1} \lambda_{i}^{\downarrow} E_{i}^{A}. \nonumber \end{align} \tag{ A.16 }$$

For the second minimisation problem in (A.14), note that since $\newcommand{\SA}{{_{A}}} \newcommand{\ket}[1]{| #1\rangle} \{\ket{\varphi_{i}^{A}}\}_{i=0}^{d\SA-1}$ and $\newcommand{\SB}{{_{B}}} \newcommand{\ket}[1]{| #1\rangle} \{\ket{\varphi_{i}^{B}}\}_{i=0}^{d\SB-1}$ are orthonormal bases, the matrix representation of U can be extended to a unitary $\newcommand{\SB}{{_{B}}} d\SB \times d\SB$ matrix $V$. Similarly, the matrix representation of $\varrho$ can be extended to a positive $\newcommand{\SB}{{_{B}}} d\SB\times d\SB$ matrix $\bar{\varrho}$ by padding it with zeroes. With this we have

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \displaystyle \min_{U, \varrho\ {\rm s.t.}\ \lambda_{\varrho}=\lambda} \tr (U \varrho U^{\dagger} H\SB)\geqslant \min_{V, \bar{\varrho}\ {\rm s.t.}\ \lambda_{\bar{\varrho}}=(\lambda,0,\dots,0)} \tr (V \bar{\varrho} V^{\dagger} H\SB). \end{align} \tag{ A.17 }$$

Again the passive state with spectrum $(\lambda,0,\dots,0)$ solves the right-hand side of inequality (A.17) and we obtain

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \displaystyle \min_{V, \bar{\varrho}\ {\rm s.t.}\ \lambda_{\bar{\varrho}}=(\lambda,0,\dots,0)} \tr (V \bar{\varrho} V^{\dagger} H\SB) = \sum_{i=0}^{d-1} \lambda_{i}^{\downarrow} E_{i}^{B}. \nonumber \end{align} \tag{ A.18 }$$

This solution is in fact also an attainable solution of the left-hand side of inequality (A.17). Hence

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \displaystyle \min_{U, \varrho\ {\rm s.t.}\ \lambda_{\varrho}=\lambda} \tr (U \varrho U^{\dagger} H\SB)= \sum_{i=0}^{d-1} \lambda_{i}^{\downarrow} E_{i}^{B}. \nonumber \end{align} \tag{ A.19 }$$

We have therefore reduced the minimisation problem of (A.14) to solving

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \min_{\lambda} \sum_{i=0}^{d-1} \lambda_{i}^{\downarrow} (E_{i}^{A}+ E_{i}^{B}),\ {\rm s.t.}\ H(\lambda)=\kappa,\ \sum_{i} \lambda_{i}^{\downarrow}=1. \nonumber \end{align} \tag{ A.20 }$$

This problem, in turn, can be solved by means of Lagrange multipliers, which yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lambda_{i}^{\downarrow}= \frac{{\rm e}^{-\beta (E_{i}^{A}+ E_{i}^{B})}}{\sum\nolimits_{i=0}^{d-1} {\rm e}^{-\beta (E_{i}^{A}+ E_{i}^{B})}}, \nonumber \end{align} \tag{ A.21 }$$

which is precisely what is delivered by the solution $(\varrho_{\mathrm{opt}}(\beta(\kappa)),U_{\mathrm{opt}})$.

We can further check that for every $\kappa \in [0, \ln(d)]$ there exists a unique $\beta \in [0,\infty]$ such that $H(\lambda_{\varrho_{\mathrm{opt}}})=S(\varrho_{\mathrm{opt}}(\beta))=\kappa$, i.e. that the notation $\beta(\kappa)$ is well defined and can be understood as a function. This can be seen from

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S(\varrho_{\mathrm{opt}}(0)) = \ln(d), \nonumber \end{align} \tag{ A.22 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S(\varrho_{\mathrm{opt}}(\infty)) = 0, \nonumber \end{align} \tag{ A.23 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{\rm d}{{\rm d}\beta} S(\varrho_{\mathrm{opt}}(\beta)) = \beta (\tr(\varrho \tilde{H})^{2}- \tr(\varrho \tilde{H}^2)) < 0 , \; \forall \beta \in (0, \infty). \nonumber \end{align} \tag{ A.24 }$$

Strictly speaking, the last line is not valid when $\tilde{H} \propto \mathbb{1}$, but in that case one can straightforwardly check that $(U_{\mathrm{opt}}, \frac{\mathbb{1}}{d})$ solves our original problem (A.8) for any allowed c. We hence (tacitly) discard it from the start.

Having established this fact about $(\varrho_{\mathrm{opt}}(\beta(\kappa)),U_{\mathrm{opt}})$, let $(\varrho,U)$ be such that $S(\varrho)=\kappa$. Then $\varrho$ has some spectrum $\{\mu_0,\dots, \mu_{d-1}\}$. From the above we thus have

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle &\tr (\varrho H\SA) + \tr (U \varrho U^{\dagger} H\SB) \,\geqslant\, \sum_{i=0}^{d-1} \mu_i^{\downarrow} E_{i}^{A}+\sum_{i=0}^{d-1} \mu_i^{\downarrow} E_{i}^{B}\nonumber \\ &\geqslant\, \sum_{i=0}^{d-1} \frac{{\rm e}^{-\beta(\kappa) (E_{i}^{A}+ E_{i}^{B})}}{\sum\nolimits_{i=0}^{d-1} {\rm e}^{-\beta(\kappa)(E_{i}^{A}+E_{i}^{B})}} (E_{i}^{A}+E_{i}^{B}) \nonumber \\ & = \tr(\varrho_{\mathrm{opt}}(\beta(\kappa)) H\SA) + \tr (U_{\mathrm{opt}} \varrho_{\mathrm{opt}}(\beta(\kappa)) U_{\mathrm{opt}}^{\dagger} H\SB), \nonumber \end{align} \tag{ A.25 }$$

which proves our claim, because $S(\varrho_{\mathrm{opt}}(\beta(\kappa)))=\kappa$. □

Now let us define the function

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle f(\kappa)= \tr(\varrho_{\mathrm{opt}}(\beta(\kappa)) H\SA) + \tr (U_{\mathrm{opt}} \varrho_{\mathrm{opt}}(\beta(\kappa)) U_{\mathrm{opt}}^{\dagger} H\SB). \nonumber \end{align} \tag{ A.26 }$$

We can then establish the following proposition.

Proposition A.2. If $f(\kappa_1) < f(\kappa_2)$ then $\kappa_1 < \kappa_2$.

The above is saying that if f  has an inverse then that inverse is strictly monotonically increasing. This is indeed how the proof proceeds.

Proof. First, note that we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\rm d}{{\rm d}\kappa} f(\kappa) = \frac{\rm d}{{\rm d}\beta} f(\beta) \frac{{\rm d}\beta}{{\rm d}\kappa} = \frac{\rm d}{{\rm d}\beta} f(\beta) \left(\frac{\rm d}{{\rm d}\beta} S(\varrho_{\mathrm{opt}}(\beta))\right)^{-1}. \nonumber \end{align} \tag{ A.27 }$$

Further, we have already seen that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{\rm d}{{\rm d}\beta} S(\varrho_{\mathrm{opt}}(\beta)) = \beta (\tr(\varrho \tilde{H})^{2}- \tr(\varrho \tilde{H}^2)) < 0. \nonumber \end{align} \tag{ A.28 }$$

Similarly one calculates

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{\rm d}{{\rm d}\beta} f(\beta) = (\tr(\varrho \tilde{H})^{2}- \tr(\varrho \tilde{H}^2)) < 0. \nonumber \end{align} \tag{ A.29 }$$

So for $\kappa \in (0, \ln(d))$ we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\rm d}{{\rm d}\kappa} f(\kappa) = \frac{1}{\beta} > 0. \nonumber \end{align} \tag{ A.30 }$$

This proves that f  is strictly monotonically increasing. It therefore has an inverse f ⁻¹ that is also strictly monotonically increasing which proves our claim. □

We are now ready to solve the maximisation in (A.8).

Proposition A.3. There is a (unique) $\beta$ such that $(\varrho_{\mathrm{opt}}(\beta), U_{\mathrm{opt}})$ is a solution of (A.8).

Proof. First, if $c \geqslant \frac{1}{d} \sum\nolimits_{i=0}^{d-1} E_{i}^{A}+E_{i}^{B}$, then $(\varrho_{\mathrm{opt}}(0),U_{\mathrm{opt}})$ is the sought after solution, because

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H((\frac{1}{d}, \dots, \frac{1}{d})) = S(\varrho_{\mathrm{opt}}(0)) \,\geqslant\, S(\varrho)\ \forall\, \varrho, \nonumber \end{align} \tag{ A.31 }$$

and $\newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \tr(\varrho_{\mathrm{opt}}(0) \tilde{H}) \leqslant c$. Second, if $c \leqslant \frac{1}{d} \sum\nolimits_{i=0}^{d-1} E_{i}^{A}+E_{i}^{B}$ then

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\rm d}{{\rm d}\beta} f(\beta) < 0 \nonumber \end{align} \tag{ A.32 }$$

ensures that there exists a (unique) $\beta$, say $\bar{\beta}$, such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm Tr}} \newcommand{\re}{{\rm Re}} \displaystyle f(\bar{\beta}) = \tr (\varrho_{\mathrm{opt}} (\bar{\beta}) \tilde{H}) = c. \nonumber \end{align} \tag{ A.33 }$$

Let us now prove that $(\varrho_{\mathrm{opt}}(\bar{\beta}), U_{\mathrm{opt}})$ is the desired solution in this case. Consider a pair $(\varrho_1,U_1)$ such that

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tr (\varrho_1 H\SA)+ \tr (U_1 \varrho_1 U_1^{\dagger} H\SB) \,\leqslant\, c\, = \tr (\varrho_{\mathrm{opt}} (\bar{\beta}) \tilde{H}). \nonumber \end{align} \tag{ A.34 }$$

Let us further consider $\varrho_{\mathrm{opt}}(\beta(\kappa))$, the solution of (A.9) and choose $\kappa=S(\varrho_1)$. Then it holds that

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\tr}{{\rm Tr}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \displaystyle \tr(\varrho_{\mathrm{opt}}(\beta(\kappa)) \tilde{H}) &\leqslant\, \tr (\varrho_1 H\SA)+ \tr (U_1 \varrho U_1^{\dagger} H\SB) \nonumber \\ &\leqslant\, \tr (\varrho_{\mathrm{opt}} (\bar{\beta}) \tilde{H}), \nonumber \end{align} \tag{ A.35 }$$

which, using proposition A.2, implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S(\varrho_{\mathrm{opt}}(\beta(\kappa))) \,\leqslant\, S(\varrho_{\mathrm{opt}} (\bar{\beta})). \nonumber \end{align} \tag{ A.36 }$$

But $S(\varrho_{\mathrm{opt}}(\beta(\kappa))) = S(\varrho_1)$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S(\varrho_1) \,\leqslant\, S(\varrho_{\mathrm{opt}} (\bar{\beta})). \nonumber \end{align} \tag{ A.37 }$$

□

Here we discuss a protocol that can be applied to n copies of the initial state, i.e. $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \tau\SA (\beta){}^{\otimes n} \otimes \tau\SB (\beta){}^{\otimes n}$, to increase the temperature of the marginals in cases where STUs for single copies exist for small temperature differences, but are not attainable for larger temperature differences. For n copies, the protocol consists of n consecutive steps. In each step, a fixed STU achieving some finite temperature increase is applied to different LCSs that correspond to particular pairs of subsystems, such that for $i,j=1,\dots,n$, each subsystem A_i interacts with one and only one subsystem B_j, and no subsystem interacts again with a subsystem it has previously interacted with. This ensures that all marginals are left diagonal and thermal after each step, and leads to a step-wise increase in the temperature of the marginals. In particular, in the j th step the unitary is applied to the subsystems A_i and $B_{i + j -1 \, {\rm mod}(n)}$. This construction guarantees that the tensor product structure of thermal states is preserved for each pair, i.e. $\newcommand{\SA}{{_{A}}} \newcommand{\SB}{{_{B}}} \tau\SA (\beta_j) \otimes \tau\SB (\beta_j)$ where $\beta_j$ denotes the effective local inverse temperature after the j th step. Furthermore, this means that in each step the STU can be applied, as the marginals are in a thermal state with the same temperature and moreover, product to each other. While this protocol ensures that the desired structures are preserved, the necessary conditions for the protocol to achieve arbitrary final temperatures are still part of ongoing research. In general the question remains, whether it is possible to reach any arbitrary temperature difference, $\beta'< \beta$, within finitely many steps. This discussion contrasts the proof of the existence of STUs in the asymptotic case, discussed in section 2.2.3.

To illustrate the protocol, let us illustrate the case of n  =  4 copies here. In the following, each dot reprsents a subsystem and the application of STUs is denoted by lines connnecting pairs of subsystems. Initially one is confronted with the following situation.

In the first step, we connect subsystems corresponding to the same copy, i.e. A_i and B_i.

As specified, we are now in the situation where we have slightly decreased the inverse temperature of the marginals to $\tau(\beta_1)$. We now apply a unitary to the subsystems A_i and $B_{i +1 \, {\rm mod}(n)}$, further decreasing the inverse temperature of the marginals.

In the third step, the subsystems A_i and $B_{i + 2 \, {\rm mod}(n)}$ are connected.

Lastly, we connect A_i and $B_{i + 3 \, {\rm mod}(n)}$ as prescribed.

We have now used the four copies of the initial state, to increase the correlations between the two sides step by step, while retaining the product structure of the thermal marginals, arriving at the final reduced states $\rho_{A/B}=\tau(\beta_4){}^{\otimes 4}$ for the subsystems A and B.

In this section, we will first show that for any $3\times 3$ doubly-stochastic matrix M, there exists another $3\times 3$ doubly-stochastic matrix $\tilde{M}$ such that $M (\mathbb{1}+\Pi)=(\mathbb{1}+\Pi) \tilde{M}$, a statement we called lemma 1 in the main text.

Proof of lemma 1. In general, all doubly-stochastic matrices can be written as convex combinations of permutation matrices as the following

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M=\sum_{i=1}^{d!} \alpha_i\, \Pi^{(i)}, \label{eq:gen. doubly stoch} \nonumber \end{align} \tag{ A.38 }$$

where $\Pi^{(i)}$ indicate permutation matrices in dimension d, and $\sum_i \alpha_i =1$. Since we need to show that the statement is true for any doubly-stochastic matrix M, it is sufficient to prove that the statement is true for M being a permutation matrix. In the following we will just focus on dimension 3. In this dimension there are $3!=6$ permutation matrices where we collect the cyclic permutations $\Pi^{(1)}_{\mathrm{C}}=\mathbb{1}$, $\Pi^{(2)}_{\mathrm{C}}=\Pi$, $\Pi^{(3)}_{\mathrm{C}}=\Pi^{2}$, and anticyclic permutations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Pi^{(1)}_{\mathrm{AC}} = \left(\begin{array}{@{}ccc@{}} 1 & 0 & 0 \nonumber \\ 0 & 0 & 1 \nonumber \\ 0 & 1 & 0 \end{array}\right),\ \Pi^{(2)}_{\mathrm{AC}} = \left(\begin{array}{@{}ccc@{}} 0 & 1 & 0 \nonumber \\ 1 & 0 & 0 \nonumber \\ 0 & 0 & 1 \end{array}\right),\ \Pi^{(3)}_{\mathrm{AC}} = \left(\begin{array}{@{}ccc@{}} 0 & 0 & 1 \nonumber \\ 0 & 1 & 0 \nonumber \\ 1 & 0 & 0 \end{array}\right),\ \nonumber \end{align} \tag{ A.39 }$$

in a particular order, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Pi^{(1)} & = \Pi^{(1)}_{\mathrm{C}},\ \Pi^{(2)} = \Pi^{(1)}_{\mathrm{AC}},\ \Pi^{(3)} = \Pi^{(3)}_{\mathrm{C}},\nonumber \\ \Pi^{(4)} & = \Pi^{(3)}_{\mathrm{AC}},\ \Pi^{(5)} = \Pi^{(2)}_{\mathrm{C}},\ \Pi^{(6)} = \Pi^{(2)}_{\mathrm{AC}} \label{permutationdim3dim} \nonumber \end{align} \tag{ A.40 }$$

where $\Pi^{(1)}$, $\Pi^{(5)}$, and $\Pi^{(3)}$ trivially commute with $\Pi^{(5)}=\Pi$. It is further straightforward to show that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Pi^{(6)} (\mathbb{1}+\Pi^{(5)}) &= (\mathbb{1}+\Pi^{(5)}) \Pi^{(2)}\nonumber \\ \Pi^{(2)} (\mathbb{1}+\Pi^{(5)}) &= (\mathbb{1}+\Pi^{(5)}) \Pi^{(4)} \label{eq:comm. relation 3}\nonumber \\ \Pi^{(4)} (\mathbb{1}+\Pi^{(5)}) &= (\mathbb{1}+\Pi^{(5)}) \Pi^{(6)}\nonumber . \nonumber \end{align} \tag{ A.41 }$$

For any $3\times 3$ doubly-stochastic matrix M given in the form of equation (A.38), by using equation (A.41), one may thus find a $3\times 3$ doubly-stochastic matrix $\tilde{M}$ of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{M} &=\alpha_1\, \Pi^{(1)}+\alpha_3\, \Pi^{(3)}+\alpha_5\, \Pi^{(5)}\nonumber \\ &\quad +\alpha_2\, \Pi^{(4)}+\alpha_4\, \Pi^{(6)}+\alpha_6\, \Pi^{(2)}, \nonumber \end{align} \tag{ A.42 }$$

which satisfies the required condition. □

One of the approaches to investigating the existence of the unitaries mentioned in question 1 is to generalise the 'majorised marginals approach' of section 3.1. This approach appears to be promising because of its simplicity and utility. For such a generalisation to higher dimensions to be successful, equation (57) would demand to prove the following claim.

Claim A.1. For any $d \times d$ doubly-stochastic matrix M, with $d\geqslant 4$, there exists a doubly-stochastic matrix $\tilde{M}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M (\mathbb{1}+\Pi^i)=(\mathbb{1}+\Pi^i) \tilde{M}_i~ \forall i \in \{1, \dots , \lfloor \frac{d}{2} \rfloor\}. \label{eq:gen comm M} \nonumber \end{align} \tag{ A.43 }$$

Note that by permuting the chosen basis vectors in a particular manner, $(\mathbb{1}+\Pi^i)$ can always be written in the form of $(\mathbb{1}+\Pi)$. Additionally, we know that such a transformation of the basis transforms any doubly-stochastic matrix to another doubly-stochastic matrix. Since equation (A.43) should be true for any doubly-stochastic matrix M, these two facts imply that equation (A.43) can be reduced to showing that for any doubly-stochastic $d\times d$ matrix M there exists a doubly-stochastic matrix $\tilde{M}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M (\mathbb{1}+\Pi) = (\mathbb{1}+\Pi) \tilde{M}. \label{eq:swapping relation} \nonumber \end{align} \tag{ A.44 }$$

Using a counterexample, we show that claim A.1 does not hold in general for dimensions $d\geqslant 4$. In particular, we construct a counterexample in dimension 4.

A.5.1. Counterexample to claim A.1.

Let

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M=\left(\begin{array}{@{}cccc@{}} 1&0&0&0 \nonumber \\ 0&0&0&1 \nonumber \\ 0&1&0&0 \nonumber \\ 0&0&1&0. \end{array}\right){.} \label{m-temp} \nonumber \end{align} \tag{ A.45 }$$

As we will now show, there exists no doubly-stochastic $\tilde{M}$ such that $M (\mathbb{1} + \Pi)= (\mathbb{1}+\Pi) \tilde{M}$. To do so, we first calculate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M(\mathbb{1}+\Pi)=\left(\begin{array}{@{}cccc@{}}1&0&0&1 \nonumber \\ 0&0&1&1 \nonumber \\ 1&1&0&0 \nonumber \\ 0&1&1&0 \end{array}\right). \nonumber \end{align} \tag{ A.46 }$$

We are then interested in determining whether this can be equal to $(\mathbb{1}+\Pi) \tilde{M}$ for some $\tilde{M}$, with elements $m_{ij} \geqslant 0$, for $i,j =0,\dots,3$ such that each column and row sums to 1. From the first row of $(\mathbb{1}+\Pi)\tilde{M}$, we obtain $m_{{00}}+m_{{30}}=1$, which implies $m_{{10}}=m_{{20}}=0$, $m_{{01}}=m_{{31}}=0$, and $m_{{02}}=m_{{32}}=0$, and further we obtain $m_{{03}}+m_{{33}}=1$, which implies $m_{{13}}=m_{{23}}=0$. Thus, $\tilde{M}$ must be of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{M}=\left(\begin{array}{@{}cccc@{}} m_{00}&0&0&m_{03} \nonumber \\ 0&m_{11}&m_{12}&0 \nonumber \\ 0&m_{21}&m_{22}&0 \nonumber \\ m_{30}&0&0&m_{33}\end{array}\right). \nonumber \end{align} \tag{ A.47 }$$

Then, from the second row of $(\mathbb{1}+\Pi)\tilde{M}$ we obtain m₀₀  =  0, implying m₃₀  =  1, as well as m₁₁  =  0, which implies m₂₁  =  1. Moreover, we have m₁₂  =  1, which means m₂₂  =  0, while m₀₃  =  1 suggests m₃₃  =  0. We hence have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{M}=\left(\begin{array}{@{}cccc@{}} 0&0&0&1 \nonumber \\ 0&0&1&0 \nonumber \\ 0&1&0&0 \nonumber \\ 1&0&0&0\end{array}\right). \nonumber \end{align} \tag{ A.48 }$$

But since $((\mathbb{1}+\Pi) \tilde{M})_{31}=0 \neq (M (\mathbb{1}+\Pi))_{31}=0$ we arrive at a contradiction.

Since lemma 1 does not hold in general, one can try to relax it to a weaker statement which is still suitable for our purposes. One way to do so is to demand that an equivalent of equation (A.44) holds only when applied to (certain) vectors, in the spirit of the observation that although the set of doubly-stochastic matrices does not coincide with that of unistochastic ones, given a vector $\newcommand{\re}{{\rm Re}} \mathbf {v}$ with nonnegative components and a doubly-stochastic matrix M, there is yet a unistochastic M_U such that $\newcommand{\re}{{\rm Re}} M \mathbf {v} = M_U \mathbf {v}$. Thus, we investigate whether the following statement is true or not:

Claim A.2. Given a vector $\newcommand{\re}{{\rm Re}} \mathbf {v}$ of nonnegative components and a doubly-stochastic matrix M, there exists a doubly-stochastic matrix $\tilde{M}_{v}$, which may depend on $\newcommand{\re}{{\rm Re}} \mathbf {v}$, such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle M (\mathbb{1}+\Pi) \mathbf {v}=(\mathbb{1}+\Pi) \tilde{M}_{v} \mathbf {v}. \nonumber \end{align} \tag{ A.49 }$$

The relaxation being clearly that now $\tilde{M}$ is allowed to depend on $\newcommand{\re}{{\rm Re}} \mathbf {v}$. Unfortunately, the previous counterexample carries over to claim A.2, as we shall see.

A.5.2. Counterexample to claim A.2.

Let us consider M as in equation (A.45) and let $\newcommand{\re}{{\rm Re}} \mathbf {v}=(1,0,0,0){}^{{T}}$. Then we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:Monv} M (\mathbb{1}+\Pi) \mathbf {v} = (1,\,0,\,1,\,0)^{T}. \nonumber \end{align} \tag{ A.50 }$$

We then want to know if this equals $\newcommand{\re}{{\rm Re}} (\mathbb{1}+\Pi) \tilde{M}_{v} \mathbf {v}$ for some $\tilde{M}_{v}$ with elements $m_{ij} \geqslant 0 , \, i, j ={0},\dots,{3}$ such that each column and row sums to 1. We obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle (\mathbb{1}+\Pi) \tilde{M}_{v} \mathbf {v} = \left(\begin{array}{@{}c@{}} m_{{00}}+m_{{30}} \nonumber \\ m_{{00}}+m_{{10}} \nonumber \\ m_{{10}}+m_{{20}} \nonumber \\ m_{{20}}+m_{{30}} \end{array}\right). \nonumber \end{align} \tag{ A.51 }$$

Comparing this with equation (A.50), we have $m_{00}+m_{30}=m_{10}+m_{20}$ and $m_{20} + m_{30} = m_{00} + m_{10}$, which yields $m_{30}=m_{10}$, and concurrently $m_{00} + m_{30}=0$ and $m_{00} + m_{30}=1$, which means we arrive at a contradiction.

Here, via counterexample, we show that the following claim (discussed in section 3.2) does not hold in general:

Claim A.3. The relation $\newcommand{\re}{{\rm Re}} \frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{\mathbf {b}_i}{\Vert \mathbf {b}_i \Vert}$ holds $\forall$ $d\geqslant2$ and $\beta' \leqslant \beta$.

A.6.1. Counterexample to claim A.3.

Here, we present the proof of lemma 2 for any d-dimensional system, which ensures that condition (i) holds in general.

proof To prove that equation (48) holds, we need to show that for any $\beta> \beta^{\prime}$, $g(\beta):= \Vert {\bf q} \Vert$ is greater than $\newcommand{\re}{{\rm Re}} g(\beta^{\prime})= \Vert \mathbf {a} \Vert$ in d-dimensional systems. To achieve this, one may use the positivity of the first derivative of the function $g(\beta)$. We therefore calculate

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \partial_\beta \Vert {\bf q} \Vert &=\partial_\beta \left(\frac{\sum\nolimits_i {\rm e}^{-2\beta E_i}}{\sum\nolimits_{m,n} {\rm e}^{-\beta (E_m+E_n)}} \right) \nonumber \\ & = \frac{1}{z(\beta)^4}\sum_{i,m,n}(E_m+E_n-2E_i){\rm e}^{-\beta (2E_i+E_m+E_n)}\nonumber \\ & = \frac{2}{z(\beta)^3}\sum_{i,m}(E_m-E_i){\rm e}^{-\beta (2E_i+E_m)}\nonumber \\ & = \frac{2}{z(\beta)^3}\sum_{m > i}\big[(E_m-E_i){\rm e}^{-\beta (2E_i+E_m)}\nonumber \\ &\ \ \ \ \ + (E_i-E_m){\rm e}^{-\beta (E_i+2E_m)}\big]\nonumber \\ & = \frac{2}{z(\beta)^3}\sum_{m > i}(E_m-E_i)({\rm e}^{-\beta (2E_i+E_m)}-{\rm e}^{-\beta (E_i+2E_m)}).\nonumber \end{align} \tag{ A.52 }$$

Without loss of generality, we have ordered the energy eigenvalues in increasing order, $E_{i+1}\geqslant E_i$. Then, for any pair $(m,i)$ with m  >  i, we know that $({\rm e}^{-\beta (2E_i+E_m)}-{\rm e}^{-\beta (E_i+2E_m)})$ is nonnegative, implying $\partial_\beta \Vert {\bf q} \Vert \geqslant 0$. Also note that the last inequality is strict unless $E_i=E_m$ for all m  >  i implying $H \propto \mathbb{1}$ for which the problem is already solved. For all practical purposes, the inequality is therefore strict. □

proof Using our convention of LCSs, see section 2.2.1, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\bf r}_i:= \sum_{j=0}^{d-1} p_{j\, j+i} \,\mathbf {e}_{j}\ \ \ {\rm and}\ \ \ \mathbf {b}_i:= \sum_{j=0}^{d-1} p_{j\, j+i}^\prime \,\mathbf {e}_{j}, \nonumber \end{align} \tag{ A.53 }$$

where $p_{j\, j+i}(\beta):={\rm e}^{-\beta (E_j+E_{j+i})}/\mathcal{Z}(\beta){}^2$ and $p^\prime_{j\, j+i}:=p_{j\, j+i}(\beta')$. We then calculate

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{{\bf r}_i}{\Vert {\bf r}_i \Vert} &= \sum_{j=0}^{d-1} \frac{p_{j,j+i}}{\sum\nolimits_{k=0}^{d-1} p_{k\,k+i}} \,\mathbf {e}_{j} = \sum_{j=0}^{d-1} \frac{{\rm e}^{-\beta (E_j+E_{j+i})}}{\sum\nolimits_{k=0}^{d-1} {\rm e}^{-\beta (E_k+E_{k+i})}} \,\mathbf {e}_{j} \nonumber \\ &= \sum_{j=0}^{d-1} \frac{{\rm e}^{-\beta A_j}}{\sum\nolimits_{k=0}^{d-1} {\rm e}^{-\beta A_k}} \,\mathbf {e}_{j}, \nonumber \end{align} \tag{ A.54 }$$

where $A_j:= E_j+E_{j+i}$ for all $j=0,\dots, d-1$. The last expression is nothing else than the vectorised form of a thermal state at inverse temperature $\beta$ with respect to the Hamiltonian $\newcommand{\ket}[1]{| #1\rangle} \newcommand{\bra}[1]{\langle #1 |} \sum\nolimits_{j=0}^{d-1} A_j \ket{j}\!\bra{j}$. Since the vector of diagonal entries of any thermal state majorises the vectorised diagonal of any thermal state (with respect to the same Hamiltonian) with lower inverse temperature $\beta'< \beta$, which is nothing else than $\newcommand{\re}{{\rm Re}} \mathbf {b}_i/\Vert \mathbf {b}_i \Vert$, equation (49) holds, which concludes the proof. □

In this section, we will show how one can prove the existence of the STUs if condition (i) and the stronger version of condition (ii) from section 3.2 hold.

condition (i) ensures that there exist doubly-stochastic matrices $M_{r_i}$ which transform ${\bf r}_i/ \Vert {\bf r}_i\Vert$ to $\newcommand{\re}{{\rm Re}} \mathbf {b}_i/ \Vert \mathbf {b}_i\Vert$ and also ${\bf q}/ \Vert {\bf q}\Vert$ to $\newcommand{\re}{{\rm Re}} \mathbf {a}/ \Vert \mathbf {a} \Vert$, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle M_{r_i} \frac{{\bf r}_i}{\Vert {\bf r}_i\Vert}=\frac{\mathbf {b}_i}{\Vert \mathbf {b}_i\Vert} ~~\mathrm{for} ~ i=0,\dots,k. \label{eq:rimajbiddd} \nonumber \end{align} \tag{ A.55 }$$

Using equations (21) and (A.55), the marginals then become

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}} = M_q{\bf q} + \sum_{i=1}^{k} (\lfloor\frac{2i}{d}\rfloor+1)^{-1} (\mathbb{1}+ \Pi^i) \frac{\Vert {\bf r}_i\Vert}{\Vert \mathbf {b}_i\Vert}\mathbf {b}_i . \label{eq:EqualEvoMarginal Cprime11} \nonumber \end{align} \tag{ A.56 }$$

To reach a thermal state with higher temperature, one can compensate the required norm of $\newcommand{\re}{{\rm Re}} (\mathbb {1}+\Pi^i)\,\mathbf {b}_i$ by extra norm from the vector ${\bf q}$. To achieve that, we need to have $\newcommand{\re}{{\rm Re}} \Vert {\bf q}\Vert \geqslant \Vert \mathbf {a} \Vert$, which is ensured by condition (i). In addition, another requirement for shifting the norm is the existence of doubly-stochastic matrices $M_{q\to2 {b}_i}$ to transform ${\bf q}/ \Vert {\bf q}\Vert$ to $\newcommand{\re}{{\rm Re}} (\mathbb {1}+\Pi^i)\,\mathbf {b}_i/2\Vert \mathbf {b}_i\Vert$. Due to the HLP theorem, condition (ii) ensures that such matrices exist. Therefore, M_q can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_q= \alpha_0 \,M_{q\to a}+ \sum_{i=1}^{k} \alpha_i\, M_{q\to 2 {b}_i} , \label{eq: cov. com. M0} \nonumber \end{align} \tag{ A.57 }$$

where $M_{q\to a}$ and the $M_{q\to {b}_i}$ are doubly-stochastic matrices that map ${\bf q}$ to $\newcommand{\re}{{\rm Re}} \mathbf {a}$ and ${\bf q}$ to each of the $\newcommand{\re}{{\rm Re}} (\mathbb {1}+\Pi^i)\,\mathbf {b}_i/2\Vert \mathbf {b}_i\Vert$, respectively. Applying M_q to the vector ${\bf q}$ in equation (A.56), we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}} = \alpha_0\,\frac{\Vert {\bf q} \Vert}{\Vert \mathbf {a} \Vert} \,\mathbf {a}+ \sum_{i=1}^{k} (\alpha_i\,\frac{\Vert {\bf q} \Vert}{2\Vert \mathbf {b}_i \Vert}+(\lfloor\frac{2i}{d}\rfloor+1)^{-1}\frac{\Vert {\bf r}_i\Vert}{\Vert \mathbf {b}_i\Vert}) (\mathbb {1}+\Pi) \mathbf {b}_i. \label{eq:trans. 4d pass 2} \nonumber \end{align} \tag{ A.58 }$$

In order to obtain a thermal state at inverse temperature $\beta'\leqslant \beta$, one must choose the coefficients $\alpha_i$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:alphai} \alpha_0=\frac{\Vert \mathbf {a}\Vert}{\Vert {\bf q}\Vert},~\alpha_i=2(\lfloor\frac{2i}{d}\rfloor+1)^{-1}\frac{(\Vert \mathbf {b}_i\Vert -\Vert {\bf r}_i\Vert)}{\Vert {\bf q}\Vert}. \nonumber \end{align} \tag{ A.59 }$$

Using the fact that convex combinations of doubly-stochastic matrices are doubly-stochastic matrices, to ensure that M_q is indeed a doubly-stochastic matrix, one needs to show that the coefficients $\alpha_i$ as given by equation (A.59) are positive and sum to one. The inequality $\newcommand{\re}{{\rm Re}} \Vert \mathbf {b}_i\Vert \geqslant \Vert {\bf r}_i\Vert$ in condition (ii) ensures their positivity. Furthermore by using that $\newcommand{\re}{{\rm Re}} \Vert {\bf p} \Vert = \Vert {\tilde{\mathbf {p}}} \Vert$ we have that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \Vert {\bf p}\Vert =\Vert {\bf q}\Vert + 2\sum_{i=1}^{k} (\lfloor\frac{2i}{d}\rfloor+1)^{-1} {\Vert {\bf r}_i\Vert}=\Vert \mathbf {a}\Vert + 2\sum_{i=1}^{k} (\lfloor\frac{2i}{d}\rfloor+1)^{-1} {\Vert \mathbf {b}_i\Vert}, \nonumber \end{align} \tag{ A.60 }$$

which shows that $\alpha_0 +\sum\nolimits_{i=1}^{k} \alpha_i=1$, concluding the proof. □

Here we present detailed proofs of lemmas 3 and 4.

Proof of lemma 3. To prove $\newcommand{\re}{{\rm Re}} \frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{(\mathbb{1} + \Pi) b}{2 \Vert \mathbf {b} \Vert}$, we can employ the following two majorisation relations:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{(\mathbb{1} + \Pi) {\bf r}}{2 \Vert \mathbf {r} \Vert}, \label{eq:qmajrapp} \nonumber \end{align} \tag{ A.61 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{(\mathbb{1} + \Pi) {\bf r}}{2 \Vert \mathbf {r} \Vert} \succ \frac{(\mathbb{1}+\Pi) \mathbf {b}}{2 \Vert \mathbf {b} \Vert}. \label{eq:2rmaj2bapp} \nonumber \end{align} \tag{ A.62 }$$

If we show that these relations are true for any $\beta' \leqslant \beta$, we can say that our statement is proven.

To prove the relation in inequality (A.61), we need to show that the greatest/smallest entry of ${\bf q}/\Vert {\bf q} \Vert$ is greater/smaller than the greatest/smallest entry of $\newcommand{\re}{{\rm Re}} \frac{\mathbb{1}+\Pi}{2} \frac{{\bf r}}{\Vert \mathbf {r} \Vert}$. That is, we first need to check that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{00}}{p_{00}+p_{11}+p_{22}} \geqslant \frac{p_{01}+p_{20}}{2(\,p_{01}+p_{12}+p_{20})}. \nonumber \end{align} \tag{ A.63 }$$

Disregarding the trivial case where p ₀  =  1, and $p_{1}=p_{2}=0$, this inequality can be transformed to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &{\rm e}^{-\beta E_1} + 2 {\rm e}^{-\beta E_1} {\rm e}^{- \beta E_2} + {\rm e}^{-\beta E_2}\nonumber \\ &\ \ \geqslant ({\rm e}^{- \beta E_1}+ {\rm e}^{-\beta E_2}) ({\rm e}^{- 2 \beta E_1}+ {\rm e}^{-2 \beta E_2}). \label{eq:mystery inequality} \nonumber \end{align} \tag{ A.64 }$$

This inequality indeed holds, since the left-hand side is larger than (or equal to) $({\rm e}^{-\beta E_1}+{\rm e}^{-\beta E_2}){}^2$, which, in turn, is larger or equal to the righ-hand side of inequality (A.64).

Then, second, we need to check the inequality

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{22}}{p_{00}+p_{11}+p_{22}} {\leqslant}\frac{p_{20}+p_{12}}{2 (\,p_{01}+p_{12}+p_{20})}, \nonumber \end{align} \tag{ A.65 }$$

where the relevant case (i.e. for $p_{1}\neq1$ such that at least $p_{1}\neq0$) can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & {\rm e}^{-2 \beta E_2} ({\rm e}^{-\beta (E_1+E_2)} + 2 {\rm e}^{-\beta E_1} + {\rm e}^{-\beta E_2})\nonumber \\ & \ \ \leqslant ({\rm e}^{-\beta E_2}+ {\rm e}^{-\beta (E_1+E_2)}) (1+{\rm e}^{- 2 \beta E_1}). \label{eq:mystery inequality 2} \nonumber \end{align} \tag{ A.66 }$$

This second inequality holds since the right-hand side of inequality (A.66) is larger than (or equal to) $({\rm e}^{-\beta E_2}+{\rm e}^{-\beta (E_1+E_2)}){}^2$, which, in turn, is larger or equal to the righ-hand side of inequality (A.66).

To prove equation (A.62), we show that the greatest entry of $\newcommand{\re}{{\rm Re}} \frac{(\mathbb{1}+\Pi){\bf r}}{2 \Vert \mathbf {r} \Vert}$ is monotonically increasing with $\beta$, and that its smallest entry is monotonically decreasing with $\beta$. For the former we calculate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \partial_{\beta} &\frac{p_{01}+p_{20}}{2(\,p_{01}+p_{12}+p_{20})}= \partial_{\beta} \frac{{\rm e}^{-\beta E_1}+{\rm e}^{-\beta E_2}}{2({\rm e}^{- \beta E_1}+ {\rm e}^{-\beta (E_1+E_2)}+{\rm e}^{-\beta E_2})}\nonumber \\ & = \frac{E_2 {\rm e}^{-\beta (2E_1+E_2)}+E_1 {\rm e}^{-\beta (E_1+2 E_2)}}{2({\rm e}^{-\beta E_1}+{\rm e}^{-\beta (E_1+E_2)}+{\rm e}^{-\beta E_2})^2}\geqslant 0. \nonumber \end{align} \tag{ A.67 }$$

And for the latter

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \partial_{\beta} &\frac{p_{20}+p_{21}}{2(\,p_{01}+p_{12}+p_{20})}= \partial_{\beta} \frac{{\rm e}^{-\beta E_2}+{\rm e}^{-\beta (E_1+E_2)}}{2({\rm e}^{- \beta E_1}+ {\rm e}^{-\beta (E_1+E_2)}+{\rm e}^{-\beta E_2})}\nonumber \\ & = \frac{(E_1-E_2) {\rm e}^{-\beta (E_1+E_2)}-E_2 {\rm e}^{-\beta (2E_1+ E_2)}}{2({\rm e}^{-\beta E_1}+{\rm e}^{-\beta (E_1+E_2)}+{\rm e}^{-\beta E_2})^2}\leqslant 0. \nonumber \end{align} \tag{ A.68 }$$

□

proof We prove this lemma under the stated condition on the gap structure of the local Hamiltonians, i.e. $\delta_{i}{\geqslant} \delta_{i+1}$. Here we need to show that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{\mathbb{1}+\Pi}{2} \frac{\mathbf {b}_1}{\Vert \mathbf {b}_1 \Vert}{,}\label{qmajb14d} \nonumber \end{align} \tag{ A.69 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{\mathbb{1}+\Pi^2}{2} \frac{\mathbf {b}_2}{\Vert \mathbf {b}_2 \Vert}= \frac{\mathbf {b}_2}{\Vert \mathbf {b}_2 \Vert}\label{qmajb24d}. \nonumber \end{align} \tag{ A.70 }$$

To do so, we argue that the majorisation relations $\frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{\mathbb{1}+\Pi^i}{2} \frac{{\bf r}_i}{\Vert {\bf r}_i \Vert}$ and $\newcommand{\re}{{\rm Re}} \frac{\mathbb{1}+\Pi^i}{2} \frac{{\bf r}_i}{\Vert {\bf r}_i \Vert}\succ \frac{\mathbb{1}+\Pi^i}{2} \frac{\mathbf {b}_i}{\Vert \mathbf {b}_i \Vert}$, for $i=1,2{,}$ hold.

Let us start with equation (A.69). The condition $\delta_{i} \geqslant \delta_{i+1}$ sets the ordering of $\frac{\mathbb{1}+\Pi}{2} \frac{{{\bf r}_{1}}}{\Vert {{\bf r}_{1}} \Vert}$. Still, some ambiguity remains as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p_{01}+p_{30} \geqslant p_{12}+p_{01} \geqslant p_{30}+p_{23} \geqslant p_{23}+p_{12}. \nonumber \end{align} \tag{ A.71 }$$

Thus, we need to prove the following three inequalities:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{00}}{p_{00}+p_{11}+p_{22}+p_{33}} \geqslant \frac{p_{01}+p_{30}}{2 (\,p_{01}+p_{12}+p_{23}+p_{30})}{,}\label{item a4} \nonumber \end{align} \tag{ A.72a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{00}+p_{11}}{p_{00}+p_{11}+p_{22}+p_{33}} \geqslant \frac{p_{01}+p_{30}+p_{12}+p_{01}}{2 (\,p_{01}+p_{12}+p_{23}+p_{30})}{,}\label{item b4} \nonumber \end{align} \tag{ A.72b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{33}}{p_{00}+p_{11}+p_{22}+p_{33}} \leqslant \frac{p_{32}+p_{12}}{2 (\,p_{01}+p_{12}+p_{23}+p_{30})}\label{item c4}. \nonumber \end{align} \tag{ A.72c }$$

Inequality (A.72a) can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{1}{1+{\rm e}^{-2 \beta E_1}+{\rm e}^{-2 \beta E_2}+{\rm e}^{-2 \beta E_3}} \geqslant \frac{1}{2 \left(1+ \frac{p_{12} + p_{23}}{p_{01}+p_{30}}\right)}, \nonumber \end{align} \tag{ A.73 }$$

which can further be turned into the inequality

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1+1+{\rm e}^{-\beta E_2}+{\rm e}^{-\beta E_2} \geqslant 1+{\rm e}^{-2 \beta E_1}+{\rm e}^{-2 \beta E_2}+{\rm e}^{-2 \beta E_3}, \nonumber \end{align} \tag{ A.74 }$$

which holds because $E_1 \geqslant 0$ and $E_3 \geqslant E_2$.

Similarly, inequality (A.72b) can be recast as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{1}{1+\frac{p_{22}+p_{33}}{p_{00}+p_{11}}} \geqslant \frac{1}{2} \left(\frac{1}{1+\frac{p_2}{p_0}}+\frac{1}{1+\frac{p_3}{p_1}} \right) \nonumber \end{align} \tag{ A.75 }$$

which implies the inequality

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f\left(\frac{p_2^2+p_3^2}{p_0^2+p_1^2}\right) \geqslant \frac{1}{2} \left[\,f\left(\frac{p_2}{p_0}\right)+f\left(\frac{p_3}{p_1}\right)\right], \nonumber \end{align} \tag{ A.76 }$$

where we have introduced $f(x) := 1/(1+x)$. We can see that $f'(x) \leqslant 0$ and $\frac{p_2^2+p_3^2}{p_0^2+p_1^2} \leqslant \frac{p_2}{p_0}$, and in the case of $E_3 \leqslant E_1 +E_2$ also $\frac{p_2}{p_0} \leqslant \frac{p_3}{p_1}$. Thus,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f\left(\frac{p_3}{p_1} \right) \leqslant f \left(\frac{p_2}{p_0} \right) \leqslant f \left(\frac{p_2^2 +p_3^2}{p_0^2+p_1^2} \right), \nonumber \end{align} \tag{ A.77 }$$

from which

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{1}{2} \left[\,f\left(\frac{p_3}{p_1} \right) + f \left(\frac{p_2}{p_0} \right)\right] \leqslant f \left(\frac{p_2^2 +p_3^2}{p_0^2+p_1^2} \right) \nonumber \end{align} \tag{ A.78 }$$

follows. This proves inequality (A.72b).

Inequality (A.72c) is equivalent to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{1}{{\rm e}^{2 \beta E_3}+{\rm e}^{2 \beta (E_3-E_1)}+ {\rm e}^{2 \beta (E_3-E_2)}+1} \leqslant \frac{1}{2 \left(1 + \frac{p_{01}+p_{30}}{p_{32}+p_{12}} \right)}, \nonumber \end{align} \tag{ A.79 }$$

which can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1 + 2\frac{p_0}{p_2} \leqslant {\rm e}^{2 \beta E_3}+{\rm e}^{2 \beta (E_3-E_1)}+ {\rm e}^{2 \beta (E_3-E_2)}, \nonumber \end{align} \tag{ A.80 }$$

which is true since the right-hand side is larger or equal ${\rm e}^{2 \beta E_2}+2$, which, in turn, is larger than (or equal to) the left-hand side, because $(1-{\rm e}^{\beta E_2}){}^2 \geqslant 0$. This proves inequality (A.72c).

Thus far, we have shown that $\frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{\mathbb{1}+\Pi}{2} \frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert}$. To complete the first part of the proof, we also need to show that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \frac{\mathbb{1}+\Pi}{2} \frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert} \succ \frac{\mathbb{1}+\Pi}{2} \frac{\mathbf {b}_1}{\Vert \mathbf {b}_1 \Vert}. \label{eq-tmp-a} \nonumber \end{align} \tag{ A.81 }$$

To do so, we first prove that the derivatives with respect to $\beta$ of the two largest elements of $\frac{\mathbb{1}+\Pi}{2} \frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert}$ are positive and those of the two smallest elements are negative. Let us represent the elements of $\frac{\mathbb{1}+\Pi}{2} \frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert}$ as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\mathbb{1}+\Pi}{2} \frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert}= \bigl(\,f_{{0}}(\beta), f_{{1}}(\beta), f_{{2}}(\beta), f_{{3}}(\beta) \bigr)^{T}. \nonumber \end{align} \tag{ A.82 }$$

Note that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{{0}}(\beta)&= \frac{p_{01} + p_{30}}{2(\,p_{01}+p_{12}+p_{23}+p_{30})}=\frac{1}{2} \frac{p_0(\,p_1+p_3)}{(\,p_0+p_2)(\,p_1+p_3)}\nonumber \\ &= \frac{1}{2} \frac{1}{1+{\rm e}^{-\beta E_2}}{,} \nonumber \end{align} \tag{ A.83 }$$

which yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:f1deri} f'_{{0}}(\beta)= \frac{1}{2} \frac{E_2\,{\rm e}^{-\beta E_2}}{(1+{\rm e}^{-\beta E_2})^2} \geqslant 0. \nonumber \end{align} \tag{ A.84 }$$

Similarly, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:fderi} f_{{1}}(\beta)&= \frac{1}{2} \frac{1}{1+{\rm e}^{-\beta (E_3-E_1)}} , \,f'_{{1}}(\beta)= \frac{(E_3-E_1){\rm e}^{-\beta (E_3-E_1)}}{2(1+{\rm e}^{-\beta (E_3-E_1)})^2} \geqslant 0{,}\nonumber \\ f_{{2}}(\beta)&= \frac{1}{2} \frac{1}{1+{\rm e}^{\beta E_2}} ,\,f'_{{2}}(\beta)= \frac{-E_2\,{\rm e}^{\beta E_2}}{2(1+{\rm e}^{\beta E_2})^2} \leqslant 0{,}\nonumber \\ f_{{3}}(\beta)&= \frac{1}{2} \frac{1}{1+{\rm e}^{\beta (E_3-E_1)}} ,\,f'_{{3}}(\beta)= \frac{-(E_3-E_1){\rm e}^{\beta (E_3-E_1)}}{2(1+{\rm e}^{\beta (E_3-E_1)})^2} \leqslant 0. \nonumber \end{align} \tag{ A.85 }$$

In the regime $E_3 \leqslant E_1+E_2$, the components are ordered as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{{0}}(\beta) \geqslant f_{{4}}(\beta) \geqslant f_{{3}}(\beta) \geqslant f_{{2}} (\beta), \nonumber \end{align} \tag{ A.86 }$$

and thus one needs to show for $\beta' \leqslant \beta$ that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{{0}}(\beta) &\geqslant f_{{0}}(\beta'){,}\nonumber \\ f_{{0}}(\beta)+f_{{1}}(\beta) &\geqslant f_{{0}}(\beta')+f_{{1}}(\beta'){,}\nonumber \\ f_{{2}}(\beta) &\leqslant f_{{2}}(\beta'){,}\nonumber \\ f_{{2}}(\beta)+f_{{3}}(\beta) &\leqslant f_{{2}}(\beta')+f_{{3}}(\beta').\nonumber \end{align} \tag{ A.87 }$$

These relations follow straightforwardly from equations (A.84) and (A.85), noting that a function f _i(x) that fulfils $f'_{i}(x)\leqslant 0$ ($f'_{i}(x)\geqslant 0$) the relation $f_{i}(x_1)\leqslant f_{i}(x_{2})$ ($f_{i}(x_1)\geqslant f_{i}(x_{2})$) is also satisfied for $x_{1}\geqslant x_{2}$.

In a similar fashion, in the regime $E_3 \geqslant E_1+E_2$ the ordering of the components is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{{1}}(\beta) \geqslant f_{{0}}(\beta) \geqslant f_{{2}}(\beta) \geqslant f_{{3}} (\beta), \nonumber \end{align} \tag{ A.88 }$$

and thus one needs to show for $\beta' \leqslant \beta$ that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{{1}}(\beta) &\geqslant f_{{1}}(\beta'){,}\nonumber \\ f_{{0}}(\beta)+f_{{1}}(\beta) &\geqslant f_{{0}}(\beta')+f_{{1}}(\beta'){,}\nonumber \\ f_{{3}}(\beta) &\leqslant f_{{3}}(\beta'){,}\nonumber \\ f_{{2}}(\beta)+f_{{3}}(\beta) &\leqslant f_{{2}}(\beta')+f_{{3}}(\beta').\nonumber \end{align} \tag{ A.89 }$$

These relations also follow straightforwardly from equations (A.84) and (A.85).

We next turn our attention to equation (A.70). Because $p_{02} \geqslant p_{13}$, we need to prove the following inequalities:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{00}}{p_{00}+p_{11}+p_{22}+p_{33}} \geqslant \frac{p_{02}}{2 (\,p_{02}+p_{13})}{,}\label{item aa4} \nonumber \end{align} \tag{ A.90a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{00}+p_{11}}{p_{00}+p_{11}+p_{22}+p_{33}} \geqslant \frac{2p_{02}}{2 (\,p_{02}+p_{13})}{,}\label{item bb4} \nonumber \end{align} \tag{ A.90b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{33}}{p_{00}+p_{11}+p_{22}+p_{33}} \leqslant \frac{p_{13}}{2 (\,p_{02}+p_{13})}\label{item cc4}. \nonumber \end{align} \tag{ A.90c }$$

Inequality (A.90a) can be simplified to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{0}}{p_{00}+p_{11}+p_{22}+p_{33}} \geqslant \frac{p_{2}}{2 (\,p_{02}+p_{13})}, \nonumber \end{align} \tag{ A.91 }$$

where we can use the inequalities $p_{0}p_{02}\geqslant p_{2}p_{11}$ and $p_{0}p_{13}\geqslant p_2 p_{33}$ to arrive at the condition $p_{0}p_{13}\geqslant p_{2}p_{22}$. Due to our assumption on the energy eigenvalues, we have $p_{03}\geqslant p_{22}$, which proves inequality (A.90a).

Inequality (A.90b) can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{00}+p_{11}}{p_{00}+p_{11}+p_{22}+p_{33}} \geqslant \frac{2p_{02}}{2 (\,p_{02}+p_{13})}, \nonumber \end{align} \tag{ A.92 }$$

which implies $p_{13}(\,p_{00}+p_{11}) \geqslant p_{02}(\,p_{22}+p_{33})$. Moreover, because $p_{00}+p_{11}\geqslant 2 p_{01}$, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p_{13} (\,p_{00}+p_{11}) &\geqslant 2 p_{01} p_{13}\geqslant p_{01} p_{13}+p_{02} p_{33}\nonumber \\ &\geqslant p_{02} p_{22} + p_{02} p_{33}, \nonumber \end{align} \tag{ A.93 }$$

where in the last step the condition $\delta_i \geqslant \delta_{i+1}$ was used. This proves inequality (A.90b).

Inequality (A.90c) can be recast as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{3}}{p_{00}+p_{11}+p_{22}+p_{33}} \leqslant \frac{p_{1}}{2 (\,p_{02}+p_{13})}, \nonumber \end{align} \tag{ A.94 }$$

from which we arrive at

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 2p_3 p_{02}+ p_3p_{13} &\leqslant p_1 (\,p_{01}+p_{10}+p_{22})\nonumber \\ &\leqslant p_1 (\,p_{00}+p_{11}+p_{22}), \nonumber \end{align} \tag{ A.95 }$$

which is true for any energy spectrum. □

Unfortunately, equation (54) in lemma 4 fails when $\delta_i \ngeq \delta_{i+1}$.

A.9.1. Counterexample to equation (54) for $\boldsymbol{\delta_i \,{<}\, \delta_{i+1}}$.

The relation $\frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{\mathbb{1}+\Pi}{2} \frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert}$ fails in the regime $\delta_i \,{<}\, \delta_{i+1}$. In that regime, the greatest element of $\frac{\mathbb{1}+\Pi}{2} \frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert}$ is $\frac{p_{12}+p_{01}}{2(\,p_{01}+p_{12}+p_{23}+p_{30})}$ and so

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{00}}{p_{00}+p_{11}+p_{22}+p_{33}} \geqslant \frac{p_{12}+p_{01}}{2(\,p_{01}+p_{12}+p_{23}+p_{30})} \nonumber \end{align} \tag{ A.96 }$$

must hold in order for the relation to be true. But in the limit $E_3 \rightarrow \infty, E_1 \rightarrow 0, E_2 \rightarrow 0$ the above inequality becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_{00}}{3 p_{00}} \geqslant \frac{2 p_{00}}{2(\,p_{00}+p_{00})}\Leftrightarrow \frac{1}{3} \geqslant \frac{1}{2}, \nonumber \end{align} \tag{ A.97 }$$

which is obviously invalid. The relation $\frac{{\bf q}}{\Vert {\bf q} \Vert} \succ \frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert}$ also fails sometimes. Looking at the limit $E_3 \rightarrow \infty$, $\frac{{\bf q}}{\Vert {\bf q} \Vert} \rightarrow (\,p_{00}, p_{11}, p_{22}, 0)$ and $\frac{{\bf r}_1}{\Vert {\bf r}_1 \Vert} \rightarrow (\,p_{02}, 0 , p_{20}, 0)$. For the majorisation relation to hold, we would need $p_{22} \leqslant 0$, which is clearly not true in general.

Proof of equation (55) in lemma 4.

proof We need to prove that $\newcommand{\re}{{\rm Re}} \Vert {\bf r}_i\Vert \leqslant \Vert \mathbf {b}_i\Vert$, for $i=1,2$, for d  =  4. To do so, we again use the first derivative of $\Vert {\bf r}_i\Vert$ and show that it is always negative. From

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Vert {\bf r}_i\Vert =\sum_{j}p_j p_{j+i}=\frac{\sum\nolimits_j {\rm e}^{-\beta(E_j+E_{j+i})}}{\sum\nolimits_{m,n}{\rm e}^{-\beta(E_m+E_n)}}{,} \nonumber \end{align} \tag{ A.98 }$$

the partial derivative with respect to $\beta$ reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\partial_\beta \Vert {\bf r}_i\Vert =\sum_{j,m,n}(E_m+E_n-E_j-E_{j+i})\frac{{\rm e}^{-\beta(E_m+E_n+E_j+E_{j+i})}}{z(\beta)^4}\nonumber \\ &=\frac{1}{z(\beta)^3}\sum_{j,m}(2E_m-E_j-E_{j+i}){\rm e}^{-\beta(E_m+E_j+E_{j+i})}. \label{rjprime} \nonumber \end{align} \tag{ A.99 }$$

Now we need to show that for any $i\neq 0$, $\partial_\beta \Vert {{\bf r}_{i}}\Vert$ is negative. For i  =  1 we find

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \partial_\beta \Vert {\bf r}_1\Vert &=\frac{-1}{z(\beta)^3}\big[E_1({\rm e}^{-\beta E_1}+{\rm e}^{-\beta(E_1+E_2)}-{\rm e}^{-2\beta E_1})+(E_2-E_1)(-{\rm e}^{-\beta(E_1+E_2)}+{\rm e}^{-\beta(2E_1+E_2)}\nonumber \\ &\quad +{\rm e}^{-\beta(E_1+E_2+E_3)}-{\rm e}^{-\beta(E_1+2E_2)})+(E_3-E_2)({\rm e}^{-\beta(E_2+E_3)}-{\rm e}^{-\beta(E_1+E_2+E_3)}+{\rm e}^{-\beta(2E_2+E_3)}-{\rm e}^{-\beta(E_2+2E_3)})\nonumber \\ &\quad +E_3({\rm e}^{-\beta(E_2+E_3)}+{\rm e}^{-\beta E_3}-{\rm e}^{-\beta(E_1+ E_3)}-{\rm e}^{-2\beta E_3})\big]\leqslant 0. \label{r1prime} \nonumber \end{align} \tag{ A.100 }$$

Since we have ordered the energy eigenvalues in the increasing order, it is straightforward to show that $\partial_\beta \Vert {\bf r}_1\Vert$ is always nonpositive for any set of energy eigenvalues in dimension 4.

Now similarly, we can show that $\newcommand{\re}{{\rm Re}} \Vert {\bf r}_2\Vert \leqslant \Vert \mathbf {b}_2\Vert$. By employing equation (A.99), $\partial_\beta \Vert {\bf r}_2\Vert$ is obtained as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \partial_\beta \Vert {\bf r}_2 \Vert &=\frac{-2}{z(\beta)^3}\big[E_1({\rm e}^{-\beta E_2}+2{\rm e}^{-\beta (E_1 +E_3)}-{\rm e}^{-\beta (E_1+ E_2)}-{\rm e}^{-2 \beta E_2}-{\rm e}^{-\beta (E_2+E_3)})+(E_2-E_1)({\rm e}^{-\beta E_2}+{\rm e}^{-\beta (E_1 +E_3)}\nonumber \\ &\quad +{\rm e}^{-\beta (E_1+ E_2)}+{\rm e}^{-\beta (2E_1+ E_3)}-{\rm e}^{-2 \beta E_2}-{\rm e}^{-\beta (E_1+ E_2+E_3)}-{\rm e}^{-\beta (E_2+E_3)}-{\rm e}^{-\beta (E_1+2E_3)})\nonumber \\ &\quad +(E_3-E_2)({\rm e}^{-\beta (E_1 +E_3)}+{\rm e}^{-\beta (2E_1+ E_3)}+{\rm e}^{-\beta (E_1+ E_2+E_3)}-2{\rm e}^{-\beta (E_2+E_3)}-{\rm e}^{-\beta (E_1+2E_3)})\big]. \nonumber \end{align} \tag{ A.101 }$$

The above expression is also nonpositive under the assumption $\delta_i\,{\geqslant}\, \delta_{i+1}$. □

In this section, we will show that the condition (II) mentioned in section 3.3 is true for every choice of the set $\{E_i\}_{i=0}^{d-1}$ and any initial inverse temperature.

Expressing a general point in the curve ${\bf p}(\beta')=(x_0(\beta'),\dots,x_{d-2}(\beta'),-1)$ as a linear combination of the vertices $\newcommand{\re}{{\rm Re}} \mathbf {v}_i$ in condition (I) we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(\begin{array}{@{}c@{}}x_0(\beta') \nonumber \\ x_1(\beta') \nonumber \\ \vdots \nonumber \\ x_{d-2}(\beta') \nonumber \\ -1 \end{array}\right)=\left(\begin{array}{@{}ccccc@{}} x_0(\beta) & 0 & 0 &\dots & 0 \nonumber \\ x_1(\beta) & x_1(\beta) & 0 &\dots & 0 \nonumber \\ \vdots & \ddots & x_2(\beta) &\ddots & 0 \nonumber \\ x_{d-2}(\beta) & \dots & \ddots &\ddots & 0 \nonumber \\ -1 & -1 & -1 &\dots & -1 \end{array}\right) \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_1 \nonumber \\ \vdots \nonumber \\ a_{d-2} \nonumber \\ a_{d-1} \end{array}\right), \nonumber \end{align} \tag{ A.102 }$$

where the $(a_0,\dots,a_{d-1})$ are coefficients given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle a_i = \frac{x_i(\beta')}{x_i(\beta)} - \frac{x_{i-1}(\beta')}{x_{i-1}(\beta)} = \frac{x_i(\beta')}{x_{i-1}(\beta)} \left(\frac{x_{i-1}(\beta)}{x_i(\beta)} - \frac{x_{i-1}(\beta')}{x_i(\beta')} \right) , \nonumber \end{align} \tag{ A.103 }$$

which are positive for all i and satisfy $\sum\nolimits_i a_i=1$ (and thus below 1 for all i) if the following condition holds

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\rm d}{{\rm d} \beta}\Bigl(\frac{x_{m}}{x_{m+1}} (\beta) \Bigr) =\frac 1 {x_{m+1}}\,\Bigl(\frac{\partial x_{m}}{\partial \beta} - \frac{x_{m}}{x_{m+1}}\frac{\partial x_{m+1}}{\partial \beta}\Bigr)\,\geqslant\,0\ \forall m , \nonumber \end{align} \tag{ A.104 }$$

which means that the function $x_i/x_{i+1} (\beta)$ is monotonically decreasing with $\beta$. This in turn means that the final point can be reached as a convex combination of the vertices.

With a convenient relabeling of the index we can write $x_{m+1}/x_{m}(\beta)$ in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{x_{m}}{x_{m-1}} &=\frac{(m+1){\rm e}^{-\beta E_{m+1}}-\sum\nolimits_{i=0}^{m}{\rm e}^{-\beta E_i}}{me^{-\beta E_{m}}-\sum\nolimits_{i=0}^{m-1}{\rm e}^{-\beta E_i}}\nonumber \\ &=\frac{(m+1){\rm e}^{-\beta E_{m+1}}-me^{-\beta E_{m}}+me^{-\beta E_{m}}-\sum\nolimits_{i=0}^{m-1}{\rm e}^{-\beta E_i}-{\rm e}^{-\beta E_{m}}}{me^{-\beta E_{m}}-\sum\nolimits_{i=0}^{m-1}{\rm e}^{-\beta E_i}}\nonumber \\ &= 1+\frac{(m+1)({\rm e}^{-\beta E_{m+1}}-{\rm e}^{-\beta E_{m}})}{me^{-\beta E_{m}}-\sum\nolimits_{i=0}^{m-1}{\rm e}^{-\beta E_i}}. \nonumber \end{align} \tag{ A.105 }$$

To prove our claim (condition (II)), we need to show that for any $\beta\geqslant0$ the function $x_{m}/x_{m-1}$ is monotonically decreasing with $\beta$. That is, we have to show that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R(\beta) = \frac{\rm d}{{\rm d}\beta}\Bigl(\frac{x_{m}}{x_{m-1}} (\beta) \Bigr) \,\leqslant\, 0. \label{ineq R} \nonumber \end{align} \tag{ A.106 }$$

The derivative takes the form

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle R(\beta) &=\frac{m+1}{x_{m-1}^2} \Bigl[-(E_{m+1}{\rm e}^{-\beta E_{m+1}}-E_{m}{\rm e}^{-\beta E_{m}})(me^{-\beta E_{m}}-\sum_{i=0}^{m-1}{\rm e}^{-\beta E_i})+({\rm e}^{-\beta E_{m+1}}-{\rm e}^{-\beta E_{m}})(m\,E_{m}{\rm e}^{-\beta E_{m}}-\sum_{i=0}^{m-1}E_{i}{\rm e}^{-\beta E_i})\Bigr]\nonumber \\ &=\frac{m+1}{x_{m-1}^2}\Bigl(\sum_{i=0}^{m-1}[(E_{m+1}-E_i){\rm e}^{-\beta(E_{m+1}+E_i)}-(E_{m}-E_i){\rm e}^{-\beta(E_{m}+E_i)}]-m (E_{m+1}-E_m){\rm e}^{-\beta(E_{m+1}+E_m)}\Bigr)\nonumber \\ &=\frac{m+1}{x_{m-1}^2}\Bigl(\sum_{i=0}^{m-1}[(E_{m+1}-E_i){\rm e}^{-\beta(E_{m+1}+E_i)}-(E_{m}-E_i){\rm e}^{-\beta(E_{m}+E_i)}]-m (E_{m+1}-E_i+E_i-E_m){\rm e}^{-\beta(E_{m+1}+E_m)}\Bigr)\nonumber \\ &=\frac{m+1}{x_{m-1}^2}\sum_{i=0}^{m-1}\bigl[(E_{m+1}-E_i)({\rm e}^{-\beta(E_{m+1}+E_i)}-{\rm e}^{-\beta(E_{m+1}+E_m)})-(E_{m}-E_i)({\rm e}^{-\beta(E_{m}+E_i)}-{\rm e}^{-\beta(E_{m+1}+E_m)})\bigr]\nonumber \\ &=\frac{m+1}{x_{m-1}^2}{\rm e}^{-\beta(E_{m+1}+E_m)}\sum_{i=0}^{m-1}\bigl[(E_{m+1}-E_i)({\rm e}^{\beta(E_{m}-E_i)}-1)-(E_{m}-E_i)({\rm e}^{\beta(E_{m+1}-E_i)}-1)\bigr]. \nonumber \end{align} \tag{ A.107 }$$

To prove the nonpositivity of $R(\beta)$ we need to show that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (E_{m+1}-E_i)({\rm e}^{\beta(E_{m}-E_i)}-1)-(E_{m}-E_i)({\rm e}^{\beta(E_{m+1}-E_i)}-1)\leqslant 0 \ \ \ {\rm for}\ \ E_i\leqslant E_m \leqslant E_{m+1}. \label{eq: negativity R} \nonumber \end{align} \tag{ A.108 }$$

To do so, we define $y_m:=\beta(E_{m}-E_i)$, and write inequality (A.108) in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h(y_{m+1}) = \frac{y_{m+1}}{{\rm e}^{y_{m+1}}-1}\,\leqslant\, \frac{y_{m}}{{\rm e}^{y_{m}}-1}. \nonumber \end{align} \tag{ A.109 }$$

This inequality is satisfied for all $y_{m+1}\geqslant y_{m}$ if the function $h(y)$ is monotonically decreasing. To see that this is the case, we calculate the derivative

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial h(y)}{\partial y} = \frac{(1-y){\rm e}^{y}-1}{({\rm e}^{y}-1)^2}, \nonumber \end{align} \tag{ A.110 }$$

which can be seen to be smaller or equal to zero since $(1-y){\rm e}^{y}$ is a monotonically decreasing function for $y\geqslant0$, i.e. $\frac{\partial}{\partial y}(1-y){\rm e}^{y}=-y {\rm e}^{y}\leqslant0$ and hence has its maximum value of 1 at y   =  0, confirming that $(1-y){\rm e}^{y}-1\leqslant0$ and that $h(y)$ is monotonically decreasing.

A.10.1. Explicit partial derivatives in dimension 3.

To show explicitly that the curve $x(\beta)$ is monotonically increasing in d  =  3, we first calculate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial x}{\partial\beta} = - \mathcal{Z}^{-2}f_{x}(\beta), \label{eq:partial x wrt beta1} \nonumber \end{align} \tag{ A.111a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial y}{\partial\beta} = - 3\mathcal{Z}^{-2}f_{y}(\beta), \label{eq:partial y wrt beta1} \nonumber \end{align} \tag{ A.111b }$$

where the functions $f_{x}(\beta)$ and $f_{y}(\beta)$ are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{x}(\beta) = E_{1}{\rm e}^{-\beta E_{1}}(2+{\rm e}^{-\beta E_{2}})+E_{2}{\rm e}^{-\beta E_{2}}(1-{\rm e}^{-\beta E_{1}})\,\geqslant\,0, \label{eq:fx1} \nonumber \end{align} \tag{ A.112a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{y}(\beta) = (E_{2}-E_{1}){\rm e}^{-\beta (E_{1}+E_{2})}+E_{2}{\rm e}^{-\beta E_{2}}\,\geqslant\,0. \label{eq:fy1} \nonumber \end{align} \tag{ A.112b }$$

With this, we can then evaluate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial y}{\partial x} = \left(\frac{\partial x}{\partial\beta}\right)^{-1}\frac{\partial y}{\partial\beta} = \frac{3f_{y}(\beta)}{f_{x}(\beta)}\,\geqslant\,0, \label{eq:partial y wrt x1} \nonumber \end{align} \tag{ A.113 }$$

which allows us to confirm that the curve segment lies above (with respect to the coordinates x and y  in the plane of the polytope) the line connecting $\newcommand{\bi}{\boldsymbol} \bigl(x(\beta),y(\beta)\bigr)$ and $\newcommand{\bi}{\boldsymbol} \bigl(\tilde{x},\tilde{y}\bigr)$.

To show that this curve is a convex function, we need to show that $\partial^{2}y/\partial x^{2}\geqslant0$, and calculate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial^{2}y}{\partial x^{2}} & = \left(\frac{\partial x}{\partial\beta}\right)^{-1}\frac{\partial}{\partial\beta}\,\left[\left(\frac{\partial x}{\partial\beta}\right)^{-1}\frac{\partial y}{\partial\beta}\right]\nonumber \\ & = \left(\frac{\partial x}{\partial\beta}\right)^{-2}\left\{\frac{\partial^{2}y}{\partial \beta^{2}}-\left(\frac{\partial y}{\partial \beta}\right)\left(\frac{\partial x}{\partial\beta}\right)^{-1}\left(\frac{\partial^{2}x}{\partial \beta^{2}}\right)\right\}. \label{eq:2nd partial y wrt x} \nonumber \end{align} \tag{ A.114 }$$

The second partial derivatives with respect to $\beta$ are found to be

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial x^{2}}{\partial\beta^{2}} = \mathcal{Z}^{-3}\bigl[2\frac{\partial\mathcal{Z}}{\partial\beta}f_{x}-\mathcal{Z}\frac{\partial f_{x}}{\partial\beta}\bigr], \label{eq:2nd partial x wrt beta} \nonumber \end{align} \tag{ A.115a }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial y^{2}}{\partial\beta^{2}} = 3\mathcal{Z}^{-3}\bigl[2\frac{\partial\mathcal{Z}}{\partial\beta}f_{y}-\mathcal{Z}\frac{\partial f_{y}}{\partial\beta}\bigr]. \label{eq:2nd partial y wrt beta} \nonumber \end{align} \tag{ A.115b }$$

Combining this with equation (A.114) we find

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial^{2}y}{\partial x^{2}} = 3\,\mathcal{Z}^{2}f_{x}^{-3}\bigl(\,f_{y}\frac{\partial f_{x}}{\partial\beta}-f_{x}\frac{\partial f_{y}}{\partial\beta}\bigr). \label{eq:2nd partial y wrt x rewritten} \nonumber \end{align} \tag{ A.116 }$$

The derivatives of $f_{x}(\beta)$ and $f_{y}(\beta)$ are

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial f_{x}}{\partial\beta} = \bigl(E_{2}^{2}-E_{1}^{2}\bigr){\rm e}^{-\beta (E_{1}+E_{2})}-2E_{1}^{2}{\rm e}^{-\beta E_{1}}-E_{2}^{2}{\rm e}^{-\beta E_{2}},\label{eq:derivative of fx} \nonumber \end{align} \tag{ A.117a }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial f_{y}}{\partial\beta} = - \bigl(E_{2}^{2}-E_{1}^{2}\bigr){\rm e}^{-\beta (E_{1}+E_{2})} - E_{2}^{2}{\rm e}^{-\beta E_{2}}, \label{eq:derivative of fy} \nonumber \end{align} \tag{ A.117b }$$

which we insert into equation (A.116) to arrive at

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial^{2}y}{\partial x^{2}} = 6\,\mathcal{Z}^{3}f_{x}^{-3} E_{1}E_{2}\bigl(E_{2}-E_{1}\bigr){\rm e}^{-\beta(E_{1}+E_{2})}\,\geqslant\,0, \label{eq:2nd partial y wrt x rewritten final} \nonumber \end{align} \tag{ A.118 }$$

where we have made use of the fact that $f_{x}\geqslant0$ (since $1\geqslant {\rm e}^{-\beta E_{1}}$), along with $E_{1},E_{2},\mathcal{Z}\geqslant 0$ and $E_{2}\geqslant E_{1}$.

A.10.2. Explicit partial derivatives in dimension 4.

In this appendix, we present the explicit expressions for the first and second partial derivatives. The derivatives of the coordinates $x(\beta)$, $y(\beta)$ and $z(\beta)$ along the curve of thermal states can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial x}{\partial\beta} = - \mathcal{Z}^{-2}f_{x}(\beta), \label{eq:partial x wrt beta} \nonumber \end{align} \tag{ A.119a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial y}{\partial\beta} = - \mathcal{Z}^{-2}f_{y}(\beta), \label{eq:partial y wrt beta} \nonumber \end{align} \tag{ A.119b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial z}{\partial\beta} = - \mathcal{Z}^{-2}f_{z}(\beta), \label{eq:partial z wrt beta} \nonumber \end{align} \tag{ A.119c }$$

where the functions $f_{x}(\beta)$, $f_{y}(\beta)$, and $f_{z}(\beta)$ are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{x}(\beta) & = E_{1}{\rm e}^{-\beta E_{1}}(2+{\rm e}^{-\beta E_{2}}+{\rm e}^{-\beta E_{3}})\nonumber \\ & \ \ +(E_{2}{\rm e}^{-\beta E_{2}}+E_{3}{\rm e}^{-\beta E_{3}})(1-{\rm e}^{-\beta E_{1}})\,\geqslant\,0, \label{eq:fx} \nonumber \end{align} \tag{ A.120a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{y}(\beta) & = 3(E_{2}-E_{1}){\rm e}^{-\beta (E_{1}+E_{2})}+3E_{2}{\rm e}^{-\beta E_{2}}\nonumber \\ & \ \ +E_{3}{\rm e}^{-\beta E_{3}}(1+{\rm e}^{-\beta E_{1}}-2{\rm e}^{-\beta E_{2}})\nonumber \\ & \ \ +{\rm e}^{-\beta E_{3}}(2E_{2}{\rm e}^{-\beta E_{2}}-E_{1}{\rm e}^{-\beta E_{1}})\geqslant0,\label{eq:fy} \nonumber \end{align} \tag{ A.120b }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle f_{z}(\beta) & = 4{\rm e}^{-\beta E_{3}}\bigl[E_{3}+{\rm e}^{-\beta E_{1}}(E_{3}-E_{1})\nonumber \\ & \ \ \ +{\rm e}^{-\beta E_{2}}(E_{3}-E_{2})\bigr]\geqslant0, \nonumber \end{align} \tag{ A.120c }$$

where $f_{y}\geqslant0$ follows, since $E_{3}\geqslant E_{2}$ and $(1+{\rm e}^{-\beta E_{1}}-2{\rm e}^{-\beta E_{2}})\geqslant0$ which implies that the terms in f _y proportional to ${\rm e}^{-\beta E_{3}}$ can be bounded by $\newcommand{\bi}{\boldsymbol} E_{3}{\rm e}^{-\beta E_{3}}\bigl(E_{2}+{\rm e}^{-\beta E_{1}}(E_{2}-E_{1})\bigr)\geqslant0$. Since f _x, f _y, and f _z are nonnegative, all partial derivatives with respect to $\beta$ are nonpositive, but

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial y}{\partial x} = \left(\frac{\partial x}{\partial\beta}\right)^{-1}\frac{\partial y}{\partial\beta} = \frac{f_{y}(\beta)}{f_{x}(\beta)}\,\geqslant\,0, \label{eq:partial y wrt x} \nonumber \end{align} \tag{ A.121 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial z}{\partial x} = \left(\frac{\partial x}{\partial\beta}\right)^{-1}\frac{\partial z}{\partial\beta} = \frac{f_{z}(\beta)}{f_{x}(\beta)}\,\geqslant\,0, \label{eq:partial z wrt x} \nonumber \end{align} \tag{ A.122 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial z}{\partial y} = \left(\frac{\partial y}{\partial\beta}\right)^{-1}\frac{\partial z}{\partial\beta} = \frac{f_{z}(\beta)}{f_{y}(\beta)}\,\geqslant\,0. \label{eq:partial z wrt y} \nonumber \end{align} \tag{ A.123 }$$

For the second partial derivatives, we have

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial^{2}y}{\partial x^{2}} & = \left(\frac{\partial x}{\partial\beta}\right)^{-1}\frac{\partial}{\partial\beta}\,\left[\left(\frac{\partial x}{\partial\beta}\right)^{-1}\frac{\partial y}{\partial\beta}\right]\nonumber \\ & = \left(\frac{\partial x}{\partial\beta}\right)^{-2}\left\{\frac{\partial^{2}y}{\partial \beta^{2}}-\left(\frac{\partial y}{\partial \beta}\right)\left(\frac{\partial x}{\partial\beta}\right)^{-1}\left(\frac{\partial^{2}x}{\partial \beta^{2}}\right)\right\}\nonumber \\ & = \mathcal{Z}^{2}f_{x}^{{-3}}\bigl(\,f_{y}\frac{\partial f_{x}}{\partial\beta}-f_{x}\frac{\partial f_{y}}{\partial\beta}\bigr), \label{eq:2nd partial y wrt x alternative} \nonumber \end{align} \tag{ A.124 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial^{2}z}{\partial x^{2}} = \mathcal{Z}^{2}f_{x}^{-3}\bigl(\,f_{z}\frac{\partial f_{x}}{\partial\beta}-f_{x}\frac{\partial f_{z}}{\partial\beta}\bigr), \nonumber \end{align} \tag{ A.125 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial^{2}z}{\partial y^{2}} = \mathcal{Z}^{2}f_{y}^{-3}\bigl(\,f_{z}\frac{\partial f_{y}}{\partial\beta}-f_{y}\frac{\partial f_{z}}{\partial\beta}\bigr). \nonumber \end{align} \tag{ A.126 }$$

We thus have to calculate the derivatives of the functions f _x, f _y, and f _z with respect to $\beta$, which evaluate to

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial f_{x}}{\partial \beta} & = -E_{1}^{2}{\rm e}^{-\beta E_{1}}\bigl(2+{\rm e}^{-\beta E_{2}}+{\rm e}^{-\beta E_{3}}\bigr) \nonumber \\ & \,\,\, \ -\bigl(E_{2}^{2}{\rm e}^{-\beta E_{2}}+E_{3}^{2}{\rm e}^{-\beta E_{3}}\bigr)\bigl(1-{\rm e}^{-\beta E_{1}}\bigr), \label{eq:partial fx wrt beta} \nonumber \end{align} \tag{ A.127 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial f_{y}}{\partial \beta} & = E_{1}^{2}{\rm e}^{-\beta E_{1}} \bigl(3{\rm e}^{-\beta E_{2}}+{\rm e}^{-\beta E_{3}}\bigr)\nonumber \\ & \ \ \ -E_{2}^{2}{\rm e}^{-\beta E_{2}} \bigl(3+3{\rm e}^{-\beta E_{1}}+2{\rm e}^{-\beta E_{3}}\bigr)\nonumber \\ & \ \ \,\,\,-E_{3}^{2}{\rm e}^{-\beta E_{3}} \bigl(1+{\rm e}^{-\beta E_{1}}-2{\rm e}^{-\beta E_{2}}\bigr), \label{eq:partial fy wrt beta} \nonumber \end{align} \tag{ A.128 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial f_{z}}{\partial \beta} & = 4{\rm e}^{-\beta E_{3}} \bigl[E_{1}^{2}{\rm e}^{-\beta E_{1}}+E_{2}^{2}{\rm e}^{-\beta E_{2}} \nonumber \\ & \ \ \ - E_{3}^{2}\bigl(1+{\rm e}^{-\beta E_{1}}+{\rm e}^{-\beta E_{2}}\bigr)\bigr]. \label{eq:partial fz wrt beta} \nonumber \end{align} \tag{ A.129 }$$

Then, since $\mathcal{Z}^{2}\geqslant0$ and $f_{x},f_{y},f_{z}\geqslant0$, we just have to evaluate

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle f_{y}\frac{\partial f_{x}}{\partial\beta}-f_{x}\frac{\partial f_{y}}{\partial\beta} & = 2\bigl[E_{1}E_{2}{\rm e}^{-\beta(E_{1}+E_{2})}(E_{2}-E_{1}) + E_{1}E_{3}{\rm e}^{-\beta(E_{1}+E_{3})}(E_{3}-E_{1})\bigr]\bigl(3+{\rm e}^{-\beta E_{3}}\bigr)\mathcal{Z}\nonumber \\ & \ \ \ + 2E_{2}E_{3}{\rm e}^{-\beta(E_{2}+E_{3})}(E_{3}-E_{2})\bigl(1+{\rm e}^{-\beta E_{1}}\bigr)\bigl(2+2{\rm e}^{-\beta E_{1}}-{\rm e}^{-\beta E_{2}}+{\rm e}^{-\beta E_{3}}\bigr) \,\geqslant\,0, \nonumber \end{align} \tag{ A.130a }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle f_{z}\frac{\partial f_{x}}{\partial\beta}-f_{x}\frac{\partial f_{z}}{\partial\beta} & = 4\mathcal{Z}{\rm e}^{-\beta E_{3}}\,\bigl[E_{1}{\rm e}^{-\beta E_{1}}(E_{2}-E_{1}) \bigl({\rm e}^{-\beta E_{2}}(E_{3}-E-{2})+2E_{3}\bigr)\nonumber \\ & \ \ \ +(E_{3}-E_{2})E_{3}\bigl(E_{1}{\rm e}^{-\beta(E_{1}+E_{2})}+2E_{1}{\rm e}^{-\beta E_{1}}+E_{2}{\rm e}^{-\beta E_{2}}(1-{\rm e}^{-\beta E_{1}})\bigr) \bigr] \,\geqslant\,0, \nonumber \end{align} \tag{ A.130b }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle f_{z}\frac{\partial f_{y}}{\partial\beta}-f_{y}\frac{\partial f_{z}}{\partial\beta} = 12\mathcal{Z}{\rm e}^{-\beta (E_{2}+E_{3})}\bigl(E_{3}-E_{2}\bigr)\,\bigl[E_{1}^{2}{\rm e}^{-\beta E_{1}}+E_{3}\bigl(E_{2}-E_{1}\bigr){\rm e}^{-\beta E_{1}}+E_{2}\bigl(E_{3}-{\rm e}^{-\beta E_{1}}E_{1}\bigr) \bigr] \,\geqslant\,0. \nonumber \end{align} \tag{ A.130c }$$

For local dimension d  =  3, the transformations we consider here to ensure symmetric marginals imply that the vectors collecting the diagonal entries of the final state marginals are of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}}={\tilde{\mathbf {p}}}\SA={\tilde{\mathbf {p}}}\SB = M_q \mathbf{q} + (\mathbb{1}+ \Pi)M_{r} \mathbf{r}, \label{eq:EqualEvoMarginal d3} \nonumber \end{align} \tag{ A.131 }$$

where $\mathbf{q}=(q_{i})$ with $q_{i}=p_{i}^{2}$ and $\newcommand{\bi}{\boldsymbol} (\mathbb{1}+ \Pi)\mathbf{r}=\bigl(\,p_{i}(1-p_{i})\bigr)$. This means that the set of reachable states corresponds to the convex hull of the points $\tilde{\mathbf{p}}^{(i,j)}:=\Pi^{(i)}\mathbf{q}+(\mathbb{1}+ \Pi)\Pi^{(\,j)}\mathbf{r}$ for $i,j\in\{1,\ldots,3!\}$. However, in our case the problem can be further reduced to that of showing that the points $\mathbf{p}(\beta'<\beta)$ lie within a restricted part of the polytope, i.e. a triangle between the points $\mathbf{p}(\beta)$, $\mathbf{p}(\beta\rightarrow0)$, and $\tilde{\mathbf{p}}^{(6,1)}=\Pi^{(6)}\mathbf{q}+(\mathbb{1}+ \Pi)\Pi^{(1)}\mathbf{r}$, where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Pi^{(1)}=\mathbb{1},\ \ \Pi^{(6)} = \left(\begin{array}{@{}ccc@{}} 0 & 1 & 0 \nonumber \\ 1 & 0 & 0 \nonumber \\ 0 & 0 & 1 \end{array}\right),\ \ \Pi = \left(\begin{array}{@{}ccc@{}} 0 & 0 & 1 \nonumber \\ 1 & 0 & 0 \nonumber \\ 0 & 1 & 0 \end{array}\right). . \nonumber \end{align} \tag{ A.132 }$$

Note that the point ${\bf p}(\beta\rightarrow0)$ trivially lies in the polytope, since it can be obtained from an equally weighted convex combination of all cyclic permutations obtained for $\Pi^{(i)}=\Pi^{(\,j)}=\mathbb{1},\Pi,\Pi^{2}$. Then we switch from the coordinates $\{\,p_{0},p_{1},p_{2}\}$ to the coordinates $\{x,y,z\}$, given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x = - p_{0} + p_{1}, \nonumber \end{align} \tag{ A.133a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y = - p_{0} - p_{1} + 2p_{2}, \nonumber \end{align} \tag{ A.133b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z = -p_{0} - p_{1} - p_{2} = -1. \nonumber \end{align} \tag{ A.133c }$$

The plane containing the polytope is defined by z  =  −1 such that $p_{2}=1-p_{0}-p_{1}$, and we therefore have $y=2-3p_{0}-3p_{1}$. The corners $(\,p_{0},p_{1},p_{2})=(1,0,0)$, $(0,1,0)$, and $(0,0,1)$ of the outer triangle containing all (diagonal) qutrit states are mapped to $(x,y)=(-1,-1)$, $(1,-1)$, and $(0,2)$, respectively, see figure 1. Moreover, the curve of thermal states is parameterised as $\newcommand{\bi}{\boldsymbol} \bigl(x(\beta),y(\beta)\bigr)$ with $x(\beta)=-p_{0}(\beta)+p_{1}(\beta)$ and $y(\beta)=2-3p_{0}(\beta)-3p_{1}(\beta)$, and connects the points $(x,y)=(-1,-1)$ corresponding to $\mathbf{p}(\beta\rightarrow\infty)$ and $(x,y)=(0,0)$ corresponding to $\mathbf{p}(\beta\rightarrow0)$.

Next, we note that for any given thermal state with coordinates $\newcommand{\bi}{\boldsymbol} \bigl(x(\beta),y(\beta)\bigr)$, the coordinates $\newcommand{\bi}{\boldsymbol} \bigl(x=\tilde{x},y=\tilde{y}\bigr)$ of the corresponding point $\tilde{{\bf p}}^{(6,1)}$ are given by

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \tilde{x} = (\,p_{0}-p_{1})\bigl[2(\,p_{0}+p_{1})-1\bigr]\,\geqslant\,0, \nonumber \end{align} \tag{ A.134a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{y} = 2 - 3p_{0} - 3p_{1} = y(\beta), \nonumber \end{align} \tag{ A.134b }$$

and consequently this point always lies on the same height $\tilde{y}=y(\beta)$ as the point corresponding to ${\bf p}(\beta)$, but to the right, i.e. $\tilde{x}\geqslant0\geqslant x(\beta)$ of any of the points corresponding to thermal states (in particular $\mathbf{p}(\beta\rightarrow0)$). This shows that the point with coordinates $(x_{0},x_{1})=(0,y(\beta))$ is included in the polytope, since it can be obtained as a convex combination of ${\bf p}(\beta)$ and $\tilde{{\bf p}}^{(6,1)}$.

Moreover, evaluating the partial derivatives, we have found $(\partial y(\beta)/\partial x(\beta)\geqslant0)$ and $(\partial^{2} y/\partial x^{2}\geqslant0)$, see appendix A.10.2, i.e. the curve described by $\newcommand{\bi}{\boldsymbol} \bigl(x(\beta),y(\beta)\bigr)$ is a monotonically increasing $(\partial y(\beta)/\partial x(\beta)\geqslant0)$ and convex $(\partial^{2} y/\partial x^{2}\geqslant0)$ function, implying that the curve segment corresponding to all points $\newcommand{\bi}{\boldsymbol} \bigl(x(\beta'),y(\beta')\bigr)$ with $\beta'\leqslant\beta$ is contained within the triangle $\mathbf{p}(\beta)$, $\tilde{\mathbf{p}}^{(6,1)}$, and $\mathbf{p}(\beta\rightarrow0)$.

For local dimension d  =  4, we have marginal vectors of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\SB}{{_{B}}} \newcommand{\SA}{{_{A}}} \newcommand{\re}{{\rm Re}} \displaystyle {\tilde{\mathbf {p}}}={\tilde{\mathbf {p}}}\SA={\tilde{\mathbf {p}}}\SB = M_q \mathbf{q} + (\mathbb{1}+ \Pi)M_{r_{1}} \mathbf{r}_{1} + \frac{1}{2}(\mathbb{1}+ \Pi^{2})M_{r_{2}} \mathbf{r}_{2}. \label{eq:EqualEvoMarginal d4} \nonumber \end{align} \tag{ A.135 }$$

The set of reachable states corresponds to the convex hull of the points $\tilde{\mathbf{p}}^{(i,j,k)}:=\Pi^{(i)}\mathbf{q}+(\mathbb{1}+ \Pi)\Pi^{(\,j)}\mathbf{r}_{1}+\frac{1}{2}(\mathbb{1}+ \Pi^{2})\Pi^{(k)}\mathbf{r}_{2}$ for $i,j,k\in\{1,\ldots,4!\}$. Switching from the coordinates $\{\,p_{0},p_{1},p_{2},p_{3}\}$ to the new coordinates $\{x,y,z,w\}$, given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x = - p_{0} + p_{1}, \nonumber \end{align} \tag{ A.136a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y = - p_{0} - p_{1} + 2p_{2}, \nonumber \end{align} \tag{ A.136b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z = - p_{0} - p_{1} - p_{2} + 3p_{3}, \nonumber \end{align} \tag{ A.136c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w = - p_{0} - p_{1} - p_{2} - p_{3} = -1. \nonumber \end{align} \tag{ A.136d }$$

We here show that the points $(x,y,z)=(0,y(\beta),z(\beta))$ and $(0,0,z(\beta))$ are included. To achieve this, we identify a set of 5 points that all lie in the plane of constant $z=z(\beta)$, and show that the two-dimensional polytope spanned by these points encloses both $(0,y(\beta),z(\beta))$ and $(0,0,z(\beta))$, see figure A1. These five points are $A=\tilde{\mathbf{p}}^{(1,1,1)}$, $B=\tilde{\mathbf{p}}^{(7,1,1)}$, $C=\tilde{\mathbf{p}}^{(13,7,1)}$, $D=\tilde{\mathbf{p}}^{(9,7,1)}$, and $E=\tilde{\mathbf{p}}^{(9,1,1)}$, where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Pi & = \left(\begin{array}{@{}cccc@{}} 0 & 0 & 0 & 1 \nonumber \\ 1 & 0 & 0 & 0 \nonumber \\ 0 & 1 & 0 & 0 \nonumber \\ 0 & 0 & 1 & 0 \end{array}\right),\ \Pi^{(7)} = \left(\begin{array}{@{}cccc@{}} 0 & 1 & 0 & 0 \nonumber \\ 1 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 1 & 0 \nonumber \\ 0 & 0 & 0 & 1 \end{array}\right),\nonumber \\ \Pi^{(9)} & = \left(\begin{array}{@{}cccc@{}} 0 & 1 & 0 & 0 \nonumber \\ 0 & 0 & 1 & 0 \nonumber \\ 1 & 0 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 1 \end{array}\right),\ \Pi^{(13)} = \left(\begin{array}{@{}cccc@{}} 0 & 0 & 1 & 0 \nonumber \\ 1 & 0 & 0 & 0 \nonumber \\ 0 & 1 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 1 \end{array}\right). \nonumber \end{align} \tag{ A.137 }$$

Zoom In Zoom Out Reset image size

The point B can be seen to have the same y -coordinate as the starting point A, $y(B)=y(A)=y(\beta)$, while

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle x(B) = \bigl(\,p_{0}(\beta) - p_{1}(\beta)\bigr)\bigl(2p_{0}(\beta) + 2p_{1}(\beta)-1\bigr)\,\geqslant\,0, \nonumber \end{align} \tag{ A.138 }$$

since $p_{0}(\beta)\geqslant p_{1}(\beta)$ and $p_{0}(\beta)+p_{1}(\beta)\geqslant\frac{1}{2}$ for any Hamiltonian and any temperature. This means the point $(0,y(\beta),z(\beta))$ is contained in polytope since it can be obtained as a convex combination of A and B.

Finally, to show that the point $\tilde{O}=(0,0,z(\beta))$ is contained within the polytope, we show that the points C, D and E satisfy, $x(C)\geqslant0$, $y(D)\geqslant0$ and $x(E)\leqslant 0$, while $(A\tilde{O}\times AE)_{z}\geqslant0$, $(E\tilde{O}\times ED)_{z}\geqslant0$ and $(D\tilde{O}\times DC)_{z}\geqslant0$, where for any two points X and Y we use the notation XY  =  Y  −  X. In other words, the closed path connecting the five points A, B, C, D and E encircles the point $\tilde{O}$, meaning that $\tilde{O}$ can be written as a convex combination of these points.

In a slight abuse of notation, we now use the shorthand $\newcommand{\e}{{\rm e}} p_{i}\equiv p_{i}(\beta)$ for the thermal state diagonal components to express the relevant coordinates for the points C, D and E as

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle x(C) & = \bigl(\,p_{0} - p_{1}\bigr) \bigl(2 p_{0} + 2p_{1} + p_{2} - 1\bigr)\nonumber \\ & \ \,\,\,\,\, + \bigl(\,p_{1} - p_{2}\bigr) \bigl(\,p_{0} + p_{1} + p_{2}\bigr)\,\geqslant\,0, \nonumber \end{align} \tag{ A.139a }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle y(D) & = 3\bigl(\,p_{0} - p_{1}\bigr) \bigl(\,p_{0} + p_{1} + p_{2} - \frac{1}{3}\bigr)\nonumber \\ & \ \,\,\,\, + 3\bigl(\,p_{1} - p_{2}\bigr) \bigl(\,p_{0} + p_{1} + p_{2}-\frac{2}{3}\bigr)\,\geqslant\,0, \nonumber \end{align} \tag{ A.139b }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle x(E) & = -\bigl[\bigl(\,p_{0} - p_{1}\bigr)\bigl(1 - p_{0} - p_{1}\bigr)\nonumber \\ & \ \ \ \ + \bigl(\,p_{1} - p_{2}\bigr) \bigl(\,p_{1} + p_{2}\bigr)\bigr]\,\leqslant\,0, \nonumber \end{align} \tag{ A.139c }$$

while the z-components of the relevant cross products are given explicitly below. In appendix A.10.2 we have further shown that the partial derivatives $\partial y/\partial x$, $\partial z/\partial x$, and $\partial z/\partial y$ as well as second derivatives $\partial^{2} y/\partial x^{2}$, $\partial^{2} z/\partial x^{2}$, and $\partial^{2} z/\partial y^{2}$ are nonnegative along the curve of thermal states, meaning that condition (I) and (II) are satisfied for d  =  4.

A.12.1. Relevant cross products.

Finally, the z-components of the relevant cross products are given

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle (A\tilde{O}\times AE)_{z} = & \bigl(\,p_{0} - p_{1}\bigr)^{2} \bigl(2 p_{0} - 2p_{1} + 3p_{2}\bigr) \nonumber \\ & + \bigl(\,p_{0} - p_{1}\bigr)\bigl(\,p_{1} - p_{2}\bigr) \bigl(\,p_{0} - p_{1} + 4p_{2}\bigr) + 2\bigl(\,p_{1} - p_{2}\bigr)^{2} \bigl(\,p_{1} + p_{2}\bigr) \,\geqslant\,0, \nonumber \end{align} \tag{ A.140a }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle (E\tilde{O}\times ED)_{z} = 2p_{1}\bigl(\,p_{0}-p_{2}\bigr) \Bigl[\bigl(\,p_{0}-p_{1}\bigr) \bigl(1-p_{1}+p_{2}\bigr) + \bigl(\,p_{1}-p_{2}\bigr)\bigl(\,p_{1}+2p_{2}-p_{3}\bigr)\Bigr]\,\geqslant\,0, \nonumber \end{align} \tag{ A.140b }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle (D\tilde{O}\times DC)_{z} & = \bigl(\,p_{0} - p_{1}\bigr)^{2} \bigl(\,p_{0} + p_{2}\bigr)\bigl((3 p_{0} + 3 p_{1} - 1) + 3 (\,p_{1} - p_{3})\bigr) \nonumber \\ & \ \ + \bigl(\,p_{0} - p_{1}\bigr)\bigl(\,p_{1} - p_{2}\bigr) \bigl(3 (\,p_{0} - p_{1}) (\,p_{0} + p_{1}) + 6 p_{0} (\,p_{2} - p_{3}) + 2 p_{2} (2 p_{0} + 5 p_{1} + 2 p_{2} - p_{3})\bigr) \nonumber \\ &\ \ + \bigl(\,p_{1} - p_{2}\bigr)^{2} \bigl(\bigl[\,p_{0} - 3 p_{3} (1 - p_{3}) + 6 p_{1} p_{2}\bigr] + (\,p_{1} - p_{2}) + 3 p_{2}^{2} + 2 (\,p_{0} - p_{1}) (1 + 3 p_{0})\bigr) \,\geqslant\,0, \nonumber \end{align} \tag{ A.140c }$$

where the last inequality follows since the term in angled brackets in the last line can be shown to nonnegative using $6p_{1}p_{2}\geqslant 6p_{3}^{2}$ and $3 p_{3} (1 - 3p_{3})\leqslant \frac{1}{4}\leqslant p_{0}$.

Here we present a possible recursive generalisation of the geometric approach, based on re-expressing condition (I) in terms of majorisation relations rather than trying to prove it in a full geometric way.

To illustrate this method we refer back to the proof of theorem 3 for the case d  =  4. We have seen that the point $\newcommand{\re}{{\rm Re}} \mathbf {v}_1 =(0,x_1(\beta),x_2(\beta),-1)$ can be reached with a transformation of the type $\newcommand{\re}{{\rm Re}} \mathbf {v}_1 =M_{q_1}{\bf q}+(\mathbb {1}+\Pi)\mathbb {1}{\bf r}_1 +\frac{1}{2} (\mathbb {1}+\Pi^2)\mathbb {1} {\bf r}_2$, where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{q_1}=\left(\begin{array}{@{}cccc@{}} m & 1-m & 0 & 0 \nonumber \\ 1-m & m & 0 & 0 \nonumber \\ 0 & 0 & 1 & 0 \nonumber \\ 0 & 0 & 0 & 1 \end{array}\right) , \nonumber \end{align} \tag{ A.141 }$$

with $m=1-1/2(\,p_0+p_1)$, which is doubly-stochastic since $p_0+p_1\geqslant 1/2$ holds for d  =  4.

To prove that $\newcommand{\re}{{\rm Re}} \mathbf {v}_2=(0,0,x_2(\beta),-1)$ can be reached, we also considered $M_{r_2}=\mathbb {1}$, which turns out to be a big simplification, due to the following reasoning. Assuming $M_{r_2}=\mathbb {1}$ we can rewrite the general transformation that (potentially) reaches $\newcommand{\re}{{\rm Re}} \mathbf {v}_2$ as

$$\begin{align} \newcommand{\e}{{\rm e}} mathbf \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:geomproof1} \mathbf {v}_2+{\bf q} - {\bf p} - (\mathbb {1}+\Pi)(M_{r_1}-\mathbb {1}) {\bf r}_1 = M_q {\bf q} , \nonumber \end{align} \tag{ A.142 }$$

where we used the fact that $\frac 1 2 (\mathbb {1}+\Pi^2){\bf r}_2 = {\bf p}-{\bf q} -(\mathbb {1}+\Pi) {\bf r}_1$. Thus, with this simplification of assuming $M_{r_2}=\mathbb {1}$, the problem reduces to showing that there exist some pair $(M_{q},M_{r_1})$ such that equation (A.142) holds. In particular, we use the fact that $(1+\Pi)(M_{r_1}-\mathbb {1})$ has norm equal to zero, which in turn implies that the sum of vectors on the left hand side of equation (A.142) has the same norm as ${\bf q}$, since both $\newcommand{\re}{{\rm Re}} \mathbf {v}_2$ and ${\bf p}$ have norm 1 (i.e. they are probabilities). Note also that the vector $(\mathbb {1}+\Pi)(M_{r_1}-\mathbb {1})$ by itself need not be ordered with decreasing components, as well as the whole expression in the left hand side of equation (A.142). More precisely, in this case the greatest component of expression $\newcommand{\re}{{\rm Re}} \mathbf {v}_2+{\bf q} - {\bf p}$ is the third.

To prove the statement, then, we can show that $\newcommand{\re}{{\rm Re}} \mathbf {v}_2+{\bf q} - {\bf p} - (\mathbb {1}+\Pi)(M_{r_1}-\mathbb {1}) {\bf r}_1$ is majorised by ${\bf q}$, or, more precisely, we have to find a proper stochastic matrix $M_{r_1}$ such that the above majorisation holds. This is achieved again with a matrix of the form $M_{r_1}=M_{q_1}$ that leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:geomproof2} (\mathbb {1}+\Pi)(M_{r_1}-\mathbb {1}) {\bf r}_1=(-a,0,a,0) , \nonumber \end{align} \tag{ A.143 }$$

with $a=(m-1)\,p_1(\,p_0-p_2)$. Then, the first condition for the majorisation relation to hold is that the third component of the left hand side of equation (A.142) is smaller than the first component of ${\bf q}$, i.e. that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{p_0+p_1-2p_2}{3}+p_2^2+(m-1)\,p_1(\,p_0-p_2) \leqslant p_0^2 , \nonumber \end{align} \tag{ A.144 }$$

which implies that we have to choose

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1-m= \max \{0, \frac{p_0+p_1-2p_2-3(\,p_0^2-p_2^2)}{3p_1(\,p_0-p_2)} \} , \nonumber \end{align} \tag{ A.145 }$$

such that the majorisation relation is correctly satisfied, since the other conditions are also trivially satisfied.

The last step to check is that with this choice of the parameter 1  −  m, $M_{r_1}$ is indeed a doubly-stochastic matrix. This is true because

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:1md4} 1-m \leqslant 1 \iff p_0+p_1-2p_2-3(\,p_0^2-p_2^2) \leqslant 3p_1(\,p_0-p_2) \nonumber \end{align} \tag{ A.146 }$$

also holds, as well as $1-m\geqslant 0$, which holds by construction.

For higher dimensions, we can try a sort of recursive reasoning, that we work out here for the case d  =  5. Let us hypotesise that the three vertices $\newcommand{\re}{{\rm Re}} (\mathbf {v}_1, \mathbf {v}_2, \mathbf {v}_{3})$ can be found with a set of transformations $\{((M_q)_i, (M_{r_1})_i , \mathbb {1})\}_{i=1,2,3}$, i.e. that are such that the last stochastic matrix is always the identity $(M_{r_{2}})_i=\mathbb {1}$ for all $i=1,2,3$. This allows us to decompose the last vector as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\mathbb{1}+ \Pi^2){\bf r}_2 = {\bf p}- {\bf q} - (\mathbb{1}+ \Pi){\bf r}_1 , \nonumber \end{align} \tag{ A.147 }$$

that is ${\bf p} - {\bf q}$ minus a sum of terms that have zero norm. Then, we rewrite the condition for a vertex $v_i$ to be reached as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:geomproof3} \mathbf {v}_{i}- {\bf p} + {\bf q} - (\mathbb{1}+ \Pi)((M_{r_1})_i-\mathbb {1}) {\bf r}_1 = (M_q)_i{\bf q} , \nonumber \end{align} \tag{ A.148 }$$

which is satisfied if and only if the left hand side is majorised by ${\bf q}$. So now, the task would be to find suitable doubly-stochastic matrices $(M_{r_1})_i$ such that the left hand side is majorised by ${\bf q}$ for each of the $\newcommand{\re}{{\rm Re}} \mathbf {v}_i$. Let us then start by considering $\newcommand{\re}{{\rm Re}} \mathbf {v}_1$, which in the probability coordinates is $\newcommand{\re}{{\rm Re}} \mathbf {v}_1=((\,p_0+p_1)/2,(\,p_0+p_1)/2,p_2,p_3,p_4)$. In this case, the greatest component of $\newcommand{\re}{{\rm Re}} \mathbf {v}_{i}- {\bf p} + {\bf q}$ is the second, namely $\frac{p_0-p_1} 2 +p_1^2$. Unlike lower dimensions, however, in this case we have that $\frac{p_0-p_1} 2 +p_1^2\leqslant p_0^2 \iff p_0+p_1 \geqslant 1/2$ does not always hold, meaning that this time $(M_{r_1})_1=\mathbb {1}$ does not suffice. Thus, in general we have to find $(M_{r_1})_1$ such that the second component of the vector $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {r}}}_1:=(\mathbb{1}+ \Pi)((M_{r_1})_1-\mathbb {1}) {\bf r}_1$ is negative. The simplest would be if $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {r}}}_1=(-a,a,0,0,0)$ for some positive a. However, it can be readily checked that this is not possible. Nevertheless, we can look for vectors, for example, of the form $\newcommand{\re}{{\rm Re}} {\tilde{\mathbf {r}}}_1=(a_1,a_2,-a_1-a_2,0,0)$, for positive a₂ such that we can satisfy $\frac{p_0-p_1} 2 +p_1^2-a_2\leqslant p_0^2$ by choosing

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle a_2 = \max \{0, (\,p_0-p_1)(\frac 1 2-p_0-p_1)\} . \nonumber \end{align} \tag{ A.149 }$$

Afterwards, the next condition for the majorisation to hold is either

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p_0^2 +p_1^2 - a_1 -a_2 \leqslant p_0^2+p_1^2 , \nonumber \end{align} \tag{ A.150 }$$

or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:lastappcond3} \frac{p_0-p_1} 2 +p_1^2+p_2^2+a_1\leqslant p_0^2+p_1^2 . \nonumber \end{align} \tag{ A.151 }$$

While the first is trivially satisfied if a₁ and a₂ are both positive, the second needs to be checked. After this, all other conditions for the majorisation will be trivially satisfied. This can be achieved with the following matrix:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (M_{r_1})_1 =\left(\begin{array}{@{}ccccc@{}} 1-m_1 & m_1 & 0 & 0 & 0 \nonumber \\ 0 & 1-m_1 & m_1 & 0 & 0 \nonumber \\ m_1 & 0 & 1-m_1 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 1-m_2 & m_2 \nonumber \\ 0 & 0 & 0 & m_2 & 1-m_2 \end{array}\right) , \nonumber \end{align} \tag{ A.152 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle a_2&=m_2p_4(\,p_0-p_3) = m_1(\,p_0p_1-p_2p_3) , \nonumber \\ a_1&=m_1p_2(\,p_1-p_3) , \nonumber \end{align} \tag{ A.153 }$$

and we have to ensure the following conditions (only in the case $1/2\geqslant p_0+p_1$), coming from $m_1\leqslant 1$ and $m_2\leqslant 1$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:1md5II} (\,p_0-p_1)(\frac 1 2-p_0-p_1)/(\,p_0p_1-p_2p_3) \leqslant 1 , \nonumber \\ (\,p_0-p_1)(\frac 1 2-p_0-p_1)/(\,p_4(\,p_0-p_3)) \leqslant 1 , \nonumber \end{align} \tag{ A.154 }$$

plus we have to guarantee that equation (A.151) is satisfied, which leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:lastappcond31} (\,p_0-p_1)(\frac 1 2-p_0-p_1)\,p_2(\,p_1-p_3)/(\,p_0p_1-p_2p_3) \leqslant p_0^2-p_2^2-\frac{p_0-p_1} 2 , \nonumber \end{align} \tag{ A.155 }$$

since $a_1=(\,p_0-p_1)(\frac 1 2-p_0-p_1)\,p_2(\,p_1-p_3)/(\,p_0p_1-p_2p_3)$. The above conditions are nontrivial, but still satisfied, as it is shown below in appendix A.13.1.

Let us now come to $\newcommand{\re}{{\rm Re}} \mathbf {v}_2$. For this case we can try to re-use the result of d  =  4, and show that the majorisation relation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq:geomproof4} \mathbf {v}_{2}- {\bf p} + {\bf q} + (\mathbb{1}+ \Pi)((M_{r_1})_2-\mathbb {1}) {\bf r}_1 = M_{q_2}{\bf q} , \nonumber \end{align} \tag{ A.156 }$$

holds for

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (M_{r_1})_2=\left(\begin{array}{@{}ccccc@{}} m & 1-m & 0 & 0 & 0 \nonumber \\ 1-m & m & 0 & 0 & 0 \nonumber \\ 0 & 0 & 1 & 0 & 0 \nonumber \\ 0 & 0 & 0 & 1 & 0 \nonumber \\ 0 & 0 & 0 & 0 & 1 \end{array}\right) , \nonumber \end{align} \tag{ A.157 }$$

with the appropriate coefficient 1  −  m, which is exactly the same as in the previous case, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1-m= \max \{0, \frac{p_0+p_1-2p_2-3(\,p_0^2-p_2^2)}{3p_1(\,p_0-p_2)} \} , \nonumber \end{align} \tag{ A.158 }$$

and again we have to ensure that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:reld5d4} 1-m \leqslant 1 \iff p_0+p_1-2p_2-3(\,p_0^2-p_2^2) \leqslant 3p_1(\,p_0-p_2) . \nonumber \end{align} \tag{ A.159 }$$

However, while equation (A.159) holds alwasy true for d  =  4, this is no longer the case for d  =  5. Still, there is a range of parameters $\{E_i\}_{i=1}^5$ where it holds (see also appendix A.13.1 below).

Last, for the point $\newcommand{\re}{{\rm Re}} \mathbf {v}_{2}$ we can consider a matrix similar to the above, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (M_{r_1})_3=\left(\begin{array}{@{}ccccc@{}} 1 & 0 & 0 & 0 & 0 \nonumber \\ 0 & m & 1-m & 0 & 0 \nonumber \\ 0 & 1-m & m & 0 & 0 \nonumber \\ 0 & 0 & 0 & 1 & 0 \nonumber \\ 0 & 0 & 0 & 0 & 1 \end{array}\right) , \nonumber \end{align} \tag{ A.160 }$$

but now the condition to be satisfied is $\frac{p_0+p_1+p_2+p_3} 4 -p_3+p_3^2-(1-m)\,p_2(\,p_1-p_3)\leqslant p_0^2$ which leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1-m = \max \{0, \frac{p_0+p_1+p_2-3p_3-4(\,p_0^2-p_3^2)}{4p_2(\,p_1-p_3)} \} . \nonumber \end{align} \tag{ A.161 }$$

The requirement $1-m \leqslant 1$ thus translates to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:proofd5last} p_0+p_1+p_2-3p_3-4(\,p_0^2-p_3^2) \leqslant 4p_2(\,p_1-p_3) . \nonumber \end{align} \tag{ A.162 }$$

While this last relation might not hold in general, it can be satisfied in some appropriately chosen range of the parameters $\{E_i\}_{i=1}^4$. Thus, in order to fully prove the d  =  5 case with this method, one would need to find other possible $(M_{r_1})_2$ and $(M_{r_1})_3$ for the respective complementary parameter regions. Nevertheless, this can possibly still be achieved while keeping $M_{r_2}=\mathbb {1}$.

In conclusion, we have seen a possible strategy to prove that all the vertices listed in Condition (I) can be reached with the transformations preserving equal marginals. This exploits the geometry of such achievable thermal marginals and tries to avoid the difficulties encountered in the 'passing on the norm' approach, by considering the assumption that it is always possible to choose one of the matrices $M_{r_i}$ as the identity. This assumption is in fact supported by the geometric proofs in dimensions d  =  3 and d  =  4.

A.13.1. Validity of doubly-stochasticity of ${M_{r_i}}$ for d  =  5.

Here we study the validity of the relations that guarantee that the matrices $M_{r_i}$ needed for the proof of d  =  5 are indeed doubly-stochastic. We start with equation (A.146), which was also needed for the d  =  4 case. To prove equation (A.146) we can see that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p_0+p_1-2p_2-3(\,p_0^2-p_2^2)\leqslant 3p_1(\,p_0-p_2) . \nonumber \end{align} \tag{ A.163 }$$

The above relation in fact holds since we can rewrite it as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\,p_0-p_2)(1-3(\,p_0+p_1+p_2))+(\,p_1-p_2)\leqslant 0 , \nonumber \end{align} \tag{ A.164 }$$

and we can see that indeed

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\,p_0-p_2)(1-3(\,p_0+p_1+p_2))+(\,p_1-p_2) \nonumber \\ \leqslant (\,p_0-p_2)(2-3(\,p_0+p_1+p_2)) \leqslant 0 , \nonumber \end{align} \tag{ A.165 }$$

since $p_0+p_1+p_2 \geqslant 2/3$ holds for d  =  4, but not always d  =  5. In fact, for d  =  5 we can have points close to $p_0=p_1=p_2=1/5$.

To prove equation (A.154), we have to prove that, whenever $\frac 1 2 \geqslant p_0+p_1$ it holds that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\,p_0-p_1)(\frac 1 2-p_0-p_1) \leqslant p_4(\,p_0-p_3) \leqslant p_0p_1-p_2p_3 . \nonumber \end{align} \tag{ A.166 }$$

This is indeed true, since we can rewrite it as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\,p_0-p_1)(\frac 1 2-p_0-p_1) - p_4(\,p_0-p_3) \leqslant (\,p_0-p_1)(\frac 1 2-p_0-p_1-p_4) \leqslant 0 , \nonumber \end{align} \tag{ A.167 }$$

where the last relation holds whenever $\frac 1 2 \geqslant p_0+p_1$ because we can write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac 1 2-p_0-p_1-p_4 = p_2+p_3 - \frac 1 2 \leqslant p_0+p_1 - \frac 1 2 \leqslant 0 . \nonumber \end{align} \tag{ A.168 }$$

Let us know consider equation (A.155). Since $p_2(\,p_1-p_3)\leqslant (\,p_0p_1-p_2p_3)$ we can rewrite it as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\,p_0-p_1)(\frac 1 2-p_0-p_1) \leqslant p_0^2-p_2^2-\frac{p_0-p_1} 2 , \nonumber \end{align} \tag{ A.169 }$$

which is always satisfied since $p_0-p_1 \leqslant 2(\,p_0^2-p_1^2) \leqslant 2p_0^2-p_2^2-p_1^2$.

Finally, let us consider equation (A.162). We can rewrite the expression as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\,p_0-p_3)(1-4(\,p_0+p_3))+(\,p_1-p_3)(1-4p_2)+p_2-p_3 \leqslant 0 , \nonumber \end{align} \tag{ A.170 }$$

and again, we can see that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & (\,p_0-p_3)(3-4(\,p_0+p_2+p_3)) \nonumber \\ & \leqslant (\,p_0-p_3)(1-4(\,p_0+p_3))+(\,p_1-p_3)(1-4p_2)+p_2-p_3 \nonumber \end{align} \tag{ A.171 }$$

and in particular both the above expressions are surely negative whenever

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p_0+p_2+p_3 \geqslant \frac{3}{4}, \nonumber \end{align} \tag{ A.172 }$$

which is not always true, but it is so in a certain range of parameters.