# Entanglement generation from athermality

A. de Oliveira Junior∗ Center for Macroscopic Quantum States bigQ, Department
of Physics, Technical University of Denmark, Fysikvej 307, 2800 Kgs. Lyngby,
Denmark  Jeongrak Son∗ School of Physical and Mathematical Sciences, Nanyang
Technological University, 637371, Singapore  Jakub Czartowski  Doctoral School
of Exact and Natural Sciences, Jagiellonian University, 30-348 Kraków, Poland
Faculty of Physics, Astronomy and Applied Computer Science, Jagiellonian
University, 30-348 Kraków, Poland  Nelly H. Y. Ng  School of Physical and
Mathematical Sciences, Nanyang Technological University, 637371, Singapore

###### Abstract

We investigate the thermodynamic constraints on the pivotal task of
entanglement generation using out-of-equilibrium states through a model-
independent framework with minimal assumptions. We establish a necessary and
sufficient condition for a thermal process to generate bipartite qubit
entanglement, starting from an initially separable state. Consequently, we
identify the set of system states that _cannot_ be entangled, when no external
work is invested. In the regime of infinite temperature, we analytically
construct this set; while for finite temperature, we provide a simple
criterion to verify whether any given initial state is or is not entanglable.
Furthermore, we provide an explicit construction of the _future thermal cone
of entanglement_ — the set of entangled states that a given separable state
can thermodynamically evolve to. We offer a detailed discussion on the
properties of this cone, focusing on the interplay between entanglement and
its volumetric properties. We conclude with several key remarks on the
generation of entanglement beyond two-qubit systems, and discuss its dynamics
in the presence of dissipation.

##  I Introduction

00footnotetext: ∗ These authors contributed equally to this work.

uantum dynamics and thermodynamics have been burgeoning areas of study,
particularly in the domain of quantum information science [1, 2, 3, 4, 5].
Upon moving from the macroscopic world to the quantum realm, it is natural to
ask how thermodynamics affects quantum information processing [6, 7, 8, 9]. A
thermodynamics-inspired paradigm for quantifying quantum correlations has led
to the conclusion that one cannot extract the maximum work from an entangled
state when it is shared by two separated parties [6, 10]. Conversely,
correlations shared between quantum systems allow for optimal work storage
[9]. In tandem with these inquiries, the fundamental energetic limitations on
creating both classical and quantum correlations have been thoroughly
investigated [11, 8, 12, 11, 13], assuming closed system dynamics.

The task of correlation engineering becomes increasingly interesting, in the
more realistic scenario where quantum systems interact with a thermal
environment, exchanging energy and entropy. Counterintuitively, it has been
shown that steady-state entanglement can be generated [14, 15, 16, 17] even
when no external agents are involved [18]. In particular, one can generate a
maximally entangled steady-state between two qubits deterministically using a
third ancillary qubit [19]. These results suggest a paradigm shift where
certain dissipative processes, rather than being detrimental to quantum
correlations, might harness environment interactions for thermodynamically
stable entanglement generation. As a result, a natural question arises – _what
are the thermodynamic limitations for quantum correlation and entanglement
generation?_

A powerful toolkit for answering this question is the _resource-theoretic_
approach [20, 21]. One may initially think that the entanglement resource
theory of local operations and classical communication (LOCC) should be the
relevant framework for the entanglement manipulation problem. However, LOCC
operations are naturally motivated only after entanglement has been
distributed; if one is instead at the preparation phase in a local lab, prior
to entanglement distribution, then it is more reasonable to consider a
resource theory that contains entangling operations. If entangling operations
are allowed, the main limiting factor must come from quantum control and
energetics [22, 23, 24, 25], which leads us to the resource theory of quantum
thermodynamics. This model-independent framework allows one to investigate
general thermodynamic processes by relying solely on the assumptions of total
energy conservation and thermality of environments. Its broad scope has
allowed us to establish bounds on the size of the bath required for the
thermalization of a many-body system [26], bridge the gap between Markovian
and non-Markovian thermal processes [27, 28], and uncover optimal heat-bath
algorithmic cooling strategies [29], among other significant results (see [30,
31] for recent reviews).

In this work, we systematically probe the capability of using thermal
processes to generate entanglement when given access to some initial state. In
the particular case of two qubits, we derive necessary and sufficient
conditions for the existence of a thermal process that transforms a given
separable initial state into a final state that is entangled. These conditions
depend only on the energetic populations of the two-qubit system. A key
insight that led to this finding, is the observation that a generic entangling
strategy can be decomposed into two independent procedures: one that is
explicitly achieved by open system dynamics (interaction with heat bath that
results in changes of energetic population), and another closed unitary
operation, i.e. an optimal rotation within the global energy subspace. The
consistency of our results is confirmed by the problem addressed in Refs. [11,
13], where the authors investigate two-qubit mixed states and determine
unconstrained global unitary operations that generate maximal entanglement. We
approach the considered problem with additional restrictions on the alllowed
operations, motivated by cornerstone findings on thermodynamic resource
theories. To be specific, we assume that entangling operations must fall into
the set of thermal operations [22, 32].

Using the above framework we provide an explicit construction of the set of
separable states that can become entangled. The degree of entanglement in this
set can be assessed either by the volume of the set or by their negativity. In
this regard, our results show that a Bell pair can always be thermodynamically
produced at any temperature regime, provided that the system is initialized in
its pure excited state. We also discuss the generation of entanglement in the
presence of dissipation, examining the behavior of the volume of the future
thermal cone of entanglement under Markovian thermal processes, and comparing
it with the closed dynamics case. Finally, we offer some remarks on how our
results can be extended to systems beyond two qubits.

The paper is organized as follows: We start by motivating how thermal
resources can be used to generate entanglement in Section II, using an
experimentally relevant example involving two spins interacting with a cavity
mode. After this warm-up, we set the scene in Section III by recalling the
resource-theoretic approach to quantum thermodynamics. In Section IV, we
rigorously define the main quantifiers used for entanglement generation under
thermal operations in arbitrary dimensions. Section V focuses on the two-qubit
system, and we derive necessary and sufficient conditions for entanglement
preparation, and proceeds to quantify the entanglement generated and discuss
its applicability. We discuss in Section VI how the task of entanglement
generation changes for continuous-time thermal processes and high-dimensional
systems, and conclude with an outlook in Section. VII.

##  II Motivating examples

Let us begin with a warm-up example that gives a flavor of the main
investigation. Consider two qubits $\mathsf{A}$ and $\mathsf{B}$, with the
former prepared in the ground state and the latter in the excited state. The
joint state of the composite system is given by
$\rho_{\mathsf{AB}}=\outerproduct{01}{01}$. Furthermore, assume that each
qubit is described by the Hamiltonian
$H_{0}=E\left|1\right\rangle\\!\\!\left\langle 1\right|$, the composite
$\mathsf{AB}$ is described by a Hamiltonian
$H=H_{0}\otimes\mathbbm{1}+\mathbbm{1}\otimes H_{0}$. This Hamiltonian gives
rise to the degenerate energy subspace
$\mathcal{V}_{E}=\mathrm{span}\\{\ket{01},\ket{10}\\}$, associated with energy
$E$. Then any unitary acting non-trivially only on $\mathcal{V}_{E}$ is
energy-preserving, including

$U=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\\\
1&-1\end{pmatrix}\oplus\mathbbm{1}_{\mathcal{H}\setminus\mathcal{V}_{E}},$ (1)

that maps product states $\ket{01}$ and $\ket{10}$ into Bell states
$\ket{\phi^{+}}$ and $\ket{\phi^{-}}$, where
$\ket{\phi^{\pm}}=\frac{1}{\sqrt{2}}(\ket{01}\pm\ket{10})$. In particular,
when the initial state is either $\ket{01}$ or $\ket{10}$, this transformation
can be implemented with Jaynes-Cummings interactions only [33]

$\displaystyle H_{\mathrm{JC}}=\sigma^{+}_{\textsf{A(B)}}\otimes
a_{\textsf{R}}+\sigma^{-}_{\textsf{A(B)}}\otimes a^{\dagger}_{\textsf{R}},$
(2)

albeit indirectly by introducing a mediating bosonic mode R [34]. We use the
notation $\sigma^{+}=\left|1\right\rangle\\!\\!\left\langle 0\right|$ and
$\sigma^{-}=\left|0\right\rangle\\!\\!\left\langle 1\right|$ for each qubit A
or B, while $a_{\textsf{R}}$ and $a_{\textsf{R}}^{\dagger}$ are the
annihilation and creation operators for the bosonic mode.

Figure 1: Negativity as a function of $\beta E$. Preprocessing the spin system
by coupling to a thermal bath opens up the possibility of generating
entanglement through energy-preserving unitaries that encode the laws of
thermodynamics. We plotted the maximal negativity [35] [defined in Eq. (31)]
obtainable by a two-step process _i)_ transforming one of the spins via
Jaynes-Cummings interaction with a thermal bosonic mode and _ii_) applying the
entangling unitary, Eq. (1). The $x$-axis is the product of inverse
environment temperature $\beta$ and energy gap $E$. The black solid curve
represents the resulting negativity from a state initially prepared in the
ground state, while the purple dashed curve is that of the initially fully
excited state.

Suppose now that both qubits are prepared in the ground state,
$\rho_{\textsf{AB}}=\outerproduct{00}{00}$. Then, it is _impossible_ to
entangle them using Eq. (1), or any other energy-preserving unitary acting on
the composite $\mathsf{AB}$, as $\rho_{\textsf{AB}}$ has no support in
$\mathcal{V}_{E}$. However, in the presence of a thermal bath, the populations
of the composite system can be modified, thereby enabling the generation of
entanglement. A simple yet experimentally relevant example illustrating this
effect involves a two-step protocol. First, assume access to a bosonic mode
with a resonant Hamiltonian
$H_{\textsf{R}}=Ea_{\textsf{R}}^{\dagger}a_{\textsf{R}}$, thermalized with
respect to an inverse temperature $\beta$. Qubit B interacts with the bosonic
mode for a finite time, which creates non-zero support of the composite state
in $\mathcal{V}_{\textsf{E}}$. Second, apply the unitary given in Eq. (1).
This two-step protocol transforms $\rho_{\textsf{AB}}$ into

$\sigma_{\textsf{AB}}=f(\beta E)\outerproduct{00}{00}+[1-f(\beta
E)]\outerproduct{\phi^{+}}{\phi^{+}}.$ (3)

The ground state population $f(\beta E)$ depends on the interaction
Hamiltonian and the duration of the interaction, in addition to the
temperature and the energy gap. Let us assume that we apply Jaynes-Cummings
interaction in Eq. (2) and optimize the time duration to minimize $f(\beta
E)$. It turns out that this optimized quantity decreases monotonically as
$\beta E$ decreases. Hence, by additionally allowing the initial state
$\rho_{\textsf{AB}}=\outerproduct{00}{00}$ to interact with a thermal bath,
one can generate entangled states; higher amounts of entanglement can be
observed at higher temperatures as well as with smaller energy gaps. The same
protocol can be applied when both spins are excited. However, for this initial
state, contrary to the previous example, $f(\beta E)$ is the smallest – and
thus the final state is more entangled – when the ambient temperature is
lower. In Fig. 1, we quantify the optimal amount of entanglement produced in
these two cases, measured via the negativity, as a function of the inverse
temperature of the bosonic mode.

In the above discussion, we have constructed specific thermal processes that
precondition the composite state into entanglable ones, by using a simple
bosonic mode and the Jaynes-Cummings interaction. A natural question follows:
_what is the ultimate bound of the entanglement generation, for the most
generic, all-encompassing set of thermal processes?_

##  III Setting the scene

We investigate entanglement generation in an $N$-partite system whose initial
state and the non-interacting Hamiltonian are given as

$\rho\in\mathrm{conv}\quantity(\bigotimes_{i=1}^{N}\rho_{i})\quad\text{and}\quad
H=\sum_{i=1}^{N}\mathbbm{1}^{\otimes(i-1)}\otimes
H_{i}\otimes\mathbbm{1}^{\otimes(N-i)}.$ (4)

The assumption that the Hamiltonian is non-interacting is natural, and
excludes the trivial scenario of entangling simply via a fully thermalizing
process. The Hilbert space of the composite system is denoted by
$\mathcal{H}=\bigotimes_{i=1}^{N}\mathcal{H}_{i}$. The energy level structure
creates a natural division of the full Hilbert space into energy-degenerate
subspaces $\mathcal{V}_{E}$,

$\mathcal{H}=\bigoplus_{E}\mathcal{V}_{E},$ (5)

where
$\mathcal{V}_{E}=\mathrm{span}\left\\{\ket{\psi}\in\mathcal{H}:H\ket{\psi}=E\ket{\psi}\right\\}$.
Note that energy-preserving unitaries acting on $\mathcal{H}$ are those that
can be written as $U=\bigoplus_{E}U_{E}$, where each unitary $U_{E}$ acts on
$\mathcal{V}_{E}$. In other words, they do not mix states from different
energy subspaces – an essential feature for subsequent considerations.

###  III.1 Thermal operations

The multipartite system interacts with the thermal environment that is
initialized in a Gibbs state

$\gamma_{\textsf{R}}=\frac{e^{-\beta H_{\textsf{R}}}}{\Tr(e^{-\beta
H_{\textsf{R}}})},$ (6)

with a Hamiltonian $H_{\textsf{R}}$ and an inverse temperature $\beta$. We
make minimal assumptions on the joint system-bath dynamics: the composite
system is closed and thus evolves unitarily; in addition, this unitary
evolution is energy-preserving. Formally, this leads to a set of operations
known as _thermal operations_ (TOs) [22], which transform the state $\rho$ as
follows:

$\mathcal{E}(\rho)=\Tr_{\textsf{R}}\left[U\left(\rho\otimes\gamma_{\textsf{R}}\right)U^{\dagger}\right],$
(7)

where $U$ is a joint unitary that commutes with the total Hamiltonian of the
system and the bath

$[U,H\otimes\mathbbm{1}_{\textsf{R}}+\mathbbm{1}\otimes H_{\textsf{R}}]=0,$
(8)

and the environment Hamiltonian $H_{\textsf{R}}$ can be chosen arbitrarily.
Given that there are no further constraints on $U$, the system and the bath
can develop arbitrarily strong correlations (including entanglement).

###  III.2 State transformation conditions

