Minimal energy cost of entanglement extraction

TL;DR

This paper calculates the lowest possible energy increase needed to extract entanglement from the ground state of quadratic bosonic or fermionic systems, providing a universal lower bound applicable to many physical systems.

Contribution

It introduces a protocol-independent method to determine the minimal energy cost of entanglement extraction in quadratic systems, establishing a fundamental lower bound.

Findings

Derived the minimal energy increase $\, ext{Δ}E_{ ext{min}}$ for entanglement extraction.

Constructed modes that achieve this minimal energy cost.

Applied results to various physical systems with examples.

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 $\Delta E_{\mathrm{min}}$ 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.

Equations

$$v^a\equiv {\left(q_1(v),p_1(v),\cdots ,q_N(v),p_N(v)\right)}^{⊺}\in {\mathbb{R}}^{2N},$$

$${\Omega }^{ab}$$

$${\mathbf{\Omega }}_2=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right).$$

$${\mathcal{O}}_w:V\to \mathbb{R}:v^a\mapsto w_a{v}^a.$$

$${\mathcal{O}}_h:V\to \mathbb{R}:v^a\mapsto \frac{1}{2}{h}_{ab}{v}^a{v}^b.$$

$$[{\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}1$$

$${\hat{\xi }}^a\equiv ({\hat{q}}_1,\cdots ,{\hat{q}}_N,{\hat{p}}_1,\cdots ,{\hat{p}}_N{)}^{⊺}$$

$${\hat{\mathcal{O}}}_w=w_a{\hat{\xi }}^a\text{and}{\hat{\mathcal{O}}}_h=\frac{1}{2}{h}_{ab}{\hat{\xi }}^a{\hat{\xi }}^b.$$

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

$${\mathcal{O}}_t:v\mapsto \frac{1}{n!}{t}_{a_1\cdots a_n}{\xi }^{a_1}(v)\cdots {\xi }^{a_n}(v)$$

$${\xi }^a\equiv (q_1,p_1,\cdots ,q_N,p_N{)}^{⊺}.$$

$$\{{\xi }^a,{\xi }^b\}={\Omega }^{ab}.$$

$$\{{\mathcal{O}}_w,{\mathcal{O}}_u\}=w_a{u}_b\{{\xi }^a,{\xi }^b\}=w_a{u}_b{\Omega }^{ab}.$$

$$z^a$$

$$G^{ab}$$

$$G^{ab}{\Omega }_{bc}^{-1}{G}^{cd}{\Omega }_{de}^{-1}=-{\delta }^a{}_e.$$

$$J^a{}_b=-G^{ac}{\Omega }_{cb}^{-1}.$$

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

$$\frac{1}{2}({\delta }^a{}_b+\mathrm{i}{J}^a{}_b)({\hat{\xi }}^b-z^b)|G,z⟩=0.$$

$$C_n^{a_1\cdots a_n}=⟨G,z|({\hat{\xi }}^{a_1}-z^{a_1})\cdots ({\hat{\xi }}^{a_n}-z^{a_n})|G,z⟩.$$

$$C_2^{ab}=\frac{1}{2}(G^{ab}+\mathrm{i}{\Omega }^{ab}).$$

$$C_n^{a_1\cdots a_n}$$

$$=C_2^{a_1{a}_2}\cdots C_2^{a_{n-1}{a}_n}+\cdots ,$$

$$C_4^{abcd}=C_2^{ab}{C}_2^{cd}+C_2^{ac}{C}_2^{bd}+C_2^{ad}{C}_2^{bc}.$$

$${\Omega }^{ab}\equiv \left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)\text{and}{G}^{ab}\equiv \left(\begin{array}{cc}1/\omega & 0\\ 0& \omega \end{array}\right).$$

$$J^a{}_b=-G^{ac}{\Omega }_{cb}^{-1}\equiv \left(\begin{array}{cc}0& 1/\omega \\ -\omega & 0\end{array}\right).$$

$$\frac{1}{2}({\delta }^a{}_b+\mathrm{i}{J}^a{}_b){\hat{\xi }}^b\equiv \frac{1}{2}\left(\begin{array}{c}\hat{q}+\frac{\mathrm{i}}{\omega }\hat{p}\\ \hat{p}-\mathrm{i}\omega \hat{q}\end{array}\right)=\left(\begin{array}{c}\sqrt{\frac{1}{2\omega }}\hat{a}\\ -\mathrm{i}\sqrt{\frac{\omega }{2}}\hat{a}\end{array}\right),$$

