Thermodynamic cost of creating correlations

3.1.1. Bipartite systems

Let us start by considering the case of a two qudit system, i.e. two d-level systems. Correlations between the two systems (which shall be referred to as A₁ and A₂) can naturally be measured using the quantum mutual information $I({{A}_{1}}:{{A}_{2}})$

$$\begin{eqnarray}&&I({{A}_{1}}:{{A}_{2}})=S({{A}_{1}})+S({{A}_{2}})-S({{A}_{1}}{{A}_{2}}),\end{eqnarray} \tag{ 4 }$$

where $S(X)=-{\rm Tr}({{\rho }_{X}}{\rm log} {{\rho }_{X}})$ is the von Neumann entropy of system X.

The goal is then to find the the optimal unitary operation U such that ${{\rho }_{{\rm f}}}=U{{\tau }_{\beta }}\otimes {{\tau }_{\beta }}{{U}^{\dagger }}$ has the maximal possible mutual information. Note first that initially $I({{A}_{1}}:{{A}_{2}})=0$, as the initial state factorizes. Thus, to create correlations, one must find a global unitary that increases the local entropies $S({{A}_{i}})$ of ${{\rho }_{{\rm f}}}$, since the total entropy $S({{A}_{1}}{{A}_{2}})=2S({{\tau }_{\beta }})$ cannot change. Since for a d-level system the local entropy is upper bounded by $S({{A}_{i}})\leqslant {\rm log} d$, the maximal possible mutual information is upper bounded by

$$\begin{eqnarray}&&I({{A}_{1}}:{{A}_{2}})\leqslant 2[{\rm log} d-S({{\tau }_{\beta }})].\end{eqnarray} \tag{ 5 }$$

This bound can always be achieved, by making use of the following protocol, which amounts to rotating from the energy eigenbasis to the generalized Bell basis, i.e. to a basis of maximally entangled qudit states. In more detail, for all d one can define the unitary operators

$$\begin{eqnarray}&&X=\mathop{\sum }\limits_{m}|m+1{\rm mod} d\rangle \langle m|,\qquad Z=\mathop{\sum }\limits_{m}{{\omega }^{m}}|m\rangle \langle m|,\end{eqnarray} \tag{ 6 }$$

with $\omega ={{{\rm e}}^{2\pi {\rm i}/d}}$ as generalizations of the (qubit) Pauli operators ${{\sigma }_{x}}$ and ${{\sigma }_{z}}$. The Bell basis $\{|{{\phi }_{ij}}\rangle {{\}}_{ij}}$ is then given by

$$\begin{eqnarray}&&|{{\phi }_{ij}}\rangle ={{Z}^{i}}\otimes {{X}^{j}}|\phi \rangle ,\end{eqnarray} \tag{ 7 }$$

where $|\phi \rangle =\frac{1}{\sqrt{d}}\mathop{\sum }\limits_{i}|ii\rangle$. Finally, we consider the operation given by

$$\begin{eqnarray}&&U=\mathop{\sum }\limits_{ij}|{{\phi }_{ij}}\rangle \langle ij|.\end{eqnarray} \tag{ 8 }$$

Since the initial state is a mixture of energy eigenstates, ${{\rho }_{{\rm f}}}$ is a mixture of Bell states. Finally, since these all have maximally mixed marginals, i.e. ${\rm T}{{{\rm r}}_{{{A}_{k}}}}(|{{\phi }_{ij}}\rangle \langle {{\phi }_{ij}}|)={\mathbb{1}}/d$, the bound (5) is achieved [28]. We end by noting that the maximally mixed state ${\mathbb{1}}/d$ corresponds to the infinite-temperature thermal state ${{\tau }_{0}}$. We shall see in the second half of the paper that when one has in addition a constraint on the energy, the optimal protocol produces thermal marginals, only there at lower temperatures.

Finally, we note that for all finite initial temperatures $\beta \ne 0$ the mutual information that can be created between the two systems is non-zero, i.e. that one can produce correlations between them at arbitrary finite temperatures.

3.1.2. Multipartite systems

In the multipartite setting one can generalize the notion of mutual information by considering the difference between the sum of local entropies and the total entropy of the global system. That is, for a collection of n systems ${{A}_{1}},\ldots ,{{A}_{n}}$, we define the multipartite mutual information as

$$\begin{eqnarray}&&I({{A}_{i}}:\cdots :{{A}_{n}})=\mathop{\sum }\limits_{i=1}^{n}S({{A}_{i}})-S({{A}_{1}}\cdots {{A}_{n}}),\end{eqnarray} \tag{ 9 }$$

which vanishes only when the global system is a product state. Again, since the total entropy is conserved, to maximize this quantity one must maximize the sum of final local entropies after the protocol. The analogous upper bound

$$\begin{eqnarray}&&I({{A}_{i}}:\cdots :{{A}_{n}})\leqslant n({\rm log} d-S({{\tau }_{\beta }})),\end{eqnarray} \tag{ 10 }$$

is seen to hold, and can again be achieved by rotating the energy eigenbasis to a basis of generalized GHZ states. Namely, one can define the basis $\{|\phi _{{{i}_{1}}\cdots {{i}_{n}}}^{n}\rangle {{\}}_{{{i}_{1}}\cdots {{i}_{n}}}}$ by

$$\begin{eqnarray}&&|\phi _{{{i}_{1}}\cdots {{i}_{n}}}^{n}\rangle ={{Z}^{{{i}_{1}}}}\otimes {{X}^{{{i}_{2}}}}\otimes \cdots \otimes {{X}^{{{i}_{n}}}}|{{\phi }^{n}}\rangle ,\end{eqnarray} \tag{ 11 }$$

where $|{{\phi }^{n}}\rangle =\frac{1}{\sqrt{d}}\mathop{\sum }\limits_{i}|i{{\rangle }^{\otimes n}}$ and the operation U given by

$$\begin{eqnarray}&&U=\mathop{\sum }\limits_{{{i}_{i}},\ldots ,{{i}_{n}}}|{{\phi }_{{{i}_{1}}\cdots {{i}_{n}}}}\rangle \langle {{i}_{1}}\cdots {{i}_{n}}|.\end{eqnarray} \tag{ 12 }$$

Again, since the final state of the global system is a mixture of generalized GHZ states, all of which have maximally mixed marginals ${\rm T}{{{\rm r}}_{{{{\bar{A}}}_{k}}}}(|{{\phi }_{ij}}\rangle \langle {{\phi }_{ij}}|)={\mathbb{1}}/d$ (where ${{\bar{A}}_{k}}$ denotes tracing over all systems except A_k) the bound is seen to be saturated. Finally, as long as the initial temperature is not infinite $\beta \ne 0$, then the bound is non-zero, and a finite amount of correlation can be created.

Having seen in the previous section that it is possible to create correlations in a multipartite thermal state starting at arbitrary temperatures in a mathematically easy fashion, we now move on to the more interesting question of creating entanglement. We will first look at the case of bipartite systems, where there is a single notion of entanglement, before moving on to multipartite systems, where there are a number of inequivalent notions of entanglement that we will study. In all cases we will restrict ourselves to the study of qubits.

3.2.1. Bipartite systems

We shall start our study of the bipartite case with the simplest possible scenario, involving two qubits. Although there is only a single notion of entanglement, one can nevertheless define many inequivalent measures of entanglement. Here for concreteness we will focus on the concurrence [29], which for pure states is the linear entropy of the reduced state of one party

$$\begin{eqnarray}&&C(\psi )=\sqrt{2(1-{\rm tr}(\rho _{A}^{2}))},\end{eqnarray} \tag{ 13 }$$

where ${{\rho }_{A}}={\rm t}{{{\rm r}}_{{\rm B}}}|\psi \rangle \langle \psi {{|}_{A{\rm B}}}$, and is extended to mixed states via the convex-roof construction

$$\begin{eqnarray}&&C(\rho )={\rm inf} \mathop{\sum }\limits_{i}{{p}_{i}}C({{\psi }_{i}}),\end{eqnarray} \tag{ 14 }$$

where the infimum is taken over all pure state decompositions $\rho ={{\sum }_{i}}{{p}_{i}}|{{\psi }_{i}}\rangle \langle {{\psi }_{i}}|$. The concurrence is important as for qubits the convex roof can be analytically calculated and the entanglement of formation [29] is functionally related to it.

