# Thermal Area Law in Long-Range Interacting Systems

Donghoon Kim<EMAIL_ADDRESS>Analytical quantum complexity RIKEN Hakubi
Research Team, RIKEN Cluster for Pioneering Research (CPR)/ RIKEN Center for
Quantum Computing (RQC), Wako, Saitama 351-0198, Japan  Tomotaka Kuwahara
<EMAIL_ADDRESS>Analytical quantum complexity RIKEN Hakubi Research
Team, RIKEN Cluster for Pioneering Research (CPR)/ RIKEN Center for Quantum
Computing (RQC), Wako, Saitama 351-0198, Japan  PRESTO, Japan Science and
Technology (JST), Kawaguchi, Saitama 332-0012, Japan  Keiji Saito
<EMAIL_ADDRESS>Department of Physics, Kyoto University,
Kyoto 606-8502, Japan

###### Abstract

The area law of the bipartite information measure characterizes one of the
most fundamental aspects of quantum many-body physics. In thermal equilibrium,
the area law for the mutual information universally holds at arbitrary
temperatures as long as the systems have short-range interactions. In systems
with power-law decaying interactions, $r^{-\alpha}$ ($r$: distance),
conditions for the thermal area law are elusive. In this work, we aim to
clarify the optimal condition $\alpha>\alpha_{c}$ such that the thermal area
law universally holds. A standard approach to considering the conditions is to
focus on the magnitude of the boundary interaction between two subsystems.
However, we find here that the thermal area law is more robust than this
conventional argument suggests. We show the optimal threshold for the thermal
area law by $\alpha_{c}=(D+1)/2$ ($D$: the spatial dimension of the lattice),
assuming a power-law decay of the clustering for the bipartite correlations.
Remarkably, this condition encompasses even the thermodynamically unstable
regimes $\alpha<D$. We verify this condition numerically, finding that it is
qualitatively accurate for both integrable and non-integrable systems.
Unconditional proof of the thermal area law is possible by developing the
power-law clustering theorem for $\alpha>D$ above a threshold temperature.
Furthermore, the numerical calculation for the logarithmic negativity shows
that the same criterion $\alpha>(D+1)/2$ applies to the thermal area law for
quantum entanglement.

Introduction.— Quantum correlation and entanglement play pivotal roles in
understanding quantum many-body systems from an information-theoretic
standpoint. They help in identifying the exotic quantum phases [1, 2, 3] and
serve as crucial resources in quantum information processing [4]. One hallmark
of quantum many-body systems is the area law for the ground state [5], which
posits that the correlation and entanglement between two subsystems are
constrained by the surface area of their interface. This principle has
undergone extensive verification through both analytical and numerical
approaches in a wide array of scenarios [6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16]. The area law not only elucidates the inherent complexities of quantum
systems but also significantly enhances their simulatability using classical
computers [17, 18].

The quantum Gibbs state describing thermal equilibrium obeys an analogous
principle, known as the thermal area law, as long as the systems have short-
range interactions [19, 20, 21, 22, 23]. This law is reflected by the behavior
of mutual information, which contains both classical and quantum correlations,
as illustrated by the following inequality:

$\displaystyle\mathcal{I}(A:B)$ $\displaystyle\leq 2\beta\|H_{\partial
A}\|\,,$ (1)

where $\mathcal{I}$ represents the mutual information of the Gibbs state
between two subsystems $A$ and $B$ at an inverse temperature $\beta$. Here,
$\|H_{\partial A}\|$ denotes the operator norm of the boundary interaction
Hamiltonian, defined as the sum of interactions whose supports intersect with
the boundary $\partial A$ of $A$. In the short-range interacting case, the
norm of these boundary interactions scales with the boundary area $|\partial
A|$, thereby affirming the thermal area law. This law is not only relevant to
mutual information but also extends to a wider range of bipartite information
measures [20, 21, 22]. We also note that a recent study has further clarified
that genuine quantum correlations cannot be as large as the system size for
short-range interacting systems [24].

Figure 1: Schematic of the results. The thermal area law is satisfied for
$\alpha>(D+1)/2$ (indicated by a red region). The thermal area law is much
more robust than the relation (1) suggests (a blue region). The key physics
behind the area law is the power-law clustering of bipartite correlations.

In contrast to short-range interactions, the general aspects of locality in
long-range interacting systems remain less understood. Long-range interactions
here refer to interactions that decay as $r^{-\alpha}$ with distance $r$ [25],
which are ubiquitous in nature such as in atomic, molecular, and optical
systems [26, 27, 28, 29, 30, 31]. These interactions lead to non-trivial
phenomena absent in short-range systems [32, 33, 34, 35, 36, 37, 38, 39, 40].
Key unresolved issues include (i) identifying the critical threshold
$\alpha_{\rm c}$ above which any systems have the thermal area law, (ii)
understanding the underlying physics, and (iii) distinguishing between
classical and quantum correlations. This paper addresses these questions as
fundamental physics. In a practical aspect, this direction may lead to finding
an efficiency-guaranteed algorithm for simulating the long-range quantum Gibbs
state. Should the fundamental physics be governed solely by the magnitude of
the boundary interaction, then the inequality (1) could already be considered
optimal; that is, the condition of the thermal area law is given by
$\alpha>D+1$ (or $\alpha_{c}=D+1$), which is required to ensure $\|H_{\partial
A}\|\propto|\partial A|$ as in [39]. Indeed, the threshold given by the finite
boundary interaction recurs across various contexts in physics, indicating
that physical properties undergo qualitative changes at this threshold in
various many-body systems [41, 42, 43, 44, 45, 39].

However, in this work, we discover that the underlying physics of the thermal
area law is influenced by a different characteristic: the clustering property,
which is expected to remain robust in any noncritical phases. Note that the
clustering property has been widely adopted as a basic and natural assumption
in discussing various quantum properties such as entanglement area law [14],
quantum Gibbs sampling [46, 47], the eigenstate thermalization [48, 49], etc.
Consequently, we show that the thermal area law for mutual information
possesses a greater resilience than traditionally anticipated. Utilizing the
power-law clustering properties in bipartite correlations, we analytically
find that the sufficient condition for the thermal area law in the mutual
information is given $\alpha>(D+1)/2$. Regarding the critical threshold
$\alpha_{c}$ above which any systems satisfy the thermal area law, our
numerical analysis concludes that $\alpha_{c}=(D+1)/2$, since below
$\alpha_{c}$ we provide explicit examples violating the area law. Remarkably,
this condition encompasses $\alpha<D$, where the thermodynamic extensivity is
violated [25]. In this regime, most analytical techniques break down, and
universal behaviors remain largely unexplored. While the numerical
investigation suggests that the clustering property robustly holds even for
thermodynamic nonextensive regimes, we demonstrate that the power-law
clustering theorem can be rigorously derived for $\alpha>D$ by identifying a
suitable temperature regime. Consequently, the thermal area law
unconditionally holds for $\alpha>D$. Our key results are depicted
schematically in Fig. 1.

In addition, we investigate whether differences exist between quantum and
classical correlations in the context of the thermal area law. By employing
the logarithmic negativity proposed by Shapourian, Shiozaki, and Ryu for
bilinear fermionic systems [50], we demonstrate that this criterion equally
governs the thermal area law for quantum entanglement.

Setup.— We consider a quantum system with $N$ qudits or fermions on the
$D$-dimensional lattice. Let $\Lambda$ be the set of all sites on the lattice
of $N$ sites. We consider the Hamiltonian $H$ in the $k$-local form:

$\displaystyle H=\sum_{Z:|Z|\leq k}h_{Z}\,,$ (2)

where $Z$ is a set of interaction sites and $|Z|$ denotes the number of
interacting sites (i.e., $|Z|$-body interaction). We consider the long-range
interactions satisfying the following condition

$\displaystyle
J_{i,i^{\prime}}:=\sum_{Z:Z\ni\\{i,i^{\prime}\\}}\left\|h_{Z}\right\|\leq{g\over(1+d_{i,i^{\prime}})^{\alpha}}\,.$
(3)

Here, $J_{i,i^{\prime}}$ is the maximum norm of the local Hamiltonian
including the sites $i$ and $i^{\prime}$, and $g$ is some constant. The symbol
$d_{i,i^{\prime}}$ stands for the Manhattan distance between the sites $i$ and
$i^{\prime}$. The parameter $\alpha$ is the index of power-law decay in the
interaction. We will provide examples of the Hamiltonian in (9) and (10)
later. Note that the Hamiltonian does not contain the Kac factor for the
thermodynamically nonextensive regime $\alpha<D$, since such a factor does not
exist in realistic situations [28, 51].

The quantum Gibbs state $\rho_{\beta}$ at an inverse temperature $\beta$ is
defined as $\rho_{\beta}:=e^{-\beta H}/{\rm tr}(e^{-\beta H})$. The mutual
information of $\rho_{\beta}$ between two regimes $A$ and $B$ is defined as

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B):=S(\rho_{\beta}^{A})+S(\rho_{\beta}^{B})-S(\rho_{\beta}^{AB}),$
(4)

where $\rho_{\beta}^{X}$ is the reduced density matrix of $\rho_{\beta}$ on
the subset $X\subseteq\Lambda$, and $S(\rho_{\beta}^{X}):=-{\rm
tr}(\rho_{\beta}^{X}\log\rho_{\beta}^{X})$ is the von Neumann entropy. We
define the correlation function of two observables $O_{A}$ and $O_{B}$ as

$\displaystyle{\rm Cor}_{\rho_{\beta}}(O_{A},O_{B}):={\rm
tr}(\rho_{\beta}O_{A}O_{B})-{\rm tr}(\rho_{\beta}O_{A})\cdot{\rm
tr}(\rho_{\beta}O_{B}).$ (5)

Main theorem.— We find that the thermal area law of the mutual information
holds above some power-law threshold with the assumption of clustering.

Theorem 1: Let us assume that the correlation function of the quantum Gibbs
state $\rho_{\beta}$ between two arbitrary operators $O_{i}$ and
$O_{i^{\prime}}$ supported on the sites $i$ and $i^{\prime}$ satisfies the
following power-law clustering property:

$\displaystyle{\rm
Cor}_{\rho_{\beta}}(O_{i},O_{i^{\prime}})\leq\frac{C}{d_{i,i^{\prime}}^{\alpha}}\|O_{i}\|\cdot\|O_{i^{\prime}}\|,$
(6)

where $C$ is an $\mathcal{O}(1)$ constant. Then, for the bipartition $A,B$ of
$D$-dimensional lattice $\Lambda$ $(A\cup B=\Lambda)$ and $\alpha>(D+1)/2$,
the mutual information is upper bounded by

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)\leq\text{const.}\beta|\partial
A|\,.$ (7)

We postpone the proof to the end of the paper (the details are in the
supplementary material (SM) [52]). Below, we focus on physically crucial
points. The theorem ensures that the thermal area law holds for the regime
$\alpha>(D+1)/2$. This condition encompasses a broader scope compared to the
condition found in previously established thermal area law (1). It is known
that (e.g., Ref. [39]) the operator norm $\|H_{\partial A}\|$ in (1) is upper-
bounded by the boundary area between $A$ and $B$ when $\alpha>D+1$ for generic
interacting systems (the condition for the bilinear systems is
$\alpha>D/2+1$). Hence, the condition given in the theorem means that the
thermal area law is much more robust than expected from a simple argument on
the norm of the boundary interaction. Remarkably, the condition partially
includes the thermodynamically nonextensive regime.

Power-law clustering and unconditional thermal area law.— Theorem 1 suggests
that the crucial physics behind the thermal area law is not only from the
magnitude of the boundary interaction but also from the clustering property.
Note that the clustering property is expected to be one of the robust physical
properties in non-critical thermal phases, rendering it an inherently natural
assumption [14]. We here demonstrate that the power-law clustering property
can be rigorously proven for $\alpha>D$ above a temperature threshold as shown
in the following theorem 2.

Theorem 2: Under the general Hamiltonian (2) with (3) for the regime
$\alpha>D$, the following clustering property holds for the temperatures above
$\beta_{c}^{-1}$:

$\displaystyle{\rm Cor}_{\rho_{\beta}}(O_{X},O_{Y})$ $\displaystyle\leq
C\|O_{X}\|\cdot\|O_{Y}\|{|X||Y|\,e^{{(|X|+|Y|)/k}}\over d_{X,Y}^{\alpha}}\,,$
(8)

