# Minimal energy cost of entanglement extraction

Lucas Hackl  Max Planck Institute of Quantum Optics, Hans-Kopfermann-Str. 1,
85748 Garching, Germany  Munich Center for Quantum Science and Technology,
Schellingstraße 4, D-80799 München, Germany<EMAIL_ADDRESS>Robert H.
Jonsson  QMATH, Department of Mathematical Sciences, University of Copenhagen,
Universitetsparken 5, 2100 Copenhagen, Denmark<EMAIL_ADDRESS>
(August 27, 2024)

###### Abstract

We compute the minimal energy cost for extracting entanglement from the ground
state of a bosonic or fermionic quadratic system. Specifically, we find the
minimal energy increase in the system resulting from replacing an entangled
pair of modes, sharing entanglement entropy $\Delta S$, by a product state,
and we show how to construct modes achieving this minimal energy cost. Thus,
we obtain a protocol independent lower bound on the extraction of pure state
entanglement from quadratic systems. Due to their generality, our results
apply to a large range of physical systems, as we discuss with examples.

##  1 Introduction

Entanglement is a hallmark of quantum theory: On a fundamental level, its very
existence has deep implications for our understanding of nature, _e.g._ , as
pertaining to its realism and localism [1, 2, 3, 4, 5, 6]. From an applied
point of view, entanglement is a key resource in quantum technologies and
quantum information processing [7, 8, 9, 10]. And moreover, studying
entanglement structures provided deep insights in various fields ranging from
quantum many body systems, over quantum field theory to quantum gravity [11,
12, 13, 14, 15, 16, 17].

The fact that the ground states of local Hamiltonians are generally highly
entangled gives rise to the idea that such systems could serve as a source of
entanglement, _e.g._ , for applications in quantum information. In particular,
the idea of extracting entanglement from the ground state of a system has been
explored in the context of the vacuum state of a quantum field. Starting with
the observation that the entanglement present in the vacuum fluctuations can
violate Bell’s inequalities [18, 19, 20, 21, 22], various works studied the
ability of the quantum field vacuum to entangle systems which couple locally
to the quantum field. In particular, the Unruh-DeWitt particle detector model
has been employed to investigate aspects ranging from the impact of
relativistic motion [23], spacetime curvature and topology [24, 25, 26], the
nature of field states and interactions [27, 28, 29, 30, 31, 32], and
potential implementations [33, 34, 35, 36].

Previous works on the relation between entanglement and energy content of a
system focused mainly on the energy cost of creating entanglement and
correlations, or looked at the conversion of one into the other, _e.g._ , in
the context of quantum thermodynamics [37, 38, 39, 40, 41, 42, 43]. Also it
was shown that it is possible to use the entanglement present in the vacuum of
a quantum field to teleport energy from one spacetime region to another. In
quantum energy teleportation one observer performs a local measurement of the
field, classically communicates the result to another observer who then
through local unitary interactions is able to extract energy from the field
[44, 45, 46, 47].

However, the systematic study of the energy cost of entanglement extraction
from complex quantum systems was only initiated recently by Beny et al. [48].
When entanglement is extracted from a system, its state is changed and thus
its energy content may change as well. In particular, the extraction of
entanglement from a system’s _unique_ ground state is inevitably associated
with an energy increase injected into the system when changing its state.

This interplay of entanglement extraction and energy cost is not only relevant
with respect to potential implementations and quantum thermodynamics but, in
fact, also connects to quantum field theory on curved spacetimes. For example,
Jacobson [49] recently showed that the semi-classical Einstein equations are
equivalent to the property of the vacuum of a field to locally maximize the
entanglement entropy. Attempting to increase a sphere’s entanglement by
varying the field state induces an increase in energy density which curves
spacetime and shrinks the sphere’s surface, thus reducing the entanglement
contributions proportional to the surface area, and exactly cancelling out the
attempted increase of entanglement.

In this article, we present a first result in the direction initiated by [48]
which applies to large classes of systems of physical interest: We study the
extraction of entanglement from the ground state of bosonic and fermionic
modes governed by a quadratic Hamiltonian. This also includes spin systems
that can be mapped to quadratic fermionic Hamiltonians via the Jordan-Wigner
transform. In particular, we find the minimal energy cost for the extraction
of mode pairs in a pure entangled state.

An important tool to our analysis is the partner mode formula, which was
developed a in series of recent papers [50, 51, 52], mostly in the context of
bosonic quantum field theory. We generalize this construction, such that it
applies to any bosonic or fermionic pure Gaussian state. This is achieved by
phrasing the construction in terms of the complex structure associated with a
pure Gaussian state [53, 54], thus illuminating the geometric origin of the
partner mode construction.

The partner mode construction is based on the well-known property of pure
Gaussian states to factorize into pairs of entangled mode pairs [55] (see also
[56]). If a system of many modes is in an overall pure Gaussian state, then
for every mode in that system, which is not in a pure state on its own, there
exists a unique partner mode, with which it shares all its entanglement. Thus,
the mode and its partner are in a pure entangled state, that is in a product
state with the rest of the system.

This property is the motivation for extracting entanglement by unitarily
swapping the entangled state of such partner modes with an unentangled state
of target modes in some laboratory system. Thus, the system acts as a source
of an entangled state that is traded with an unentangled state from the target
system. This framework was introduced in [51] where the relationship between
energy cost and entanglement extraction was discussed for a harmonic chain.

In this article, we generalize this framework to arbitrary systems of bosonic
or fermionic modes with quadratic Hamiltonians. We find the minimal energy
cost for the extraction of partner mode pairs from their ground state as
analytical functions of the Hamiltonian’s spectrum. Furthermore, we show how
to explicitly construct partner modes that achieve this minimal energy cost.

The structure of this article is as follows: In Section 2, we review bosonic
Gaussian states and quadratic Hamiltonians. In Section 3, we present our
approach of extracting entanglement from partner modes and derive an
analytical formula for the minimal energy cost. The following two Sections 4
and 5 mirror the previous two for fermionic systems, _i.e._ , we first review
fermionic Gaussian states and quadratic Hamiltonians in Section 4 and then
derive the respective minimal energy cost for fermions in Section 5. Finally,
Section 6 discusses applications of our results for both bosonic and fermionic
models, before we conclude with a discussion in Section 7. Appendix A
rederives the pairwise correlation structure of bosonic and fermionic Gaussian
states in the newly developed language of linear complex structures. Appendix
B reviews relevant results in linear algebra on how the spectrum of a linear
map changes under restrictions.

The article is structured to be largely self-contained. To the experienced
reader Sections 2 and 4 will mostly serve to establish notation for bosons and
fermions respectively, and to present a covariant way of computing partner
modes.

##  2 Quadratic bosonic systems and their partner mode construction

The partner mode construction is a central tool in our derivation of the
minimal energy of entanglement extraction. In this section, we introduce the
partner mode formula for arbitrary pure Gaussian states of bosonic modes
(specifically in Section 2.5). To be largely self-contained and to establish
notation, we begin with a review of quadratic bosonic systems. We follow the
conventions of [57, 58, 59], while other comprehensive reviews include [60,
14, 61, 62, 63].

###  2.1 Phase space and observables

For $N$ bosonic degrees of freedom, we consider a classical phase space
$V\simeq\mathbb{R}^{2N}$ and its dual $V^{*}\simeq\mathbb{R}^{2N}$. We denote
vectors $v^{a}\in V$ with an upper index and dual vectors $w_{a}\in V^{*}$
with a lower index. As phase space, both $V$ and $V^{*}$ are equipped with
symplectic forms, namely $\Omega^{ab}:V^{*}\times V^{*}\to\mathbb{R}$ on
$V^{*}$ and its inverse $\Omega_{ab}^{-1}:V\times V\to\mathbb{R}$ satisfying
$\Omega^{ac}\Omega_{cb}^{-1}=\delta^{a}{}_{b}$.

We can characterize a vector $v^{a}\in V$ by its coordinate values with
respect to a basis as

$\displaystyle
v^{a}\equiv\left(q_{1}(v),p_{1}(v),\cdots,q_{N}(v),p_{N}(v)\right)^{\intercal}\in\mathbb{R}^{2N}\,,$
(1)

where “$\equiv$” indicates that $v^{a}$ is represented by the column vector
with respect to the specified basis.

Typically, we choose coordinate functions $q_{i},p_{i}\in C^{\infty}(V)$ which
are linear, _i.e._ , $q_{i},p_{i}\in V^{*}$, and which consist of conjugate
variables, _i.e._ , satisfying $\Omega(q_{i},p_{j})=\delta_{ij}$.

With respect to such a basis, the symplectic form $\Omega^{ab}$ takes the
standard Darboux form, _i.e._ , is represented by the block diagonal matrix

$\displaystyle\Omega^{ab}$
$\displaystyle\equiv\mathbf{\Omega}=\left(\begin{array}[]{ccc}\mathbf{\Omega}_{2}&&0\\\
&\ddots&\\\ 0&&\mathbf{\Omega}_{2}\end{array}\right)\,,$ (5)

where we use boldface symbols $\mathbf{\Omega}$ for the matrix representation
with respect to a specific basis, and define the $2\times 2$-matrix

$\displaystyle\mathbf{\Omega}_{2}=\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}\,.$
(6)

In the classical theory, observables correspond to smooth functions on the
phase space $V$, _i.e._ , functions in $C^{\infty}(V)$. In particular, a
covector $w\in V^{*}$ defines a linear observable through

$\displaystyle\mathcal{O}_{w}:V\to\mathbb{R}:v^{a}\mapsto w_{a}v^{a}\,.$ (7)

Similarly, a bilinear form $h_{ab}$ on $V$ describes uniquely a quadratic
observable

$\displaystyle\mathcal{O}_{h}:V\to\mathbb{R}:v^{a}\mapsto\frac{1}{2}h_{ab}v^{a}v^{b}\,.$
(8)

In quantum theory, observables correspond to linear operators on a Hilbert
space $\mathcal{H}$. To define linear and multi-linear observables it is
convenient to introduce the operator-valued vector $\hat{\xi}^{a}$ referred to
as quantization map. It is subject to the constraint

$\displaystyle[\hat{\xi}^{a},\hat{\xi}^{b}]=\hat{\xi}^{a}\hat{\xi}^{b}-\hat{\xi}^{b}\hat{\xi}^{a}=\\{\xi^{a},\xi^{b}\\}=\mathrm{i}\Omega^{ab}\mathbb{1}$
(9)

where $\mathbb{1}$ represents the identity operator on $\mathcal{H}$. With
respect to our initial basis where $\Omega$ takes the form of (5),
$\hat{\xi}^{a}$ is represented as

$\displaystyle\hat{\xi}^{a}\equiv(\hat{q}_{1},\cdots,\hat{q}_{N},\hat{p}_{1},\cdots,\hat{p}_{N})^{\intercal}$
(10)

with the familiar position and momentum operators. From here, we can construct
the well-known creation and annihilation operators
$\hat{a}_{i}=\frac{1}{\sqrt{2}}(\hat{q}_{i}-\mathrm{i}\hat{p}_{i})$ and
$\hat{a}_{i}^{\dagger}=\frac{1}{\sqrt{2}}(\hat{q}_{i}+\mathrm{i}\hat{p}_{i})$.

The quantum observables associated to the linear form $w_{a}$ and quadratic
form $h_{ab}$ are then

$\displaystyle\hat{\mathcal{O}}_{w}=w_{a}\hat{\xi}^{a}\quad\text{and}\quad\hat{\mathcal{O}}_{h}=\frac{1}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}\,.$
(11)

In contrast to classical observables, a quantum observable constructed from a
general tensor $t_{a_{1}\cdots a_{n}}$ depends on the tensor’s non-symmetric
part since the operators $\hat{q}_{i}$ and $\hat{p}_{i}$ do not commute.

Note that in the classical theory, the quantization map reduces to the vector-
valued function

$\displaystyle\xi^{a}:V\to V:v^{a}\mapsto\xi^{a}(v)=v^{a}\,,$ (12)

which is, in fact, just the identity map on $V$. It can be used to represent
linear and multi-linear observables, _e.g._ , through

$\displaystyle\mathcal{O}_{t}:v\mapsto\frac{1}{n!}t_{a_{1}\cdots
a_{n}}\xi^{a_{1}}(v)\cdots\xi^{a_{n}}(v)$ (13)

for a general tensor $t_{a_{1}\cdots a_{n}}$. In terms of a linear basis of
$V$ it is represented as the tupel of functions

$\displaystyle\xi^{a}\equiv(q_{1},p_{1},\cdots,q_{N},p_{N})^{\intercal}.$ (14)

Because each component of $\xi^{a}$, _i.e._ , $q_{i},p_{i}:V\to\mathbb{R}$, is
a function on phase space, the vector valued function encodes the Poisson
brackets of the classical theory through

$\displaystyle\\{\xi^{a},\xi^{b}\\}=\Omega^{ab}\,.$ (15)

This is particularly convenient and allows for a clean notation when
expressing Poisson brackets of observables abstractly. For instance, for two
linear observables $\mathcal{O}_{w}=w_{a}\xi^{a}$ and
$\mathcal{O}_{u}=u_{a}\xi^{a}$, we find immediately

$\displaystyle\\{\mathcal{O}_{w},\mathcal{O}_{u}\\}=w_{a}u_{b}\\{\xi^{a},\xi^{b}\\}=w_{a}u_{b}\Omega^{ab}\,.$
(16)

###  2.2 Bosonic Gaussian states

Bosonic Gaussian states are an important class of quantum states because they
can be described through powerful analytic methods [60]. They appear as ground
states of quadratic Hamiltonians and their wave function representation is
given by multi-dimensional complex Gaussian distributions. When representing
them as quasi-probability distributions on classical phase space, they are the
only states with positive Wigner functions [64, 65]. In particular, their
Wigner functions are multi-dimensional real Gaussians again. In the context of
Bogoliubov transformations, Gaussian states are often identified as the vacuum
associated to a set of annihilation operators satisfying bosonic commutation
relations. Finally, their $n$-point correlation functions are completely
determined from their $1$-point and 2-point correlation function via Wick’s
theorem [66].

We label a Gaussian state $\ket{G,z}$ by a displacement vector $z^{a}$ and a
covariance matrix $G^{ab}$

$\displaystyle z^{a}$ $\displaystyle=\bra{G,z}\hat{\xi}^{a}\ket{G,z}\,,$ (17)
$\displaystyle G^{ab}$
$\displaystyle=\bra{G,z}\hat{\xi}^{a}\hat{\xi}^{b}+\hat{\xi}^{b}\hat{\xi}^{a}\ket{G,z}-2z^{a}z^{b}\,.$
(18)

The covariance matrix $G^{ab}$ is a positive definite symmetric bilinear form
on the dual phase space $V^{*}$, _i.e._ , $G$ equips $V^{*}$ with an inner
product. In particular, for each covariance matrix there exists a basis of the
phase space $V$ with respect to which the covariance matrix is represented by
the identity matrix $G^{ab}\equiv\mathbb{1}$, while $\Omega$ takes the
standard form (5).

Furthermore, the covariance matrix of a pure Gaussian state and the inverse
symplectic form satisfy the condition

$\displaystyle
G^{ab}\Omega^{-1}_{bc}G^{cd}\Omega^{-1}_{de}=-\delta^{a}{}_{e}\,.$ (19)

In fact, we can use this property to define a linear map on $V$

$\displaystyle J^{a}{}_{b}=-G^{ac}\Omega^{-1}_{cb}\,.$ (20)

This map is referred to as linear complex structure on $V$ because it
satisfies

$\displaystyle J^{2}=(G\Omega^{-1})^{2}=-\mathbb{1}\,,$ (21)

_i.e._ , it acts like the multiplication with the complex number $\mathrm{i}$
on the real vector space. For a given covariance matrix and displacement
vector, the Gaussian state $\ket{G,z}\in\mathcal{H}$ is then fully specified
by the requirement that

$\displaystyle\frac{1}{2}(\delta^{a}{}_{b}+\mathrm{i}J^{a}{}_{b})(\hat{\xi}^{b}-z^{b})\ket{G,z}=0\,.$
(22)

Here, the first term corresponds to a projector onto the space spanned by
those annihilation operators that annihilate $\ket{G,z}$. Describing Gaussian
states in terms of linear complex structures was first introduced to describe
unitarily inequivalent Fock space representations [53, 54, 63], but since then
has been used to study entanglement production [57, 58], typical entanglement
of energy eigenstates [67, 68, 69] and to explore circuit complexity in free
field theories [68, 70].

Wick’s theorem provides an efficient way to compute arbitrary $n$-point
functions

$\displaystyle C^{a_{1}\cdots
a_{n}}_{n}=\bra{G,z}(\hat{\xi}^{a_{1}}-z^{a_{1}})\cdots(\hat{\xi}^{a_{n}}-z^{a_{n}})\ket{G,z}.$
(23)

where the subtraction of $z^{a}$ simplifies later expressions. Clearly,
knowing $C^{a_{1}\cdots a_{n}}_{n}$ enables us to compute arbitrary
correlations even without $z^{a}$. With the definitions of $G^{ab}$ and
$z^{a}$, and the commutation relations
$[\hat{\xi}^{a},\hat{\xi}^{b}]=\mathrm{i}\Omega^{ab}$, we find the 2-point
function to be

$\displaystyle C^{ab}_{2}=\frac{1}{2}(G^{ab}+\mathrm{i}\Omega^{ab})\,.$ (24)

With this expression, we can state Wick’s theorem compactly as

$\displaystyle C^{a_{1}\cdots a_{n}}_{n}$
$\displaystyle=\sum\text{(contractions of $C_{2}^{ab}$)}$ (25)
$\displaystyle=C^{a_{1}a_{2}}_{2}\cdots C^{a_{n-1}a_{n}}_{2}+\cdots\,,$ (26)

where the sum over all contraction refers to writing $C^{a_{1}\cdots
a_{n}}_{n}$ as sum of products of 2-point functions $C_{2}^{ab}$ with all
possible inequivalent assignments of ordered indices $a_{i}$. This implies
that all odd $n$-point functions vanish, because contractions require an even
number of indices. Therefore, after the 2-point function, the next non-trivial
correlation function is the four-point function given by

$\displaystyle
C_{4}^{abcd}=C_{2}^{ab}C_{2}^{cd}+C_{2}^{ac}C_{2}^{bd}+C_{2}^{ad}C_{2}^{bc}\,.$
(27)

Note the importance of ordering the indices on each 2-point function
$C^{a_{i}a_{j}}_{2}$ in the same way as on $C_{n}^{a_{1}\cdots a_{n}}$, _i.e._
, we require $i<j$.

##### Example.

In order to illustrate our methods and to fix conventions, let us look at the
ground state $\ket{G,z}$ of the harmonic oscillator with Hamiltonian
$\hat{H}=\frac{1}{2}(\hat{p}^{2}+\omega^{2}\hat{q}^{2})$. Clearly, we have
$z^{a}=0$, because the ground state has zero expectation value of position
$\hat{q}$ and momentum $\hat{p}.$ We represent everything in the basis
$\hat{\xi}^{a}\equiv(\hat{q},\hat{p})^{\intercal}$. Here, we have

$\displaystyle\Omega^{ab}\equiv\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}\,\,\,\text{and}\,\,\,G^{ab}\equiv\begin{pmatrix}1/\omega&0\\\
0&\omega\end{pmatrix}\,.$ (28)

We can use $\Omega^{ab}$ and $G^{ab}$ to compute the matrix representation of
the linear complex structure