$$\hat{H}=\frac{1}{2}{h}_{ab}{\hat{\xi }}^a{\hat{\xi }}^b+f_a{\hat{\xi }}^a,$$

$$E=\frac{1}{4}{h}_{ab}{G}^{ab}+\frac{1}{2}{h}_{ab}{z}^a{z}^b+f_a{z}^a.$$

$$z^a=-(h^{-1}{)}^{ab}{f}_b,$$

$$G^{ab}=-K^a{}_c|K^{-1}{|}^c{}_d{\Omega }^{db},$$

Full text

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

With this, the projection of equation (22) reads

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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})$

[TABLE]

Applying formula (34) yields

[TABLE]

which agrees with the result from (28). Here

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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]

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

[TABLE]

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

[TABLE]

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:

**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

[TABLE]

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

[TABLE]

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.

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where the individual blocks are given by

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

The energy variance of the final system state is

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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$.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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.

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Wick’s theorem applies in almost the same way

[TABLE]

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

[TABLE]

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

[TABLE]

With this, the fermionic projector becomes

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

which gives via (197) rise to

[TABLE]

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

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

then its partner mode is given by

[TABLE]

This means that the partner mode is in standard form itself

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

where

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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)

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

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]

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

with the partner modes’ squeezing parameter $r$.

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

[TABLE]

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

[TABLE]

Similarly, we find the Hamiltonian matrix to be

[TABLE]

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:

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.

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

[TABLE]

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.,

[TABLE]

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.

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 within the subspace of only the two lowest energy eigenmodes of the system. Hence their time evolution is restricted to this subspace. In general, for partner modes constructed from a range of energy eigenmodes, a more complex time evolution is expected. This would affect how the extraction from such source modes can be practically implemented, as well as it affects how the energy injected during the extraction process dissipates inside the source system after the extraction process. Indeed, the prospect that such complex dynamics potentially could be exploited in the design of entanglement extraction protocols motivates further study of the dynamics and time evolution of partner mode pairs.

7.1 Outlook: Extraction from arbitrary modes

In this work, we assumed that the system modes which are swapped onto the target modes can be selected freely from the source system. This led us to choose a pair of partner modes which is the only way to obtain a pure final state in the target modes. Hence, we derived the energy cost of replacing the state of arbitrary partner modes in a system’s ground state with the product state of their restricted ground states. In particular, we established an upper bound to this cost in terms of the systems largest excitation energy.

In practical implementations of entanglement extraction protocols, however, the access to the source system might be limited. For example, the modes could be restricted to be localized in a given subregion of the system, or to lie in the span of a limited set of modes. In particular, such restriction could make it impossible to access pairs of partner modes. It could also be interesting to consider randomly chosen modes, see, e.g., [80].

This raises the important question what our results imply for more general extraction protocols that are not based on partner modes. We expect that the minimal energy cost for entanglement extraction from partner modes, that we derived here, also constitutes a lower bound for the minimal energy cost for swapping pairs of arbitrary single modes $A$ and $B$ out of the system that are not necessarily partners.

If the modes $A$ and $B$ are not partner modes then other entanglement monotones than the entanglement entropy $S(\rho_{A})$ or $S(\rho_{B})$ are necessary to measure the entanglement between them, such as logarithmic negativity or entanglement of purification need to be used instead. This is because the modes $A$ and $B$ share entanglement with the remainder of the system, such that their joint state $\rho_{AB}$, i.e., the state of the subsystem $A\oplus B$, is mixed.

Whichever mixed state entanglement $\mathcal{N}_{AB}$ measure we employ, we may assume that the lower entanglement entropy of the two modes $S_{\min}=\min(S(\rho_{A}),S(\rho_{B}))$ implies an upper bound on $\mathcal{N}_{AB}$ because the entanglement between a mode $A$ and its complement is larger than the entanglement between modes $A$ and $B$. In particular, this means that a pair of partner modes sharing entanglement $S_{\min}$ is at least as valuable as the two-mode state $\rho_{AB}$.

However, we expect that the minimal energy cost for the extraction of a partner mode pair $A$ and $\bar{A}$ with entanglement $S_{\min}$ is always lower than the energy cost of extracting modes $A$ and $B$. In fact, this can be proven for the case of a Hamiltonian with degenerate excitation energies and we have numerical evidence for the general case [81].

In this way, the derived minimal energy cost from partner mode extraction actually applies much more broadly to essentially any entanglement monotone (of Gaussian states) and any two target modes $A$ and $B$. Note, however, that this minimal energy cost will only be saturated if we are actually using partner modes, which are in this sense optimal. Where this is prevented by additional constraints, e.g., locality assumptions, we are likely to underestimate the actually required energy cost.