We focus on _energy-incoherent_ states $\rho$, which means that each partition
is (block) diagonal in the energy eigenbasis, resulting in an output state
that is also (block) diagonal. Mathematically, the multipartite state $\rho$
commutes with the Hamiltonian, i.e., $[\rho,H]=0$. Under this assumption, the
state $\rho$ can be fully represented by a probability vector $\boldsymbol{p}$
consisting of the eigenvalues of $\rho$, which correspond to populations in
the energy eigenbasis, $\rho=\sum_{i}p_{i}\outerproduct{E_{i}}{E_{i}}$. We
henceforth directly refer to $\boldsymbol{p}=(p_{1},...,p_{d})$ as energy-
incoherent states, which live within the probability simplex

$\Delta_{d}=\quantity{\boldsymbol{p}=(p_{1},...,p_{d})\in\mathbbm{R}_{\geq
0}^{d}:\sum_{i}p_{i}=1}.$ (9)

Analogously, the thermal equilibrium state of the system, Eq. (6) with
$H_{\textsf{R}}$ replaced by $H$, can be represented by a vector of thermal
populations, denoted as $\boldsymbol{\gamma}$.

When considering transformations between different energy-incoherent states
under thermal operations, an important aspect of a state is its associated
$\beta$-order – a permutation which arranges populations, rescaled with
respect to the Gibbs state, in a non-increasing order.

######  Definition 1 ($\beta$-ordering).

Suppose that population vectors $\boldsymbol{p}$ and $\boldsymbol{\gamma}$,
corresponding to a $d$-dimensional system state and its Gibbs state with
respect to the ambient temperature $1/\beta$, are given. The $\beta$-ordering
of $\boldsymbol{p}$ is given as a permutation
$\pi_{\boldsymbol{p}}\in\mathcal{S}_{d}$ of energy levels within the symmetric
group $\mathcal{S}_{d}$, such that the element-wise ratio between
$\boldsymbol{p}$ and $\boldsymbol{\gamma}$ is non-increasing, i.e.

$\displaystyle\frac{p_{(\pi_{\boldsymbol{p}})i}}{\gamma_{(\pi_{\boldsymbol{p}})i}}\geq\frac{p_{(\pi_{\boldsymbol{p}})j}}{\gamma_{(\pi_{\boldsymbol{p}})_{j}}},\quad\text{for
all}\quad i\leq j.$ (10)

For simplicity, we also denote a $\beta$-ordering $\pi$ as a vector
$\boldsymbol{\pi}=(\pi(1),\pi(2),\cdots,\pi(d))$ in the rest of the
manuscript. For the specific case of $d=4$, we denote the $\beta$-ordering
$(2,1,3,4)$ as $\boldsymbol{\pi}^{\star}$, since it plays a special role in
entanglement generation, which will become clear in Section IV (Theorem 2).

Finally, the set of states $\mathcal{T}_{+}(\boldsymbol{p})$ that can be
achieved via thermal operations from a given initial state $\boldsymbol{p}$ is
called the _future thermal cone_ [36, 37]. This set is a convex polytope
comprising at most $d!$ extreme points, each corresponding to different
$\beta$-orderings $\boldsymbol{\pi}$ [36, 37]. Mathematically, it is
represented by:

$\mathcal{T}_{+}(\boldsymbol{p})=\operatorname{conv}[{\boldsymbol{p}^{\boldsymbol{\pi}},\boldsymbol{\pi}\in{\cal
S}_{d}}],$ (11)

where $\boldsymbol{p}^{\boldsymbol{\pi}}$ are extreme points characterized by
thermomajorization [38] (see Lemma 5, Appendix A). Again, for $d=4$, we denote
the extreme point corresponding to the ordering $\boldsymbol{\pi}^{\star}$ as
$\boldsymbol{p}^{\star}$.

Although not explicitly stated, it is important to mention that the existence
of a thermal operation that induces the transition between incoherent states
$\boldsymbol{p}\rightarrow\boldsymbol{p}^{\prime}$ is equivalent to the
existence of a Gibbs-preserving (GP) stochastic matrix mapping between these
states [22]. Therefore, when studying the thermodynamic interconversion
problem between block diagonal states, one can focus on stochastic matrices
and probability vectors instead of CPTP maps and density matrices.
Nonetheless, our forthcoming discussion will primarily center on the concept
of thermal cones, rather than on GP matrices or thermal operations.

##  IV Entanglability under thermal operations

In this section, we introduce the main definitions that will enable us to
derive the necessary and sufficient conditions for entangling a separable
state under thermal operations. At first glance, this might appear abstract
due to its generality, but these notions will be applied to a specific case in
Section V.

We begin by defining a subset of unitary operations that preserve the energy
of the system and keep the energy level structure intact. Note that this is
motivated by capturing the thermodynamic constraints in producing entanglement
when the system is closed, i.e. when the unitary dynamics satisfy energy
conservation. We refer to this subset as _energy subspace-unitaries_ , denoted
by $\mathcal{U}_{\mathrm{subs}}$. Formally, this set is given by

$\mathcal{U}_{\mathrm{subs}}=\quantity{U_{\mathrm{subs}}=\bigoplus_{E}U_{E}},$
(12)

where $U_{E}$ acts on $\mathcal{V}_{E}$ subspace of definite energy $E$.
Consequently, any unitary operation from this set is an admissible operation
within the framework of thermal operations.

Next, we introduce an entanglement measure $\mathcal{M}$, which is non-
negative for all separable states, $\mathcal{M}(\rho_{\textsf{sep}})\geq 0$.
One such entanglement measure which we focus on is derived from the famous
Horodecki-Peres positive partial transpose (PPT) criterion [39, 40], which
states that separable states have positive partial transpose,
$\rho^{T_{\textsf{A}}}\geq 0$; using this, we define
$\mathcal{M}(\rho)\coloneqq\min_{\ket{\psi}}\expectationvalue{\rho^{T_{A}}}{\psi}=\lambda_{-}$
as the minimal eigenvalue of the partial transpose. Therefore, the quantity

$f(\boldsymbol{p})\coloneqq\min_{U\in\mathcal{U}_{\mathrm{subs}}}\mathcal{M}\quantity(U\rho
U^{\dagger})$ (13)

is a good witness for the thermodynamic capacity of creating entanglement from
$\boldsymbol{p}$ via closed, energy-conserving system dynamics.

######  Definition 2 (Subspace (non-)entanglable set).

Given a state $\boldsymbol{p}\in\Delta_{d}$, the subspace non-entanglable set
is defined as the set of states that remain PPT under subspace rotations, i.e.

$\mathbbm{NE}\coloneqq\quantity{\boldsymbol{p}\,\big{|}\,f(\boldsymbol{p})\geq
0}.$ (14)

The function $f$ acts as an entanglement witness, i.e. $f(\boldsymbol{p})\leq
0$ implies the possibility of creating entanglement. Furthermore, the
_subspace entanglable set_ is defined as the complement of $\mathbbm{NE}$,
i.e. $\mathbbm{E}\coloneqq\Delta_{d}\setminus\mathbbm{NE}$.

This definition can be extended by employing multiple entanglement measures
$\mathcal{M}_{1},\ldots,\mathcal{M}_{M}$, resulting in multiple witnesses
$f_{1},\ldots,f_{M}$, and requiring that all these functions yield positive
values.

It is important to note that, states in the set $\mathbbm{E}$ are entanglable
even when neither heat or work exchange are allowed. Notably, when thermal
resources are available, i.e. under thermal operations, the set of entanglable
states expand.

######  Definition 3 (Future thermal cone of entanglement).

The future thermal cone of entanglement of a state $\boldsymbol{p}$ is defined
as

$\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{p})\coloneqq\quantity{\boldsymbol{q}\,\big{|}\,f(\boldsymbol{q})<0,\
\boldsymbol{q}\in\mathcal{T}_{+}(\boldsymbol{p})}=\mathcal{T}_{+}(\boldsymbol{p})\cap\mathbbm{E}.$
(15)

Now we introduce the primary object of this paper: the set of states that
cannot become entangled under thermal operations.

######  Definition 4 (Thermally non-entanglable set).

Given the ambient temperature $\beta^{-1}$, the thermally non-entanglable set
is defined as

$\mathbbm{TNE}(\beta)=\quantity{\boldsymbol{p}\,\big{|}\,\min_{\boldsymbol{q}\in\mathcal{T}_{+}(\boldsymbol{p})}f(\boldsymbol{q})\geq
0},$ (16)

with the entanglement measure $\mathcal{M}$. Analogous to $\mathbbm{E}$, we
also define the set of thermally entanglable states
$\mathbbm{TE}(\beta)\coloneqq\Delta_{d}\setminus\mathbbm{TNE}(\beta)$.

By definition, $\boldsymbol{p}$ is thermally non-entanglable if and only if
its future cone of entanglement is empty, i.e.

$\boldsymbol{p}\notin\mathbbm{TNE}(\beta)~{}~{}\Longleftrightarrow~{}~{}\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{p})=\varnothing.$
(17)

We can then characterize $\mathbbm{TNE}(\beta)$ with the following
proposition. Hereafter, we drop the explicit dependence on $\beta$ in
$\mathbbm{TNE}$ and $\mathbbm{TE}$ when there is no concern for confusion.

######  Proposition 1 (Thermally non-entanglable set).

Suppose that the entanglement measure $\mathcal{M}$ is a concave function.
Then the set $\mathbbm{TNE}$ can be written as

$\mathbbm{TNE}=\quantity{\boldsymbol{p}\mid\forall_{\pi}f(\boldsymbol{p}_{\pi})\geq
0}.$ (18)

When $\beta E=0$, the above expression simplifies to

$\mathbbm{TNE}=\bigcap_{\pi\in\mathcal{S}_{d}}\quantity{\boldsymbol{p}\mid
f(\Pi\boldsymbol{p})\geq
0}=\vcentcolon\bigcap_{\pi\in\mathcal{S}_{d}}\,\Omega_{\pi}\ ,$ (19)

where $\Pi$ represents a permutation matrix corresponding to
$\pi\in\mathcal{S}_{d}$; in particular, for $\pi=\mathrm{id}$ we have
$\Omega_{\mathrm{id}}=\mathbbm{NE}$.

Definition 3 also allows us to quantify the entanglement generation power
using geometric measures. The volume of the future thermal cone of
entanglement can gauge how useful an initial state is for producing entangled
states. Given a state $\boldsymbol{p}$, this quantity is defined as

$\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{p})\coloneqq\frac{V\left[\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{p})\right]}{V(\Delta_{d})},$
(20)

where $V$ is the volume measured using the Euclidean metric. Likewise, the
volume of (non-)entanglable sets, defined in Definitions 2 and 4, are denoted
as

$\mathbf{V}_{\mathbbm{X}}\coloneqq\frac{V({\mathbbm{X}})}{V(\Delta_{d})}\quad,\quad\mathbbm{X}=\mathbbm{NE},\
\mathbbm{E},\ \mathbbm{TNE}.$ (21)

##  V Two-qubit systems

Let us now apply the general framework we developed to a specific scenario:
two-qubit entanglement. In this section, we find an efficient criterion for
determining whether a two-qubit state $\boldsymbol{p}$ is thermally non-
entanglable, which only checks if $\boldsymbol{p}^{\star}$ is subspace non-
entanglable. We also study the geometry of the future thermal cone of
entanglement and the thermally non-entanglable set for different ambient
temperatures.

###  V.1 Subspace entanglable set for two qubits

Assume that Alice and Bob share a pair of qubits that are initially prepared
in a separable state $\rho_{\textsf{AB}}$. Each of them is described by a
local Hamiltonian $H_{A}=H_{B}=E\outerproduct{1}{1}$, and they are allowed to
interact with a single heat bath via thermal operations. If there exists a
thermal operation that maps $\rho_{\textsf{AB}}$ to $\sigma_{\textsf{AB}}$,
where $[\rho_{\textsf{AB}},H]=0$ and $[\sigma_{\textsf{AB}},H]=0$, then
$\sigma_{\textsf{AB}}$ can be expressed as

$\sigma_{\textsf{AB}}=q_{1}\outerproduct{00}{00}+q_{2}\outerproduct{\psi}{\psi}+q_{3}\outerproduct{\psi_{\perp}}{\psi_{\perp}}+q_{4}\outerproduct{11}{11},$
(22)

where $\ket{\psi},\lvert\psi^{\perp}\rangle\in\mathcal{V}_{E}$ and
$\langle\psi|\psi^{\perp}\rangle=0$. This observation follows from the fact
that one can freely rotate the populations within the degenerate energy
subspace. The unitary that connects $\ket{\psi}=U\ket{10}$ and
$\ket{\psi_{\perp}}=U\ket{01}$ can be written with two parameters
$\theta,\phi\in[0,2\pi]$

$U=\begin{pmatrix}\cos\theta&\sin\theta e^{i\phi}\\\ -\sin\theta
e^{-i\phi}&\cos\theta\end{pmatrix}\oplus\mathbbm{1}_{\mathcal{H}\setminus\mathcal{V}_{E}},$
(23)

apart from an irrelevant phase factor, global to $\mathcal{V}_{E}$. Then the
final state $\sigma$ reads

$\\!\\!\sigma=\begin{pmatrix}q_{1}&0&0&0\\\
0&q_{2}\cos^{2}\theta+q_{3}\sin^{2}\theta&\frac{(q_{3}-q_{2})}{2}\sin 2\theta
e^{i\phi}&0\\\ 0&\frac{(q_{3}-q_{2})}{2}\sin 2\theta
e^{-i\phi}&q_{2}\sin^{2}\theta+q_{3}\cos^{2}\theta&0\\\
0&0&0&q_{4}\end{pmatrix}.$ (24)

For a two-qubit system, separability can be fully determined by the positive
partial transposition (PPT) criterion [39, 40]. When the partial transposition
is applied to the state $\sigma$ in Eq. (24), its smallest eigenvalue is

$\\!\\!\\!\lambda_{-}=\frac{1}{2}\left[\left(q_{1}+q_{4}\right)-\sqrt{\left[(q_{2}-q_{3})^{2}\sin^{2}2\theta+(q_{1}-q_{4})^{2}\right]}\right],$
(25)

which is minimized when $\sin^{2}2\theta=1$111This minimum value coincides
with the case where we diagonalize the energy degenerate subspace in Bell
basis, i.e., $\ket{\psi},\lvert\psi_{\perp}\rangle=\lvert\Psi^{\pm}\rangle$..
Thus, $\lambda_{-}$ can become negative if and only if

$4q_{1}q_{4}<(q_{2}-q_{3})^{2}.$ (26)

In other words, a state with population vector
$\boldsymbol{q}=(q_{1},q_{2},q_{3},q_{4})$ is subspace entanglable if and only
if Eq. (26) is true. We will refer to this inequality as the entanglement
constraint or criterion. In addition, we define an associated entanglement
measure

