# Extracting work from correlated many-body quantum systems

Raffaele Salvia<EMAIL_ADDRESS>Scuola Normale Superiore and
University of Pisa, I-56127 Pisa, Italy  Vittorio Giovannetti  NEST, Scuola
Normale Superiore and Istituto Nanoscienze-CNR, I-56127 Pisa, Italy

###### Abstract

The presence of correlations in the input state of a non-interacting many-body
quantum system can lead to an increase in the amount of work we can extract
from it under global unitary processes (ergotropy). The present work explore
such effect on translationally invariant systems relaying on the Matrix
Product Operator formalism to define a measure of how much they are
correlated. We observe that in the thermodynamic limit of large number of
sites, complete work extraction can be attained for relatively small
correlation strength (a reduction of a 2 factor in dB unit). Most importantly
such an effect appears not to be associated with the presence of quantum
correlations (e.g. entanglement) in the input state (classical correlation
suffices), and to be attainable by only using incoherent ergotropy. As a
byproduct of our analysis we also present a rigorous formulation of the
heuristic typicality argument first formulated in [Alicki and Fannes, 2013],
which gives the maximum work extractable for a set of many identical quantum
systems in the asymptotic limit.

##  I Introduction

For purely classical models, the maximum work $\Delta W^{(\max)}_{\rm{clas}}$
that can be extracted from a closed system by means of reversible, iso-
entropic cycles is equal to the difference between the input energy $E_{in}$
and the energy $E_{th}$ of the thermal equilibrium state that has the same
entropy of the initial configuration, i.e. $\Delta
W^{(\max)}_{\rm{clas}}=E_{in}-E_{th}$ Thomson (1852); Fermi  _et al._ (1955).
The corresponding quantity for a quantum system with Hamiltonian $H$,
initialised in a state $\rho$, is called  _ergotropy_ $\mathcal{E}(\rho;H)$
Allahverdyan  _et al._ (2004) and is formally defined as the maximal amount of
energy one can recover from the system by means of reversible coherent (i.e.
unitary) processes. In general $\mathcal{E}(\rho;H)$ is strictly smaller than
the threshold one would get from a naive translation of the classical optimal
term $\Delta W^{(\max)}_{\rm{clas}}$, i.e. the difference
$E(\rho;H)-\mathfrak{E}_{S(\rho)}$ between the mean energy
$E(\rho;H):=\mbox{Tr}[H\rho]$ of the input state $\rho$, and the mean energy
$\mathfrak{E}_{S(\rho)}:=\mbox{Tr}[H\omega_{\beta}]$ of the thermal Gibbs
state $\omega_{\beta}:=e^{-\beta H}/Z_{\beta}$ that has the same von Neumann
entropy $S(\rho)$. Nonetheless, thanks to the fact that $\mathcal{E}(\rho;H)$
is an extensive super-additive functional, the gap between the work
extractable at the quantum level and $E(\rho;H)-\mathfrak{E}_{S(\rho)}$ can be
progressively reduced by allowing joint operations on an increasing collection
of identical copies of the system. In particular, indicating with $H^{(N)}$
the global Hamiltonian of the $N$ copies written in terms of a sum of
independent, uniform, local terms (see Eq. (2) below), it turns out that the
ergotropy per-site ${\mathcal{E}(\rho^{\otimes N};H^{(N)})}/{N}$ of the model
is a non-decreasing function of $N$ whose asymptotic value (a quantity called
from time-to-time total ergotropy Niedenzu  _et al._ (2019)) verifies the
identity Alicki and Fannes (2013)

$\displaystyle\mathcal{E}_{\text{tot}}(\rho;H)$ $\displaystyle:=$
$\displaystyle\lim_{N\to\infty}\frac{\mathcal{E}(\rho^{\otimes
N};H^{(N)})}{N}$ (1) $\displaystyle=$ $\displaystyle
E(\rho;H)-\mathfrak{E}_{S(\rho)}\;.$

The limit (1) corresponds to the saturation of the classical threshold for
reversible, iso-entropic work extraction processes involving a quantum system.
In all non trivial cases it is explicitly attained for an infinite number of
copies, the only exceptions being associated to the cases where
$\mathcal{E}(\rho^{\otimes N};H^{(N)})=0$ for all $N$. A direct consequence of
this observation is that when operating on non correlated input quantum system
$\rho^{\otimes N}$ we can extract less work than in the classical case. Such
limitation however does not necessary hold if we allowed quantum correlations
in the input state of the system: indeed if we permit operations that act on
the entire many-body system, the classical threshold (1) can be overcame
allowing for the extraction of additional work Huber  _et al._ (2015);
Francica  _et al._ (2017) (Notice that this is the opposite of what happens if
instead we are forced to act only on individual parts of the quantum state: in
this scenario in fact correlations typically are detrimental for work
extraction see e.g. Refs. Alicki and Fannes (2013); Oppenheim  _et al._
(2002); Vitagliano  _et al._ (2018); Goold  _et al._ (2016); Bera  _et al._
(2017); Manzano  _et al._ (2018); Andolina  _et al._ (2019)). Aim of the
present article is to investigate this issue by studying the ergotropy
$\mathcal{E}(\rho^{(N)};H^{(N)})$ of correlated multipartite quantum states
$\rho^{(N)}$ of translationally invariant many-body quantum systems composed
by $N$ sites. To express how much the state $\rho^{(N)}$ is correlated, we
will use as a measure the minimum bond link rank $M$ necessary to represent
$\rho^{(N)}$ as a Matrix Product Operator Verstraete  _et al._ (2004):
indicating with $r$ the rank (number of strictly positive eigenvalues) of the
single-site density matrix $\rho$ of the model, the quantity $M$ can range
from $1$ (corresponding to factorized state $\rho^{(N)}=\rho^{\otimes N}$), to
$r^{2}$ (corresponding to the case of where $\rho^{(N)}$ is a pure GHZ
entangled state). Via constructive examples, we shall hence produce lower
bounds for the maximum ergotropy per site $\mathcal{E}(\rho^{(N)};H^{(N)})/N$
attainable in the system for fixed values of the correlation parameter $M$,
comparing them to the classical threshold limit (1) reachable when the
correlations are removed from the model.

Our main finding is to show that, while forcing $\rho^{(N)}$ to be pure by
taking $M=r^{2}$ allows one to trivially boost
$\mathcal{E}(\rho^{(N)};H^{(N)})/N$ to its natural upper bound (i.e. the
single-site mean energy $E(\rho;H)$), in the limit of large $N$ one can
asymptotically reach the same result for much smaller values of $M$. In
particular we prove that setting $M=r$ is enough, corresponding to a reduction
of a factor 2 in dB unit. Most importantly it turns out that while for
$M=r^{2}$ the saturation to $E(\rho;H)$ relays explicitly on the presence of
the quantum correlations (entanglement) that forces $\rho^{(N)}$ to be pure,
for $M=r$ the same can be attained by only exploiting classical correlations.
Our analysis shows also that this effect can be attained through operations
which do no affect the coherence of the input state Baumgratz  _et al._
(2014): using the notation introduced in Francica  _et al._ (2020); Çakmak
(2020); Touil  _et al._ (2021), this means that the asymptotic saturation of
$\mathcal{E}(\rho^{(N)};H^{(N)})/N$ to $E(\rho;H)$ we report here for $M=r$
can be obtained by just using  _incoherent ergotropy_.

We conclude mentioning that as a byproduct of our study we provide a rigorous
derivation of the heuristic typicality argument first formulated in Alicki and
Fannes (2013), which gives the maximum work extractable for a set of identical
quantum systems in the asymptotic limit of a large number of copies $N$.

The manuscript is organized as follows. In Sec. II we set the problem
introducing the model, presenting the notation, and giving the basic
definitions that will be used in the remaining of the paper. Also in an effort
to improve readability, in Sec. II.1 we present a brief technical summary of
the main results, with comments and general considerations. With Sec. III we
enter into the technical part of the manuscript: here in Sec. III.1 and III.2
we give a detailed account of the Matrix Product Operator representation of
translationally invariant states, while in Sec. III.3 and III.4 we focus on
special families of the states that, in the forthcoming sections will be
adopted to provide estimations for the attainable values of
$\mathcal{E}(\rho^{(N)};H^{(N)})/N$ for given bond link rank $M$. The proofs
of the results are presented in Sec. IV, and conclusions are drawn in Sec. V.
The paper contains also a series of technical appendices.

##  II The problem

In our analysis we shall focus on a many-body quantum system consisting on an
ordered collection $Q_{1}$, $Q_{2}$, $\cdots$, $Q_{N}$ of $N$ $d$-dimensional
sites – see Fig. 1. Indicating with ${\cal H}$ and ${\cal H}^{\otimes N}$ the
single-site and many-body Hilbert spaces, we represent with $\mathfrak{L}$ and
$\mathfrak{L}^{(N)}$ the associated spaces of the linear operators, and with
$\mathfrak{S}$ and $\mathfrak{S}^{(N)}$ the sets of the corresponding density
matrices. Assuming no interactions among the various sub-systems, we write the
joint Hamiltonian of the model as a sum of homogenous local terms

$\displaystyle H^{(N)}:=\sum_{\ell=1}^{N}H_{\ell}\;,$ (2)

with $H_{\ell}$ representing the same single-site operator $H\in\mathfrak{L}$
acting on the $\ell$-th subsystem $Q_{\ell}$. Without loss of generality in
what follows we shall take $H$ (hence also $H^{(N)}$) to be positive
semidefinite and to have zero ground state eigenvalue; in particular we shall
use the symbols $\epsilon_{j}$ to represent its eigenvalues that we order via
the inequalities,

$\displaystyle\epsilon_{d}\geq\epsilon_{d-1}\geq\cdots\geq\epsilon_{2}\geq\epsilon_{1}=0\;,$
(3)

and the symbols $|\epsilon_{i}\rangle$ to represent their associated
eigenvectors, i.e

$\displaystyle H|\epsilon_{j}\rangle=\epsilon_{j}|\epsilon_{j}\rangle.$ (4)

Given hence $\rho^{(N)}\in\mathfrak{S}^{(N)}$ a generic (possible correlated)
quantum state of the compound we write its ergotropy as Allahverdyan  _et al._
(2004)

${{\cal
E}(\rho^{(N)};H^{(N)})}:=E(\rho^{(N)};H^{(N)})-\min_{U}E(U\rho^{(N)}U^{\dagger};H^{(N)})\;,$
(5)

where the minimization is performed over all transformation $U\in{\bf U}({\cal
H}^{\otimes N})$ of the unitary set on ${\cal H}^{\otimes N}$, and where

$E(\rho^{(N)};H^{(N)}):=\Tr[\rho^{(N)}H^{(N)}]\;,$ (6)

represents the energy expectation values on $\rho^{(N)}$ which by construction
is non-negative.

Figure 1: Pictorial representation of the model: we consider a collection of
$N$ identical quantum systems $Q_{1}$, $Q_{2}$, $\cdots$, $Q_{N}$ initialized
into (possibly correlated) density matrices $\rho^{(N)}$ that are invariant
under translation of the subsystem indexes. The system Hamiltonian consists
into a sum of homogeneous local terms with no interactions. Fixing the single-
site density matrix $\rho$ of the model (see Eq. (8), the degree of intra-site
correlations of $\rho^{(N)}$ are measured in terms of the parameter $m$ which
defines the state BLR value (12): for $m=1$, $\rho^{(N)}$ corresponds to the
factorized state $\rho^{\otimes N}$, while for $m$ equal to the square of the
rank $r$ of $\rho$, we obtain the maximally entangled state
$|\mbox{GHZ}^{(N)}_{\rho}\rangle$ of Eq. (14).

In this paper we will be concerned on $\rho^{(N)}$ belonging to the special
subset $\mathfrak{S}_{T}^{(N)}$ of $\mathfrak{S}^{(N)}$ formed by density
matrices that are invariant under ciclic translations of the indexes sites,
i.e.

$\displaystyle\rho^{(N)}=T\rho^{(N)}T^{\dagger}\;,$ (7)

where $T$ is the transformation on ${\cal H}^{\otimes N}$ which for all
$\ell=1,\cdots,N$ maps $Q_{\ell}$ into $Q_{\ell\oplus 1}$ (here $\oplus$
represents the sum modulus $N$). Since by construction these states are
locally uniform, we can associate to each one of them the same single-site
reduced density matrix

$\Tr_{Q_{1}\dots Q_{\ell-1}Q_{\ell+1}\dots
Q_{N}}[\rho^{(N)}]=\rho\;,\quad\forall\ell\in[1,N]\;,$ (8)

that allows us to express the global mean energy of the system as

$E\left(\rho^{(N)};H^{(N)}\right)=NE(\rho;H)\;.$ (9)

Factorized density matrices $\rho^{(N)}=\rho^{\otimes N}$ are special
instances of $\mathfrak{S}_{T}^{(N)}$ which, due to the absence of
correlations among the various subsystems, are fully characterized by their
single-site counterpart (8). A proper representation of
$\mathfrak{S}_{T}^{(N)}$ can be obtained in terms of the _translationally
invariant matrix product operators_ (TI-MPO) formalism Verstraete  _et al._
(2004). The key observation here is that given $\\{\ket{i_{1}i_{2}\dots
i_{N}}=|i_{1}\rangle\otimes|i_{2}\rangle\otimes\cdots\otimes|i_{N}\rangle\\}$
an orthonormal basis of ${\cal H}^{\otimes N}$ constructed from the elements
of reference, single-site basis $\\{\ket{i};i=1,\cdots,d\\}$, for each
$\rho^{(N)}\in\mathfrak{S}_{T}^{(N)}$ it is possible to associate a four rank
tensor $\mathbb{A}\in\mathbb{C}^{d\times d\times M\times M}$ via the identity

$\displaystyle\rho^{(N)}$ $\displaystyle=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{i_{1}\dots
i_{N}}\bra{j_{1}\dots j_{N}}$ (10)
$\displaystyle\;\times\Tr[{A}^{i_{1},j_{1}}\dots A^{i_{N},j_{N}}]\;,$

with $\\{A^{i,j}|{i,j\in 1,\cdots d}\\}$ being $d^{2}$ matrices of dimension
$M\times M$ that we use to represent $\mathbb{A}$ – see Sec. III.1 for an
explicit derivation of this fact. In this construction the single-site
counterpart (8) of $\rho^{(N)}$ rewrites

$\rho=\sum_{i,j}\ket{i}\bra{j}\;\Tr[{A}^{i,j}{\bar{A}}^{N-1}]\;,\qquad\bar{A}:=\sum_{i}{A}^{i,i},$
(11)

while the value of the parameter $M$ can be used to gauge the strength of the
correlations among the various sites. For instance, as we shall see, setting
$M=d^{2}$ one can construct multi-site GHZ entangled pure states (see Sec.
III.3), while to get factorized states $\rho^{\otimes N}$ it is sufficient to
have $M=1$. Notice however the correspondence between $\mathfrak{S}_{T}^{(N)}$
and the set of tensors $\mathbb{A}$ is not one to one. In particular
determining whether the right-hand-side term of Eq. (10) would produce a
proper density matrix for a given $\mathbb{A}$ is a NP-hard problem, and it
becomes undecidable in the thermodynamic limit Kliesch  _et al._ (2014). Most
importantly for our purposes, multiple inequivalent choices of $\mathbb{A}$,
characterized by different values of $M$, can be assigned to each
$\rho^{(N)}\in\mathfrak{S}_{T}^{(N)}$. To remove this ambiguity we define the
_TI-MPO bond link rank_ (BLR) $\mathfrak{br}[\rho^{(N)}]$ as the minimum $M$
required to represent the state $\rho^{(N)}\in\mathfrak{S}_{T}^{(N)}$ in the
form (10) Navascues and Vertesi (2018):

$\mathfrak{br}[\rho^{(N)}]:=\min\left\\{M\geq
1\;\middle|\;\exists\;\mathbb{A}\in\mathbb{C}^{d\times d\times M\times
M}\right.\\\ \left.\mbox{s.t. Eq.~{}(\ref{MPO_rhoN}) holds true}\right\\}\;.$
(12)

This quantity is a proper functional of $\rho^{(N)}$: in particular, as
explicitly shown in Sec. III.2, it does not depend on the specific choice of
the local basis $\\{\ket{i};i=1,\cdots,d\\}$ that we use to define the TI-MPO
representation. Unfortunately determining the exact value of
$\mathfrak{br}[\rho^{(N)}]$ is a rather challenging task. For our analysis
however it will be sufficient to identify educated bounds for a special sub-
set of states.

Equipped with the above definitions we can now classify the elements of
$\mathfrak{S}_{T}^{(N)}$ in terms of their local properties and of their BLR.
Specifically, given $\rho\in\mathfrak{S}$ we first define
$\mathfrak{S}_{T}^{(N)}(\rho)$ as the subset of $\mathfrak{S}_{T}^{(N)}$
formed by states $\rho^{(N)}$ that admits $\rho$ as single-site density matrix
(8), i.e.

$\displaystyle\mathfrak{S}_{T}^{(N)}(\rho):=\Big{\\{}\rho^{(N)}\in\mathfrak{S}_{T}^{(N)}\;\big{|}\;$
$\displaystyle\qquad\Tr_{Q_{1}\dots Q_{\ell-1}Q_{\ell+1}\dots
Q_{N}}[\rho^{(N)}]=\rho\;,\;\forall\ell\in[1,N]\Big{\\}}\;,$

then we decompose $\mathfrak{S}_{T}^{(N)}(\rho)$ in terms of the associated
BLR, introducing the partitions

$\displaystyle\mathfrak{S}_{T}^{(N,m)}(\rho):=\Big{\\{}\rho^{(N)}\in\mathfrak{S}_{T}^{(N)}(\rho)\;\big{|}\;\mathfrak{br}[\rho^{(N)}]\leq
m\Big{\\}}\;.$

Notice that for $m=1$, $\mathfrak{S}_{T}^{(N,1)}(\rho)$ contains only the
factorized state $\rho^{\otimes N}$, i.e.

$\displaystyle\mathfrak{S}^{(N,1)}_{T}(\rho)=\Big{\\{}\rho^{\otimes
N}\Big{\\}}\;.$ (13)

On the contrary it can be shown that given $\mbox{r}[\rho]$ the rank of the
local density matrix $\rho$, the set
$\mathfrak{S}_{T}^{(N,\mbox{r}^{2}[\rho])}(\rho)$ includes the maximally
entangled pure state

$\displaystyle|\mbox{GHZ}^{(N)}_{\rho}\rangle:=\sum_{i=1}^{\mbox{r}[\rho]}\sqrt{\lambda_{i}}\;|\lambda_{i}\dots\lambda_{i}\rangle\;,$
(14)

with $\lambda_{i}>0$ being the non-zero eigenvalues of $\rho$ and
$|\lambda_{i}\rangle$ being the associated eigenvectors – an explicit proof of
this fact is presented in Sec. III.3. More generally for fixed $\rho$, the
$\mathfrak{S}_{T}^{(N,m)}(\rho)$’s form a family of sets of increasing size
which, in the non trivial case $\mbox{r}[\rho]\geq 2$, partition the
collection of density matrices $\rho^{(N)}$ fulfilling the local constraint
(8) into groups of increasing multi-site correlation strength, i.e.

$\mathfrak{S}_{T}^{(N,m-1)}(\rho)\subseteq\mathfrak{S}_{T}^{(N,m)}(\rho)\;,\quad\forall
m\geq 2\;,$ (15)

(for $\mbox{r}[\rho]=1$, i.e. when the local state is a pure, Eq. (15) makes
no sense as in this case only $\mathfrak{S}_{T}^{(N,m=1)}(\rho)$ exists). Our
goal is to characterize the maximum ergotropy value that can be achieved for a
given degree $m$ of the correlation. For this purpose we define the functional

$\displaystyle\mathcal{E}^{(N,m)}_{\max}(\rho;H):=\max_{\rho^{(N)}\in\mathfrak{S}_{T}^{(N,m)}(\rho)}\frac{\mathcal{E}(\rho^{(N)};H^{(N)})}{N}\;,$
(16)

that represents the maximum value of the ergotropy per-site that one can
extract from the system when it is in a joint state $\rho^{(N)}$ of
$\mathfrak{S}_{T}^{(N,m)}(\rho)$. From (15) it follows that (16) inherits the
monotonic behaviour

$\mathcal{E}^{(N,m-1)}_{\max}(\rho;H)\leq\mathcal{E}^{(N,m)}_{\max}(\rho;H)\;,\quad\forall
m\geq 2\;.$ (17)

Notice also that from (13) we trivially get

$\mathcal{E}^{(N,1)}_{\max}(\rho;H)=\frac{\mathcal{E}(\rho^{\otimes
N};H^{(N)})}{N}\;,$ (18)

that represents a lower bond for all the other
$\mathcal{E}^{(N,m)}_{\max}(\rho;H)$’s and which, in the asymptotic limit of
large $N$ approaches from below the classical threshold for the work
extraction from the single-site state $\rho$, given by the total ergotropy
function (1), i.e.

$\mathcal{E}^{(\infty,1)}_{\max}(\rho;H):=\lim_{N\rightarrow\infty}\mathcal{E}^{(N,1)}_{\max}(\rho;H)=\mathcal{E}_{\text{tot}}(\rho;H)\;.$
(19)

An upper bound for $\mathcal{E}^{(N,m)}_{\max}(\rho;H)$ can instead be easily
obtained from (157) and (9) which allows us to rewrite Eq. (16) in the
equivalent form

$\displaystyle\mathcal{E}^{(N,m)}_{\max}(\rho;H):=E(\rho;H)$ (20)
$\displaystyle\qquad\quad-\min_{\rho^{(N)}\in\mathfrak{S}_{T}^{(N,m)}(\rho)}\max_{U}\frac{E(U\rho^{(N)}U^{\dagger};H^{(N)})}{N}\;.$

Remembering then that $H^{(N)}$ is positive semidefinite, we can hence write

$\displaystyle\mathcal{E}^{(N,m)}_{\max}(\rho;H)\leq E(\rho;H)\;,\quad\forall
m\;.$ (21)

Observe next that, for fixed $N$ and $\rho$, the inequality (21) it is
certainly saturated for all $m$ larger than or equal to the square of the rank
of $\rho$, i.e.

$\mathcal{E}^{(N,m)}_{\max}(\rho;H)=E(\rho;H)\;,\quad\forall
m\geq\mbox{r}^{2}[\rho]\;.$ (22)

This is a consequence of the mononicity relation (17) and of the fact that
$\mathfrak{S}_{T}^{(N,\mbox{r}^{2}[\rho])}(\rho)$ contains at least the pure
state $|\mbox{GHZ}^{(N)}_{\rho}\rangle$ of Eq. (14) which allows one to set
equal to zero the negative term on the right-hand-side of Eq. (20) by choosing
$U$ to be unitary transformation that moves such vector into the ground state
of $H^{(N)}$ (remember that in our model the ground state eigenvalue is equal
to zero).

Determining how for fixed $N$ and $\rho$, one passes from the lower threshold
(18) to the upper threshold (22) by increasing the correlation parameter $m$,
is the focus of the present work. For this purpose, in the following sections,
we shall produce a series of lower bounds for
$\mathcal{E}^{(N,m)}_{\max}(\rho;H)$, that at least in same regimes allows one
to determine its exact value.

###  II.1 Summary of the main findings

In this section we anticipate the main results of the paper postposing the
derivations in the second part of the work.

Figure 2: Plot of the lower bounds of the asymptotic ($N\rightarrow\infty$)
per-site maximum many-body ergotropy $\mathcal{E}^{(\infty,m)}_{\max}(\rho;H)$
rescaled by $E(\rho;H)$, as a function of the BLR value $m$. Notice that in
agreement with (33), $\mathcal{E}^{(\infty,m)}_{\max}(\rho;H)/E(\rho;H)$
saturates to the optimal value $1$ already when $m$ coincide with the rank $r$
of the matrix $\rho$, a with a reduction of a factor 2 in dB unit with respect
to what one could expect from (22) that instead predicts this to happen for
$m=\mbox{r}^{2}[\rho]$. For smaller $m$ the our best estimation of the lower
bound (orange line in the plot) is the function $\mathcal{E}^{(\infty,m)}_{\bf
C}(\rho;H)/E(\rho;H)$ of Eq. (29); the blu curve is instead the lower bound
$\mathcal{E}^{(\infty,m)}_{\bf B}(\rho;H)/E(\rho;H)$ of Eq. (25). In this
example data were obtained taking $H$ to be the Hamiltonian of a quantum
system of dimension $d=40$ with randomly selected spaced level, and $\rho$ a
randomly selected, full rank ($r=d=40$) density matrix. The white, yellow, and
green stars in the plot represent, respectively, the ergotropy values achieved
on the factorized state $\rho^{\otimes N}$, the classically correlated state
$\rho_{cc}^{(N)}$ (see Eq. (105) below), and the GHZ state
$|\mbox{GHZ}^{(N)}_{\rho}\rangle$ of Eq. (14).

####  II.1.1 Lower bounds for $\mathcal{E}^{(N,m)}_{\max}(\rho;H)$ and
asymptotic saturation at $E(\rho;H)$ for $m=\mbox{r}[\rho]$

Our first observation is a lower bound for
$\mathcal{E}^{(N,m)}_{\max}(\rho;H)$ which, for $N$ sufficiently large, cover
the full spectrum of the BLR values, providing an interpolation between Eqs.
(18) and (22):

######  Proposition 1 (Preliminary bound).

Given $\eta\in]0,1[$ and a single-site density matrix $\rho$ of rank
$\mbox{r}[\rho]=r$, for $m\leq{r}^{2}$ and $N$ sufficiently large the
following inequality holds true

$\displaystyle\mathcal{E}^{(N,m)}_{\max}(\rho;H)$ (23) $\displaystyle\quad\geq
E(\rho;H)-\mathfrak{E}_{s_{\bf
B}}-(2\sqrt{2.01C\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}}})\;,$