7.2 Outlook: Mixed states

Furthermore, it will be interesting to consider extraction from states other than the ground state of the source system. Of particular relevance for practical implementations should be general mixed states and thermal states. There, we may encounter situations with vanishing or even negative energy cost. Similarly, in quantum systems with spontaneously broken symmetries, i.e., with several ground states, one might be able to extract entanglement by moving from one ground state to another.

In the context of mixed states, but also in general, it should be relevant to consider more general methods of entanglement extraction. Instead of replacing the entangled state of a subsystem by a product state in one discrete step, one could consider more general local operations on the subsystem such as the local operations and classical communication (LOCC) envisioned in the general framework of [48]. This may yield an infinitesimal relation between entropy and energy increase around a given system state.

For these and other future research directions, the advancement of the partner mode formula provided in this work should prove instrumental because it can yield better insight into the correlation structure of mixed Gaussian states. In particular, we are interested in extending the partner mode formula to mixed states. This allows, e.g., for the decomposition of a mixed Gaussian states into chains of $k$-tuples of modes where the modes in each $k$-tuple are not correlated with each other, and only correlated with modes in the preceding and subsequent tuple. In fact, this construction does not only apply to mixed states but to bi-linear forms in general. Thus, it may be used to further develop techniques such as [82] to efficiently describe the interaction of several quantum systems with a given bath of harmonic oscillators, with applications ranging from open system dynamics to fundamental quantum communication between local observers via relativistic fields [83, 84, 85, 86, 87, 88].

Acknowledgements

The authors thank Dmytro Bondarenko, Tommaso Guaita, Terry Farrelly, and Tommaso Roscilde for inspiring discussions. LH is supported by the the Max Planck Harvard Research Center for Quantum Optics and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2111 – 39081486. LH thanks the QMATH center for their hospitality during several visits related to this project. RHJ acknowledges support by ERC Advanced grant 321029 and VILLUM FONDEN via the QMATH center of excellence (grant no.10059).

Appendix A Standard forms of Gaussian states and the complex structure

The linear complex structure associated to a quantum state is a powerful method to parametrize Gaussian states. This is because most properties of the partial, mixed state of subsystems, such as the von Neumann entropy and the full entanglement spectrum, are encoded in the restriction of the complex structure to the relevant sub-phase space. Therefore, this appendix uses the complex structure to rephrase some known properties of Gaussian states in a largely basis-independent way. The key properties of the complex structure are reviewed and derived based on compatibility with the symplectic form (for bosons) or the positive-definite metric (for fermions). This appendix is in parts based on the respective appendix in [59].

Definition 1**.**

Given a finite-dimensional real vector space $V$, we refer to a linear map $J:V\to V$ as complex structure if it satisfies $J^{2}=-\mathbb{1}$.

The spectrum of a linear complex structure takes the following simple form.

Proposition 1**.**

A complex structure is diagonalizable over the complex numbers with eigenvalues that come in conjugate pairs as $\pm\mathrm{i}$ with equal multiplicity. This also implies that $V$ must be even dimensional.

Proof.

The equality $J^{2}=-\mathbb{1}$ implies that $J$ is diagonalizable and every eigenvalue squares to $-1$. Moreover, $J$ is a real map which means that all non-real eigenvalues come in complex conjugate pairs of equal multiplicity. Thus, the spectrum of $J$ is given by $\pm\mathrm{i}$ with equal multiplicity equal to half of the dimension of $V$. The last part implies that such a linear map $J$ can only exist in an even dimensional real vector space $V$. ∎

We will be particularly interested in the restriction of linear complex structures.

Proposition 2**.**

Given a complex structure $J$ over $V$ and a direct sum decomposition of $V=A\oplus B$, we have the unique decomposition

[TABLE]

*where the projections $\mathbb{P}_{A}:V\to A$ and $\mathbb{P}_{B}:V\to B$ with $\mathbb{P}_{A}+\mathbb{P}_{B}=\mathbb{1}$ provide the unique decomposition of a vector $v=a+b$ into its part in $a=p_{A}(v)\in A$ and $b=p_{B}(v)\in B$.

Given such a decomposition, the spectrum of $J^{2}_{A}$ and $J^{2}_{B}$ is the same except for the number of eigenvalues equal to $-1$. Put differently, given an eigenvalue $\lambda$ of $J^{2}_{A}$, we have $\lambda=-1$ or $J^{2}_{B}$ also has $\lambda$ as eigenvalue.*