$\displaystyle
J^{a}{}_{b}=-G^{ac}\Omega^{-1}_{cb}\equiv\begin{pmatrix}0&1/\omega\\\
-\omega&0\end{pmatrix}\,.$ (29)

With this, the projection of equation (22) reads

$\displaystyle\frac{1}{2}(\delta^{a}{}_{b}+\mathrm{i}J^{a}{}_{b})\hat{\xi}^{b}\equiv\frac{1}{2}\begin{pmatrix}\hat{q}+\frac{\mathrm{i}}{\omega}\,\hat{p}\\\
\hat{p}-\mathrm{i}\omega\hat{q}\end{pmatrix}=\begin{pmatrix}\sqrt{\frac{1}{2\omega}}\,\hat{a}\\\
-\mathrm{i}\sqrt{\frac{\omega}{2}}\hat{a}\end{pmatrix}\,,$ (30)

_i.e._ , we project onto the subspace spanned by the annihilation operator
$\hat{a}=\sqrt{\frac{\omega}{2}}(\hat{q}+\frac{\mathrm{i}}{\omega}\,\hat{p})$.

###  2.3 Quadratic bosonic Hamiltonians

The most general quadratic Hamiltonian

$\displaystyle\hat{H}=\frac{1}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}+f_{a}\hat{\xi}^{a}\,,$
(31)

contains a linear and a quadratic term, where we require $h_{ab}$ to be
positive definite and symmetric. This ensures that $\hat{H}$ is bounded from
below and its ground state is a Gaussian state $\ket{G,z}$. The covariance
matrix $G^{ab}$ and the linear displacement $z^{a}$ can be computed from
$h_{ab}$ and $f_{a}$:

  * •

Computation of $z^{a}$
The ground state minimizes the energy expectation value
$E=\bra{G,z}\hat{H}{\ket{G,z}}$ given by

$\displaystyle
E=\frac{1}{4}h_{ab}G^{ab}+\frac{1}{2}h_{ab}z^{a}z^{b}+f_{a}z^{a}\,.$ (32)

This implies $z^{a}h_{ab}+f_{b}=0$ giving

$\displaystyle z^{a}=-(h^{-1})^{ab}f_{b}\,,$ (33)

with the inverse form $(h^{-1})^{ac}h_{cb}=\delta^{a}{}_{b}$, which exists
since $h_{ab}$ is positive definite, thus non-degenerate.

  * •

Computation of $G^{ab}$
We can compute $G^{ab}$ directly from $h_{ab}$ via

$\displaystyle G^{ab}=-K^{a}{}_{c}|K^{-1}|^{c}{}_{d}\Omega^{db}\,,$ (34)

with $K^{a}{}_{b}=\Omega^{ac}h_{cb}$. $|K^{-1}|$ refers to the linear map
constructed by replacing the eigenvalues of $K^{-1}$ with their absolute
values.

We can think of these formulas as a projection from the space of quadratic
Hamiltonians onto the space of Gaussian states. Every Hamiltonian gives rise
to a unique ground state, but for each state there are plenty of Hamiltonians
that share the same ground state. This redundancy is visible in (34) where
rescaling $h_{ab}$ with an overall factor (or individually in each frequency
sector) leaves $G^{ab}$ unchanged.

The energy expectation value $E$ and its variance $\Sigma_{E}$ for Gaussian
states of the form $\ket{G^{\prime},z}$, _i.e._ , with the above $z^{a}$ but
with arbitrary $G^{\prime}$, are readily calculated as

$\displaystyle E$
$\displaystyle=\bra{G^{\prime},z}\hat{H}\ket{G^{\prime},z}=\operatorname{Tr}(hG^{\prime})/4+E_{z}\,,$
(35) $\displaystyle\Sigma_{E}$
$\displaystyle=\sqrt{\bra{G^{\prime},z}\hat{H}^{2}\ket{G^{\prime},z}-E^{2}}$
(36) $\displaystyle=\sqrt{\operatorname{Tr}(hG^{\prime}hG^{\prime}+h\Omega
h\Omega)}/4\,,$ (37)

where
$E_{z}=-\frac{1}{2}h_{ab}z^{a}z^{b}=-\frac{1}{2}(h^{-1})^{ab}f_{a}f_{b}$. This
follows from (33) to rewrite

$\displaystyle\hat{H}=\frac{1}{2}h_{ab}(\hat{\xi}-z)^{a}(\hat{\xi}-z)^{b}+E_{z}\,.$
(38)

##### Example.

Let us consider the simple harmonic oscillator
$\hat{H}=\frac{1}{2}(\hat{p}^{2}+\omega^{2}\hat{q}^{2})$ from before. We find
$f_{a}=0$ and with respect to $\hat{\xi}^{a}\equiv(\hat{q},\hat{p})$

$\displaystyle h_{ab}\equiv\begin{pmatrix}\omega^{2}&0\\\ 0&1\end{pmatrix}\,.$
(39)

Applying formula (34) yields

$\displaystyle K^{a}{}_{b}\equiv\begin{pmatrix}0&1\\\
-\omega^{2}&0\end{pmatrix}\,\,\Rightarrow\,\,G^{ab}\equiv\begin{pmatrix}1/\omega&0\\\
0&\omega\end{pmatrix}\,,$ (40)

which agrees with the result from (28). Here

$\displaystyle|K^{-1}|\equiv\begin{pmatrix}\frac{1}{\omega}&0\\\
0&\frac{1}{\omega}\end{pmatrix}\,,$ (41)

because $K^{-1}$ has eigenvalues $\pm\frac{\mathrm{i}}{\omega}$.

###  2.4 Entanglement of bosonic Gaussian states

In terms of a quantum system’s Hilbert space $\mathcal{H}$, a bi-partition of
the system corresponds to a tensor product
$\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. A bi-partition of a
system of $N$ bosonic modes, into subsystems of $N_{A}$ and $N_{B}$ modes
respectively, induces a decomposition of the phase space $V=A\oplus B$ into
two symplectic complements $A$ and $B$.

Generally a Gaussian state on the full system $\ket{G,z}$ contains
correlations, _i.e._ , entanglement between $A$ and $B$. However, there always
exist bases $(q_{1}^{A},p_{1}^{A},\cdots,q^{A}_{N_{A}},p^{A}_{N_{A}})$ of $A$
and $(q_{1}^{B},p_{1}^{B},\cdots,q_{N_{B}}^{B},p^{B}_{N_{B}})$ of $B$, such
that the covariance matrix $G^{ab}$ takes the standard form [55] re-derived in
appendix A.1

$\displaystyle
G\equiv\left(\begin{array}[]{ccc|cccccc}\mathbf{ch}_{1}&\cdots&0&\mathbf{sh}_{1}&\cdots&0&0&\cdots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\cdots&\mathbf{ch}_{N_{A}}&0&\cdots&\mathbf{sh}_{N_{A}}&0&\cdots&0\\\
\hline\cr\mathbf{sh}_{1}&\cdots&0&\mathbf{ch}_{1}&\cdots&0&0&\cdots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\cdots&\mathbf{sh}_{N_{A}}&0&\cdots&\mathbf{ch}_{N_{A}}&0&\cdots&0\\\
0&\cdots&0&0&\cdots&0&\mathbb{1}_{2}&\cdots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\cdots&0&0&\cdots&0&0&\cdots&\mathbb{1}_{2}\\\ \end{array}\right)\,,$ (51)

built from the $2\times 2$-matrices, with $r_{i}\geq 0$,

$\displaystyle\mathbf{ch}_{i}$
$\displaystyle=\left(\begin{array}[]{cc}\operatorname{ch}{2r_{i}}&0\\\
0&\operatorname{ch}{2r_{i}}\end{array}\right)\,,$ (54)
$\displaystyle\mathbf{sh}_{i}$
$\displaystyle=\left(\begin{array}[]{cc}\operatorname{sh}{2r_{i}}&0\\\
0&-\operatorname{sh}{2r_{i}}\end{array}\right)\,.$ (57)

The standard form highlights two important features of the state’s
entanglement structure. Firstly, the basis modes bring the diagonal blocks,
corresponding to $A$ and $B$ into diagonal form, _i.e._ , they provide a
normal mode decomposition of the subsystems’ partial states: The partial
states of $A$, and $B$ respectively, factorize into a product state over the
modes in each subsystem basis.

Secondly, all entanglement between $A$ and $B$ is contained in entangled mode
pairs: If $r_{i}>0$, for a mode $(q_{i}^{A},p_{i}^{A})$ in $A$, then it is
entangled with the mode $(q_{i}^{B},p_{i}^{B})$ with which it is in a two-mode
squeezed state.

The total von Neumann entropy $S_{A}(\ket{G,z})$ between $A$ and $B$ is given
by the sum of the mode pairs’ individual entropies [71, 72]

$\displaystyle
S_{A}(\ket{G,z})=\sum_{i=1}^{N_{A}}s_{b}(\operatorname{ch}{2r_{i}})$ (58)

with
$s_{b}(x)=\left(\frac{x+1}{2}\right)\log\left(\frac{x+1}{2}\right)-\left(\frac{x-1}{2}\right)\log\left(\frac{x-1}{2}\right)$.
Here, we chose the convention to compute all logarithms with respect to base
$2$.

Note, that the total entropy can be expressed in a basis invariant way as [57,
59]

$\displaystyle
S_{A}(\ket{G,z})=\mathrm{Tr}\left[\left(\frac{\mathbb{1}_{A}+\mathrm{i}[J]_{A}}{2}\right)\log\left|\frac{\mathbb{1}_{A}+\mathrm{i}[J]_{A}}{2}\right|\,\\!\right]\,,$
(59)

where $[J]_{A}$ is the restriction of the complex structure
$J^{a}{}_{b}=-G^{ac}\Omega^{-1}_{cb}$ to the symplectic subspace $A$, _i.e._ ,
we restrict it to a $2N_{A}$-by-$2N_{A}$ matrix. The absolute value should be
understood in terms of eigenvalues of $[J]_{A}$ which come in conjugate pairs
$\pm\mathrm{i}\operatorname{ch}{2r_{i}}$.

###  2.5 Bosonic partner mode construction

If we pick a single mode from a larger system in a Gaussian state $\ket{G,z}$,
the single mode is in general entangled with the rest of the modes. However,
as implied by the standard form (51), there exists a second mode, the so
called partner mode, which shares all of the first mode’s entanglement, _i.e._
, these two modes are in a product state with the rest of the system.

Recent works have developed formulae for the partner mode [50, 51, 52]. We
here recast them in terms of the complex structure $J$ of the state, thus
illuminating their structure, and generalizing them to arbitrary Gaussian
states. In Appendix A.1 we derive the partner mode construction in terms of
the complex structure.

Let us assume that the two observables $\hat{Q}_{A}=x_{a}\hat{\xi}^{a}$ and
$\hat{P}_{A}=k_{a}\hat{\xi}^{a}$ define a mode, _i.e._ ,

$\displaystyle\Omega^{ab}x_{a}k_{b}=1\,,$ (60)

and that they yield the standard form of the mode according to (51), _i.e._ ,
for some $r>0$

$\displaystyle G^{ab}x_{a}x_{b}$
$\displaystyle=G^{ab}k_{a}k_{b}=\operatorname{ch}2r\,,$ (61) $\displaystyle
G^{ab}x_{a}k_{b}$ $\displaystyle=0\,.$ (62)

Then the partner mode of this mode is given by
$\hat{Q}_{\bar{A}}=\bar{x}_{a}\hat{\xi}^{a}$ and
$\hat{P}_{\bar{A}}=\bar{k}_{a}\hat{\xi}^{a}$ with

$\displaystyle\bar{x}_{a}$
$\displaystyle=\coth(2r)x_{a}+\frac{1}{\operatorname{sh}2r}\left(J^{\intercal}\right)_{a}{}^{c}k_{c}\,,$
(63) $\displaystyle\bar{k}_{a}$
$\displaystyle=-\coth(2r)k_{a}+\frac{1}{\operatorname{sh}2r}\left(J^{\intercal}\right)_{a}{}^{c}x_{c}\,.$
(64)

####  2.5.1 Partner modes yield standard form for $G$

This definition yields a normalized mode
$\Omega^{ab}\bar{x}_{a}\bar{k}_{b}=1$, which commutes with the original mode
and is correlated exactly as predicted by (51), _i.e._ , as is straightforward
to check using the identities $(J^{\intercal})^{2}=-\mathbb{1}$ and
$G\Omega^{-1}G=-\Omega$, we have

$\displaystyle G^{ab}\bar{x}_{a}\bar{x}_{b}$
$\displaystyle=G^{ab}\bar{k}_{a}\bar{k}_{b}=\operatorname{ch}2r\,,$ (65)
$\displaystyle\ G^{ab}x_{a}\bar{x}_{b}$
$\displaystyle=-G^{ab}k_{a}\bar{k}_{b}=\operatorname{sh}2r\,,$ (66)
$\displaystyle G^{ab}k_{a}\bar{x}_{b}$
$\displaystyle=G^{ab}x_{a}\bar{k}_{b}=G^{ab}\bar{x}_{a}\bar{k}_{b}=0\,.$ (67)

Note that repeating the partner mode construction gives back the original mode
again, _i.e._ , the partner mode’s partner is the original mode.

There is a close relationship between the correlation of other modes with the
original mode, and the commutator of those other modes with the partner mode:
Let $x^{\prime}_{a}$ and $k^{\prime}_{a}$ define another mode in the system,
_i.e._ , $\Omega^{ab}x^{\prime}_{a}k^{\prime}_{b}=1$. Then

$\displaystyle\Omega^{ab}x^{\prime}_{a}\bar{x}_{b}$
$\displaystyle=\coth(2r)\Omega^{ab}x^{\prime}_{a}x_{b}+\frac{G^{ad}x^{\prime}_{a}k_{d}}{\operatorname{sh}2r}\,,$
(68) $\displaystyle\Omega^{ab}k^{\prime}_{a}\bar{x}_{b}$
$\displaystyle=\coth(2r)\Omega^{ab}k^{\prime}_{a}x_{b}+\frac{G^{ad}k^{\prime}_{a}k_{d}}{\operatorname{sh}2r}\,,$
(69) $\displaystyle\Omega^{ab}x^{\prime}_{a}\bar{k}_{b}$
$\displaystyle=-\coth(2r)\Omega^{ab}x^{\prime}_{a}k_{b}+\frac{G^{ad}x^{\prime}_{a}x_{d}}{\operatorname{sh}2r}\,,$
(70) $\displaystyle\Omega^{ab}x^{\prime}_{a}\bar{x}_{b}$
$\displaystyle=-\coth(2r)\Omega^{ab}k^{\prime}_{a}k_{b}+\frac{G^{ad}k^{\prime}_{a}x_{d}}{\operatorname{sh}2r}\,.$
(71)

This means that if the mode commutes with the original mode
$\Omega(x,x^{\prime})=\Omega(x,k^{\prime})=\Omega(k,x^{\prime})=\Omega(k,k^{\prime})=0$,
then the commutator with the partner mode is proportional to the correlation
with the original mode. In particular, it follows that with respect to a basis
of modes which contains the original mode and the partner mode as its first
two modes, the covariance matrix $G^{ab}$ takes precisely the standard form.

$\displaystyle
G\equiv\left(\begin{array}[]{cc|cc|c}\operatorname{ch}{2r}&0&\operatorname{sh}{2r}&0&0\\\
0&\operatorname{ch}{2r}&0&-\operatorname{sh}{2r}&0\\\
\hline\cr\operatorname{sh}{2r}&0&\operatorname{ch}{2r}&0&0\\\
0&-\operatorname{sh}{2r}&0&\operatorname{ch}{2r}&0\\\ \hline\cr
0&0&0&0&\mathbb{1}_{N-2}\end{array}\right)$ (77)

This shows that $\ket{G,z}$ is in a product state consisting of an entangled
two-mode state (between original modes and its partner) and another Gaussian
state describing the rest of the system. Note that the displacement vector
$z^{a}$ does not affect the entanglement structure of the state.

####  2.5.2 Unsqueezing the partner modes

From the matrix representation of $G$ above we see that partner modes are in a
pure two-mode squeezed state. Hence, they can be viewed as arising from
squeezing two modes which are in a pure product state with each other. We
obtain these two modes by acting with the inverse squeezing transformation on
the two partner modes $(x_{a}\hat{\xi}^{a},k_{a}\hat{\xi}^{a})$ and
$(\bar{x}_{a}\hat{\xi}^{a},\bar{k}_{a}\hat{\xi}^{a})$. This yields the modes
$(y_{a}\hat{\xi}^{a},l_{a}\hat{\xi}^{a})$ and
$(z_{a}\hat{\xi}^{a},m_{a}\hat{\xi}^{a})$, given by

$\displaystyle y_{a}$
$\displaystyle=\operatorname{ch}(r)x_{a}-\operatorname{sh}(r)\bar{x}_{a}=\frac{x_{a}-\left(J^{\intercal}\right)_{a}{}^{c}k_{c}}{2\operatorname{ch}r}\,,$
(78) $\displaystyle l_{a}$
$\displaystyle=\operatorname{ch}(r)k_{a}+\operatorname{sh}(r)\bar{k}_{a}=\frac{k_{a}+\left(J^{\intercal}\right)_{a}{}^{c}x_{c}}{2\operatorname{ch}r}\,,$
(79) $\displaystyle z_{a}$
$\displaystyle=-\operatorname{sh}(r)x_{a}+\operatorname{ch}(r)\bar{x}_{a}=\frac{x_{a}+\left(J^{\intercal}\right)_{a}{}^{c}k_{c}}{2\operatorname{sh}r}\,,$
(80) $\displaystyle m_{a}$
$\displaystyle=\operatorname{sh}(r)k_{a}+\operatorname{ch}(r)\bar{k}_{a}=\frac{-k_{a}+\left(J^{\intercal}\right)_{a}{}^{c}x_{c}}{2\operatorname{sh}r}\,.$
(81)

Since the inverse squeezing is a symplectic transformation these _unsqueezed_
modes form a normalized and commuting basis of the subsystem spanned by the
two partner modes. In particular, with respect to them the covariance matrix
is represented by the identity matrix $G^{ab}\equiv\mathbb{1}$.

##  3 Entanglement extraction from quadratic bosonic systems

In this section, we present the central result of this paper for bosonic
systems. We begin by defining the entanglement extraction procedure in Section
3.1. Our proof then follows three steps:

  1. 1.

Splitting the problem into two parts
In section 3.2, we show that the problem of finding the minimal energy cost
can be broken up into first solving the problem for a system of two modes, and
then optimizing how to embed these two modes used for entanglement extraction
into the larger system.

  2. 2.

Solving the two-mode problem
In Section 3.3, we solve the two mode problem, _i.e._ , we answer the
question: What is the minimal energy cost $\Delta E_{\min}$ of extracting two
partner modes with entanglement entropy $\Delta S$ from a system consisting of
two modes with a quadratic Hamiltonian? For this, we first look at the case of
a degenerate two-mode Hamiltonian in Section 3.3.1, before tackling the
general case in Section 3.3.2. The result is a general formula for the minimal
energy cost as a function of the Hamiltonian’s excitation energies
$\epsilon_{1}$ and $\epsilon_{2}$.

  3. 3.

Solving the embedding problem
In Section 3.4, we derive how to select modes from a large system that
optimize the energy cost. Not surprisingly, the lowest energy cost is achieved
only if the two modes are chosen from the subsystem spanned by the two lowest
energy eigenmodes of the Hamiltonian. Based on this insight, we construct the
optimal partner modes and the associated excitation energies in the subsystem
explicitly. Furthermore, we generalize the expressions and constructions to
scenarios where the two partner modes are restricted to lie in a different
two-mode subsystem.