Crucially, for our purposes the problem of finding the state of maximal concurrence given only its spectrum was solved in [30, 31], which is an alternative way of phrasing the problem which we are interested in here. Moreover, it was shown that the optimal protocol not only maximizes the concurrence (and therefore the entanglement of formation), but also two other important measures of entanglement, the relative entropy of entanglement, and the negativity.

The protocol of [31] is easiest understood by decomposing it into a product of two unitaries, $U={{V}_{2}}{{V}_{1}}$, where V₁ is a CNOT gate

$$\begin{eqnarray}&&{{V}_{1}}=|00\rangle \langle 00|+|01\rangle \langle 01|+|11\rangle \langle 10|+|10\rangle \langle 11|,\end{eqnarray} \tag{ 15 }$$

and V₂ is a rotation in the subspace spanned by $\{|00\rangle ,|11\rangle \}$ to maximally entangled states

$$\begin{eqnarray}&&{{V}_{2}}=|{{\phi }_{00}}\rangle \langle 00|+|01\rangle \langle 01|+|10\rangle \langle 10|+|{{\phi }_{10}}\rangle \langle 11|.\end{eqnarray} \tag{ 16 }$$

Denoting by ${{\{{{\lambda }_{i}}\}}_{i}}$ the eigenvalues of the initial state ${{\rho }_{{\rm i}}}$ arranged in non-increasing order, the concurrence of the final state ${{\rho }_{{\rm f}}}={{V}_{2}}{{V}_{1}}{{\rho }_{{\rm i}}}V_{1}^{\dagger }V_{2}^{\dagger }$ is given by

$$\begin{eqnarray}&&C={\rm max} (0,{{\lambda }_{1}}-{{\lambda }_{3}}-2\sqrt{{{\lambda }_{2}}{{\lambda }_{4}}}).\end{eqnarray} \tag{ 17 }$$

Applied to the case at hand, with ${{\rho }_{{\rm i}}}={{\tau }_{\beta }}\otimes {{\tau }_{\beta }}$ we finally obtain

$$\begin{eqnarray}&&{{C}_{{\rm max} }}={\rm max} (0,2{{p}^{2}}-p-2(1-p)\sqrt{p(1-p)}).\end{eqnarray} \tag{ 18 }$$

It follows therefore, that unlike when considering correlations, there is now a threshold temperature, ${{k}_{{\rm B}}}{{T}_{{\rm max} }}/E\approx 1.19$ (or equivalently a threshold ground-state population ${{p}_{{\rm min} }}\approx 0.698$), such that for all $T\geqslant {{T}_{{\rm max} }}$ (or $p\leqslant {{p}_{{\rm min} }}$) no entanglement can be created between the two qubits.

3.2.2. Multipartite systems

We now switch our attention to the multipartite setting. Here we will see that the limiting temperature ${{T}_{{\rm max} }}$ below which one can create entanglement can be increased when several thermal qubits are jointly processed. Essentially, as more qubits are available, the global system contains larger energy gaps and thus subspaces with higher purity, which can then potentially be more easily entangled. To illustrate this point we can look at the example of two thermal qubits.The renormalized subspace spanned by $\{|00\rangle ,|11\rangle \}$ constitutes a qubit with higher purity that the single thermal qubit and in multipartite systems we can make use of this by finding virtual two qubit subspaces with increased purity. In the following we make this intuition precise by studying the dependence of ${{T}_{{\rm max} }}$ on the number of qubits n. At the same time, we study several classes of entanglement that naturally appear in the multipartite case including its strongest form: GME.

Entanglement in all bipartitions.

To start our discussion, we consider the case of n qubits and a straightforward generalization of the above two-qubit protocol. That is, we consider a rotation in the $|0{{\rangle }^{\otimes n}}$, $|1{{\rangle }^{\otimes n}}$ subspace, of the form (16),

$$\begin{eqnarray}&&U=|{{\phi }^{n}}\rangle \langle 0{{|}^{\otimes n}}+|{{\phi }^{{{n}^{\prime }}}}\rangle \langle 1{{|}^{\otimes n}}+{\mathbb{1}}-(|0\rangle {{\langle 0|)}^{\otimes n}}-(|1\rangle {{\langle 1|)}^{\otimes n}},\end{eqnarray} \tag{ 19 }$$

where $|{{\phi }^{{{n}^{\prime }}}}\rangle =|\phi _{10\cdots 0}^{n}\rangle$. For a given bipartition $j|n-j$ (i.e. a partition of j qubits versus $n-j$ qubits), the concurrence in the final state ${{\rho }_{{\rm f}}}$ can be conveniently lower bounded using the relation [32]

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} C\geqslant 2\left( \left| {{\langle 0{{|}^{\otimes n}}{{\rho }_{{\rm f}}}|1\rangle }^{\otimes n}} \right| \right. \\ \quad -\left. \sqrt{{{\langle 0{{|}^{\otimes j}}{{\langle 1{{|}^{\otimes (n-j)}}{{\rho }_{{\rm f}}}|0\rangle }^{\otimes j}}|1\rangle }^{\otimes (n-j)}}}\sqrt{{{\langle 1{{|}^{\otimes j}}{{\langle 0{{|}^{\otimes (n-j)}}{{\rho }_{{\rm f}}}|1\rangle }^{\otimes j}}|0\rangle }^{\otimes (n-j)}}} \right) \\ \end{array}\end{eqnarray} \tag{ 20 }$$

and due to the simple form of ${{\rho }_{{\rm f}}}$, these bounds are in fact tight [33]. Evaluating explicitly, we then obtain

$$\begin{eqnarray}&&C={{\lambda }_{0}}-{{\lambda }_{n}}-2\sqrt{{{\lambda }_{j}}{{\lambda }_{n-j}}},{{\lambda }_{j}}={{\langle 0{{|}^{\otimes (n-j)}}{{\langle 1{{|}^{\otimes j}}{{\rho }_{{\rm i}}}|0\rangle }^{\otimes (n-j)}}|1\rangle }^{\otimes j}}.\end{eqnarray} \tag{ 21 }$$

which is independent of the bipartition, and given by

$$\begin{eqnarray}&&C={{p}^{n}}-{{(1-p)}^{n}}-2{{p}^{n/2}}{{(1-p)}^{n/2}}.\end{eqnarray} \tag{ 22 }$$

By demanding $C\gt 0$, we can characterize the smallest p, and thus the largest T, that allows for entanglement to be created simultaneously across all bipartitions, as a function of n. We find a linear scaling in n for this critical temperature $T_{{\rm E}}^{({\rm all}\,{\rm bip}.)}$,

$$\begin{eqnarray}&&\frac{{{k}_{{\rm B}}}T_{{\rm E}}^{({\rm all}\,{\rm bip}.)}}{E}\geqslant \frac{n}{2{\rm ln} (1+\sqrt{2})}.\end{eqnarray} \tag{ 23 }$$

Hence it follows that for every finite temperature T, there exists a sufficiently large number of qubits n that still allows for entanglement creation across all bipartitions. We note also that if one used instead of concurrence the negativity across a bipartition, a straightforward calculation shows that the same bound is obtained.

Entanglement in a single bipartition.

The above protocol can be improved if the aim is to generate entanglement in a given single bipartition $j|n-j$. As in the two-qubit protocol, the idea is to perform a permutation of the initial diagonal elements before applying the rotation (19). From expression (21), we see that the optimal permutation is the one where ${{\lambda }_{0}}={{p}^{n}}$, ${{\lambda }_{n}}={{\lambda }_{j}}=p{{(1-p)}^{n-1}}$ and ${{\lambda }_{n-j}}={{(1-p)}^{n}}$. In such a case, we a similar analysis to above leads to the limiting temperature, which, for large n is given by

$$\begin{eqnarray}&&\frac{{{k}_{{\rm B}}}{{T}_{{\rm E}}}}{E}\gtrsim \frac{n-1/2}{{\rm ln} (3)}.\end{eqnarray} \tag{ 24 }$$

Hence the threshold temperature for the creation of bipartite entanglement using this protocol is also linear in n (for high temperatures), but improves upon the above protocol in the constants. Thus for fixed n, one can generate entanglement across a single bipartition for slightly higher temperatures.

Genuine multipartite entanglement

GME is the strongest form of entanglement in multipartite systems. While entanglement across all partitions is a necessary requirement for GME it is by no means sufficient. A state ρ is GME iff it only admits decompositions of the form