Proof.

We write $J^{2}=-\mathbb{1}$ in blocks to find

[TABLE]

Let us consider the eigenvector $a\in A$ with $J_{A}^{2}\,a=\alpha\,a$. From the first diagonal block equation, we find that $a$ must also be eigenvector of $J_{AB}J_{BA}$ with $J_{AB}J_{BA}\,a=-(1+\alpha)a$. Let us define $b:=J_{BA}\,a$. Now there are two possibilities: Either $b=0$ which implies $\alpha=-1$ or $b\neq 0$, in which case we must have $J_{BA}J_{AB}\,b=-(1+\alpha)J_{BA}\,a=-(1+\alpha)\,b$. The second diagonal block equation then implies $J_{B}^{2}\,b=\alpha\,b$. ∎

A.1 Compatible symplectic form (bosons)

For bosonic systems, we are interested in linear complex structures that are compatible with the underlying symplectic form or vice versa.

Definition 2**.**

An antisymmetric, non-degenerate bilinear form $\omega:V\times V\to\mathbb{R}$ is called compatible to a complex structure $J:V\to V$ if and only if it satisfies the following conditions:

• (I)

Skew-symmetry: The symplectic form satisfies $\omega(Jv,w)=-\omega(v,Jw)$.

• (II)

Taming: The bilinear form $g(v,w):=\omega(v,Jw)$ is positive definite.

We will refer to $\omega$ as a $J$-compatible symplectic form and to $J$ as a $\omega$-compatible complex structure.

The compatibility conditions implies the following invariance property.

Proposition 3**.**

A $J$-compatible symplectic form $\omega$ is invariant under $J$ meaning

[TABLE]

for all $v,w\in V$, thus $J\in\mathrm{Sp}(V,\omega)$. The bilinear form

[TABLE]

is symmetric and positive definite and thus a proper metric.

Proof.

Straightforward computation leads to

[TABLE]

Every linear map $M:V\to V$ that preserves the symplectic form $\omega$ in the sense of $\omega(Mv,Mw)=\omega(v,w)$ is part of the symplectic group $\mathrm{Sp}(V,\omega)$. The form $g(v,w)=\omega(v,Jw)$ is symmetric due to

[TABLE]

The taming property (II) of $\omega$ ensures that $g$ is positive definite, and thus a proper positive definite metric. ∎

The compatibility conditions characterize the possible spectra of restrictions of the linear complex structure to symplectic subspaces.

Proposition 4**.**

Let the complex structure $J$ and the symplectic form $\omega$ on a vector space $V$ be compatible, and let $V=A\oplus B$ a decomposition of $V$ into a direct sum of even dimensional symplectic complements. Then, the restricted complex structure $J_{A}$ is diagonalizable and its spectrum consists of complex conjugate pairs $\pm ic_{i}$ with $c_{i}\in[1,\infty)$. Both $J_{A}$ and $\omega_{A}$ can simultaneously be brought into the block diagonal form

[TABLE]

Proof.

The restricted linear complex structure $J_{A}$ is diagonalizable due to the fact that it squares to the diagonalizable matrix $J_{A}^{2}=-\mathbb{1}_{A}-J_{AB}J_{BA}$. Given a non-zero eigenvector $a\in A$ with $J_{A}^{2}\,a=\alpha\,a$, we have $\alpha\in(-\infty,-1]$ by the following argument: From

[TABLE]

follows $\alpha\leq 0$. Moreover, we can compute

[TABLE]

where we used in the last step that $A$ and $B$ are symplectic complements with $g(Ja,J_{BA}a)=\omega(Ja,JJ_{BA}a)=-\omega(J^{2}a,J_{BA}a)=\omega(a,J_{BA}a)=0$. We thus find $g(J_{A}a,J_{A}a)\geq g(a,a)$ which implies $\alpha\in(-\infty,-1]$.

From this, we find that $J_{A}$ must have eigenvalues appearing in complex conjugate pairs $\pm\mathrm{i}c$ with $c=\sqrt{-\alpha}$, meaning that every eigenvalue $\alpha$ of $J_{A}^{2}$ appears with even multiplicity. Since $J_{A}$ is a linear real map with complex conjugated eigenvalues $\pm\mathrm{i}c$, it can can be brought into the block diagonal form of (366). In particular, in such a basis, $J_{A}^{2}$ is diagonal. Next, we show that $\omega_{A}$ simultaneously takes the form claimed above.