where $d_{X,Y}=\min_{i\in X,j\in Y}d_{i,j}$ and $|X|$ and $|Y|$ are numbers of
sites in the regions $X$ and $Y$, respectively. $C$ is an $\mathcal{O}(1)$
constant and the threshold temperature is given by $\beta_{c}=1/(8ugk)$ with
$u=\sum_{j\in\Lambda}1/(1+d_{i,j})^{\alpha}$.

The proof is based on the cluster expansion technique. The details in the
proof are provided in the SM [52]. To the best of our knowledge, this is the
first result to establish the power-law clustering theorem in quantum long-
range interacting systems.

From this statement, the clustering property in Theorem 1 is not an assumption
but rigorously holds above a threshold temperature as long as the
thermodynamically extensive regime (i.e., $\alpha>D$) is considered. In
particular, in 1D systems, the clustering property is believed to hold at
arbitrary temperatures [53, 54]. Moreover, we stress that even in the
thermodynamically nonextensive regime ($\alpha<D$), numerical calculations for
the two models below show the validity of the clustering property.

Numerical verification.— We numerically verify our theorems in both integrable
and nonintegrable systems, establishing the bound’s tightness.

As a typical example of integrable systems, we use the following long-range
bilinear fermion system:

$\displaystyle
H=-\sum_{i,j\in\Lambda}\frac{t_{i,j}}{d_{i,j}^{\alpha}}(c_{i}^{\dagger}c_{j}+c_{j}^{\dagger}c_{i}),$
(9)

where $t_{i,j}$ is a hopping parameter of order $\mathcal{O}(1)$, and $c_{i}$
and $c_{i}^{\dagger}$ are the annihilation and creation operators of the
spinless fermion at site $i$, respectively. We consider the Hamiltonian in the
one dimension (1D) and the two dimension (2D).

We first verify the clustering property (6) in the Hamiltonian (9) [55, 56].
We analyze a 1D chain of $N=1000$ sites and a 2D square lattice of side length
$N=40$, focusing on two-point correlation functions $\langle
c_{i}^{\dagger}c_{i+r}\rangle$ in 1D and $\langle
c_{{\bm{i}}}^{\dagger}c_{{\bm{i}}+{\bm{r}}}\rangle$ in 2D, setting $\beta=2$.
To consider general bilinear systems, we assign randomly chosen values
$\in[0,1]$ to the hopping variable, and take the average over $1000$ samples.
In Figs. 2(a) and (b), we observe that $|\langle
c_{i}^{\dagger}c_{i+r}\rangle|\times r^{\alpha}$ and $|\langle
c_{{\bm{i}}}^{\dagger}c_{{\bm{i}}+{\bm{r}}}\rangle|\times|{\bm{r}}|^{\alpha}$
saturate at large distances for any $\alpha$, thus verifying the clustering
property in both dimensions. Note also that the clustering property holds even
for $\alpha<D$. While we show here the averaged data over samples, even
individual data show the clustering property (not shown). This indicates all
correlation functions adhere to this property due to the Wick theorem.

After establishing the clustering property, we affirm the thermal area law is
applicable for $\alpha>(D+1)/2$, a sufficient condition whose optimality
warrants examination. Let us investigate how tight the condition is, through
the numerical calculation with the same model. The system is partitioned into
subsystems $A$ and $B$, with $A$ including the first $N/2$ sites and $B$ the
rest in 1D. For 2D, $A$ comprises the initial $(N\times N/2)$ sites, and $B$
the remaining. The mutual information $\mathcal{I}_{\rho_{\beta}}(A:B)$ is
calculated with the same parameters from the clustering property study.
Results are presented in Figs. 2(c) for 1D and (d) for 2D. The mutual
information increases with the system size if $\alpha<1$ and it saturates to a
constant value if $\alpha>1$ in 1D. In 2D, the mutual information divided by
the boundary area $N$ grows for $\alpha<1.5$ and steadies for $\alpha>1.5$.
These numerical results support the thermal area law holds for
$\alpha>(D+1)/2$, and hence the condition is optimal.

Figure 2: 1D and 2D long-range bilinear fermions. (a) The two-point
correlation function $|\langle c_{i}^{\dagger}c_{i+r}\rangle|$ in the 1D chain
of $N=1000$ sites between the $(i=250)$th site and $(i+r)$-th site multiplied
by $r^{\alpha}$ for different $\alpha$’s. (b) The two-point correlation
function $|\langle c_{{\bm{i}}}^{\dagger}c_{{\bm{i}}+{\bm{r}}}\rangle|$ in the
2D square lattice with a side length $N=40$ between the site
${\bm{i}}=(10,10)$ and the site ${\bm{i}}+{\bm{r}}=(10+r,10+r)$ multiplied by
$|{\bm{r}}|^{\alpha}$ for different $\alpha$’s. (c) The mutual information
$\mathcal{I}_{\rho_{\beta}}(A:B)$ in the 1D chain of sites $N$ between half
($A$) of the system and the other half ($B$). (d) The mutual information
$\mathcal{I}_{\rho_{\beta}}(A:B)$ in the 2D square lattice with a side length
$N$ between half ($A$) of the system and the other half ($B$) divided by $N$.
All the figures are the average over the 1000 samples of random variables
$t_{i,j}\in[0,1]$ in Eq. (9) at $\beta=2$.

In general, bilinear systems can be special in several physical aspects not
only for the non-equilibrium properties [57] but also for the static
properties including the thermal entanglement [58, 59, 60]. However, the
clustering is expected to be robust, regardless of the (non)integrability, and
hence the thermal area law should also hold universally. We here check this
universality using the specific nonintegrable system. Let us consider the one-
dimensional long-range Heisenberg spin-1/2 chain of $N$ sites

$\displaystyle H=\sum_{1\leq i<j\leq
N}\frac{1}{d_{i,j}^{\alpha}}{\bm{S}}_{i}\cdot{\bm{S}}_{j},$ (10)

with a spin-1/2 operator ${\bm{S}}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$ at
site $i$. Exploring the clustering property of
$\langle{\bm{S}}_{i}\cdot{\bm{S}}_{i+r}\rangle$ with $N=200$ and $\beta=2$, we
employ the exponential tensor renormalization group (XTRG) algorithm, as
detailed in references [61, 62], in conjunction with the QSpace tensor library
[63, 64]. This approach enables us to construct the matrix product operator
for the Gibbs state $\rho_{\beta}$. See [65] for details. Fig. 3(a) presents
$|\langle{\bm{S}}_{i}\cdot{\bm{S}}_{i+r}\rangle|\times r^{\alpha}$, evidencing
clustering by constant upper bounds. This observation, as well as the bilinear
Hamiltonian case, strongly suggests that the clustering universally holds.

Figure 3: 1D long-range Heisenberg chain. (a) The two-point correlation
function $|\langle{\bm{S}}_{i}\cdot{\bm{S}}_{i+r}\rangle|$ multiplied by
$r^{\alpha}$ in the 1D chain with $N=200$ sites between $(i=50)$th site and
$(i+r)$th site at $\beta=2$ for different $\alpha$’s. (b) The entanglement
entropy Eq. (12) of the thermofield double state of the Gibbs state at
$\beta=2$ for different $\alpha$’s. We choose two subsystems as half of the
system and the other half and change the size $N$ of the chain. Here, we use
the XTRG algorithm and set bond dimensions ($\mathrm{SU}(2)$ multiplets) of
210 and 120 for (a) and (b), which correspond approximately 850 and 450
states, respectively.

To verify the thermal area law, we compute the entanglement entropy
($\mathcal{E}$) for the thermofield double (TFD) state of Gibbs state
$\rho_{\beta}$ as it upper bounds the mutual information [66]. We consider two
identical copies $\mathcal{H}_{L}$ and $\mathcal{H}_{R}$ of the original
Hilbert space $\mathcal{H}$. For the eigenvalues $\\{E_{n}\\}$ and eigenstates
$\\{|n\rangle\\}$ of the Hamiltonian, the TFD state of the Gibbs state
$\rho_{\beta}=\sum_{n}e^{-\beta E_{n}}/{\rm tr}(e^{-\beta H})|n\rangle\langle
n|$ is defined as

$\displaystyle|\mathrm{TFD}\rangle:=\frac{1}{\sqrt{{\rm tr}(e^{-\beta
H})}}\sum_{n}e^{-\beta E_{n}/2}|n\rangle_{L}\otimes|n\rangle_{R}.$ (11)

Partitioning the system into halves $A$ and $B$, we define $\mathcal{E}$ as
the von Neumann entropy of the reduced density matrix
$\sigma_{A_{L},A_{R}}={\rm
tr}_{B_{L},B_{R}}(|\mathrm{TFD}\rangle\langle\mathrm{TFD}|)$:

$\displaystyle\mathcal{E}:=S(\sigma_{A_{L},A_{R}})=-{\rm
tr}(\sigma_{A_{L},A_{R}}\log\sigma_{A_{L},A_{R}}),$ (12)

where $A_{L}$ and $A_{R}$ ($B_{L}$ and $B_{R}$) are the copies of the original
subsystem $A$ ($B$). Then the mutual information is upper-bounded by
$2\mathcal{E}$, i.e. $\mathcal{I}_{\rho_{\beta}}(A:B)\leq 2\mathcal{E}$ [67,
22]. We show the results of $\mathcal{E}$ in Fig. 3(b). $\mathcal{E}$ trends
upward for $\alpha<1$ and stabilizes for $\alpha>1$, supporting the thermal
area law’s optimality condition. While this discussion focuses solely on
Hamiltonian (10), it’s worth noting that similar behaviors are numerically
demonstrated even when the disorder is added into the Hamiltonian (see SM
[52]).

Thermal area law in the quantum entanglement.— The mutual information includes
both classical and quantum correlations. We now extract purely quantum
correlations from the Gibbs state to analyze the thermal area law condition.
To this end, we investigate bilinear fermionic systems (9) using Shapourian,
Shiozaki, and Ryu’s (SSR) logarithmic negativity [50] to quantify quantum
entanglement between subsystems $A$ and $B$ in a mixed state efficiently. The
SSR logarithmic negativity is formulated through the partial time-reversal
transform $R_{A}$, yielding:

$\displaystyle E_{\mathrm{SSR}}(\rho):=\log\|\rho^{R_{A}}\|_{1}\,,$ (13)

where $\|\cdot\|_{1}$ is the trace norm. Details on the SSR logarithmic
negativity for bilinear systems are in the SM [52]. We calculate this for 1D
and 2D Hamiltonian systems (9) with identical parameters and partitioning as
mutual information, showing results in Figs. 4 (a) for 1D and (b) for 2D.
These figures demonstrate that the SSR logarithmic negativity shows the same
behavior as the mutual information in Figs 2(c) and (d). Hence, the figures
indicate that under the same conditions, namely $\alpha>(D+1)/2$, quantum
entanglement likewise follows the thermal area law.

Figure 4: The SSR negativity $E_{\mathrm{SSR}}$ of long-range bilinear
fermions. (a) $E_{\mathrm{SSR}}$ of the 1D chain of $N$ sites and (b)
$E_{\mathrm{SSR}}/N$ of the 2D square lattice with a side length $N$ between
one half of the system and the other half for different $\alpha$’s. The
figures are averaged over the 1000 samples of random variables
$t_{i,j}\in[0,1]$ in Eq. (9) at $\beta=2$.

Proof of main theorem.— We here provide the outline of the proof of the main
theorem. For the sake of simplicity, we focus on the simple case where $h_{Z}$
is $h_{i,j}$ that acts only on two different points $i,j\in\Lambda$. The proof
for the general case is in the SM [52]. Let us consider the total Hamiltonian
$H=H_{A}+H_{B}+H_{\partial A}$, where $H_{A}$, $H_{B}$ are the Hamiltonian
supported only on $A,B$, respectively. The part $H_{\partial A}$ is the
interaction Hamiltonian between $A$ and $B$, given as $H_{\partial
A}=\sum_{i\in A}\sum_{j\in B}h_{i,j}$ with
$h_{i,j}=\sum_{s=1}^{d_{0}^{4}}h_{i}^{(s)}\otimes h_{j}^{(s)}$. Here, $d_{0}$
is the Hilbert dimension of the one site. Following Ref. [19], one has the
following inequality from the Gibbs variational principle:

$\displaystyle F(\rho)-F(\rho_{\beta})$
$\displaystyle=\beta^{-1}D(\rho||\rho_{\beta})\geq 0\,,$ (14)