###  3.1 Pure state entanglement extraction from bosonic modes

The general framework for entanglement extraction consists of two target
systems and the source system $S$. Initially, the total state is assumed to be
a product state

$\displaystyle\rho_{i}=\sigma_{1}\otimes\tau_{2}\otimes\ket{\Psi}\bra{\Psi}_{S}\,,$
(82)

where $\ket{\Psi}_{S}$ is the initial state of the source system, and
$\sigma_{1}$ and $\tau_{2}$ are the initial states of the first and second
target system, respectively. Without loss of generality, the interaction
between the source and target systems can be modeled by two unitary operations
$U_{1}$ and $U_{2}$, which each couple one of the targets to the source. The
final state of the two target systems is then

$\displaystyle\rho_{12}$
$\displaystyle=\operatorname{Tr}_{S}\left(U_{1}U_{2}\rho_{i}U_{2}^{\dagger}U_{1}^{\dagger}\right)\,.$
(83)

In general, the entanglement content of this state depends on the type of
source and target system and their realizable couplings. In particular, it
also depends on the resources such as the energy which is available to
implement the extraction.

The concrete type of source system we study in this section, is a system
composed of bosonic modes with a quadratic Hamiltonian. This system is
supposed to be in its ground state initially, _i.e._ , the state
$\ket{\Phi}_{S}$ is Gaussian. Consequently, its entanglement structure can be
analyzed using the formalism of partner modes introduced before. The two
target modes are assumed to be single bosonic modes.

Figure 1: General framework of entanglement extraction from a system of
bosonic modes. Two target modes get entangled by swapping their states with a
pair of entangled modes in the source system.

The model of interaction we use is that the target modes, with quadratures
$(\hat{Q}_{1},\hat{P}_{1})$ and $(\hat{Q}_{2},\hat{P}_{2})$, each couple to
one mode inside the source system, with quadratures
$(\hat{Q}_{A},\hat{P}_{A})$ and $(\hat{Q}_{B},\hat{P}_{B})$. The unitary
couplings $U_{1}$ and $U_{2}$ are assumed to swap the states of the target
modes and the source modes, _i.e._ ,

$\displaystyle\begin{split}U_{1}^{\dagger}\hat{Q}_{1}U_{1}&=\hat{Q}_{A}\,,\quad
U_{1}^{\dagger}\hat{Q}_{A}U_{1}=\hat{Q}_{1}\,,\\\
U_{1}^{\dagger}\hat{P}_{1}U_{1}&=\hat{P}_{A}\,,\quad\,\,U_{1}^{\dagger}\hat{P}_{A}U_{1}=\hat{P}_{1}\,,\end{split}$
(84)

and $U_{2}$ swaps $(\hat{Q}_{2},\hat{P}_{2})$ and $(\hat{Q}_{B},\hat{P}_{B})$,
accordingly.

The only restriction imposed on the source modes is that they commute
according to

$\displaystyle\left[{\hat{Q}_{A}},{\hat{Q}_{B}}\right]\\!=\\!\left[{\hat{Q}_{A}},{\hat{P}_{B}}\right]\\!=\\!\left[{\hat{P}_{A}},{\hat{Q}_{B}}\right]\\!=\\!\left[{\hat{P}_{A}},{\hat{P}_{B}}\right]\\!=\\!0.$
(85)

This condition ensures that the entanglement of the two target modes was
preexisting in the system, but not created by the couplings between target
modes and system [31]. Other than that, the modes can be chosen freely inside
the source systems.

In particular, the source modes can be chosen to be partner modes, _i.e._ ,
$\hat{Q}_{B}=\hat{Q}_{\bar{A}}$ and $\hat{P}_{B}=\hat{P}_{\bar{A}}$ in the
notation of Section 2.5. As discussed, this implies that the ground state of
the source system $S$ is a product state between the two modes and the rest of
the source system

$\displaystyle\ket{\Phi}_{S}=\ket{\psi}_{A\bar{A}}\otimes\ket{\phi}_{R}\,,$
(86)

where $\ket{\psi}_{A\bar{A}}$ is a two-mode squeezed state, and
$\ket{\phi}_{R}$ is the ground state of the rest of the system. When the
target modes and the source modes are swapped this puts the total system into
the state

$\displaystyle\rho_{f}=\ket{\psi}\bra{\psi}_{12}\otimes\sigma_{A}\otimes\tau_{\bar{A}}\otimes\ket{\phi}\bra{\phi}_{R}\,,$
(87)

where now the target modes are in the entangled state $\ket{\psi}$, whereas
the two source modes are in the initial states of the target modes. The rest
of the source system is unaffected by the extraction and continues to be in
its ground state $\ket{\phi_{R}}$.

Choosing partner modes in the source system is advantageous for two reasons:
Firstly, it is the only way for the target modes to end up in a pure state. If
the system modes are not partners then after the swap the two target modes are
correlated with the source system, _i.e._ , after tracing out the source
system the partial state of target modes is mixed. Secondly, choosing the
modes as partners maximizes the extracted entanglement. If the second party
chooses a different mode, then the (mixed state) entanglement between the two
modes is never larger than between the first mode and its partner. This is
because the mixed state of the two modes is obtained by tracing out part of a
larger state which included the partner mode.

Motivated by this framework and these considerations we study the following
concise question: If in the ground state of a quadratic Hamiltonian the state
of two partner modes sharing entanglement entropy $\Delta S$ is replaced by a
product state, by at least what amount $\Delta E$ does this increase the
energy expectation value of the system.

While in the following we refer to $\Delta E$ as the energy cost of
entanglement extraction, in a realistic implementation of an extraction
protocol additional energy costs are expected: For one, we do not take into
account the energy content of the target modes outside of the source system.
Furthermore, the implementation of the unitary interaction between targets and
source may require energy as well. However, our reason to focus on the energy
$\Delta E$ injected into the source system only, is that the other additional
energy costs depend on the details and the design of the specific extraction
protocol implemented.

Another interesting aspect is the dynamics of the extraction process: In
general, the source modes are not eigenmodes of the Hamiltonian, i.e., they
have a non-trivial time evolution and their physical localization in the
source system may change over time. This dynamic needs to be taken into
account in the design of the interaction of source and target modes. In the
scope of this work, we assume that the swap of target and source modes happens
within a time much shorter than the time scale of the time evolution of the
target modes. Hence we neglect the dynamics and consider the interaction of
target a The vale modes to take part instantaneously. idity of this ind
sourcdealization has to be reviewed with respect to any given implementation
of an extraction protocol. We return to a discussion of this point at the end
of the article.

###  3.2 Energy cost of entanglement extraction

The energy cost of entanglement extraction, as introduced above, is given by
the difference in expectation value of the source system’s Hamiltonian before
and after the extraction, _i.e._ ,

$\displaystyle\Delta
E=\operatorname{Tr}\left(\sigma_{A}\otimes\tau_{\bar{A}}\otimes\ket{\phi}\bra{\phi}_{R}\hat{H}-\ket{\Psi}\bra{\Psi}_{S}\hat{H}\right)\,.$
(88)

The Hamiltonian
$\hat{H}=\frac{1}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}+f_{a}\hat{\xi}^{a}$ of
the source system may be coupling all different modes of the system.
Nevertheless, for the calculation of $\Delta E$, only the part acting on the
two selected partner modes $A$ and $B=\bar{A}$ matter.

To see this, we split up the Hamiltonian into a sum of three terms

$\displaystyle\hat{H}=\hat{H}_{A\bar{A}}+\hat{H}_{R}+\hat{H}_{A\bar{A},R}$
(89)

where $\hat{H}_{A\bar{A}}$ acts only on the two partner modes, $\hat{H}_{R}$
acts only on the rest of the source system’s modes, and $\hat{H}_{A\bar{A},R}$
is the part which couples the partner modes and the rest of the system.111Put
differently, the operators $\hat{H}_{A\bar{A}}$ and $\hat{H}_{R}$ are the
restrictions of $\hat{H}$ onto the subspace spanned by the partner modes and
the rest of the system, respectively.
$\hat{H}_{A\bar{A},R}:=\hat{H}-\hat{H}_{A\bar{A}}-\hat{H}_{R}$ then contains
all terms in $\hat{H}$ that act on both $A\bar{A}$ and $R$ simultaneously.

The summand $\hat{H}_{R}$ does not contribute to the energy cost $\Delta E$,
because the rest of the system is in the same state before and after the
extraction. Also the part $\hat{H}_{A\bar{A},R}$ does not contribute to
$\Delta E$, because both the initial state (86) of the system and its final
state (87) are product states with respect to the partition $A\bar{A}\oplus
V$, hence

$\displaystyle\operatorname{Tr}\left(\sigma_{A}\otimes\tau_{\bar{A}}\otimes\ket{\phi}\bra{\phi}_{R}\hat{H}_{A\bar{A},R}\right)=0\,.$
(90)

Therefore, the energy cost of entanglement extraction is determined by
$\hat{H}_{A\bar{A}}$ acting on the partner modes, only:

$\displaystyle\Delta
E=\operatorname{Tr}\left(\sigma_{A}\otimes\tau_{\bar{A}}\hat{H}_{A\bar{A}}\right)-\operatorname{Tr}\left(\ket{\psi}\bra{\psi}_{A\bar{A}}\hat{H}_{A\bar{A}}\right)\,.$
(91)

From this follows which initial states of the target modes minimize the energy
cost: Because

$\displaystyle\operatorname{Tr}\left(\sigma_{A}\otimes\tau_{\bar{A}}\hat{H}_{A\bar{A}}\right)=\operatorname{Tr}\left(\sigma_{A}\hat{H}_{A}\right)+\operatorname{Tr}\left(\tau_{\bar{A}}\hat{H}_{\bar{A}}\right)\,,$
(92)

the initial states $\sigma$ and $\tau$ of the target modes need to be the
ground states of the restrictions $\hat{H}_{A}$ and $\hat{H}_{\bar{A}}$ of the
Hamiltonian $\hat{H}$ onto the partner modes. These states are pure Gaussian
states.

Equation (91) effectively simplifies finding the minimal energy cost to a two-
mode problem: If we know the minimal energy cost for the extraction of an
amount $\Delta S$ of entanglement from a source system consisting of two
modes, then applying this result to $\hat{H}_{A\bar{A}}$ and optimizing over
different choices of $A$ yields the minimal cost for extracting $\Delta S$
from a larger system of $N$ modes.

###  3.3 Energy cost for two bosonic modes

In this section, we analyze entanglement extraction from a source system that
consists of only two bosonic modes and has a quadratic Hamiltonian. We show
that the minimal energy cost for the extraction from a pair of partner modes
with a fixed amount of entanglement is a function only of the Hamiltonian’s
excitation energies. Our derivation is constructive in the sense that we
derive the symplectic transformation which maps the eigenmodes of the
Hamiltonian to a pair of partner modes which minimize the energy cost.

Since the Hamiltonian is quadratic, we can choose a basis
$\hat{\xi}^{a}\equiv(\hat{q}_{1},\hat{p}_{1},\hat{q}_{2},\hat{p}_{2})$ such
that the quadratic part of the Hamiltonian $h_{ab}$ and the covariance matrix
$G^{ab}$ of its ground state are represented by diagonal matrices $\mathbf{h}$
and $\mathbf{G}$ given by

$\displaystyle h_{ab}$ $\displaystyle\stackrel{{\scriptstyle
1,2}}{{\equiv}}\mathbf{h}=\left(\begin{array}[]{cccc}\epsilon_{1}&0&0&0\\\
0&\epsilon_{1}&0&0\\\ 0&0&\epsilon_{2}&0\\\
0&0&0&\epsilon_{2}\end{array}\right)\,,$ (97) $\displaystyle G^{ab}$
$\displaystyle\stackrel{{\scriptstyle
1,2}}{{\equiv}}\mathbf{G}=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&1&0&0\\\
0&0&1&0\\\ 0&0&0&1\end{array}\right)\,.$ (102)

These modes $1$ and $2$ are not entangled and, therefore, cannot be used as
modes $A$ and $\bar{A}$ for entanglement extraction. Instead, we need to find
a pair of partner modes
$(\hat{Q}_{A},\hat{P}_{A},\hat{Q}_{\bar{A}},\hat{P}_{\bar{A}})$ with respect
to which $G^{ab}$ is represented in standard form (77) given by

$\displaystyle G^{ab}\stackrel{{\scriptstyle
A,\bar{A}}}{{\equiv}}\mathbf{\tilde{G}}=\begin{pmatrix}\operatorname{ch}{2r}&0&\operatorname{sh}{2r}&0\\\
0&\operatorname{ch}{2r}&0&-\operatorname{sh}{2r}\\\
\operatorname{sh}{2r}&0&\operatorname{ch}{2r}&0\\\
0&-\operatorname{sh}{2r}&0&\operatorname{ch}{2r}\end{pmatrix}$ (103)

for a fixed $r>0$, which we choose based on the amount of entanglement that we
want to extract. There is an infinite number of choices for such $A$ and
$\bar{A}$, but we will identify those that minimize the energy cost. This is,
we find the symplectic transformations which map the energy eigenmodes $1$ and
$2$ to those partner modes $A$ and $\bar{A}$ which are optimal for
entanglement extraction.

The two-mode squeezing transformation

$\displaystyle\mathbf{M}_{r}=\left(\begin{array}[]{cccc}\operatorname{ch}{r}&0&\operatorname{sh}{r}&0\\\
0&\operatorname{ch}{r}&0&-\operatorname{sh}{r}\\\
\operatorname{sh}{r}&0&\operatorname{ch}{r}&0\\\
0&-\operatorname{sh}{r}&0&\operatorname{ch}{r}\end{array}\right)$ (108)

brings $\mathbf{G}$ into the form from (103) via

$\displaystyle\mathbf{\tilde{G}}=\mathbf{M}_{r}\mathbf{G}\mathbf{M}^{\intercal}_{r}\,,$
(109)

_i.e._ , it maps the eigenmodes to a pair of partner modes with squeezing
parameter $r$. Under the symplectic transformation $\mathbf{M}_{r}$, the
Hamiltonian matrix $\mathbf{h}$ transforms under the action of the inverse
transform $\mathbf{M}^{-1}_{r}=\mathbf{M}_{-r}$, _i.e._ ,

$\displaystyle\mathbf{\tilde{h}}=\mathbf{M}_{-r}^{\intercal}\mathbf{h}\mathbf{M}_{-r}\,.$
(110)

However, $\mathbf{M}_{r}$ is not the only transformation doing this: We can
compose the squeezing $\mathbf{M}_{r}$ with any symplectic transformation
$\mathbf{N}$ which leaves the covariance matrix invariant, _i.e._ ,
$\mathbf{G}=\mathbf{N}\mathbf{G}\mathbf{N}^{\intercal}$. This yields the most
general form of symplectic transformation
$\mathbf{T}_{r}=\mathbf{M}_{r}\mathbf{N}$ that transforms the modes $1$ and
$2$ into $A$ and $\bar{A}$ for fixed $r$. All of these choices share the same
amount of entanglement across $A$ and $\bar{A}$ because the covariance matrix
takes the same form
$\mathbf{\tilde{G}}=\mathbf{T}_{r}\mathbf{G}\mathbf{T}^{\intercal}_{r}$.
However, the Hamiltonian matrix $\mathbf{h}$ is in general not left invariant
by $\mathbf{N}$, _i.e._ ,
$\mathbf{h}\neq(\mathbf{N}^{-1})^{\intercal}\mathbf{h}\mathbf{N}^{-1}$ and
thus the energy cost will be different.

####  3.3.1 Degenerate two-mode Hamiltonian
$(\epsilon_{1}\\!=\\!\epsilon_{2})$

We first consider the case where $\epsilon:=\epsilon_{1}=\epsilon_{2}$. This
case is the simplest to solve because the matrix representations $\mathbf{h}$
and $\mathbf{G}$ are proportional to each other in the eigenmode basis and,
thus, in all bases. Indeed, any symplectic transformation $\mathbf{N}$ that
leaves $\mathbf{G}$ invariant, _i.e._ ,
$\mathbf{N}\mathbf{G}\mathbf{M}^{\intercal}=\mathbf{G}$, also leaves
$\mathbf{h}$ invariant, _i.e._ ,
$\mathbf{h}=(\mathbf{N}^{-1})^{\intercal}\mathbf{h}\mathbf{N}^{-1}$.
Therefore, for any pair of partner modes $A$ and $\bar{A}$ where the
covariance matrix takes the standard form $\mathbf{\tilde{G}}$ for fixed $r$,
we find

$\displaystyle\mathbf{\tilde{h}}=\left(\begin{array}[]{cc|cc}\epsilon\operatorname{ch}{2r}&0&-\epsilon\operatorname{sh}{2r}&0\\\
0&\epsilon\operatorname{ch}{2r}&0&\epsilon\operatorname{sh}{2r}\\\
\hline\cr-\epsilon\operatorname{sh}{2r}&0&\epsilon\operatorname{ch}{2r}&0\\\
0&\epsilon\operatorname{sh}{2r}&0&\epsilon\operatorname{ch}{2r}\end{array}\right).$
(115)

As discussed following equation (92), the entanglement extraction swaps the
state of the two modes against the product state
$\sigma_{A}\otimes\tau_{\bar{A}}$ of two Gaussian states $\sigma$ and $\tau$.
Hence, denoting $\sigma_{A}=\ket{G_{A}^{\prime}}\bra{G_{A}^{\prime}}$ and
$\tau_{\bar{A}}=\ket{G_{\bar{A}}^{\prime}}\bra{G_{\bar{A}}^{\prime}}$, with
respect to partner mode basis, the covariance matrix after the extraction
takes the block diagonal form

$\displaystyle\mathbf{\tilde{G}^{\prime}}=\left(\begin{array}[]{c|c}\mathbf{\tilde{G}}^{\prime}_{A}&0\\\\[2.84526pt]
\hline\cr\\\\[-8.53581pt]
0&\mathbf{\tilde{G}}^{\prime}_{\bar{A}}\end{array}\right).$ (119)

To minimize the energy cost, we choose the initial states of the target modes
to be the ground states of the restricted single-mode Hamiltonians, _i.e._ ,
$\mathbf{\tilde{G}}^{\prime}_{A}=\mathbf{\tilde{G}}^{\prime}_{\bar{A}}=\mathbb{1}_{2}$.
In terms of this new covariance matrix, the energy cost (91) of the
entanglement extraction is

$\displaystyle\Delta E$
$\displaystyle=\frac{1}{4}\operatorname{Tr}\left(\mathbf{\tilde{h}}^{\intercal}\mathbf{\tilde{G}^{\prime}}\right)-\frac{1}{4}\operatorname{Tr}\left(\mathbf{h}^{\intercal}{\mathbf{G}}\right)$
(120)
$\displaystyle=\frac{1}{4}\operatorname{Tr}\left(\mathbf{\tilde{h}}_{A}\mathbf{\tilde{G}}_{A}+\mathbf{\tilde{h}}_{\bar{A}}\mathbf{\tilde{G}}_{\bar{A}}\right)-\epsilon,$
(121)

where
$\mathbf{\tilde{h}}_{A}=\mathbf{\tilde{h}}_{\bar{A}}=\epsilon\operatorname{ch}{2r}\mathbb{1}_{2}$
are the diagonal blocks of $\mathbf{\tilde{h}}$ above.

With these initial states, the minimal energy cost for a pure two-mode
squeezed state with entanglement entropy $\Delta
S=s_{b}(\operatorname{ch}{2r})$ (see (58)) is