While the full $J$ satisfies $\omega(Jv,Jw)=\omega(v,w)$ for $v,w\in V$, the restricted $J_{A}$ does in general not satisfy $\omega_{A}(J_{A}a,J_{A}a^{\prime})=\omega_{A}(a,a^{\prime})$ for $a,a^{\prime}\in A$. However, we can use the relation $\omega(Jv,w)=-\omega(v,Jw)$ for $v,w\in V$ to derive the relation

[TABLE]

which is the condition for $J_{A}$ to represent a symplectic algebra element, i.e., $J_{A}\in\mathfrak{sp}(2N_{A},\mathbb{R})$. For a diagonalizable symplectic algebra element with purely imaginary eigenvalues, we can always choose a basis, e.g., see proposition 3.1.18 in [89], such that the matrix representations of $J_{A}$ and $\omega_{A}$ are given by (366). ∎

This result is central for computing the entanglement spectrum of bosonic systems. In particular, the fact that the magnitude of possible eigenvalues lies in the interval $[1,\infty)$ indicates that the entanglement entropy in bosonic systems can be arbitrarily large. This can also be related to the non-compactness of symplectic group whose group elements relate different linear complex structures.

Before we construct the standard form for the matrix representation of $J$ for a general symplectic decomposition $V=A\oplus B$, we show that if $J_{A}^{2}=-\mathbb{1}_{A}$, i.e., when the restricted complex structure defines a complex structure on the subspace $A$, then the complex structure leaves $A$ and $B$ invariant.

Proposition 5**.**

In the setup of the previous Proposition 4, if $J_{A}^{2}=-\mathbb{1}$, then $J_{BA}=0$ vanishes, furthermore $J_{B}^{2}=-\mathbb{1}$ and $J_{AB}=0$.

Proof.

Because $V=A\oplus B$ is a symplectic decomposition, for all $a\in A,b\in B$

[TABLE]

hence $\omega_{A}(J_{A}a,J_{AB}b)=-\omega_{B}(J_{BA}a,J_{B}b)$. For $J_{A}^{2}=-\mathbb{1}_{A}$, we find $J_{AB}J_{BA}=0$ due to $-\mathbb{1}=J^{2}$ and $-\mathbb{1}_{A}=J_{A}^{2}+J_{AB}J_{BA}$. To show that $J_{BA}=0$, let $b:=J_{BA}a\in B$ leading to

[TABLE]

Hence $g(b,b)=0$ which implies $b=0$ and since $a\in A$ was arbitrary, this means $J_{BA}=0$.

Since $J_{A}^{2}=-\mathbb{1}$, by Proposition 2 we have $J_{B}^{2}=-\mathbb{1}$. Hence, by repeating the argument above, we obtain $J_{AB}=0$. ∎

The derived form for $J_{A}$ above induces a standard representation of $J$ which yields the standard form of the bosonic covariance matrix presented in (51).

Proposition 6**.**

As before, let the complex structure $J$ and the symplectic form $\omega$ on a vector space $V$ be compatible, and let $V=A\oplus B$ be a decomposition of $V$ into a direct sum of even dimensional symplectic complements. Without loss of generality assume that the dimension of $A$ is lower than, or equal to the dimension of $B$. Then there exists a symplectic basis with respect to which $J$ takes the form

[TABLE]

for some $r_{i}\geq 0$, where $\mathbb{A}_{2}$ and $\mathbb{S}_{2}$ are given by

[TABLE]

Proof.

The form of the upper diagonal block, representing $J_{A}$, follows from the previous proposition, by choosing $\operatorname{ch}(2r_{i})=c_{i}$. Now for each $i\in\{1,\dots,N_{A}\}$, we denote the basis vectors bringing $J_{A}$ into this standard form by $a^{(i)}_{1},a^{(i)}_{2}\in A$, i.e.,

[TABLE]

Now, if $r_{i}>0$, define

[TABLE]

These vectors lie in the eigenspace of $J_{B}^{2}$ with eigenvalue $-\operatorname{ch}(2r_{i})^{2}$. By definition, they bring the lower left off-diagonal block, representing $J_{BA}$ into the claimed form. Using $J_{B}J_{BA}=-J_{BA}J_{A}$ which follows from (354), we find that

[TABLE]

This means that the vectors also bring $J_{B}$ into standard form on the eigenspaces of $J_{B}^{2}$ of eigenvalue $-\operatorname{ch}(2r_{i})^{2}$. They also bring $\omega$ into standard form, which follows from

[TABLE]