$f(\boldsymbol{q}):=4q_{1}q_{4}-(q_{2}-q_{3})^{2},$ (27)

which detects entanglement whenever $f(\boldsymbol{q})<0$.

Figure 2: Subspace entanglable set $\mathbbm{E}$. Geometry of the set of
states that can get entangled via energy-preserving unitaries acting on the
degenerate energy subspace. Note that the set $\mathbbm{E}$ consists of two
disjoints sets each on the opposite side of the plane $p_{2}=p_{3}$ coloured
in blue. Figure 3: Thermally non-entanglable states $\mathbbm{TNE}$ for $\beta
E=0$. States that belong to the coloured geometrical shape cannot get
entangled under thermal operations when $\beta E=0$.

Figure 2 depicts the subspace entanglable set $\mathbbm{E}$ (see Figure 10 in
Appendix B for its complement $\mathbbm{NE}$). States in $\mathbbm{E}$ can
become entangled by some energy-preserving unitary of the form Eq. (23).

The sets $\mathbbm{NE}$ and $\mathbbm{E}$ are independent of the ambient
temperature and their volume can be calculated by numerical methods or direct
integration. Remarkably, the normalised volume of $\mathbbm{E}$ is
approximately $\mathbf{V}_{\mathbbm{E}}\approx 21/62$; much less than half of
the states in $\Delta_{4}$ are subspace entanglable. In contrast, as we note
in Section V-C, the future thermal cone of entanglement of the highest excited
state $\boldsymbol{p}=(0,0,0,1)$ is precisely $\mathbbm{E}$, as this state can
achieve any other state in the probability simplex, regardless of the ambient
temperature. Hence, the full volume $\mathbf{V}_{\mathbbm{E}}$ can be achieved
via thermal operations. Futhermore, we prove the convexity of the subspace
non-entanglable set, find its extreme points, and identify a continuous family
of boundary points (Appendix B).

###  V.2 Thermally entanglable states

We now construct the set of states that can generate entanglement with the
help of a thermal environment. In the simple case of $\beta E=0$, its
complement, the set of thermally non-entanglable states $\mathbbm{TNE}$, can
be immediately constructed from Proposition 1. The resulting set is convex and
occupies approximately one third of the probability simplex [see Fig. (3)].
However, the structure of this set becomes more intricate when $\beta E\neq
0$. Nevertheless, we simplify the problem by showing that the task of
determining whether $\boldsymbol{p}\in\mathbbm{TNE}$ – defined by Equation
(16) as including a non-trivial minimization – is reduced to verifying if
$\boldsymbol{p}^{\star}\in\mathbbm{NE}$, which only requires checking a single
inequality $f(\boldsymbol{p}^{\star})\geq 0$.

######  Theorem 2 (Criterion for determining $\boldsymbol{p}\in\mathbbm{TE}$).

Let $\boldsymbol{p}$ be a two-qubit state in $\mathbbm{TE}$. Then
$\boldsymbol{p}^{\star}$, the extreme point of
$\mathcal{T}_{+}(\boldsymbol{p})$ corresponding to the $\beta$-ordering
$\boldsymbol{\pi}^{\star}=(2,1,3,4)$, is in $\mathbbm{E}$.

The full proof is detailed in Appendix C, and consists in two main steps.
First, we identify $\boldsymbol{\pi}^{\star}$ as a special $\beta$-ordering
for subspace entanglability: whenever $\boldsymbol{p}\in\mathbbm{TNE}$, there
exists a state
$\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{p})$ with the
ordering $\boldsymbol{\pi}^{\star}$ (Lemma 7). Next, we use a geometric
argument (Proposition 8) to prove that whenever there exists such
$\boldsymbol{q}$, the extreme point $\boldsymbol{p}^{\star}$ is also subspace
entanglable, i.e.
$\boldsymbol{p}^{\star}\in\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{p})$. It
is however important to note that Theorem 2 does not imply that the minimum of
$f(\boldsymbol{q})$ over $\boldsymbol{q}\in\mathcal{T}_{+}(\boldsymbol{p})$ is
always achieved when $\boldsymbol{q}=\boldsymbol{p}^{\star}$. For a simple
counter-example, consider a state
$\boldsymbol{p}=(\epsilon,1-\epsilon-\epsilon^{2},0,\epsilon^{2})$ with
$\epsilon\ll e^{-2\beta E}$. The state $\boldsymbol{p}$ has a $\beta$-ordering
$(2,1,4,3)$, yet a smaller $f$ value
$f(\boldsymbol{p})<f(\boldsymbol{p}^{\star})$ compared to
$\boldsymbol{p}^{\star}$.

Figure 4: Optimal entanglement generation regime. The upper plot illustrates
the approximated log volume of the set of states that cannot become entangled,
as a function of $\beta$. The bottom plot demonstrates how the geometry of
this set changes at specific values of $\beta$.

The interplay between the temperature $\beta^{-1}$ of the immediate
environment and the ability to entangle a system manifests in Theorem 2
through the $\beta$-ordering. We numerically quantify the volume of the
thermally non-entanglable state set $\mathbbm{TNE}$ in Fig. 4. For any finite
temperature regime $\beta>0$, determining all states that are thermally non-
entanglable is not straightforward. We devise a computationally efficient
algorithm, using Theorem 2 and the vacuum-packing-inspired procedure, that
approximates $\mathbbm{TNE}$ for given $\beta$:

  1. 1.

Generate a regular grid on the surface of the probability simplex
$\Delta_{3}$.

  2. 2.

Select a point $\boldsymbol{p}^{(0)}_{o}$ from the grid and the Gibbs state
$\boldsymbol{p}^{(0)}_{i}=\boldsymbol{\gamma}$. With this choice, we are
guaranteed to have $\boldsymbol{p}^{(0)}_{o}\in\mathbbm{TE}$ and
$\boldsymbol{p}^{(0)}_{i}\in\mathbbm{TNE}$222Except for a single case
$\boldsymbol{p}^{(0)}_{o}=\frac{1}{1+2e^{-\beta E}}(1,e^{-\beta E},e^{-\beta
E},0)$..

  3. 3.

In the $j$-th step of the procedure we define midpoint
$\boldsymbol{p}^{(j)}_{m}=\frac{1}{2}\quantity(\boldsymbol{p}^{(j)}_{i}+\boldsymbol{p}^{(j)}_{o})$
and, based on Theorem 2, decide whether $\boldsymbol{p}^{(j)}_{m}$ is an
element from the non-entanglable set. If
$\boldsymbol{p}^{(j)}_{m}\in\mathbbm{TNE}$, then
$\boldsymbol{p}^{(j+1)}_{i}=\boldsymbol{p}^{(j)}_{m}$; else
$\boldsymbol{p}^{(j+1)}_{o}=\boldsymbol{p}^{(j)}_{m}$.

  4. 4.

Repeat the bisection step $n=30$ times, achieving the precision
$\norm{\boldsymbol{p}^{(n)}_{i}-\boldsymbol{p}^{(n)}_{o}}\approx
10^{-9}\norm{\boldsymbol{p}^{(0)}_{i}-\boldsymbol{p}^{(0)}_{o}}$.

  5. 5.

Add the point $\boldsymbol{p}^{(n)}_{m}$ to the approximate set
$\widetilde{\mathbbm{TNE}}$.

  6. 6.

Repeat steps 2-5 for all the gridpoints.

The convex hull of the final set
$\operatorname{conv}(\widetilde{\mathbbm{TNE}})$ is a good approximation of
the actual non-entanglable set $\mathbbm{TNE}$.

The numerical result in Fig. 4 exhibits a monotonic behaviour, indicating that
low temperatures enhance entanglability. Indeed, this outcome is intuitively
expected. Thermal processes permit the local bypassing of strict energy-
preservation, thereby offering advantage in terms of entanglement generation.
Yet, thermal noise inevitably accompanies thermal processes, impacting the
amount of achievable entanglement. The negative impact of thermal noise is
anticipated to be minimized at lower environmental temperatures. We confirm
this first by observing that the set of non-entanglable states $\mathbbm{TNE}$
vanishes at the zero temperature limit $\beta\rightarrow\infty$. Moreover, we
can make a more rigorous observation.

######  Lemma 3.

The set $\mathbbm{TNE}(\infty)=\\{(1,0,0,0)\\}$, i.e. it contains only the
ground state.

###### Proof.

Given any other initial state, it is possible to thermalize one of the
subsytems; in the zero-temperature limit ($\beta\rightarrow\infty$), this
transforms a state $\boldsymbol{p}=(p_{1},p_{2},p_{3},p_{4})$ into
$\boldsymbol{q}=(p_{1}+p_{3},p_{2}+p_{4},0,0)$ or
$(p_{1}+p_{2},p_{3}+p_{4},0,0)$. However, in this specific scenario, it is
sufficient to assume entanglement because the last two levels are empty while
the second level is populated. Consequently, $f$ always becomes negative.
Thus, for any state $\boldsymbol{p}\neq(1,0,0,0)$, the final state
$\boldsymbol{q}$ is subspace entanglable, i.e.
$\boldsymbol{p}\in\mathbbm{TE}$. ∎

One may be tempted to view the ground state $(1,0,0,0)$ then as useless for
entanglement generation; however, this is not true. In particular, the ground
state becomes thermally entanglable whenever $\beta^{-1}\neq 0$: one can
thermalize one of the subsystems to obtain
$\boldsymbol{q}=\frac{1}{1+e^{-\beta E}}(1,e^{-\beta E},0,0)$, which is
subspace entanglable. What is observed is that as long as there is a certain
amount of athermality in the initial state (with respect to environment
temperature), we can often trade it for some amount of entanglement.

Given the above, it is then intuitive to see that given an initial state
$\boldsymbol{p}$, the range of environmental temperature that allows for
entanglement generation is very much state-dependent. In particular, we are
able to identify critical inverse temperatures $\beta_{C}$, at which a state
$\boldsymbol{p}$ starts to become thermally (non-)entanglable. In fact, we
show that $\beta_{C}$, if exists, can be directly calculated for any state by
employing Theorem 2 alongside the entanglement constraint in Eq. (26).
Interestingly, for some states two critical temeperatures
$\beta_{C_{1}}<\beta_{C_{2}}$ exist and the state becomes entanglable when
$\beta<\beta_{C_{1}}$ or $\beta>\beta_{C_{2}}$.

A particularly clear analysis can be done when the initial state themselves
are thermal, having a temperature $\beta_{\textsf{S}}\neq\beta$ in principle
different from the ambient temperature. We report the conclusions here, while
directing the curious reader to the full calculation in Appendix D.

Thermal initial states are written as

$\boldsymbol{p}=\frac{1}{Z_{\textsf{S}}}(1,e^{-\beta_{\textsf{S}}E},e^{-\beta_{\textsf{S}}E},e^{-2\beta_{\textsf{S}}E}),$
(28)

where $Z_{\textsf{S}}=(1+e^{-\beta_{\textsf{S}}E})^{2}$. When the system is
colder than the environment, $\beta_{\textsf{S}}>\beta$, the $\beta$-ordering
of the state is $(1,2,3,4)$. Finding the explicit form of the extreme point
$\boldsymbol{p}^{\star}$ and solving $f(\boldsymbol{p}^{\star})=0$, we may
obtain the critical temperature

$\displaystyle\beta_{C_{1}}\equiv\beta_{\textsf{S}}-E^{-1}\log\left[1+2e^{-\beta_{\textsf{S}}E}\left(\sqrt{e^{2\beta_{\textsf{S}}E}+1}-1\right)\right]\
.$ (29)

Whenever the ambient temperature $\beta^{-1}$ is higher than the critical
temperature $\beta_{C_{1}}^{-1}$, the thermal state $\boldsymbol{p}$ becomes
thermally entanglable. We can obtain a simple approximation
$\beta_{C_{1}}\simeq\beta_{\textsf{S}}-E^{-1}\log 3$ of Eq. (29) up to the
zeroth order of $e^{-\beta_{\textsf{S}}E}$ in the low system temperature limit
$\beta_{\textsf{S}}E\gg 1$.

On the other hand, when the system is hotter than its surrounding, i.e.
$\beta_{\textsf{S}}<\beta$, the critical temperature can be obtained from a
quartic equation and we do not have an explicit formula for it. However, in
the same limit $\beta_{\textsf{S}}E\gg 1$, we can approximate
$\beta_{C_{2}}\simeq\beta_{\textsf{S}}+E^{-1}\log 3$. In other words, if

$|\beta_{S}-\beta|\geq E^{-1}\log 3,$ (30)

then we may operate such a thermal machine to generate entanglement. See Fig.
5 for the numerical values of the critical temperatures, which clearly
displays the validity of our approximation
$\beta_{C}\simeq\beta_{\textsf{S}}\pm E^{-1}\log 3$.

As mentioned earlier, the temperature difference between the system and the
environment is a resource for entanglement generation, as it is for many other
thermodynamic tasks. This process can be conceptualized as a type of heat
engine which aims not to output work, but entanglement. For instance, consider
a 3-stroke thermal-entanglement engine whose working substance is a two-qubit
system. This engine operates between an ambient bath at temperature
$\beta^{-1}$ and a working bath with temperature
$\beta^{-1}_{W}:\absolutevalue{\beta^{-1}_{W}-\beta^{-1}_{C}}>0$. The cycle
proceeds as follows: _(i)_ the two-qubit system thermalises with the working
bath, reaching equilibrium at the working bath’s temperature $\beta_{W}^{-1}$;
_(ii)_ the system is then brought into contact with the ambient bath, which is
at temperature $\beta^{-1}$. During this stage, entanglement is generated
within the two-qubit system due to thermal interactions with the ambient bath.
To close the cycle, _(iii)_ the generated entanglement is used, and the system
is returned to its initial state, which is energy-incoherent.

Figure 5: Thermal state & entanglability. Critical temperatures of an
environment for thermal initial states, where the qubit energy gap is set to
$E=1$. The $x$-axis indicates the inverse temperature $\beta_{\textsf{S}}$ of
the initial state and the $y$-axis marks the ambient inverse temperature
$\beta$. If the pair $(\beta_{\textsf{S}},\beta)$ is in either blue or red
shaded region, the thermal initial state is entanglable. The blue and red
lines mark the hot critical (inverse) temperature and cold critical (inverse
temperature) for the environment to make the thermal initial state
entanglable.

###  V.3 Quantifiers for achievable entanglement under thermal operations