$\displaystyle\Delta E_{\mathrm{min}}$
$\displaystyle=2\epsilon\operatorname{sh}^{2}{r}\,.$ (122)

The final state of the source system after the swap is not an energy
eigenstate anymore. Its energy variance is

$\displaystyle\Sigma_{E}=\frac{\sqrt{\operatorname{Tr}\left(\mathbf{\tilde{h}}\mathbf{\tilde{G}}^{\prime}\mathbf{\tilde{h}}\mathbf{\tilde{G}}^{\prime}+\mathbf{\tilde{h}}\mathbf{\Omega}\mathbf{\tilde{h}}\mathbf{\Omega}\right)}}{4}=\frac{\epsilon\operatorname{sh}{2r}}{\sqrt{2}}\,.$
(123)

We notice that for the degenerate Hamiltonian with
$\epsilon_{1}=\epsilon_{2}$, the energy cost and energy variance is actually
the same for all partner modes with entanglement entropy $\Delta S$. In the
next section, we will see that this is different for the case
$\epsilon_{1}<\epsilon_{2}$.

####  3.3.2 General two-mode Hamiltonian ($\epsilon_{1}\\!<\\!\epsilon_{2}$)

The case of a general Hamiltonian, with $\epsilon_{1}\leq\epsilon_{2}$, is
more involved to analyze because there exist symplectic transformations which
leave the covariance matrix of the state invariant but do change the
Hamiltonian matrix and vice versa.

Thus, before we apply the squeezing transformation $\mathbf{M}_{r}$, we can
apply another symplectic transformation
$\mathbf{N}=\exp\left(\theta\mathbf{K}\right)$ with

$\displaystyle\mathbf{K}$
$\displaystyle=\left(\begin{array}[]{cccc}0&0&\cos{\phi}&-\sin{\phi}\\\
0&0&\sin{\phi}&\cos{\phi}\\\ -\cos{\phi}&-\sin{\phi}&0&0\\\
\sin{\phi}&-\cos{\phi}&0&0\end{array}\right)\,,$ (128)

which leaves $\mathbf{G}$ invariant, but possibly changes $\mathbf{h}$. Above
general form can be derived by writing $\mathbf{N}$ as a matrix exponential
$\mathbf{N}=\exp\left(\theta\mathbf{K}\right)$ of a generator $\mathbf{K}$
subject to the constraints
$\mathbf{K}\mathbf{G}+\mathbf{G}\mathbf{K}^{\intercal}=0$ (leaving
$\mathbf{G}$ invariant) and
$\mathbf{K}\mathbf{\Omega}+\mathbf{\Omega}\mathbf{K^{\intercal}}=0$ (being
symplectic).

While $\mathbf{N}=\exp\left(\theta\mathbf{K}\right)$ leaves $\mathbf{G}$
invariant, it changes $\mathbf{h}$. Consequently,
$\mathbf{T}_{r}=\mathbf{M}_{r}\mathbf{N}$ provides the most general
transformation that turns the covariance matrix from (102) into
$\mathbf{\tilde{G}}=\mathbf{T}_{r}\mathbf{G}\mathbf{T}_{r}^{\intercal}$ in the
form of (103). Under this transformation, the Hamiltonian matrix transforms to

$\displaystyle\mathbf{\tilde{h}}=(\mathbf{T}_{r}^{-1})^{\intercal}\mathbf{h}\mathbf{T}_{r}^{-1}=\left(\begin{array}[]{c|c}\mathbf{h}_{A}&\mathbf{h}_{X}\\\\[2.84526pt]
\hline\cr\\\\[-11.38109pt]
\mathbf{h}_{X}^{\intercal}&\mathbf{h}_{\bar{A}}\end{array}\right)\,,$ (132)

where the individual blocks are given by

$\displaystyle\mathbf{h}_{A}$
$\displaystyle=\left(\begin{array}[]{cc}\epsilon_{+}\operatorname{ch}{2r}-\epsilon_{-}(\cos{2\theta}+\sin{2\theta}\cos{\phi}\operatorname{sh}{2r})&-\epsilon_{-}\sin{2\theta}\sin{\phi}\operatorname{sh}{2r}\\\
-\epsilon_{-}\sin{2\theta}\sin{\phi}\operatorname{sh}{2r}&\epsilon_{+}\operatorname{ch}{2r}-\epsilon_{-}(\cos{2\theta}-\sin{2\theta}\cos{\phi}\operatorname{sh}{2r})\end{array}\right)\,,$
(135) $\displaystyle\mathbf{h}_{\bar{A}}$
$\displaystyle=\left(\begin{array}[]{cc}\epsilon_{+}\operatorname{ch}{2r}+\epsilon_{-}(\cos{2\theta}-\sin{2\theta}\cos{\phi}\operatorname{sh}{2r})&\epsilon_{-}\sin{2\theta}\sin{\phi}\operatorname{sh}{2r}\\\
\epsilon_{-}\sin{2\theta}\sin{\phi}\operatorname{sh}{2r}&\epsilon_{+}\operatorname{ch}{2r}+\epsilon_{-}(\cos{2\theta}+\sin{2\theta}\cos{\phi}\operatorname{sh}{2r})\end{array}\right)\,,$
(138) $\displaystyle\mathbf{h}_{X}$
$\displaystyle=\left(\begin{array}[]{cc}-\epsilon_{+}\operatorname{sh}{2r}+\epsilon_{-}\sin{2\theta}\cos{\phi}\operatorname{ch}{2r}&-\epsilon_{-}\sin{2\theta}\sin{\phi}\operatorname{ch}{2r}\\\
\epsilon_{-}\sin{2\theta}\sin{\phi}\operatorname{ch}{2r}&\epsilon_{+}\operatorname{sh}{2r}+\epsilon_{-}\sin{2\theta}\cos{\phi}\operatorname{ch}{2r}\end{array}\right)$
(141)

with $\epsilon_{\pm}=\frac{1}{2}(\epsilon_{2}\pm\epsilon_{1})$. The
entanglement extraction swaps the state of the two modes against the initial
state of the target modes. We denote its covariance matrix with respect $A$
and $\bar{A}$ by $\mathbf{\tilde{G}^{\prime}}$ as in (119). The energy cost is
then given by

$\displaystyle\Delta E$
$\displaystyle=\frac{1}{4}\operatorname{Tr}\left(\mathbf{\tilde{h}}_{A}\mathbf{\tilde{G}^{\prime}}_{A}+\mathbf{\tilde{h}}_{\bar{A}}\mathbf{\tilde{G}^{\prime}}_{\bar{A}}\right)-\frac{\epsilon_{1}+\epsilon_{2}}{2}\,.$
(142)

To minimize the energy cost we choose the initial states to be the ground
states of the single-mode Hamiltonians of $\mathbf{h_{A}}$ and
$\mathbf{h_{\bar{A}}}$. This yields

$\displaystyle\Delta
E=\frac{1}{2}(\epsilon_{A}+\epsilon_{\bar{A}}-\epsilon_{1}-\epsilon_{2})\,,$
(143)

where the excitation energies of the single-mode Hamiltonians are given by the
symplectic eigenvalues of $\mathbf{h}_{A}$ and $\mathbf{h}_{\bar{A}}$, _i.e._
,

$\displaystyle\epsilon_{A}$
$\displaystyle=\frac{1}{2}\sqrt{4\epsilon_{+}^{2}\operatorname{ch}^{2}{2r}+\epsilon_{-}^{2}\left(3+\cos{4\theta}-(1-\cos{4\theta})\operatorname{ch}{4r}\right)-8\epsilon_{+}\epsilon_{-}\cos{2\theta}\operatorname{ch}{2r}}\,,$
(144) $\displaystyle\epsilon_{\bar{A}}$
$\displaystyle=\frac{1}{2}\sqrt{4\epsilon_{+}^{2}\operatorname{ch}^{2}{2r}+\epsilon_{-}^{2}\left(3+\cos{4\theta}-(1-\cos{4\theta})\operatorname{ch}{4r}\right)+8\epsilon_{+}\epsilon_{-}\cos{2\theta}\operatorname{ch}{2r}}\,.$
(145)

Since the energy cost is independent of the parameter $\phi$, we only need to
compute the minimum with respect to $\theta$, in order to find partner modes
which yield the minimal energy cost of entanglement extraction. $\Delta E$
only depends on $\theta$ through $\epsilon_{A}+\epsilon_{\bar{A}}$ and we find
that the minimum is attained for $\theta=\frac{\pi}{4}+\frac{n\pi}{2}$ with
$n\in\mathbb{Z}$. This value gives
$\epsilon_{A}=\epsilon_{\bar{A}}=\frac{1}{2}\sqrt{\epsilon_{1}^{2}+\epsilon_{2}^{2}+2\epsilon_{1}\epsilon_{2}\operatorname{ch}{4r}}$.
Hence the minimal energy cost for the extraction of a pure squeezed state with
entanglement entropy $\Delta S=S(\operatorname{ch}{2r})$ is given by

$\displaystyle\Delta E_{\mathrm{min}}$
$\displaystyle=\frac{1}{2}\left(\sqrt{\epsilon_{1}^{2}+\epsilon_{2}^{2}+2\epsilon_{1}\epsilon_{2}\operatorname{ch}{4r}}-\epsilon_{1}-\epsilon_{2}\right)\,.$
(146)

The energy variance of the final system state is

$\displaystyle\Sigma_{E}$
$\displaystyle=\sqrt{\frac{\epsilon_{1}^{2}\epsilon_{2}^{2}(1+\operatorname{ch}{4r})}{\epsilon_{1}^{2}+\epsilon_{2}^{2}+2\epsilon_{1}\epsilon_{2}\operatorname{ch}{4r}}}\,\operatorname{sh}{2r}\,.$
(147)

Figure 2: Minimal energy cost $\Delta E_{\mathrm{min}}$ as a function of the
extracted entanglement entropy $\Delta S$ from a bosonic, quadratic two-mode
Hamiltonian with excitation energies $\epsilon_{1}$ and $\epsilon_{2}$.

The minimal energy cost has a couple of interesting properties, which are seen
in Figure 2. First, we note that the expression above correctly reduces to the
degenerate case for $\epsilon_{2}=\epsilon_{1}$. For large amounts of
extracted entanglement, the energy cost grows exponentially as

$\displaystyle\Delta E_{\mathrm{min}}\sim\mathrm{e}^{2r}\sim e^{\Delta
S}\quad\text{as}\quad r\to\infty\,,$ (148)

where we used the asymptotics $\Delta S\sim 2r$ as $r\to\infty$. For small
$r$, we find

$\displaystyle\Delta
E_{\mathrm{min}}\sim\frac{4\epsilon_{1}\epsilon_{2}}{\epsilon_{1}+\epsilon_{2}}r^{2}+\mathcal{O}(r^{4})\quad\text{as}\quad
r\to 0\,,$ (149)

which should be seen in the context of the asymptotics $\Delta S\sim
r^{2}\left(1-\log\left(r^{2}\right)\right)$ as $r\to 0$.

Furthermore, the energy cost is a monotonously increasing function in both
excitation energies. This means, for example that $\Delta E_{\mathrm{min}}$
increases when $\epsilon_{2}$ is increased and $\epsilon_{1}$ is kept fixed.
However, there exists an upper bound on the energy cost in this case, which is
determined by the lower excitation energy

$\displaystyle\lim_{\epsilon_{2}\to\infty}\Delta
E_{\mathrm{min}}=\epsilon_{1}\left(\operatorname{sh}2r\right)^{2}\,.$ (150)

Comparing the minimal energy cost for a non-degenerate Hamiltonian $\Delta
E_{\mathrm{min},\epsilon_{2}>\epsilon_{1}}$ with the cost for a degenerate
Hamiltonian $\Delta E_{\epsilon_{2}=\epsilon_{1}}$, we find that the ratio
between them asymptotes to

$\displaystyle\lim_{r\to\infty}\frac{\Delta
E_{\mathrm{min},\epsilon_{2}>\epsilon_{1}}}{\Delta
E_{\epsilon_{2}=\epsilon_{1}}}=\sqrt{\frac{\epsilon_{2}}{\epsilon_{1}}}\,.$
(151)

Another interesting figure of merit is to look at the ratio between a given
amount of extracted entanglement and the minimal energy cost necessary for its
extraction $\Delta S/\Delta E_{\mathrm{min}}$. Since $\Delta
S=S(\operatorname{ch}{2r})\sim(1-2\log(r))r^{2}$ as $r\to 0$, the ratio
diverges in the limit of small amounts of extracted entanglement, _i.e._ , we
have

$\displaystyle\lim_{r\to 0}\frac{\Delta S}{\Delta E_{\mathrm{min}}}=\infty\,.$
(152)

As the extracted entanglement increases the ratio is strictly decreasing, and
in the limit of large entanglement $r\to\infty$ it decays exponentially. This
means that it requires more energy to extract an amount of entanglement
$\Delta S$ from a system once, than twice time the energy required to extract
half the amount $\Delta S/2$.

Figure 3:  We show the energy variance $\Sigma_{E}$ as a function of the
extracted entanglement entropy $\Delta S$ from a bosonic, quadratic two-mode
Hamiltonian with excitation energies $\epsilon_{1}$ and $\epsilon_{2}$

###  3.4 Optimal modes in large source systems

Having understood how to minimize the energy cost of entanglement extraction
from a two-mode source system, allows us to also understand how to optimally
extract entanglement from large source systems.

The partner modes from which entanglement can be extracted for the least
possible energy cost are linear combinations of the two lowest energy
eigenmodes of the system only. No other choice of partner modes can have a
lower energy cost, because the excitation energies $\epsilon_{1}$ and
$\epsilon_{2}$ of the restricted two-mode Hamiltonian $\hat{H}_{A\bar{A}}$ are
restricted by the excitation energies of $\hat{H}$ by

$\displaystyle\omega_{1}\leq\epsilon_{1}\leq\omega_{N-1}\,,\quad\omega_{2}\leq\epsilon_{2}\leq\omega_{N}$
(153)

as discussed in Appendix B.

These modes can be constructed according to above derivation. Take the two
lowest energy eigenmodes $(\hat{q}_{1},\hat{p}_{1},\hat{q}_{2},\hat{p}_{2})$
of the system and mix them by acting on them with the transformation

$\displaystyle\exp\left(\frac{\pi}{4}\mathbf{K}\right)=\left(\begin{array}[]{cccc}\frac{1}{\sqrt{2}}&0&\frac{\cos\phi}{\sqrt{2}}&-\frac{\sin\phi}{\sqrt{2}}\\\
0&\frac{1}{\sqrt{2}}&\frac{\sin\phi}{\sqrt{2}}&\frac{\cos\phi}{\sqrt{2}}\\\
-\frac{\cos\phi}{\sqrt{2}}&-\frac{\sin\phi}{\sqrt{2}}&\frac{1}{\sqrt{2}}&0\\\
\frac{\sin\phi}{\sqrt{2}}&-\frac{\cos\phi}{\sqrt{2}}&0&\frac{1}{\sqrt{2}}\\\
\end{array}\right)$ (158)

where $\phi$ can be chosen freely. The two obtained modes need to by squeezed
by acting with $\mathbf{M}_{r}$ on them, in order to obtain the globally
optimal partner modes of the source system.

In many scenarios the globally optimal modes may not be accessible to the
experimenter, but they may be restricted to use a particular pair of partner
modes $(\hat{Q}_{A},\hat{P}_{A},\hat{Q}_{\bar{A}},\hat{P}_{\bar{A}})$ as
source modes. In this case, applying results of the previous section to the
two-mode restriction $\hat{H}_{A\bar{A}}$ yields a lower bound on the energy
cost of entanglement extraction, as discussed in Section 3.2.

However, the modes
$(\hat{Q}_{A},\hat{P}_{A},\hat{Q}_{\bar{A}},\hat{P}_{\bar{A}})$ can only
achieve the minimal energy cost $\Delta E_{\mathrm{min}}$ which the spectrum
of $\hat{H}_{A\bar{A}}$ allows for according to (146), if they exactly
correspond to the optimal modes with respect to $\hat{H}_{A\bar{A}}$ as we
constructed them in the derivation of the previous section.

To check if this is the case it is sufficient to check if $\hat{H}_{A\bar{A}}$
acts symmetrically on the two modes $A$ and $\bar{A}$ because this is exactly
what happens for the Hamiltonian matrix in (132) when
$\theta=\frac{\pi}{4}+\frac{n\pi}{2}$ is chosen optimally.

The Hamiltonian matrix of $\hat{H}_{A\bar{A}}$ with respect to the basis
$(\hat{Q}_{A},\hat{P}_{A},\hat{Q}_{\bar{A}},\hat{P}_{\bar{A}})$ is
conveniently calculated using an inner product on $V^{*}$ which is induced by
$h_{ab}$. To do this, we denote

$\displaystyle\hat{Q}_{A}$
$\displaystyle=x_{a}\hat{\xi}^{a}\,,\qquad\hat{P}_{A}=k_{a}\hat{\xi}^{a}\,,$
(159) $\displaystyle\hat{Q}_{\bar{A}}$
$\displaystyle=\bar{x}_{a}\hat{\xi}^{a},\qquad\hat{P}_{\bar{A}}=\bar{k}_{a}\hat{\xi}^{a}\,.$
(160)

and assume, without loss of generality, that the modes are in standard form as
in Section 2.5. With respect to this basis, the quadratic part of
$\hat{H}_{A\bar{A}}$ is represented by

$\displaystyle\left(\begin{array}[]{cc|cc}\check{h}(k,k)&-\check{h}(k,x)&\check{h}(k,\bar{k})&-\check{h}(k,\bar{x})\\\
-\check{h}(x,k)&\check{h}(x,x)&-\check{h}(x,\bar{k})&\check{h}(x,\bar{x})\\\
\hline\cr\check{h}(\bar{k},k)&-\check{h}(\bar{k},x)&\check{h}(\bar{k},\bar{k})&-\check{h}(\bar{k},\bar{x})\\\
-\check{h}(\bar{x},k)&\check{h}(\bar{x},x)&-\check{h}(\bar{x},\bar{k})&\check{h}(\bar{x},\bar{x})\end{array}\right)\,,$
(165)

with the inner product $\check{h}:V^{*}\times V^{*}\to\mathbb{R}$ with

$\displaystyle\check{h}(w,u)={h}_{ab}\Omega^{ac}\Omega^{bd}w_{c}u_{d}={h}_{ab}G^{ac}G^{bd}w_{c}u_{d}\,,$
(166)

induced by $h_{ab}$ of the full Hamiltonian. Denoting
$(J^{\intercal}w)_{a}=(J^{\intercal})_{a}{}^{b}r_{b}$, the inner product
satisfies

$\displaystyle\check{h}(J^{\intercal}w,J^{\intercal}u)=\check{h}(w,u)\,,$
(167)

which just implies that $J^{\intercal}$ is orthogonal with respect to
$\check{h}$ in the same way as $J$ is orthogonal with respect to $h$. In
particular, this implies $\check{h}(w,J^{\intercal}w)=0$, such that the matrix
expressions simplify

$\displaystyle\check{h}(\bar{x},\bar{x})$
$\displaystyle=\check{h}(x,x)+\Delta\,,$ (168)
$\displaystyle\check{h}(\bar{k},\bar{k})$
$\displaystyle=\check{h}(k,k)+\Delta\,,$ (169)
$\displaystyle\check{h}(\bar{x},\bar{k})$ $\displaystyle=-\check{h}(x,k)\,,$
(170) $\displaystyle\check{h}(x,\bar{x})$
$\displaystyle=\frac{\operatorname{ch}(2r)\,\check{h}(x,x)+\check{h}(x,J^{\intercal}k)}{\operatorname{sh}(2r)}\,,$
(171) $\displaystyle\check{h}(k,\bar{k})$
$\displaystyle=-\frac{\operatorname{ch}(2r)\,\check{h}(k,k)\\!+\\!\check{h}(x,J^{\intercal}k)}{\operatorname{sh}(2r)}\,,$
(172) $\displaystyle\check{h}(x,\bar{k})$
$\displaystyle=-\frac{\operatorname{ch}(2r)}{\operatorname{sh}(2r)}\,\check{h}(x,k)=-\check{h}(k,\bar{x})\,,$
(173)