where $F(\rho)$ is the nonequilibrium free energy defined as $F(\rho):={\rm
tr}(H\rho)+\beta^{-1}{\rm tr}[\rho\log(\rho)]$ and $D$ is the quantum relative
entropy defined as $D(\rho||\sigma)={\rm tr}(\rho\log\rho)-{\rm
tr}(\rho\log\sigma)$. Plugging $\rho=\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}$
into this inequality, one obtains the relation on the mutual information:

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)$
$\displaystyle\leq\beta\,\mathrm{tr}\left[(\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}-\rho_{\beta}^{AB})H_{\partial
A}\right]\,$ $\displaystyle=\beta\sum_{i\in A}\sum_{j\in
B}\sum_{s}\mathrm{Cor}_{\rho_{\beta}}(h_{i}^{(s)},h_{j}^{(s)})\,.$ (15)

The assumption of the clustering property leads to

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)\leq C\sum_{i\in A}\sum_{j\in
B}\sum_{s}\beta\|h_{i}^{(s)}\|\cdot\|h_{j}^{(s)}\|d_{i,j}^{-\alpha}.$ (16)

By the condition (3), we have

$\displaystyle\max_{s}\left(\|h_{i}^{(s)}\|,\|h_{j}^{(s)}\|\right)\leq
J_{i,j}\leq\frac{g}{d_{i,j}^{\alpha}}.$ (17)

Combining Eq. (16) and Eq. (17), we finally obtain

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)\leq\sum_{i\in A}\sum_{j\in
B}\frac{\beta d_{0}^{4}gC}{d_{i,j}^{2\alpha}}.$ (18)

It can be proven that the summation of the distance $\sum_{i\in A}\sum_{j\in
B}d_{i,j}^{-2\alpha}$ is upper bounded by the boundary area when
$2\alpha>D+1$. Therefore, Ineq. (18) leads to the desired inequality (7).

Summary and outlook.— We consider the validity of the thermal area law in the
systems with long-range interactions $r^{-\alpha}$. Under the assumption of
the clustering property, we derive the critical threshold $\alpha_{\rm
c}:=(D+1)/2$ above which any systems obey the thermal area law of the mutual
information (See Fig.1). Remarkably, the regime covers thermodynamically
nonextensive regime. This criteria potentially may allow for an efficient
representation of quantum Gibbs states exhibiting a power-law decay up to
$r^{-\alpha_{\rm c}}$. Given this criterion, it is a crucial future problem to
develop an efficiency-guaranteed algorithm for simulating long-range
interacting systems at finite temperatures, as well as constructing tensor
network states with polynomial bond dimensions [68, 69].

## References

  * Kitaev and Preskill [2006] A. Kitaev and J. Preskill, Topological Entanglement Entropy, Phys. Rev. Lett. 96, 110404 (2006).
  * Levin and Wen [2006] M. Levin and X.-G. Wen, Detecting Topological Order in a Ground State Wave Function, Phys. Rev. Lett. 96, 110405 (2006).
  * Chen  _et al._ [2010] X. Chen, Z.-C. Gu, and X.-G. Wen, Local unitary transformation, long-range quantum entanglement, wave function renormalization, and topological order, Phys. Rev. B 82, 155138 (2010).
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Rev. Mod. Phys. 82, 277 (2010).
  * Vidal  _et al._ [2003] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Entanglement in Quantum Critical Phenomena, Phys. Rev. Lett. 90, 227902 (2003).
  * Audenaert  _et al._ [2002] K. Audenaert, J. Eisert, M. B. Plenio, and R. F. Werner, Entanglement properties of the harmonic chain, Phys. Rev. A 66, 042327 (2002).
  * Plenio  _et al._ [2005] M. B. Plenio, J. Eisert, J. Dreißig, and M. Cramer, Entropy, Entanglement, and Area: Analytical Results for Harmonic Lattice Systems, Phys. Rev. Lett. 94, 060503 (2005).
  * Gioev and Klich [2006] D. Gioev and I. Klich, Entanglement Entropy of Fermions in Any Dimension and the Widom Conjecture, Phys. Rev. Lett. 96, 100503 (2006).
  * Hastings [2007] M. B. Hastings, An area law for one-dimensional quantum systems, Journal of Statistical Mechanics: Theory and Experiment 2007, P08024 (2007).
  * Masanes [2009] L. Masanes, Area law for the entropy of low-energy states, Phys. Rev. A 80, 052104 (2009).
  * Arad  _et al._ [2012] I. Arad, Z. Landau, and U. Vazirani, Improved one-dimensional area law for frustration-free systems, Phys. Rev. B 85, 195145 (2012).
  * Arad  _et al._ [2013] I. Arad, A. Kitaev, Z. Landau, and U. Vazirani, An area law and sub-exponential algorithm for 1D systems, arXiv preprint arXiv:1301.1162 (2013), arXiv:1301.1162 .
  * Brandão and Horodecki [2013] F. G. Brandão and M. Horodecki, An area law for entanglement from exponential decay of correlations, Nature Physics 9, 721 (2013).
  * Kuwahara and Saito [2020a] T. Kuwahara and K. Saito, Area law of noncritical ground states in 1D long-range interacting systems, Nature Communications 11, 4478 (2020a).
  * Anshu  _et al._ [2022] A. Anshu, I. Arad, and D. Gosset, An area law for 2d frustration-free spin systems, in  _Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing_, STOC 2022 (Association for Computing Machinery, New York, NY, USA, 2022) p. 12–18.
  * Landau  _et al._ [2015] Z. Landau, U. Vazirani, and T. Vidick, A polynomial time algorithm for the ground state of one-dimensional gapped local Hamiltonians, Nature Physics 11, 566 (2015).
  * Arad  _et al._ [2017] I. Arad, Z. Landau, U. Vazirani, and T. Vidick, Rigorous RG Algorithms and Area Laws for Low Energy Eigenstates in 1D, Communications in Mathematical Physics 356, 65 (2017).
  * Wolf  _et al._ [2008] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Area Laws in Quantum Systems: Mutual Information and Correlations, Phys. Rev. Lett. 100, 070502 (2008).
  * Sherman  _et al._ [2016] N. E. Sherman, T. Devakul, M. B. Hastings, and R. R. P. Singh, Nonzero-temperature entanglement negativity of quantum spin models: Area law, linked cluster expansions, and sudden death, Phys. Rev. E 93, 022128 (2016).
  * Berta  _et al._ [2018] M. Berta, F. G. S. L. Brandão, J. Haegeman, V. B. Scholz, and F. Verstraete, Thermal states as convex combinations of matrix product states, Phys. Rev. B 98, 235154 (2018).
  * Kuwahara  _et al._ [2021] T. Kuwahara, A. M. Alhambra, and A. Anshu, Improved Thermal Area Law and Quasilinear Time Algorithm for Quantum Gibbs States, Phys. Rev. X 11, 011047 (2021).
  * Lemm and Siebert [2023] M. Lemm and O. Siebert, Thermal Area Law for Lattice Bosons, Quantum 7, 1083 (2023).
  * Kuwahara and Saito [2022] T. Kuwahara and K. Saito, Exponential clustering of bipartite quantum entanglement at arbitrary temperatures, Physical Review X 12, 021022 (2022).
  * Defenu  _et al._ [2023] N. Defenu, T. Donner, T. Macrì, G. Pagano, S. Ruffo, and A. Trombettoni, Long-range interacting quantum systems, Rev. Mod. Phys. 95, 035002 (2023).
  * Yan  _et al._ [2013] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Observation of dipolar spin-exchange interactions with lattice-confined polar molecules, Nature 501, 521 (2013).
  * Jurcevic  _et al._ [2014] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt, and C. F. Roos, Quasiparticle engineering and entanglement propagation in a quantum many-body system, Nature 511, 202 (2014).
  * Richerme  _et al._ [2014] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C. Monroe, Non-local propagation of correlations in quantum systems with long-range interactions, Nature 511, 198 (2014).
  * Islam  _et al._ [2013] R. Islam, C. Senko, W. C. Campbell, S. Korenblit, J. Smith, A. Lee, E. E. Edwards, C.-C. J. Wang, J. K. Freericks, and C. Monroe, Emergence and Frustration of Magnetism with Variable-Range Interactions in a Quantum Simulator, Science 340, 583 (2013).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, _et al._ , Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017), article.
  * Zhang  _et al._ [2017] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator, Nature 551, 601 (2017).
  * Koffel  _et al._ [2012] T. Koffel, M. Lewenstein, and L. Tagliacozzo, Entanglement Entropy for the Long-Range Ising Chain in a Transverse Field, Phys. Rev. Lett. 109, 267203 (2012).
  * Vodola  _et al._ [2014] D. Vodola, L. Lepori, E. Ercolessi, A. V. Gorshkov, and G. Pupillo, Kitaev Chains with Long-Range Pairing, Phys. Rev. Lett. 113, 156402 (2014).
  * Gong  _et al._ [2016] Z.-X. Gong, M. F. Maghrebi, A. Hu, M. Foss-Feig, P. Richerme, C. Monroe, and A. V. Gorshkov, Kaleidoscope of quantum phases in a long-range interacting spin-1 chain, Phys. Rev. B 93, 205115 (2016).
  * Kozin and Kyriienko [2019] V. K. Kozin and O. Kyriienko, Quantum time crystals from hamiltonians with long-range interactions, Phys. Rev. Lett. 123, 210602 (2019).
  * Kuwahara and Saito [2020b] T. Kuwahara and K. Saito, Strictly Linear Light Cones in Long-Range Interacting Systems of Arbitrary Dimensions, Phys. Rev. X 10, 031010 (2020b).
  * Tran  _et al._ [2021a] M. C. Tran, A. Y. Guo, A. Deshpande, A. Lucas, and A. V. Gorshkov, Optimal State Transfer and Entanglement Generation in Power-Law Interacting Systems, Phys. Rev. X 11, 031016 (2021a).
  * Tran  _et al._ [2021b] M. C. Tran, A. Y. Guo, C. L. Baldwin, A. Ehrenberg, A. V. Gorshkov, and A. Lucas, Lieb-robinson light cone for power-law interactions, Phys. Rev. Lett. 127, 160401 (2021b).
  * Minato  _et al._ [2022] T. Minato, K. Sugimoto, T. Kuwahara, and K. Saito, Fate of measurement-induced phase transition in long-range interactions, Phys. Rev. Lett. 128, 010603 (2022).
  * Tamaki and Saito [2020] S. Tamaki and K. Saito, Energy current correlation in solvable long-range interacting systems, Physical Review E 101, 042118 (2020).
  * Dyson [1969] F. J. Dyson, Existence of a phase-transition in a one-dimensional ising ferromagnet, Communications in Mathematical Physics 12, 91 (1969).
  * Kosterlitz [1976] J. M. Kosterlitz, Phase transitions in long-range ferromagnetic chains, Phys. Rev. Lett. 37, 1577 (1976).
  * Araki [1975] H. Araki, On uniqueness of KMS states of one-dimensional quantum lattice systems, Communications in Mathematical Physics 44, 1 (1975).
  * Kishimoto [1976] A. Kishimoto, On uniqueness of KMS states of one-dimensional quantum lattice systems, Communications in Mathematical Physics 47, 167 (1976).
  * Gong  _et al._ [2017] Z.-X. Gong, M. Foss-Feig, F. G. S. L. Brandão, and A. V. Gorshkov, Entanglement Area Laws for Long-Range Interacting Systems, Phys. Rev. Lett. 119, 050501 (2017).
  * Kastoryano and Brandão [2016] M. J. Kastoryano and F. G. S. L. Brandão, Quantum Gibbs Samplers: The Commuting Case, Communications in Mathematical Physics 344, 915 (2016).
  * Brandão and Kastoryano [2019] F. G. S. L. Brandão and M. J. Kastoryano, Finite Correlation Length Implies Efficient Preparation of Quantum Thermal States, Communications in Mathematical Physics 365, 1 (2019).
  * Brandao and Cramer [2015] F. G. Brandao and M. Cramer, Equivalence of statistical mechanical ensembles for non-critical quantum systems, arXiv preprint arXiv:1502.03263 (2015), arXiv:1502.03263 .
  * Kuwahara and Saito [2020c] T. Kuwahara and K. Saito, Eigenstate Thermalization from the Clustering Property of Correlation, Phys. Rev. Lett. 124, 200604 (2020c).
  * Shapourian  _et al._ [2017] H. Shapourian, K. Shiozaki, and S. Ryu, Partial time-reversal transformation and entanglement negativity in fermionic systems, Phys. Rev. B 95, 165101 (2017).
  * Neyenhuis  _et al._ [2017] B. Neyenhuis, J. Zhang, P. W. Hess, J. Smith, A. C. Lee, P. Richerme, Z.-X. Gong, A. V. Gorshkov, and C. Monroe, Observation of prethermalization in long-range interacting spin chains, Science advances 3, e1700672 (2017).
  * [52] Supplementary material .
  * Araki [1969] H. Araki, Gibbs states of a one dimensional quantum lattice, Communications in Mathematical Physics 14, 120 (1969).
  * Kimura and Kuwahara [2024] Y. Kimura and T. Kuwahara, Clustering theorem in 1d long-range interacting systems at arbitrary temperatures (2024), arXiv:2403.11431 [quant-ph] .
  * Latorre and Riera [2009] J. I. Latorre and A. Riera, A short review on entanglement in quantum spin systems, Journal of Physics A: Mathematical and Theoretical 42, 504002 (2009).
  * Sen  _et al._ [2016] A. Sen, S. Nandy, and K. Sengupta, Entanglement generation in periodically driven integrable systems: Dynamical phase transitions and steady state, Phys. Rev. B 94, 214301 (2016).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Žnidarič  _et al._ [2008] M. Žnidarič, T. c. v. Prosen, and I. Pižorn, Complexity of thermal states in quantum spin chains, Phys. Rev. A 78, 022103 (2008).
  * Bernigau  _et al._ [2015] H. Bernigau, M. J. Kastoryano, and J. Eisert, Mutual information area laws for thermal free fermions, Journal of Statistical Mechanics: Theory and Experiment 2015, P02008 (2015).
  * Barthel [2017] T. Barthel, One-dimensional quantum systems at finite temperatures can be simulated efficiently on classical computers, arXiv preprint arXiv:1708.09349 (2017), arXiv:1708.09349 .
  * Chen  _et al._ [2018] B.-B. Chen, L. Chen, Z. Chen, W. Li, and A. Weichselbaum, Exponential Thermal Tensor Network Approach for Quantum Lattice Models, Phys. Rev. X 8, 031082 (2018).
  * Li  _et al._ [2019] H. Li, B.-B. Chen, Z. Chen, J. von Delft, A. Weichselbaum, and W. Li, Thermal tensor renormalization group simulations of square-lattice quantum spin models, Phys. Rev. B 100, 045110 (2019).
  * Weichselbaum [2012] A. Weichselbaum, Non-abelian symmetries in tensor networks: A quantum symmetry space approach, Annals of Physics 327, 2972 (2012).
  * Weichselbaum [2020] A. Weichselbaum, X-symbols for non-Abelian symmetries in tensor networks, Phys. Rev. Research 2, 023385 (2020).
  * ft [0] We implement $\mathrm{SU}(2)$ spin symmetry of Eq. (10) with the tensor library QSpace. We set the bond dimensions of the matrix product operator (MPO) to 210 for Fig. 3(a) and 120 for Fig. 3(b). These bond dimensions are associated with the 210 and 120 $\mathrm{SU}(2)$ multiplets, which are roughly equivalent to $850$ and $450$ states.
  * Cottrell  _et al._ [2019] W. Cottrell, B. Freivogel, D. M. Hofman, and S. F. Lokhande, How to build the thermofield double state, Journal of High Energy Physics 2019, 58 (2019).
  * Bagchi and Pati [2015] S. Bagchi and A. K. Pati, Monogamy, polygamy, and other properties of entanglement of purification, Phys. Rev. A 91, 042323 (2015).
  * Molnar  _et al._ [2015] A. Molnar, N. Schuch, F. Verstraete, and J. I. Cirac, Approximating Gibbs states of local Hamiltonians efficiently with projected entangled pair states, Phys. Rev. B 91, 045138 (2015).
  * Kliesch  _et al._ [2014] M. Kliesch, C. Gogolin, M. J. Kastoryano, A. Riera, and J. Eisert, Locality of Temperature, Phys. Rev. X 4, 031019 (2014).
  * Kuwahara  _et al._ [2020] T. Kuwahara, K. Kato, and F. G. S. L. Brandão, Clustering of Conditional Mutual Information for Quantum Gibbs States above a Threshold Temperature, Phys. Rev. Lett. 124, 220601 (2020).
  * Kuwahara and Saito [2020d] T. Kuwahara and K. Saito, Gaussian concentration bound and Ensemble equivalence in generic quantum many-body systems including long-range interactions, Annals of Physics 421, 168278 (2020d).