$$\begin{eqnarray}&&\rho =\mathop{\sum }\limits_{i}{{p}_{i}}|{{\phi }_{i}}\rangle \langle {{\phi }_{i}}|\end{eqnarray} \tag{ 25 }$$

where at least one $|{{\phi }_{i}}\rangle$ is entangled in every possible bipartition. It follows that a necessary but not sufficient condition for GME is that ρ itself is entangled across every bipartition. This suggests that the previously considered protocol for generating entanglement in all bipartitions is a natural candidate to gain a first insight on the maximal temperature for GME creation.

After applying the unitary (19), the state ${{\rho }_{{\rm f}}}$ is essentially a GHZ-state mixed with (diagonal) noise. For such a simple form, the techniques of [32, 35] give us necessary and sufficient conditions for the creation of GME [34], namely

$$\begin{eqnarray}&&{{\rho }_{{\rm f}}}\;{\rm is}\;{\rm GME}\ \Longleftrightarrow \ {{p}^{n}}-{{(1-p)}^{n}}-2({{2}^{n-1}}-1){{p}^{n/2}}{{(1-p)}^{n/2}}\geqslant 0.\end{eqnarray} \tag{ 26 }$$

This condition leads to a lower bound on the threshold temperature for creating GME, ${{T}_{{\rm GME}}}$, which turns out to be asymptotically independent of n, and given by

$$\begin{eqnarray}&&\frac{{{k}_{{\rm B}}}T_{{\rm GME}}^{({\rm GHZ})}}{E}\simeq \frac{1}{2{\rm ln} (2)},\end{eqnarray} \tag{ 27 }$$

where we added the suffix GHZ because the target entangled state of this protocol is a GHZ state. Moreover, as we show in the appendix, this result holds for all states whose density matrix features only diagonal and anti-diagonal elements, also known as X-states [34, 36].

Geniune multipartite entanglement II

Recall that there are many inequivalent types of multipartite entangled states and GHZ states only constitute one prominent class. In fact it is much more favorable to use protocols that target another type of entangled states, namely Dicke states [37]. An n-qubit Dicke state with k excitations is defined as:

$$\begin{eqnarray}&&|D_{k}^{n}\rangle =\frac{1}{\sqrt{\left( \frac{n}{k} \right)}}\mathop{\sum }\limits_{j}{{P}_{j}}{{\{|1\rangle }^{\otimes k}}|0{{\rangle }^{\otimes (n-k)}}\},\end{eqnarray} \tag{ 28 }$$

where ${{\sum }_{j}}{{P}_{j}}\{\}$ is a sum over all possible permutations. Besides being relevant for the theory of light–matter interaction, Dicke states are useful for various quantum information tasks [38], have been detected experimentally [39, 40] and have shown to be GME [41, 42].

By constructing a protocol that uses the state (28) as the target entangled state, we will show that the threshold temperature for generating GME is given by

$$\begin{eqnarray}&&\frac{{{k}_{{\rm B}}}{{T}_{{\rm GME}}}}{E}\geqslant \frac{n}{(k+1){\rm ln} n}+\mathcal{O}\left[ \frac{n}{{{({\rm ln} n)}^{2}}} \right].\end{eqnarray} \tag{ 29 }$$

The scaling is almost linear with n, which allows now for the creation of GME for an arbitrarily high temperature $T\lt \infty$, by considering a sufficient number of qubits n. Moreover, the simple witnesses used to demonstrate GME are known not to be optimal in many cases [45] indicating the improvements might even be possible. Note that this result is quite counter-intuitive, as the complexity of the task we consider, entangling all qubits, increases with n. Furthermore, it is in stark contrast with the results obtained above for the GHZ class, and thus indicates that different types of entanglement behave in a very different manner.

Let us now sketch the idea of the protocol for creating Dicke type entanglement; all details are in appendix A.2. As in the previous cases, the protocol consists of two steps: a permutation of the diagonal elements followed by a rotation to maximally entangled states (in this case to Dicke states). The permutation first moves the largest eigenvalue, pⁿ, plus the small eigenvalues, ${{p}^{k}}{{(1-p)}^{n-k}}$, into the degenerate subspace of energy kE, thus purifying the subspace. It also moves other small eigenvalues⁵ into the subspaces of $k-1$ and $k+1$ excitations, as this is favorable for the considered entanglement witness [43]. Now, in the degenerate subspace of k excitations, the state with the biggest population pⁿ is rotated to the Dicke state (28). In order for the transformation to be unitary, the rest of the energy eigenvectors of the subspace are rotated to the set of orthonormal states

$$\begin{eqnarray}&&|d_{k,l}^{n}\rangle =\frac{1}{{{N}_{k}}}\mathop{\sum }\limits_{j}{{{\rm e}}^{{\rm i}\frac{2\pi lj}{{{N}_{k}}}}}{{P}_{j}}{{\{|1\rangle }^{\otimes k}}|0{{\rangle }^{\otimes n-k}}\},{{N}_{k}}=\sqrt{\left( \frac{n}{k} \right)},\end{eqnarray} \tag{ 30 }$$

with $j=\{1,\ldots ,{{N}_{k}}-1\}$. This concludes the protocol leading to (29) (see the appendix for detailed computations).

The fact that the creation of Dicke type GME is so much more favorable can be understood intuitively by recalling that Dicke states are in general much more robust to noise compared to GHZ states [43, 44]. Notice also from (29) that it is most favorable to create entanglement in the first excited subspace, where the Dicke state becomes the well-known W state.

Upper bounds and discussion

So far, we have investigated explicit protocols, which allowed us to place lower bounds on the threshold temperature that still allows for the creation of entanglement. To study the limitations imposed by a thermal background it is essential to also find upper bounds on the maximal temperature. For that purpose, a first approach is to use results on the geometry of quantum states. In particular, it is known that the maximally mixed state is always surrounded by a ball of finite size that contains only separable states, and it is possible to place lower bounds on the radius of such a ball [46, 47]. By applying these results we obtain an upper bound that scales exponentially with n. Therefore, there is an exponential gap between lower and upper bounds, thus making this approach essentially useless for large n.

The results from [46, 47] are useful for any state, as long as it is sufficiently close to the identity, whereas here we are concerned with a very particular form of states, namely those states with a thermal spectrum. This information can be used to obtain better upper bounds. Indeed, the following theorem was proven in [48]: let $\rho \in {{\mathcal{H}}_{2}}\otimes {{\mathcal{H}}_{d}}$ have eigenvalues ${{\lambda }_{1}}\geqslant {{\lambda }_{2}}\geqslant ...\geqslant {{\lambda }_{2d}}$, then

$$\begin{eqnarray}&&U\rho {{U}^{\dagger }}\;{\rm is}\;{\rm separable}\;\forall {\rm U}\;\ \Longleftrightarrow \ {{\lambda }_{1}}-{{\lambda }_{2d-1}}-2\sqrt{{{\lambda }_{2d-2}}{{\lambda }_{2d}}}\leqslant 0.\end{eqnarray} \tag{ 31 }$$

By taking $d={{2}^{n-1}}$, this criterion applies to any qubit versus qudit bipartition of the n-qubit thermal system we considered. Furthermore, notice that this condition amounts to calculating the concurrence in a specific $4\times 4$ subspace, which happens to be exactly the purest one we used in the protocol leading to (24). Hence that protocol is optimal for generating entanglement in any qubit versus qudit bipartition. While the possibility to obtain a better ${{T}_{{\rm E}}}$ in a qudit versus qudit bipartition remains open, this criterion does yield upper bounds for $T_{{\rm E}}^{({\rm all}\,{\rm bip}.)}$ and ${{T}_{{\rm GME}}}$ ⁶ , obtaining

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \frac{{{k}_{{\rm B}}}T_{{\rm E}}^{({\rm all}\,{\rm bip}.)}}{E}\leqslant \frac{n-1}{{\rm ln} 3}, \\ \quad \frac{{{k}_{{\rm B}}}{{T}_{{\rm GME}}}}{E}\leqslant \frac{n-1}{{\rm ln} 3}. \\ \end{array}\end{eqnarray} \tag{ 32 }$$

Therefore we obtain upper bounds on (23) and (29) that also scale linearly with n, showing that this scaling between the maximal temperature and the number of qubits is a fundamental property, and that our protocols perform close to optimal for entanglement and GME generation at high temperatures. The results are summarized in figure 1.

Zoom In Zoom Out Reset image size