where

$\displaystyle\Delta$
$\displaystyle=\frac{\check{h}(x,x)+\check{h}(k,k)+2\operatorname{ch}(2r)\,\check{h}(x,J^{\intercal}k)}{\operatorname{sh}^{2}(2r)}\,.$
(174)

For $\Delta=0$, the Hamiltonian matrix is symmetric with respect to the two
partner modes. This indicates that the partner modes are in fact optimal modes
inside the two-mode subspace on which $\hat{H}_{A\bar{A}}$ acts and they
achieve the minimal energy cost. Why this is the case is further illuminated
by the following analysis of the case where $\Delta$ does not vanish

To analyze the general case of $\Delta\neq 0$ it is useful to reverse-engineer
the derivation of the previous section. There, we showed that all partner
modes in the subspace of $\hat{H}_{A\bar{A}}$ are connected to the eigenmodes
of $\hat{H}_{A\bar{A}}$ through a transformation
$\mathbf{T}_{r}=\mathbf{M}_{r}\mathbf{N}_{\theta}$ which is the product of an
orthogonal transformation
$\mathbf{N}_{\theta}=\exp\left(\theta\mathbf{K}\right)$ and a squeezing
operation $\mathbf{M}_{r}$.

If we apply the first orthogonal transformation $\mathbf{N}_{\theta}$ to the
energy eigenmodes we obtain a pair of modes which still diagonalize the
covariance matrix $\mathbf{G}$. This pair of modes we also obtain by acting
with the inverse squeezing transformation $\mathbf{M}_{-r}$ on the partner
modes, _i.e._ , they are the unsqueezed modes $(y,l,z,m)$ as defined in (78).

Now we are able to express the matrix representing $\hat{H}_{A\bar{A}}$ in two
different ways. On the one hand, we obtain it starting from its diagonal form
$\mathbf{h}$ in (97) for the energy eigenmodes, and on the other hand we can
express its entries using the form $\check{h}$ and the expressions (78) for
the unsqueezed modes.

$\displaystyle\begin{split}&\left(\mathbf{N}_{\theta}^{-1}\right)^{\intercal}\mathbf{h}\mathbf{N}_{\theta}^{-1}=\\\
&\left(\begin{array}[]{cc|cc}\check{h}(l,l)&-\check{h}(l,y)&\check{h}(l,m)&-\check{h}(l,z)\\\
-\check{h}(y,l)&\check{h}(y,y)&-\check{h}(y,m)&\check{h}(y,z)\\\
\hline\cr\check{h}(m,l)&-\check{h}(m,y)&\check{h}(m,m)&-\check{h}(m,z)\\\
-\check{h}(z,l)&\check{h}(z,y)&-\check{h}(z,m)&\check{h}(z,z)\end{array}\right)\,.\end{split}$
(175)

From this equation the values of the parameters $\phi$ and $\theta$ which
yield the transformation $\mathbf{N}$ from the eigenmodes of
$\hat{H}_{A\bar{A}}$ to the unsqueezed modes of $A$ and $\bar{A}$ as well as
the spectrum of $\hat{H}_{A\bar{A}}$ can be determined. The latter is given by

$\displaystyle\epsilon_{+}$
$\displaystyle=\frac{\check{h}(x,x)+\check{h}(k,k)+\Delta}{2\operatorname{ch}{2r}}\,,$
(176) $\displaystyle\epsilon_{-}$
$\displaystyle=\sqrt{\frac{4\check{h}(x,k)^{2}+\left(\check{h}(x,x)-\check{h}(k,k)\right)^{2}}{4\operatorname{sh}^{2}(2r)}+\frac{\Delta^{2}}{4}}\,,$
(177)

with $\epsilon_{\pm}=\frac{1}{2}(\epsilon_{2}\pm\epsilon_{1})$ as before. The
parameter values are determined by

$\displaystyle\tan\phi$
$\displaystyle=\frac{2\check{h}(x,k)}{\check{h}(x,x)-\check{h}(k,k)}\,,$ (178)
$\displaystyle\tan 2\theta$
$\displaystyle=\frac{\sqrt{4\epsilon_{-}^{2}-\Delta^{2}}}{\Delta}\,,$ (179)

or, alternatively, $\Delta=2\epsilon_{-}\cos 2\theta$.

With these parameters, we can calculate the minimal energy cost $\Delta E$ of
entanglement extraction from a pair of given partner modes
$(x,k\bar{,}x,\bar{k})$ using (143). The energy cost lies in the interval

$\displaystyle\Delta E_{\theta=\pi/4}\leq\Delta E\leq\Delta E_{\theta=0}\,,$
(180)

which is bounded by the minimal energy cost $\Delta E_{\min}$, achieved by
$\theta=\frac{\pi}{4}+\frac{n\pi}{2}$,

$\displaystyle\Delta E_{\theta=\pi/4}=\Delta
E_{\mathrm{min}}=\sqrt{\epsilon_{+}^{2}\operatorname{ch}^{2}{2r}-\epsilon_{-}^{2}\operatorname{sh}^{2}2r}$
(181)

and the energy cost for a pair of partner modes obtained with the least
favorable choice for the parameter $\theta=0+\frac{n\pi}{2}$, which is

$\displaystyle\Delta E_{\theta=0}=2\epsilon_{+}\operatorname{sh}^{2}{r}\,.$
(182)

In fact, the last equality implies an upper bound on the energy cost of
entanglement extraction with arbitrary partner modes of the system. Due to
(153), we have

$\displaystyle\Delta E_{\theta=0}\leq
2\omega_{N}\operatorname{sh}^{2}r=:\Delta E_{\mathrm{max}}\,.$ (183)

This means that for any pair of partner modes in the system the energy cost
$\Delta E$ (see (91) and (92)) of swapping their entangled state against the
product of their one-mode restricted ground states is upper bounded by the
minimal cost for a degenerate Hamiltonian (see (122)) with only the system’s
largest excitation energy $\omega_{N}$ in its spectrum.

Finally, we note that the analysis above allows us to construct from a given
pair of partner modes
$(\hat{Q}_{A},\hat{P}_{A},\hat{Q}_{\bar{A}},\hat{P}_{\bar{A}})$ the pair of
modes which are optimal for extracting entanglement from the two-mode subspace
of $\hat{H}_{A\bar{A}}$. For this we only need to correct for the mismatch
between $\mathbf{N}_{\theta}$ corresponding to the given modes, and
$\mathbf{N}_{\theta=\pi/4}$ which yields the optimal pair. So to obtain the
ideal modes we only need to act with the transformation

$\displaystyle\mathbf{M}_{r}\exp\left((\frac{\pi}{4}-\theta)\mathbf{K}\right)\mathbf{M}_{-r}$
(184)

on the modes $(\hat{Q}_{A},\hat{P}_{A},\hat{Q}_{\bar{A}},\hat{P}_{\bar{A}})$.

##  4 Quadratic fermionic systems and their partner mode construction

We review the most important definitions of quadratic fermionic systems and
fermionic Gaussian states. More detailed reviews can be found in [14, 61, 62],
but we follow mostly the conventions introduced in [67, 59]. In particular,
defining fermionic Gaussian states using linear complex structures emphasizes
the similarities to the bosonic case. The key difference between bosonic and
fermionic Gaussian states is that $\Omega$ and $G$ exchange their roles. From
a group theoretic perspective, this implies that the symplectic group
$\mathrm{Sp}(2N,\mathbb{R})$ preserving $\Omega$ (governing commutation
relations) is now replaced by the orthogonal group $\mathrm{O}(2N)$ preserving
a positive definite metric $G$ (governing anti commutation relations).

###  4.1 Classical and quantum theory

As for bosons, we can formally define a classical phase space
$V\simeq\mathbb{R}^{2N}$ and its dual $V^{*}\simeq\mathbb{R}^{2N}$ for a
system with $N$ fermionic degrees of freedom. In contrast to bosons, $V$ and
$V^{*}$ are now equipped with a positive-definite bilinear form
$G^{ab}:V^{*}\times V^{*}\to\mathbb{R}$ and its inverse $G^{-1}_{ab}:V\times
V\to\mathbb{R}$ satisfying $G^{ac}G^{-1}_{cb}=\delta^{a}{}_{b}$.

We can introduce an orthonormal basis of linear observables $q_{i},p_{i}\in
V^{*}\subset C^{\infty}(V)$, such that $G(p_{i},p_{j})=0$ and
$G(q_{i},q_{j})=G(p_{i},p_{j})=\delta_{ij}$. As in the bosonic case, $q_{i}$
and $p_{i}$ yield a component representation of a vector $v^{a}\in V$ as

$\displaystyle v^{a}\equiv(q_{1}(v),p_{1}(v),\dots,q_{N}(v),p_{N}(v))\,.$
(185)

With respect to this basis, the bilinear form $G^{ab}$ takes the standard form
$G^{ab}\equiv\mathbb{1}$. We use $\xi^{a}:V\to V:v^{a}\mapsto\xi^{a}(v)=v^{a}$
in the same way as for bosons. However, the fermionic classical observables
are not smooth functions on $V$, but rather the Grassmann algebra generated by
$\xi^{a}$, _i.e._ , a formal series of anti-symmetrized powers in $\xi^{a}$.

When we quantized the bosonic system, we used the classical Poisson brackets
encoded in $\Omega^{ab}$. For fermions, we implement the canonical anti-
commutation relations

$\displaystyle[\hat{\xi}^{a},\hat{\xi}^{b}]_{+}=\hat{\xi}^{a}\hat{\xi}^{b}+\hat{\xi}^{b}\hat{\xi}^{a}=\\{\xi^{a},\xi^{b}\\}_{+}=G^{ab}\,,$
(186)

where $\\{\xi^{a},\xi^{b}\\}_{+}=G^{ab}$ represents the classical fermionic
Poisson brackets.

###  4.2 Fermionic Gaussian states

Fermionic Gaussian states appear under various names in the literature,
ranging from fermionic vacua to Slater determinant states. There is no
Gaussian wave function representations of them, but they share many properties
with their bosonic counterparts. They form the set of ground states of
quadratic Hamiltonians and also—in contrast to bosonic Gaussian states—their
higher excited eigenstates. In the context of fermionic Bogoliubov
transformations, they appear as vacua of annihilation operators satisfying
fermionic anti commutation relations. Finally, their $n$-point correlation
functions are completely determined from their 2-point correlation function.

We label a fermionic Gaussian state $\ket{\Omega}$ by its covariance matrix

$\displaystyle\Omega^{ab}=\bra{\Omega}\hat{\xi}^{a}\hat{\xi}^{b}-\hat{\xi}^{b}\hat{\xi}^{a}\ket{\Omega}\,,$
(187)

which is antisymmetric, rather than symmetric as in the bosonic case. Note
that there is no linear displacement for physical fermionic Gaussian states,
_i.e._ , we require $\bra{\Omega}\hat{\xi}^{a}\ket{\Omega}=0$. We can use
$\Omega$ to construct a linear complex structure

$\displaystyle J^{a}{}_{b}=\Omega^{ac}G^{-1}_{cb}=-G^{ac}\Omega^{-1}_{bc}\,,$
(188)

that satisfies $J^{2}=-\mathbb{1}$ just as for bosons. Moreover, we have the
condition
$\frac{1}{2}(\delta^{a}{}_{b}+\mathrm{i}J^{a}{}_{b})\hat{\xi}^{b}\ket{\Omega}=0$,
analogously to (22).

In contrast to bosons, there is no displacement vector $z^{a}$ for fermions.
Hence, fermionic $n$-point functions are directly defined as

$\displaystyle C_{n}^{a_{1}\cdots
a_{n}}=\bra{\Omega}\hat{\xi}^{a_{1}}\dots\hat{\xi}^{a_{n}}\ket{\Omega}\,.$
(189)

Apart from this, the definition of the two-point function is unchanged:

$\displaystyle C_{2}^{ab}=\frac{1}{2}(G^{ab}+\mathrm{i}\Omega^{ab})\,.$ (190)

Wick’s theorem applies in almost the same way

$\displaystyle C^{a_{1}\cdots a_{n}}_{n}$
$\displaystyle=\sum\text{(contractions of $C_{2}^{ab}$)}$ (191)
$\displaystyle=\sigma_{(a_{1},\dots,a_{n})}C^{a_{1}a_{2}}_{2}\cdots
C^{a_{n-1}a_{n}}_{2}+\cdots\,,$ (192)

however, each contraction enters with a positive or negative sign, depending
on the parity $\sigma_{(a_{1},\dots,a_{n})}\in\\{-1,1\\}$ of its index
ordering. This subtlety needs to be taken into account carefully when applying
Wick’s theorem to fermions.

##### Example.

Looking at the ground state $\ket{\Omega}$ of the fermionic oscillator
$\hat{H}=\frac{\mathrm{i}}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}=\mathrm{i}\hat{q}\hat{p}$.
The covariance matrix $\Omega^{ab}$ and the metric $G^{ab}$ with respect to
the basis $\hat{\xi}^{a}\equiv(\hat{q},\hat{p})^{\intercal}$ are given by

$\displaystyle\Omega^{ab}\equiv\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}\quad\text{and}\quad G^{ab}\equiv\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\,.$ (193)

Just as in the bosonic case, we can construct the linear complex structure

$\displaystyle J^{a}{}_{b}=\Omega^{ac}G^{-1}_{cb}\equiv\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}\,.$ (194)

With this, the fermionic projector becomes

$\displaystyle\frac{1}{2}(\delta^{a}{}_{b}+\mathrm{i}J^{a}{}_{b})\hat{\xi}^{b}\equiv\frac{1}{2}\begin{pmatrix}\hat{q}+\mathrm{i}\hat{p}\\\
\hat{p}-\mathrm{i}\hat{q}\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}\hat{a}\\\
-\mathrm{i}\hat{a}\end{pmatrix}\,,$ (195)

_i.e._ , the condition
$\frac{1}{2}(\delta^{a}{}_{b}+\mathrm{i}J^{a}{}_{b})\hat{\xi}^{b}\ket{\Omega}=0$
is the same as for bosons. The projector selects the subspace spanned by
complex linear combinations of annihilation operators
$\hat{a}=\frac{1}{\sqrt{2}}(\hat{q}+\mathrm{i}\hat{p})$.

###  4.3 Quadratic fermionic Hamiltonians

We are considering the most general quadratic fermionic Hamiltonian given by

$\displaystyle\hat{H}=\frac{\mathrm{i}}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}\,,$
(196)

where $h_{ab}$ is antisymmetric and non-degenerate. Note that $\hat{H}$ is
always bounded from below.

The ground state of $\hat{H}$ is given by a Gaussian state $\ket{\Omega}$. As
for bosons, we can compute the ground state covariance matrix from $h_{ab}$:

  * •

Computation of $\Omega^{ab}$
Similar to (86), we compute

$\displaystyle\Omega^{ab}=K^{a}{}_{c}|K^{-1}|^{c}{}_{d}G^{db}$ (197)

with $K^{a}{}_{b}=G^{ac}h_{cb}$.

The energy expectation value
$E=\bra{\Omega^{\prime}}\hat{H}\ket{\Omega^{\prime}}$ and its variance
$\Sigma_{E}$ for an arbitrary Gaussian state $\ket{\Omega^{\prime}}$ can be
efficiently computed as

$\displaystyle E$
$\displaystyle=-\frac{1}{4}h_{ab}{\Omega^{\prime}}^{ab}=-\frac{\operatorname{Tr}(h{\Omega^{\prime}}^{\intercal})}{4}\,,$
(198) $\displaystyle\Sigma_{E}$
$\displaystyle=\frac{\sqrt{-\operatorname{Tr}(hGhG+h{\Omega^{\prime}}h{\Omega^{\prime}})}}{4}\,.$
(199)

Again, we note that the symplectic form $\Omega^{ab}$ and the positive
definite metric $G^{ab}$ switched their roles with respect to the bosonic
formulae in Section 2.

##### Example.

We consider the harmonic oscillator of a single fermionic mode
$\hat{H}=\frac{\mathrm{i}}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}=\mathrm{i}\hat{q}\hat{p}$.
The bilinear form $h$ with respect to $\hat{\xi}^{a}\equiv(\hat{q},\hat{p})$
is

$\displaystyle h_{ab}\equiv\begin{pmatrix}0&\omega\\\
-\omega&0\end{pmatrix}\,,$ (200)

which gives via (197) rise to

$\displaystyle K^{a}{}_{b}\equiv\begin{pmatrix}0&\omega\\\
-\omega&0\end{pmatrix}\,\,\Rightarrow\,\,\Omega^{ab}\equiv\begin{pmatrix}0&1\\\
-1&0\end{pmatrix}\,,$ (201)

where we used $|K^{-1}|=\mathbb{1}/\omega$, just like in the bosonic case. The
ground state energy is given by $E=\bra{\Omega}\hat{H}\ket{\Omega}=-\omega/2$
and expressing the Hamiltonian in creation and annihilation operators gives
$\hat{H}=\frac{\omega}{2}(\hat{a}^{\dagger}\hat{a}-1)$.

###  4.4 Entanglement of fermionic Gaussian states

Given a fermionic Gaussian state $\ket{\Omega}$ and a system decomposition
into two even dimensional orthogonal complements $V=A\oplus B$, we can always
choose an orthonormal basis
$(q_{1}^{A},p_{1}^{A},\dots,p_{N_{A}}^{A},q^{A}_{N_{A}},q_{1}^{B},p_{1}^{B},\dots,p_{N_{B}}^{B},q^{B}_{N_{B}})$,
such that the covariance matrix takes the standard form derived in appendix
A.2

$\displaystyle\Omega\equiv\left(\begin{array}[]{ccc|cccccc}\mathbf{cos}_{1}&\cdots&0&\mathbf{sin}_{1}&\cdots&0&0&\cdots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\cdots&\mathbf{cos}_{N_{A}}&0&\cdots&\mathbf{sin}_{N_{A}}&0&\cdots&0\\\
\hline\cr-\mathbf{sin}_{1}&\cdots&0&\mathbf{cos}_{1}&\cdots&0&0&\cdots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\cdots&-\mathbf{sin}_{N_{A}}&0&\cdots&\mathbf{cos}_{N_{A}}&0&\cdots&0\\\
0&\cdots&0&0&\cdots&0&\mathbb{A}_{2}&\cdots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\cdots&0&0&\cdots&0&0&\cdots&\mathbb{A}_{2}\end{array}\right)\,,$ (211)

where $\mathbb{A}_{2}$, $\mathbf{cos}_{i}$ and $\mathbf{sin}_{i}$ are given by

$\displaystyle\mathbb{A}_{2}$ $\displaystyle=\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right)\,,$ (214) $\displaystyle\mathbf{cos}_{i}$
$\displaystyle=\left(\begin{array}[]{cc}0&\cos{2r_{i}}\\\
-\cos{2r_{i}}&0\end{array}\right)\,,$ (217) $\displaystyle\mathbf{sin}_{i}$
$\displaystyle=\left(\begin{array}[]{cc}0&\sin{2r_{i}}\\\
\sin{2r_{i}}&0\end{array}\right)\,.$ (220)