## acknowledgments

D. K. thanks M. Kim for helpful discussions about QSpace. D. K. and T. K.
acknowledge the Hakubi projects of RIKEN. T. K. was supported by Japan Science
and Technology Agency Precursory Research for Embryonic Science and Technology
(Grant No. JPMJPR2116). K. S. was supported by JSPS Grants-in-Aid for
Scientific Research (No. JP19H05603, No. JP19H05791 and No. JP23H01099).

Supplementary Material for “Thermal Area Law in Long-Range Interacting
Systems”
Donghoon Kim1, Tomotaka Kuwahara1,2, Keiji Saito3

1Analytical quantum complexity RIKEN Hakubi Research Team, RIKEN Cluster for
Pioneering Research (CPR)/ RIKEN Center for Quantum Computing (RQC), Wako,
Saitama 351-0198, Japan
2PRESTO, Japan Science and Technology (JST), Kawaguchi, Saitama 332-0012,
Japan
3Department of Physics, Kyoto University, Kyoto 606-8502, Japan

###### Contents

  1. S.I Setup
  2. S.II Proof of Main Theorem in a General Setting
  3. S.III Power-law clustering theorem for $\alpha>D$
  4. S.IV Exact Diagonalization Results for Long-Range Heisenberg Chains
  5. S.V Numerical Detail of Shapourian-Shiozaki-Ryu Negativity
    1. S.V.1 Numerical Procedure

##  S.I Setup

We consider a quantum system, where we use $\Lambda$ to represent the set of
sites. Given an arbitrary subset $X$ in $\Lambda$, $X\subseteq\Lambda$, we
denote the cardinality of the set $X$, i.e., the number of sites contained in
$X$, by $|X|$. For any two subsets $X,Y$ in $\Lambda$, $X,Y\subseteq\Lambda$,
we define the distance between the two sets, $d_{X,Y}$, as the shortest path
length of the graph that connects $X$ and $Y$. We remark that when $X$
intersects with $Y$, $X\cap Y\neq\emptyset$, the distance between the sets $X$
and $Y$ is zero: $d_{X,Y}=0$. We denote the complementary set of $X$ as
$X^{\rm c}$, namely $X^{\rm c}:=\Lambda\setminus X$, and the surface subset of
$X$ as $\partial X$, i.e. $\partial X:=\\{i\in X|d_{i,X^{\rm c}}=1\\}$,
respectively. Also, for an arbitrary operator $O$, we denote the support of
$O$ by ${\rm Supp}(O)$. We often describe the support of the operator
explicitly by adding an Index of lower right subscript such as $O_{X}$ with
${\rm Supp}(O_{X})=X$.

We consider a $k$-local Hamiltonian $H$ with long-range interactions:

$\displaystyle H=\sum_{Z:|Z|\leq k}h_{Z}\,,$ (S.0)

where $Z$ is a subset indicating the interacting sites, and $h_{Z}$ is the
local Hamiltonian. Note that $|Z|$ is the number of the interacting sites, and
hence $h_{Z}$ means $Z$-body interaction. We assume that the interaction form
satisfies the following property

$\displaystyle
J_{i,i^{\prime}}:=\sum_{Z:Z\ni\\{i,i^{\prime}\\}}\left\|h_{Z}\right\|\leq{g\over(1+d_{i,i^{\prime}})^{\alpha}}.$
(S.1)

Here, $J_{i,i^{\prime}}$ is the maximum operator norm of the local
Hamiltonians containing the sites $i$ and $i^{\prime}$.

We consider the quantum Gibbs state $\rho_{\beta}$ with a fixed inverse
temperature $\beta$ defined as follows:

$\displaystyle\rho_{\beta}:=e^{-\beta H}/Z,\quad Z={\rm tr}(e^{-\beta H}).$
(S.2)

The mutual information of $\rho_{\beta}$ between two subsets $A$ and $B$ is
defined as

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B):=S(\rho_{\beta}^{A})+S(\rho_{\beta}^{B})-S(\rho_{\beta}^{AB}),$
(S.3)

where $\rho_{\beta}^{X}$ is the density matrix of $\rho_{\beta}$ for the
regime $X\in\\{A,B,AB\\}$ and $S(\rho_{\beta}^{X}):=-{\rm
tr}(\rho_{\beta}^{X}\log\rho_{\beta}^{X})$ is the von Neumann entropy. The
standard correlation function of the density matrix $\rho$ between two
observables $O_{A}$ and $O_{B}$ is defined as

$\displaystyle{\rm Cor}_{\rho}(O_{A},O_{B}):={\rm tr}(\rho O_{A}O_{B})-{\rm
tr}(\rho O_{A})\cdot{\rm tr}(\rho O_{B}).$ (S.4)

##  S.II Proof of Main Theorem in a General Setting

######  Theorem 1.

Let us assume that the correlation function of the quantum Gibbs state
$\rho_{\beta}$ between two arbitrary operators $O_{X}$ and $O_{Y}$ supported
on the subsets $X$ and $Y$ in $\Lambda$ satisfies the following power-law
clustering property:

$\displaystyle{\rm
Cor}_{\rho_{\beta}}(O_{X},O_{Y})\leq\frac{C}{d_{X,Y}^{\alpha}}\|O_{X}\|\cdot\|O_{Y}\|,$
(S.5)

where $C$ is an $\mathcal{O}(1)$ constant. Then, for the bipartition $A,B$ of
$D$-dimensional lattice $\Lambda$ ($A\cup B=\Lambda$) and $\alpha>(D+1)/2$,
the mutual information is upper bounded by

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)\leq\beta\cdot C\cdot|\partial
A|\cdot\text{const}.$ (S.6)

###### Proof.

We partition the lattice $\Lambda$ into two subsets $A,B\subset\Lambda$ and
write the Hamiltonian as $H=H_{A}+H_{B}+H_{\partial A}$. Here, $H_{A}$,
$H_{B}$ are the Hamiltonian supported only on $A,B$, respectively, and
$H_{\partial A}$ is the interaction between $A$ and $B$, i.e.,

$\displaystyle H_{\partial A}=\sum_{\begin{subarray}{c}Z:Z\cap
A\neq\emptyset\\\ Z\cap B\neq\emptyset\end{subarray}}h_{Z}.$ (S.7)

Using the non-negativity of the quantum relative entropy
$D(\rho\|\sigma):={\rm tr}(\rho(\log\rho-\log\sigma))$ yields

$\displaystyle D(\rho\|\rho_{\beta})=\beta{\rm tr}(\rho H)-S(\rho)+\log
Z=\beta F_{\beta}(\rho)-\beta F_{\beta}(\rho_{\beta})\geq 0,$ (S.8)

for an arbitrary density matrix $\rho$. Here, $F_{\beta}(\rho)$ is the free
energy defined as

$\displaystyle F_{\beta}(\rho):={\rm tr}(\rho H)-\beta^{-1}S(\rho),$ (S.9)

and it becomes $F_{\beta}(\rho_{\beta})=-\beta^{-1}\log Z$ at equilibrium.
Substituting $\rho=\rho_{\beta}^{A}\otimes\rho_{\beta}^{B}$ and using
$\mathcal{I}_{\rho_{\beta}}(A:B)=S(\rho_{\beta}^{A}\otimes\rho_{\beta}^{B})-S(\rho_{\beta}^{AB})$
and ${\rm tr}((H_{A}+H_{B})\rho_{\beta}^{A}\otimes\rho_{\beta}^{B})={\rm
tr}((H_{A}+H_{B})\rho_{\beta}^{AB})$, we obtain

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)\leq\beta\cdot{\rm
tr}((\rho_{A}\otimes\rho_{B}-\rho_{AB})H_{\partial A}).$ (S.10)

For the $k$-local Hamiltonian, each local term $h_{Z}$ can be expanded by sum
of the tensor products of two operators supported on $Z_{A}=Z\cap A$ and
$Z_{B}=Z\cap B$:

$\displaystyle h_{Z}=\sum_{s=1}^{d_{0}^{2k}}h_{Z_{A}}^{(s)}\otimes
h_{Z_{B}}^{(s)}.$ (S.11)

Equations (S.7), (S.10) and (S.11) lead to

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)\leq\sum_{\begin{subarray}{c}Z:Z\cap
A\neq\emptyset\\\ Z\cap
B\neq\emptyset\end{subarray}}\sum_{s=1}^{d_{0}^{2k}}\beta\cdot{\rm
tr}((\rho_{A}\otimes\rho_{B}-\rho_{AB})h_{Z_{A}}^{(s)}\otimes
h_{Z_{B}}^{(s)})=\sum_{\begin{subarray}{c}Z:Z\cap A\neq\emptyset\\\ Z\cap
B\neq\emptyset\end{subarray}}\sum_{s=1}^{d_{0}^{2k}}\beta\cdot\mathrm{Cor}_{\rho_{AB}}(h_{Z_{A}}^{(s)},h_{Z_{B}}^{(s)}).$
(S.12)

The assumption of the clustering property gives

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)$
$\displaystyle\leq\sum_{\begin{subarray}{c}Z:Z\cap A\neq\emptyset\\\ Z\cap
B\neq\emptyset\end{subarray}}\sum_{s=1}^{d_{0}^{2k}}\beta\cdot\|h_{Z_{A}}^{(s)}\|\|h_{Z_{B}}^{(s)}\|\frac{C}{d_{Z_{A},Z_{B}}^{\alpha}}$
(S.13) $\displaystyle\leq\sum_{i\in A}\sum_{j\in
B}\sum_{\begin{subarray}{c}Z:Z_{A}\ni i\\\ Z_{B}\ni
j\end{subarray}}\sum_{s=1}^{d_{0}^{2k}}\beta\cdot\|h_{Z_{A}}^{(s)}\|\|h_{Z_{B}}^{(s)}\|\frac{C}{d_{i,j}^{\alpha}}$
(S.14) $\displaystyle\leq\sum_{i\in A}\sum_{j\in
B}\sum_{\begin{subarray}{c}Z:Z_{A}\ni i\\\ Z_{B}\ni j\end{subarray}}\beta\cdot
d_{0}^{2k}\max_{s}\|h_{Z_{A}}^{(s)}\|\|h_{Z_{B}}^{(s)}\|\frac{C}{d_{i,j}^{\alpha}}.$
(S.15)

In the second inequality, we have used the fact that
$d_{Z_{A},Z_{B}}^{-\alpha}$ can be expressed as $d_{i,j}^{-\alpha}$ for some
$i\in Z_{A}$ and $j\in Z_{B}$. In the third inequality, we has used
$\sum_{s=1}^{d_{0}^{2k}}\|h_{Z_{A}}^{(s)}\|\|h_{Z_{B}}^{(s)}\|\leq
d_{0}^{2k}\max_{s}\|h_{Z_{A}}^{(s)}\|\|h_{Z_{B}}^{(s)}\|$. By the condition
Eq. (S.1) of the long-range interaction, we obtain

$\displaystyle\sum_{\begin{subarray}{c}Z:Z_{A}\ni i\\\ Z_{B}\ni
j\end{subarray}}\max_{s\in\mathbb{Z}_{[1,d_{0}^{2k}]}}\|h_{Z_{A}}^{(s)}\|\|h_{Z_{B}}^{(s)}\|\leq\sum_{Z:Z\ni\\{i,j\\}}\|h_{Z}\|\leq\frac{g}{d_{i,j}^{\alpha}}.$
(S.16)

Combining Eq. (S.15) and Eq. (S.16),

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)\leq\sum_{i\in A}\sum_{j\in
B}\frac{\beta d_{0}^{2k}gC}{d_{i,j}^{2\alpha}}.$ (S.17)

It was proven that $\sum_{i\in A}\sum_{j\in B}d_{i,j}^{-2\alpha}$ is upper
bounded by the boundary area $|\partial A|$ when $2\alpha>D+1$ [39]. With Eq.
(S.17), we finally obtain

$\displaystyle\mathcal{I}_{\rho_{\beta}}(A:B)\leq\beta\cdot C\cdot|\partial
A|\cdot\text{const},$ (S.18)

for $2\alpha>D+1$. ∎

##  S.III Power-law clustering theorem for $\alpha>D$

We consider the regime $\alpha>D$, so that the extensivity property holds as

$\displaystyle\left\|\widehat{H_{L}}\right\|\leq
g|L|,\quad\widehat{H_{L}}=\sum_{Z:Z\cap L\neq\emptyset}h_{Z}=H-H_{L^{\rm c}}$
(S.19)

for an arbitrary subset $L\subseteq\Lambda$. Then, the clustering theorem for
the bipartite correlations is rigorously proven as

######  Theorem 2.

The following relation holds for long-range interacting systems for the
temperatures above $\beta_{c}^{-1}$:

$\displaystyle{\rm Cor}_{\rho_{\beta}}(O_{X},O_{Y})$ $\displaystyle\leq
C\|O_{X}\|\cdot\|O_{Y}\|{|X||Y|\,e^{{(|X|+|Y|)/k}}\over d_{X,Y}^{\alpha}}\,,$
(S.20)

where $C$ is an $\mathcal{O}(1)$ constant and the threshold temperature is
given by $\beta_{c}=W(1/ue)/(2gk)$ with the Lambert $W$ function (i.e., the
inverse function of $xe^{x}$) and
$u=\sum_{j\in\Lambda}1/(1+d_{i,j})^{\alpha}$. We can lower-bound $\beta_{c}$
as $\beta_{c}\geq 1/(8ugk)$ by using

$W(1/ue)\leq 1/(1+ue)\leq 1/(4u),$

where we use $u$ is always larger than $1$ from the definition.

For the proof of the clustering theorem, we utilize the cluster expansion [70,
71, 69]. First of all, we describe the correlation function ${\rm
Cor}_{\rho}(O_{X},O_{Y})$ can be written as

$\displaystyle{\rm Cor}_{\rho_{\beta}}(O_{X},O_{Y})=\frac{1}{{\cal Z}^{2}}{\rm
tr}\left(e^{-\beta H^{(+)}}O_{X}^{(0)}O_{Y}^{(1)}\right),$ (S.21)

where we define the operators $O^{(+)}$, $O^{(0)}$ and $O^{(1)}$ as

$\displaystyle O^{(+)}=O\otimes\hat{1}+\hat{1}\otimes O,\quad
O^{(0)}=O\otimes\hat{1},\quad O^{(1)}=O\otimes\hat{1}-\hat{1}\otimes O\,,$
(S.22)

by considering a copy of the original Hilbert space. The quantity ${\cal Z}$
is a normalization factor ${\cal Z}^{2}={\rm tr}(e^{-\beta H^{(+)}})$($=[{\rm
tr}(e^{-\beta H})]^{2}$ in the original Hilbert space). It is simple to
observe that the following lemma holds:

######  Lemma 1.

For an arbitrary set of operators $\\{O_{Z_{j}}\\}_{j=1}^{m}$, the following
relation holds:

$\displaystyle{\rm tr}\left(O_{Z_{1}}^{(+)}O_{Z_{2}}^{(+)}\cdots
O_{Z_{m}}^{(+)}O_{X}^{(0)}O_{Y}^{(1)}\right)=0$ (S.23)

when the subsets $\\{Z_{1},Z_{2},\ldots,Z_{m}\\}$ cannot connect the subsets
$X$ and $Y$ [see Fig. 5 (a)], i.e.,

$\displaystyle\left(X\cup Z_{i_{1}}\cup Z_{i_{2}}\cup\cdots\cup
Z_{i_{s}}\right)\cap\left(Y\cup Z_{i_{s+1}}\cup Z_{i_{s+2}}\cup\cdots\cup
Z_{i_{m}}\right)=\emptyset$ (S.24)

with $\\{i_{1},i_{2},\ldots,i_{m}\\}=\\{1,2,\ldots,m\\}$.

Proof of Lemma 1. Without loss of generality, we let
$\\{i_{1},i_{2},\ldots,i_{s}\\}=\\{1,2,\ldots,s\\}$ and
$\\{1,2,\ldots,s\\}\setminus\\{i_{1},i_{2},\ldots,i_{m}\\}=\\{s+1,s+2,\ldots,m\\}$.
We then denote

$\displaystyle O_{Z_{1}}^{(+)}O_{Z_{2}}^{(+)}\cdots
O_{Z_{s}}^{(+)}=W_{1},\quad O_{Z_{s+1}}^{(+)}O_{Z_{s+2}}^{(+)}\cdots
O_{Z_{m}}^{(+)}=W_{2},$ (S.25)

with

$\displaystyle{\rm Supp}(X\cup W_{1})\cap{\rm Supp}(Y\cup W_{2})=\emptyset,$
(S.26)

where the second equation is equivalent to the condition (S.24). We then
decompose as

${\rm tr}\left(O^{(+)}_{Z_{1}}O^{(+)}_{Z_{2}}\ldots
O^{(+)}_{Z_{m}}O^{(0)}_{X}O^{(1)}_{Y}\right)={\rm
tr}\left(W_{1}O^{(0)}_{X}\right){\rm tr}\left(W_{2}O^{(1)}_{Y}\right).$ (S.27)

Here, the operators
$O^{(+)}_{Z_{i_{k}}}=O_{Z_{i_{k}}}\otimes\hat{1}+\hat{1}\otimes O_{Z_{i_{k}}}$
are symmetric for the exchange of the two Hilbert space while the operator
$O^{(1)}_{Y}=O_{Y}\otimes\hat{1}-\hat{1}\otimes O_{Y}$ is antisymmetric, which
yields ${\rm tr}\left(W_{2}O^{(1)}_{Y}\right)=0$. We thus prove the main
equation (S.23). $\square$

[ End of Proof of Lemma 1]

---
Figure 5: (a): Schematic picture of unconnected subsets. (b): Connected
subsets. (c): Decomposition of the shortest subsets and the rest of the
subsets.

By using Lemma 1, the function ${\rm Cor}_{\rho_{\beta}}(O_{X},O_{Y})$ can be
expanded as

$\displaystyle{\rm Cor}_{\rho_{\beta}}=\frac{1}{{\cal
Z}^{2}}\sum_{m=0}^{\infty}\sum_{Z_{1},Z_{2},\ldots,Z_{m}:{\rm
connected}}\frac{(-\beta)^{m}}{m!}{\rm
tr}\left(h_{Z_{1}}^{(+)}h_{Z_{2}}^{(+)}\cdots
h_{Z_{m}}^{(+)}O_{X}^{(0)}O_{Y}^{(1)}\right),$ (S.28)

where the summation $\sum_{Z_{1},Z_{2},\ldots,Z_{m}:{\rm connected}}$ is taken
over the collections $\\{X,Z_{1},Z_{2},\ldots,Z_{m},Y\\}$ that connect $X$ and
$Y$. See Fig. 5 (b). Note that $\\{Z_{s}\\}_{s=1}^{m}$ can be the same, e.g.,
$Z_{1}=Z_{2}$. Let us define the operator $\rho_{\rm cl}$ as

$\displaystyle\rho_{{\rm cl}}=\frac{1}{{\cal
Z}^{2}}\sum_{m=0}^{\infty}\sum_{Z_{1},Z_{2},\ldots,Z_{m}:{\rm
connected}}\frac{(-\beta)^{m}}{m!}h_{Z_{1}}^{(+)}h_{Z_{2}}^{(+)}\cdots
h_{Z_{m}}^{(+)}.$ (S.29)

Using the operator $\rho_{{\rm cl}}$, the correlation ${\rm
Cor}_{\rho_{\beta}}(O_{X},O_{Y})$ has the following upper bound:

$\displaystyle{\rm Cor}_{\rho_{\beta}}(O_{X},O_{Y})\leq
2\left\|O_{X}\right\|\cdot\left\|O_{Y}\right\|\cdot\left\|\rho_{{\rm
cl}}\right\|_{1}$ (S.30)

To count the summation of $\sum_{Z_{1},Z_{2},\ldots,Z_{m}:{\rm connected}}$,
we first decompose as

$\displaystyle\\{Z_{1},Z_{2},\ldots,Z_{m}\\}=w_{\rm cl}\oplus w_{\rm cl}^{\rm
c},\quad w_{\rm cl}^{\rm c}:=\\{Z_{1},Z_{2},\ldots,Z_{m}\\}\setminus w_{\rm
cl},$ (S.31)