The problem of the attainable entanglement in the unitary orbit of mixed states has been considered in the context of nuclear magnetic resonance (see [23] and references therein). The best protocol in [23] obtains precisely the scaling (23), improving on protocols based on algorithmic cooling [50] and effective pure states [49]. Our result (24) provides a tighter bound on the minimal temperature required for entanglement generation, and the upper bound derived from [48] gives evidence that it is tight⁷ . We also studied the minimal temperature for GME, finding a surprising positive scaling with the number of qubits. Our results thus provide bounds on the number of required qubits to generate entanglement and GME at finite temperature, while showing that in the asymptotic limit generation of entanglement and GME is possible at any temperature.

In analogy with the previous section, let us start by considering the case of a two qudit system, i.e. two d-level systems. The goal is now to maximize $I({{A}_{1}}:{{A}_{2}})$, as defined in (4), over all global unitaries constrained by $W\leqslant \Delta E$.

Note first that initially $I({{A}_{1}}:{{A}_{2}})=0$, as the initial state factorizes. Now, to create correlations, we must apply a global unitary that will increase the local entropies $S({{A}_{i}})$ of ${{\rho }_{{\rm f}}}$, since the total entropy $S({{A}_{1}}{{A}_{2}})=2S({{\tau }_{\beta }})$ will clearly not change. Recalling that the thermal state maximizes the entropy of a system with fixed average energy (see, for example, [27]), we find that

$$\begin{eqnarray}&&{{I}_{\Delta E}}\leqslant 2\left[ S({{\tau }_{{{\beta }^{\prime }}}})-S({{\tau }_{\beta }}) \right],\end{eqnarray} \tag{ 33 }$$

where $\beta ^{\prime}$ is chosen such that $\Delta E={\rm Tr}[{{H}_{{\rm tot}}}(\tau _{{{\beta }^{\prime }}}^{\otimes 2}-\tau _{\beta }^{\otimes 2})]$. Hence in order to obtain correlations at minimal energy cost, one should look for a protocol such that the local states of ${{\rho }_{f}}$ are thermal states at equal temperature. That is, the optimal unitary ${{U}^{*}}$ satisfies

$$\begin{eqnarray}&&{\rm T}{{{\rm r}}_{{{A}_{1}}}}({{U}^{*}}{{\rho }_{i}}{{U}^{*\dagger }})={\rm T}{{{\rm r}}_{{{A}_{2}}}}({{U}^{*}}{{\rho }_{i}}{{U}^{*\dagger }})={{\tau }_{{{\beta }^{\prime }}}}.\end{eqnarray} \tag{ 34 }$$

This unitary effectively heats up the system locally, while the global system preserves its entropy. In the appendix (first section) we construct ${{U}^{*}}$, thus reaching the bound (33), for Hamiltonians with equally spaced energy levels and for arbitrary Hamiltonians if the temperature difference is big enough. In figure 2 we illustrate our results for two qubits and for various values of ${{k}_{{\rm B}}}T/E$. Finally, notice that expression (33) recovers the case of maximal correlations, (33), in the limit $\beta ^{\prime} \to 0$, with a corresponding work cost

$$\begin{eqnarray}&&W=2\left( \frac{1}{d}{\rm Tr}H-\frac{1}{\mathcal{Z}}{\rm Tr}H{{e}^{-\beta H}} \right).\end{eqnarray} \tag{ 35 }$$

Zoom In Zoom Out Reset image size

These results are easily extendible to the multipartite case. The generalized mutual information (9) is maximized (for a given energy cost) by those unitaries that satisfy (34) for every local state.

4.2.1. Bipartite systems

Next we derive the minimal cost of creating entanglement for the simplest case of two qubits. Consider first the case T = 0, i.e. $\tau =|0\rangle \langle 0|$. If the state is pure, entanglement can be measured by the entropy of entanglement, which is simply given by the local entropy of the state. The problem at hand is thus equivalent to the maximization of the mutual information, so the same reasoning can be used here⁸ . In particular, the optimal unitary, ${{U}^{*}}$ in (34), can be generated by a rotation in the $|00\rangle$, $|11\rangle$ subspace. From this we find the relation

$$\begin{eqnarray}&&C=\sqrt{\frac{\Delta E}{E}\left( 2-\frac{\Delta E}{E} \right)}.\end{eqnarray} \tag{ 36 }$$

Moving to non-zero temperature, finding the optimal unitary is no longer straightforward. Nevertheless the problem can be attacked from two directions. First, we maximize C numerically, with respect to all possible unitaries, for a given cost W. Second, we use an ansatz protocol, inspired by the optimal unitaries to achieve C_max in (17). These unitaries have the form of first rotating in the subspace of $|10\rangle$ and $|11\rangle$, followed by rotating in the subspace of $|00\rangle$ and $|11\rangle$. Our ansatz is to optimise over such unitaries, now a much simpler optimization over the two unknown angles (one for each rotation). The results are presented in figure 3, where the solid line shows the result of the full optimization and the dashed line shows the results of the ansatz. We see that when there is no restriction on the amount of available energy W, then our ansatz protocol performs optimally. However, this is not the case when W is limited. Note that the amount of energy required to reach ${{C}_{{\rm max} }}$ is decreasing as T increases, shown in inset (a), where we also see that for low temperatures (${{k}_{{\rm B}}}T/E\lesssim 0.1$), we can generate essentially one Bell state of two qubits, i.e. ${{C}_{{\rm max} }}\simeq 1$. Moreover, for any $T\gt 0$, there is a minimal amount of energy required for generating entanglement, shown in inset (b). This is because some energy is always needed to leave the set of separable states,when starting from thermal product states.

Zoom In Zoom Out Reset image size

4.2.2. Multipartite systems

Quantification and characterization of multipartite entanglement is still a highly active field of research (see e.g. [26]). The main challenge is a consistent quantification of multipartite entanglement in operational terms. It seems that this task may not be as easy as in bipartite systems where in the LOCC paradigm entanglement can be quantified by a unique resource. Here we circumvent this issue by studying a measure independent question: what is the energy cost of transforming a thermal state into an entangled one, either GME or entangled in all bipartitions.

The work cost associated to the unitary (19) is easily computed to be

$$\begin{eqnarray}&&W=\frac{nE(1-{{{\rm e}}^{-\beta En}})}{2{{(1+{{{\rm e}}^{-\beta E}})}^{n}}}.\end{eqnarray} \tag{ 37 }$$

By inserting ${{T}^{{\rm all}\,{\rm bip}.}}$ in (37), one obtains that the cost to leave the separable set (for this particular protocol) is exactly

$$\begin{eqnarray}&&{{W}^{{\rm sep}.}}=nE\frac{1+\sqrt{2}}{{{\left( {{\left( 1+\sqrt{2} \right)}^{2/n}}+1 \right)}^{n}}}\end{eqnarray} \tag{ 38 }$$

which is exponentially small in n. This shows that having more qubits not only opens the possibility to generate entanglement at a higher temperature, but also reduces the energy cost of leaving the separable set. An exponential decrease of the work cost with n is also found for the other protocols for GME generation in the multipartite setting (see appendices). The reason behind this behavior is that the considered protocols only act on particular subspaces, whose population becomes negligible in the limit of large n. This also implies that the amount of generated entanglement decreases with the number of qubits. Interestingly, in the multipartite setting, even a small amount of entanglement might be enough to obtain a substantial quantum advantage. In particular, in the field of quantum computation, for a computational speed up (in pure states) entanglement is required across every bipartition [51, 52], but the actual amount can be polynomially small in the system's size [53] ⁹ . Therefore, our protocols for multipartite entanglement generation are not only interesting from a fundamental point of view—as they set fundamental bounds on the maximal temperatures—but might find applications in the field of quantum computation.

Finding protocols that generate a substantial amount of entanglement and GME at high temperatures remains as an important future direction, as this would give a bigger resistance to noise and is important for other applications of GME, such as metrology [54].

In this appendix we prove that for equal energy spacings (${{E}_{k}}=k{{E}_{0}}$) any local temperature $T^{\prime} \gt T$ lies in the unitary orbit of thermal product states at temperature T. Furthermore this proof provides a constructive protocol that in this context reaches any amount of mutual information at minimal possible energy cost. Furthermore for general Hamiltonians this protocols provides a means to reach any sufficiently larger T' (for an exact condition see below) at minimal energy cost. First let us adopt the following convenient notation for the eigenvalue distribution of the marginal ${{\vec{p}}_{A}}:={\rm diag}({{\rho }_{A}})$, which is sufficient for our purposes as the thermal states will always just be diagonal in the energy eigenbasis. The general idea of the protocol that follows is that the global unitary should induce a doubly stochastic transformation M on the marginal probability vector, while ensuring that no coherences are created in any subsystem. First we decompose the marginal vectors as