This decomposition is analogous to the bosonic one in (51), where the mode
$(q_{i}^{A},p_{i}^{A})$ is entangled with the mode $(q_{i}^{B},p_{i}^{B})$.
The parameters $r_{i}$ can always be chosen to lie in the interval
$[0,\pi/4]$. Again, the amount of entanglement is encoded in the parameter
$r_{i}$, but now given by [73]

$\displaystyle S_{A}(\ket{\Omega})=\sum^{N_{A}}_{i=1}s_{f}(\cos{2r_{i}})$
(221)

with
$s_{f}(x)=-\frac{1+x}{2}\log\left(\frac{1+x}{2}\right)-\frac{1-x}{2}\log\left(\frac{1-x}{2}\right)$.
As for bosons, there exists a compact formula to express the entanglement
entropy in terms of the linear complex structure $J$, namely

$\displaystyle\textstyle
S_{A}(\ket{G,z})=-\mathrm{Tr}\left[\left(\frac{\mathbb{1}_{A}+\mathrm{i}[J]_{A}}{2}\right)\log\left(\frac{\mathbb{1}_{A}+\mathrm{i}[J]_{A}}{2}\right)\,\\!\right]\,,$
(222)

which closely resembles (59), but where the modulus is replaced by an overall
minus sign.

###  4.5 Fermionic partner mode construction

The standard form of the fermionic covariance matrix introduced above implies
that also for a system of fermionic modes in an overall pure state, each
single mode has a unique partner with which it is entangled.

In fact the partner mode for bosonic modes can be almost directly carried over
to the fermionic case. Only the hypergeometric functions need to be replaced
by their trigonometric analogues. If $(x_{a},k_{a})$ define a mode in standard
form

$\displaystyle G^{ab}x_{a}x_{b}=G^{ab}k_{a}k_{b}=1\,,$ (223) $\displaystyle
G^{ab}x_{a}k_{b}=0\,,$ (224) $\displaystyle\Omega^{ab}x_{a}k_{b}=\cos 2r\,,$
(225)

then its partner mode is given by

$\displaystyle\bar{x}_{a}$ $\displaystyle=\cot(2r)x_{a}+\frac{1}{\sin
2r}\left(J^{\intercal}\right)_{a}^{c}k_{c}\,,$ (226)
$\displaystyle\bar{k}_{a}$ $\displaystyle=-\cot(2r)k_{a}+\frac{1}{\sin
2r}\left(J^{\intercal}\right)_{a}^{c}x_{c}\,.$ (227)

This means that the partner mode is in standard form itself

$\displaystyle G^{ab}\bar{x}_{a}\bar{x}_{b}=G^{ab}\bar{k}_{a}\bar{k}_{b}=1\,,$
(228) $\displaystyle G^{ab}\bar{x}_{a}\bar{k}_{b}=0\,,$ (229)
$\displaystyle\Omega^{ab}\bar{x}_{a}\bar{k}_{k}=\cos 2r\,,$ (230)

and the correlations with the original mode are as required by the standard
form

$\displaystyle
G^{ab}x_{a}\bar{x}_{b}=G^{ab}k_{a}\bar{k}_{b}=G^{ab}x_{a}\bar{k}_{b}=0\,$
(231)
$\displaystyle\Omega^{ab}x_{a}\bar{k}_{b}=\Omega^{ab}k_{a}\bar{x}_{b}=\sin
2r\,.$ (232)

Just as in the bosonic case, also for fermionic modes the partner of the
partner mode is the original mode itself. Also, we observe that any other mode
$(x^{\prime},k^{\prime})$ which anti-commutes with $(x,k)$, has an anti-
commutatator with the partner mode proportional to the covariance with the
first mode, _i.e._ ,

$\displaystyle G^{ab}x^{\prime}_{a}\bar{x}_{b}=\cot
2r\,G^{ab}x^{\prime}_{a}x_{b}-\frac{1}{\sin
2r}\Omega^{ab}x^{\prime}_{a}k_{b}\,.$ (233)

This means that any mode which anticommutes with both partner modes is not
correlated with them. In particular, a basis of modes which contains the two
partner modes brings the covariance matrix of the state into block diagonal
form

$\displaystyle\Omega\equiv\left(\begin{array}[]{c|c}\mathbf{\tilde{\Omega}}&0\\\
\hline\cr 0&\mathbf{\Omega}_{R}\end{array}\right)$ (236)

where

$\displaystyle\mathbf{\tilde{\Omega}}=\left(\begin{array}[]{cccc}0&\cos
2r&0&\sin 2r\\\ -\cos 2r&0&\sin 2r\\\ 0&-\sin 2r&0&\cos 2r\\\ -\sin 2r&0&-\cos
2r\end{array}\right)$ (241)

takes the standard form of (211).

Also fermionic partner modes can be unsqueezed to find a pair of modes $(y,l)$
and $(z,m)$ which are in a pure product state and from which $(x,k)$ and
$(\bar{x},\bar{k})$ are obtained by squeezing. These modes are given by:

$\displaystyle y_{a}$
$\displaystyle=\cos(r)x_{a}-\sin(r)\bar{x}_{a}=\frac{x_{a}-(J^{\intercal})_{a}{}^{b}k_{b}}{2\cos
r}\,,$ (242) $\displaystyle l_{a}$
$\displaystyle=\cos(r)k_{a}+\sin(r)\bar{k}_{a}=\frac{k_{a}+(J^{\intercal})_{a}{}^{b}x_{b}}{2\cos
r}\,,$ (243) $\displaystyle z_{a}$
$\displaystyle=\sin(r)x_{a}+\cos(r)\bar{x}_{a}=\frac{x_{a}+(J^{\intercal})_{a}{}^{b}k_{b}}{2\sin
r}\,,$ (244) $\displaystyle m_{a}$
$\displaystyle=-\sin(r)k_{a}+\cos(r)\bar{k}_{a}=\frac{-k_{a}+(J^{\intercal})_{a}{}^{b}x_{b}}{2\sin
r}\,.$ (245)

##  5 Entanglement extraction from quadratic fermionic systems

In this section, we analyze the minimal energy cost of pure state entanglement
extraction from the ground state of a system of fermionic modes with a
quadratic Hamiltonian. Our approach is motivated by the same framework as
discussed for bosonic modes in Section 3. Again, after reviewing the
extraction procedure in section 5.1, our proof follows exactly the same three
steps, split over the sections 5.2, 5.3 and 5.4.

###  5.1 Pure state entanglement extraction from bosonic modes

We assume that the source system consists of fermionic modes which are in the
ground state $\ket{\Psi}_{S}$ of an arbitrary quadratic Hamiltonian $\hat{H}$,
which is Gaussian. Within the source system, we select two anti-commuting
source modes from which we want to extract entanglement. For this we assume
that it is possible to swap the state of the two source modes onto a pair of
target modes outside of the system. We can always use Gaussian evolution,
_i.e._ , quadratic coupling terms, to swap the entangled Gaussian state in the
source system with the unentangled Gaussian state in the target system. The
target modes can be prepared in an arbitrary product state. Hence the joint
initial state of target modes and source system is of the form

$\displaystyle\rho_{i}=\sigma_{1}\otimes\tau_{2}\otimes\ket{\Psi}\bra{\Psi}_{S}.$
(246)

Just as in the bosonic case, also for fermionic modes the partner mode
structure of entanglement of Gaussian states implies that the only way to
extract a pure entangled state in this way is to chose the two source modes to
be partner modes which we denote by $A$ and $\bar{A}$. This is also the
optimal choice in the sense that any entanglement measure of one mode with
with a second non-partner mode is never larger than of the mode with its
partner.

The ground state of the source system is a product state
$\ket{\Psi}=\ket{\psi}_{A\bar{A}}\otimes\ket{\phi}_{R}$ between the partner
modes and the remainder $R$ of the system. After the entanglement extraction,
_i.e._ , after the states of the partner modes and the target modes have been
swapped, the source system is in the state

$\displaystyle\rho^{\prime}_{S}=\sigma_{A}\otimes\tau_{\bar{A}}\otimes\ket{\phi}\bra{\phi}_{R}$
(247)

where now the partner modes are in the product state given by the initial
state of the target modes.

###  5.2 Energy cost of entanglement extraction

By the minimal energy cost of entanglement extraction we are referring to the
minimum value by which the expectation value of $\hat{H}$ is increased, _i.e._
, by how much energy on average is injected into the system when the modes are
swapped. Since we have used the same notation, the discussion and formulae
from Section 3.2 apply equally to the fermionic case, as considered here. In
particular, the energy cost for the extraction of a fermionic pair of partner
modes is given by

$\displaystyle\Delta
E=\operatorname{Tr}\left(\sigma_{A}\otimes\tau_{\bar{A}}\hat{H}_{A\bar{A}}-\ket{\psi}\bra{\psi}_{A\bar{A}}\hat{H}_{A\bar{A}}\right)\,.$
(248)

It results from the restriction $\hat{H}_{A\bar{A}}$ of the full Hamiltonian
onto the space spanned by the two partner modes. Thus, again the problem of
finding the minimal energy cost reduces to a two-mode problem.

Also as before, we find that the final energy expectation value of the system
is

$\displaystyle\operatorname{Tr}\left(\sigma_{A}\otimes\tau_{\bar{A}}\hat{H}_{A\bar{A}}\right)=\operatorname{Tr}\left(\sigma_{A}\hat{H}_{A}\right)+\operatorname{Tr}\left(\sigma_{\bar{A}}\hat{H}_{\bar{A}}\right)$
(249)

which implies that in order to minimize the energy cost, the target modes need
to be initialized in the ground states of the single-mode restrictions
$\hat{H}_{A}$ and $\hat{H}_{\bar{A}}$ of $\hat{H}$ onto the individual partner
modes. Also for fermionic modes these states are Gaussian states, as discussed
in the previous section. Thus, we can perform our entire analysis within the
framework of Gaussian states.

###  5.3 Energy cost for two fermionic modes

This section presents the fermionic analogue to 3.3. We analyze entanglement
extraction from a source system that consists of only two fermionic modes and
has a quadratic Hamiltonian.

Since the Hamiltonian is quadratic, we can choose a basis
$\hat{\xi}^{a}\equiv(\hat{q}_{1},\hat{p}_{1},\hat{q}_{2},\hat{p}_{2})$ of
energy eigenmodes such that the quadratic part of the Hamiltonian $h_{ab}$ and
the covariance matrix $\Omega^{ab}$ of its ground state are both represented
by block diagonal matrices $\mathbf{h}$ and $\mathbf{\Omega}$ given by

$\displaystyle h_{ab}$ $\displaystyle\stackrel{{\scriptstyle
1,2}}{{\equiv}}\mathbf{h}=\left(\begin{array}[]{cccc}0&\epsilon_{1}&0&0\\\
-\epsilon_{1}&0&0&0\\\ 0&0&0&\epsilon_{2}\\\
0&0&-\epsilon_{2}&0\end{array}\right)\,,$ (254) $\displaystyle\Omega^{ab}$
$\displaystyle\stackrel{{\scriptstyle
1,2}}{{\equiv}}\mathbf{\Omega}=\left(\begin{array}[]{cccc}0&1&0&0\\\
-1&0&0&0\\\ 0&0&0&1\\\ 0&0&-1&0\end{array}\right)\,.$ (259)

The energy eigenmodes are not entangled and, therefore, cannot be used for
entanglement extraction. Instead, we need to find a pair of partner modes
$(\hat{Q}_{A},\hat{P}_{A},\hat{Q}_{\bar{A}},\hat{P}_{\bar{A}})$ with respect
to which $\Omega^{ab}$ is represented by the matrix corresponding to the
standard form (241)

$\displaystyle\Omega^{ab}\stackrel{{\scriptstyle
A,\bar{A}}}{{\equiv}}\mathbf{\tilde{\Omega}}$ (260)

for some fixed $r\in[0,\pi/4]$, which depends on the amount of entanglement we
want to extract. Among all choices of partner modes $A$ and $\bar{A}$, we want
to identify those that minimize the energy cost associated to the entanglement
extraction. Consequently, we will search for the orthogonal transformation
which maps the energy eigenmodes $1$ and $2$ to those optimal partner modes
$A$ and $\bar{A}$ for entanglement extraction.

The two-mode squeezing transformation

$\displaystyle\mathbf{M}_{r}=\left(\begin{array}[]{cccc}\cos{r}&0&\sin{r}&0\\\
0&\cos{r}&0&-\sin{r}\\\ -\sin{r}&0&\cos{r}&0\\\
0&\sin{r}&0&\cos{r}\end{array}\right)$ (265)

brings $\mathbf{\Omega}$ into the form from (241) via

$\displaystyle\mathbf{\tilde{\Omega}}=\mathbf{M}_{r}\mathbf{\Omega}\mathbf{M}^{\intercal}_{r}\,,$
(266)

_i.e._ , it maps the eigenmodes to a pair of partner modes with squeezing
parameter $r$. Under the orthogonal transformation $\mathbf{M}_{r}$, the
Hamiltonian matrix $\mathbf{h}$ transforms under the action of the inverse
transform $\mathbf{M}^{-1}_{r}=\mathbf{M}_{-r}$, _i.e._ ,

$\displaystyle\mathbf{\tilde{h}}=\mathbf{M}_{-r}^{\intercal}\mathbf{h}\mathbf{M}_{-r}\,.$
(267)

However, $\mathbf{M}_{r}$ is not the only transformation doing this: We can
compose the squeezing $\mathbf{M}_{r}$ with any orthogonal transformation
$\mathbf{N}$ which leaves the covariance matrix invariant, _i.e._ ,
$\mathbf{\Omega}=\mathbf{N}\mathbf{\Omega}\mathbf{N}^{\intercal}$, to obtain
the most general transformation $\mathbf{T}_{r}=\mathbf{M}_{r}\mathbf{N}$ that
transforms the modes $1$ and $2$ into $A$ and $\bar{A}$ for fixed $r$. All of
these choices share the same amount of entanglement across $A$ and $\bar{A}$
because the covariance matrix $\Omega$ takes the same form
$\mathbf{\tilde{\Omega}}=\mathbf{T}_{r}\mathbf{\Omega}\mathbf{T}_{r}^{\intercal}$
with respect to them. However, the Hamiltonian matrix $\mathbf{h}$ is in
general not left invariant by $\mathbf{N}$, _i.e._ ,
$\mathbf{h}\neq(\mathbf{N}^{-1})^{\intercal}\mathbf{h}\mathbf{N}^{-1}$ and
thus the energy cost can be different.

####  5.3.1 Degenerate two-mode Hamiltonian
$(\epsilon_{1}\\!=\\!\epsilon_{2})$

We first consider the case where $\epsilon:=\epsilon_{1}=\epsilon_{2}$. This
case is the simplest to solve because the matrix representations $\mathbf{h}$
and $\mathbf{\Omega}$ are proportional to each other in the eigenmode basis
and, thus, in all bases. Indeed, any orthogonal transformation $\mathbf{N}$
that leaves $\mathbf{\Omega}$ invariant, _i.e._ ,
$\mathbf{N}\mathbf{\Omega}\mathbf{M}^{\intercal}=\mathbf{\Omega}$, also leaves
$\mathbf{h}$ invariant, _i.e._ ,
$\mathbf{h}=(\mathbf{N}^{-1})^{\intercal}\mathbf{h}\mathbf{N}^{-1}$.
Therefore, for any pair of partner modes $A$ and $\bar{A}$ where the
covariance matrix takes the standard form $\mathbf{\tilde{\Omega}}$ for fixed
$r$, we find

$\displaystyle\mathbf{\tilde{h}}=\left(\begin{array}[]{cc|cc}0&\epsilon\cos{2r}&0&-\epsilon\sin{2r}\\\
-\epsilon\cos{2r}&0&-\epsilon\sin{2r}&0\\\ \hline\cr
0&\epsilon\sin{2r}&0&\epsilon\cos{2r}\\\
\epsilon\sin{2r}&0&-\epsilon\cos{2r}&0\end{array}\right).$ (272)

As discussed before, the entanglement extraction swaps the state of the chosen
partner modes against a product state $\sigma_{\bar{A}}\otimes\tau_{\bar{A}}$.
This brings the covariance matrix into block diagonal form with respect to the
partner modes:

$\displaystyle\mathbf{\tilde{\Omega}}^{\prime}=\left(\begin{array}[]{c|c}\mathbf{\tilde{\Omega}}_{A}&0\\\\[2.84526pt]
\hline\cr\\\\[-8.53581pt]
0&\mathbf{\tilde{\Omega}}_{\bar{A}}\end{array}\right).$ (276)

In terms of this new covariance matrix, the energy cost (248) of the
entanglement extraction is

$\displaystyle\textstyle\Delta E$
$\displaystyle=-\frac{1}{4}\operatorname{Tr}\left(\mathbf{\tilde{h}}\,\mathbf{\tilde{\Omega}}^{\intercal}\right)+\frac{1}{4}\operatorname{Tr}\left(\mathbf{h}\,{\mathbf{\Omega}}^{\intercal}\right)$
(277)
$\displaystyle=\epsilon-\frac{1}{4}\operatorname{Tr}\left(\mathbf{\tilde{h}}_{A}\mathbf{\tilde{\Omega}}_{A}+\mathbf{\tilde{h}}_{\bar{A}}\mathbf{\tilde{\Omega}}_{\bar{A}}\right)\,,$
(278)

where
$\mathbf{\tilde{h}_{A}}=\mathbf{\tilde{h}_{\bar{A}}}=\epsilon\cos{2r}\,\mathbb{A}_{2}$
are the Hamiltonian matrices of the restriction of $\hat{H}$ onto the single
partner modes respectively, _i.e._ , the diagonal blocks of
$\mathbf{\tilde{h}}$ in (272). To minimize $\Delta E$, the target modes need
to be prepared in the ground state of $\mathbf{\tilde{h}_{A}}$ and
$\mathbf{\tilde{h}_{\bar{A}}}$, _i.e._ ,
$\mathbf{\tilde{\Omega}_{A}}^{\prime}=\mathbf{\tilde{\Omega}_{\bar{A}}}^{\prime}=\mathbb{A}_{2}$.

With these initial states, the minimal energy cost for the extraction of a
pure two-mode squeezed state with entanglement entropy $\Delta S=S_{f}(\cos
2r)$ is

$\displaystyle\Delta E$ $\displaystyle=2\epsilon\sin^{2}{r}\,.$ (279)

In addition to this elevation of the energy expectation value, the final state
of the source system has an energy variance of

$\displaystyle\textstyle\Sigma_{E}$
$\displaystyle=\frac{\sqrt{-\operatorname{Tr}\left(\mathbf{\tilde{h}}\mathbf{\tilde{G}}\mathbf{\tilde{h}}\mathbf{\tilde{G}}+\mathbf{\tilde{h}}\mathbf{\tilde{\Omega}}^{\prime}\mathbf{\tilde{h}}\mathbf{\tilde{\Omega}}^{\prime}\right)}}{4}$
(280) $\displaystyle=\frac{\epsilon}{\sqrt{2}}\sin{2r}\,.$ (281)

####  5.3.2 General two-mode Hamiltonian
($\epsilon_{1}\\!\leq\\!\epsilon_{2}$)

The case of a general Hamiltonian, with $\epsilon_{1}\leq\epsilon_{2}$, is
slightly more involved than the degenerate case. Interestingly, there exists a
critical value of extracted entanglement $\Delta S$ below which the minimal
energy cost is independent of the choice of partner modes. Above the critical
amount, it is optimal to chose partner modes obtained from squeezing the
eigenmodes of $\hat{H}$ directly.