where $\epsilon_{d}$ is maximum eigenvalue (3) of the single-site Hamiltonian
$H$, $\mathfrak{E}_{s}$ is the mean energy of a single-site thermal Gibbs
state with entropy $s\geq 0$ (see Eq. (189) in the Appendix),

$\displaystyle s_{\bf
B}:=\ln\left\lceil{r/\sqrt{m}}\;\right\rceil+N^{-1}\ln\sqrt{m}\;,$ (24)

while finally $C$ and $\alpha$ are finite, positive quantities which depend
upon $\rho$ and $H$ only.

The derivation of this result is reported in Sec. IV.1.1 by focusing on a
special sub-class of TI-MPO states (the block-wise purified density matrices
$\rho^{(N)}_{{\bf B}}$ introduced in Sec. III.4.2). Taking the
$N\rightarrow\infty$ limit for $\eta$ assigned and then sending the latter to
zero, from Eq. (27) we can extrapolate the following simplified inequality

$\mathcal{E}^{(\infty,m)}_{\max}(\rho;H):=\lim_{N\rightarrow\infty}\mathcal{E}^{(N,m)}_{\max}(\rho;H)\geq\mathcal{E}^{(\infty,m)}_{\bf{B}}(\rho;H)\;,$
(25)

which applies for all $m\leq\mbox{r}^{2}[\rho]$ with

$\displaystyle\mathcal{E}^{(\infty,m)}_{\bf{B}}(\rho;H)$ $\displaystyle:=$
$\displaystyle
E(\rho;H)-\mathfrak{E}_{\ln\left\lceil{\mbox{r}[\rho]/\sqrt{m}}\;\right\rceil}\;.$
(26)

Notice that since $\mathfrak{E}_{s}=0$ when $s=0$, for $m=\mbox{r}^{2}[\rho]$,
$\mathcal{E}^{(\infty,m)}_{\bf{B}}(\rho;H)$ reduces to $E(\rho;H)$ and Eq.
(25) leads to (22).

Our next result is a refinement of Proposition 1. To begin with we observe
that an improved version of Eq. (25) can be found in the low BLR regime
$m\leq\mbox{r}[\rho]$:

######  Proposition 2 (Low BLR regime).

Given $\eta\in]0,1[$ and a single-site density matrix $\rho$ of rank $r$, for
$m\leq{r}$ and $N$ sufficiently large the following inequality holds true

$\displaystyle\mathcal{E}^{(N,m)}_{\max}(\rho;H)$ (27) $\displaystyle\quad\geq
E(\rho;H)-\mathfrak{E}_{s_{\bf
C}}-(2\sqrt{2.01C\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}}})\;,$

where $\epsilon_{d}$, $C$, and $\alpha$ are as in Proposition 1 while now

$\displaystyle s_{\bf C}:=\ln\left\lceil{r/{m}}\;\right\rceil+N^{-1}\ln{m}\;.$
(28)

The derivation of this result is reported in Sec. IV.1.2: interestingly enough
this is done by focusing on a special class of TI-MPO states (i.e. the set
$\rho^{(N)}_{{\bf C}}$ of Sec. III.4.3) which represent classically correlated
(not entangled) density matrices of the system which can be chosen to be
diagonal in the energy eigenbasis. Similarly to what done for Proposition 1
the analysis simplifies when taking the proper $N\rightarrow\infty$ limit. In
this case from Eq. (27) we get the inequality

$\mathcal{E}^{(\infty,m)}_{\max}(\rho;H)\geq\mathcal{E}^{(\infty,m)}_{\bf{C}}(\rho;H)\;,\qquad\forall
m\leq\mbox{r}[\rho]\;,$ (29)

which for

$\displaystyle\mathcal{E}^{(\infty,m)}_{\bf{C}}(\rho;H)$ $\displaystyle:=$
$\displaystyle
E(\rho;H)-\mathfrak{E}_{\ln\left\lceil{\mbox{r}[\rho]/{m}}\right\rceil}\;,$
(30)

represents a clear improvement with respect to the constraint imposed by Eq.
(26) – see Fig. 2. Notice in particular that since
$\mathcal{E}^{(\infty,\mbox{r}[\rho])}_{\bf{C}}(\rho;H)=E(\rho;H)$, Eq. (29)
predicts a saturation of the upper bound (21) already for $m=\mbox{r}[\rho]$,
i.e. for values of $m$ which are much smaller than those suggested by (22).
Even more interesting is the fact that such asymptotic saturation is attained
with states (the density matrices $\rho^{(N)}_{{\bf C}}$) which, as already
mentioned, contain no quantum correlations and which can be chosen to the
diagonal with respect to the energy eigenbasis. The first property of the
$\rho^{(N)}_{{\bf C}}$’s implies that, at least in the asymptotic regime of
infinitely many sites, classical correlations are sufficient to enable full
work extraction from the system. The second property instead tell us that this
can be done via incoherent ergotropy Francica  _et al._ (2020); Çakmak (2020);
Touil  _et al._ (2021).

To clarify the asymptotic attainability of the upper bound (21) for large $m$
we present a final statement that applies for BLR values that larger than or
equal to $\mbox{r}[\rho]$:

######  Proposition 3 (High BLR regime).

Given $N$ integer and a single-site density matrix $\rho$ of rank $r$, for
$m\in\\{r,\cdots,r^{2}\\}$ the following inequality holds

$\mathcal{E}^{(N,m)}_{\max}(\rho;H)\geq E(\rho;H)-\frac{\Delta}{N}\;\;,$ (31)

with $\Delta$ being a positive constant term which depends upon the spectra of
$\rho$ and $H$.

The proof of this result which at variance with Propositions 1 and 2 does not
require to have $N$ large, is given in Sec. IV.2 by focusing on the special
class of (non pure) quantum states $\rho_{\bf A}^{(N)}$ that correspond to
proper mixtures of multi-sites GHZ states (see Sec. III.4.1). In that section
the quantity $\Delta$ appearing in Eq. (31) is also identified with the gap
between the single-site mean energy, and the single-site ergotropy, i.e.

$\displaystyle\Delta:={E}(\rho,H)-{\cal E}(\rho,H)\;.$ (32)

A direct consequence of Propositions 3 is the universal identity

$\displaystyle\mathcal{E}^{(\infty,m)}_{\max}(\rho;H)=E(\rho;H)\;,\quad
m\geq\mbox{r}[\rho]\;,$ (33)

which clearly superseds (25) in the high BRL regime, confirming the saturation
of the bound (21) observed in Eq. (29) for $m=\mbox{r}[\rho]$.

As evident e.g. from Fig. 2, Propositions 2 and 3 provide our best estimations
for $\mathcal{E}^{(N,m)}_{\max}(\rho;H)$. For large enough $N$ and $m\geq r$
they are optimal since lead to the exact evaluation of
$\mathcal{E}^{(\infty,m)}_{\max}(\rho;H)$. For small $m$ on the contrary the
constraint posed by Eq. (27) is certainly suboptimal: to see this observe for
instance that for $m=1$ we get

$\displaystyle\mathcal{E}^{(\infty,1)}_{\bf{C}}(\rho;H)=E(\rho;H)-\mathfrak{E}_{\ln\left\lceil{\mbox{r}[\rho]}\right\rceil}\;,$
(34)

which is clearly not larger than the exact value (19) due the fact that
$\mathfrak{E}_{\ln\left\lceil{\mbox{r}[\rho]}\right\rceil}$ is always larger
than or equal to the thermal energy $\mathfrak{E}_{S(\rho)}$ associated with
the Gibbs state that is iso-entropic with $\rho$ (same considerations apply of
course for $\mathcal{E}^{(\infty,1)}_{\bf{B}}(\rho;H)$). An improvement for
small $m$ can however be obtained by invoking an heuristic (not rigorous)
argument that we shall discuss in Remark 3 of Sec. IV.1.2. This suggests that
for $m\ll r$ it should be possible to identify a special sub-class of states
$\rho^{(N)}_{{\bf C}}$ we employed in the derivation of Proposition 2 which
allow us to replace $s_{\bf C}$ in Eq. (27) with the improved value

$\displaystyle\left.s_{\bf C}\right|_{\rm(heu)}:=S(\rho)-\frac{N-1}{N}\ln
m\;,$ (35)

and (34) with

$\left.\mathcal{E}^{(\infty,1)}_{\bf{C}}(\rho;H)\right|_{\rm(heu)}=E(\rho;H)-\mathfrak{E}_{S(\rho)}=\mathcal{E}_{\text{tot}}(\rho;H)\;.$
(36)

####  II.1.2 Lower bound for $\frac{\mathcal{E}(\rho^{\otimes N};H^{(N)})}{N}$

As mentioned in the introduction a secondary result our analysis is to set on
rigorous ground of the typicality argument first formulated by Alicki and
Fannes Alicki and Fannes (2013) for the maximum work extractable for a set of
many identical quantum systems. Specifically we show that

######  Corollary 1.

Given $\eta\in]0,1[$ and a single-site density matrix $\rho$, for $N$
sufficiently large the following inequality holds true

$\frac{{\cal E}(\rho^{\otimes
N};H^{(N)})}{N}\geq\mathcal{E}_{\text{tot}}(\rho;H)-(2\sqrt{2.01C\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}}})\;,$
(37)

with $\epsilon_{d}$, $C$, and $\alpha$ as in Proposition 1.

The proof of this result can be seen as a special instance of Proposition 1
and it is given in Sec. IV.3.

##  III TI-MPO representation of translationally invariant states

This section is dedicated to clarify some technical aspects of the TI-MPO
representation and to present explicit examples of elements of
$\mathfrak{S}_{T}^{(N,m)}(\rho)$ to be used to construct our lower bounds for
$\mathcal{E}^{(N,m)}_{\max}(\rho;H)$. First of all in Sec. III.1 we formally
show that any translational invariant state of the $N$ sites, admits a TI-MPO
representation (10) for a proper choice of the tensor $\mathbb{A}$; then in
Sec. III.2 we show that the definition of the BLR $\mathfrak{br}[\rho^{(N)}]$
defined in Eq. (12) is unaffected by the choice of the local basis
$\\{\ket{i};i=1,\cdots,d\\}$ entering in the TI-MPO representation (10); in
Sec. III.3 we prove instead that the GHZ-like states of the form (14) belong
to the set $\mathfrak{S}_{T}^{(N,\mbox{r}^{2}[\rho])}(\rho)$; finally in Sec.
III.4 we provide an explicit TI-MPO representation for three different
families of correlated states in $\mathfrak{S}_{T}^{(N)}(\rho)$.

###  III.1 Existence of TI-MPO representation

In order to show that the elements $\rho^{(N)}$ of $\mathfrak{S}_{T}^{(N)}$
can be expressed in the TI-MPO representation of Eq. (10), let start from a
not necessarily TI - MPO representation of such an element, i.e.

$\displaystyle\rho^{(N)}$ $\displaystyle=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{i_{1}\dots
i_{N}}\bra{j_{1}\dots j_{N}}$ (38)
$\displaystyle\;\times\Tr[{A}_{1}^{i_{1},j_{1}}\dots A_{N}^{i_{N},j_{N}}]\;,$

where for $k=1,\cdots,N$, $\\{A_{k}^{i,j}|{i,j\in 1,\cdots d}\\}$ are a set of
$d^{2}$, $\bar{M}\times\bar{M}$ matrices associated to the $k$-th site, which
always exists for $\bar{M}$ sufficiently large Parker  _et al._ (2020). Now
exploiting the translational invariance property (7) and the ciclicity of the
trace, it follows that for all $\ell=0,\cdots,N-1$ we can also write

$\displaystyle\rho^{(N)}$ $\displaystyle=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{i_{1}\dots
i_{N}}\bra{j_{1}\dots j_{N}}$ (39)
$\displaystyle\;\times\Tr[{A}_{1\oplus\ell}^{i_{1},j_{1}}\dots
A_{N\oplus\ell}^{i_{N},j_{N}}]\;,$

with $\oplus$ representing the sum modulus $N$. Hence we get

$\displaystyle\rho^{(N)}$ $\displaystyle=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{i_{1}\dots
i_{N}}\bra{j_{1}\dots j_{N}}$ (40)
$\displaystyle\;\times\frac{1}{N}\sum_{\ell=0}^{N-1}\Tr[{A}_{1\oplus\ell}^{i_{1},j_{1}}\dots
A_{N\oplus\ell}^{i_{N},j_{N}}]\;,$