where $w_{\rm cl}$ is taken such that $w_{\rm cl}\oplus\\{X,Y\\}$ are
connected to each other (See Fig. 5 (b)), in other words,

$\displaystyle Z\cap Z^{\prime}=\emptyset\quad\textrm{for}\quad Z\in w_{\rm
cl},\quad Z^{\prime}\in w^{\rm c}_{\rm cl}.$ (S.32)

We define all the sets of $w_{\rm cl}$ and $w_{\rm cl}^{\rm c}$ as
$\mathcal{G}_{\rm cl}$ and $\mathcal{G}^{\rm c}_{\rm cl}$, respectively. We
also define the total support of the $w_{\rm cl}$ and $w_{\rm cl}^{\rm c}$ as
$V_{w_{\rm cl}}$ and $V_{w_{\rm cl}^{\rm c}}$, respectively. We then obtain

$\displaystyle\rho_{\rm cl}$ $\displaystyle={1\over{\cal
Z}^{2}}\sum_{m=0}^{\infty}{(-\beta)^{m}\over
m!}\sum_{s=0}^{m}{m!\over(m-s)!s!}\sum_{\\{Z_{1},\cdots,Z_{s}\\}=w_{\rm
cl}\in\mathcal{G}_{\rm
cl}}\\!\\!\\!\\!\tilde{h}(Z_{1},\cdots,Z_{s})\sum_{w_{\rm cl}^{\rm
c}\in\mathcal{G}^{\rm c}_{\rm cl}}\sum_{Z_{s+1},\cdots,Z_{m}\atop\in w_{\rm
cl}^{\rm c}}\\!\\!\\!\\!\\!\\!h^{(+)}(Z_{s+1})h^{(+)}(Z_{s+2})\cdots
h^{(+)}(Z_{m})$ $\displaystyle={1\over{\cal
Z}^{2}}\sum_{s=0}^{\infty}{(-\beta)^{s}\over s!}\sum_{w_{\rm
cl}\in\mathcal{G}_{\rm cl}}\tilde{h}(Z_{1},\cdots,Z_{s})e^{-\beta
H^{(+)}_{V_{w^{\rm c}_{\rm cl}}}}.$ (S.33)

Here, we have defined

$\displaystyle\tilde{h}(Z_{1},\cdots,Z_{s})$
$\displaystyle=\sum_{P}h^{(+)}(Z_{P_{1}})\cdots h^{(+)}(Z_{P_{s}})\,,$ (S.34)

where $P$ implies taking all combinations of different elements for $Z$. In
particular, when $Z$ are all different, it means taking a permutation. From
the expression (S.33), we obtain

$\displaystyle\|{\rho_{{\rm cl}}}\|_{1}$ $\displaystyle\leq{1\over{\cal
Z}^{2}}\sum_{s=0}^{\infty}{\beta^{s}\over s!}\sum_{w_{\rm
cl}\in\mathcal{G}_{\rm cl}}\|\tilde{h}(Z_{1},\cdots,Z_{s})\|\,{\rm
tr}\Bigl{[}e^{-\beta H^{(+)}_{V_{w^{\rm c}_{\rm
cl}}}}\Bigr{]}\leq\sum_{s=0}^{\infty}{\beta^{s}\over s!}\sum_{w_{\rm
cl}\in\mathcal{G}_{\rm cl}}e^{2\beta g|V_{w_{\rm
cl}}|}\|\tilde{h}(Z_{1},\cdots,Z_{s})\|\,.$ (S.35)

Here, we use the Golden-Thompson inequality as follows:

$\displaystyle{\rm tr}\Bigl{[}e^{-\beta H^{(+)}_{V_{w^{\rm c}_{\rm
cl}}}}\Bigr{]}$ $\displaystyle={\rm
tr}\Bigl{[}e^{-\beta\left(H^{(+)}+H^{(+)}_{V_{w^{\rm c}_{\rm
cl}}}-H^{(+)}\right)}\Bigr{]}\leq{\rm tr}\Bigl{[}e^{-\beta
H^{(+)}}e^{-\beta\left(H^{(+)}_{V_{w^{\rm c}_{\rm
cl}}}-H^{(+)}\right)}\Bigr{]}\leq{\cal Z}^{2}e^{2\beta\|H_{V^{\rm c}_{w_{\rm
cl}}}-H\|}$ $\displaystyle\leq{\cal Z}^{2}e^{2\beta g|V_{w_{\rm
cl}}|}\leq{\cal Z}^{2}e^{2\beta gks}$ (S.36)

where we use the condition (S.19) in the last inequality. From the
$k$-locality of the Hamiltonian, we have $|V_{w_{\rm cl}}|\leq k|w_{\rm
cl}|=ks$.

We divide each element in the set $w_{\rm cl}$ into the two subsets $w_{\rm
cl}^{\rm s}$ and $w_{\rm cl}^{\rm r}$, where $w_{\rm cl}^{\rm s}$ is a set
representing the shortest path with the length $s^{\prime}$ and $w_{\rm
cl}^{\rm r}$ is a set for the rest part with the length $(s-s^{\prime})$. See
Fig. 5 (c). We thus obtain

$\displaystyle\|\rho_{{\rm cl}}\|_{1}$
$\displaystyle\leq\sum_{s=0}^{\infty}{(\beta e^{2\beta gk})^{s}\over
s!}\sum_{s^{\prime}=1}^{s}{s!\over
s^{\prime}!(s-s^{\prime})!}\sum_{\\{Z_{1},\cdots,Z_{s^{\prime}}\\}=w_{\rm
cl}^{\rm
s}}\|\tilde{h}(Z_{1},\cdots,Z_{s^{\prime}})\|^{\prime}\sum_{\\{Z_{s^{\prime}+1},\cdots,Z_{s}\\}=w_{\rm
cl}^{\rm r}}\|\tilde{h}(Z_{s^{\prime}+1},\cdots,Z_{s})\|^{\prime}$
$\displaystyle=\sum_{s=0}^{\infty}(\beta e^{2\beta
gk})^{s}\sum_{s^{\prime}=1}^{s}\sum_{w_{\rm cl}^{\rm
s}}{\|\tilde{h}(Z_{1},\cdots,Z_{s^{\prime}})\|^{\prime}\over
s^{\prime}!}\sum_{w_{\rm cl}^{\rm
r}}{\|\tilde{h}(Z_{s^{\prime}+1},\cdots,Z_{s})\|^{\prime}\over(s-s^{\prime})!}\,,$
(S.37)

where

$\displaystyle\|\tilde{h}(Z_{1},\cdots,Z_{s^{\prime}})\|^{\prime}$
$\displaystyle:=\sum_{P}\|h_{Z_{P_{1}}}^{(+)}\|\cdots\|h_{Z_{P_{s^{\prime}}}}^{(+)}\|\leq
2^{s^{\prime}}\sum_{P}\|h_{Z_{P_{1}}}\|\cdots\|h_{Z_{P_{s^{\prime}}}}\|\,.$
(S.38)

Similarly to (S.38), $\|\tilde{h}(Z_{s^{\prime}+1},\cdots,Z_{s})\|^{\prime}$
is defined. Below, we evaluate the terms on $w_{\rm cl}^{\rm s}$ and $w_{\rm
cl}^{\rm r}$ in (S.37). To this end, the following lemma is useful.

######  Lemma 2.

For $\alpha>D$, the quantity $J_{i,j}$ defined in (S.1) satisfies the
following inequality

$\displaystyle\sum_{j\in\Lambda}J_{i,j}J_{j,k}\leq{g^{2}u\over(1+d_{i,k})^{\alpha}}\,,$
(S.39)

where $u=2^{\alpha}\sum_{j\in\Lambda}(1+d_{i,j})^{-\alpha}$. The iterative use
of (S.39) leads to
$\left[{\bm{J}}^{\ell}\right]_{i,k}\leq{g^{\ell}u^{\ell-1}/(1+d_{i,k})^{\alpha}}$.

Proof of Lemma 2.

$\displaystyle\sum_{j\in\Lambda}J_{i,j}J_{j,k}$
$\displaystyle=g^{2}\sum_{j\in\Lambda}{1\over(1+d_{i,j})^{\alpha}}{1\over(1+d_{j,k})^{\alpha}}\leq{g^{2}\over(1+d_{i,k})^{\alpha}}\sum_{j\in\Lambda}{(2+d_{i,j}+d_{j,k})^{\alpha}\over(1+d_{i,j})^{\alpha}}{1\over(1+d_{j,k})^{\alpha}}.$
(S.40)

Note $(x+y)^{\alpha}\leq 2^{\alpha-1}(x^{\alpha}+y^{\alpha})$ for $\alpha>1$.
Then we have

$\displaystyle\sum_{j\in\Lambda}J_{i,j}J_{j,k}$
$\displaystyle\leq{g^{2}2^{\alpha-1}\over(1+d_{i,k})^{\alpha}}\sum_{j\in\Lambda}{1\over(1+d_{i,j})^{\alpha}}+{1\over(1+d_{j,k})^{\alpha}}={g^{2}u\over(1+d_{i,k})^{\alpha}}.$
(S.41)

The finiteness of $u$ is guaranteed by the condition $\alpha>D$.

[ End of Proof of Lemma 2]

Let us first consider the terms on $w_{\rm cl}^{\rm s}$ in (S.37).

$\displaystyle\sum_{w_{\rm cl}^{\rm
s}}{\|\tilde{h}(Z_{1},\cdots,Z_{s^{\prime}})\|^{\prime}\over s^{\prime}!}$
$\displaystyle\leq{2^{s^{\prime}}\over s^{\prime}!}\sum_{i_{1}\in
X}\sum_{i_{2}\in\Lambda}\sum_{Z_{1}\ni\\{i_{1},i_{2}\\}}\sum_{i_{3}\in\Lambda}\sum_{Z_{2}\ni\\{i_{2},i_{3}\\}}\cdots\sum_{i_{s^{\prime}}\in\Lambda}\sum_{i_{s^{\prime}+1}\in
Y}\sum_{Z_{s^{\prime}}\ni\\{i_{s^{\prime}},i_{s^{\prime}+1}\\}}\|h_{Z_{1}}\|\cdots\|h_{Z_{s^{\prime}}}\|$
$\displaystyle={2^{s^{\prime}}\over s^{\prime}!}\sum_{i_{1}\in
X}\sum_{i_{2}\in\Lambda}\cdots\sum_{i_{s^{\prime}}\in\Lambda}\sum_{i_{s^{\prime}+1}\in
Y}J_{i_{1},i_{2}}\cdots
J_{i_{s^{\prime}},i_{s^{\prime}+1}}\leq{u^{-1}(2gu)^{s^{\prime}}\over
s^{\prime}!}|X|\,|Y|{1\over(1+d_{X,Y})^{\alpha}}\,.$ (S.42)

We next consider the terms of $w_{\rm cl}^{\rm r}$ in (S.37). We define
$n:=s-s^{\prime}$ and $L:=X\cup Y\cup\omega_{\rm cl}^{\rm s}$ for a
simplicity.

$\displaystyle{1\over n!}\sum_{w_{\rm cl}^{\rm
r}}\|\tilde{h}(Z_{s^{\prime}+1}\cdots,Z_{s})\|^{\prime}\leq{2^{n}\over
n!}\sum_{i_{1}\in L}\sum_{Z_{1}\ni i_{1}}\sum_{i_{2}\in L\cup
Z_{1}}\sum_{Z_{2}\ni i_{2}}\sum_{i_{3}\in L\cup Z_{1}\cup Z_{2}}\sum_{Z_{3}\ni
i_{3}}\cdots\sum_{i_{n}\in L\cup Z_{1}\cup\cdots\cup Z_{n-1}}\sum_{Z_{n}\ni
i_{n}}\|h_{Z_{1}}\|\cdots\|h_{Z_{n}}\|$ $\displaystyle\leq{(2gu)^{n}\over
n!}\sum_{i_{1}\in L}\sum_{i_{2}\in L\cup Z_{1}}\sum_{i_{3}\in L\cup Z_{1}\cup
Z_{2}}\cdots\sum_{i_{n}\in L\cup Z_{1}\cup\cdots\cup
Z_{n-1}}\leq{(2gu)^{n}\over n!}|L|(|L|+k)\cdots(|L|+(n-1)k)\,$
$\displaystyle\leq{(2gu)^{n}\over(n/e)^{n}}|L|(|L|+k)\cdots(|L|+(n-1)k)\leq(2guek)^{n}\bigl{(}1+{|L|/k\over
n}\bigr{)}^{n}$
$\displaystyle\leq(2guek)^{n}e^{|L|/k}\leq(2guek)^{n}e^{(|X|+|Y|+s^{\prime}k)/k}\,,$
(S.43)