In the previous sections, we focused on determining whether a state is
entanglable or not using thermal resources. Here, we are additionally
interested in the extent of entanglement that thermal operations can produce,
and how this varies across different temperature regimes. Two figures of merit
can be used to quantify this extent: the _thermal volume of entanglement_ ,
$\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{p})$ as defined in Eq.
(20) and the maximum negativity [35] that the initial state can produce. In
particular, the thermal volume of entanglement is a thermodynamic monotone, in
other words,
$\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{p})\geq\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{q})$
for all $\boldsymbol{q}\in\mathcal{T}_{+}(\boldsymbol{p})$. Further insights
on this quantity can be developed by considering the geometry of the future
thermal cone of entanglement $\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{p})$.
Recall in Definition 3, $\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{p})$ is
the intersection of two sets: $\mathbbm{E}$ and
$\mathcal{T}_{+}(\boldsymbol{p})$. We also know that $\mathbbm{E}$ is non-
convex and not connected, hence
$\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{p})$ is also non-convex and not
connected [see bottom panel of Fig. (6)]. However, as the temperature
decreases, two disjoint parts grow in sizes and approach each other, becoming
almost connected; this is a natural consequence of the future thermal cone
approaching the plane $p_{4}=0$, which is on the boundary of $\mathbbm{E}$.

Figure 6: Regions of entanglement. The evolution of the future and
entanglement thermal cones as a function of $\beta$ reveals interesting
behavior. For a bipartite product state
$\boldsymbol{p}=(0.5,0.5)\otimes(0.24,0.76)$ with an identical energy gap
$E=1$, we observe that at $\beta=0$, the state cannot be entangled via thermal
operations. However, for $\beta\geq\beta_{C}\approx 0.21$, the state becomes
entanglable. As $\beta$ increases, the two disjoint parts of
$\mathcal{T}^{\mathrm{ENT}}(\boldsymbol{p})$ become closer, until at
$\beta\to\infty$, they are separated only by a plane $p_{2}=p_{3}$.

The change in volume $\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{p})$
with respect to $\beta$ is, however, heavily dependent on the initial state.
To see this, consider the energy eigenstates, i.e. four vertices of the
probability simplex $\Delta_{4}$ as described explicitly in the caption of
Figure 7. For the ground state $\boldsymbol{p}_{1}$,
$\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{p}_{1})$ decreases with
$\beta$ [black solid curve, top panel of Fig. 7], because the future thermal
cone $\mathcal{T}_{+}(\boldsymbol{p}_{1})$ shrinks as the Gibbs state $\gamma$
approaches the ground state as $\beta$ increases. On the other hand, for the
maximally excited state $\boldsymbol{p}_{4}$,
$\mathcal{T}_{+}(\boldsymbol{p}_{4})$ remains constant [red dot curve, top
panel of Fig. 7], since the future thermal cone encompasses the entire
probability simplex, regardless of the ambient temperature. In this case, the
future thermal cone of entanglement is precisely $\mathbbm{E}$. Finally, for
both $\boldsymbol{p}_{2}$ and $\boldsymbol{p}_{3}$, the thermal volume of
entanglement decreases with $\beta$, albeit slower than that of
$\boldsymbol{p}_{1}$.

Figure 7: Volume and negativity as a function of $\beta$. For all sharp states
$\boldsymbol{p}_{1}=(1,0,0,0)$, $\boldsymbol{p}_{2}=(0,1,0,0)$,
$\boldsymbol{p}_{3}=(0,0,1,0)$, and $\boldsymbol{p}_{4}=(0,0,0,1)$, we plot
the thermal volume of entanglement (top), and the largest negativity as a
function of the inverse temperature. Among all states, two are distinct: when
both qubits start in the excited state $\boldsymbol{p}_{4}$ and when they
start in the ground state (black solid curve) $\boldsymbol{p}_{1}$. Any other
permutation of the initial state characterizes a different active state, and
their volume is represented by the blue dotted curve on the top panel. Note
that this plot also represents the normalized volume of subspace entanglable
states, since the thermal cone of $\boldsymbol{p}_{4}$ is the full set space.
Moreover, although the maximum achievable negativity of $\boldsymbol{p}_{4}$
is the same as, e.g., $(0,1,0,0)$, it is much more powerful in terms of state
preparation, which can be useful since entanglement is in general irreversibly
manipulated. In the bottom panel, the maximum negativity of these states
coincides with the maximally active state.

The second figure of merit is negativity [35],

$\mathcal{N}(\rho)=\frac{\norm{\rho^{T_{\textsf{A}}}}_{1}-1}{2}=\sum_{i}\frac{|\lambda_{i}|-\lambda_{i}}{2},$
(31)

where $\lambda_{i}$ represents all the eigenvalues of $\rho^{T_{\textsf{A}}}$.
Equivalently, $\mathcal{N}(\rho)=\sum_{j}|{\lambda}^{-}_{j}|$, where
${\lambda}^{-}_{j}$ are negative eigenvalues of $\rho^{T_{\textsf{A}}}$.
Therefore, if the system is qubit-qubit or qubit-qutrit, non-zero negativity
is necessary and sufficient for the state being entangled; for any larger
systems, there exist entangled states with zero negativity [41, 42]. In the
qubit-qubit system we are examining, given a population vector
$\boldsymbol{q}=\mathrm{eig}(\sigma)$, the maximal negativity can be obtained
from Eq. (25) as

$\mathcal{N}(\sigma)=\frac{1}{2}\quantity[\sqrt{(q_{1}-q_{4})^{2}+(q_{2}-q_{3})^{2}}-(q_{1}+q_{4})],$
(32)

when $f(\boldsymbol{q})<0$, which has also been expected from Ref. [11] in a
different context. Observe that Eq. (32) is related to Eq. (27) via

$\displaystyle\mathcal{N}(\sigma)=\frac{1}{2}\quantity[\sqrt{(q_{1}+q_{4})^{2}-f(\boldsymbol{q})}-(q_{1}+q_{4})],$
(33)

and monotonically decreases as $f(\boldsymbol{q})$ increases. For
$f(\boldsymbol{q})\geq 0$, the negativity is always zero.

The maximum negativity over all probability simplex $\Delta_{4}$ occurs when
$\boldsymbol{q}=(0,1,0,0)$ or $\boldsymbol{q}=(0,0,1,0)$ and the corresponding
state are the Bell states. If the system is initially in one of the four
vertices of $\Delta_{4}$, it is possible to prepare Bell states via thermal
operations when $\beta=0$. However, as $\beta\geq 0$, the ground state
$\boldsymbol{p}_{1}$ loses its ability to produce a Bell state, while three
other pure states remain capable. The maximum negativity that
$\boldsymbol{p}_{1}$ can produce decreases as $\beta$ increases.

##  VI Additional discussions

###  VI.1 Generation of entanglement in the presence of dissipation

So far, our discussion revolved around the assumption that any energy-
preserving system-bath dynamics is accessible. However, one can also consider
the scenario where the system undergoes an open dynamics described by a
Lindblad master equation [43, 44, 45],

$\frac{d\rho(t)}{dt}=-i\quantity[H,\rho(t)]+\mathcal{L}_{t}\quantity[\rho(t)],$
(34)

where $[\cdot,\cdot]$ denotes the commutator and $\mathcal{L}_{t}$ is the
Lindbladian having the general form

$\mathcal{L}_{t}(\rho)=\sum_{i}r_{i}(t)\left[L_{i}(t)\rho
L_{i}(t)^{\dagger}-\frac{1}{2}\Bigl{\\{}L_{i}(t)^{\dagger}L_{i}(t),\rho\Bigl{\\}}\,\\!\right]\\!,\\!$
(35)

with $\\{\cdot,\cdot\\}$ denoting the anticommutator, $L_{i}(t)$ being jump
operators, and $r_{i}(t)\geq 0$ being non-negative jump rates. Furthermore,
assuming that the quantum system interacts with a large heat bath implies that
the system Gibbs state is a stationary solution of the dynamics,
$\mathcal{L}_{t}(\gamma)=0$ and that the Lindbladian $\mathcal{L}_{t}$
commutes with the generator of the Hamiltonian dynamics $-i[H,\cdot]$ for all
times $t$.

Figure 8: Entanglement generation via Markovian thermal process. The upper
plot illustrates the volume of the future thermal cone of entanglement as a
function of $\beta$ for a bipartite system
$\rho_{\textsf{A}B}=\outerproduct{11}{11}$ undergoing a Markovian thermal
process, depicted by the solid green curve. The dashed orange line represents
the volume of the future thermal cone of entanglement for the same state, but
under thermal operations. In the lower plot, the evolution of the Markovian
future thermal cone of entanglement (green) with respect to $\beta$ is shown.
Here, the orange shape depicts the future thermal cone under TO, while the
internal lighter color with green vertices represents the future cone under
MTP.

The dynamics generated by master equations that satisfy two aforementioned
properties lead to what is known as a Markovian thermal process, which is
known to be a strict subset of thermal operations [46]. The differences
between them arise from our degree of control on the bath. Nonetheless, the
entanglement criterion given in Eq. (26) remains valid for Markovian thermal
processes. Consequently, we can study the entanglement generation in the
presence of dissipation using the same criterion function $f$.

The additional restriction stemming from the Markovianity of the Lindblad
master equations is the most evident when considering the highest excited
state $\boldsymbol{p}=(0,0,0,1)$ as the initial state. As analysed in Section
V-C, the volume of the future thermal cone of entanglement
$\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{p})$ defined in Eq. (20)
remains maximal for all inverse temperatures $\beta$, when we can utilize any
thermal operation to generating entanglement. In Fig. 8a, we present how the
volume $\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{p})$ varies as a
function of $\beta$, when only Markovian thermal processes are allowed. For
any $\beta$, the volume
$\mathbf{V}^{\mathrm{ENT}}_{\mathrm{TO}}(\boldsymbol{p})>0$ indicates that the
entanglement generation via dissipative dynamics is possible. However, the
difference between the entanglement generation via thermal operations is the
most pronounced in the high-temperature limit, i.e. when $\beta$ is small.
Furthermore, in the limit $\beta\rightarrow\infty$, Markovian thermal
processes are as capable as thermal operations for generating entanglement.
The detailed behaviour of the future thermal cone of entanglement under MTP is
depicted in Fig. 8b.

###  VI.2 Remarks on higher-dimensional systems

The entanglement structure for multipartite systems with more than two parties
is notoriously complex. Multipartite entanglement can be classified using
equivalence classes defined with invertible local transformations. Two states
are considered to have the same class of entanglement if one can be
transformed into the other via LOCC with a non-zero probability [47]. Even in
the simplest three-qubit systems, one encounters two distinct classes of
entanglement: GHZ-like entanglement and W-like entanglement [47, 48],
represented by the states

$\displaystyle\ket{\mathrm{GHZ}}=\frac{\ket{000}+\ket{111}}{\sqrt{2}},\ \
\ket{\mathrm{W}}=\frac{\ket{001}+\ket{010}+\ket{100}}{\sqrt{3}}.$ (36)

The $n$-qubit extensions of two classes exist. In particular, $n$-qubit GHZ
states are defined as

$\ket{\mathrm{GHZ}_{n}}=\frac{1}{\sqrt{2}}\left(\underbrace{\ket{0\ldots
0}}_{n}+\underbrace{\ket{1\ldots 1}}_{n}\right),$ (37)

an equal superpositions of all $n$ qubits being in the ground state and in the
excited state. The concept of W-like entanglement is similarly generalized to
$n$-qubit Dicke states

$\ket{D_{m}}=\binom{n}{m}^{-1/2}\sum_{\sigma\in\mathcal{S}_{n}}P_{\sigma}\quantity(\underbrace{\ket{0\ldots
0}}_{n-m}\otimes\underbrace{\ket{1\ldots 1}}_{m})\ ,$ (38)

where the sum goes over all permutations $\sigma$ on $n$ elements and
$P_{\sigma}$ are the respective matricial representations. These states form a
basis for the symmetric subspace of an $n$-qubit system.

Remarkably, different energy level structures give rise to different
entanglement structures that can be achieved using thermal operations. Let us
assume that the Hamiltonian of the $i$-th qubit is
$H_{i}=E_{i}\left|1\right\rangle\\!\\!\left\langle 1\right|$. First, if all
the energy gaps are equal, i.e. $E_{i}=E$, then the $n$-qubit Hilbert space is
divided into energy-degenerate subspaces:

$\mathcal{H}=\bigoplus_{i=0}^{n}\mathcal{V}_{iE},$ (39)

each of which spanned by states with a fixed number of excitations $i$ and
thus the energy equal to $i\cdot E$. With this energy-level structure, Dicke
states $\ket{D_{i}}$, or W-like entanglement in the three-qubit case, can be
generated via thermal operations with nonzero components within the energy-
degenerate subspaces.

On the other hand, if we set $E_{n}=\sum_{i=1}^{n-1}E_{i}$, we find that the
degenerate subspace corresponding to this energy

$\mathcal{V}_{E_{n}}=\mathrm{span}\quantity(|\underbrace{1\ldots
1}_{n-1}0\rangle,\,|\underbrace{0\ldots 0}_{n-1}1\rangle)$ (40)

contains states possessing GHZ-like entanglement. Specifically, in the case of
three qubits, this configuration leads to GHZ states up to a NOT operation on
the third qubit. Similar arguments can be easily extended to systems with
larger local dimensions.

Concerning systems with larger local dimensionality, we note that PPT
criterion, based on negativity of partial transpose of a state, is a complete
entanglement measure not only for two-qubit systems, but also for qubit-qutrit
systems as well [40]. In particular, we may consider a system in an energy-
incoherent state $\rho$ described by population
$\boldsymbol{p}=\quantity(p_{1},\ldots,p_{6})$ and subject to local
hamiltonians $H_{1}=\outerproduct{1}{1}$ and
$H_{2}=\outerproduct{1}{1}+2\outerproduct{2}{2}$. Based on similar arguments
as for the qubit-qubit systems, if any of the following two conditions,

$\displaystyle 4\,p_{1}\min(p_{4},p_{5})$ $\displaystyle<(p_{2}-p_{3})^{2},$
(41) $\displaystyle 4p_{6}\min(p_{2},p_{3})$
$\displaystyle<(p_{4}-p_{5})^{2},$

is satisfied, one of the eigenvalues of the partial transpose $\rho^{T_{A}}$
can become negative under energy-preserving unitary operations, and therefore,
the state is entanglable. Similar expressions can be easily put forward for
qubit-qudit systems.

It should be highlighted, however, that in general, complete entanglement
measures for multipartite systems beyond $2\times 2$ and $2\times 3$ do not
exist, unlike the one based on the PPT criterion, which is effective for two-
qubit and two-qutrit cases [39, 40], which is reflected in the appearance of
so-called bound entanglement [49]. Nevertheless, there exist partial
entanglement measures, such as the 3-tangle for three-qubit systems [50], can
be used to provide an outer approximation of the entanglable and non-
entanglable sets in such scenarios.