where we used that $J_{AB}J_{BA}a_{i}=\operatorname{sh}^{2}{2r_{i}}a_{i}$ which follows from $-\mathbb{1}=J_{A}^{2}+J_{AB}J_{BA}$ in (354). This relation we also use to confirm that the defined vectors bring the lower off-diagonal block, representing $J_{BA}$ into the claimed standard form:

[TABLE]

In the case where $r_{i}=0$, Proposition 5 ensures that the matrix entries of off-diagonal blocks vanish. Hence in the subspaces of $B$ corresponding to $r_{i}=0$, and in the remaining dimensions of $B$, we can choose any basis which brings $J_{B}$ in its standard form. Thus, we obtain the claimed matrix representation for $J$. Finally, this result implies the standard form (51) of the bosonic covariance matrix $G$ by applying $G=-J\Omega$ with $\Omega\equiv\mathbf{\Omega}$ from (5). ∎

To conclude this subsection, let us comment on the connection to the standard form of the bosonic covariance matrix and partner formula of Section 2. The symplectic form $\omega$ and the metric $g$ used in this section act on the phase space $V$. In fact, they are the inverses of $\Omega^{ab}$ and $G^{ab}$ used in the main body of the article which act on $V^{*}$, the co-vector space of the phase space.

With respect to a basis in which $J$ is represented by the matrix form which we just derived, $G$ is then represented, due to equation (20), by the matrix $\mathbf{G}=-\mathbf{J}\mathbf{\Omega}$ which is exactly as in (51).

Each pair of vectors $a_{1}^{(i)},a_{1}^{(i)}\in V$, appearing in the proof above, is dual to a pair of covectors $x,k\in V^{*}$ which defines a mode. From the proof it follows that the mode’s partner mode is defined by the covectors $(\bar{x},\bar{k})$ given by

[TABLE]

which is precisely the partner formula (63).

A.2 Compatible metric (fermions)

For fermionic systems the analysis of the complex structure is very similar. Here, we study the properties of a linear complex structure which is compatible with a positive definite metric.

Definition 3**.**

A positive definite symmetric bilinear form $g:V\times V\to\mathbb{R}$ is called compatible to a complex structure $J:V\to V$ if and only if it satisfies

[TABLE]

We will refer to $g$ as a $J$-compatible metric and to $J$ as a $g$-compatible complex structure.

The compatibility conditions implies the following invariance property.

Proposition 7**.**

A $J$-compatible metric $g$ is invariant under $J$ meaning $g(Jv,Jw)=g(v,w)$ for all $v,w\in V$ which implies $J\in O(V,g)$. Moreover, we can define the antisymmetric bilinear form $\omega(v,w)=g(Jv,w)$ which is a well-defined symplectic form on $V$.

Proof.

Straightforward computation gives

[TABLE]

Every linear map $M:V\to V$ that preserves metric $g$ in the sense of $g(Mv,Mw)=g(v,w)$ is an element of the orthogonal group $O(V,g)$. A symplectic form on a vector space $V$ is required to be (a) antisymmetric and (b) non-degenerate:

• (a)

Antisymmetry: We compute $\omega(v,w)=g(Jv,w)=g(w,Jv)=-g(Jw,v)=-\omega(w,v)$.

• (b)

Non-Degeneracy: We need to show that for every non-zero vector $v$, there is at least one vector $w$, such that $\omega(v,w)\neq 0$. This can easily be done by choosing $w=Jv$, where we find $\omega(v,w)=g(Jv,Jv)>0$.

This proves that every compatible every pair of a complex structure $J$ and a metric $g$ defines a symplectic form $\omega$. ∎

Again, the compatibility condition characterizes the possible spectra of restrictions of the linear complex structure to even dimensional subspaces.

Proposition 8**.**

Given a compatible pair of a complex structure $J$ and a metric $g$, we can decompose the vector space $V$ into a direct sum of orthogonal sum of orthogonal complements $V=A\oplus B$. In this case, the restricted complex structure $J_{A}$ is diagonalizable and its spectrum consists of complex conjugate pairs $\pm ic_{i}$ with $c_{i}\in[0,1]$. If $A$ is odd-dimensional, $J_{A}$ also has the eigenvalue [math]. Provided that we choose a decomposition into even-dimensional orthogonal complements with dimensions $2N_{A}$, we can bring $J_{A}$ and $g_{A}$ simultaneously in the following block diagonal form:

[TABLE]

Proof.

The restricted complex structure $J_{A}$ is itself anti-hermitian with respect to the metric $g_{A}$ because we have for $a_{1},a_{2}\in A$

[TABLE]