In the non-degenerate case the _orthogonal_ transformations which map the
energy eigenmodes to any pair of partner modes is of the general form
$\mathbf{T}_{r}=\mathbf{M}_{r}\mathbf{N}$. This means that before the
squeezing transformation $\mathbf{M}_{r}$ we can apply a transformation
$\mathbf{N}$ which can change $\mathbf{h}$ but leaves $\mathbf{\Omega}$
invariant. To find all possible forms of $\mathbf{N}$ we write it as a matrix
exponential $\mathbf{N}=\exp\left(\theta\mathbf{K}\right)$ of the generator
$\mathbf{K}$ subject to the constraints
$\mathbf{K}\mathbf{G}+\mathbf{G}\mathbf{K}^{\intercal}=0$ (leaving
$\mathbf{G}$ invariant, _i.e._ , being orthogonal) and
$\mathbf{K}\mathbf{\Omega}+\mathbf{\Omega}\mathbf{K^{\intercal}}=0$ (leaving
$\mathbf{\Omega}$ invariant, _i.e._ , being symplectic).222These are indeed
the same conditions as for bosons. Mathematically speaking, the intersection
of the orthogonal algebra $\mathrm{so}(2N)$ and the symplectic algebra
$\mathrm{sp}(2N)$ is given by the unitary algebra $\mathrm{u}(N)$. From this,
we find the most general form

$\displaystyle\mathbf{K}$
$\displaystyle=\left(\begin{array}[]{cccc}0&0&\cos{\phi}&-\sin{\phi}\\\
0&0&\sin{\phi}&\cos{\phi}\\\ -\cos{\phi}&-\sin{\phi}&0&0\\\
\sin{\phi}&-\cos{\phi}&0&0\end{array}\right)\,,$ (286)

giving rise to $\mathbf{N}=\exp(\theta\mathbf{K})$.

Consequently, $\mathbf{T}_{r}=\mathbf{M}_{r}\mathbf{N}$ provides the most
general transformation which brings the covariance matrix
$\mathbf{\tilde{\Omega}}=\mathbf{T}_{r}\mathbf{\Omega}\mathbf{T}_{r}^{\intercal}$
in the form of (241). Such a transformation brings the Hamiltonian matrix into
the general form

$\displaystyle\mathbf{\tilde{h}}=(\mathbf{T}_{r}^{-1})^{\intercal}\mathbf{h}\mathbf{T}_{r}^{-1}=\left(\begin{array}[]{c|c}\mathbf{h}_{A}&\mathbf{h}_{X}\\\\[2.84526pt]
\hline\cr\\\\[-11.38109pt]
-\mathbf{h}_{X}^{\intercal}&\mathbf{h}_{\bar{A}}\end{array}\right)\,,$ (290)

with blocks that read, using
$\epsilon_{\pm}=\frac{1}{2}(\epsilon_{2}\pm\epsilon_{1})$,

$\displaystyle\mathbf{h}_{A}$
$\displaystyle=\begin{pmatrix}0&\epsilon_{+}\cos{2r}-\epsilon_{-}\cos{2\theta}\\\
-\epsilon_{+}\cos{2r}+\epsilon_{-}\cos{2\theta}&0\end{pmatrix}\,,$ (291)
$\displaystyle\mathbf{h}_{\bar{A}}$
$\displaystyle=\begin{pmatrix}0&\epsilon_{+}\cos{2r}+\epsilon_{-}\cos{2\theta}\\\
-\epsilon_{+}\cos{2r}-\epsilon_{-}\cos{2\theta}&0\end{pmatrix}\,,$ (292)
$\displaystyle\mathbf{h}_{X}$ $\displaystyle=\begin{pmatrix}\epsilon_{-}\sin
2\theta\sin\phi&\epsilon_{+}\sin 2r+\epsilon_{-}\cos\phi\sin 2\theta\\\
\epsilon_{+}\sin 2r-\epsilon_{-}\cos\phi\sin 2\theta&\epsilon_{-}\sin
2\theta\sin\phi\end{pmatrix}\,.$ (293)

The excitation energies of the restricted one-mode Hamiltonian $\hat{H}_{A}$
and $\hat{H}_{\bar{A}}$ are thus

$\displaystyle\epsilon_{A}=|\epsilon_{+}\cos{2r}-\epsilon_{-}\cos{2\theta}|\,,$
(294)
$\displaystyle\epsilon_{\bar{A}}=|\epsilon_{+}\cos{2r}+\epsilon_{-}\cos{2\theta}|\,.$
(295)

As discussed above, the entanglement extraction brings the partner modes $A$
and $\bar{A}$ into a pure Gaussian product state. Since $A$ and $\bar{A}$ are
fermionic modes, there are only four different product states. Their
covariance matrices are

$\displaystyle\mathbb{A}^{(k,l)}=\left(\begin{array}[]{c|c}(-1)^{k}\mathbb{A}_{2}&0\\\\[2.84526pt]
\hline\cr\\\\[-8.53581pt] 0&(-1)^{l}\mathbb{A}_{2}\end{array}\right)\,,$ (299)

with $k,l\in\\{0,1\\}$. Which of these states minimizes the energy cost,
_i.e._ , is the product of the ground states of $\mathbf{h}_{A}$ and
$\mathbf{h}_{\bar{A}}$, depends on spectrum of the Hamiltonian and the
squeezing parameter $r$. This is because the energy cost arising from the
different $\mathbb{A}^{(k,l)}$ is