##  VII Outlook

In this paper, we have investigated the interplay between entanglement and
thermodynamics. The minimal thermodynamic assumptions we made, namely that the
system-bath is closed and evolves via an energy-preserving unitary, constrain
the final entanglement achievable from given initial states. By applying the
PPT criterion and constructing future thermal cones, we derived necessary and
sufficient conditions for two-qubit initial states that can be transformed
into entangled states within thermal operations framework. It turns out that
any separable state, apart from a cluster of states around the thermal state,
can become entangled. In other words, our investigation demonstrates that
entanglement generation is possible, given certain degree of athermality in
the initial state. Furthermore, we captured how the entanglement generating
capability varies with the ambient temperature, by analysing the volume of
entanglable initial state set, and found that low temperature allows more
states to become entanglable. Nevertheless, when the initial state is fixed,
we discovered that lowering temperature is not always beneficial and that more
detailed considerations on the geometry of thermal cones are needed.

There are several promising avenues for extending and generalizing the results
presented in this paper. A concrete and direct generalization consists of
using an auxiliary system, a catalyst, to facilitate a process that would not
occur spontaneously; see recent reviews [51, 52]. The idea of using a catalyst
in our context arises from Ref. [53], where the authors demonstrate that
catalysis enables the creation of otherwise inaccessible types of correlations
by accessing those locked in non-degenerate energy subspaces. As we have
discussed, the entanglement generation capacity depends on the structure of
degenerate energy subspaces. Employing a catalyst not only increases the
number of degenerate energy subspaces but also lift some of the dynamical
restrictions imposed by energy conservation. In this regard, the initial steps
demonstrating the utility of generating entanglement via thermal operations
with the aid of a strict catalyst were presented in Ref. [54]. Furthermore, in
Appendix E, we provide an example where correlated catalysis enables the
generation of entanglement from an initial state that cannot be entangled
without a catalyst. This raises the intriguing question of how one can
systematically study the generation of entanglement with the assistance of a
catalyst.

In this work, we only considered entanglement generation from initially
energy-incoherent states. When we focus on the entanglement generated between
degenerate energy levels, coherence between different energy levels cannot
contribute. Nevertheless, coherence can still nontrivially affect the
entanglement, including those between energy levels with different energy
values. Currently, tools and methods available for analysing coherent state
transformations under thermal operations are rather limited [36, 55, 56]. For
two-qubit systems, existing techniques to characterize single qubit coherent
state transformations [55, 36] may be adapted to study the entanglement
generation. This direction of study will extend our results to the realm where
three resources – entanglement, athermality, and coherence – compete and
cooperate with each other.

We make a few final remarks on existing literature, to which our work can
connect. First, it would be interesting to study the conventional quantifiers
of entanglement such as entanglement cost and distillable entanglement [1]
which are originally defined via LOCC, and evaluate them with regard to
thermodynamic constraints instead. In particular, the gap between these
quantities would then quantify the amount of irreversibility when it comes to
the local preparation of entanglement. Lastly, entanglement in continuous-
variable (CV) systems is relatively less explored. Extending our analysis to
CV systems would be extremely fruitful, considering their widespread use in
practical experiments [57] involving optomechanics [58, 59], and
superconducting setups [60]. For instance, we may use the toolkits developed
in the study of Gaussian thermal operations [61] for an analysis similar to
what have been done in the current work, and aim to explicate CV entanglement
results, such as the recent no-go result [62] for the CV steady-state
entanglement generation.

###### Acknowledgements.

We thank Kamil Korzekwa for valuable discussions that significantly
contributed to the execution of this project, and also Patryk Lipka-Bartosik
for very fruitful comments on the first version of the manuscript. AOJ
acknowledges financial support from VILLUM FONDEN through a research grant
(40864). JS and NN are supported by the start-up grant for Nanyang Assistant
Professorship of Nanyang Technological University, Singapore, awarded to Nelly
H. Y. Ng. JCz acknowledges financial support by NCN PRELUDIUM BIS no.
DEC-2019/35/O/ST2/01049.

## References

  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Kosloff [2013] R. Kosloff, Quantum thermodynamics: A dynamical viewpoint, Entropy 15, 2100 (2013).
  * Goold  _et al._ [2016] J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, The role of quantum information in thermodynamics—a topical review, J. Phys. A Math. Theor. 49, 143001 (2016).
  * Vinjanampathy and Anders [2016] S. Vinjanampathy and J. Anders, Quantum thermodynamics, Contemporary Physics 57, 545 (2016).
  * Deffner and Campbell [2019] S. Deffner and S. Campbell,  _Quantum Thermodynamics: An introduction to the thermodynamics of quantum information_ (Morgan & Claypool Publishers, 2019).
  * Oppenheim  _et al._ [2002] J. Oppenheim, M. Horodecki, P. Horodecki, and R. Horodecki, Thermodynamical approach to quantifying quantum correlations, Phys. Rev. Lett. 89, 180402 (2002).
  * Hovhannisyan  _et al._ [2013] K. V. Hovhannisyan, M. Perarnau-Llobet, M. Huber, and A. Acín, Entanglement generation is not necessary for optimal work extraction, Phys. Rev. Lett. 111, 240401 (2013).
  * Huber  _et al._ [2015] M. Huber, M. Perarnau-Llobet, K. V. Hovhannisyan, P. Skrzypczyk, C. Klöckl, N. Brunner, and A. Acín, Thermodynamic cost of creating correlations, New J. Phys 17, 065008 (2015).
  * Perarnau-Llobet  _et al._ [2015] M. Perarnau-Llobet, K. V. Hovhannisyan, M. Huber, P. Skrzypczyk, N. Brunner, and A. Acín, Extractable work from correlations, Phys. Rev. X 5, 041011 (2015).
  * Alicki and Fannes [2013] R. Alicki and M. Fannes, Entanglement boost for extractable work from ensembles of quantum batteries, Phys. Rev. E 87, 042123 (2013).
  * Verstraete  _et al._ [2001] F. Verstraete, K. Audenaert, and B. De Moor, Maximally entangled mixed states of two qubits, Phys. Rev. A 64, 012316 (2001).
  * Friis  _et al._ [2016] N. Friis, M. Huber, and M. Perarnau-Llobet, Energetics of correlations in interacting systems, Phys. Rev. E 93, 042135 (2016).
  * Ishizaka and Hiroshima [2000] S. Ishizaka and T. Hiroshima, Maximally entangled mixed states under nonlocal unitary operations in two qubits, Phys. Rev. A 62, 022310 (2000).
  * Plenio and Huelga [2002] M. B. Plenio and S. F. Huelga, Entangled light from white noise, Phys. Rev. Lett. 88, 197901 (2002).
  * Hartmann  _et al._ [2006] L. Hartmann, W. Dür, and H.-J. Briegel, Steady-state entanglement in open and noisy quantum systems, Phys. Rev. A 74, 052304 (2006).
  * Hartmann  _et al._ [2007] L. Hartmann, W. Dür, and H. Briegel, Entanglement and its dynamics in open, dissipative systems, New J. Phys 9, 230 (2007).
  * Bellomo and Antezza [2013] B. Bellomo and M. Antezza, Steady entanglement out of thermal equilibrium, Europhysics Letters 104, 10006 (2013).
  * Brask  _et al._ [2015] J. B. Brask, G. Haack, N. Brunner, and M. Huber, Autonomous quantum thermal machine for generating steady-state entanglement, New J. Phys 17, 113029 (2015).
  * Khandelwal  _et al._ [2024] S. Khandelwal, B. Annby-Andersson, G. F. Diotallevi, A. Wacker, and A. Tavakoli, Maximal steady-state entanglement in autonomous quantum thermal machines, arXiv:2401.01776 (2024).
  * Gour  _et al._ [2015] G. Gour, M. P. Müller, V. Narasimhachar, R. W. Spekkens, and N. Yunger Halpern, The resource theory of informational nonequilibrium in thermodynamics, Phys. Rep 583, 1–58 (2015).
  * Brandão  _et al._ [2013] F. G. S. L. Brandão, M. Horodecki, J. Oppenheim, J. M. Renes, and R. W. Spekkens, Resource theory of quantum states out of thermal equilibrium, Phys. Rev. Lett. 111, 250404 (2013).
  * Janzing  _et al._ [2000] D. Janzing, P. Wocjan, R. Zeier, R. Geiss, and T. Beth, Thermodynamic cost of reliability and low temperatures: Tightening landauer’s principle and the second law, Int. J. Theor. Phys. 39, 2717 (2000).
  * Horodecki and Oppenheim [2013] M. Horodecki and J. Oppenheim, (quantumness in the context of) resource theories, Int. J. Mod. Phys. B 27, 1345019 (2013).
  * Brandão  _et al._ [2015] F. G. S. L. Brandão, M. Horodecki, N. H. Y. Ng, J. Oppenheim, and S. Wehner, The second laws of quantum thermodynamics, Proc. Natl. Acad. Sci. U.S.A. 112, 3275 (2015).
  * Lostaglio [2019a] M. Lostaglio, An introductory review of the resource theory approach to thermodynamics, Rep. Prog. Phys. 82, 114001 (2019a).
  * Sparaciari  _et al._ [2021] C. Sparaciari, M. Goihl, P. Boes, J. Eisert, and N. H. Y. Ng, Bounding the resources for thermalizing many-body localized systems, Commun. Phys. 4, 3 (2021).
  * Czartowski  _et al._ [2023] J. Czartowski, A. de Oliveira Junior, and K. Korzekwa, Thermal recall: Memory-assisted markovian thermal processes, PRX Quantum 4, 040304 (2023).
  * Son and Ng [2023] J. Son and N. H. Y. Ng, A hierarchy of thermal processes collapses under catalysis, arXiv:2303.13020 (2023).
  * Alhambra  _et al._ [2019] Á. M. Alhambra, M. Lostaglio, and C. Perry, Heat-bath algorithmic cooling with optimal thermalization strategies, Quantum 3, 188 (2019).
  * Lostaglio [2019b] M. Lostaglio, An introductory review of the resource theory approach to thermodynamics, Rep. Prog. Phys 82, 114001 (2019b).
  * Torun  _et al._ [2023] G. Torun, O. Pusuluk, and Ö. E. Müstecaplıoğlu, A compendious review of majorization-based resource theories: quantum information and quantum thermodynamics, Turk. J. Phys 47, 141 (2023).
  * Horodecki and Oppenheim [2013] M. Horodecki and J. Oppenheim, Fundamental limitations for quantum and nanoscale thermodynamics, Nat. Commun. 4, 2059 (2013).
  * Jaynes and Cummings [1963] E. T. Jaynes and F. W. Cummings, Comparison of quantum and semiclassical radiation theories with application to the beam maser, Proc. IEEE 51, 89 (1963).
  * Cirac and Zoller [1994] J. I. Cirac and P. Zoller, Preparation of macroscopic superpositions in many-atom systems, Phys. Rev. A 50, R2799 (1994).
  * Vidal and Werner [2002] G. Vidal and R. F. Werner, Computable measure of entanglement, Phys. Rev. A 65, 032314 (2002).
  * Korzekwa [2017] K. Korzekwa, Structure of the thermodynamic arrow of time in classical and quantum theories, Phys. Rev. A 95, 052318 (2017).
  * de Oliveira Junior  _et al._ [2022] A. de Oliveira Junior, J. Czartowski, K. Życzkowski, and K. Korzekwa, Geometric structure of thermal cones, Phys. Rev. E 106, 064109 (2022).
  * Lostaglio  _et al._ [2018] M. Lostaglio, Á. M. Alhambra, and C. Perry, Elementary Thermal Operations, Quantum 2, 52 (2018).
  * Peres [1996] A. Peres, Separability criterion for density matrices, Phys. Rev. Lett. 77, 1413 (1996).
  * Horodecki  _et al._ [1996] M. Horodecki, P. Horodecki, and R. Horodecki, Separability of mixed states: necessary and sufficient conditions, Phys. Lett. A 223, 1 (1996).
  * Horodecki  _et al._ [1998] M. Horodecki, P. Horodecki, and R. Horodecki, Mixed-state entanglement and distillation: Is there a “bound” entanglement in nature?, Phys. Rev. Lett. 80, 5239 (1998).
  * Chen  _et al._ [2017] L. Chen, Y. Yang, and W.-S. Tang, Schmidt number of bipartite and multipartite states under local projections, Quantum Inf Process 16, 1 (2017).
  * Kossakowski [1972] A. Kossakowski, On quantum statistical mechanics of non-hamiltonian systems, Rep. Math. Phys. 3, 247 (1972).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of n-level systems, J. Math. Phys. 17, 821 (1976).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * Lostaglio and Korzekwa [2022] M. Lostaglio and K. Korzekwa, Continuous thermomajorization and a complete set of laws for Markovian thermal processes, Phys. Rev. A 106, 012426 (2022).
  * Dür  _et al._ [2000] W. Dür, G. Vidal, and J. I. Cirac, Three qubits can be entangled in two inequivalent ways, Phys. Rev. A 62, 062314 (2000).
  * Acín  _et al._ [2001] A. Acín, D. Bruß, M. Lewenstein, and A. Sanpera, Classification of mixed three-qubit states, Phys. Rev. Lett. 87, 040401 (2001).
  * Gaida and Kleinmann [2023] M. Gaida and M. Kleinmann, Seven definitions of bipartite bound entanglement, Journal of Physics A: Mathematical and Theoretical 56, 385302 (2023).
  * Coffman  _et al._ [2000] V. Coffman, J. Kundu, and W. K. Wootters, Distributed entanglement, Phys. Rev. A 61, 052306 (2000).
  * Datta  _et al._ [2023] C. Datta, T. V. Kondra, M. Miller, and A. Streltsov, Catalysis of entanglement and other quantum resources, Rep. Prog. Phys 86, 116002 (2023).
  * Lipka-Bartosik  _et al._ [2023a] P. Lipka-Bartosik, H. Wilming, and N. H. Y. Ng, Catalysis in quantum information theory, arXiv:2306.0079 (2023a).
  * Lipka-Bartosik  _et al._ [2023b] P. Lipka-Bartosik, G. F. Diotallevi, and P. Bakhshinezhad, Fundamental limits on anomalous energy flows in correlated quantum systems, arXiv:2307.03828 (2023b).
  * Czartowski and de Oliveira Junior [2024] J. Czartowski and A. de Oliveira Junior, We can catalyse it for you thermally: Catalytic transformations for thermal operations, arXiv:2403.XXXXX (2024).
  * Lostaglio  _et al._ [2015] M. Lostaglio, K. Korzekwa, D. Jennings, and T. Rudolph, Quantum coherence, time-translation symmetry, and thermodynamics, Phys. Rev. X 5, 021001 (2015).
  * Gour  _et al._ [2018] G. Gour, D. Jennings, F. Buscemi, R. Duan, and I. Marvian, Quantum majorization and a complete set of entropic conditions for quantum thermodynamics, Nat. Commun. 9 (2018).
  * Andersen  _et al._ [2010] U. Andersen, G. Leuchs, and C. Silberhorn, Continuous-variable quantum information processing, Laser Photonics Rev 4, 337 (2010).
  * Rocheleau  _et al._ [2009] T. Rocheleau, T. Ndukum, C. Macklin, J. B. Hertzberg, A. A. Clerk, and K. C. Schwab, Preparation and detection of a mechanical resonator near the ground state of motion, Nature 463, 72–75 (2009).
  * Teufel  _et al._ [2011] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W. Simmonds, Sideband cooling of micromechanical motion to the quantum ground state, Nature 475, 359–363 (2011).
  * Reiter  _et al._ [2013] F. Reiter, L. Tornberg, G. Johansson, and A. S. Sørensen, Steady-state entanglement of two superconducting qubits engineered by dissipation, Phys. Rev. A 88, 032317 (2013).
  * Serafini  _et al._ [2020] A. Serafini, M. Lostaglio, S. Longden, U. Shackerley-Bennett, C.-Y. Hsieh, and G. Adesso, Gaussian thermal operations and the limits of algorithmic cooling, Phys. Rev. Lett. 124, 010602 (2020).
  * Longstaff  _et al._ [2023] B. Longstaff, M. G. Jabbour, and J. B. Brask, Impossibility of bosonic autonomous entanglement engines in the weak-coupling limit, Phys. Rev. A 108, 032209 (2023).
  * Mazurek and Horodecki [2018] P. Mazurek and M. Horodecki, Decomposability and convex structure of thermal processes, New J. Phys. 20, 053040 (2018).
  * Mazurek [2019] P. Mazurek, Thermal processes and state achievability, Phys. Rev. A 99, 042110 (2019).
  * Son and Ng [2024] J. Son and N. H. Y. Ng, Catalysis in action via elementary thermal operations, New Journal of Physics (2024).