Therefore, $J_{A}$ is normal with respect to $g_{A}$. This ensures that $J_{A}$ is diagonalizable (with a complete set of eigenvectors) and so is $J_{A}^{2}$. Given a non-zero eigenvector $a\in A$ with $J_{A}^{2}\,a=\alpha\,a$, we have $\alpha\in[-1,0]$ due to the following argument. Let us compute

[TABLE]

This already implies that $\alpha\leq 0$. Moreover, we can also compute

[TABLE]

where we used in the last step that $A$ and $B$ are orthogonal which elimates crossing terms. This equation implies the inequality $g(a,a)\geq g(J_{A}a,J_{A}a)=|\alpha|g(a,a)$ leading to $\alpha\in[-1,0]$. From here, we find that $J_{A}$ must have eigenvalues appearing in complex conjugate pairs $\pm\mathrm{i}c$ with $c=\sqrt{-\alpha}$ meaning that also every eigenvalue $\alpha$ of $J_{A}^{2}$ appears with even multiplicity unless $\alpha=0$.

While the full $J$ satisfies $g(Jv,Jw)=g(v,w)$ for $v,w\in V$, the restricted $J_{A}$ does in general not satisfy $g_{A}(J_{A}a,J_{A}a^{\prime})=\omega_{A}(a,a^{\prime})$ for $a,a^{\prime}\in A$. However, we can use the relation $g(Jv,w)=-g(v,Jw)$ for $v,w\in V$ to derive the relation

[TABLE]

which is well-known to be the condition for $J_{A}$ to represent an orthogonal algebra element, i.e., $J_{A}\in\mathfrak{so}(2N_{A})$. An orthogonal algebra element is antisymmetric with respect to $g_{A}$ and has purely imaginary eigenvalues. It is well-known that we can always choose an orthonormal basis, such that the matrix representations of $J_{A}$ and $g_{A}$ are given by (401). ∎

Proposition 9**.**

In the setup of Proposition 8, if $J_{A}^{2}=-\mathbb{1}$ this implies $J_{BA}=0$, furthermore $J_{B}^{2}=-\mathbb{1}$ and $J_{AB}=0$.

Proof.

For any two vectors $a,a^{\prime}\in A$, we have

[TABLE]

which is shown using $g(v,Jw)=-g(Jv,w)$ and the fact that $A$ and $B$ are orthogonal subspaces:

[TABLE]

If $J_{A}^{2}=-\mathbb{1}$, then from (354) follows $J_{AB}J_{BA}=-\mathbb{1}-J_{A}^{2}=0$. This however means that for all $a\in A$ we have

[TABLE]

Since $g$ is non-degenerate this implies $J_{BA}=0$. By Proposition 2, $J_{B}^{2}=-\mathbb{1}$ follows from $J_{A}^{2}=-\mathbb{1}$, then repeating the argument above shows $J_{AB}=0$. ∎

This result is central for computing the entanglement spectrum of fermionic systems. In particular, the fact that the magnitude of possible eigenvalues lies in the interval $[0,1]$ implies that the entanglement entropy of fermionic systems is bounded by $\log{2}$ per fermionic degree of freedom.

Proposition 10**.**

As before, let the complex structure $J$ and the metric $g$ on a vector space $V$ be compatible, and let $V=A\oplus B$ be a decomposition of $V$ into even-dimensional orthogonal complements. Without loss of generality, assume that the dimension of $A$ is lower than, or equal to the dimension of $B$. Then there exists an orthonormal basis with respect to which $J$ takes the form

[TABLE]

with $r_{i}\in[0,\frac{\pi}{4}]$ and the two-by-two matrices $\mathbb{A}_{2}$ and $\mathbb{S}_{2}$ defined in (386).

Proof.

By Proposition 8 the upper left diagonal block, representing $J_{A}$, can be brought into the claimed form by choosing $\cos 2r_{i}=c_{i}$. Let $a^{(i)}_{1},\,a^{(i)}_{2}$ denote the pair of basis vectors spanning the space corresponding to $r_{i}$. Then, if $r_{i}>0$, define the vectors

[TABLE]

These vectors are orthonormal

[TABLE]

where we used $J_{AB}J_{BA}a=(-1+\cos^{2}2r_{i})a=-\sin^{2}(2r_{i})a$ which follows from (354). Similarly, one also shows that $g(b^{(i)}_{n},b^{(j)}_{m})=0$ if $i\neq j$. By definition, these vectors bring the lower left diagonal block, representing $J_{BA}$ into the claimed standard form. To see that the vectors also bring the diagonal blocks of $J_{B}$ into the claimed form, we calculate