which is of the form Eq. (10) by identifying $M=\bar{M}N$ and taking the
matrices ${A}^{i,j}$ as

$\displaystyle
A^{i,j}=\frac{1}{{N}^{1/N}}\sum_{\ell=1}^{N}{A}_{\ell}^{i,j}\otimes|v_{\ell}\rangle\langle
v_{\ell\oplus 1}|\;,$ (41)

with $\\{|v_{\ell}\rangle;\ell=1,\cdots,N-1\\}$ being orthonormal vectors of
an auxiliary $N$-dimensional vector space (hereafter, in the writing of the
MPO matrices we adopt the Dirac notation).

###  III.2 Invariance of the BLR with respect to the choice of the local basis

Here we show that the choice of the local basis $\\{|i\rangle\\}$ entering the
TI-MPO representation (10) does not affect the definition BLR (12) of a state
$\rho^{(N)}\in\mathfrak{S}_{T}^{(N)}$. To see this consider
$\\{\ket{\phi_{i}};i=1,\cdots,d\\}$ a new basis of ${\cal H}$ connected to
$\\{\ket{i};i=1,\cdots,d\\}$ via a unitary transformation $U$, i.e.

$\displaystyle\ket{\phi_{i}}=U\ket{i}\;,\forall i=1,\cdots,d\;.$ (42)

Replacing this into (10) we get

$\displaystyle\rho^{(N)}$ $\displaystyle=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}U^{{\dagger}\otimes
N}\ket{\phi_{i_{1}}\dots\phi_{i_{N}}}\bra{\phi_{j_{1}}\dots\phi_{j_{N}}}U^{\otimes
N}$ (43) $\displaystyle\;\times\Tr[{A}^{i_{1},j_{1}}\dots A^{i_{N},j_{N}}]\;,$
$\displaystyle=$
$\displaystyle\sum_{i^{\prime}_{1},j^{\prime}_{1},\dots,i^{\prime}_{N},j^{\prime}_{N}}\ket{\phi_{i^{\prime}_{1}}\dots\phi_{i^{\prime}_{N}}}\bra{\phi_{j^{\prime}_{1}}\dots\phi_{j^{\prime}_{N}}}$
$\displaystyle\;\times\Tr[{B}^{i^{\prime}_{1},j^{\prime}_{1}}\dots
B^{i^{\prime}_{N},j^{\prime}_{N}}]\;,$

where for $i^{\prime},j^{\prime}=1,\cdots,d$

$\displaystyle{B}^{i^{\prime},j^{\prime}}=\sum_{i,j=1}^{d}\langle\phi_{i^{\prime}}|U^{\dagger}|\phi_{i}\rangle{A}^{i,j}\langle\phi_{j}|U^{\dagger}|\phi_{j^{\prime}}\rangle$
(44)

forming a new set of $d^{2}$, $M\times M$ matrices. Equation (44) creates a
one-to-one correspondence between the TI-MPO representations of $\rho^{(N)}$
constructed with the basis $\\{\ket{\phi_{i}};i=1,\cdots,d\\}$ with bond link
$M$, and the TI-MPO representations associated with the basis
$\\{\ket{i};i=1,\cdots,d\\}$ with the same bond link: accordingly the BLR (12)
computed with respect to those two basis will produce the same value. Notice
in particular that this allows us to identify $\ket{i}$ appearing in (10) with
the eigenvectors $|\lambda_{i}\rangle$ of the local density matrix $\rho$ of
Eq. (8), a choice that, via Eq. (11), will allow us to identify the associated
eigenvalues as

$\displaystyle\Tr[{A}^{i,j}{\bar{A}}^{N-1}]=\lambda_{i}\;\delta_{i,j}\;.$ (45)

###  III.3 TI-MPO representation for GHZ-like states

Consider $\rho$ a single-site density matrix of rank $\mbox{r}[\rho]$

$\displaystyle\rho=\sum_{i=1}^{\mbox{r}[\rho]}\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|\;,$
(46)

with $\\{|\lambda_{i}\rangle,i=1,\cdots,d\\}$ its eigenvectors, and
$\\{\lambda_{i};i=1,\cdots,\mbox{r}[\rho]\\}$ the corresponding non-zero
eigenvalues (for $i\geq\mbox{r}[\rho]+1$ of course
$\rho|\lambda_{i}\rangle=0$). Consider also the associated GHZ-like state (14)
which is explicitly translationally invariant. Here we prove that
$|\mbox{GHZ}^{(N)}_{\rho}\rangle$ belongs to the set
$\mathfrak{S}_{T}^{(N,\mbox{r}^{2}[\rho])}(\rho)$, i.e.

$\displaystyle|\mbox{GHZ}^{(N)}_{\rho}\rangle\langle\mbox{GHZ}^{(N)}_{\rho}|\in\mathfrak{S}_{T}^{(N,\mbox{r}^{2}[\rho])}(\rho)\;,$
(47)

or equivalently that its BLR (12) is smaller or equal than
$\mbox{r}^{2}[\rho]$, i.e.

$\displaystyle\mathfrak{br}\left[|\mbox{GHZ}^{(N)}_{\rho}\rangle\langle\mbox{GHZ}^{(N)}_{\rho}|\right]\leq\mbox{r}^{2}[\rho]\;.$
(48)

To show this fact it is sufficient to exhibit a TI-MPO representation (10) of
$|\mbox{GHZ}^{(N)}_{\rho}\rangle\langle\mbox{GHZ}^{(N)}_{\rho}|$ with a tensor
$\mathbb{A}\in\mathbb{C}^{d\times d\times M\times M}$ characterized by
$M=\mbox{r}^{2}[\rho]$. Exploiting the observation of Sec. III.2 we shall
focus on representations which uses the eigenvectors
$\\{|\lambda_{i}\rangle;i=1,\cdots,d\\}$ of $\rho$ as local basis
$\\{|i\rangle;i=1,\cdots,d\\}$. Then we notice that the solution is provided
by the identity

$\displaystyle|\mbox{GHZ}^{(N)}_{\rho}\rangle\langle\mbox{GHZ}^{(N)}_{\rho}|$
$\displaystyle=$
$\displaystyle\sum_{i,j=1}^{\mbox{r}[\rho]}\sqrt{\lambda_{i}\lambda_{j}}\;\ket{\lambda_{i}\dots\lambda_{i}}\bra{\lambda_{j}\dots\lambda_{j}}$
(49) $\displaystyle=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{\lambda_{i_{1}}\dots\lambda_{i_{N}}}\bra{\lambda_{j_{1}}\dots\lambda_{j_{N}}}$
$\displaystyle\;\times\Tr[{A}^{i_{1},j_{1}}\dots A^{i_{N},j_{N}}]\;,$

where defining $\\{|v_{i}\rangle;i=1,\cdots,\mbox{r}[\rho]\\}$ an orthonormal
basis of an auxiliary vector space ${\cal H}_{A}$ of dimension
$\mbox{r}[\rho]$, we took ${A}^{i,j}$ to be the
$\mbox{r}^{2}[\rho]\times\mbox{r}^{2}[\rho]$ matrices