##  Appendix A Future thermal cone

Figure 9: Thermomajorisation $\&$ future thermal cone. For two states
$\boldsymbol{p}=(0.7,0.2,0.1)$ and $\boldsymbol{q}=(0.6,0.2,0.2)$ with
$\beta=0.5$, (a) we plot their thermomajorization curves
$\mathcal{L}_{\boldsymbol{p}}(x)$ and $\mathcal{L}_{\boldsymbol{q}}(x)$,
respectively. Since $\boldsymbol{p}$ thermomajorizes $\boldsymbol{q}$, it
implies that one can transform $\boldsymbol{p}$ into $\boldsymbol{q}$ under
thermal operations. This can also be observed by examining their (b) future
thermal cones, where $\boldsymbol{q}$ lies within the future thermal cone of
$\boldsymbol{p}$. The Gibbs state $\boldsymbol{\gamma}$ is depicted by a black
star $\star$.

We briefly review state transition rules under thermal operations and the
notion of the future thermal cone, a set of state that can be transformed from
a given initial state with thermal operations. First, we introduce
thermomajorization relations, which fully characterize the convertibility
between two states via thermal operations. The thermomajorization curves,
associated to each energy-incoherent state, are defined as the following:

######  Definition 5 (Thermomajorization curve).

Suppose that the Gibbs state $\boldsymbol{\gamma}\in\Delta_{d}$ with respect
to the system Hamiltonian $H$ and the ambient inverse temperature $\beta$ is
given. For each state $\boldsymbol{p}\in\Delta_{d}$, we define the
thermomajorisation curve of $\boldsymbol{p}$ as a piecewise-linear function
$\mathcal{L}_{\boldsymbol{p}}:[0,1]\to[0,1]$ composed of linear segments
connecting the point $(0,0)$ and the elbow points
$\\{(\sum_{i=1}^{k}\gamma_{\pi_{\boldsymbol{p}}^{-1}(i)},~{}\sum_{i=1}^{k}p_{\pi_{\boldsymbol{p}}^{-1}(i)})\\}_{k=1}^{d}$.
Here $\pi_{\boldsymbol{p}}$ is the $\beta$-ordering of $\boldsymbol{p}$, which
imposes
$p_{\pi_{\boldsymbol{p}}^{-1}(i)}/\gamma_{\pi_{\boldsymbol{p}}^{-1}(i)}\geq
p_{\pi_{\boldsymbol{p}}^{-1}(j)}/\gamma_{\pi_{\boldsymbol{p}}^{-1}(j)}$ for
all $i\leq j$.

Note that thermomajorization curves are always concave and non-decreasing by
definition. For the Gibbs state $\boldsymbol{\gamma}$, the thermomajorization
curve $\mathcal{L}_{\boldsymbol{\gamma}}$ is a straight line from $(0,0)$ to
$(1,1)$. From the concavity of the curve and the fact that any
$\mathcal{L}_{\boldsymbol{p}}$ passes points $(0,0)$ and $(1,1)$, any
thermomajorization curve is above the curve for the Gibbs state. i.e.
$\mathcal{L}_{\boldsymbol{p}}(x)\geq\mathcal{L}_{\boldsymbol{\gamma}}(x)$ for
all $x\in[0,1]$. See Fig. 9 (a) for examples of thermomajorization curves.

Using thermomajorization curves, thermomajorization relations between two
states can be established.

######  Definition 6 (Thermomajorization relation).

The thermomajorization relation is a preorder between two energy-incoherent
states. If two states $\boldsymbol{p},\boldsymbol{q}\in\Delta_{d}$ have
corresponding thermomajorization curves $\mathcal{L}_{\boldsymbol{p}}$ and
$\mathcal{L}_{\boldsymbol{q}}$, such that
$\mathcal{L}_{\boldsymbol{p}}(x)\geq\mathcal{L}_{\boldsymbol{q}}(x)$ for all
$x\in[0,1]$, we say $\boldsymbol{p}$ thermomajorizes $\boldsymbol{q}$, or more
simply, $\boldsymbol{p}\succ_{\beta}\boldsymbol{q}$.

Since the thermomajorization relation gives a preorder, not a total order,
there exist pairs of states $\boldsymbol{p}$ and $\boldsymbol{q}$ such that
neither $\boldsymbol{p}\succ_{\beta}\boldsymbol{q}$ nor
$\boldsymbol{q}\succ_{\beta}\boldsymbol{p}$ holds. Such pairs are dubbed
incomparable. Reamrkably, thermomajorization relations coincide with the state
convertibility under thermal operations.

######  Theorem 4 (Theorem 2 of Ref. [32]).

Suppose that $\boldsymbol{p},\boldsymbol{q}\in\Delta_{d}$ describe two energy-
incoherent states of a system with Hamiltonian $H$ with access to thermal
environments with the inverse temperature $\beta$. If
$\boldsymbol{p}\succ_{\beta}\boldsymbol{q}$, than $\boldsymbol{p}$ can be
transformed to $\boldsymbol{q}$ via some thermal operation.

In particular, the Gibbs state $\boldsymbol{\gamma}$ is thermomajorized by any
other state and thus any state can be converted to the Gibbs state.
Conversely, the most excited state with $\boldsymbol{p}=(0,\cdots,0,1)$
thermomajorizes all other states in $\Delta_{d}$ and any state can be
converted from it.

Theorem 4 provides an efficient way of checking whether a transition from one
state to another is feasible. However, it is also possible to characterize the
entire set of states that can be attained from an initial state via thermal
operations. We refer to such set with an initial state $\boldsymbol{p}$ as the
future thermal cone of $\boldsymbol{p}$ and denote it with a symbol
$\mathcal{T}_{+}(\boldsymbol{p})$.

######  Lemma 5 (Lemma 12 of Ref. [38]).

Suppose $\boldsymbol{p}\in\Delta_{d}$. The future thermal cone
$\mathcal{T}_{+}(\boldsymbol{p})$ is a convex polytope with at most $d!$
extreme points each corresponding to permutations $\boldsymbol{\pi}\in{\cal
S}_{d}$. The extreme point $\boldsymbol{p}^{\boldsymbol{\pi}}$ corresponding
to the permutation $\boldsymbol{\pi}$ can be identified by first constructing
its following thermomajorization curve
$\mathcal{L}_{\boldsymbol{p}^{\boldsymbol{\pi}}}$:

  1. 1.

The $x$-coordinates of the $\mathcal{L}_{\boldsymbol{p}^{\boldsymbol{\pi}}}$
elbow points are $x_{k}=\sum_{i=1}^{k}\gamma_{\pi^{-1}(i)}$ for
$k\in\left\\{1,\dots,d\right\\}$.

  2. 2.

The $y$-coordinates of the $\mathcal{L}_{\boldsymbol{p}^{\boldsymbol{\pi}}}$
elbow points are $y_{k}=\mathcal{L}_{\boldsymbol{p}}(x_{k})$ for
$k\in\left\\{1,\dots,d\right\\}$.

Therefore, $\mathcal{L}_{\boldsymbol{p}^{\boldsymbol{\pi}}}$ is the piecewise-
linear curve that connects $(0,0)$ and elbow points
$\\{(x_{k},y_{k})\\}_{k=1}^{d}$, which all lie on the curve
$\mathcal{L}_{\boldsymbol{p}}$.

The thermomajorization relation between the initial state $\boldsymbol{p}$ and
extreme states $\boldsymbol{p}^{\boldsymbol{\pi}}$ are also known as tight-
thermomajorization, since all elbow points of
$\mathcal{L}_{\boldsymbol{p}^{\boldsymbol{\pi}}}$ are on the curbe
$\mathcal{L}_{\boldsymbol{p}}$. Furthermore, by construction, each state
$\boldsymbol{p}^{\boldsymbol{\pi}}$ has the $\beta$-ordering
$\boldsymbol{\pi}$. Lemma 5 enables us to characterize the future thermal cone
of $\boldsymbol{p}$ by constructing states
$\boldsymbol{p}^{\,\boldsymbol{\pi}}$ for each
$\boldsymbol{\pi}\in\mathcal{S}_{d}$, and by taking their convex hull. Here,
we provide the construction of ${\cal T}_{+}(\boldsymbol{p})$ as a simple
corollary of Lemma 5:

######  Theorem 6 (Future thermal cone).

The future thermal cone of a $d$-dimensional energy-incoherent state
$\boldsymbol{p}$ is given by

${\cal
T}_{+}(\boldsymbol{p})=\operatorname{conv}[\\{\boldsymbol{p}^{\boldsymbol{\pi}},\boldsymbol{\pi}\in{\cal
S}_{d}\\}]\,.$ (42)

See Fig. 9 (b) for example future thermal cones of two three-level states. In
this example, $\boldsymbol{p}\succ_{\beta}\boldsymbol{q}$ and ${\cal
T}_{+}(\boldsymbol{q})\subset{\cal T}_{+}(\boldsymbol{p})$; these two
relations are indeed equivalent in general. The equivalence directly follows
from the transitivity of the preorder $\succ_{\beta}$, i.e. if
$\boldsymbol{p}\succ_{\beta}\boldsymbol{q}$ and
$\boldsymbol{q}\succ_{\beta}\boldsymbol{r}$, it is always true that
$\boldsymbol{p}\succ_{\beta}\boldsymbol{r}$; hence, $\boldsymbol{r}\in{\cal
T}_{+}(\boldsymbol{q})$ implies $\boldsymbol{r}\in{\cal
T}_{+}(\boldsymbol{p})$.

It is worth mentioning that $\mathcal{T}^{\,\beta}_{+}(\boldsymbol{p})$ can
also be constructed by identifying the entire set of extremal Gibbs-preserving
stochastic matrices, as discussed in [38, 63, 64]. However, if we are only
interested in the future thermal cones, following Lemma 5 and Theorem 6 is
more efficient.

##  Appendix B Convexity of the subspace non-entanglable set $\mathbbm{NE}$

Figure 10: Subspace non-entanglable set $\mathbbm{NE}$. Geometry of the set of
qubit-qubit states that cannot get entangled using energy-preserving
unitaries. The extreme points of this set include the trivial ones, given by
$(1,0,0,0),(0,0,0,1)$ and $(0,\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}},0)$, as
well as a continuous family of boundary points specified by the solution to
the quadratic equation $-(p_{2}-p_{3})^{2}+4p_{1}(1-p_{1}-p_{2}-p_{3})=0$ when
solved for $p_{3}$.

For a bipartite two-qubit state, we now show that the subspace non-entanglable
set $\mathbbm{NE}$ is convex. The set $\mathbbm{NE}$ is defined (in Definition
2) as the set of states satisfying $f(\boldsymbol{p})\geq 0$, where $f$ is a
function given in Eq. (27):

$f(\boldsymbol{p})=4p_{1}(1-p_{1}-p_{2}-p_{3})-(p_{2}-p_{3})^{2}.$ (43)

To show that the set of states defined by $f(\boldsymbol{p})\geq 0$ is convex,
we need to demonstrate that the function $f$ is, counterintuitively, concave.
Our task is to verify that for any two arbitrary points
$P_{1}=(x_{1},y_{1},z_{1})$ and $P_{2}=(x_{2},y_{2},z_{2})$, the following
condition holds:

$f[\lambda(x_{1},y_{1},z_{1})+(1-\lambda)(x_{2},y_{2},z_{3})]\geq\lambda
f(x_{1},y_{1},z_{1})+(1-\lambda)f(x_{2},y_{2},z_{2}),$ (44)

for $\lambda\in[0,1]$. To make the calculation less cumbersome, we can prove
convexity by considering two arbitrary boundary points. To do this, we first
find the roots of Eq. (43) with respect to $p_{3}$:

$f(p_{1},p_{2},p_{3})=-p^{2}_{3}+2p_{3}(p_{2}-2p_{1})-[p_{2}(4p_{1}+p_{2})-4p_{1}(1-p_{1})]=0,$
(45)

whose solutions are simply given by

$\displaystyle p^{(\pm)}_{3}=-2p_{1}+p_{2}\pm 2\sqrt{p_{1}-2p_{1}p_{2}},$ (46)

By focusing on the boundary conditions, we can simplify the analysis
concerning convexity. Specifically, when considering two arbitrary boundary
points, the right-hand side of Eq. (44) is zero. This reduces the problem to
verifying whether the left-hand side of this equation is less than or equal to
zero. Given that there are two solutions for Eq. (43), we need to demonstrate
that all possible combinations of boundary points
$P_{i}=(x_{1},y_{1},p^{(\pm)}_{3})$ for $i\\{1,2\\}$, the following
inequalities are satisfied:

$\displaystyle
f[\lambda(x_{1},y_{1},p^{(\pm)}_{3})+(1-\lambda)(x_{2},y_{2},p^{(\pm)}_{3})]$
$\displaystyle\geq 0,$ (47) $\displaystyle
f[\lambda(x_{1},y_{1},p^{(\pm)}_{3})+(1-\lambda)(x_{2},y_{2},p^{(\mp)}_{3})]$
$\displaystyle\geq 0,$ (48)

Furthermore, we must impose the constrain $y_{i}\leq 1/2$ for $i\in\\{1,2\\}$
to ensure that the values within the square root in (46) are real. By
rearranging and rewriting the inequalities above, we observe that both
inequalities in (47), as well as inequalities (48), are equivalent.
Consequently, our analysis boils down to demonstrating the validity of the
following two inequalities:

$\displaystyle 4\lambda(1-\lambda)[x_{1}(1-2y_{2})+x_{2}(1-2y_{1})\pm
2\sqrt{(x_{1}-2x_{1}y_{1})(x_{2}-2x_{2}y_{2})}]$ $\displaystyle\geq 0.$ (49)

Since the function $4\lambda(1-\lambda)$ is concave and positive on
$\lambda\in[0,1]$ interval, we only need to prove that the second term is
positive. By setting $a=\sqrt{x_{1}(1-2y_{2})}$ and $b=\sqrt{x_{2}(1-2y_{1}}$
we see that the entire expression is reduced to

$\displaystyle 4\lambda(1-\lambda)(a\pm b)^{2}$ $\displaystyle\geq 0,$ (50)

which completes the proof of concavity of function $f$.

##  Appendix C Proof of Theorem 2

We first establish the existence of a subspace entanglable $\beta$-ordering
$(2,1,3,4)$ state for all thermally entanglable states. Yet, the subspace
entanglability of the extreme point $\boldsymbol{p}^{\star}$ is not guaranteed
from Lemma 7 yet.

######  Lemma 7 (Critical $\beta$-ordering).

For any state $\boldsymbol{p}\in\mathbbm{TE}$, there exists a state
$\boldsymbol{q}$, such that
$\boldsymbol{q}\in\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{p})$ and the
$\beta$-order $\boldsymbol{\pi}_{\boldsymbol{q}}=(2,1,3,4)$.

###### Proof.

Suppose that the Lemma is true for all $\boldsymbol{p}\in\mathbbm{E}$. Then
for any state $\boldsymbol{r}\in\mathbbm{TE}$, there exists
$\boldsymbol{p}\in\mathbbm{E}$ and
$\boldsymbol{q}\in\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{p})\subset\mathcal{T}_{+}^{\mathrm{ENT}}(\boldsymbol{r})$
with $\boldsymbol{\pi}_{\boldsymbol{q}}=(2,1,3,4)$. Hence, it is sufficient to
prove the Lemma for $\boldsymbol{p}\in\mathbbm{E}$.

Without loss of generality, assume that $p_{2}\geq p_{3}$. From the assumption
$\boldsymbol{p}\in\mathbbm{E}$, we have $f(\boldsymbol{p})<0$ for $f$ defined
in Eq. (27). Now, consider another state $\boldsymbol{r}$ such that
$\boldsymbol{r}\in\mathcal{T}_{+}(\boldsymbol{p})$ and

$\displaystyle\boldsymbol{r}=\left(p_{1}-\delta_{1},p_{2}+\sum_{i=1,3,4}\delta_{i},p_{3}-\delta_{3},p_{4}-\delta_{4}\right),$
(51)

for some $\delta_{1},\delta_{3},\delta_{4}\geq 0$. We can observe that
$f(\boldsymbol{r})$ is decreasing monotonically as
$\delta_{1},\delta_{3},\delta_{4}$ increases. If
$(\boldsymbol{\pi}_{\boldsymbol{p}})_{1}\neq 2$, there always exists a state
$\boldsymbol{r}$ in the form of Eq. (51) and
$(\boldsymbol{\pi}_{\boldsymbol{r}})_{1}=2$. From the monotonicity of
$f(\boldsymbol{r})$, this state
$\boldsymbol{r}\in\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{p})$. We prove
the rest of the Lemma by showing there exists
$\boldsymbol{q}\in\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{r})$ and use
$\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{r})\subset\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{p})$.
Hence, we assume that $(\boldsymbol{\pi}_{\boldsymbol{p}})_{1}=2$.

Now, suppose that the state $\boldsymbol{p}$ has $p_{4}\geq p_{1}e^{-2\beta
E}$. Then a state $\boldsymbol{r}$ such that $r_{2}=p_{2}$, $r_{3}=p_{3}$, and

$\displaystyle r_{1}=\left(1-e^{-2\beta E}\right)p_{1}+p_{4},\quad
r_{4}=e^{-2\beta E}p_{1},$ (52)

is in $\mathcal{T}_{+}(\boldsymbol{p})$ and satisfies $r_{1}e^{-2\beta E}\geq
r_{4}$. From direct calculation,

$\displaystyle f(\boldsymbol{r})=f(\boldsymbol{p})-4\left(1-e^{-2\beta
E}\right)p_{1}\left(p_{4}-e^{-2\beta E}p_{1}\right)\leq f(\boldsymbol{p}),$
(53)

i.e. $\boldsymbol{r}\in\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{p})$
whenever $\boldsymbol{p}\in\mathbbm{E}$.

Therefore, one can always find
$\boldsymbol{q}\in\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{p})$ whose
$\beta$-ordering is either _(i)_ $\pi_{\boldsymbol{q}}=(2,3,1,4)$, _(ii)_
$\pi_{\boldsymbol{q}}=(2,1,4,3)$, or _(iii)_ $\pi_{\boldsymbol{q}}=(2,1,3,4)$.
We will show that cases _(i)_ and _(ii)_ again imply the existence of
$\boldsymbol{q}\in\mathcal{T}^{\mathrm{ENT}}_{+}(\boldsymbol{p})$
corresponding to the case _(iii)_ , which proves the Lemma.

Case _(i)_ : $\pi_{\boldsymbol{p}}=(2,3,1,4)$ with a strict ordering
$p_{3}>p_{1}e^{-\beta E}$. It is always possible to find
$\boldsymbol{p}^{\prime}\in\mathcal{T}_{+}(\boldsymbol{p})$, such that
$\pi_{\boldsymbol{p}^{\prime}}=(2,3,1,4)$ and $f(\boldsymbol{p}^{\prime})=0$
by a convex combination
$\boldsymbol{p}^{\prime}=(1-\lambda)\boldsymbol{p}+\lambda\boldsymbol{\gamma}$,
where $\lambda\in[0,1]$ and $\gamma$ is the Gibbs state. Now assume that
$\boldsymbol{q}\in\mathcal{T}_{+}(\boldsymbol{p}^{\prime})$ is a state with
$\pi_{\boldsymbol{q}}=(2,1,3,4)$ and

$\displaystyle\boldsymbol{q}=\left(p^{\prime}_{1}+\delta,p^{\prime}_{2},p^{\prime}_{3}-\delta,p^{\prime}_{4}\right),$
(54)

with some $\delta>0$. Such a state always exists since
$p^{\prime}_{3}>p^{\prime}_{1}e^{-\beta E}$. The entanglement witness gives

$\displaystyle f(\boldsymbol{q})$
$\displaystyle=4\left(p^{\prime}_{1}+\delta\right)\left(p^{\prime}_{4}\right)-\left(p^{\prime}_{2}-p^{\prime}_{3}+\delta\right)^{2}$
$\displaystyle=-\delta^{2}+2\left(p^{\prime}_{3}+2p^{\prime}_{4}-p^{\prime}_{2}\right)\delta,$
(55)

where we used $f(\boldsymbol{p}^{\prime})=0$ for the second equality. If
$p^{\prime}_{1}=0$, we have $p^{\prime}_{4}=0$ from the $\beta$-order of
$\boldsymbol{p}^{\prime}$ and $p^{\prime}_{2}=p^{\prime}_{3}$ from
$f(\boldsymbol{p}^{\prime})=0$; hence, $f(\boldsymbol{q})<0$. When
$p^{\prime}_{1}\neq 0$, we use
$p^{\prime}_{4}=(p^{\prime}_{2}-p^{\prime}_{3})^{2}/(4p^{\prime}_{1})$ to get

$\displaystyle
f(\boldsymbol{q})=-\delta^{2}-2\left(p^{\prime}_{2}-p^{\prime}_{3}\right)\left(1-\frac{p^{\prime}_{2}-p^{\prime}_{3}}{2p^{\prime}_{1}}\right)\delta.$
(56)

Recall $4(p^{\prime}_{1})^{2}\geq
4p^{\prime}_{1}p^{\prime}_{4}=(p^{\prime}_{2}-p^{\prime}_{3})^{2}$ implying
$1-\frac{p^{\prime}_{2}-p^{\prime}_{3}}{2p^{\prime}_{1}}\geq 0$ and thus
$f(\boldsymbol{q})<0$.

Case _(ii)_ : $\pi_{\boldsymbol{p}}=(2,1,4,3)$. There exists
$\boldsymbol{p}^{\prime}\in\mathcal{T}_{+}(\boldsymbol{p})$, such that
$\pi_{\boldsymbol{p}^{\prime}}=(2,1,4,3)$ and $f(\boldsymbol{p}^{\prime})=0$
and also $\boldsymbol{q}\in\mathcal{T}_{+}(\boldsymbol{p}^{\prime})$, such
that

$\displaystyle\boldsymbol{q}=\left(p^{\prime}_{1},p^{\prime}_{2},p^{\prime}_{3}+\delta,p^{\prime}_{4}-\delta\right),$
(57)

with $\pi_{\boldsymbol{q}}=(2,1,3,4)$. Observe that

$\displaystyle f(\boldsymbol{q})$
$\displaystyle=4\left(p^{\prime}_{1}\right)\left(p^{\prime}_{4}-\delta\right)-\left(p^{\prime}_{2}-p^{\prime}_{3}-\delta\right)^{2}$
$\displaystyle=-\delta^{2}-2\left(2p^{\prime}_{1}-p^{\prime}_{2}+p^{\prime}_{3}\right)\delta<0$
(58)

for any $\delta>0$. ∎

Now we use geometric arguments to prove the remaining part of Theorem 2. To do
this, we define generalized boundary sets.

######  Definition 7 ($d$-dimensional boundaries).

Let $X$ be a convex set and $C_{d}(X)$ be a set of all $(d+1)$-dimensional
convex subsets $Y\subset X$. $d$-dimensional boundary of $X$ is defined as

$\displaystyle{\operatorname{Bd}}_{d}(X)\coloneqq\left\\{x\,\big{|}\,\not\exists
Y\in C_{d}(X),\ \text{such that}\ x\in\mathrm{int}(Y)\right\\}.$ (59)

Intuitively, ${\operatorname{Bd}}_{d}(X)$ cannot be written as a convex
combination of $d+2$ extreme points of $X$ with nonzero weights; hence,
${\operatorname{Bd}}_{0}(X)$ is a set of extreme points of $X$, while
${\operatorname{Bd}}_{1}(X)$ is a set of points of the edges of $X$.
Furthermore, if $X$ is $D$-dimensional, ${\operatorname{Bd}}_{D-1}(X)$ reduces
to the conventional boundary, $\partial X$

Note that ${\operatorname{Bd}}_{d}(X)\subset{\operatorname{Bd}}_{d+1}(X)$ for
any $d\geq 0$: if it is impossible write a state $x$ as a convex combination
of any $d+2$ extreme points, it is also impossible to write it as a convex
combination of any $d+3$ extreme points.

Using this definition, we establish a helpful proposition.

######  Proposition 8.

Suppose that a convex set $A$ is a subset of a convex set $B$ and their
boundaries interesect, i.e. $A\subset B$ and $\partial A\cap\partial
B\neq\emptyset$. Then, for all $d$-dimensional boundaries,

$\displaystyle
A\cap{\operatorname{Bd}}_{d}(B)\subset{\operatorname{Bd}}_{d}(A).$ (60)

If we have additional condition
$[{\operatorname{Bd}}_{d}(A)\setminus{\operatorname{Bd}}_{d-1}(A)]\cap{\operatorname{Bd}}_{d}(B)=\emptyset$,
for some $d$, we find

$\displaystyle
A\cap{\operatorname{Bd}}_{d}(B)\subset{\operatorname{Bd}}_{d-1}(A).$ (61)

###### Proof.

For the first part of the proposition, suppose that there exists a point
$x\notin{\operatorname{Bd}}_{d}(A)$ such that $x\in
A\cap{\operatorname{Bd}}_{d}(B)$. Then there exists a convex subset $Y$ of $A$
having $d+2$ extreme points and containing $x$ in $\mathrm{int}(Y)$. Yet, this
set $Y$ is also a convex subset of $B$ having the same property, which
contradicts $x\in A\cap{\operatorname{Bd}}_{d}(B)$.

The second part of the proposition follows directly from combining Eq. (60)
and the additional assumption
$[{\operatorname{Bd}}_{d}(A)\setminus{\operatorname{Bd}}_{d-1}(A)]\cap{\operatorname{Bd}}_{d}(B)=\emptyset$.
∎

Now we proceed to prove Theorem 2, which states that
$\boldsymbol{p}^{\star}\in\mathbbm{E}$ whenever
$\boldsymbol{p}\in\mathbbm{TE}$. First, consider mixtures of the initial state
$\boldsymbol{p}$ and the Gibbs state $\boldsymbol{\gamma}$ that can be written
as $\boldsymbol{p}^{(x)}=(1-x)\boldsymbol{p}+x\boldsymbol{\gamma}$ for some
$x\in[0,1]$. Future thermal cones between these states have a hierarchy
$\mathcal{T}_{+}(\boldsymbol{p}^{(y)})\subset\mathcal{T}_{+}(\boldsymbol{p}^{(x)})$
whenever $y\geq x$. Then there exists a coefficient $a\in(0,1)$, such that

$\displaystyle\mathcal{T}_{+}\left(\boldsymbol{p}^{(a)}\right)\subset\mathbbm{NE},\quad\partial\left(\mathcal{T}_{+}\left(\boldsymbol{p}^{(a)}\right)\right)\cap\partial\left(\mathbbm{NE}\right)\neq\emptyset,$
(62)