$$\begin{eqnarray}&&{\rm diag}({{\rho }_{A}})=\mathop{\sum }\limits_{i=0}^{d-1}{{\vec{p}}_{i}}\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&{\rm diag}({{\rho }_{B}})=\mathop{\sum }\limits_{i=0}^{d-1}{{\Pi }^{i}}{{\vec{p}}_{i}}\end{eqnarray} \tag{ A.2 }$$

with $p_{i}^{j}=\langle j|\otimes \langle j+i|\rho |j\rangle \otimes |j+i\rangle ={{{\rm e}}^{-\beta ({{E}_{j}}+{{E}_{j+1}})}}/{{\mathcal{Z}}^{2}}$ and $\Pi ={{\sum }_{k}}|k\rangle \langle k+1|$. It will be useful to consider again the following generalized Bell states

$$\begin{eqnarray}&&|{{B}_{i,j}}\rangle =\mathop{\sum }\limits_{k=0}^{d-1}{{\omega }^{ki}}|k\rangle \otimes |k+j\rangle ,\end{eqnarray} \tag{ A.3 }$$

with $\omega ={{{\rm e}}^{2\pi {\rm i}/d}}$. Now it is straightforward to see that rotating in the subspaces spanned by ${{S}_{i}}={\rm span}(\{{{B}_{0,i}},{{B}_{1,i}},\cdots ,{{B}_{d-1,i}}\})$ ensures that every diagonal element that can be created by these rotations is being traced over. Unitarity ensures that rotations in the subspaces S_i induce a doubly stochastic transformation of the diagonal part of the density matrix in this subspace. Now all that is left is to observe that the probabilities in the decomposition of the subsystems correspond exactly to the rotations in the subspaces spanned by the maximally entangled states defined before, i.e.

$$\begin{eqnarray}&&{\rm diag}({\rm T}{{{\rm r}}_{B}}(U\rho {{U}^{\dagger }}))=\mathop{\sum }\limits_{i}{{M}_{i}}{{\vec{p}}_{i}}\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&{\rm diag}({\rm T}{{{\rm r}}_{A}}(U\rho {{U}^{\dagger }}))=\mathop{\sum }\limits_{i}{{\Pi }^{i}}({{M}_{i}}{{\vec{p}}_{i}}),\end{eqnarray} \tag{ A.5 }$$

where each M_i is a doubly stochastic matrix. Now we can use the symmetry of the initial state, i.e. ${{p}_{ij}}={{p}_{ji}}$, and define a target doubly stochastic matrix M that should describe the transformation of both marginals. Since the vector ${{\vec{p}}_{0}}$ is equal for both marginals it is evident that ${{M}_{0}}=M$ already takes the first part out of the picture without restricting the generality of stochastic transformations. In general if every doubly stochastic matrix is equal, i.e. ${{M}_{i}}=M$ and commutes with all ${{\Pi }^{i}}$, i.e. is a circulant matrix, it is evident (due to $\sum _{i=0}^{d-1}{{\vec{p}}_{i}}=\sum _{i=0}^{d-1}{{\Pi }^{i}}{{\vec{p}}_{i}}$) that both subsystems' probability vector will just be transformed by M. In other words it is easily achievable to transform the subsystems probability distribution by any doubly stochastic matrix that commutes with all cyclic permutations, i.e. a circulant matrix.

The final question is thus whether circulant doubly stochastic transformations of the form $T={{\sum }_{i}}{{\alpha }_{i}}{{\Pi }^{i}}$ are sufficient to reach any temperature, i.e. $T\vec{p}(\beta )=\vec{p}({{\beta }^{\prime }})\forall \beta \geqslant \beta ^{\prime}$? Obviously starting from $\beta =\infty$ one can reach all temperatures via choosing ${{\alpha }_{i}}={{p}_{i}}({{\beta }^{\prime }})$ and from any β one can reach the infinite temperature distribution via ${{\alpha }_{i}}=\frac{1}{d}$ for all i.

To address this question to its fullest extent we will first construct a general convex sum of cyclic permutations that achieves this general task and then check whether all coefficients are positive. We require that

$$\begin{eqnarray}&&\mathop{\sum }\limits_{i}{{\alpha }_{i}}{{\Pi }^{x}}{{p}_{i}}(\beta )={{p}_{x}}({{\beta }^{\prime }}).\end{eqnarray} \tag{ A.6 }$$

For sake of simplicity we will first define $\alpha _{i}^{\prime }={{\alpha }_{i}}\frac{Z}{Z^{\prime} }$ such that the condition is simplified to

$$\begin{eqnarray}&&\mathop{\sum }\limits_{i}\alpha _{i}^{\prime }{{{\rm e}}^{-\beta {{E}_{i+x}}}}={{{\rm e}}^{-{{\beta }^{\prime }}{{E}_{x}}}}.\end{eqnarray} \tag{ A.7 }$$

A set of $\alpha _{k}^{\prime }$ solving this equation system is given by

$$\begin{eqnarray}&&\alpha _{k}^{\prime }=\frac{{{\epsilon }_{{{i}_{1}}{{i}_{2}}(...){{i}_{n}}}}{{{\rm e}}^{-\beta {{E}_{{{i}_{1}}}}}}{{{\rm e}}^{-\beta {{E}_{{{i}_{2}}-1}}}}(\cdots ){{{\rm e}}^{-{{\beta }^{\prime }}{{E}_{{{i}_{n}}-n-1-k}}}}}{{{\epsilon }_{{{i}_{1}}{{i}_{2}}(...){{i}_{n}}}}{{{\rm e}}^{-\beta {{E}_{{{i}_{1}}}}}}{{{\rm e}}^{-\beta {{E}_{{{i}_{2}}-1}}}}(\cdots ){{{\rm e}}^{-\beta {{E}_{{{i}_{n}}-n-1}}}}}.\end{eqnarray} \tag{ A.8 }$$

From this explicit form we can easily find negative coefficients and thus prove that circulant matrices are insufficient to reach any arbitrarily higher temperature. On the other hand we immediately see from very simple geometric considerations that for a sufficiently high difference in temperatures $\Delta T=T^{\prime} -T$ circulant matrices are always sufficient. Since the original probability vector is linearly independent from all its cyclic permutations and all of them are equally far in Euclidean distance from the infinite temperature distribution we can study the convex cone with ${{\Pi }^{i}}{{\vec{p}}_{A}}$ as extremal rays. Since all of the extremal rays share the same distance to the center ray (infinite temperature), we know that a sufficient condition for circulant matrices to achieve the higher temperature Boltzmann distribution is simply given by the minimal distance from the central ray to all faces of the cone. This is always easily calculable for any energy distribution and gives a sufficient condition on $\Delta T$ for this protocol to work.

Furthermore we can use the explicit solution to find Hamiltonians for which this protocol always works. One important example is equal energy spacing between the different levels, i.e. ${{E}_{k}}=k{{E}_{0}}$. In this case the explicit solution for the ${{\alpha }_{k}}$ is given as

$$\begin{eqnarray}&&{{\alpha }_{k}}=\frac{1-Z^{\prime} p_{1}^{\prime }}{1-Z{{p}_{1}}}\left( {{\delta }_{0,k}}+\frac{Z^{\prime} p_{1}^{\prime }-Z{{p}_{1}}}{1-Z^{\prime} p_{1}^{\prime }}p_{k+d{{\delta }_{0,k}}-1}^{\prime } \right)\end{eqnarray} \tag{ A.9 }$$

which is positive for all k due to the fact that ${{p}_{0}}\gt p_{0}^{\prime }$, i.e. we have derived a protocol that delivers the maximally possible amount of mutual information at minimum energy costs for all Hamiltonians with equal energy spacing (and thus qubits as a special case).

Given an n-qubit thermal state $\Omega =\tau _{\beta }^{\otimes n}$ we here find the $n\gg 1$ asymptotic behaviour of the maximal temperature T_GME that allows to unitarily create GME in the ensemble with the W-state protocol and also calculate the energy cost of the protocol. Here ${{\tau }_{\beta }}={\rm diag}(p,vp)$, where $p=1/(1+v)$ and $v={{{\rm e}}^{-\beta E}}$ is the Boltzmann weight. If the eigenvectors corresponding to the first excited level of the total Hamiltonian are $\{|w_{i}^{(1)}\rangle \}_{i=1}^{n}$ and the ones corresponding to the second excited level are $\{|w_{a}^{(2)}\rangle \}_{a=1}^{n(n-1)/2}$, then the measure we use has the form [55]