$\displaystyle{A}^{i,j}:=\left\\{\begin{array}[]{ll}\left({\lambda_{i}\lambda_{j}}\right)^{\frac{1}{2N}}\;|v_{i}\rangle\langle
v_{i}|\otimes|v_{j}\rangle\langle v_{j}|\;,&\forall
i,j\leq\mbox{r}[\rho]\;,\\\ \\\ 0\;,&\mbox{otherwise.}\end{array}\right.$ (53)

We remark that the argument we presented here doesn’t ensure that
$\mathfrak{br}[|\mbox{GHZ}^{(N)}_{\rho}\rangle\langle\mbox{GHZ}^{(N)}_{\rho}|]$
coincides with $\mbox{r}^{2}[\rho]$, but (as implied by Eq. (48)) only that
the former is an upper bound for the latter. In other words, formally
speaking, we cannot exclude that
$|\mbox{GHZ}^{(N)}_{\rho}\rangle\langle\mbox{GHZ}^{(N)}_{\rho}|$ belongs to
$\mathfrak{S}_{T}^{(N,m)}(\rho)$ for same $m<\mbox{r}^{2}[\rho]$ (we
conjecture however that this is not the case).

###  III.4 TI-MPO representation of three families of correlated states of
$\mathfrak{S}_{T}^{(N)}(\rho)$

In this section we present an explicit TI-MPO representation for three
different families of states in $\mathfrak{S}_{T}^{(N)}(\rho)$ which, as
pictorially depicted in Fig. 3, allows us span various BLR values. All these
examples are constructed starting from the set of the non-null eigenvalues of
the single-site density matrix $\rho$ of the model, i.e.

$\displaystyle{\mathfrak{K}}[\rho]:=\\{\lambda_{i};i=1,\cdots,\mbox{r}[\rho]\\}\;,$
(54)

with $\lambda_{i}$ and $\mbox{r}[\rho]$ defined as in Eq. (46), and dividing
it into smaller, non-empty subsets. Specifically, given $L$ integer greater
than or equal to $1$ and no larger than $\mbox{r}[\rho]$, we consider a
collection

$\displaystyle{\mathfrak{P}}:=\\{{\mathfrak{K}}_{1},{\mathfrak{K}}_{2},\cdots,{\mathfrak{K}}_{L}\\}\;,$
(55)

of $L$ non-overlapping, non-null subsets of ${\mathfrak{K}}[\rho]$, which
provides a partition of such set, i.e.

$\displaystyle{\mathfrak{K}}_{\ell}\cap{\mathfrak{K}}_{\ell^{\prime}}$
$\displaystyle=$ $\displaystyle\O\;,\qquad\forall\ell\neq\ell^{\prime}\;,$
(56) $\displaystyle\bigcup_{\ell=1}^{M}{\mathfrak{K}}_{\ell}$ $\displaystyle=$
$\displaystyle{\mathfrak{K}}[\rho]\;.$ (57)

In what follow we shall use the symbol $\\#{\mathfrak{K}}_{\ell}$ to represent
the cardinality of ${\mathfrak{K}}_{\ell}$ and indicate with $\\#_{\max}$ and
$|{\mathfrak{P}}|^{(p)}$ their maximum value and the sum of their $p$-powers,
i.e. the quantities

$\displaystyle\\#_{\max}$ $\displaystyle:=$
$\displaystyle\max_{\ell=1,\cdots,L}\left(\\#{\mathfrak{K}}_{\ell}\right)\;,$
(58) $\displaystyle|{\mathfrak{P}}|^{(p)}$ $\displaystyle:=$
$\displaystyle\sum_{\ell=1}^{L}(\\#{\mathfrak{K}}_{\ell})^{p}\;,$ (59)

which are related by the inequality

$\displaystyle L\leq|{\mathfrak{P}}|^{(p)}\leq(\\#_{\max})^{p}L\;.$ (60)

By construction the $\\#{\mathfrak{K}}_{\ell}$’s are all greater than or equal
to $1$ and sum up to the rank of $\rho$, i.e.

$\displaystyle\sum_{\ell=1}^{L}\\#{\mathfrak{K}}_{\ell}=\mbox{r}[\rho]\;,$
(61)

which implies

$\displaystyle\frac{\mbox{r}[\rho]}{L}\leq$ $\displaystyle\\#_{\max}$
$\displaystyle\leq\mbox{r}[\rho]-L+1\;,$ (62)
$\displaystyle\mbox{r}[\rho]\leq$ $\displaystyle|{\mathfrak{P}}|^{(p)}$
$\displaystyle\leq\mbox{r}^{p}[\rho]\;.$ (63)

Notice that if $L$ is an exact divisor of $\mbox{r}[\rho]$ one can force all
the subsets $\mathfrak{K}_{\ell}$ to have the same cardinality, i.e.

$\displaystyle\\#{\mathfrak{K}}_{\ell}={\mbox{r}[\rho]}/{L}\;,\qquad\forall\ell\in
1,\cdots,L\;,$ (64)

corresponding to an exact saturation of the lower bound in (62): when this
happens one gets $|{\mathfrak{P}}|^{(p)}=\mbox{r}^{p}[\rho]/{L^{p-1}}$, and we
say that the associated ${\mathfrak{P}}$ is an uniform partition of
${\mathfrak{K}}[\rho]$. In case $L$ is not an exact divisor of
$\mbox{r}[\rho]$ we say instead that ${\mathfrak{P}}$ is an almost uniform
partition of ${\mathfrak{K}}[\rho]$ if

$\displaystyle\\#_{\max}=\left\lceil{\mbox{r}[\rho]/L}\right\rceil\;,$ (65)

leading to
$|{\mathfrak{P}}|^{(p)}\leq(\left\lceil{\mbox{r}[\rho]/L}\right\rceil)^{p}L$.

Next we introduce a labelling for the elements of each individual subsets of
${\mathfrak{P}}$; specifically, given $\ell=1,\cdots,L$, we indicate with
$\lambda_{k,\ell}$ the $k$-th element of ${\mathfrak{K}}_{\ell}$, so that

$\displaystyle{\mathfrak{K}}_{\ell}=\\{\lambda_{k,\ell};k=1,\cdots,\\#{\mathfrak{K}}_{\ell}\\}\;.$
(66)

Notice that, thanks to (56) and (57), for each $i=1,\cdots,\mbox{r}[\rho]$
there is a unique choice of the indexes $\ell=1,\cdots L$ and
$k=1,\cdots,\\#{\mathfrak{K}}_{\ell}$ that identifies the $i$-th eigenvalue of
$\rho$ as the $k$ element of the subset ${\mathfrak{K}}_{\ell}$: in the
following we shall exploit this one-to-one correspondence indicating with the
symbols $k^{[i]}$ and $\ell^{[i]}$ those special values via the mapping

$\displaystyle i\longmapsto
k^{[i]},\ell^{[i]}\quad\mbox{s.t.}\quad\lambda_{i}=\lambda_{k^{[i]},\ell^{[i]}}\;.$
(67)

Finally for all $\ell=1,\cdots,L$, we introduce the quantity

$\displaystyle\mathscr{S}_{\ell}:=\sum_{k=1}^{\\#{\mathfrak{K}}_{\ell}}\lambda_{k,\ell}\;,$
(68)

to gauge the statistical weight of the set $\mathfrak{K}_{\ell}$: these terms
are clearly positive and fulfils the normalization condition

$\displaystyle\sum_{\ell=1}^{L}\mathscr{S}_{\ell}$ $\displaystyle=$
$\displaystyle 1\;.$ (69)

Figure 3: Graphical representation of the families of states of $\mathfrak{S}_{T}^{(N)}(\rho)$ defined in Sec. III.4. Specifically by the black trajectory represents the set of convex convolution of GHZ-like states $\rho^{(N)}_{\bf A}$ defined in Eq. (75) of Sec. III.4.1, which connect the pure state $|\mbox{GHZ}^{(N)}_{\rho}\rangle$ of Eq. (14) (green star in the figure) with its classically correlated counterpart $\rho^{(N)}_{\rm cc}$ defined in Eq. (77) and represented by the yellow star in the picture; the blue trajectory is instead the set of block-wise purified states $\rho^{(N)}_{\bf B}$ defined in Eq. (96) of Sec. III.4.2, which interpolates from $|\mbox{GHZ}^{(N)}_{\rho}\rangle$ to the fully factorized element $\rho^{\otimes N}$ (white star); finally the orange trajectory represents the set of classically correlated states $\rho^{(N)}_{\bf C}$ defined in Eq. (105) of Sec. III.4.3, which interpolates from $\rho^{\otimes N}$ to $\rho^{(N)}_{\rm cc}$. The grey lines represent the subsets $\mathfrak{S}_{T}^{(N,m)}(\rho)$ of Eq. (II): the larger is the set represent the higher is the associated values of the BLR $\mathfrak{br}[\rho^{(N)}]$ of the corresponding states. | $\rho^{(N)}_{\bf A}$ | $\rho^{(N)}_{\bf B}$ | $\rho^{(N)}_{\bf C}$
---|---|---|---
rank |  $L$ | $\left(\\#_{\max}\right)^{N}$ | $|{\mathfrak{P}}|^{(N)}$
BLR | $\leq|{\mathfrak{P}}|^{(2)}$ | $\leq L^{2}$ |  $\leq L$
Table 1: Rank and BLR (12) for the states $\rho^{(N)}_{\bf A}$,
$\rho^{(N)}_{\bf B}$, and $\rho^{(N)}_{\bf C}$ of
$\mathfrak{S}_{T}^{(N)}(\rho)$ defined in Sec. III.4. The integer parameter
$L=1,\cdots,\mbox{r}[\rho]$ counts the number of elements of the partition
${\mathfrak{P}}$ which we used to divide the set of positive eigenvalues
${\mathfrak{K}}[\rho]$ of the single-side state (see Eq. (55)); $\\#_{\max}$
is the maximum among the cardinalities of the elements of ${\mathfrak{P}}$
(see Eq. (58)); while for $p$ integer, $|{\mathfrak{P}}|^{(p)}$ is the sum of
the $p$-powers of the cardinalities (see Eq. (59)). For uniform partitions we
get $\\#_{\max}=\mbox{r}[\rho]/L$ and
$|{\mathfrak{P}}|^{(p)}=\mbox{r}^{p}[\rho]/L^{p-1}$ with $\mbox{r}[\rho]$
being the rank of the single-site density matrix.

####  III.4.1 TI-MPO representation for convex convolutions of non overlapping
GHZ-like states

The first family of states of $\mathfrak{S}_{T}^{(N)}(\rho)$ we consider is a
generalization of the GHZ-like construction introduced in Sec. III.3 (this set
will be used to prove Proposition 3). In particular for element
$\mathfrak{K}_{\ell}$ of the the partition (55) we define a corresponding GHZ-
like state via the identity

$\displaystyle|\mbox{GHZ}^{(N)}_{\rho_{\ell}}\rangle:=\frac{1}{\sqrt{\mathscr{S}_{\ell}}}\sum_{k=1}^{\\#{\mathfrak{K}}_{\ell}}\sqrt{\lambda_{k,\ell}}|\lambda_{k,\ell}\dots\lambda_{k,\ell}\rangle\;.$
(70)

This a translational invariant state whose associated single-site local
density matrix (8) is given by

$\displaystyle\rho_{\ell}:=\frac{1}{\mathscr{S}_{\ell}}\sum_{k=1}^{\\#{\mathfrak{K}}_{\ell}}\lambda_{k,\ell}|\lambda_{k,\ell}\rangle\langle_{k,\ell}\lambda_{k,\ell}|\;,$
(71)

which has a rank $\mbox{r}[\rho_{\ell}]$ that corresponds to the cardinality
of ${\mathfrak{K}}_{\ell}$, i.e.

$\displaystyle\mbox{r}[\rho_{\ell}]=\\#{\mathfrak{K}}_{\ell}\;.$ (72)

From Eq. (57) and (69) it follows that the vectors (70) are orthonormal, i.e.

$\displaystyle\langle\mbox{GHZ}^{(N)}_{\rho_{\ell}}|\mbox{GHZ}^{(N)}_{\rho_{\ell^{\prime}}}\rangle=\delta_{\ell,\ell^{\prime}}\;,$
(73)

and provide a decomposition of the state (14) via the identity

$\displaystyle|\mbox{GHZ}^{(N)}_{\rho}\rangle=\sum_{\ell=1}^{L}\sqrt{\mathscr{S}_{\ell}}\;|\mbox{GHZ}^{(N)}_{\rho_{\ell}}\rangle\;.$
(74)

Furthermore by direct inspection one can easily verify that the convex
convolution of rank $L$ obtained by mixing the vectors
$|\mbox{GHZ}^{(N)}_{\rho_{\ell}}\rangle$ with statistical weight
$\mathscr{S}_{\ell}$, i.e. the density matrix

$\displaystyle\rho^{(N)}_{\bf
A}:=\sum_{\ell=1}^{L}\mathscr{S}_{\ell}\;|\mbox{GHZ}^{(N)}_{\rho_{\ell}}\rangle\langle\mbox{GHZ}^{(N)}_{\rho_{\ell}}|\;,$
(75)

is an element of $\mathfrak{S}_{T}^{(N)}(\rho)$, i.e. it admits $\rho$ of Eq.
(46) as single-site reduced density matrix (ultimately this is a consequence
of the identity $\sum_{\ell=1}^{L}\mathscr{S}_{\ell}\rho_{\ell}=\rho$), and
has a rank equal to $L$ as reported in Table 1. Notice in particular that as
$L$ varies from the extremal cases $L=1$ and $L=\mbox{r}[\rho]$ the states
(75) interpolate from the GHZ-like configuration to its classically correlated
counterpart, i.e.

$\displaystyle\rho^{(N)}_{{\bf A}}\Big{|}_{L=1}$ $\displaystyle=$
$\displaystyle|\mbox{GHZ}^{(N)}_{\rho}\rangle\langle\mbox{GHZ}^{(N)}_{\rho}|\;,$
(76) $\displaystyle\rho^{(N)}_{{\bf A}}\Big{|}_{L=\mbox{r}[\rho]}$
$\displaystyle=$ $\displaystyle\rho^{(N)}_{\rm
cc}:=\sum_{i=1}^{\mbox{r}[\rho]}\lambda_{i}\ket{\lambda_{i}\dots\lambda_{i}}\bra{\lambda_{i}\dots\lambda_{i}},$
(77)

(see Fig. 3).

We now show that the BLR of $\rho^{(N)}_{\bf A}$ fulfils the inequality

$\displaystyle\mathfrak{br}[\rho^{(N)}_{\bf
A}]\leq|{\mathfrak{P}}|^{(2)}=\sum_{\ell=1}^{L}(\\#{\mathfrak{K}}_{\ell})^{2}\;,$
(78)

or equivalently that

$\displaystyle\rho^{(N)}_{\bf
A}\in\mathfrak{S}_{T}^{(N,|{\mathfrak{P}}|^{(2)})}(\rho)\;.$ (79)

For the trivial choice $L=1$ where ${\mathfrak{P}}$ contains only
${\mathfrak{K}}[\rho]$ as unique element, $\rho^{(N)}_{\bf A}$ corresponds to
the GHZ-like state (14) which purifies $\rho$,
$|{\mathfrak{P}}|^{(2)}=\mbox{r}^{2}[\rho]$, and Eq. (79) reduces to the
property (47) which we proved in the previous section. To show that (79) holds
also for $L>1$ we adopt a similar scheme and present an explicit TI-MPO
decomposition

$\displaystyle\rho^{(N)}_{\bf A}$ $\displaystyle=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{\lambda_{i_{1}}\dots\lambda_{i_{N}}}\bra{\lambda_{j_{1}}\dots\lambda_{j_{N}}}$
(80) $\displaystyle\;\times\Tr[{A}^{i_{1},j_{1}}_{{\bf A}}\dots
A^{i_{N},j_{N}}_{{\bf A}}]\;,$

that employs the eigenvectors $\\{|\lambda_{i}\rangle;i=1,\cdots,d\\}$ of
$\rho$ as local basis and which explicitly uses matrices
$|{\mathfrak{P}}|^{(2)}\times|{\mathfrak{P}}|^{(2)}$ matrices ${A}^{i,j}_{{\bf
A}}$. For this purpose observe first that given $\ell=1,\cdots,L$, the state
$|\mbox{GHZ}^{(N)}_{\rho_{\ell}}\rangle$ admits a TI-MPO representation

$\displaystyle|\mbox{GHZ}^{(N)}_{\rho_{\ell}}\rangle\langle\mbox{GHZ}^{(N)}_{\rho_{\ell}}|$
$\displaystyle=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{\lambda_{i_{1}}\dots\lambda_{i_{N}}}\bra{\lambda_{j_{1}}\dots\lambda_{j_{N}}}$
(81) $\displaystyle\;\times\Tr[{A}_{\ell}^{i_{1},j_{1}}\dots
A_{\ell}^{i_{N},j_{N}}]\;,$

with

${A}_{\ell}^{i,j}:=\left\\{\begin{array}[]{ll}\left({\frac{\lambda_{i}\lambda_{j}}{\mathscr{S}^{2}_{\ell}}}\right)^{\frac{1}{2N}}\\!\\!|v_{k^{[i]}}^{(\ell)}\rangle\langle
v_{k^{[i]}}^{(\ell)}|\otimes|v_{k^{[j]}}^{(\ell)}\rangle\langle
v_{k^{[j]}}^{(\ell)}|,&\forall\lambda_{i},\lambda_{j}\in{\mathfrak{K}}_{\ell},\\\
\\\ \\\ \\\ 0\;,&\mbox{otherwise,}\end{array}\right.$ (82)

where $k^{[i]}$ and $k^{[j]}$ are defined by the mapping (67), and where
$\\{|v_{k}^{(\ell)}\rangle;k=1,\cdots,\\#{\mathfrak{K}}_{\ell}\\}$ form an
orthonormal basis of an auxiliary vector space ${\cal H}_{\ell}$ of dimension
$\mbox{r}[\rho_{\ell}]$ (notice that the $A_{\ell}^{i,j}$’s operate on two
copies of ${\cal H}_{\ell}$, i.e. on ${\cal H}_{\ell}^{\otimes 2}$, and are
hence $(\\#{\mathfrak{K}}_{\ell})^{2}\times(\\#{\mathfrak{K}}_{\ell})^{2}$).
Equation (80) can now be obtained by identifying the ${A}^{i,j}_{{\bf A}}$’s
as proper direct sums over the index $\ell$ of the matrices (82).
Specifically, let first observe that the auxiliary space ${\cal H}_{{\bf
A}}:=\bigoplus_{\ell=1}^{L}{\cal H}_{\ell}^{\otimes 2}$ admits
$\bigcup_{\ell=1}^{L}\\{|v_{k}^{(\ell)}\rangle\otimes|v_{k^{\prime}}^{(\ell)}\rangle;k,k^{\prime}=1,\cdots,\\#{\mathfrak{K}}_{\ell}\\}$
as an orthonormal set and has dimension

$\displaystyle\mbox{dim}[{\cal H}_{{\bf
A}}]=\sum_{\ell=1}^{L}\mbox{dim}^{2}[{\cal
H}_{\ell}]=\sum_{\ell=1}^{L}(\\#{\mathfrak{K}}_{\ell})^{2}=|{\mathfrak{P}}|^{(2)}\;.$
(83)

On such space we then introduce the matrices

$\displaystyle{A}^{i,j}_{{\bf A}}$ $\displaystyle:=$
$\displaystyle\bigoplus_{\ell=1}^{L}{\mathscr{S}_{\ell}}^{\tfrac{1}{N}}\;{A}_{\ell}^{i,j}\;,$
(84)

which can be equivalently expressed as

$\displaystyle{A}^{i,j}_{{\bf A}}$ $\displaystyle=$
$\displaystyle(\sqrt{\lambda_{i}\lambda_{j}})^{1/N}\;\delta_{\ell^{[i]},\ell^{[j]}}$
$\displaystyle\times|v_{k^{[i]}}^{(\ell^{[i]})}\rangle\langle
v_{k^{[i]}}^{(\ell^{[i]})}|\otimes|v_{k^{[j]}}^{(\ell^{[j]})}\rangle\langle
v_{k^{[j]}}^{(\ell^{[j]})}|\;,$

with $\ell^{[i]}$ and $\ell^{[j]}$ defined by the mapping (67) and with
$\delta_{\ell^{[i]},\ell^{[j]}}$ being the Kronecker delta symbol that forces
$\ell^{[i]}=\ell^{[j]}$. The identity (80) finally follows by observing that

$\Tr[{A}^{i_{1},j_{1}}_{{\bf A}}\dots A^{i_{N},j_{N}}_{{\bf
A}}]=\sum_{\ell=1}^{L}\mathscr{S}_{\ell}\Tr[{A}^{i_{1},j_{1}}_{\ell}\dots
A^{i_{N},j_{N}}_{\ell}]\;,$ (86)

and from (81) and (75).

####  III.4.2 TI-MPO representation for block-wise purified density matrices

The second family of elements of $\mathfrak{S}_{T}^{(N)}(\rho)$ we consider
allow us to connect the completely uncorrelated $\rho^{\otimes N}$ to the pure
GHZ-like state configuration $|\mbox{GHZ}^{(N)}_{\rho}\rangle$ – see
trajectory $\rho^{(N)}_{{\bf B}}$ in Fig. 3. We dub these states block-wise
purified density matrices and define them by means of an explicit TI-MPO
representation. They will be used in the derivation of Proposition 1.

The starting point of the analysis is again the partition (55) of the non-null
eigenvalues of $\rho$ which satisfies the properties (56) and (57). For the
sake of simplicity we report here the explicit derivation for the special case
in which $L$ is an exact divisor of $\mbox{r}[\rho]$ and ${\mathfrak{P}}$ is
an uniform partition of ${\mathfrak{K}}[\rho]$, so that Eq. (64) holds true:
the extension of this construction to the general case is given in Appendix C.
Introduce next an auxiliary Hilbert space ${\cal H}_{\bf B}$ of dimension $L$,
with an orthonormal basis $\left\\{\ket{v_{\ell}};\ell=1,\cdots,L\right\\}$
and replace the matrices $A^{ij}_{\bf{A}}$ of Eq. (III.4.1) with the following
set of $L^{2}\times L^{2}$ matrices:

$\displaystyle A^{ij}_{\bf{B}}$ $\displaystyle:=$
$\displaystyle\sqrt{\lambda_{i}\lambda_{j}}\left(\mathscr{S}_{{\ell^{[i]}}}\mathscr{S}_{{\ell^{[j]}}}\right)^{\frac{1-N}{2N}}\;\delta_{k^{[i]},k^{[j]}}$
(87)
$\displaystyle\times\ket{v_{\ell^{[i]}}}\bra{v_{\ell^{[i]}}}\otimes\ket{v_{\ell^{[j]}}}\bra{v_{\ell^{[j]}}}\;,$

where $k^{[i]}$, $\ell^{[i]}$ and $\mathscr{S}_{\ell}$ defined as in Eqs. (67)
and (68), respectively. We now define the block-wise purified density operator
$\rho^{(N)}_{{\bf B}}$ via the TI-MPO representation

$\displaystyle\rho^{(N)}_{{\bf B}}$ $\displaystyle:=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{\lambda_{i_{1}}\dots\lambda_{i_{N}}}\bra{\lambda_{j_{1}}\dots\lambda_{j_{N}}}$
(88) $\displaystyle\;\times\Tr[{A}^{i_{1},j_{1}}_{{\bf B}}\dots
A^{i_{N},j_{N}}_{{\bf B}}]\;.$

By construction we get that $\rho^{(N)}_{{\bf B}}$ has non-zero entries

$\displaystyle\bra{\lambda_{j_{1}}\lambda_{j_{2}}\dots\lambda_{j_{N}}}\rho^{(N)}_{{\bf
B}}\ket{\lambda_{i_{1}}\lambda_{i_{2}}\dots\lambda_{i_{N}}}$ (89)
$\displaystyle\qquad\qquad\qquad=\sqrt{\frac{\lambda_{i_{1}}\lambda_{i_{2}}\cdots\lambda_{i_{N}}}{\mathscr{S}_{\ell^{[i_{1}]}}^{N-1}}}\sqrt{\frac{\lambda_{j_{1}}\lambda_{j_{2}}\cdots\lambda_{j_{N}}}{\mathscr{S}_{\ell^{[j_{1}]}}^{N-1}}}\;,$

if the following conditions are true:

  * •

$\ell^{[i_{1}]}=\ell^{[i_{2}]}=\cdots=\ell^{[i_{N}]}$,

  * •

$\ell^{[j_{1}]}=\ell^{[j_{2}]}=\cdots=\ell^{[j_{N}]}$,

  * •

$\left(k^{[i_{1}]},k^{[i_{2}]},\cdots
k^{[i_{N}]}\right)=\left(k^{[j_{1}]},k^{[j_{2}]},\dots k^{[j_{N}]}\right)$;

while

$\bra{\lambda_{j_{1}}\lambda_{j_{2}}\dots\lambda_{j_{N}}}\rho^{(N)}_{{\bf
B}}\ket{\lambda_{i_{1}}\lambda_{i_{2}}\dots\lambda_{i_{N}}}=0\;,$ (90)

otherwise. Making use of Eq. (45) and the fact that in the present case we
have $\bar{A}_{\bf
B}=\sum_{\ell=1}^{L}\mathscr{S}_{{\ell}}^{\frac{1}{N}}\ket{v_{\ell}}\bra{v_{\ell}}\otimes\ket{v_{\ell}}\bra{v_{\ell}}$,
it is not difficult to verify that $\rho^{(N)}_{{\bf B}}$ satisfies the
partial trace condition (8): hence we can claim that $\rho^{(N)}_{{\bf B}}$ is
an element of $\mathfrak{S}_{T}^{(N,L^{2})}(\rho)$, i.e.

$\displaystyle\rho^{(N)}_{{\bf B}}\in\mathfrak{S}_{T}^{(N,L^{2})}(\rho)\;.$
(91)

A more explicit form for the density matrix $\rho^{(N)}_{{\bf B}}$ can be
obtained exploiting correspondence (67) and the fact that, thanks to the
choice of working with uniform partitions, the index $k$ of
${\lambda}_{k_{a},\ell}$ run from $1$ to $\mbox{r}[\rho]/L$ irrespectively
from the value of $\ell$. This leads to

$\rho^{(N)}_{{\bf
B}}=\sum_{\vec{k}}\sum_{\ell,\ell^{\prime}=1}^{L}\sqrt{\tfrac{{\lambda}^{(N)}_{\vec{k},\ell}}{\mathscr{S}_{\ell}^{N-1}}}\sqrt{\tfrac{{\lambda}^{(N)}_{\vec{k},\ell^{\prime}}}{\mathscr{S}_{\ell^{\prime}}^{N-1}}}\ket{{\lambda}^{(N)}_{\vec{k},\ell}}\bra{{\lambda}^{(N)}_{\vec{k},\ell^{\prime}}}\;,$
(92)

where the first sum runs over the $N$-uple $\vec{k}:=\left(k_{1},k_{2},\dots
k_{N}\right)$, and where we used the short-hand notation

$\displaystyle{\lambda}^{(N)}_{\vec{k},\ell}$ $\displaystyle:=$
$\displaystyle\prod_{a=1}^{N}{\lambda}_{k_{a},\ell}\;,$ (93)
$\displaystyle\ket{\lambda^{(N)}_{\vec{k},\ell}}$ $\displaystyle:=$
$\displaystyle\ket{\lambda_{k_{1},\ell}\;\lambda_{k_{2},\ell}\;\dots\;\lambda_{k_{N},\ell}}\;.$
(94)

Figure 4: Example of a block-wise purified state $\rho^{(N)}_{{\bf B}}$ of Eq.
(96) with $N=2$, $\mbox{r}[\rho]=r=12$, $L=3$, $M=4$. The colored dots are the
entries of the $r^{N}\times r^{N}$ density matrix which are different from
zero. The dots connected by lines constitute a bloc of rank one.

A close inspection of Eq. (92) reveals that the matrix $\rho^{(N)}_{{\bf B}}$
can be divided in a collection of uncoupled blocs, each one identified by a
value of the $N$-uple $\vec{k}$, each having rank one, and with non-zero
eigenvalue given by

$\Lambda^{(N)}_{\vec{k}}:=\sum_{\ell=1}^{L}\mathscr{S}_{\ell}\left(\frac{\lambda^{(N)}_{\vec{k},\ell}}{\mathscr{S}_{\ell}^{N}}\right)=\sum_{\ell=1}^{L}\frac{\lambda^{(N)}_{\vec{k},\ell}}{\mathscr{S}_{\ell}^{N-1}}\;,$
(95)

(see Fig. 4). Specifically we get

$\displaystyle\rho^{(N)}_{{\bf
B}}=\sum_{\vec{k}}\Lambda_{\vec{k}}^{(N)}\;|\Psi^{(N)}_{\vec{k}}\rangle\langle\Psi^{(N)}_{\vec{k}}|\;,$
(96)

with

$|\Psi^{(N)}_{\vec{k}}\rangle:=\tfrac{1}{\sqrt{\Lambda_{\vec{k}}^{(N)}}}\sum_{\ell=1}^{L}\sqrt{\tfrac{\lambda^{(N)}_{\vec{k},\ell}}{\mathscr{S}_{\ell}^{N-1}}}\;\ket{\lambda^{(N)}_{\vec{k},\ell}}\;,$
(97)

being orthonormal elements of $\mathfrak{S}^{(N)}$. It is worth stressing that
(96) is properly normalized thanks to the fact that

$\displaystyle\sum_{\vec{k}}\Lambda_{\vec{k}}^{(N)}$ $\displaystyle=$
$\displaystyle\sum_{\vec{k}}\sum_{\ell=1}^{L}\mathscr{S}_{\ell}^{1-N}\left(\prod_{a=1}^{N}\lambda_{k_{a},\ell}\right)$
$\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\mathscr{S}_{\ell}^{1-N}\left(\sum_{k=1}^{\\#{\mathfrak{K}}_{\ell}}\lambda_{k,\ell}\right)^{N}=\sum_{\ell=1}^{L}\mathscr{S}_{\ell}=1\;.$

At variance with the states (75) defined in the previous section, the
individual eigenvectors $|\Psi^{(N)}_{\vec{k}}\rangle$ of $\rho^{(N)}_{{\bf
B}}$ are in general not translationally invariant, while of course (96) obeys
to such symmetry. The rank of $\rho^{(N)}_{{\bf B}}$ is equal to the number of
$N$-uple $\vec{k}$: using the assumption (64) this leads to

$\displaystyle\mbox{rank}[\rho^{(N)}_{{\bf B}}]=({\mbox{r}[\rho]}/{L})^{N}\;,$
(99)

corresponding to a reduction of a factor $L^{-N}$ with respect to
$\rho^{\otimes N}$, which instead has rank $\mbox{r}[\rho]^{N}$. More
generally, as shown in Appendix C, when $L$ in not an exact divisor of
$\mbox{r}[\rho]$, Eq. (99) gets replaced by

$\displaystyle\mbox{rank}[\rho^{(N)}_{{\bf
B}}]=\left(\\#_{\max}\right)^{N}\;,$ (100)

with $\\#_{\max}$ the maximum cardinalities of the of the elements of
${\mathfrak{P}}$. In the extremal cases where $L=1$ and $L=\mbox{r}[\rho]$ we
get

$\displaystyle\rho^{(N)}_{{\bf B}}\Big{|}_{L=1}$ $\displaystyle=$
$\displaystyle\rho^{\otimes N}\;,$ (101) $\displaystyle\rho^{(N)}_{{\bf
B}}\Big{|}_{L=\mbox{r}[\rho]}$ $\displaystyle=$
$\displaystyle|\mbox{GHZ}^{(N)}_{\rho}\rangle\langle\mbox{GHZ}^{(N)}_{\rho}|\;,$
(102)

with the last expression implying that the present TI-MPO representation
assigns to $|\mbox{GHZ}^{(N)}_{\rho}\rangle$ the same BLR value as the the TI-
MPO representation of Sec. III.4.1 – indeed in both cases predict the state to
be an element of $\mathfrak{S}_{T}^{(N,\mbox{r}^{2}[\rho])}(\rho)$.

####  III.4.3 TI-MPO representation for classically correlated density
matrices

Our final example of TI-MPO quantum states is formed by a family of
classically correlated density matrices represented by the elements
$\rho^{(N)}_{{\bf C}}$ of Fig. 3, which connect the completely uncorrelated
$\rho^{\otimes N}$ to the classical correlated state $\rho^{(N)}_{\rm cc}$ of
Eq. (77). As summarized in Table 1 this family exhibits behaviours in terms of
rank and BLR which is almost complementary with respect to those of the family
$\rho^{(N)}_{{\bf A}}$ of Sec. III.4.1. They will be used to derive
Proposition 2.

Starting again from the partition ${\mathfrak{P}}$ of Eq. (55) we now consider
the following TI-MPO operator

$\displaystyle\rho^{(N)}_{{\bf C}}$ $\displaystyle:=$
$\displaystyle\sum_{i_{1},j_{1},\dots,i_{N},j_{N}}\ket{\lambda_{i_{1}}\dots\lambda_{i_{N}}}\bra{\lambda_{j_{1}}\dots\lambda_{j_{N}}}$
(103) $\displaystyle\;\times\Tr[{A}^{i_{1},j_{1}}_{{\bf C}}\dots
A^{i_{N},j_{N}}_{{\bf C}}]\;,$

with matrices

$A^{ij}_{{\bf
C}}:=\lambda_{i}\mathscr{S}_{\ell^{[i]}}^{\frac{1-N}{N}}\;\delta_{i,j}\;\ket{{v}_{\ell^{[i]}}}\bra{{v}_{\ell^{[i]}}}\;,$
(104)

operating on an auxiliary Hilbert space ${\cal H}_{\bf C}$ of dimension $L$
and characterized by an orthonormal basis
$\left\\{\ket{{v}_{\ell}};\ell=1,\cdots,L\right\\}$ (in the above expressions
$\ell^{[i]}$ and $\mathscr{S}_{\ell}$ are defined as in Eqs. (67) and (68),
respectively). By direct inspection one can verify that these states respect
the partial trace condition (8) and admit the following diagonal form

$\displaystyle\rho^{(N)}_{{\bf C}}$ $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\mathscr{S}_{\ell}^{1-N}\left(\sum_{k=1}^{\\#\mathfrak{K}_{\ell}}\lambda_{k,\ell}\ket{\lambda_{k,\ell}}\bra{\lambda_{k,\ell}}\right)^{\otimes
N}$ (105) $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\sum_{\vec{k}\in\\{1,\cdots,\\#\mathfrak{K}_{\ell}\\}^{N}}\Lambda_{\vec{k},\ell}^{(N)}|\lambda_{\vec{k},\ell}^{(N)}\rangle\langle\lambda_{\vec{k},\ell}^{(N)}|\;,$

with eigenvectors $|\lambda_{\vec{k},\ell}^{(N)}\rangle$ defined as in Eq.
(94), and associated eigenvalues given by

$\displaystyle\Lambda_{\vec{k},\ell}^{(N)}:=\mathscr{S}_{\ell}\left(\frac{\lambda^{(N)}_{\vec{k},\ell}}{\mathscr{S}_{\ell}^{N}}\right)=\mathscr{S}_{\ell}^{1-N}\lambda_{\vec{k},\ell}^{(N)}\;.$
(106)

Figure 5: Example of the state $\rho^{(N)}_{{\bf C}}$ defined in (105) with
$N=2$, $\mbox{r}[\rho]=r=12$, $L=4$, ${M}=3$. The colored dots are the entries
of the $r^{N}\times r^{N}$ density matrix which are different from zero.

We stress that the states $\rho^{(N)}_{{\bf C}}$ are explicitly separable with
respect to all possible partitions of the sites, and are diagonal in the same
basis of the tensor product state $\rho^{\otimes N}$. In particular, by
choosing $\rho$ to be diagonal in the energy eigenbasis of $H$, we can force
$\rho^{(N)}_{{\bf C}}$ to be diagonal in the eigenbasis of $H^{(N)}$. By
construction it also follows that the density matrix $\rho^{(N)}_{{\bf C}}$
has BLR value that is upper bounded by $L$, so that

$\displaystyle\rho^{(N)}_{{\bf C}}\in\mathfrak{S}_{T}^{(N,L)}(\rho)\;,$ (107)

(remember that for the case discussed in Sec. III.4.1, $L$ measured the rank
of the state $\rho^{(N)}_{{\bf A}}$). The rank of $\rho^{(N)}_{{\bf C}}$ is
instead equal to the total number of terms
$|\lambda_{\vec{k},\ell}^{(N)}\rangle$ entering (105), i.e.

$\displaystyle\mbox{rank}[\rho^{(N)}_{{\bf
C}}]=\sum_{\ell=1}^{L}\left(\\#\mathfrak{K}_{\ell}\right)^{N}=:|{\mathfrak{P}}|^{(N)}\;,$
(108)

where in the last identity we invoked Eq. (60) – for $\rho^{(N)}_{{\bf A}}$,
the term $|{\mathfrak{P}}|^{(2)}$ was instead an upper bound for the BLR. In
particular for uniform partitions Eq. (108) corresponds to have

$\displaystyle\mbox{rank}[\rho^{(N)}_{{\bf
C}}]\Big{|}_{{\mathfrak{P}}_{\mbox{\tiny{unif}}}}=L({\mbox{r}[\rho]}/{L})^{N}\;,$
(109)

with a reduction of a factor $L^{-N+1}$ with respect to the rank of the
completely uncorrelated state $\rho^{\otimes N}$.

We conclude by observing that in the extremal cases where $L=1$ and
$L=\mbox{r}[\rho]$ we get the tensor product and the the fully correlated
classical state respectively, i.e.

$\displaystyle\rho^{(N)}_{{\bf C}}\Big{|}_{L=1}=\rho^{\otimes
N}\;,\qquad\rho^{(N)}_{{\bf C}}\Big{|}_{L=\mbox{r}[\rho]}=\rho^{(N)}_{\rm
cc}\;.$ (110)

Notice finally that the TI-MPO representation associated with the family
$\rho^{(N)}_{{\bf C}}$ assigns to $\rho^{(N)}_{\rm cc}$ the same BLR value as
the the TI-MPO representation of the family $\rho^{(N)}_{{\bf A}}$ of Sec.
III.4.1 – indeed in both cases predict the state to be an element of
$\mathfrak{S}_{T}^{(N,\mbox{r}[\rho])}(\rho)$.

##  IV Derivations

This section is dedicated to derive the results anticipated in Sec. II.1.

###  IV.1 Proof of Propositions 1 and 2

Recalling the definition $\mathfrak{E}_{s}$ of the thermal energy associated
with the single-site entropy $s$ given in Eq. (189) and of the corresponding
effective heat capacity $\mathfrak{C}_{s}$, we start by proving the following
statement:

######  Proposition 4.

For any single-site entropy value $s_{0}>0$ there exist
$\eta_{\star}\in\;]0,1[$ and a function $N_{\star}(\eta)$ mapping
$\eta\in\;]0,\eta_{\star}[$ into $\mathbb{R}^{+}$, such that given
$\eta\in\;]0,\eta_{\star}[$ and $N\geq N_{\star}(\eta)$, any state
$\rho^{(N)}\in\mathfrak{S}_{T}^{(N)}$ admitting a eigenvalue subset
${\mathfrak{T}}^{(N)}$ that has cardinality and associated total population
satisfying the constraints

$\displaystyle\\#{\mathfrak{T}}^{(N)}$ $\displaystyle\leq$ $\displaystyle
e^{N\left(s_{0}+\eta\right)}\;,$ (111)
$\displaystyle\sum_{\lambda_{i}^{(N)}\in{\mathfrak{T}}^{(N)}}\lambda_{i}^{(N)}$
$\displaystyle\geq$ $\displaystyle 1-e^{-\frac{2N\eta^{2}}{\alpha^{2}}}\;,$
(112)

for some given $\alpha>0$, fulfils the following inequality

$\displaystyle\frac{{\cal E}(\rho^{(N)};H^{(N)})}{N}$ (113)
$\displaystyle\quad\geq
E(\rho;H)-\mathfrak{E}_{s_{0}}-(2\sqrt{2.01\mathfrak{C}_{s_{0}}\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}}})\;,$

where $\epsilon_{d}$ is the largest eigenvalue of the single-site Hamiltonian
$H$.

Proof:– Given $s_{0}>0$ and $\alpha>0$, assume that $\rho^{(N)}$ is an element
of $\mathfrak{S}_{T}^{(N)}$ which admits an subset ${\mathfrak{T}}^{(N)}$ of
eigenvalues fulfilling the conditions (111) and (112) for some $\eta>0$ whose
value will be determined in the following. Since the ergotropy, as defined in
Eq. (157), is a maximum over the set of unitary transformation, to prove
statement of the proposition it would be sufficient to show an unitary
transformation $U_{\star}\in{\bf U}({\cal H}^{\otimes N})$ such that the
associated mean energy difference

$\displaystyle
W_{U_{\star}}(\rho^{(N)};H^{(N)}):=E(\rho^{(N)};H^{(N)})-E(U_{\star}\rho^{(N)}U^{\dagger}_{\star};H^{(N)})$

is greater than the right-hand side of (113). As suitable candidates for
$U_{\star}$ we consider a unitary transformation which exchanges the
eigenvectors associated with the elements of the set ${\mathfrak{T}}^{(N)}$
with a subset of the eigenvectors corresponding to the elements of the sets
${\mathfrak{A}_{\xi}^{(N)}}$ of Lemma 2 associated to the same value of
$s_{0}$ and to a proper choice of $\xi$. This construction might be regarded
as a more general case of the one proposed in Perarnau-Llobet and Uzdin
(2019), where it is applied for reducing the fluctuations in the extracted
work - while we are interested in the mean values
$W_{U_{*}}\left(\rho^{\otimes N};H^{(N)}\right)$. The transformation we are
targeting clearly can be identified if in ${\mathfrak{T}}^{(N)}$ there are
more elements than in ${\mathfrak{A}^{(N)}}$, i.e.

$\\#{\mathfrak{A}_{\xi}^{(N)}}\geq\\#{\mathfrak{T}}^{(N)}\;.$ (115)

The hypothesis (112) gives us the upper bound for the cardinality of
${\mathfrak{T}}^{(N)}$, while Lemma 2 provides the lower bound (201) for the
cardinality of ${\mathfrak{A}_{\xi}^{(N)}}$ providing that the selected $\xi$
and $N$ fulfils the constraints

$\displaystyle\xi\leq\xi_{*}\;,\qquad N\geq N_{*}(\xi)\;,$ (116)

with $\xi_{*}$ and $N_{*}(\xi)$ dependent upon $s_{0}$ and the structure of
single-site Hamiltonian – see Remark 4 below Lemma 2 in Appendix B.1.
Combining these facts, we can see that the requirement (115) can be satisfied
by imposing

$\eta=\frac{\xi^{2}}{2.01\mathfrak{C}_{s_{0}}}\;,$ (117)

which we write incorporating the assumption (216). Notice now that with this
choice all the eigenvalues of $\rho^{(N)}$ belonging to ${\mathfrak{T}}^{(N)}$
will be associated with energy eigenvalues $\epsilon_{\vec{i}}$ that can be
upper bounded via the inequality (200); the remaining one instead can be
bounded by the maximum eigenvalue of $H^{(N)}$, i.e. $N\epsilon_{d}$.
Accordingly we can write

$\displaystyle E(U_{\star}\rho^{(N)}U_{\star}^{\dagger};H^{(N)})$
$\displaystyle\quad\leq
N(\mathfrak{E}_{s_{0}}+2\xi)\left(\sum_{\lambda_{i}^{(N)}\in{\mathfrak{T}}^{(N)}}\lambda_{i}\right)+N\epsilon_{d}\left(\sum_{\lambda_{i}^{(N)}\notin{\mathfrak{T}}^{(N)}}\lambda_{i}^{(N)}\right)$
$\displaystyle\quad\leq
N(\mathfrak{E}_{s_{0}}+2\xi)+N\epsilon_{d}\left(1-\sum_{\lambda_{i}^{(N)}\in{\mathfrak{T}}^{(N)}}\lambda_{i}\right)$
$\displaystyle\quad\leq
N(\mathfrak{E}_{s_{0}}+2\xi)+N\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}}}$
$\displaystyle\quad=N(\mathfrak{E}_{s_{0}}+2\sqrt{2.01\mathfrak{C}_{s_{0}}\eta})+N\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}}}\;,$

where in the last passage we used Eq. (112). From this we can then obtain the
thesis

$\displaystyle\frac{{\cal E}(\rho^{(N)};H^{(N)})}{N}\geq
W_{U_{\star}}(\rho^{(N)};H^{(N)})$ (118) $\displaystyle\quad\geq
E(\rho;H)-\mathfrak{E}_{s_{0}}-(2.01\sqrt{2\mathfrak{C}_{s_{0}}\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}}})\;,$

which according to (117) and (116) is valid for $\eta>0$ satisfying the
constraint

$\displaystyle\eta\leq\eta_{\star}:=\min\left\\{\frac{\xi_{*}^{2}}{2.01\mathfrak{C}_{s_{0}}},1\right\\}\;,$
(119)

and for $N$ integers such that

$\displaystyle N\geq
N_{\star}(\eta):=N_{*}(\sqrt{2.01\mathfrak{C}_{s_{0}}\eta})\;.$ (120)

Notice in particular that thanks to (215) we can fix the functional dependence
of $N_{\star}(\eta)$ as

$\displaystyle N_{\star}(\eta)={K_{\star}}/{\eta}\;,$ (121)

with the constant term

$\displaystyle K_{\star}:=\frac{K_{*}}{2.01\mathfrak{C}_{s_{0}}}\;,$ (122)

depending upon $s_{0}$ and $H$. $\square$

Remark 1: We stress that the coefficient $2.01$ appearing on the right-hand-
side of Eq. (113) is a byproduct of the choice (216): generalization for
arbitrary $\zeta\in]0,1[$ can be obtained by simply replacing $2.01$ with
$2/\zeta$.

Remark 2: The thermal capacity $\mathfrak{C}_{s_{0}}$ is always finite. In
fact, as shown in Ref.Correa  _et al._ (2015), it can be bounded by a constant
which depends only on the dimension $d$ of the Hilbert space. Therefore we can
always rewrite (113) in the weaker form

$\displaystyle\frac{{\cal E}(\rho^{(N)};H^{(N)})}{N}$ (123)
$\displaystyle\quad\geq
E(\rho;H)-\mathfrak{E}_{s_{0}}-(2\sqrt{2.01\mathfrak{C}_{\max}\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}}})\;,$

with the definition

$\displaystyle\mathfrak{C}_{\max}:=\max_{0<s<\ln\mbox{r}[\rho]}\mathfrak{C}_{s}\;.$
(124)

####  IV.1.1 Proof of Proposition 1

The result follows by showing that family of density matrix $\rho^{(N)}_{{\bf
B}}$ introduced in Sec. III.4.2 and Appendix C fulfils the hypotheses of
Proposition 4. In particular we shall focus on those cases where the partition
${\mathfrak{P}}$ of Eq. (55) is almost uniform, so that Eq. (65) holds true.
Under this condition for every $\ell\in\\{1,\cdots,L\\}$, let us introduce the
quantities

$s[\mathfrak{K}_{\ell}]:=-\sum_{k=1}^{\\#\mathfrak{K}_{\ell}}\frac{\lambda_{k,\ell}}{\mathscr{S}_{\ell}}\ln{\frac{\lambda_{k,\ell}}{\mathscr{S}_{\ell}}}\,,$
(125)

which represents the Shannon entropy of the probability distribution
$\\{\lambda_{k,\ell}/\mathscr{S}_{\ell}\\}_{k=1,\dots,\\#\mathfrak{K}_{\ell}}$.
Since each $s[\mathfrak{K}_{\ell}]$ is the sum of $\\#\mathfrak{K}_{\ell}$
terms, from (58) and (65) we have the trivial bound

$s[\mathfrak{K}_{\ell}]\leq\ln(\\#_{\max})=\ln\left\lceil{\mbox{r}[\rho]/L}\right\rceil\;.$
(126)

Identifying then the probability distribution $\mathfrak{X}$ of Lemma 1 with
$\\{\lambda_{k,\ell}/\mathscr{S}_{\ell}\\}_{k=1,\dots,\\#\mathfrak{K}_{\ell}}$,
it follows that the set

${\mathfrak{X}}_{\ell}^{(N)}:=\left\\{\vec{k}\;\middle|\;\mathscr{S}^{-N}_{\ell}\prod_{a=1}^{N}\lambda_{k_{a},\ell}\geq
e^{-N(s[\mathfrak{K}_{\ell}]+\eta)}\right\\}\;,$ (127)

has cardinality bounded by

$\\#{\mathfrak{X}}_{\ell}^{(N)}\leq e^{N(s[\mathfrak{K}_{\ell}]+\eta)}\;,$
(128)

and satisfy the inequality

$\sum_{\vec{k}\in{\mathfrak{X}}_{\ell}^{(N)}}\prod_{a=1}^{N}\lambda_{k_{a},\ell}\geq\mathscr{S}^{N}_{\ell}\left(1-e^{-\frac{2N\eta^{2}}{\alpha^{2}[\mathfrak{K}_{\ell}]}}\right)\;,$
(129)

with

$\alpha[\mathfrak{K}_{\ell}]:=\max_{k=1,\cdots,\\#\mathfrak{K}_{\ell}}\ln\lambda_{k,\ell}-\min_{k=1,\cdots,\\#\mathfrak{K}_{\ell}}\ln\lambda_{k,\ell}\;.$
(130)

Notice that (130) is smaller than or equal to the logarithmic spectral ratio
$\alpha(\rho)$ of the single-site density matrix, i.e. the quantity

$\displaystyle\alpha(\rho):=\ln\left(\lambda_{\max}/\lambda_{\min}\right)\;,$
(131)

with $\lambda_{\max}$ and $\lambda_{\min}$ being respectively the maximum and
the minimum positive eigenvalues of $\rho$. We can hence replace (130) with
the inequality

$\sum_{\vec{k}\in{\mathfrak{X}}_{\ell}^{(N)}}\prod_{a=1}^{N}\lambda_{k_{a},\ell}\geq\mathscr{S}^{N}_{\ell}\left(1-e^{-\frac{2N\eta^{2}}{\alpha^{2}(\rho)}}\right)\;.$
(132)

Observe next that the union of the all the sets ${\mathfrak{X}}_{\ell}^{(N)}$,

${\mathfrak{X}}^{(N)}:=\bigcup_{\ell=1}^{L}{\mathfrak{X}}_{\ell}^{(N)}\;,$
(133)

has cardinality bounded by

$\\#{\mathfrak{X}}^{(N)}\leq\sum_{\ell=1}^{L}\\#{\mathfrak{X}}^{(N)}_{\ell}\leq\sum_{\ell=1}^{L}e^{N(s[\mathfrak{K}_{\ell}]+\eta)}\leq
Le^{N(\ln\left\lceil{\mbox{r}[\rho]/L}\right\rceil+\eta)}\;,$ (134)

where in the third inequality we used (128) and in the last passage we invoked
(126).

Let us recall that the eigenvalues of $\rho^{(N)}_{{\bf B}}$ are the
quantities $\Lambda^{(N)}_{\vec{k}}$ reported in Eq. (95) and consider the
subset ${\mathfrak{T}}_{\bf B}^{(N)}$ of such values characterized by
$\vec{k}$ vectors belonging to ${\mathfrak{X}}^{(N)}$, i.e.

$\displaystyle{\mathfrak{T}}^{(N)}_{\bf
B}:=\left\\{\Lambda^{(N)}_{\vec{k}}\big{|}\vec{k}\in{\mathfrak{X}}^{(N)}\right\\}\;.$
(135)

On one hand, using Eq. (134) we can hence write

$\\#{\mathfrak{T}}^{(N)}_{\bf B}=\\#{\mathfrak{X}}^{(N)}\leq
e^{N(\ln\left\lceil{\mbox{r}[\rho]/L}\right\rceil+N^{-1}\ln L+\eta)}\;,$ (136)

while, on the other hand, using (133) and (132) we can infer that the total
population associated with such subset is at least

$\displaystyle\sum_{\Lambda^{(N)}_{\vec{k}}\in{\mathfrak{T}}^{(N)}_{\bf
B}}\Lambda^{(N)}_{\vec{k}}$ $\displaystyle=$
$\displaystyle\sum_{\vec{k}\in{\mathfrak{X}}^{(N)}}\sum_{\ell=1}^{L}\mathscr{S}_{\ell}^{1-N}\prod_{a=1}^{N}\lambda_{k_{a},\ell}$
$\displaystyle\geq$
$\displaystyle\sum_{\ell=1}^{L}\sum_{\vec{k}\in{\mathfrak{X}}_{\ell}^{(N)}}\mathscr{S}_{\ell}^{1-N}\prod_{a=1}^{N}\lambda_{k_{a},\ell}$
$\displaystyle\geq$
$\displaystyle\left(1-e^{-\frac{2N\eta^{2}}{\alpha^{2}(\rho)}}\right)\sum_{\ell=1}^{L}\mathscr{S}_{\ell}=1-e^{-\frac{2N\eta^{2}}{\alpha^{2}(\rho)}}\;.$

Equations (134) and (IV.1.1) certify that the set ${\mathfrak{T}}_{\bf
B}^{(N)}$ fulfils the hypotheses (112) and (111) of Proposition 4, with

$\displaystyle s_{0}=\ln\left\lceil{\mbox{r}[\rho]/L}\right\rceil+N^{-1}\ln
L\;.$ (138)

Therefore we can conclude there exist $\eta_{\star}\in\;]0,1[$ and a function
$N_{\star}(\eta)$ mapping $\eta\in\;]0,\eta_{\star}[$ into $\mathbb{R}^{+}$,
such that given $\eta\in\;]0,\eta_{\star}[$ and $N\geq N_{\star}(\eta)$, we
have

$\displaystyle\frac{{\cal E}(\rho^{(N)}_{{\bf B}};H^{(N)})}{N}$ (139)
$\displaystyle\quad\geq E(\rho;H)-\mathfrak{E}_{s_{\bf
B}}-(2\sqrt{2.01\mathfrak{C}_{\max}\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}(\rho)}})\;,$

where we used the weaker version (123) of (113) discussed in Remark 2, and
recall Eq. (91) to set $m=L^{2}$ to identify $s_{0}$ of Eq. (138) with $s_{\bf
B}$ of Eq. (24). Proposition 1 now finally follows by identifying the
constants $C$ and $\alpha$ of Eq. (23) with $\mathfrak{C}_{\max}$ and
$\alpha(\rho)$ respectively, and observing that for each $m$ and $\eta$, by
choosing $N$ sufficiently large Eq. (139) is a trivial lower bound for the
$\mathcal{E}^{(N,m)}_{\max}(\rho;H)$. $\square$

####  IV.1.2 Proof of Proposition 2

The proof exploits again Proposition 4 and closely mimics the one we presented
for Proposition 1, the only difference being that this time we replace the
states $\rho^{(N)}_{{\bf B}}$ with the family of classically correlated states
$\rho^{(N)}_{{\bf C}}$ defined in Sec. III.4.2 (orange trajectory of Fig. 3)
under the hypothesis that it is generated by an almost uniform partition
${\mathfrak{P}}$ (see Appendix C), so that (65) is true.

To prove Proposition 2, we need to show that the state $\rho^{(N)}_{{\bf C}}$
has a subset of eigenvalues ${\mathfrak{T}}_{\bf C}^{(N)}$ which satisfies
(112) and (111). Observe that in this case the eigenvalues of
$\rho^{(N)}_{{\bf C}}$ are given by the expressions
$\Lambda_{\vec{k},\ell}^{(N)}$ of Eq. (106) which are labelled by $\ell$ and
$\vec{k}$ (instead the eigenvalues of $\rho^{(N)}_{{\bf B}}$ where identified
only by the vectors $\vec{k}$). Accordingly we define ${\mathfrak{T}}_{\bf
C}^{(N)}$ as

${\mathfrak{T}}_{\bf C}^{(N)}:=\bigcup_{\ell=1}^{L}{\mathfrak{T}}_{{\bf
C},\ell}^{(N)}\;,$ (140)

where for $\ell\in\\{1,\cdots,L\\}$ we take

$\displaystyle{\mathfrak{T}}^{(N)}_{{\bf
C},\ell}:=\left\\{\Lambda^{(N)}_{\vec{k},\ell}\big{|}\vec{k}\in{\mathfrak{X}}_{\ell}^{(N)}\right\\}\;,$
(141)

with ${\mathfrak{X}}_{\ell}^{(N)}$ defined as in Eq. (127). Accordingly we
have

$\displaystyle\\#{\mathfrak{T}}^{(N)}_{\bf
C}=\sum_{\ell=1}^{L}\\#{\mathfrak{T}}^{(N)}_{{\bf C},\ell}$ $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\\#{\mathfrak{X}}_{\ell}^{(N)}$
$\displaystyle\leq$ $\displaystyle
e^{N(\ln\left\lceil{\mbox{r}[\rho]/L}\right\rceil+N^{-1}\ln L+\eta)}\;,$

where the first identity follows from the fact that the sets
${\mathfrak{T}}^{(N)}_{{\bf C},\ell}$ are disjoint, and where in the last
inequality we invoked (134). Furthermore we observe

$\displaystyle\sum_{\Lambda^{(N)}_{\vec{k},\ell}\in{\mathfrak{T}}^{(N)}_{\bf
C}}\Lambda^{(N)}_{\vec{k},\ell}$ $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\sum_{\Lambda^{(N)}_{\vec{k},\ell}\in{\mathfrak{T}}^{(N)}_{{\bf
C},\ell}}\Lambda^{(N)}_{\vec{k},\ell}$ (143) $\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\sum_{{\vec{k}}\in{\mathfrak{X}}^{(N)}_{\ell}}\Lambda^{(N)}_{\vec{k},\ell}$
$\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\sum_{{\vec{k}}\in{\mathfrak{X}}^{(N)}_{\ell}}\mathscr{S}_{\ell}^{1-N}\lambda_{\vec{k},\ell}^{(N)}$
$\displaystyle\geq$ $\displaystyle
1-e^{-\frac{2N\eta^{2}}{\alpha^{2}(\rho)}}\;,$

where the last passage follows directly from (IV.1.1). We can therefore apply
Proposition 4 with $s_{0}$ as in Eq. (138) obtaining that there exist
$\eta_{\star}\in\;]0,1[$ and a function $N_{\star}(\eta)$ mapping
$\eta\in\;]0,\eta_{\star}[$ into $\mathbb{R}^{+}$, such that given
$\eta\in\;]0,\eta_{\star}[$ and $N\geq N_{\star}(\eta)$, we have

$\displaystyle\frac{{\cal E}(\rho^{(N)}_{{\bf C}};H^{(N)})}{N}$ (144)
$\displaystyle\quad\geq E(\rho;H)-\mathfrak{E}_{s_{\bf
C}}-(2\sqrt{2.01\mathfrak{C}_{\max}\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}(\rho)}})\;,$

where we recalled (107) to set $m=L$ and transform $s_{0}$ into $s_{\bf C}$ of
Eq. (28). Proposition 2 now finally follows by identifying again the constants
$C$ and $\alpha$ with $\mathfrak{C}_{\max}$ and $\alpha(\rho)$ respectively,
and by observing that for each $m$ and $\eta$, by choosing $N$ sufficiently
large Eq. (144) is a trivial lower bound for the
$\mathcal{E}^{(N,m)}_{\max}(\rho;H)$. $\square$

Remark 3: We now present the heuristic argument in support of the fact that,
as mentioned in Sec. II.1, it is reasonable to think that the functional
dependency of $s_{\bf C}$ of Eq. (28) of Proposition 2 can be improved for
small levels of correlation. To see this observe that the following identity
holds true

$\displaystyle\sum_{\ell=1}^{L}\mathscr{S}_{\ell}\;s[\mathfrak{K}_{\ell}]$
$\displaystyle=$
$\displaystyle-\sum_{\ell=1}^{L}\sum_{k=1}^{\\#\mathfrak{K}_{\ell}}{\lambda_{k,\ell}}\ln{\frac{\lambda_{k,\ell}}{\mathscr{S}_{\ell}}}$
(145) $\displaystyle=$
$\displaystyle-\sum_{\ell=1}^{L}\sum_{k=1}^{\\#\mathfrak{K}_{\ell}}{\lambda_{k,\ell}}\ln{{\lambda_{k,\ell}}}+\sum_{\ell=1}^{L}\sum_{k=1}^{\\#\mathfrak{K}_{\ell}}\lambda_{k,\ell}\ln{\mathscr{S}_{\ell}}$
$\displaystyle=$ $\displaystyle
S(\rho)+\sum_{\ell=1}^{L}\mathscr{S}_{\ell}\ln{\mathscr{S}_{\ell}}\;.$

Because of (69), equation (145) can be see as a weighted mean of the entropic
quantities $s[\mathfrak{K}_{\ell}]$ when they are averaged with weights
$\mathscr{S}_{\ell}$. Therefore, letting

$s_{\max}:=\max_{\ell}s[\mathfrak{K}_{\ell}]\;,$ (146)

it is always true that

$s_{\max}\geq
S(\rho)+\sum_{\ell=1}^{L}\mathscr{S}_{\ell}\ln{\mathscr{S}_{\ell}}\geq
S(\rho)-\ln L\;,$ (147)

with the inequalities (147) becoming equalities in the case in which all the
partial sums $\mathscr{S}_{\ell}$ were equal, i.e. if

$\mathscr{S}_{\ell}={1}/{L}\;,\qquad\forall\ell\in\\{1,\cdots,L\\}\;.$ (148)

The values of $\mathscr{S}_{\ell}$ are determined by the spectrum
$\\{\lambda_{i};i=1,\cdots,\mbox{r}[\rho]\\}$ of the state $\rho$ and by our
choice of the partition ${\mathfrak{P}}$ of Eq. (55). When
$L\ll\mbox{r}[\rho]$, and if the eigenvalues of $\rho$ are sufficiently evenly
distributed, we expect it to be possible to group them in $L$ subsets such
that the total populations $\mathscr{S}_{\ell}$ in each subset are
approximately equal in order to fulfil (148) with good approximation.
Accordingly we expect that in such conditions we can write

$s_{\max}\simeq S(\rho)-\ln L\;.$ (149)

Thanks to this we can now replace (126), with

$s[\mathfrak{K}_{\ell}]\leq s_{\max}\simeq S(\rho)-\ln L\;,$ (150)

and hence (134) with

$\displaystyle\\#{\mathfrak{X}}^{(N)}\leq\sum_{\ell=1}^{L}\\#{\mathfrak{X}}^{(N)}_{\ell}$
$\displaystyle\leq$
$\displaystyle\sum_{\ell=1}^{L}e^{N(s[\mathfrak{K}_{\ell}]+\eta)}$ (151)
$\displaystyle\lesssim$ $\displaystyle Le^{N(S(\rho)-\ln L+\eta)}\;,$

and (IV.1.2) with

$\displaystyle\\#{\mathfrak{T}}^{(N)}_{\bf C}\lesssim
e^{N(S(\rho)-\frac{N-1}{N}\ln L+\eta)}\;,$ (152)

Following the final passages of the proof we thus arrive to the conclusion
that Eq. (144) holds with $s_{\bf C}$ replaced by the term $s_{\bf
C}|_{\rm(heu)}$ of Eq. (35).

###  IV.2 Proof of Proposition 3

To prove Proposition 3 we focus on the family of states $\rho^{(N)}_{\bf A}$
introduced in Sec. III.4.1.

From Eq. (63) it follows that the quantity $|{\mathfrak{P}}|^{(2)}$ fulfils
the inequality

$\displaystyle\mbox{r}[\rho]\leq$ $\displaystyle|{\mathfrak{P}}|^{(2)}$
$\displaystyle\leq\mbox{r}^{2}[\rho]\;,$ (153)

with the lower and upper value being attained respectively by fixing the
number $L$ of elements of the partition choosing the partition
${\mathfrak{P}}$ equal to $\mbox{r}[\rho]$ and $1$. Therefore for each given

$\displaystyle m\in\\{\mbox{r}[\rho],\cdots,\mbox{r}^{2}[\rho]\\}\;,$ (154)

we can identify a special partition ${\mathfrak{P}}$ (55), such that

$\displaystyle|{\mathfrak{P}}|^{(2)}\leq m\;.$ (155)

Invoking then Eq. (79), we can claim that given $m$ as in Eq. (154) the
associated the density matrix $\rho^{(N)}_{\bf A}$ of Eq. (75) is an element
of $\mathfrak{S}_{T}^{(N,m)}(\rho)$, so that the following lower bound holds

$\displaystyle\mathcal{E}^{(N,m)}_{\max}(\rho;H)$ $\displaystyle\geq$
$\displaystyle\frac{\mathcal{E}(\rho^{(N)}_{\bf A};H^{(N)})}{N}$
$\displaystyle=$ $\displaystyle
E(\rho;H)-\frac{\min_{U\in{\bf{U}}(d^{N})}E(U\rho^{(N)}_{\bf
A}U^{\dagger};H^{(N)})}{N}$

see Eq. (16) and (157). We remind that the minimization on the right-hand-side
of (157) can be explicitly performed producing a closed expression in terms of
the spectra of $H^{(N)}$ and $\rho_{\bf A}^{(N)}$ Mirsky (1975); Allahverdyan
_et al._ (2004). This yields to simplified formula

$\min_{U\in{\bf{U}}(d^{N})}E(U\rho^{(N)}_{\bf
A}U^{\dagger};H^{(N)})=\sum_{j=1}^{d^{N}}\epsilon_{j}^{(N)}\lambda^{(N,\downarrow)}_{j}\;,$
(157)

where $\epsilon_{j}^{(N)}$ are the eigenvalues of $H^{(N)}$ that, as in the
case of $H$ we organize in increasing order, i.e.

$\displaystyle\epsilon_{j+1}^{(N)}\geq\epsilon_{j}^{(N)}\;,$ (158)

while $\lambda_{j}^{(N,\downarrow)}$ are the eigenvalues of $\rho_{\bf
A}^{(N)}$, which instead, as indicated by the arrow, are assumed to arranged
in decreasing order, i.e.

$\displaystyle\lambda_{j+1}^{(N,\downarrow)}\geq\lambda_{j}^{(N,\downarrow)}\;.$
(159)

A closed look at Eq. (75) reveals that these last quantities can be written as

$\displaystyle\lambda_{j}^{(N,\downarrow)}=\begin{dcases}\mathscr{S}_{j}^{(\downarrow)}&\mbox{
if }j\leq L\;;\\\ 0&\mbox{otherwise}\;,\end{dcases}$ (160)

where $\mathscr{S}_{j}^{(\downarrow)}$ are the partial sums (68) rearranged in
decreasing order, i.e.
$\mathscr{S}_{1}^{(\downarrow)}\geq\mathscr{S}_{2}^{(\downarrow)}\geq\dots\geq\mathscr{S}_{L}^{(\downarrow)}$.
Therefore we can write

$\displaystyle\min_{U\in{\bf{U}}(d^{N})}E(U\rho^{(N)}_{\bf
A}U^{\dagger};H^{(N)})$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{L}\epsilon_{j}^{(N)}\mathscr{S}_{j}^{(\downarrow)}\;.$
(161)

Now we notice that the structure of the eigenvalues of the Hamiltonian
$H^{(N)}$, given by (197), warrants that

$\displaystyle\epsilon_{j}^{(N)}\leq\epsilon_{j}\qquad\forall j\in[1,d]\;.$
(162)

Since by construction $L\leq d$, the inequality (162) is also true for
$j\in[1,L]$, and we can use it in (161) obtaining

$\displaystyle\min_{U\in{\bf{U}}(d^{N})}E(U\rho^{(N)}_{\bf
A}U^{\dagger};H^{(N)})\leq\sum_{j=1}^{L}\epsilon_{j}\mathscr{S}_{j}^{(\downarrow)}\;.$
(163)

Since the eigenvalues $\\{\mathscr{S}_{j}\\}_{j=1\dots L}$ are, by definition
(68), partial sums of the eigenvalues $\\{\lambda_{j}\\}_{j=1\dots d}$ of the
single-site density matrix $\rho$, we can write (see the appendix B of
Alimuddin  _et al._ (2020)):

$\min_{U\in{\bf{U}}(d^{N})}E(U\rho^{(N)}_{\bf
A}U^{\dagger};H^{(N)})\leq\sum_{j=1}^{L}\epsilon_{j}\mathscr{S}_{j}^{(\downarrow)}\leq\sum_{j=1}^{d}\epsilon_{j}\lambda_{j}^{(\downarrow)}\;.$
(164)

Noticing finally that

$\displaystyle\sum_{j=1}^{d}\epsilon_{j}\lambda_{j}^{(\downarrow)}$
$\displaystyle=$ $\displaystyle\min_{U\in{\bf{U}}(d)}E(U\rho U^{\dagger};H)$
(165) $\displaystyle=$ $\displaystyle{E}(\rho,H)-{\cal E}(\rho,H)\;,$

we finally arrive at (31) replacing (164) into (IV.2). $\square$

###  IV.3 Proof of Corollary 1

In view of the identity (18), Corollary 1 can be seen as a refinement of
Proposition 1 for $m=1$. Indeed we can derive the statement following the same
passages of Sec. IV.1.1 and observing that $\rho^{\otimes N}$ is the unique
element of family $\rho_{\bf B}^{(N)}$ we get when setting $L=1$ (see Eq.
(102)). In this case Eqs. (125) and (126) get replaced by
$s[\mathfrak{K}_{1}]=S(\rho)$ which in turn allow us to replace Eqs. (134) and
(136) with

$\\#{\mathfrak{T}}^{(N)}_{\bf B}=\\#{\mathfrak{X}}^{(N)}\leq
e^{N(S(\rho)+\eta)}\;.$ (166)

Invoking hence (IV.1.1) that still remains valid, we can conclude that now the
set ${\mathfrak{T}}_{\bf B}^{(N)}$ fulfils the hypotheses (112) and (111) of
Proposition 4, with $s_{0}=S(\rho)$ hence leading to

$\displaystyle\frac{{\cal E}(\rho^{\otimes N};H^{(N)})}{N}$ (167)
$\displaystyle\quad\geq
E(\rho;H)-\mathfrak{E}_{S(\rho)}-(2\sqrt{2.01\mathfrak{C}_{{S(\rho)}}\eta}+\epsilon_{d}e^{-\frac{2N\eta^{2}}{\alpha^{2}(\rho)}})\;,$

that corresponds to (37) by identifying $C$ and $\alpha$ with
$\mathfrak{C}_{{S(\rho)}}$ and $\alpha(\rho)$ respectively. $\square$

##  V Conclusions

We derived some analytic lower bound for the work that, in the best case, can
be extracted with unitary transformations from a translationally invariant
correlated many-body system, using as a measure of correlation the minimum
bond link rank (BLR) necessary to represent the state as a matrix product
operator. When the number $N$ of copies of the system is finite, non-classical
correlations are required to extract as work the full energy of the system.
However, in the macroscopic limit $N\to\infty$, we found that this quantum
feature disappears, and that it is possible to create many-body states with
classically correlated state which have a relative low BRL (equal at most to
to the rank $r$ of the local state, out of a maximum of $r^{2}$). Our bounds
do not depend on the entropy $S(\rho)$ on the state, so they are worse for
states of low entropy. However, heuristic consideration suggest that, at least
for small correlations strengths, the bounds can be improved with an explicit
dependence on $S(\rho)$.

We conjecture that, for $N\geq 3$, the BLR of the family of states employed in
our analysis is equal to the upper bounds that we found by explicit
construction. If true, this could allow to derive also upper bound for the
ergotropy of a translationally invariant state with a given correlation
strength. Another possible improvement of our work could be repeating the
analysis with a measure of correlation more sophisticated the TI-MPO bond link
rank, like some form of correlation entropy Schindler  _et al._ (2020);
Rolandi and Wilming (2020).

###### Acknowledgements.

We would like to thank Giacomo De Palma for comments and discussions. This
work is supported by MIUR (Ministero dell’Istruzione, dell’Università e della
Ricerca) via project PRIN 2017 Taming complexity via Quantum Strategies a
Hybrid Integrated Photonic approach (QUSHIP) Id. 2017SRNBRK.

##  Appendix A Chernoff inequality

If we extract $N$ times a random variable $X\in[a,b]$, the Chernoff bound
Chernoff (1952) (or equivalently Hoeffding’s inequality Hoeffding (1963)) tell
us that

$\mathbb{P}\left(\sum_{i=1}^{N}X_{i}\leq
N(\mathbb{E}[\mathcal{V}]+\varepsilon)\right)\geq
1-e^{-\frac{2N\varepsilon^{2}}{(b-a)^{2}}}\;.$ (168)

In this paper we will apply the following useful consequence of the Chernoff
inequality:

######  Lemma 1.

Let $\mathfrak{X}=\\{x_{i}\\}_{i}$ be a finite collection of positive real
numbers $x_{i}\in]0,1]$, such that $\sum_{i}x_{i}=1$, and $-\sum_{i}x_{i}\ln
x_{i}=s_{0}$. Let $\mathfrak{X}^{\otimes N}$ denote the set of $N$-ples
$\vec{x}:=(x_{i_{1}},\dots x_{i_{N}})$. Then, for any $\eta>0$, the subset
${\mathfrak{X}}^{(N)}\subseteq\mathfrak{X}^{\otimes N}$ defined by

${\mathfrak{X}}^{(N)}:=\left\\{\vec{x}=(x_{i_{1}},\dots,x_{i_{N}})\middle|\prod_{k=1}^{N}x_{i_{k}}\geq
e^{-N(s_{0}+\eta)}\right\\}\;,$ (169)

has cardinality

$\displaystyle\\#{\mathfrak{X}}^{(N)}\leq e^{N(s_{0}+\eta)}\;,$ (170)

and satisfies the property

$\displaystyle\sum_{\vec{x}\in{\mathfrak{X}^{(N)}}}\prod_{k=1}^{N}x_{i_{k}}\geq
1-e^{-2N\eta^{2}/\alpha^{2}}\;,$ (171)

where

$\alpha:=\max_{x\in\mathfrak{X}}\ln x-\min_{x\in\mathfrak{X}}\ln
x=\ln\left(x_{max}/x_{min}\right)\;.$ (172)

Proof: The bound (170) on the cardinality follows trivially from the
definition (169) and the fact that

$\displaystyle\sum_{\vec{x}\in{\mathfrak{X}^{(N)}}}\prod_{k=1}^{N}x_{i_{k}}\leq\sum_{\vec{x}\in\mathfrak{X}^{\otimes
N}}\prod_{k=1}^{N}x_{i_{k}}=1\;.$ (173)

To prove (171) it is sufficient to notice that the quantity
$-\ln\prod_{k=1}^{N}x_{i_{k}}=-\sum_{k=1}^{N}\ln x_{i_{k}}$ can be regarded as
the sum of $N$ extractions of the random variable $X:=-\ln x$, that is, of the
variable that with probability $x_{i}$ takes the value $X_{i}=-\ln x_{i}$.
Then the thesis (171) follows straightforwardly from (168). $\square$

##  Appendix B Gibbs states

The Gibbs states of a single-site of our model are the density matrices

$\omega_{\beta}:=\frac{e^{-\beta H}}{Z_{\beta}}\,\;,\quad
Z_{\beta}:=\Tr[e^{-\beta H}]\;,$ (174)

where the parameter $\beta\geq 0$ can be called the _inverse temperature_ of
the system, in analogy with the classical case. They are diagonal in the
eigenbasis $\\{\ket{\epsilon_{i}};i\in[1,d]\\}$ of $H$,

$\omega_{\beta}=\sum_{i=1}^{d}\hat{\lambda}_{i}(\beta)\ket{\epsilon_{i}}\bra{\epsilon_{i}}\;,$
(175)

with population given by

$\displaystyle\hat{\lambda}_{i}(\beta):=Z^{-1}_{\beta}e^{-\beta\epsilon_{i}}\;.$
(176)

The quantity $Z_{\beta}$ is usually called the _partition function_ of the
system and allow us to establish a natural correspondences between the mean
energy of $\omega_{\beta}$ (a quantity that we shall refer to as the
equilibrium energy of the model), its entropy, and the inverse temperature
$\beta$. Specifically we have

$\displaystyle E_{\beta}$ $\displaystyle:=$ $\displaystyle
E(\omega_{\beta};H)=\Tr[\omega_{\beta}H]=-\frac{\partial}{\partial\beta}\ln
Z_{\beta}\;,$ (177) $\displaystyle S_{\beta}$ $\displaystyle:=$ $\displaystyle
S(\omega_{\beta})=-\Tr[\omega_{\beta}\ln\omega_{\beta}]=\beta E_{\beta}+\ln
Z_{\beta}$ (178) $\displaystyle=$
$\displaystyle-\beta\frac{\partial}{\partial\beta}\ln Z_{\beta}+\ln
Z_{\beta}\;,$

which lead to the identity

$\displaystyle\frac{\partial E_{\beta}}{\partial S_{\beta}}=1/\beta\;.$ (179)

The first derivative with respect to $\beta$ of the equilibrium energy (177)
define the heat capacity functional of the model, specifically

$\displaystyle C_{\beta}:=-\frac{\partial
E_{\beta}}{\partial\beta}=\frac{\partial^{2}}{\partial^{2}\beta}\ln
Z_{\beta}\;,$ (180)

which enters in the following Taylor expansions formulas

$\displaystyle E_{\beta^{\prime}}$ $\displaystyle=$ $\displaystyle
E_{\beta}-(\beta^{\prime}-\beta)C_{\beta}+\mathcal{O}\left((\beta^{\prime}-\beta)^{2}\right)\;,$
(181) $\displaystyle\ln Z_{\beta^{\prime}}$ $\displaystyle=$ $\displaystyle\ln
Z_{\beta}-(\beta^{\prime}-\beta)E_{\beta}+\frac{1}{2}(\beta^{\prime}-\beta)^{2}C_{\beta}$
(182)
$\displaystyle\qquad+\mathcal{O}\left((\beta^{\prime}-\beta)^{3}\right)\;,$

(the minus sign in Eq. (180) accounts for the fact that that $\beta$ is an
inverse temperature). The functional $\ln Z_{\beta}$ can be shown to be
decreasing and convex implying the inequality

$\displaystyle\ln Z_{\beta^{\prime}}$ $\displaystyle>$ $\displaystyle\ln
Z_{\beta}+(\beta^{\prime}-\beta)\frac{\partial\ln Z_{\beta}}{\partial\beta}$
(183) $\displaystyle=$ $\displaystyle\ln
Z_{\beta}-(\beta^{\prime}-\beta)E_{\beta}\;,$

valid for $\beta^{\prime}<\beta$. This property also ensures also the
positivity of $C_{\beta}$ which in turns implies that both $E_{\beta}$ and
$S_{\beta}$ are monotonically decreasing functions of $\beta$ in agreement
with Eq. (179). Exploiting these one-to-one correspondences with $\beta$, we
can naturally associate to each entropy value $s$ a thermal energy value
$\mathfrak{E}_{s}$, a heat capacity $\mathfrak{C}_{s}$, and an inverse
temperature $\beta_{s}$ via the identities

$\displaystyle\left\\{\begin{array}[]{l}\mathfrak{E}_{s}:=E_{\beta}\;,\\\ \\\
\mathfrak{C}_{s}:=C_{\beta}\;,\\\ \\\
\beta_{s}:=\beta\;,\end{array}\right.\iff s=S(\omega_{\beta})\;.$ (189)

More generally, given $\rho\in\mathfrak{S}$ a generic single-site state, we
define its associated thermal energy $\mathfrak{E}(\rho)$, an effective heat
capacity $\mathfrak{C}(\rho)$, and inverse effective temperature $\beta(\rho)$
via the identities

$\displaystyle\left\\{\begin{array}[]{l}\mathfrak{E}(\rho):=E_{\beta}\;,\\\
\\\ \mathfrak{C}(\rho):=C_{\beta}\;,\\\ \\\
\beta(\rho):=\beta\;,\end{array}\right.\iff S(\rho)=S(\omega_{\beta})\;.$
(195)

Thanks to these definitions we can rewrite the corresponding total ergotropy
(1) as

$\displaystyle\mathcal{E}_{\text{tot}}(\rho;H)$ $\displaystyle=$
$\displaystyle E(\rho;H)-\mathfrak{E}(\rho)\;.$ (196)

We conclude by noticing that the above construction can be trivially
generalized to the non-interacting $N$ sites model. Specifically from (2) it
follows that the associated Gibbs configurations are tensor product of single
sites Gibbs states, i.e. $\omega^{\otimes N}_{\beta}$. Introducing
$\ket{\epsilon_{\vec{i}}}:=\ket{\epsilon_{i_{1}}}\otimes\dots\otimes\ket{\epsilon_{i_{N}}}$
the eigenvectors of $H^{(N)}$ we then observe that the following identity hold

$\displaystyle
H^{(N)}\ket{\epsilon_{\vec{i}}}=\epsilon_{\vec{i}}\ket{\epsilon_{\vec{i}}}\;,\qquad\omega^{\otimes
N}_{\beta}\ket{\epsilon_{\vec{i}}}=\hat{\lambda}_{\vec{i}}\ket{\epsilon_{\vec{i}}}\;,$
(197)

where

$\displaystyle\epsilon_{\vec{i}}:=\sum_{k=1}^{N}\epsilon_{i_{k}}\;,\qquad\hat{\lambda}_{\vec{i}}(\beta):=\prod_{k=1}^{N}\hat{\lambda}_{i_{k}}(\beta)\;,$
(198)

which from Eq. (175) imply

$\ln\hat{\lambda}_{\vec{i}}(\beta)=N\ln Z_{\beta}+\beta\epsilon_{\vec{i}}\;.$
(199)

###  B.1 A useful Lemma

In this section we provide an estimation of the number of eiegenvalues of
$H^{(N)}$ with energy just above a given threshold linked via Eq. (189) to
single-site entropy values. Specifically, given $s_{0}>0$, $N$ integer, and
$\xi>0$ define

$\displaystyle{\mathfrak{A}_{\xi}^{(N)}}:=\left\\{\epsilon_{\vec{i}}\middle|N\mathfrak{E}_{s_{0}}\leq\epsilon_{\vec{i}}\leq
N(\mathfrak{E}_{s_{0}}+2\xi)\right\\}\;,$ (200)

the subset of the eigenvalues of the Hamiltonian $H^{(N)}$ whose energy share
per site is $2\xi$-close to $\mathfrak{E}_{s_{0}}$. Then the following
property holds:

######  Lemma 2.

Given $\zeta<1$ a constant strictly smaller than 1, for all $s_{0}>0$ there
exists $\xi_{*}\in\;]0,1[$ such that for all $\xi\in\;]0,\xi_{*}[$, we can
identify $N_{*}(\xi)$ integer such that for all $N\geq N_{*}(\xi)$ the
cardinality of ${\mathfrak{A}_{\xi}^{(N)}}$ is bounded by the inequality

$\displaystyle\\#{\mathfrak{A}_{\xi}^{(N)}}\geq
e^{N\left(s_{0}+\frac{\zeta\xi^{2}}{{2}\mathfrak{C}_{s_{0}}}\right)}\;.$ (201)

Proof:– We remind that according to the notation introduced in Eq. (189)
$\mathfrak{E}_{s_{0}}$ is the mean energy of $E_{\beta_{s_{0}}}$ of a single-
site Gibbs state $\omega_{\beta_{s_{0}}}$ with entropy $s_{0}$,
${\beta_{s_{0}}}$ being the corresponding inverse temperature defined as in
Eq. (195). For $s_{0}$ assigned and $\xi>0$ sufficiently small consider the
inverse temperature $\beta^{\prime}$ defined by the identity

$\displaystyle
E_{\beta^{\prime}}=E_{\beta_{s_{0}}}+\xi=\mathfrak{E}_{s_{0}}+\xi\;,$ (202)

which, by construction is slightly smaller than $\beta_{s_{0}}$. By virtue of
Eq. (200) and (199) we have

$\displaystyle\epsilon_{\vec{i}}\in{\mathfrak{A}_{\xi}^{(N)}}$
$\displaystyle\iff$
$\displaystyle\left\lvert\epsilon_{\vec{i}}-NE_{\beta^{\prime}}\right\rvert\leq
N\xi$ (203) $\displaystyle\iff$
$\displaystyle\left\lvert-\ln\hat{\lambda}_{\vec{i}}(\beta^{\prime})-NS_{\beta^{\prime}}\right\rvert\leq
N\beta^{\prime}\xi\;,$

which in particular implies

$\displaystyle\epsilon_{\vec{i}}\in{\mathfrak{A}_{\xi}^{(N)}}\Longrightarrow\hat{\lambda}_{\vec{i}}(\beta^{\prime})\geq
e^{-N(S_{\beta^{\prime}}+\beta^{\prime}\xi)}\;.$ (204)

Remember then that the expected values of
$-\ln\hat{\lambda}_{\vec{i}}(\beta^{\prime})$ for the state $\omega^{\otimes
N}_{\beta^{\prime}}$ are:

$\displaystyle\mathbb{E}[-\ln\hat{\lambda}_{\vec{i}}(\beta^{\prime})]$
$\displaystyle=$ $\displaystyle-\Tr[\omega^{\otimes
N}_{\beta^{\prime}}\ln\omega^{\otimes
N}_{\beta^{\prime}}]=NS_{\beta^{\prime}}\;.$ (205)

Identifying hence the variables $x_{i}$ of Lemma 1 of Appendix A with the
population $\hat{\lambda}_{i}(\beta^{\prime})$ of the Gibbs state
$\omega_{\beta^{\prime}}$, and $s_{0}$ with the associated entropy
$S_{\beta^{\prime}}$, from Eq. (171) we can claim that in the state
$\omega_{\beta^{\prime}}^{\otimes N}$, the set ${\mathfrak{A}_{\xi}^{(N)}}$
hosts a total population of at least

$\sum_{\vec{i}|\epsilon_{\vec{i}}\in{\mathfrak{A}_{\xi}^{(N)}}}\hat{\lambda}_{\vec{i}}(\beta^{\prime})\geq
1-e^{-\frac{2N\xi^{2}}{\epsilon^{2}_{d}}}\;,$ (206)

where according to Eq. (3) $\epsilon_{d}$ is the maximum eigenvalue of $H$.
Equation (203) also implies that

$\epsilon_{\vec{i}}\in{\mathfrak{A}_{\xi}^{(N)}}\Longrightarrow\hat{\lambda}_{\vec{i}}(\beta^{\prime})\leq
e^{-N(S_{\beta^{\prime}}-\beta^{\prime}\xi)}\;.$ (207)

From (206) and (207), it then follows that the cardinality of the set of
eigenvalues ${\mathfrak{A}_{\xi}^{(N)}}$ is at least

$\\#{\mathfrak{A}_{\xi}^{(N)}}\geq
e^{N(S_{\beta^{\prime}}-\beta^{\prime}\xi)}\left(1-e^{-\frac{2N\xi^{2}}{\epsilon^{2}_{d}}}\right)\;.$
(208)

From the second order expansion (182) we get

$\displaystyle e^{N(S_{\beta^{\prime}}-\beta^{\prime}\xi)}$ $\displaystyle=$
$\displaystyle
e^{N\left(s_{0}+\frac{1}{2}\mathfrak{C}_{s_{0}}(\beta_{0}-\beta^{\prime})^{2}+\mathcal{O}\left((\beta_{0}-\beta^{\prime})^{3}\right)\right)}$
(209) $\displaystyle=$ $\displaystyle
e^{N\left(s_{0}+\frac{\xi^{2}}{2\mathfrak{C}_{s_{0}}}+\mathcal{O}\left(\xi^{3}\right)\right)}\;,$

where in the second identity we invoked (181). The term in
$\mathcal{O}\left(\xi^{3}\right)$ could be positive or negative, depending on
the sign of $\frac{d{C}(\beta)}{d\beta}$ for $\beta=\beta_{s_{0}}$. Notice
however that the inequality (183) assures that

$\displaystyle e^{N(S_{\beta^{\prime}}-\beta^{\prime}\xi)}$ $\displaystyle=$
$\displaystyle e^{N\left(\ln
Z_{\beta^{\prime}}+\beta^{\prime}\mathfrak{E}_{s_{0}}\right)}$
$\displaystyle>$ $\displaystyle e^{N\left[\ln
Z_{\beta_{0}}+(\beta_{0}-\beta^{\prime})\mathfrak{E}_{s_{0}}+\beta^{\prime}\mathfrak{E}_{s_{0}}\right]}=e^{Ns_{0}},$

hence implying that
$\frac{\xi^{2}}{2\mathfrak{C}_{s_{0}}}+\mathcal{O}\left(\xi^{3}\right)$ in
(209) is globally positive. Replacing (209) into (208), we now get

$\displaystyle\\#{\mathfrak{A}_{\xi}^{(N)}}$ $\displaystyle\geq$
$\displaystyle\exp[N\left(s_{0}+\tfrac{\xi^{2}}{2\mathfrak{C}_{s_{0}}}+\mathcal{O}\left(\xi^{3}\right)\right)$
(211)
$\displaystyle\qquad+\ln\left(1-\exp[-\tfrac{2N\xi^{2}}{\epsilon^{2}_{d}}]\right)\Big{]}\;.$

The thesis then follows by noticing that for fixed $\xi$ in the limit of large
$N$ the exponent on the right-hand-side of (211) approaches
$N\left(s_{0}+\frac{\xi^{2}}{2\mathfrak{C}_{s_{0}}}+\mathcal{O}\left(\xi^{3}\right)\right)$
which for enough $\xi$ small can be forced to be larger than
$N\left(s_{0}+\zeta\frac{\xi^{2}}{\mathfrak{C}_{s_{0}}}\right)$. To see this
explicitly observe for instance that for each fixed $\xi>0$ there exists
$N_{*}(\xi)$ integer such that for all $N\geq N_{*}(\xi)$ we can ensure

$\displaystyle\frac{\ln(1-\exp[-\tfrac{2N\xi^{2}}{\epsilon^{2}_{d}}])}{N}\geq-\left(\frac{1-\zeta}{2}\right)\frac{\xi^{2}}{2\mathfrak{C}_{s_{0}}}\;,$
(212)

a condition that allows us to replace (211) with

$\\#{\mathfrak{A}_{\xi}^{(N)}}\geq
e^{N\left(s_{0}+\frac{1+\zeta}{2}\tfrac{\xi^{2}}{2\mathfrak{C}_{s_{0}}}+\mathcal{O}\left(\xi^{3}\right)\right)}\;.$
(213)

Now observe that taking $\xi>0$ smaller than some critical value
$\xi^{(1)}_{*}$ which depends upon $s_{0}$ and $\zeta$, we can impose

$\displaystyle\frac{1+\zeta}{2}\frac{\xi^{2}}{2\mathfrak{C}_{s_{0}}}+\mathcal{O}\left(\xi^{3}\right)\geq\zeta\frac{\xi^{2}}{2\mathfrak{C}_{s_{0}}}\;,$
(214)

hence transforming (213) into (201). We conclude noticing that the value
$\xi_{*}$ is obtained by taking the smallest among $\xi^{(1)}_{*}$ and the
threshold $\epsilon_{d}-\mathfrak{E}_{s_{0}}$ needed to ensure that the
effective inverse temperature $\beta^{\prime}$ introduced in Eq. (202) is
properly defined, i.e.
$\xi_{*}:=\min\\{\xi^{(1)}_{*},\epsilon_{d}-\mathfrak{E}_{s_{0}},1\\}$.
$\square$

Remark 4: It is worth stressing that the parameter $\xi_{*}$ introduced in
Lemma2 is a function of $s_{0}$, $\zeta$, and, due to the presence of
$\mathfrak{C}_{s_{0}}$ and $\mathfrak{E}_{s_{0}}$ in Eqs. (211) and (202), of
the Hamiltonian $H$, i.e. $\xi_{*}=\xi_{*}(s_{0},\zeta,H)$; similar
considerations hold also for $N_{*}(\xi)$ which, besides depending upon
$\xi\leq\xi_{*}$, it is also a function of $s_{0}$, $\zeta$, and $H$, i.e.
$N_{*}(\xi)=N_{*}(\xi,s_{0},\zeta,H)$. In particular, as discussed in Remark 5
below, a not necessarily optimal choice of $N_{*}(\xi,s_{0},\zeta,H)$ which is
however sufficient to ensure (212), is

$\displaystyle N_{*}(\xi)={K_{*}}/{\xi^{2}}\;,$ (215)

with $K_{*}=K_{*}(s_{0},\zeta,H)$ a factor that depends upon $s_{0}$, $\zeta$,
and $H$ – see Eq. (222) for details. Notice also that in our analysis, the
choice of the parameter $\zeta$ in $]0,1[$ is free, however the higher we take
it the larger becomes $N_{*}(\xi)$ reducing the range of $N$ for which Eq.
(201) applies. In an effort to reduce the number of parameters, in the
remaining of the paper we shall fix such constant equal to

$\displaystyle\zeta=2/2.01=0.995025\;,$ (216)

we stress however that none of the results that follow depend crucially on
such a choice.

Remark 5: In order to find an estimation of $N_{*}(\xi,s_{0},H)$ so that Eq.
(212) holds for all $N\geq N_{*}(\xi,s_{0},H)$, let us rewrite such inequality
as

$\displaystyle\frac{\ln(1-e^{-x})}{x}\geq-\left(\frac{1-\zeta}{4}\right)\frac{\epsilon^{2}_{d}}{\mathfrak{C}_{s_{0}}}\;,$
(217)

with $x:=\tfrac{2N\xi^{2}}{\epsilon^{2}_{d}}$. Notice then that for

$\displaystyle x\geq\ln(\frac{e}{e-1})\;,$ (218)

we can write

$\displaystyle\frac{\ln(1-e^{-x})}{x}\geq-\frac{1}{x}\;,$ (219)

whose right-hand-side is larger than the right-hand-side of (217) for

$\displaystyle
x\geq\left(\frac{4}{1-\zeta}\right)\frac{\mathfrak{C}_{s_{0}}}{\epsilon^{2}_{d}}\;.$
(220)

Therefore enforcing $x$ to fulfil both (218) and (220), i.e. imposing

$\displaystyle N\geq\frac{K_{*}(s_{0},\zeta,H)}{\xi^{2}}\;,$ (221)

with

$K_{*}(s_{0},\zeta,H):=\frac{1}{2}\max\left\\{{\epsilon_{d}^{2}}\ln(\frac{e}{e-1}),\frac{4{\mathfrak{C}_{s_{0}}}}{1-\zeta}\right\\}\;,$
(222)

we can ensure that (217) (i.e. (212)) applies, hence proving Eq. (215). As
mentioned in the main text this choice for $N_{*}(\xi,s_{0},\zeta,H)$ is
arguably not optimal as it relay on the correct but drastic simplification
(219).

##  Appendix C Generalization of $\rho^{(N)}_{{\bf B}}$ to non-uniform
partitions

Here we generalize the construction of Sec. III.4.2 to the case in which the
partition ${\mathfrak{P}}$ is not necessarily uniform. Also in this case we
introduce $\rho^{(N)}_{{\bf B}}$ via the expression (88) with the matrices
$L^{2}\times L^{2}$ matrices $A^{ij}_{\bf{B}}$ expressed as in Eq. (87). In
this case however we observe that for $i\neq j$, depending on the
cardinalities of the associated subsets ${\mathfrak{K}}_{\ell^{[i]}}$ and
${\mathfrak{K}}_{\ell^{[j]}}$, the indexes $k^{[i]}$ and $k^{[j]}$ can run
over sets of different sizes: the rectangular matrix
$\delta_{k^{[i]},k^{[j]}}$ appearing in (87) should hence be interpret as the
natural generalization of of Kronecker delta which take same values of the
latter on the common subsets of the indices and which is zero everywhere else;
this prescription ensures that all the identities from Eqs. (III.4.2) to (90)
still hold. In deriving the equivalent formula of (92) we need however some
extra precaution. Again the problem is related to the fact that if the
partition ${\mathfrak{P}}$ is not uniform then the subsets
${\mathfrak{K}}_{\ell}$ have different cardinalities. For compensate for this
fact, we extend these sets adding extra zero elements to push their effective
cardinality to the maximum value $\\#_{\max}$; formally this is obtained by
replacing ${{\mathfrak{K}}}_{\ell}$ with the new set

$\displaystyle\tilde{{\mathfrak{K}}}_{\ell}:=\\{\tilde{\lambda}_{k,\ell};k=1,\cdots,\\#_{\max}\\}\;,$
(223)

with

$\displaystyle\tilde{\lambda}_{k,\ell}:=\left\\{\begin{array}[]{ll}{\lambda}_{k,\ell}\;,&\quad\forall
k=1,\cdots,\\#{\mathfrak{K}}_{\ell}\;,\\\ \\\ 0\;,&\quad\forall
k=\\#{{\mathfrak{K}}}_{\ell}+1,\cdots,\\#_{\max}\;.\end{array}\right.$ (227)

Similarly we define a new set of orthonormal single-site vectors
$\\{|\tilde{\lambda}_{k,\ell}\rangle;k=1,\cdots,\\#_{\max}\\}$ with the
prescription that the first $\\#{\mathfrak{K}}_{\ell}$ elements fulfil the
condition

$\displaystyle|\tilde{\lambda}_{k,\ell}\rangle=|{\lambda}_{k,\ell}\rangle\;,\qquad\forall
k=1,\cdots,\\#{\mathfrak{K}}_{\ell}\;,$ (228)

while the remaining $\\#_{\max}-{\mathfrak{K}}_{\ell}$ can be chosen freely.
Given then an $N$-uple $\vec{k}=\left(k_{1},k_{2},\dots k_{N}\right)$ where
each component can now assume up to $\\#_{\max}$ distinct values, we introduce
the quantities

$\displaystyle\tilde{\lambda}^{(N)}_{\vec{k},\ell}$ $\displaystyle:=$
$\displaystyle\prod_{a=1}^{N}\tilde{\lambda}_{k_{a},\ell}\;,$ (229)

and the $N$ site states

$\displaystyle\ket{\tilde{\lambda}^{(N)}_{\vec{k},\ell}}$ $\displaystyle:=$
$\displaystyle\ket{\tilde{\lambda}_{k_{1},\ell}\;\tilde{\lambda}_{k_{2},\ell}\;\dots\;\tilde{\lambda}_{k_{N},\ell}}\;.$
(230)

Notice that the $\tilde{\lambda}^{(N)}_{\vec{k},\ell}$ correspond to the
positive terms of Eq. (95) if $k_{a}\leq\\#{\mathfrak{K}}_{\ell}$ for all
$a=1,\cdots,N$, and are instead equal to zero otherwise; notice also that for
each chosen $\vec{k}=\left(k_{1},k_{2},\dots k_{N}\right)$ we can ensure that
there exists at least one value of $\ell$ such that
$\tilde{\lambda}^{(N)}_{\vec{k},\ell}>0$ (this follows from the fact that
since $\\#_{\max}$ is the greatest of all $\\#{\mathfrak{K}}_{\ell}$, there is
at least one value of $\ell$ for which
$\tilde{{\mathfrak{K}}}_{\ell}={\mathfrak{K}}_{\ell}$): therefore we can
conclude that the terms

$\tilde{\Lambda}^{(N)}_{\vec{k}}:=\sum_{\ell=1}^{L}\frac{\tilde{\lambda}^{(N)}_{\vec{k},\ell}}{\mathscr{S}_{\ell}^{N-1}}\;,$
(231)

are not null for all $\vec{k}$.

With the help of these definitions can now replace (92) and (96) with

$\displaystyle\rho^{(N)}_{{\bf B}}$ $\displaystyle=$
$\displaystyle\sum_{\vec{k}}\sum_{\ell,\ell^{\prime}=1}^{L}\sqrt{\tfrac{\tilde{\lambda}^{(N)}_{\vec{k},\ell}}{\mathscr{S}_{\ell}^{N-1}}}\sqrt{\tfrac{\tilde{\lambda}^{(N)}_{\vec{k},\ell^{\prime}}}{\mathscr{S}_{\ell^{\prime}}^{N-1}}}\ket{\tilde{\lambda}^{(N)}_{\vec{k},\ell}}\bra{\tilde{\lambda}^{(N)}_{\vec{k},\ell^{\prime}}}$
(232) $\displaystyle=$
$\displaystyle\sum_{\vec{k}}\tilde{\Lambda}_{\vec{k}}^{(N)}\;|\tilde{\Psi}^{(N)}_{\vec{k}}\rangle\langle\tilde{\Psi}^{(N)}_{\vec{k}}|\;,$

with

$|\tilde{\Psi}^{(N)}_{\vec{k}}\rangle:=\tfrac{1}{\sqrt{\tilde{\Lambda}_{\vec{k}}^{(N)}}}\sum_{\ell=1}^{L}\sqrt{\tfrac{\tilde{\lambda}^{(N)}_{\vec{k},\ell}}{\mathscr{S}_{\ell}^{N-1}}}\;\ket{\tilde{\lambda}^{(N)}_{\vec{k},\ell}}\;,$
(233)

being orthonormal elements of $\mathfrak{S}^{(N)}$. It is worth stressing that
also in this case the state (232) is properly normalized thanks to the fact
that

$\displaystyle\sum_{\vec{k}}\tilde{\Lambda}_{\vec{k}}^{(N)}$ $\displaystyle=$
$\displaystyle\sum_{\vec{k}}\sum_{\ell=1}^{L}\mathscr{S}_{\ell}^{1-N}\left(\prod_{a=1}^{N}\tilde{\lambda}_{k_{a},\ell}\right)$
$\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\mathscr{S}_{\ell}^{1-N}\left(\sum_{k=1}^{\\#_{\max}}\tilde{\lambda}_{k,\ell}\right)^{N}$
$\displaystyle=$
$\displaystyle\sum_{\ell=1}^{L}\mathscr{S}_{\ell}^{1-N}\left(\sum_{k=1}^{\\#{\mathfrak{K}}_{\ell}}\lambda_{k,\ell}\right)^{N}=\sum_{\ell=1}^{L}\mathscr{S}_{\ell}=1\;.$

As anticipated in the main text, while the BLR $\mathfrak{br}[\rho^{(N)}_{{\bf
B}}]$ is upper bounded by $L^{2}$ the rank of the state $\rho^{(N)}_{{\bf B}}$
is given by the total number of $N$-uple $\vec{k}$, leading to Eq. (100).

## References

  * Thomson (1852) W. Thomson, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 4, 8 (1852).
  * Fermi  _et al._ (1955) E. Fermi, P. Pasta, S. Ulam, and M. Tsingou,  _Studies of the nonlinear problems_, Tech. Rep. (1955).
  * Allahverdyan  _et al._ (2004) A. E. Allahverdyan, R. Balian, and T. M. Nieuwenhuizen, Europhysics Letters (EPL) 67, 565 (2004).
  * Niedenzu  _et al._ (2019) W. Niedenzu, M. Huber, and E. Boukobza, Quantum 3, 195 (2019).
  * Alicki and Fannes (2013) R. Alicki and M. Fannes, Phys. Rev. E 87, 042123 (2013).
  * Huber  _et al._ (2015) M. Huber, M. Perarnau-Llobet, K. V. Hovhannisyan, P. Skrzypczyk, C. Klöckl, N. Brunner, and A. Acín, New J. Phys. 17, 065008 (2015).
  * Francica  _et al._ (2017) G. Francica, J. Goold, F. Plastina, and M. Paternostro, npj Quantum Information 3 (2017), 10.1038/s41534-017-0012-8.
  * Oppenheim  _et al._ (2002) J. Oppenheim, M. Horodecki, P. Horodecki, and R. Horodecki, Phys. Rev. Lett. 89 (2002), 10.1103/physrevlett.89.180402.
  * Vitagliano  _et al._ (2018) G. Vitagliano, C. Klöckl, M. Huber, and N. Friis, in  _Fundamental Theories of Physics_ (Springer International Publishing, 2018) pp. 731–750.
  * Goold  _et al._ (2016) J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, J. Phys. A: Math. Th.l 49, 143001 (2016).
  * Bera  _et al._ (2017) M. N. Bera, A. Riera, M. Lewenstein, and A. Winter, Nat. Comm. 8 (2017), 10.1038/s41467-017-02370-x.
  * Manzano  _et al._ (2018) G. Manzano, F. Plastina, and R. Zambrini, Phys. Rev. Lett. 121 (2018), 10.1103/physrevlett.121.120602.
  * Andolina  _et al._ (2019) G. M. Andolina, M. Keck, A. Mari, M. Campisi, V. Giovannetti, and M. Polini, Phys. Rev. Lett. 122 (2019), 10.1103/physrevlett.122.047702.
  * Verstraete  _et al._ (2004) F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
  * Baumgratz  _et al._ (2014) T. Baumgratz, M. Cramer, and M. B. Plenio, Phys. Rev. Lett. 113, 140401 (2014).
  * Francica  _et al._ (2020) G. Francica, F. C. Binder, G. Guarnieri, M. T. Mitchison, J. Goold, and F. Plastina, Phys. Rev. Lett. 125, 180603 (2020).
  * Çakmak (2020) B. Çakmak, Phys. Rev. E 102, 042111 (2020).
  * Touil  _et al._ (2021) A. Touil, B. Çakmak, and S. Deffner, “Second law of thermodynamics for quantum correlations,” (2021), arXiv:2102.13606 .
  * Kliesch  _et al._ (2014) M. Kliesch, D. Gross, and J. Eisert, Phys. Rev. Lett. 113 (2014), 10.1103/physrevlett.113.160503.
  * Navascues and Vertesi (2018) M. Navascues and T. Vertesi, Quantum 2, 50 (2018).
  * Parker  _et al._ (2020) D. E. Parker, X. Cao, and M. P. Zaletel, Phys. Rev. B 102 (2020), 10.1103/physrevb.102.035147.
  * Perarnau-Llobet and Uzdin (2019) M. Perarnau-Llobet and R. Uzdin, New J. Phys. 21, 083023 (2019).
  * Correa  _et al._ (2015) L. A. Correa, M. Mehboudi, G. Adesso, and A. Sanpera, Phys. Rev. Lett. 114 (2015), 10.1103/physrevlett.114.220405.
  * Mirsky (1975) L. Mirsky, Monatshefte für Mathematik 79, 303 (1975).
  * Alimuddin  _et al._ (2020) M. Alimuddin, T. Guha, and P. Parashar, Phys. Rev. E 102 (2020), 10.1103/physreve.102.012145.
  * Schindler  _et al._ (2020) J. Schindler, D. Šafránek, and A. Aguirre, Phys. Rev. A 102 (2020), 10.1103/physreva.102.052407.
  * Rolandi and Wilming (2020) A. Rolandi and H. Wilming, “Extensive Rényi entropies in matrix product states,” (2020), arXiv:2008.11764 .
  * Chernoff (1952) H. Chernoff, The Annals of Mathematical Statistics 23, 493 (1952).
  * Hoeffding (1963) W. Hoeffding, J. Am. Stat. Assoc. 58, 13 (1963).