where we have used $\sum_{Z\ni i}\|h_{Z}\|<gu$, $(n/e)^{n}<n!$,
$(1+x/n)^{n}<e^{x}$ and $|L|<|X|+|Y|+s^{\prime}k$.

Finally we sum over $s^{\prime}$ to get the following expression

$\displaystyle\|\rho_{\rm cl}\|_{1}$ $\displaystyle\leq
u^{-1}|X|\,|Y|{e^{(|X|+|Y|)/k}\over(1+d_{X,Y})^{\alpha}}\sum_{s=0}^{\infty}(2\beta
gkeu\,e^{2\beta gk})^{s}\sum_{s^{\prime}=1}^{s}{(1/k)^{s^{\prime}}\over
s^{\prime}!}$
$\displaystyle\leq(e^{1/k}/u)|X|\,|Y|{e^{(|X|+|Y|)/k}\over(1+d_{X,Y})^{\alpha}}{1\over
1-2\beta gkeu\,e^{2\beta gk}}\,.$ (S.44)

Here, we use
$\sum_{s^{\prime}=1}^{s}{(1/k)^{s^{\prime}}/s^{\prime}!}<e^{1/k}$. This
completes the proof of the Theorem 2.

##  S.IV Exact Diagonalization Results for Long-Range Heisenberg Chains

Let’s examine a more general 1D long-range Heisenberg chain described by the
Hamiltonian

$\displaystyle H=\sum_{1\leq i<j\leq
N}\frac{a_{i,j}}{d_{i,j}^{\alpha}}{\bm{S}}_{i}\cdot{\bm{S}}_{j},$ (S.45)

where $a_{i,j}$ is a parameter of order $\mathcal{O}(1)$, $N$ is the number of
sites, and ${\bm{S}}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$ is a spin-1/2
operator at site $i$. To support our argument that the thermal area holds for
$\alpha>(D+1)/2$ and this criterion is optimal, we compute the mutual
information $\mathcal{I}_{\rho_{\beta}}(A:B)$ between subsystems $A$ and $B$,
with $A$ representing the first $N/2$ sites and $B$ the remaining $N/2$ sites,
by using the exact diagonalization method. We randomly assign $a_{i,j}$ from
the interval $[0,1]$, and then compute the average over 2000 samples. Fig. 6
shows that the mutual information increases with the system size when
$\alpha<1$, while it does not exhibit an increase for $\alpha>1$.

Figure 6: Exact diagonalization results: the averaged mutual information
$\mathcal{I}_{\rho_{\beta}}(A:B)$ over 2000 samples of the long-range
Heisenberg chain (S.45) with random numbers $a_{i,j}$ assigned from the
interval $[0,1]$. Here $N$ is the number of sites, and $A$ and $B$ are the
first $N/2$ sites and the remaining $N/2$ sites, respectively.

##  S.V Numerical Detail of Shapourian-Shiozaki-Ryu Negativity

We consider $n$ spinless fermionic particles with the basis states $|0\rangle$
and $|1\rangle=c^{\dagger}|0\rangle$. The determinant of the $n\times n$
matrix $\mathbf{X}=(x_{i,j})$ is defined as

$\displaystyle\det\mathbf{X}=\sum_{\sigma\in
S_{n}}\mathrm{sgn}(\sigma)x_{1,\sigma(1)}\cdots x_{n,\sigma(n)},$ (S.46)

and the Pfaffian of the $2n\times 2n$ skew-symmetric matrix
$\mathbf{Y}=(y_{i,j})$, i.e. $\mathbf{Y}^{\mathrm{T}}=-\mathbf{Y}$, is defined
as

$\displaystyle\mathrm{Pf}[\mathbf{Y}]=\frac{1}{2^{n}n!}\sum_{\sigma\in
S_{2n}}\mathrm{sgn}(\sigma)y_{\sigma(1)\sigma(2)}\cdots
y_{\sigma(2n-1),\sigma(2n)}.$ (S.47)

Here, $S_{n}$ and $\mathrm{sgn}(\sigma)$ denotes the symmetric group and the
sign of the permutation $\sigma$, respectively. For an arbitrary $2n\times 2n$
matrix $\mathbf{Z}$, we use the following expression

$\displaystyle\mathbf{Z}=\begin{pmatrix}[\mathbf{Z}]^{(1,1)}&[\mathbf{Z}]^{(1,2)}\\\
[\mathbf{Z}]^{(2,1)}&[\mathbf{Z}]^{(2,2)}\end{pmatrix}$ (S.48)

to describe the $n\times n$ block matrix $[\mathbf{Z}]^{(a,b)}$ with $a$-th
row and $b$-th column ($a,b\in\\{1,2\\}$). We denote $\mathds{1}_{n}$ as an
$n\times n$ identity matrix.

We use the Grassmann variables $\\{\xi_{j},\overline{\xi}_{j}\\}$ to describe
the fermionic systems with the following notation:

$\displaystyle\xi$ $\displaystyle\equiv(\xi_{1},\cdots,\xi_{n}),$ (S.49)
$\displaystyle(\xi,\overline{\xi})$
$\displaystyle\equiv(\xi_{1},\cdots,\xi_{n},\overline{\xi}_{1},\cdots,\overline{\xi}_{n}),$
(S.50)

with the differentials

$\displaystyle d\xi$ $\displaystyle\equiv d\xi_{n}\cdots d\xi_{1},$ (S.51)
$\displaystyle d\overline{\xi}d\xi$ $\displaystyle\equiv
d\overline{\xi}_{n}\cdots d\overline{\xi}_{1}d\xi_{n}\cdots d\xi_{1},$ (S.52)
$\displaystyle d(\overline{\xi},\xi)$ $\displaystyle\equiv
d\overline{\xi}_{1}d\xi_{1}\cdots
d\overline{\xi}_{n}d\xi_{n}=(-1)^{n(n-1)/2}d\overline{\xi}d\xi.$ (S.53)

We define the state with the Grassmann variable as

$\displaystyle|\xi_{j}\rangle$
$\displaystyle\equiv|0\rangle_{j}-\xi_{j}|1\rangle_{j},$ (S.54)
$\displaystyle|\overline{\xi}_{j}\rangle$ $\displaystyle\equiv\langle
0|_{j}-\overline{\xi}_{j}\langle 1|_{j},$ (S.55)
$\displaystyle|\xi\rangle\langle\overline{\xi}|$
$\displaystyle=|\xi_{1}\rangle\langle\overline{\xi}_{1}|\otimes\cdots\otimes|\xi_{n}\rangle\langle\overline{\xi}_{n}|.$
(S.56)

The definition naturally leads to
$\langle\overline{\xi}_{j}|\xi_{j}\rangle=e^{\overline{\xi}_{j}\xi_{j}}$. We
also use the notation $(\chi,\overline{\chi})$ to describe another Grassmann
variables.

In terms of the Grassmann variables, the trace of an operator $O$ is given by

$\displaystyle\mathrm{Tr}[O]=\int\langle-\overline{\xi}|O|\xi\rangle
e^{-\sum_{j}\overline{\xi}_{j}\xi_{j}}\,d(\overline{\xi},\xi).$ (S.57)

We consider the following form of the density operator:

$\displaystyle\rho=\frac{1}{Z_{\rho}}\int
d\overline{\xi}d\xi\,|\xi\rangle\langle\overline{\xi}|e^{(\xi,\overline{\xi})\frac{1}{2}\bm{\Gamma}(\xi,\overline{\xi})^{\mathrm{T}}+\sum_{j}\overline{\xi}_{j}\xi_{j}}.$
(S.58)

Imposing $\mathrm{Tr}[\rho]=1$ with Eqs. (S.57) and (S.58), we have
$Z_{\rho}=\mathrm{Pf}[\bm{\Gamma}]$. Therefore, the density matrix is written
as

$\displaystyle\rho=\frac{1}{\mathrm{Pf}[\bm{\Gamma}]}\int
d\overline{\xi}d\xi\,|\xi\rangle\langle\overline{\xi}|e^{(\xi,\overline{\xi})\frac{1}{2}\bm{\Gamma}(\xi,\overline{\xi})^{\mathrm{T}}+\sum_{j}\overline{\xi}_{j}\xi_{j}}.$
(S.59)

Using Eq.(S.57), we can relate the two-point correlations functions of the
fermionic creation and annihilation operators to the matrix $\bm{\Gamma}$ by

$\displaystyle\bm{\Gamma}^{-1}=\begin{pmatrix}[\bm{\Gamma}^{-1}]^{(1,1)}&[\bm{\Gamma}^{-1}]^{(1,2)}\\\
[\bm{\Gamma}^{-1}]^{(2,1)}&[\bm{\Gamma}^{-1}]^{(2,2)}\end{pmatrix}=\begin{pmatrix}\langle
c_{j}c_{i}\rangle_{\rho}&-\langle c_{j}^{\dagger}c_{i}\rangle_{\rho}\\\
\langle c_{i}^{\dagger}c_{j}\rangle_{\rho}&\langle
c_{j}^{\dagger}c_{i}^{\dagger}\rangle_{\rho}\end{pmatrix}.$ (S.60)

In other words, we identify the density matrix with the covaraince matrix
$\bm{\Gamma}$.

Now we move onto the Shapourian-Shiozaki-Ryu (SSR) negativity. We partition
the total system $\Lambda$ into $A$ and $B$ with $|A|=n_{A}$ and $|B|=n_{B}$.
To define the SSR negativity between $A$ and $B$, we need the partial time-
reversal transformation $R_{A}$ on the subsystem $A$. The time-reversal
transformation maps $|\xi\rangle\langle\overline{\xi}|$ to
$|i\overline{\xi}\rangle\langle i\xi|$ in the coherent basis. Therefore, we
have

$\displaystyle(|\\{\xi_{j}\\}_{j\in A},\\{\xi_{j}\\}_{j\in
B}\rangle\langle\\{\overline{\chi}_{j}\\}_{j\in
A},\\{\overline{\chi}_{j}\\}_{j\in
B}|)^{R_{A}}=|\\{i\overline{\chi}_{j}\\}_{j\in A},\\{\xi_{j}\\}_{j\in
B}\rangle\langle\\{i\xi_{j}\\}_{j\in A},\\{\overline{\chi}_{j}\\}_{j\in B}|,$
(S.61)

and therefore the density matrix after the partial time-reversal
transformation becomes

$\displaystyle\rho^{R_{A}}=\frac{1}{\mathrm{Pf}[\bm{\Gamma}]}\int
d\overline{\xi}\,d\xi\,|i\overline{\xi}_{A},\xi_{B}\rangle\langle
i\xi_{A},\overline{\xi}_{B}|e^{(\xi,\overline{\xi})\frac{1}{2}\bm{\Gamma}(\xi,\overline{\xi})^{\mathrm{T}}+\sum_{j}\overline{\xi}_{j}\xi_{j}}.$
(S.62)

Now, we rearrange the Grassmann variables. Let the vector
$(\xi^{\prime},\overline{\xi}^{\prime})\equiv(\xi_{A}^{\prime},\xi_{B}^{\prime},\overline{\xi}_{A}^{\prime},\overline{\xi}_{B}^{\prime})=(-i\overline{\xi}_{A},\xi_{B},-i\xi_{A},\overline{\xi}_{B})$,
then we represent it by introducing the matrix $\mathbf{T}$:

$\displaystyle(\xi^{\prime},\overline{\xi}^{\prime})^{\mathrm{T}}=\mathbf{T}(\xi,\overline{\xi})^{\mathrm{T}},\qquad\mathbf{T}=\begin{pmatrix}0&0&-i\mathds{1}_{n_{A}}&0\\\
0&\mathds{1}_{n_{B}}&0&0\\\ -i\mathds{1}_{n_{A}}&0&0&0\\\
0&0&0&\mathds{1}_{n_{B}}\end{pmatrix}.$ (S.63)