$$\begin{eqnarray}&&\mathcal{E}[\Omega ]\;=\;\mathop{\sum }\limits_{i\ne j}\;|{{\Omega }_{ij}}|\;-2\;\sqrt{{{\Omega }_{00}}}\mathop{\sum }\limits_{a}\;\sqrt{{{\Omega }_{aa}}}-\;(n\;-2)\;\mathop{\sum }\limits_{i}\;{{\Omega }_{ii}},\end{eqnarray} \tag{ A.10 }$$

where ${{\Omega }_{ij}}=\langle w_{i}^{(1)}|\Omega |w_{j}^{(1)}\rangle$ and ${{\Omega }_{ab}}=\langle w_{a}^{(2)}|\Omega |w_{b}^{(2)}\rangle$.

In short, the W-state protocol is the maximization of $\mathcal{E}$ over all such unitary operations that generate non-diagonal elements only in the eigensubspace of the first excited level (which we denote by ${{\mathcal{W}}_{1}}$). These unitaries can be represented as $U\Pi$, where Π is a permutation operation on the initial state and U is a general $n\times n$ unitary living in ${{\mathcal{W}}_{1}}$. As this representation suggests, we divide the optimization procedure in two steps: (i) maximization over Us for a given Π, and (ii) maximization over Πs. After Π acts, the state becomes ${{\Omega }^{\Pi }}=\Pi \rho \Pi$ and its projection on ${{\mathcal{W}}_{1}}$ we denote by ${{\omega }^{\Pi }}$. Now, the operation U will act only on ${{\omega }^{\Pi }}$ and take it to ${{\omega }^{\prime }}=U{{\omega }^{\Pi }}{{U}^{\dagger }}$ and since U is unitary, the traces of $\omega ^{\prime}$ and ${{\omega }^{\Pi }}$ will be the same. Therefore, we can rewrite (A.10) as

$$\begin{eqnarray}&&\mathcal{E}[U\Pi \Omega \Pi {{U}^{\dagger }}]\;=\;\mathop{\sum }\limits_{i\ne j}|\omega _{ij}^{\prime }|\;-\;2\sqrt{\Omega _{00}^{\Pi }}\mathop{\sum }\limits_{\alpha }\sqrt{\Omega _{\alpha \alpha }^{\Pi }}-(n-2)\mathop{\sum }\limits_{i}\omega _{ii}^{\Pi }.\end{eqnarray} \tag{ A.11 }$$

This shows that the maximization of $\mathcal{E}$ over U is reduced to the maximization of ${{\sum }_{i\ne j}}|\omega _{ij}^{\prime }|$ over U. To find this maximum, we first observe that due to the unitarity of U, ${\rm Tr}\left( {{({{\omega }^{\prime }})}^{2}} \right)={\rm Tr}\left( {{({{\omega }^{\Pi }})}^{2}} \right)$; whence

$$\begin{eqnarray}&&\mathop{\sum }\limits_{i\lt j}|\omega _{ij}^{\prime }{{|}^{2}}=\frac{{\rm Tr}\left( {{({{\omega }^{\Pi }})}^{2}} \right)-{{\sum }_{i}}{{(\omega _{i}^{\prime })}^{2}}}{2}.\end{eqnarray} \tag{ A.12 }$$

We now relax for the moment the constraint that $\omega ^{\prime}$ and ${{\omega }^{\Pi }}$ are unitarily connected and only require that ${\rm Tr}(\omega ^{\prime} )={\rm Tr}\left( {{\omega }^{\Pi }} \right)\equiv \alpha$ and ${\rm Tr}\left( {{({{\omega }^{\prime }})}^{2}} \right)={\rm Tr}\left( {{({{\omega }^{\Pi }})}^{2}} \right)\equiv {{\alpha }^{2}}\lambda$. Here we again divide the optimization in two steps: (1) maximize ${{\sum }_{i\lt j}}|\omega _{ij}^{\prime }|$ with ${{\sum }_{i\lt j}}|\omega _{ij}^{\prime }{{|}^{2}}$ fixed and (2) maximize the latter. Now we notice that

• (1)  
  The maximum is reached for $|\omega _{ij}^{\prime }|=|\omega _{i^{\prime} j^{\prime} }^{\prime }|$ and therefore ${\rm max} {{\sum }_{i\lt j}}|\omega _{ij}^{\prime }|=\sqrt{\frac{n(n-1)}{2}{{\sum }_{i\lt j}}|\omega _{ij}^{\prime }{{|}^{2}}}$.
• (2)  
  From (A.12), the maximum for ${{\sum }_{i\lt j}}|\omega _{ij}^{\prime }{{|}^{2}}$ is reached when ${{\sum }_{i}}{{(\omega _{i}^{\prime })}^{2}}$ is minimal. Since ${{\sum }_{i}}\omega _{i}^{\prime }=\alpha$ is fixed, the minimum for ${{\sum }_{i}}{{(\omega _{i}^{\prime })}^{2}}$ is reached when all $\omega _{i}^{\prime }=\frac{\alpha }{n}$, i.e., ${\rm max} {{\sum }_{i\lt j}}|\omega _{ij}^{\prime }{{|}^{2}}={{\alpha }^{2}}(\lambda -1/n)/2$.

Finally

$$\begin{eqnarray}&&{\rm max} {{\sum }_{i\ne j}}|\omega _{ij}^{\prime }|=\alpha \sqrt{n(n-1)(\lambda -1/n)},\end{eqnarray} \tag{ A.13 }$$

and on this maximum, $\omega ^{\prime}$ has the following form:

$$\begin{eqnarray}\alpha \;\left( \;\begin{array}{ccccccccccccccc} \frac{1}{n} & \;{{{\rm e}}^{{{\phi }_{12}}}}\;\sqrt{\;\frac{\lambda -1/n}{n(n-1)}} & \;\cdots \; & {{{\rm e}}^{{{\phi }_{1n}}}}\;\sqrt{\;\frac{\lambda -1/n}{n(n-1)}} \\ {{{\rm e}}^{{{\phi }_{21}}}}\;\sqrt{\;\frac{\lambda -1/n}{n(n-1)}} & \;\frac{1}{n} & \;\cdots \; & {{{\rm e}}^{{{\phi }_{2n}}}}\;\sqrt{\;\frac{\lambda -1/n}{n(n-1)}} \\ {} & \; & \;\vdots \; & {} \\ {{{\rm e}}^{{{\phi }_{{{N}_{1}}}}}}\;\sqrt{\;\frac{\lambda -1/n}{n(n-1)}} & \;{{{\rm e}}^{{{\phi }_{{{N}_{2}}}}}}\;\sqrt{\;\frac{\lambda -1/n}{n(n-1)}} & \;\cdots \; & \frac{1}{n} \\ \end{array}\; \right).\end{eqnarray} \tag{ A.14 }$$

Obviously, being obtained in less restrictive conditions, (A.13) upper-bounds the sought $\mathop{{\rm max} }\limits_{U}{{\sum }_{i\ne j}}|\omega _{ij}^{\prime }|$. Nevertheless, one can prove, that for suitably chosen $\{{{\phi }_{ij}}\}$ the matrix in (A.14) can always be unitarily reached from ${{\omega }^{\Pi }}$. The proof is slightly more involved and is conducted by first proving the statement for n = 3 by explicitly calculating the corresponding phases (only one phase is necessary to adjust there) and then proving the statement by induction for any n.

Now, having done the maximization over U, we turn to finding the Π with largest

$$\begin{eqnarray}&&{{\mathcal{E}}^{\Pi }}=\alpha \left( \sqrt{n(n-1)(\lambda -1/n)}-n+2 \right)-2\sqrt{\Omega _{00}^{\Pi }}\mathop{\sum }\limits_{\alpha }\sqrt{\Omega _{\alpha \alpha }^{\Pi }},\end{eqnarray} \tag{ A.15 }$$