$\displaystyle\Delta
E=\epsilon_{+}-\frac{1}{4}{\operatorname{Tr}\left(\mathbf{\tilde{h}}\,\mathbb{A}^{(k,l)}\right)}$
$\displaystyle=\left\\{\begin{array}[]{cl}\epsilon_{+}(1-\cos 2r),&\text{if
}(k,l)=(0,0)\\\ \epsilon_{+}-\epsilon_{-}\cos 2\theta,&\text{if
}(k,l)=(1,0)\\\ \epsilon_{+}+\epsilon_{-}\cos 2\theta,&\text{if
}(k,l)=(0,1)\\\ \epsilon_{+}(1+\cos 2r),&\text{if
}(k,l)=(1,1)\end{array}\right.\,.$ (304)

The fact that for two of the states the energy cost is independent of $\theta$
versus it being independent of $r$ for the other two states, has an
interesting consequence: It means that for low values of $r$, _i.e._ , for low
amounts of entanglement, the minimal energy cost is independent of $\theta$
and given by $\Delta E_{\text{min}}=2\epsilon_{+}\sin^{2}r$ which is achieved
by the state $\mathbb{A}^{(0,0)}$. However, as soon as $\cos
2r\leq\frac{\epsilon_{-}}{\epsilon_{+}}$ is stisfied, it is more beneficial to
bring the modes into the final state $\mathbb{A}^{(1,0)}$ which results in an
energy cost $\Delta E=\epsilon_{+}-\epsilon_{-}\cos 2\theta$. In this case,
the energy cost is minimized by the parameter value $\theta=n\pi$ with
$n\in\mathbb{Z}$, _i.e._ , if the partner modes are obtained by squeezing
directly the energy eigenmodes.

In summary, the minimal energy cost for the extraction of a pair of partner
modes sharing entanglement entropy $\Delta S=S_{f}\left(\cos 2r\right)$ from a
fermionic quadratic Hamiltonian with excitation energies
$\epsilon_{1}\leq\epsilon_{2}$, shown in Figure 4, is

$\displaystyle\Delta
E_{\mathrm{min}}=\begin{cases}(\epsilon_{1}+\epsilon_{2})\sin^{2}{r},&\text{if
}\cos{2r}\geq\frac{\epsilon_{2}-\epsilon_{1}}{\epsilon_{2}+\epsilon_{1}}\\\
\epsilon_{1},&\text{if
}\cos{2r}\leq\frac{\epsilon_{2}-\epsilon_{1}}{\epsilon_{2}+\epsilon_{1}}\end{cases}\,.$
(305)

In particular, the minimal energy cost is bounded from above by the lower
excitation energy of the Hamiltonian $\epsilon_{1}$. This is because above the
critical squeezing parameter with $\cos 2r\leq\epsilon_{-}/\epsilon_{+}$, the
product state $\mathbb{A}^{(1,0)}$ in which the partner modes are left after
the extraction is exactly the first excited state of the Hamiltonian
$\hat{H}$. In fact, the covariance matrix $\mathbb{A}^{(0,1)}$ is left
invariant by the squeezing operation
$\mathbf{M}_{r}\mathbb{A}^{(1,0)}\mathbf{M}_{r}^{\intercal}$, _i.e._ , the
first excited state of $\hat{H}$ has the same covariance matrix with respect
to the energy eigenmodes and with respect to the partner modes $A$ and
$\bar{A}$ which are obtained by squeezing the energy eigenmodes.

Figure 4: Minimal energy cost $\Delta E_{\mathrm{min}}$ of extracting
entanglement entropy $\Delta S$ from a fermionic, quadratic two-mode
Hamiltonian with excitation energies $\epsilon_{1}$ and $\epsilon_{2}$.

As a consequence, in the regime $\cos{2r}\leq\epsilon_{-}/\epsilon_{+}$, the
energy variance $\Sigma_{E}$ of the final system state vanishes because the
system is in an energy eigenstate. For lower amounts of entanglement, where
this is not the case, the variance of the state which minimizes the energy
cost depends on the squeezing parameter. It is plotted in Figure 5 and reads

$\displaystyle\Sigma_{E}$
$\displaystyle=\left\\{\begin{array}[]{lcl}\frac{\epsilon_{+}}{\sqrt{2}}\sin{2r},&&\text{if
}\cos{2r}>\frac{\epsilon_{2}-\epsilon_{1}}{\epsilon_{2}+\epsilon_{1}}\\\
0,&&\text{if
}\cos{2r}\leq\frac{\epsilon_{2}-\epsilon_{1}}{\epsilon_{2}+\epsilon_{1}}\end{array}\right.\,,$
(308)

Note that for
$\cos{2r}>\frac{\epsilon_{2}-\epsilon_{1}}{\epsilon_{2}+\epsilon_{1}}$, it is
possible to leave the source system in an energy eigenstate such that
$\Sigma_{E}=0$. The energy cost for this is $\epsilon_{1}$.

Another relevant figure of merit is to look at the ratio between a given
amount of extracted entanglement and its minimal energy cost $\Delta S/\Delta
E$. Since $\Delta S=S_{f}(\cos{2r})\sim(1-2\log{r})\,r^{2}$ as $r\to 0$, the
ratio diverges as

$\displaystyle\lim_{r\to 0}\frac{\Delta S}{\Delta E_{\text{min}}}=\infty$
(309)

for small amounts of extracted entanglement, just as we found for bosons in
(152). As the extracted entanglement increases the ratio is strictly
decreasing and reaches its minimum $\log(2)/\epsilon_{1}$ for $r=\pi/4$.
Similar to the bosonic case, we thus find that the price per entanglement bit
increases as we are trying extract a larger total amount.

Figure 5: Energy variance $\Sigma_{E}$ of the state minimizing the energy cost
of extracting partner modes with entanglement entropy $\Delta S$ from a
fermionic, quadratic two-mode Hamiltonian with excitation energies
$\epsilon_{1}$ and $\epsilon_{2}$.

###  5.4 Choosing optimal partner modes

The analysis of the previous sections explain how the increase of energy
expectation value of the source system can be minimized when extracting a pair
of partner modes sharing a given amount of entanglement: To keep the energy
cost as low as possible the modes need to be constructed from the two lowest
energy eigenmodes of the Hamiltonian. Specifically, the partner modes should
be obtained by directly squeezing the two energy eigenmodes, such that always
the lowest possible energy cost according to (305) is achieved. Again, as
discussed in Appendix B, no other choice of partner modes yields a lower
energy cost because the energy spectrum of the restricted Hamiltonian is
bounded by the spectrum of the full Hamiltonian.

For a general pair of partner modes $A$ and $\bar{A}$ the minimal energy cost
is determined by the spectrum of the restricted two-mode Hamiltonian
$\hat{H}_{A\bar{A}}$, which enters the formula (305). In contrast to the
fermionic case, an arbitrarily chosen pair of partner modes will in fact
achieve this minimal energy price if the amount of extracted entanglement is
below the critical value, _i.e._ , if
$\cos{2r}\geq\frac{\epsilon_{-}}{\epsilon_{+}}$. If the squeezing is above
this value, according to the analysis of the previous section only the partner
modes obtained by directly squeezing the energy eigenmodes of
$\hat{H}_{A\bar{A}}$ do achieve the possible minimal energy cost of $\Delta
E_{\text{min}}=\epsilon_{1}$.

To perform this analysis for a given choice of partner modes, it is necessary
to express the Hamiltonian matrix of $\hat{H}_{A\bar{A}}$ with respect to the
partner mode operators
$\left(\hat{Q}_{A},\hat{P}_{A},\hat{Q}_{\bar{A}},\hat{P}_{\bar{A}}\right)$. As
before, we denote these mode operators as

$\displaystyle\hat{Q}_{A}$
$\displaystyle=x_{a}\hat{\xi}^{a}\,,\qquad\hat{P}_{A}=k_{a}\hat{\xi}^{a}\,,$
(310) $\displaystyle\hat{Q}_{\bar{A}}$
$\displaystyle=\bar{x}_{a}\hat{\xi}^{a},\qquad\hat{P}_{\bar{A}}=\bar{k}_{a}\hat{\xi}^{a}\,.$
(311)

In this basis, the Hamiltonian matrix representing the restriction
$\hat{H}_{A\bar{A}}$ of a general quadratic Hamiltonian
$\hat{H}=\frac{\mathrm{i}}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}$ is given by

$\displaystyle\left(\begin{array}[]{cc|cc}0&\check{h}(x,k)&\check{h}(x,\bar{x})&\check{h}(x,\bar{k})\\\
\check{h}(k,x)&0&\check{h}(k,\bar{x})&\check{h}(k,\bar{k})\\\
\hline\cr\check{h}(\bar{x},x)&\check{h}(\bar{x},k)&0&\check{h}(\bar{x},\bar{k})\\\
\check{h}(\bar{k},x)&\check{h}(\bar{k},k)&\check{h}(\bar{k},\bar{x})&0\end{array}\right)$
(316)

where we now define $\check{h}:V^{*}\times V^{*}\to\mathbb{R}$ by

$\displaystyle\check{h}(v,w)=G^{ab}G^{cd}v_{a}h_{bc}w_{d}\,,$ (317)

which for fermionic modes yields a bi-linear anti-symmetric form:

$\displaystyle\check{h}(v,w)=-\check{h}(w,v)\,.$ (318)

Using the fact that $K$ as defined in (197) commutes with $J$, it can be shown
that

$\displaystyle\check{h}(J^{\intercal}v,J^{\intercal}w)=\check{h}(v,w).$ (319)

Which allows to calculate the terms invovling the partner mode
$(\bar{x},\bar{k})$ as

$\displaystyle\check{h}(x,\bar{x})=\check{h}(k,\bar{k})=\frac{\check{h}(x,J^{\intercal}k)}{\sin
2r}\,,$ (320) $\displaystyle\check{h}(x,\bar{k})=-\cot
2r\,\check{h}(x,k)+\frac{\check{h}(x,J^{\intercal}x)}{\sin 2r}\,,$ (321)
$\displaystyle\check{h}(k,\bar{x})=-\cot
2r\,\check{h}(x,k)+\frac{\check{h}(k,J^{\intercal}k)}{\sin 2r}\,,$ (322)
$\displaystyle\check{h}(\bar{x},\bar{k})=\frac{\cot 2r}{\sin
2r}\left(\check{h}(x,J^{\intercal}x)+\check{h}(k,J^{\intercal}k)\right)$
$\displaystyle\qquad\qquad-(1+2\cot^{2}(2r))\check{h}(x,k)\,.$ (323)

Then, by equating the obtained matrix with (290), we obtain the spectrum of
$\hat{H}_{A\bar{A}}$ with
$\epsilon_{\pm}=\frac{1}{2}(\epsilon_{2}\pm\epsilon_{1})$ as

$\displaystyle\epsilon_{+}^{2}$
$\displaystyle\\!=\\!\frac{\left(\check{h}(\bar{x},\bar{k})+\check{h}(x,k)\right)^{2}+\left(\check{h}(x,\bar{k})+\check{h}(k,\bar{x})\right)^{2}}{4},$
(324) $\displaystyle\epsilon_{-}^{2}$
$\displaystyle\\!=\\!\frac{\left(\check{h}(\bar{x},\bar{k})-\check{h}(x,k)\right)^{2}+\left(\check{h}(x,\bar{k})-\check{h}(k,\bar{x})\right)^{2}}{4}$
$\displaystyle\qquad+\frac{\check{h}(x,\bar{x})^{2}}{4}\,.$ (325)

The two parameters $\phi$ and $\theta$ in the transformation $\mathbf{T}_{r}$
(mapping the eigenmodes of $\hat{H}_{A\bar{A}}$ to the partner modes $A$ and
$\bar{A}$) are determined by

$\displaystyle\tan\phi=\frac{2\check{h}(x,\bar{x})}{\check{h}(x,\bar{k})-\check{h}(k,\bar{x})}\,,$
(326) $\displaystyle\cos
2\theta=\frac{\check{h}(\bar{x},\bar{k})-\check{h}(x,k)}{2\epsilon_{-}}\,.$
(327)

With this information on the spectrum and the parameters it is possible to
check if a given pair of partner modes minimizes the energy cost: If we find
that the extracted amount of entanglement is low enough such that
$\cos{2r}\geq\frac{\epsilon_{-}}{\epsilon_{+}}$ then the pair achieves the
minimal energy cost of $\Delta E_{\text{min}}=2\epsilon_{+}\sin^{2}r$
independent of the value of $\theta$. However, beyond this regime the minimal
energy cost depends on the value of $\theta$ and according to (5.3.2) can be
lowered to

$\displaystyle\Delta E_{\theta}=\epsilon_{+}-\epsilon_{-}|\cos 2\theta|\,,$
(328)

_i.e._ , only if $\theta\in\\{0,\pi/2,\pi,\dots\\}$ the minimal energy cost of
$\Delta E_{\text{min}}=\epsilon_{+}-\epsilon_{-}=\epsilon_{1}$ is achieved by
the chosen partner modes $A$ and $\bar{A}$.

Finally, also for fermionic modes we can give an upper bound on the energy
cost of entanglement extraction from arbitrary partner modes of the system.
Combining (305) and (153) (which as shown in Appendix B applies to fermionic
Hamiltonians in the same form), we find that the energy cost to replace the
state of an arbitrary pair of partner modes by the product states of their
one-mode restricted ground states

$\displaystyle\Delta E\leq 2\epsilon_{+}\sin^{2}r\leq
2\,\omega_{N}\sin^{2}r=:\Delta E_{\mathrm{max}}$ (329)

is again upper bounded by the energy cost for a degenerate Hamiltonian with
only the system’s largest excitation energy $\omega_{N}$ in its spectrum.

##  6 Applications

In this section, we apply the previously derived minimal energy cost to
concrete physical models. The main goal is to present a proof of concept that
our general results are easily related to concrete scenarios of extracting
entanglement from modes that are accessible in a physical model. In
particular, we present examples for both bosonic and fermionic systems, the
latter being also related to spin systems via Jordan-Wigner transformation.

###  6.1 Hamiltonian of dilute Boson gas

We consider the Hamiltonian of the form

$\displaystyle\hat{H}=\sum_{k\neq
0}\omega_{k}a_{k}^{\dagger}a_{k}+\gamma_{k}\left(a_{k}^{\dagger}a_{-k}^{\dagger}+a_{k}a_{-k}\right)\,,$
(330)

with $\omega_{k}=\omega_{-k}$ and $\gamma_{k}=\gamma_{-k}$, which is well
known to describe a weakly interacting dilute Boson gas (see, _e.g._ , [74]).
The sum runs over all modes with non-zero momentum, which have very low
occupation numbers, but it excludes the zero momentum mode, which is
macroscopically occupied. In this way, the Hamiltonian describes the
interaction of excitations in the higher modes with the condensate of
particles in the zero-mode.

Typically, the squeezing term in the Hamiltonian is much smaller than the
number operator term, _i.e._ , $\gamma_{k}\ll\omega_{k}$. However, as long as
$\gamma_{k}<\omega_{k}/2$ the Hamiltonian is bounded from below, as we shall
see below.

Inserting
$a_{k}=\frac{1}{\sqrt{2}}\left(\hat{q}_{k}+\mathrm{i}\hat{p}_{k}\right)$, we
can rewrite the Hamiltonian as

$\displaystyle\hat{H}$ $\displaystyle=\sum_{k>0}\hat{H}_{k}-\sum_{k\neq
0}\frac{\omega_{k}}{2}$ (331)

with

$\displaystyle\hat{H}_{k}$
$\displaystyle=\frac{\omega_{k}}{2}\left(\hat{q}_{k}^{2}+\hat{p}_{k}^{2}+\hat{q}_{-k}^{2}+\hat{p}_{-k}^{2}\right)$
(332)
$\displaystyle\qquad\qquad+2\gamma_{k}\left(\hat{q}_{k}\hat{q}_{-k}-\hat{p}_{k}\hat{p}_{-k}\right)\,.$
(333)

It is evident, that modes with opposite momentum, _i.e._ , pairs of the form
$(\hat{q}_{k},\hat{p}_{k},\hat{q}_{-k},\hat{p}_{-k})$ are partner modes in the
ground state of $\hat{H}$ because the blocks $\hat{H}_{k}$ couple them
pairwise. Since we assumed $\omega_{k}=\omega_{-k}$ each block $\hat{H}_{k}$
is a degenerate two-mode Hamiltonian. In fact, the Hamiltonian matrix
corresponding to $\hat{H}_{k}$ is

$\displaystyle\mathbf{h}_{k}$
$\displaystyle=\begin{pmatrix}\omega_{k}&0&2\gamma_{k}&0\\\
0&\omega_{k}&0&-2\gamma_{k}\\\ 2\gamma_{k}&0&\omega_{k}&0\\\
0&-2\gamma_{k}&0&\omega_{k}\end{pmatrix}\,,$ (334)

which matches (115), with the squeezing parameter $r$ and the excitation
energy of $\hat{H}_{k}$ being

$\displaystyle\epsilon\operatorname{ch}{2r}=\omega_{k}\,,\quad\epsilon\operatorname{sh}{2r}=-2\gamma_{k}\,,$
$\displaystyle\Leftrightarrow\epsilon=\sqrt{\omega_{k}^{2}-4\gamma_{k}^{2}},\quad\operatorname{ch}{2r}=\frac{1}{\sqrt{1-4\gamma_{k}^{2}/\omega_{k}^{2}}}\,.$
(335)

This means that in order to extract pure state entanglement from the ground
state of the Hamiltonian $\hat{H}$, one needs to select two modes with
opposite momentum $k$ and $-k$. Following (58), these modes share an
entanglement entropy of

$\displaystyle\Delta
S=s_{b}\left(\left(1-4\gamma_{k}^{2}/\omega_{k}^{2}\right)^{-\frac{1}{2}}\right)\,.$
(336)

The minimal energy cost of entanglement extraction from these modes is,
according to (122),

$\displaystyle\Delta
E=\epsilon\left(\operatorname{ch}{2r}-1\right)=\omega_{k}\left(1-\sqrt{1-4\gamma_{k}^{2}/\omega_{k}^{2}}\right)\,.$
(337)

Figure 6: Ratio of energy cost to extracted entanglement from the ground state
of a dilute weakly interacting Boson gas (Hamiltonian of (330)) as a function
of $\gamma_{k}/\omega_{k}$.

It is interesting to see that while the entropy between the partner modes
diverges as $\gamma_{k}\to\omega_{k}/2$, the energy cost is upper bounded and
approaches $\omega_{k}$. At first this may seem in contradiction to the
results of Section 3, but it is due to the excitation energy of $\hat{H}_{k}$
vanishing, $\epsilon\to 0$, in this limit. In particular, this leads to an
interesting dependency of the ratio $E_{\mathrm{min}}/\Delta S$ between energy
cost and extracted entanglement on the ratio $\gamma_{k}/\omega_{k}$ of the
Hamiltoian parameters, seen in Figure 6. The energy that needs to be invested
per bit of extracted entanglement falls off drastically both for very small
amounts of entanglement, $\gamma_{k}\to 0$, and for extremely large amounts of
entanglement as $\gamma_{k}\to\omega_{k}/2$. In between the highest amount of
energy per extracted entanglement entropy is required at
$\gamma_{k}/\omega_{k}\approx 0.389368$, where the entanglement entropy is
$\Delta S\approx 1.00679$ and the energy cost $E_{\mathrm{min}}\approx
0.372648\,\omega_{k}$.

An interesting feature of the entanglement structure of this example system is
that it exhibits partner modes which may be realistically accessible in
experiment. In other systems it seems much more challenging to access both
modes in a pair of partner modes equally: For example, in the coupled harmonic
chain generally (at least) one of the partner modes are delocalized over the
complete system. The Hamiltonian (330), however, couples and entangles modes
of counter-propagating momentum symmetrically such that both partner modes
should be equally accessible in experiment.

###  6.2 XY spin model

Figure 7: Dispersion relation $\varepsilon_{k}$ of the XY model for different
values of $h/J$ and $\gamma$.

The XY model presents a prime example of a spin model that can be mapped to a
fermionic quadratic Hamiltonian via the Jordan-Wigner transformation [75]. We
consider the Hamiltonian with an external field $h$ given by [76]

$\displaystyle\begin{split}\hat{H}_{\mathrm{XY}}&=-\frac{J}{2}\sum^{N}_{j=1}\left[(1+\gamma)\hat{S}^{\mathrm{X}}_{j}\hat{S}^{\mathrm{X}}_{j+1}\right.\\\
&\qquad\quad\left.+(1-\gamma)\hat{S}^{\mathrm{Y}}_{j}\hat{S}^{\mathrm{Y}}_{j+1}\right]-h\hat{S}^{\mathrm{Z}}_{j}\,,\end{split}$
(338)

where $\hat{S}_{j}^{\mathrm{X}}$, $\hat{S}_{j}^{\mathrm{Y}}$ and
$\hat{S}_{j}^{\mathrm{Z}}$ represent the spin-$\frac{1}{2}$ operators on site
$j$. Ignoring boundary terms that are negligible in the thermodynamic limit,
the XY model can be mapped to the quadratic fermionic Hamiltonian given by
[77, 78]

$\displaystyle\begin{split}\hat{H}_{\mathrm{XY}}=\sum_{k\in\mathcal{K}}\varepsilon_{\kappa}\
\left(\hat{\eta}_{\kappa}^{\dagger}\hat{\eta}_{\kappa}-\frac{1}{2}\right)\qquad\text{with}\qquad\\\
\varepsilon_{\kappa}\\!=\\!\sqrt{h^{2}\\!+2hJ\cos(\kappa)\\!+J^{2}\\!+(\gamma^{2}-1)J^{2}\sin(\kappa)^{2}}\,,\end{split}$
(339)

where we defined $\mathcal{K}=\\{2\pi k/N|0\leq k<N\\}$. Figure 7 shows the
dispersion relation $\epsilon_{\kappa}$ as a function of $\kappa$ for various
parameter choices. The fermionic annihilation operators $\hat{\eta}_{\kappa}$
diagonalize the Hamiltonian and are related to on-site annihilation operators

$\displaystyle\hat{f}_{j}=\frac{1}{\sqrt{N}}\sum_{\kappa\in\mathcal{K}}e^{-\mathrm{i}\kappa}(u_{\kappa}\eta_{\kappa}-v_{\kappa}\eta_{-\kappa}^{\dagger})$
(340)

on site $j$, where the Bogoliubov coefficients $u_{\kappa}$ and $v_{\kappa}$
are given by

$\displaystyle\begin{split}u_{\kappa}=\frac{\varepsilon_{\kappa}+a_{\kappa}}{\sqrt{2\varepsilon_{\kappa}(\varepsilon_{\kappa}+a_{\kappa})}}\,,\quad
v_{\kappa}=\frac{\mathrm{i}b_{\kappa}}{\sqrt{2\varepsilon_{\kappa}(\varepsilon_{\kappa}+a_{\kappa})}}\,,\\\
a_{\kappa}=-J\cos(\kappa)-h\,,\quad\text{and}\quad
b_{\kappa}=\gamma\,J\sin(\kappa)\,.\end{split}$ (341)

Let us emphasize that the Jordan-Wigner transformation relating the spin
Hamiltonian (338) with the fermion Hamiltonian (339) is non-local. However, it
preserve bi-partite entanglement [15] between localized regions, _i.e._ , the
entanglement entropy associated to a single fermionic site is equal to the
respective entanglement entropy of the same site in the spin Hamiltonian
(338). (The relationship for general non-local subsystems is more subtle.)
Also, we note that the Hamiltonian in (339) is a paradigmatic model in its own
right: Choosing $\gamma=1$, we obtain the transverse field Ising model [79].

To find the minimal excitation energy, _i.e._ , the lowest eigenvalue of
$\hat{H}_{XY}$, we need to find the minimum of the dispersion relation
$\varepsilon_{\kappa}$. We find that the condition
$\cos\kappa=\frac{h}{J(\gamma^{2}-1)}$ ensures that $\varepsilon^{2}_{\kappa}$
is minimal, which yields the lowest eigenvalue

$\displaystyle\epsilon_{\min}=\min_{\kappa}\varepsilon_{\kappa}=\sqrt{\gamma^{2}\left(J^{2}+\frac{h^{2}}{\gamma^{2}-1}\right)}\,.$
(342)

The model becomes gapless for $\gamma=0$ and for $h=J\sqrt{1-\gamma^{2}}$. If
the number of sites $N$ is large enough, there are sufficiently many
excitation energies close to the minimum that we can consider the spectrum as
flat. This means, that by extracting a pair of partner modes from the subspace
of these lowest energy modes, we can achieve the minimal energy cost

$\displaystyle\Delta E_{\mathrm{min}}=2\epsilon_{\min}\sin^{2}{r}\,,$ (343)

with the partner modes’ squeezing parameter $r$.

The ground state $\ket{\Omega}$ of this model is encoded by the complex
structure [69]

$\displaystyle
J\equiv\mathbf{J}=-\frac{\mathrm{i}}{N}\sum_{{\kappa}\in\mathcal{K}}\left(\begin{array}[]{c|c}|u_{\kappa}|^{2}e^{\mathrm{i}{\kappa}(j-l)}-|v_{\kappa}|^{2}e^{\mathrm{i}{\kappa}(l-j)}&u^{*}_{\kappa}v^{*}_{\kappa}\left(e^{\mathrm{i}{\kappa}(l-j)}-e^{\mathrm{i}{\kappa}(j-l)}\right)\\\
\hline\cr
u_{\kappa}v_{\kappa}\left(e^{\mathrm{i}{\kappa}(l-j)}-e^{\mathrm{i}{\kappa}(j-l)}\right)&-|u_{\kappa}|^{2}e^{\mathrm{i}{\kappa}(l-j)}+|v_{\kappa}|^{2}e^{\mathrm{i}{\kappa}(j-l)}\end{array}\right)\,,$
(346)

here expressed with respect to the local basis
$\hat{\xi}^{a}\equiv(\hat{f}_{1}^{\dagger},\dots,\hat{f}_{N}^{\dagger},\hat{f}_{1},\dots,\hat{f}_{N})$.
We can change to a hermitian basis defined by
$\hat{q}_{i}=\frac{1}{\sqrt{2}}(\hat{f}_{i}^{\dagger}+\hat{f}_{i})$ and
$\hat{p}_{i}=\frac{\mathrm{i}}{\sqrt{2}}(\hat{f}_{i}^{\dagger}-\hat{f}_{i})$,
such that $G\equiv\mathbb{1}$. In this basis, the only non-vanishing entries
of the covariance matrix are given by

$\displaystyle\begin{split}\Omega(q_{j},p_{l})&=\frac{1}{L}\sum_{\kappa\in\mathcal{K}}\Big{[}(|v_{\kappa}|^{2}-|u_{\kappa}|^{2})\cos{\kappa(j-l)}\\\
&\qquad\qquad\,\,+2\,\mathrm{Im}(u_{k}v_{k})\sin{\kappa(j-l)}\Big{]}\,.\end{split}$
(347)

Similarly, we find the Hamiltonian matrix to be

$\displaystyle\check{h}(q_{j},p_{l})=h\delta_{jl}+\frac{J+\gamma}{2}\delta_{j+1,l}+\frac{1-\gamma}{2}\delta_{j-1,l}\,,$
(348)

where we have periodicity $N+1\equiv 1$.

We can now explore how the energy cost of randomly chosen partner modes relate
to the minimal energy cost $\Delta E_{\mathrm{min}}$ and the upper bound
$\Delta E_{\mathrm{max}}$ established above. To this end, we choose a Haar
random mode $(x_{a},k_{a})$ spanning a subsystem $A\subset V$ and then apply
the following steps:

  1. 1.

Compute entanglement entropy $S_{A}(\ket{\Omega})$.

  2. 2.

Compute partner mode $(\bar{x},\bar{k})$.

  3. 3.

Compute the restriction $\hat{H}_{A\bar{A}}$.

  4. 4.

Compute the energy cost $\Delta E$ for replacing the partner mode state by the
ground state $\ket{\Omega_{A}^{\prime}}\otimes\ket{\Omega^{\prime}_{\bar{A}}}$
of $\hat{H}_{A}+\hat{H}_{\bar{A}}$.

Selecting a random mode with respect to the Haar measure of the orthogonal
group can be accomplished by generating a random $2N$-by-$2N$ matrix with
Gaussian distributed entries and orthonormalize its column vectors with
respect to the inner product induced by $G^{ab}$. From here, we can select the
first two column vectors. This is equivalent to directly generating two
vectors and orthonormalizing them.

Figure 8: Copmarison of lower and upper bounds for the energy cost of
entanglement extraction with a sample of $10,000$ Haar randomly generated
partner modes. We consider the XY model with $N=10$, and parameter choices
$(h=.5,\gamma=.5)$ and $(h=0,\gamma=.8)$.

As seen in Figure 8, the Haar random sample is accurately bounded from below
by our predicted minimal energy cost. The upper bound on the energy cost is
due to the fact that the spectrum of the XY model is also bounded from above
by

$\displaystyle\epsilon_{\max}=\max_{\kappa}\epsilon_{\kappa}=\sqrt{J^{2}+h^{2}+2hJ}\,.$
(349)

We see that our bounds become tighter for narrower excitation spectra. As seen
from Figure 7 our spectrum becomes flat for $h\to 0$ and $\gamma\to 1$. This
is reflected in Figure 8 for $(h=0,\gamma=.8)$.

Our analysis of the energy cost for the XY model did not take into account
locality in real space so far. In particular, the optimal partner modes that
accomplish our minimal energy cost are spanned by the de-localized energy
eigenmodes with lowest energy values. It is therefore interesting, to
investigate the energy cost in a more physical scenario where one mode is
localized on a given site, while its partner mode is localized in the
complement of the site. In this framework the choice of modes is completely
fixed, hence we can scan through the $(h,\gamma)$ parameter space of the XY
model to compare the energy cost for this single site extraction to the global
lower and upper bounds.

In Figure 9, we show the energy cost of entanglement extraction for a single-
site mode and its partner. We scan through the parameter space along three
different paths of points, which are constrained to yield the same lowest
excitation energy $\epsilon_{\min}$ for the Hamiltonian $\hat{H}_{XY}$, _i.e._
,

$\displaystyle
h=\pm\frac{1}{\gamma}{\sqrt{\left(1-\gamma^{2}\right)\left(\gamma^{2}-\epsilon_{\min}^{2}\right)}}\,.$
(350)

We notice that the entanglement entropy of a single site with the rest of the
system is close to maximal, _i.e._ , $\Delta S\approx 1$. This is not
surprising, as it is known [67] that every local site becomes maximally
entanglement in the thermodynamic limit $N\to\infty$. Furthermore, we see that
for $\epsilon_{\min}=.9$, the actual energy cost from the single site
extraction is comparably close to our lower bound. This indicates that for
certain parameter choices our minimal energy cost could actually be achieved
in this concrete physical scenario.

Figure 9: We consider the XY model with $N=10$ and choose the model parameters
$h/J$ and $\gamma$, such that the minimal excitation energy $\epsilon_{\min}$
is in $\\{.5,.7,.9\\}$ as indicated in the upper panel. For each point in the
parameter space, we then compute the energy cost of entanglement extraction
between a single site mode and its partner. We compare these sample points
with our lower bounds for the given energy.

##  7 Discussion and outlook

In this work, we have found the minimal energy cost of extracting a pair of
entangled partner modes from the ground state of a quadratic Hamiltonian, and
we have characterized the modes which achieve this minimal energy cost.

The minimal energy cost is linear in the spectrum of the Hamiltonian
restricted to the two-mode space spanned by the extracted partner modes.
Hence, as one may intuitively expect, the globally lowest energy cost possible
is achieved by partner modes constructed from the two lowest energy eigenmodes
of the Hamiltonian. In some applications, access to such less localized energy
eigenmodes may physically not be possible. However, our analysis of the XY
spin model showed instances of modes that are localized on a single site but
still came relatively close to the global minimum, with their delocalized
partner mode.

Interestingly, the minimal energy cost of entanglement extraction can be used
to derive an upper bound on a Hamiltonian’s spectral gap: If it is possible to
extract a pair of partner modes sharing entanglement $S$ from the ground
state, while only increasing the system’s energy by an amount $\Delta E$, then
the Hamiltonian’s lowest eigenvalue is at most as high as the value $\epsilon$
implied by the minimal energy cost relation derived above.

Apart from the question of the global minimal energy cost, as dictated by the
spectrum of the full Hamiltonian, it is meaningful to ask whether a given pair
of partner modes minimize the energy cost locally, _i.e._ , whether they
achieve the minimal energy cost allowed for by the spectrum of the Hamiltonian
restricted to their two-mode subspace. If not, then there exist symplectic
transformations acting only on the given partner modes which yield another
pair of partner modes achieving the local minimum.

We showed that for bosonic modes there are only specific pairs of partner
modes which achieve the local minimal energy cost. These modes are obtained by
squeezing the odd and even superpositions of the eigenmodes of the restricted
Hamiltonian. However, for fermionic modes, if the amount of extracted
entanglement lies below a critical value, all partner modes achieve the local
minimal energy cost. Only above the critical value, specific partner modes are
required which are obtained by squeezing the eigenmodes of the restricted
Hamiltonian. The latter is due to the fact that above the critical value, the
lowest excited energy eigenstate of the two-mode subsystem then coincides with
the product of the restricted ground states of the partner modes.

An intriguing question for future research will be to analyze the dynamics of
the extraction process due to the time evolution of the source modes under the
source system’s Hamiltonian. In certain cases this dynamic may be relatively
simple. For example, the partner modes which minimize the energy cost lie