[TABLE]

Hence, by extending the set of all these vectors with an orthonormal set of vectors which brings the eigenspace of $J_{B}^{2}$ with eigenvalue $-1$ into the standard form, we obtain an orthonormal basis of $B$ which brings the matrix representing $J$ into the claimed form.

Finally, this result implies the standard form (220) of the fermionic covariance matrix $\Omega$ by applying $\Omega=JG$ with $G\equiv\mathbb{1}$. ∎

The propositions of these section imply both the normal form of the covariance matrix for fermions (211), as well as the fermionic partner mode. In fact, the matrix $\mathbf{\Omega}$ representing the fermionic covariance matrix $\Omega^{ab}$ and the matrix $\mathbf{J}$ representing the complex structure $J^{a}{}_{B}$ are in fact identical. This is due to $\mathbf{\Omega}=\mathbf{J}\mathbf{G}$ where the metric $G$, which is the inverse of the metric $g$ used in this subsection, is represented by the identity matrix $\mathbf{G}=\mathbb{1}$. The partner formula for fermions (226) follows by the same argument as discussed above for bosons.

Appendix B Spectra of restricted Hamiltonian

Here, we show how the spectrum of the restriction of a Hamiltonian to a subspace of partner modes, is related to the spectrum of the unrestricted Hamiltonian operator, as a consequence of the Courant–Fischer–Weyl min-max principle.

Proposition 11**.**

Let $\hat{H}$ be a quadratic Hamiltonian on a (bosonic or fermionic) system of $N$ modes, and denote the excitation energies of $\hat{H}$ by $0<\omega_{1}\leq\omega_{2}\leq\dots\leq\omega_{N}$. Assuming that $\hat{H}_{A}$ is the restriction of $\hat{H}$ onto a subsystem $A$ spanned by $N_{A}<N$ modes such that the ground state $\ket{0}$ of $\hat{H}$ is a product state

[TABLE]

between $A$ and its complement $B$, the excitation energies $\epsilon_{1}\leq\dots\leq\epsilon_{N_{A}}$ of $\hat{H}_{A}$ fulfill

[TABLE]

Proof.

We note that the excitation energies of a quadratic Hamiltonian $\hat{H}=\frac{1}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}+f_{a}\hat{\xi}^{a}$ are given by the eigenvalues of the linear map $L^{a}{}_{c}=G^{ab}h_{bc}$ obtained by composing the covariance matrix $G^{ab}$ of the Hamiltonian’s ground state with the Hamiltonian bilinear form $h_{ab}$. This is easily seen by observing that with respect to a basis of energy eigenmodes, where $G^{ab}\equiv\mathbb{1}$ is represented by the identity matrix, both $L$ and $h$ are represented by the same diagonal matrix.

In general, in a basis which is adapted to the bipartition of the $N$ modes into subsystems $A$ and $B$, the map $L$ may not be represented by a diagonal matrix. However, since we assume the ground state of $\hat{H}$ to be a product state between $A$ and $B$, the covariance matrix takes block diagonal form, and we find that $L$ is represented by the matrix

[TABLE]

The min-max principle now warrants that the eigenvalues $\omega_{1}\leq\dots\leq\omega_{N}$ of $\mathbf{L}$ (with double multiplicity) bound the eigenvalues of $\mathbf{G}_{A}\mathbf{h}_{A}$, which we denote by $\alpha_{1}\leq\dots\leq\alpha_{2N_{A}}$, as

[TABLE]

The $\alpha_{i}$ with double multiplicity represent the excitation energies of $\hat{H}_{A}$. We can see this from the ground state of $\hat{H}$, that is a product state between $A$ and $B$ such that $\ket{\psi}_{A}$ must be the ground state of $\hat{H}_{A}$ and its covariance matrix is represented by the matrix $\mathbf{G}_{A}$. Therefore, the spectrum of $L_{A}=\mathbf{G}_{A}\mathbf{h}_{A}$ yields the excitation energies of $\hat{H}_{A}$.

The equivalent statement for a fermionic Hamiltonian $\hat{H}=\frac{\mathrm{i}}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}$ follows from the same line of argument, but the linear map $L$ is now defined by $L^{a}{}_{c}=-\Omega^{ab}h_{bc}$, where $\Omega^{ab}$ represents the covariance matrix of the ground state of the fermionic Hamiltonian $\hat{H}=\frac{\mathrm{i}}{2}h_{ab}\hat{\xi}^{a}\hat{\xi}^{b}$. ∎