where we have plugged (A.13) in (A.11). The quantity λ is defined above as the sum of the squares of the normalized elements of ${{\omega }^{\Pi }}$. Therefore, it is never bigger than 1 which implies that in the $n\to \infty$ limit, ${{\mathcal{E}}^{\Pi }}$ in (A.15) will be non-negative only if $\lambda \to 1$. On the other hand, choosing a bigger α and smaller elements in the eigensubspace of the second excited level (which we denote by ${{\mathcal{W}}_{2}}$) and on the ground state will also make ${{\mathcal{E}}^{\Pi }}$ bigger. To fulfil all this we choose Π so that it takes the smallest element of Ω, ${{p}^{n}}{{v}^{n}}$, to the ground state, the biggest one, pⁿ, to ${{\mathcal{W}}_{1}}$. The rest of $(n-1)$ elements in ${{\mathcal{W}}_{1}}$ are chosen so that they are significantly smaller than pⁿ. We will take them to be all equal (so that they keep α as big as possible) and to be ${{p}^{n}}{{v}^{n-k}}$ with some k that will be discussed later on. Also, we will choose the elements in ${{\mathcal{W}}_{2}}$ to be ${{p}^{n}}{{v}^{n-m}}$ with some m that is small and independent of n. At this point we do not know which exact choice of k and m will maximize ${{\mathcal{E}}^{\Pi }}$, but fortunately the existing information about them is enough to deduce the asymptotic behavior we need.

With above described Π we have

$$\begin{eqnarray}&&\lambda =\frac{1+(n-1){{v}^{2(n-k)}}}{{{\left( 1+(n-1){{v}^{n-k}} \right)}^{2}}}.\end{eqnarray} \tag{ A.16 }$$

So, to have $\lambda \to 1$, $n{{v}^{n-k}}$ must $\to \;0$. With this condition and some algebraic manipulations employing Taylor expansions, we arrive at the following asymptotic expansion:

$$\begin{eqnarray}&&\frac{{{\mathcal{E}}^{\Pi }}}{{{p}^{n}}}=1-{{n}^{2}}{{v}^{n-k-m/2}}{{C}_{n}}+\mathcal{O}\left[ {{n}^{3}}{{v}^{2(n-k)}} \right],\end{eqnarray} \tag{ A.17 }$$

where ${{C}_{n}}={{(1-1/n)}^{2}}\left( {{v}^{m/2}}+{{v}^{k}}n/(n-1) \right)$ and is always $\mathcal{O}[1]$ since $v\lt 1$ and k and m are positive. With this, we rewrite (A.17) as

$$\begin{eqnarray}&&\frac{{{\mathcal{E}}^{\Pi }}}{{{p}^{n}}}=1-{{n}^{2}}{{v}^{n-k-m/2}}{{C}_{n}}\left( 1+\mathcal{O}\left[ n{{v}^{n-k}} \right] \right).\end{eqnarray} \tag{ A.18 }$$

Having in mind that $n{{v}^{n-k}}\to 0$ and explicitly indicating the dependence of k and T (and hence v) on n we obtain from (A.18) the asymptotic condition of the positivity of ${{\mathcal{E}}^{\Pi }}$ in the following form:

$$\begin{eqnarray}&&1\geqslant {{n}^{2}}v_{n}^{n-{{k}_{n}}-m/2}{{C}_{n}}={{{\rm e}}^{{\rm ln} {{C}_{n}}+2{\rm ln} n-\frac{n-{{k}_{n}}-m/2}{{{T}_{{\rm GME}}}/E}}}.\end{eqnarray} \tag{ A.19 }$$

From formula (A.19) it is now obvious that to maximize T_GME, k_n has to be as small as possible. So, whatever the k_n and m delivering the maximum are, they are finite numbers independent of n. Therefore,

$$\begin{eqnarray}&&T_{{\rm GME}}^{{\rm max} }=\frac{nE}{2{\rm ln} n}+\mathcal{O}\left[ \frac{nE}{{{({\rm ln} n)}^{2}}} \right].\end{eqnarray} \tag{ A.20 }$$

Finally the energy input required for such a scaling can simply be calculated from the prior permutations Π alone, as all subsequent rotations are performed in a degenerate subspace. Adding the cost of all the permutations above gives the rather cumbersome formula for the energy cost of the W-state protocol as

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} W=E{{(1-{{{\rm e}}^{-\beta E}})}^{-n}}\left( (n-1)({{{\rm e}}^{-\beta E}}-{{{\rm e}}^{-\beta E(n-1)}})+(1-{{{\rm e}}^{-\beta E}})+n{{{\rm e}}^{-\beta E(n-3)}} \right. \\ \qquad -\left. n{{{\rm e}}^{-\beta En}}+({{n}^{2}}-n)({{{\rm e}}^{-2\beta E}}-{{{\rm e}}^{-(n-2)\beta E}})+3({{{\rm e}}^{-\beta En}}-{{{\rm e}}^{-\beta E(n-3)}}) \right) \\ \end{array}\end{eqnarray} \tag{ A.21 }$$

which while seemingly complicated due to the numerous required permutations still remains exponentially small in n for any $T\gt 0$.

The W-state is but an element of a larger set of Dicke states. Correspondingly, our W-state protocol can be straightforwardly generalized to Dicke state protocols. First, let us introduce the m excitation Dicke states for n qubits:

$$\begin{eqnarray}&&|{{D}_{m}}\rangle =\frac{1}{\sqrt{(\begin{array}{ccccccccccccccc} n \\ m \\ \end{array})}}\mathop{\sum }\limits_{\{\alpha \}}|\{\alpha \}\rangle ,\end{eqnarray} \tag{ A.22 }$$

where $\{\alpha \}{\rm s}$ are the subsets of $\{i\}_{i=1}^{n}$ consisting of m elements, $|\{\alpha \}\rangle ={{\otimes }_{i\in \{\alpha \}}}|1{{\rangle }_{i}}{{\otimes }_{j\in \{j\slashed{\in }\{\alpha \}\}}}|0{{\rangle }_{j}}$, and the summation runs over all $(\begin{array}{ccccccccccccccc} n \\ m \\ \end{array})$ possible $\{\alpha \}$ s. Accordingly, the Dicke state protocol is the one when one is allowed to create non-diagonal elements only in ${{\mathcal{D}}_{m}}$—the subspace spanned by $|\{\alpha \}\rangle$ s. In that case, the GME witness is as follows [55]:

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} {{\mathcal{E}}_{m}}[\Omega ]=\mathop{\sum }\limits_{\{\gamma \}}\left( |\langle \{\alpha \}|\Omega |\{\beta \}\rangle |-\sqrt{\langle \{\alpha \}|\otimes \langle \{\beta \}|{{\Pi }_{\{\alpha \}}}\Omega \otimes \Omega {{\Pi }_{\{\alpha \}}}|\{\alpha \}\rangle \otimes |\{\beta \}\rangle } \right) \\ \qquad \qquad {\mkern 1mu} -n(n-m-1)\mathop{\sum }\limits_{\{\alpha \}}\langle \{\alpha \}|\Omega |\{\alpha \}\rangle , \\ \end{array}\end{eqnarray} \tag{ A.23 }$$

where the set $\{\gamma \}$ is the collection of all possible $(\{\alpha \},\{\beta \})$ with $\{\beta \}\in {{\mathcal{D}}_{m}}$ and such that the intersection $\{\alpha \}\cap \{\beta \}$ contains $m-1$ elements. As is straightforward to check, $\{\gamma \}$ has $m(n-m)(\begin{array}{ccccccccccccccc} n \\ m \\ \end{array})$ elements. ${{\Pi }_{\{\alpha \}}}$ is a permutation operator which, acted on some $|0...1...0...1\rangle \otimes |1...0...0...1\rangle$, swaps the parts of the vectors corresponding to $\{\alpha \}$ so that it takes first vector to ${{\mathcal{D}}_{m-1}}$ and the second one to ${{\mathcal{D}}_{m+1}}$; e.g., ${{\Pi }_{\{2,3\}}}|01100\rangle \otimes |11000\rangle \;=\;$ $|01000\rangle \otimes |11100\rangle$ (see [55] for more detailed explanations).