Here, $\mathds{1}_{n_{A}}$ ($\mathds{1}_{n_{B}}$) is an $n_{A}\times n_{A}$
($n_{B}\times n_{B}$) identity matrix. With this matrix, a partial transformed
density matrix is rewritten as follows:

$\displaystyle\rho^{R_{A}}=\frac{1}{\mathrm{Pf}[\bm{\Gamma}]}\int
d\overline{\xi}d\xi\,|\xi\rangle\langle\overline{\xi}|e^{(\xi,\overline{\xi})\frac{1}{2}\mathbf{S}^{\prime}(\xi,\overline{\xi})^{\mathrm{T}}},$
(S.64)

with a newly defined matrix

$\displaystyle\mathbf{S}^{\prime}:=\mathbf{T}\Bigg{[}\bm{\Gamma}+\begin{pmatrix}0&-\mathds{1}_{n}\\\
\mathds{1}_{n}&0\end{pmatrix}\Bigg{]}\mathbf{T}.$ (S.65)

The next step is to consider the product
$(\rho^{R_{A}})^{\dagger}\rho^{R_{A}}$. We introduce the matrix
$\mathbf{S}^{\prime\prime}$ for the operator $(\rho^{R_{A}})^{\dagger}$:

$\displaystyle\mathbf{S}^{\prime\prime}=\begin{pmatrix}([\mathbf{S}^{\prime}]^{(2,2)})^{\dagger}&([\mathbf{S}^{\prime}]^{(1,2)})^{\dagger}\\\
([\mathbf{S}^{\prime}]^{(2,1)})^{\dagger}&([\mathbf{S}^{\prime}]^{(1,1)})^{\dagger}\end{pmatrix}.$
(S.66)

This leads to

$\displaystyle(\rho^{R_{A}})^{\dagger}$
$\displaystyle=\frac{1}{\mathrm{Pf}[\bm{\Gamma}]}\int
d\overline{\xi}d\xi\,|\xi\rangle\langle\overline{\xi}|e^{(\xi,\overline{\xi})\frac{1}{2}(\mathbf{S}^{\prime\prime})(\xi,\overline{\xi})^{\mathrm{T}}}.$
(S.67)

After a straightforward computation, we find the following expression for the
product $(\rho^{R_{A}})^{\dagger}\rho^{R_{A}}$:

$\displaystyle(\rho^{R_{A}})^{\dagger}\rho^{R_{A}}$
$\displaystyle=\frac{(-1)^{n^{2}}\mathrm{Pf}[\mathbf{B}]}{\mathrm{Pf}[\bm{\Gamma}]^{2}}\int
d\overline{\xi}d\xi\,|\xi\rangle\langle\overline{\xi}|e^{(\xi,\overline{\xi})\frac{1}{2}\bm{\Gamma}^{\prime}(\xi,\overline{\xi})^{\mathrm{T}}+\sum_{j}\overline{\xi}_{j}\xi_{j}},$
(S.68)

with introducing two matrices

$\displaystyle\mathbf{B}$
$\displaystyle=\begin{pmatrix}[\mathbf{S}^{\prime}]^{(1,1)}&-\mathds{1}_{n}\\\
\mathds{1}_{n}&[\mathbf{S}^{\prime\prime}]^{(2,2)}\end{pmatrix},$ (S.69)
$\displaystyle\bm{\Gamma^{\prime}}$
$\displaystyle=-\begin{pmatrix}0&[\mathbf{S}^{\prime\prime}]^{(1,2)}\\\
[\mathbf{S}^{\prime}]^{(2,1)}&0\end{pmatrix}\mathbf{B}^{-1}\begin{pmatrix}0&[\mathbf{S}^{\prime}]^{(1,2)}\\\
[\mathbf{S}^{\prime\prime}]^{(2,1)}&0\end{pmatrix}+\begin{pmatrix}[\mathbf{S}^{\prime\prime}]^{(1,1)}&0\\\
0&[\mathbf{S}^{\prime}]^{(2,2)}\end{pmatrix}-\begin{pmatrix}0&-\mathds{1}_{n}\\\
\mathds{1}_{n}&0\end{pmatrix}.$ (S.70)

Then the SSR negativity is

$\displaystyle
E_{\mathrm{SSR}}=\log\|\rho^{R_{A}}\|_{1}=\log\mathrm{Tr}\Big{[}\sqrt{(\rho^{R_{A}})^{\dagger}\rho^{R_{A}}}\Big{]}.$
(S.71)

###  S.V.1 Numerical Procedure

Based on the above calculations, we present the following procedure to
numerically calculate the SSR negativity.

  * 1.

We first compute the covariance matrix $\bm{\Gamma}$ from Eq. (S.60). Then we
calculate $\mathbf{S}^{\prime},\mathbf{S}^{\prime\prime},$ and $\mathbf{B}$
from Eqs. (S.65), (S.66), and (S.69) to construct $\bm{\Gamma}^{\prime}$ by
Eq. (S.70).

  * 2.

Note the relation

$\displaystyle(\rho^{R_{A}})^{\dagger}\rho^{R_{A}}$
$\displaystyle=\frac{(-1)^{n^{2}}\mathrm{Pf}[\mathbf{B}]\mathrm{Pf}[\bm{\Gamma}^{\prime}]}{\mathrm{Pf}[\bm{\Gamma}]^{2}}\rho^{\prime\prime},$
(S.72) $\displaystyle\rho^{\prime\prime}$
$\displaystyle=\frac{1}{\mathrm{Pf}[\bm{\Gamma}^{\prime}]}\int
d\overline{\xi}d\xi\,|\xi\rangle\langle\overline{\xi}|e^{(\xi,\overline{\xi})\frac{1}{2}\bm{\Gamma}^{\prime}(\xi,\overline{\xi})^{\mathrm{T}}+\sum_{j}\overline{\xi}_{j}\xi_{j}}.$
(S.73)

From $\bm{\Gamma}^{\prime}$, using Eq. (S.60) again, we define the covariance
matrix $\mathbf{C}$ of $\rho^{\prime\prime}$:

$\displaystyle\mathbf{C}=\begin{pmatrix}\langle
c_{i}^{\dagger}c_{j}\rangle_{\rho^{\prime\prime}}&\langle
c_{i}^{\dagger}c_{j}^{\dagger}\rangle_{\rho^{\prime\prime}}\\\ \langle
c_{i}c_{j}\rangle_{\rho^{\prime\prime}}&\langle
c_{i}c_{j}^{\dagger}\rangle_{\rho^{\prime\prime}}\end{pmatrix}=\begin{pmatrix}[(\bm{\Gamma}^{\prime})^{-1}]^{(2,1)}&\Big{(}[(\bm{\Gamma}^{\prime})^{-1}]^{(2,2)}\Big{)}^{\mathrm{T}}\\\
\Big{(}[(\bm{\Gamma}^{\prime})^{-1}]^{(1,1)}\Big{)}^{\mathrm{T}}&\mathds{1}_{n}+[(\bm{\Gamma}^{\prime})^{-1}]^{(1,2)}\end{pmatrix}.$
(S.74)

  * 3.

We choose the unitary transformation $U$ to define another fermionic
annihilation operators $\mathbf{d}=(d_{1},\cdots,d_{n})$:

$\displaystyle\begin{pmatrix}\mathbf{d}^{\dagger}\\\
\mathbf{d}\end{pmatrix}=U\begin{pmatrix}\mathbf{c}^{\dagger}\\\
\mathbf{c}\end{pmatrix},$ (S.75)

so that it diagonalizes $\rho^{\prime\prime}$ as

$\displaystyle\rho^{\prime\prime}=\frac{e^{-\sum_{i=1}^{n}\epsilon_{i}^{\prime\prime}d_{i}^{\dagger}d_{i}}}{Z^{\prime\prime}},$
(S.76)

where
$Z^{\prime\prime}=\mathrm{Tr}[\rho^{\prime\prime}]=\prod_{i=1}^{n}(1+e^{-\epsilon_{i}^{\prime\prime}})$.
Equation. (S.76) yields

$\displaystyle\begin{pmatrix}\langle
d_{i}^{\dagger}d_{j}\rangle_{\rho^{\prime\prime}}&\langle
d_{i}^{\dagger}d_{j}^{\dagger}\rangle_{\rho^{\prime\prime}}\\\ \langle
d_{i}d_{j}\rangle_{\rho^{\prime\prime}}&\langle
d_{i}d_{j}^{\dagger}\rangle_{\rho^{\prime\prime}}\end{pmatrix}=\begin{pmatrix}\frac{1}{e^{\epsilon_{i}^{\prime\prime}}+1}\delta_{i,j}&0\\\
0&\frac{e^{\epsilon_{i}^{\prime\prime}}}{e^{\epsilon_{i}^{\prime\prime}}+1}\delta_{i,j}\end{pmatrix}=U\begin{pmatrix}\langle
c_{i}^{\dagger}c_{j}\rangle_{\rho^{\prime\prime}}&\langle
c_{i}^{\dagger}c_{j}^{\dagger}\rangle_{\rho^{\prime\prime}}\\\ \langle
c_{i}c_{j}\rangle_{\rho^{\prime\prime}}&\langle
c_{i}c_{j}^{\dagger}\rangle_{\rho^{\prime\prime}}\end{pmatrix}U^{\dagger}.$
(S.77)

It means that we can obtain $\epsilon_{i}^{\prime\prime}$ by diagonalizing the
covariance matrix $\mathbf{C}$:

$\displaystyle\mathbf{C}$ $\displaystyle=U^{\dagger}\mathrm{diag}(\langle
d_{i}^{\dagger}d_{i}\rangle_{\rho^{\prime\prime}},\langle
d_{i}d_{i}^{\dagger}\rangle_{\rho^{\prime\prime}})U,$ (S.78)
$\displaystyle\langle d_{i}^{\dagger}d_{i}\rangle_{\rho^{\prime\prime}}$
$\displaystyle=\frac{1}{e^{\epsilon_{i}^{\prime\prime}}+1}.$ (S.79)

  * 4.

From the following identity

$\displaystyle\mathrm{Tr}\Big{[}\sqrt{(\rho^{R_{A}})^{\dagger}\rho^{R_{A}}}\Big{]}$
$\displaystyle=\sqrt{\frac{(-1)^{n^{2}}\mathrm{Pf}[\mathbf{B}]\mathrm{Pf}[\bm{\Gamma}^{\prime}]}{\mathrm{Pf}[\bm{\Gamma}]^{2}}}\mathrm{Tr}\sqrt{\rho^{\prime\prime}}=\sqrt{\frac{\mathrm{Pf}[\mathbf{B}]\mathrm{Pf}[\bm{\Gamma}^{\prime}]}{\mathrm{Pf}[\bm{\Gamma}]^{2}}}\prod_{i=1}^{n}\frac{1+e^{-\frac{\epsilon_{i}^{\prime\prime}}{2}}}{\sqrt{1+e^{-\epsilon_{i}^{\prime\prime}}}},$
(S.80)

we get the following expression of the SSR negativity with computable
quantities:

$\displaystyle E_{\mathrm{SSR}}$
$\displaystyle=\log\mathrm{Tr}\Big{[}\sqrt{(\rho^{R_{A}})^{\dagger}\rho^{R_{A}}}\Big{]}$
(S.81)
$\displaystyle=\frac{1}{2}\log\Bigg{[}\frac{(-1)^{n^{2}}\mathrm{Pf}[\mathbf{B}]\mathrm{Pf}[\bm{\Gamma}^{\prime}]}{\mathrm{Pf}[\bm{\Gamma}]^{2}}\Bigg{]}+\sum_{i=1}^{n}\Bigg{[}\log\Big{(}1+e^{-\frac{\epsilon_{i}^{\prime\prime}}{2}}\Big{)}-\frac{1}{2}\log(1+e^{-\epsilon_{i}^{\prime\prime}})\Bigg{]}$
(S.82)
$\displaystyle=\frac{1}{2}\log\Bigg{[}\frac{(-1)^{n^{2}}\mathrm{Pf}[\mathbf{B}]\mathrm{Pf}[\bm{\Gamma}^{\prime}]}{\mathrm{Pf}[\bm{\Gamma}]^{2}}\Bigg{]}+\sum_{i=1}^{n}\log\Big{(}\sqrt{\langle
d_{i}^{\dagger}d_{i}\rangle_{\rho^{\prime\prime}}}+\sqrt{\langle
d_{i}d_{i}^{\dagger}\rangle_{\rho^{\prime\prime}}}\Big{)}.$ (S.83)