since $\mathcal{T}_{+}(\boldsymbol{p})\not\subset\mathbbm{NE}$, while
$\mathcal{T}_{+}(\boldsymbol{\gamma})\subset\mathrm{int}(\mathbbm{NE})$. The
Gibbs state $\boldsymbol{\gamma}$ is always full-rank and so is
$\boldsymbol{p}^{(a)}$, i.e. $\boldsymbol{p}^{(a)}\in\mathrm{int}(\Delta_{4})$
and $\mathcal{T}_{+}(\boldsymbol{p})\subset\mathrm{int}(\Delta_{4})$.

Now we show that
$[{\operatorname{Bd}}_{1}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)}))\setminus{\operatorname{Bd}}_{0}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)}))]\cap\partial(\mathbbm{NE})=\emptyset$.
Let $\boldsymbol{q}$ be a state on the edge of
$\mathcal{T}_{+}(\boldsymbol{p}^{(a)})$ but not on the vertex of it, i.e.
$\boldsymbol{q}\in[{\operatorname{Bd}}_{1}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)}))\setminus{\operatorname{Bd}}_{0}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)}))]$.
Since $\mathbbm{NE}$ is convex, if $\boldsymbol{q}\in\partial(\mathbbm{NE})$,
the entire edge containing $\boldsymbol{q}$ should be included within
$\partial(\mathbbm{NE})$. Recall that $\mathcal{T}_{+}(\boldsymbol{p}^{(a)})$
is a polyhedron whose edges are defined as segments of lines on which two of
the populations are fixed, e.g.
$\boldsymbol{q}=(p^{(a)}_{1},p^{(a)}_{2},t,1-p^{(a)}_{1}-p^{(a)}_{2}-t)$ and
for some $t\in(t_{1},t_{2})$. Then the entire edge can be written as a set of
states

$\displaystyle\left\\{\left(p^{(a)}_{1},p^{(a)}_{2},t,1-p^{(a)}_{1}-p^{(a)}_{2}-t\right)\,\big{|}\,t\in[t_{1},t_{2}]\right\\}\subset\partial(\mathbbm{NE}).$
(63)

Yet, it is impossible to satisfy
$f(p^{(a)}_{1},p^{(a)}_{2},t,1-p^{(a)}_{1}-p^{(a)}_{2}-t)=0$ for a finite
range $t\in[t_{1},t_{2}]$, which contradicts Eq. (63). The same can be proven
for any two fixed populations. Likewise,
$[{\operatorname{Bd}}_{2}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)}))\setminus{\operatorname{Bd}}_{1}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)}))]\cap\partial(\mathbbm{NE})=\emptyset$:
if any point on
$[{\operatorname{Bd}}_{2}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)}))\setminus{\operatorname{Bd}}_{1}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)}))]$
is included in $\partial(\mathbbm{NE})$, at least one entire edge of
$\mathcal{T}_{+}(\boldsymbol{p}^{(a)})$ should also be included, which is not
possible.

Using Proposition 8 twice with $A=\mathcal{T}_{+}(\boldsymbol{p}^{(a)})$ and
$B=\mathbbm{NE}$, we have

$\displaystyle\mathcal{T}_{+}(\boldsymbol{p}^{(a)})\cap\partial(\mathbbm{NE})=\mathcal{T}_{+}(\boldsymbol{p}^{(a)})\cap\partial(\mathbbm{NE})\subset{\operatorname{Bd}}_{0}(\mathcal{T}_{+}(\boldsymbol{p}^{(a)})),$
(64)

i.e. only extreme points of $\mathcal{T}_{+}(\boldsymbol{p}^{(a)})$ can be
included in the set $\partial(\mathbbm{NE})$.

We finally use Lemma 7. Since $\boldsymbol{p}^{(a)}\in\partial(\mathbbm{TE})$,
there must exist $\boldsymbol{q}\in\mathcal{T}_{+}(\boldsymbol{p}^{(a)})$ with
$\beta$-ordering $(2,1,3,4)$ that is in
$\partial(\mathbbm{E})=\partial(\mathbbm{NE})$. However, only extreme points
of $\mathcal{T}_{+}(\boldsymbol{p}^{(a)})$ can be in $\partial(\mathbbm{NE})$,
implying that the $(2,1,3,4)$ extreme point

$\displaystyle(1-a)\boldsymbol{p}^{\star}+a\boldsymbol{\gamma}\in\partial(\mathbbm{NE}).$
(65)

From $\boldsymbol{\gamma}\in\mathrm{int}(\mathbbm{NE})$, we obtain
$\boldsymbol{p}^{\star}\in\mathbbm{E}$, which concludes the proof.

##  Appendix D Thermal initial states

In this section, we derive the critical temperatures above or below which the
initially thermal state become entanglable. We denote the inverse temperature
of the system as $\beta_{\textsf{S}}$ and that of the environment as $\beta$.
To simplify the notation, we use Boltzmann weights
$\Delta_{\textsf{S}}\coloneqq e^{-\beta_{\textsf{S}}E}$, $\Delta\coloneqq
e^{-\beta E}$, , and partition functions
$Z_{\textsf{S}}\coloneqq(1+\Delta_{\textsf{S}})^{2}$,
$Z\coloneqq(1+\Delta)^{2}$. Then the initial state can be written as

$\displaystyle\boldsymbol{p}=\frac{1}{Z_{\textsf{S}}}\left(1,\Delta_{\textsf{S}},\Delta_{\textsf{S}},\Delta_{\textsf{S}}^{2}\right).$
(66)

When $\beta_{S}>\beta$, i.e. when the system is cooler than the environment,
the $\beta$-ordering of $\boldsymbol{p}$ is $(1,2,3,4)$. The extreme point
$\boldsymbol{p}^{\star}$ of $\mathcal{T}_{+}(\boldsymbol{p})$, corresponding
to the ordering $\boldsymbol{\pi}^{\star}=(2,1,3,4)$ can be easily obtained by
maximally exchanging populations (referred to as the $\beta$-swapping in Refs.
[38, 27, 65]) between first two levels $\ket{00}$ and $\ket{01}$,

$\displaystyle\boldsymbol{p}^{\star}=\frac{1}{Z_{\textsf{S}}}\left(1-\Delta+\Delta_{\textsf{S}},\Delta,\Delta_{\textsf{S}},\Delta_{\textsf{S}}^{2}\right).$
(67)

We want the entanglement measure to be negative, i.e.

$\displaystyle
f(\boldsymbol{p}^{\star})=\frac{1}{Z_{\textsf{S}}^{2}}\left(4\Delta_{\textsf{S}}^{2}\left[1-(\Delta-\Delta_{\textsf{S}})\right]-(\Delta-\Delta_{\textsf{S}})^{2}\right)<0.$
(68)

Solving this inequality, we obtain

$\displaystyle\Delta>\Delta_{C_{1}}\coloneqq\Delta_{\textsf{S}}\quantity(1+2\sqrt{1+\Delta_{\textsf{S}}^{2}}-2\Delta_{\textsf{S}}),$
(69)

equivalent to Eq. (29) in Section V-B.

When the system temperature is low, that is, when $\Delta_{\textsf{S}}\ll 1$,
Eq. (69) can be written with the Taylor expansion

$\displaystyle\Delta_{C_{1}}=3\Delta_{\textsf{S}}+O(\Delta_{\textsf{S}}^{2}),$
(70)

yielding the approximation $\beta_{C_{1}}E\simeq\beta_{\textsf{S}}E-\log 3$,
up to the leading order.

Now, consider the other case where $\beta_{\textsf{S}}<\beta$; the system is
now hotter than the surrounding environment, and the $\beta$-order of the
initial state $\mathbf{p}$ is now $(4,3,2,1)$. Thermomajorization relations
reveal that the extreme point $\boldsymbol{p}^{\star}$ now has two distinct
forms depending on the environment temperature $\beta$. If
$1>\Delta+\Delta^{2}$, that is, $\beta E\geq-\log[\frac{\sqrt{5}-1}{2}]$, the
extreme state is

$\displaystyle\boldsymbol{p}^{\star}=\frac{1}{Z_{\textsf{S}}}\left(1+(\Delta_{\textsf{S}}-\Delta)(1+\Delta),\Delta_{\textsf{S}}(1+\Delta_{\textsf{S}}-\Delta),\Delta,\Delta^{2}\right).$
(71)

On the other hand, when $\beta E<-\log[\frac{\sqrt{5}-1}{2}]$,

$\displaystyle\boldsymbol{p}^{\star}=\frac{1}{Z_{\textsf{S}}}\left(\frac{\Delta_{\textsf{S}}}{\Delta},\Delta_{\textsf{S}}(1+\Delta_{\textsf{S}}-\Delta),1+\Delta_{\textsf{S}}+\Delta_{\textsf{S}}\Delta-\Delta-\frac{\Delta_{\textsf{S}}}{\Delta},\Delta^{2}\right).$
(72)

For the latter case Eq. (72), the entanglement measure is evaluated as

$\displaystyle
f(\boldsymbol{p}^{\star})=\frac{1}{Z_{\textsf{S}}^{2}\Delta^{2}}\left(4\Delta_{\textsf{S}}\Delta^{3}-(\Delta_{\textsf{S}}-\Delta)^{2}(1+\Delta(\Delta_{\textsf{S}}-\Delta))^{2}\right).$
(73)

However, in the range where Eq. (72) is valid,
$1\geq\Delta_{\textsf{S}}>\Delta>\frac{\sqrt{5}-1}{2}$, entanglement measure
is always positive and there is no critical inverse temperature
$\beta_{C_{2}}$. For the former case Eq. (71) with
$1\geq\Delta_{\textsf{S}}>\Delta$ and $\frac{\sqrt{5}-1}{2}\geq\Delta>0$,

$\displaystyle
f(\boldsymbol{p}^{\star})=\frac{1}{Z_{\textsf{S}}^{2}}\left[-(\Delta_{\textsf{S}}-\Delta)^{2}(1+\Delta_{\textsf{S}})^{2}+4\Delta^{2}\left(1+\Delta_{\textsf{S}}-(1-\Delta_{\textsf{S}})\Delta+\Delta^{2}\right)\right]$
(74)

and the solution for $f(\boldsymbol{p}^{\star})=0$ can be obtained by solving
this quartic equation. Unfortunately, there is no simple explicit solution for
this equation. To approximate the critical value, let us assume that
$\Delta_{\textsf{S}}\ll 1$, i.e. $\beta_{\textsf{S}}E\gg 1$. Ignoring the
$Z_{\textsf{S}}$ factor in Eq. (74) and taking only the leading order terms of
$\Delta$ and $\Delta_{\textsf{S}}$, we get the approximation

$\displaystyle
Z_{\textsf{S}}^{2}f(\boldsymbol{p}^{\star})=3\Delta^{2}+2\Delta\Delta_{\textsf{S}}-\Delta_{\textsf{S}}^{2}+O(\Delta_{\textsf{S}}^{3}),$
(75)

and the solution of $f(\boldsymbol{p}^{\star})=0$ is approximately
$\Delta\simeq\frac{\Delta_{\textsf{S}}}{3}$. In other words,
$\beta_{C_{2}}E\simeq\beta_{\textsf{S}}E+\log 3$.

##  Appendix E Catalytic advantages on entanglement generation - a toy example

We now construct a simple toy example to demonstrate that if a catalyst is
allowed, entanglement can be generated from states that cannot achieve this
under thermal operations without a catalyst. For simplicity, we assume
$\beta=0$ and consider a bipartite state

$\rho_{\textsf{AB}}=\frac{2}{5}\outerproduct{00}{00}+\frac{1}{4}\outerproduct{01}{01}+\frac{33}{100}\outerproduct{10}{10}+\frac{1}{50}\outerproduct{11}{11},$
(76)

with local Hamiltonians $H_{\textsf{X}}=\outerproduct{1}{1}$ for
$\textsf{X}\in\\{\textsf{A},\textsf{B}\\}$. The initial bipartite state
$\rho_{\textsf{AB}}$ belongs to the thermally non-entanglable set, i.e.,
$\boldsymbol{p}=\mathrm{eig}(\rho)\in\mathbbm{TNE}$. Hence, at $\beta=0$,
there is no thermal operation capable of entangling $\rho_{\textsf{AB}}$.

Let us now consider a two-dimensional catalyst, prepared in a state

$\omega_{\textsf{C}}=\frac{73}{100}\outerproduct{0}{0}+\frac{27}{100}\outerproduct{1}{1},$
(77)

and described described by a Hamiltonian $H_{\textsf{C}}=\outerproduct{1}{1}$,
such that the joint Hamiltonian of the composite system ABC is given by

$H_{\textsf{ABC}}=\pNiceMatrix[margin]\Block[fill=boxcolor!,rounded-
corners]{1-1}{}0&0000000\\\ 0\Block[fill=block2c!,rounded-
corners]{3-3}{}1000000\\\ 00100000\\\ 00010000\\\
0000\Block[fill=block3c!,rounded-corners]{3-3}{}2000\\\ 00000200\\\
00000020\\\ 0000000\Block[fill=block4c!,rounded-corners]{1-1}{}3\\\ .$ (78)

where the colors in each block highlight the degenerate energy subspaces. We
assume that the total system evolves according to the following energy-
preserving unitary transformation

$U=\pNiceMatrix[margin]\Block[fill=boxcolor!,rounded-
corners]{1-1}{}1&0000000\\\ 0\Block[fill=block2c!,rounded-
corners]{3-3}{}1000000\\\ 00010000\\\ 00100000\\\
0000\Block[fill=block3c!,rounded-corners]{3-3}{}0100\\\ 00000010\\\
00001000\\\ 0000000\Block[fill=block4c!,rounded-corners]{1-1}{}1\\\ ,$ (79)

which maps the composite system to
$\sigma_{\textsf{ABC}}:=U(\rho_{\textsf{AB}}\otimes\omega_{\textsf{C}})U^{\dagger}$,
where the reduced states of AB become

$\displaystyle\sigma_{\textsf{AB}}=\frac{949}{2000}\outerproduct{00}{00}+\frac{613}{5000}\outerproduct{01}{01}+\frac{771}{2500}\outerproduct{10}{10}+\frac{189}{2000}\outerproduct{11}{11}.$
(80)

while the catalyst returns exactly to the same state as it started
$\sigma_{\textsf{C}}=\omega_{\textsf{C}}$.

Therefore, one can easily verify that the final state no longer belongs to
$\mathbbm{TNE}$, i.e.,
$\boldsymbol{q}:=\operatorname{eig}(\sigma)\notin\mathbbm{TNE}$. This implies
that a thermal operation capable of mapping $\sigma_{\textsf{AB}}$ into an
entangled state can be constructed.