As above, the idea is to maximize ${{\mathcal{E}}_{m}}[U\Pi \Omega \Pi {{U}^{\dagger }}]$ over all unitaries U acting in ${{\mathcal{D}}_{m}}$ and permutations Π. Again, for a fixed Π one has to maximize ${{\sum }_{\{\gamma \}}}|\langle \{\alpha \}|\Omega |\{\beta \}\rangle |$, but since $\{\gamma \}$ does not run over all non-diagonal elements the form (A.14) may not necessarily be the optimal one. Nevertheless, since finding maximum of the sum of absolute values of the part of non-diagonal elements of a matrix appears to be a formidable task, we will use the form (A.14) as an ansatz. In what follows we will show that the asymptotic behavior for $T_{{\rm GME}}^{{\rm max} }$ following from this ansatz is very close to the optimal one. As in the previous case, the permutation delivering the optimal asymptotics will be the one that puts pⁿ and ${{p}^{n}}{{v}^{n-k}}$ (with $k\geqslant m$ finite but sufficiently big) in ${{\mathcal{D}}_{m}}$ and fills ${{\mathcal{D}}_{m-1}}$ and ${{\mathcal{D}}_{m+1}}$ with some ${{p}^{n}}{{v}^{n-l}}$ (with sufficiently big and finite l). With this, after simple manipulations we arrive at

$$\begin{eqnarray}&&\mathcal{E}_{m}^{\Pi }={{p}^{n}}m\left( 1+((\begin{array}{ccccccccccccccc} n \\ m \\ \end{array})-1){{v}^{n-k}}-(\begin{array}{ccccccccccccccc} n \\ m \\ \end{array})(n-m){{v}^{n-k}}-(\begin{array}{ccccccccccccccc} n \\ m \\ \end{array})(n-m){{v}^{n-l}} \right.+\left. \mathcal{O}\left( {{n}^{2m}}{{v}^{2n}} \right) \right).\end{eqnarray} \tag{ A.24 }$$

So the condition for the presence of GME, $\mathcal{E}_{m}^{\Pi }\geqslant 0$ reduces to¹⁰

$$\begin{eqnarray}&&{{{\rm e}}^{(m+1){\rm ln} n-\frac{nE}{T}}}\leqslant 1\end{eqnarray} \tag{ A.25 }$$

whence we obtain

$$\begin{eqnarray}&&T\sim \frac{nE}{(m+1){\rm ln} n}\end{eqnarray} \tag{ A.26 }$$

implying that

$$\begin{eqnarray}&&T_{{\rm GME}}^{{\rm max} }\geqslant \frac{nE}{(m+1){\rm ln} n}\end{eqnarray} \tag{ A.27 }$$

for large n.

Now, returning to the question of how close to the optimal this scaling is, let us observe that the maximum for ${{\sum }_{\{\gamma \}}}|\langle \{\alpha \}|\Omega |\{\beta \}\rangle |$ is given by $\sqrt{(\lambda -1/N)N(N-1)}$ with N in this case being $(\begin{array}{ccccccccccccccc} n \\ m \\ \end{array})$ (see (A.14) and the reasoning preceding it). This value for ${{\sum }_{\{\gamma \}}}|\langle \{\alpha \}|\Omega |\{\beta \}\rangle |$ is not necessarily unitarily achievable from the initial diagonal state but is clearly an upper bound for it. Proceeding as above with this ansatz we obtain

$$\begin{eqnarray}&&T_{{\rm GME}}^{{\rm max} }\leqslant \frac{nE}{m{\rm ln} n}\end{eqnarray} \tag{ A.28 }$$

showing that the initial ansatz (A.14) is quite reasonable and that in any case $T_{{\rm GME}}^{{\rm max} }=\mathcal{O}\left( \frac{n}{{\rm ln} n} \right)$ for all m.

Given a set of n thermal qubits, $\tau _{\beta }^{\otimes n}$, in this section we study the limitations for entanglement creation with unitary transformations of the form:

$$\begin{eqnarray}&&\tau _{\beta }^{\otimes n}\to \hat{X}\end{eqnarray} \tag{ A.29 }$$

where in the energy eigenbasis $\hat{X}$ takes the form

$$\begin{eqnarray}\hat{X}=\left( \begin{array}{ccccccccccccccc} {{a}_{1}} & {} & {} & {} & {} & {} & {} & {{z}_{1}} \\ {} & {{a}_{2}} & {} & {} & {} & {} & {{z}_{2}} & {} \\ {} & {} & ... & {} & {} & ... & {} & {} \\ {} & {} & {} & {{a}_{n}} & {{z}_{n}} & {} & {} & {} \\ {} & {} & {} & z_{n}^{*} & {{b}_{n}} & {} & {} & {} \\ {} & {} & ... & {} & {} & ... & {} & {} \\ {} & z_{2}^{*} & {} & {} & {} & {} & {{b}_{2}} & {} \\ z_{1}^{*} & {} & {} & {} & {} & {} & {} & {{b}_{1}} \\ \end{array} \right),\end{eqnarray} \tag{ A.30 }$$

with $n={{2}^{N-1}}$, $|{{z}_{i}}|\leqslant \sqrt{{{a}_{i}}{{b}_{i}}}$ and ${{\sum }_{i}}({{a}_{i}}+{{b}_{i}})=1$ to ensure that $\hat{X}$ is positive and normalized (see [32] for details). A relevant example of an X-like matrix is

$$\begin{eqnarray}&&\rho =|{\rm GHZ}\rangle \langle {\rm GHZ}|+\frac{{\mathbb{1}}}{{{2}^{n}}},\end{eqnarray} \tag{ A.31 }$$

where $|{\rm GHZ}\rangle =\frac{1}{\sqrt{2}}(|0...0\rangle +|1...1\rangle )$. As shown in [32], the GME n-qubit states of the form (A.30) can be computed by the genuine multipartite concurrence

$$\begin{eqnarray}&&{{C}_{{\rm GM}}}=2{\rm max} \{0,|{{z}_{i}}|-{{w}_{i}}\},\ \;i=0,1,\ldots ,n,\end{eqnarray} \tag{ A.32 }$$

where ${{w}_{i}}=\mathop{\sum }\limits_{j\ne i}^{n}\sqrt{{{a}_{j}}{{b}_{j}}}$.

We wish to maximize (A.32) over all U acting on (A.29). The initial state, $\tau _{\beta }^{\otimes n}$, has no off-diagonal terms in the energy eigenbasis. It is then advantageous to apply a unitary operation that only generates one off-diagonal term. Indeed, creating off-diagonal terms results into a stochastic transformation of the diagonal terms, thus increasing the w_i term in (A.32) while the $|{{z}_{i}}|$ term depends only on the largest off-diagonal term. On the other hand, given two diagonal elements a_i, b_i of $\tau _{\beta }^{\otimes n}$, the biggest off-diagonal term that can be generated by a unitary operation is $|{{a}_{i}}-{{b}_{i}}|/2$, which is obtained by a rotation to the corresponding Bell states. Therefore, the optimal protocol can be thought as a combination of:

• (i)  
  Rotate two diagonal elements to Bell states in order to maximize $|{{z}_{i}}|$ in (A.32).
• (ii)  
  Permute the rest of diagonal elements to minimize w_i in (A.32). This is implemented by setting the elements in decreasing order (in ${{w}_{i}}=\sum _{j\ne i}^{n}\sqrt{{{a}_{j}}{{b}_{j}}}$, having product of biggest with smallest, second biggest with second smallest, etc).

Step 1 is optimized by acting on the ground state and the most excited state. On the other hand, the thermal state is already ordered to optimize Step 2. The first step leads to $|{{z}_{i}}|=(1-{{{\rm e}}^{-n\beta \epsilon }})/{{\mathcal{Z}}^{n}}$; and since ${{a}_{i}}{{b}_{i}}={{{\rm e}}^{-\beta \epsilon n}}$, we obtain that ${{w}_{i}}=({{2}^{n}}-2){{{\rm e}}^{-\beta \epsilon n/2}}/{{\mathcal{Z}}^{n}}$. In the limit of large n, one easily obtains that ${{k}_{{\rm B}}}{{T}_{{\rm GME}}}/\epsilon \simeq \frac{1}{2{\rm ln} (2)}$.

The previous optimization was done in two steps (first maximizing z_i and then minimizing w_i). Arguably this is not the optimal approach, as doing a bit worse in step 1 can have a global benefit. While this being true, one can easily convince themself that the differences are of $\mathcal{O}(1/n)$, and thus essentially rotating the ground state with a very excited state and then optimally permuting the rest of diagonal elements, will always lead to ${{k}_{{\rm B}}}{{T}_{{\rm GME}}}/\epsilon \simeq \frac{1}{2{\rm ln} (2)}$. Any other unitary creating X-states from thermal states can not perform better.