# Mode Entanglement in Fermionic and Bosonic Harmonium

Jan Ole Ernst  Felix Tennie<EMAIL_ADDRESS>University of
Oxford, Department of Physics, Clarendon Laboratory, Parks Road, Oxford OX13PU

(August 28, 2024)

###### Abstract

Mode entanglement in many-body quantum systems is an active area of research.
It provides crucial insight into the suitability of many-body systems for
quantum information processing tasks. Local super-selection rules must be
taken into account when assessing the amount of physically accessible
entanglement. This requires amending well-established entanglement measures by
incorporating local parity and local particle number constraints. In this
paper, we report on mode entanglement present in the analytically solvable
system of $N$-Harmonium. To the knowledge of the authors, this is the first
analytic study of the physically accessible mode and mode-mode entanglement of
an interacting many-body system in a continuous state space. We find that
super-selection rules dramatically reduce the amount of physically accessible
entanglement, which vanishes entirely in some cases. Our results strongly
suggest the need to re-evaluate intra and inter-mode entanglement in other
fermionic and bosonic systems.

###### pacs:

Valid PACS appear here

##  I Introduction

Entanglement is “the characteristic trait of quantum physics” [1] and despite
recent advances still remains somewhat elusive. It provides insights into the
behaviour of physical systems and forms a cornerstone for the investigation of
many-body quantum physics [2]. Therefore, describing and quantifying
entanglement in many-body systems is necessary for understanding and
developing adequate models. Furthermore, entanglement is a key resource for
quantum information processing, particularly for the realisation of quantum
computers [3]. It is widely conjectured that systems exhibiting only a limited
amount of entanglement could be emulated by existing noisy quantum simulation
and computation technologies [4].

Fermionic mode entanglement plays a crucial role in quantum chemistry [5, 6].
Moreover, following advances in single electron manipulation [7], fermionic
mode entanglement offers a potential resource for quantum teleportation
schemes and information processing tasks [8, 9, 10, 11]. Mode entanglement is
also highly relevant in systems of cold bosonic atoms [12], and in a
continuous variable setting, error correction codes show promise for the
realisation of fault-tolerant quantum computing [13].

The assessment and quantification of entanglement in systems of
indistinguishable particles, i.e. fermions and bosons, has been the subject of
many controversial debates [14, 15, 16, 17]. While some states, such as the
two-particle state $(|01\rangle\pm|10\rangle)$, appear to contain entanglement
in a ‘first-quantisation’ framework, it turns out that this form of
entanglement is generally not accessible by physical measurements.

Recently, Ding et al. [5] worked out an approach that consistently generalises
well-known entanglement concepts and measures to systems of indistinguishable
particles. They then evaluated these new adapted measures in analytically
solvable low dimensional fermionic model systems, and provided a numerical
study of molecular systems by truncating their infinite dimensional active
spaces.

Despite these numerical results, no analytic study on the generalised
entanglement measures introduced by Ding et al. for interacting many-body
systems with a continuous state space has been conducted so far. In most cases
it is not possible to determine an analytic form of the energy eigenstates of
an interacting many-body Hamiltonian. However, the system of Harmonium is a
notable exception. Harmonium is a model system of harmonically interacting
particles confined by a harmonic trapping potential. To leading order, it
resembles the Hamiltonians of systems such as quantum dots [18, 19, 20, 21,
22, 23, 24]). In this paper we investigate the amount of physically accessible
entanglement in Harmonium. In particular, we evaluate mode and mode-mode
entanglement of the Harmonium ground state. Remarkably, the analytic form of
the ground state permits an analytic evaluation of the generalised
entanglement measures. We find that the amount of physically accessible
entanglement is significantly smaller than the apparent amount of so-called
‘fluffy-bunny’ entanglement [25].

This paper is organised as follows: in section II we introduce correlation and
entanglement measures and review the concepts of particle and mode-reduced
density operators. We also describe the model system of Harmonium and discuss
the role and status of (local) super-selection rules. In section III we
present our results of mode and mode-mode entanglement in fermionic and
bosonic Harmonium. In section IV we provide a discussion and outlook.

##  II Background

###  II.1 Correlation and Entanglement Measures

The state of any quantum system can be described by a positive semi-definite
density operator $\rho=\sum_{i}p_{i}\ket{\psi_{i}}\bra{\psi_{i}}$ with unit
trace $\sum_{i}p_{i}=1$ [26]. Physical observables are associated with
operators $\hat{O}$ which form an algebra $\pazocal{A}(\pazocal{H})$. The
expectation value of any observable $\hat{O}$ with respect to $\rho$ reads:

$\Braket{\hat{O}}_{\rho}=tr(\rho\hat{O}).$ (1)

Provided that the total Hilbert space can be expressed as a product of two
Hilbert spaces $\pazocal{H}=\pazocal{H}_{A}\otimes\pazocal{H}_{B}$, and that
the set of physical observables form an algebra which correspondingly
decomposes as
$\pazocal{A}(\pazocal{H})=\pazocal{A}(\pazocal{H}_{A})\otimes\pazocal{A}(\pazocal{H}_{B})$,
the reduced density operator $\rho_{A}=tr_{B}(\rho)$ of a particular subsystem
$A$ is uniquely determined by the following condition on the observables
$\hat{O}_{A}$ [26]:

$\displaystyle\Braket{\hat{O}_{A}}_{\rho_{A}}=\Braket{\hat{O}_{A}\otimes\mathbb{1}_{B}}_{\rho},\forall\hat{O}_{A}\in\pazocal{A}(\pazocal{H}_{A}).$
(2)

Within the set of quantum states, one distinguishes uncorrelated states,
correlated states, separable states and entangled states. The notion of
correlation and entanglement crucially relies on the product structure of the
algebra of observables. For an uncorrelated state, all expectation values
factorise into expectation values of local measurements on the subsystems:

######  Definition II.1 (Correlated states)

A state represented by density operator $\rho$ is uncorrelated iff:
$\Braket{\hat{O}_{A}\otimes\hat{O}_{B}}_{\rho}=\Braket{\hat{O}_{A}}_{\rho_{A}}\Braket{\hat{O}_{B}}_{\rho_{B}},\hat{O}_{A/B}\in\pazocal{A}(\pazocal{H}_{A/B})$.
The set of uncorrelated states is denoted by ${S_{unc}}$. A state is
correlated iff its density operator $\rho\not\in S_{unc}$.

The notion of separable states allows one to distinguish classical
correlations and quantum entanglement. The set of all classical mixtures of
uncorrelated states forms the set of separable states. States which are not
separable are called entangled. More formally:

######  Definition II.2 (Separable and entangled states)

A bipartite state $\rho$ is separable iff it can be represented as:
$\rho=\sum_{i}p_{i}(\rho_{A}^{i}\otimes\rho_{B}^{i})$ for a probability
distribution ${p_{i}}$ and density operators $\rho_{A/B}\in\pazocal{H}_{A/B}$.
Denote the set of these states by ${S_{sep}}$. A state $\rho$ is called
entangled iff $\rho\not\in S_{sep}$.

Pure, uncorrelated states are necessarily separable. Geometrically, the set of
separable states is the convex hull of the set ${S_{unc}}$ as illustrated in
Fig. 1.111The separability of a state can be assessed by a number of criteria.
In this paper, we shall use the PPT criterion [27] which is a necessary and
sufficient condition for separability in two-dimensional subsystems and
provides a lower bound for systems in larger state spaces.

Classical correlations are quantified by the Shannon entropy. Its quantum
analogue is the von Neumann entropy:

$S(\rho):=-tr[\rho\ln(\rho)]=-\sum_{i}\lambda_{i}\ln(\lambda_{i}),$ (3)

where $\lambda_{i}$ are the eigenvalues of $\rho$ [26]. The quantum relative
entropy is a measure which captures our ability to distinguish two quantum
states $\rho$ and $\sigma$:

$S(\rho\|\sigma):=tr(\rho\;ln(\rho)-\rho\;ln(\sigma)).$ (4)

The quantum mutual information $I(\rho)$ is a widely accepted measure of the
total correlation [28, 5] between two subsystems $A$ and $B$:

$I(\rho)\equiv
I(A:B):=S(\rho_{A})+S(\rho_{B})-S(\rho)=S(\rho\|\rho_{A}\otimes\rho_{B}).$ (5)

$\rho$$\sigma*$$E$$I$$C$${{S_{u}}_{n}}_{c}$Entangled StatesSeparable
StatesUncorrelatedStates Figure 1: Visual representation of sets of quantum
states and correlation measures [5].

Geometrically, the quantum mutual information describes the distance to the
closest uncorrelated state [5] $I(\rho)=min_{\sigma\in
S_{unc}}S(\rho\|\sigma)$. We consider the relative entropy of entanglement,
$E$ [29] as our measure of entanglement. $E$ is defined as the quantum
relative entropy of the state $\rho$ and the closest separable state
$\sigma*\in{S_{sep}}$:

$E(\rho)=\min_{\sigma\in S_{sep}}S(\rho\|\sigma)=S(\rho\|\sigma*),$ (6)

$E$ is strictly positive if and only if $\rho$ is non-separable, i.e.
entangled [29]. Geometrically, the relative entropy of entanglement describes
the distance to the closest separable state $\sigma*$, cf. the dotted red line
in Fig. 1. The relative entropy of entanglement for a bipartite system in a
pure state simplifies to $E(\rho)=S(\rho_{A/B})$.

The classical correlation of a state $\rho$ is defined as the quantum relative
entropy of the separable state $\sigma*$ and the closest uncorrelated state
$\rho_{A}\otimes\rho_{B}$ [28]:

$C(\rho)=S(\sigma*\|\rho_{A}\otimes\rho_{B}).$ (7)

The classical correlation can be visualised by the geometrical distance
between $\sigma*$ and the closest uncorrelated state. In general, the amount
of entanglement and classical correlation does not sum up to the total
correlation $(I\neq E+C)$ as mixed quantum states give rise to correlation
described by the quantum discord [30].

###  II.2 $N$-Harmonium

The Hamiltonian of the $N$-Harmonium model system reads:

$H=\sum_{i=1}^{N}(\frac{{{\vec{p}}_{i}}^{\,2}}{2m}+\frac{m}{2}\vec{x}_{i}^{T}\Omega\vec{x}_{i})\\\
+\frac{K}{2}\sum_{1\leq i<j\leq N}(\vec{x}_{i}-\vec{x}_{j})^{2},$ (8)

where $\Omega=diag(\omega^{1},\cdots,\omega^{n})$ parametrises the harmonic
trap strength and $K$ parametrises the particle-particle interactions;
$(\vec{x}_{1},\ldots,\vec{x}_{N})$ and $(\vec{p}_{1},\ldots,\vec{p}_{N})$
denote particle positions and their respective momenta. Coupling can be
attractive $K<0$ or repulsive $K>0$. One requires $\frac{-\Omega^{2}}{mN}<K$
for the existence of bound states. The Hamiltonian can be separated and
diagonalised by a suitable coordinate transformation (cf. Ref. [31] for
further details).

The bosonic ground state has the profile of a Gaussian and reads:

$\Psi^{(b)}(\vec{x}_{1},\dots,\vec{x}_{N})=\mathscr{N}e^{\vec{X}^{T}\bm{B}\vec{X}-\sum_{i=1}^{N}(\vec{x}^{T}_{i}\bm{A}\vec{x}_{i})}$
(9)

with $\vec{X}\equiv\frac{1}{N}(\vec{x}_{1}+\ldots\vec{x}_{N})$ representing
the centre-of-mass vector, and constants
$\bm{A}\equiv\mbox{diag}(A^{(1)},\ldots,A^{(n)})$, where
$A^{\alpha}\equiv\frac{1}{2(\tilde{l}^{(\alpha)})^{2}}$ and
$\bm{B}\equiv\mbox{diag}(B^{(1)},\ldots,B^{(n)})$, where
$B^{(\alpha)}\equiv\frac{N}{2}\left(\frac{1}{(\tilde{l}^{(\alpha)})^{2}}-\frac{1}{(l^{(\alpha)})^{2}}\right)$,
$\tilde{\omega}^{(\alpha)}\equiv\sqrt{(\omega^{(\alpha)})^{2}+\frac{NK}{m}}$,
$\tilde{l}^{(\alpha)}\equiv\sqrt{\frac{\hbar}{m\tilde{\omega}^{(\alpha)}}}$,
$l^{(\alpha)}\equiv\sqrt{\frac{\hbar}{m\omega^{(\alpha)}}}$, and $\mathscr{N}$
represents a normalisation constant.

The fermionic ground state is found to equate to a product of a Slater
determinant of Hermite functions and an exponential factor of the centre-of-
mass coordinate. This represents a product of a determinant of polynomials
(the Vandermonde determinant, which is anti-symmetric under the exchange of
particles) and the bosonic ground state (cf. Ref. [32]). More explicitly, the
$N$-fermion ground state $\Psi^{(f)}$ of the Harmonium model is given by:

$\Psi^{(f)}(\vec{x}_{1},\dots,\vec{x}_{N})=\mathscr{N}\,\cdot\left|\begin{array}[]{c
c
c}\phi_{\bm{\mu}_{1}}^{(\mathbf{\tilde{\mathbf{l}}})}(\vec{x}_{1})&\cdots&\phi_{\bm{\mu}_{1}}^{(\mathbf{\tilde{\mathbf{l}}})}(\vec{x}_{N})\\\
\vdots&&\vdots\\\
\phi_{\bm{\mu}_{N}}^{(\mathbf{\tilde{\mathbf{l}}})}(\vec{x}_{1})&\cdots&\phi_{\bm{\mu}_{N}}^{(\mathbf{\tilde{\mathbf{l}}})}(\vec{x}_{N})\end{array}\right|\cdot
e^{\vec{X}^{T}\bm{B}\vec{X}}\,.$ (10)

Here, the $\phi_{\bm{\mu}_{k}}^{(\mathbf{\tilde{\mathbf{l}}})}$ represent
products of Hermite functions of the relevant spatial dimensions:
$\phi_{\bm{\mu}}^{(\mathbf{\tilde{\mathbf{l}})}}(\vec{x}_{i})=\prod_{\alpha=1}^{n}\phi_{\mu^{\alpha}}^{(\tilde{l}^{\alpha})}(x^{\alpha}_{i})$,
where $\phi_{\mu^{\alpha}}^{(\tilde{l}^{\alpha})}(x^{\alpha}_{i})$ are Hermite
functions of degree $\mu^{\alpha}$ and natural length scale
$\tilde{l}^{\alpha}$.

The quantum number vectors $\bm{\mu_{1}},\cdots,\bm{\mu_{N}}$ are determined
by a form of an exclusion principle while minimising the total energy. In this
instance, one should consider the spin degrees of freedom to be effectively
frozen out by imposing a strong external magnetic field which aligns them.
Consequently, in spin-aligned Harmonium no pair of $\bm{\mu_{i}}$ can be
identical. Note that $l^{(\alpha)}$ is an intrinsic length scale and one can
set $l^{(1)}=1$ without loss of generality. The lengths $\tilde{l}^{(\alpha)}$
are then determined by the dimensionless coupling constants:

$\kappa^{(\alpha)}=\left(\frac{l^{(\alpha)}}{\tilde{l}^{(\alpha)}}\right)^{4}-1.$
(11)

$\kappa^{(\alpha)}$ is bounded, as one requires $\kappa^{(\alpha)}>-1$ for
bound state solutions. We can thus express the bosonic and fermionic ground
states as a function of the dimensionless constant $\kappa^{(\alpha})$.

###  II.3 Tracing out particles and modes

In this section, we outline the concepts of tracing out particles and modes;
further details can be found in Ref. [33]. It should be noted that mode
reduced density operators fundamentally differ from particle reduced density
operators.

For a state of fixed particle number, such as the $N$-Harmonium ground state,
the $p$-particle reduced density operator reads:

$\displaystyle\rho^{(p)}:=$
$\displaystyle\int_{\mathscr{R}}\psi(x_{1},\dots,x_{p},x_{p+1},\dots,x_{N})$
(12) $\displaystyle\cdot$
$\displaystyle\psi^{*}(x_{1}^{\prime},\dots,x_{p}^{\prime},x_{p+1},\dots,x_{N})dx_{p+1}\dots
dx_{N}.$

Particle reduced density operators are well defined in systems of both
distinguishable and indistinguishable particles. They correspond to
observables of fixed particle number $p$.

Observations on a reduced set of modes give rise to the notion of mode-reduced
density operators. In systems of many indistinguishable particles, the total
Hilbert space of (anti-)symmetrised states does not have an apparent product
structure. Let $\\{{\ket{\phi_{i}^{1}}}_{i=1}^{\pazocal{D}}\\}$ form a basis
of modes of the single-particle Hilbert space $\pazocal{H}^{(1)}$.
Observations confined to a subset of modes
$\\{{\ket{\phi_{i}^{(1)}}}_{i=1}^{\pazocal{D}^{A}}\\}$ are related to a
decomposition of the single-particle Hilbert space into an orthogonal sum:

$\pazocal{H}^{(1)}=\pazocal{H}^{A}\oplus\pazocal{H}^{B}=sp\\{{\ket{\phi_{i}^{(1)}}}_{i=1}^{\pazocal{D}^{A}}\\}\;\oplus\;sp\\{{\ket{\phi_{i}^{(1)}}}_{i=\pazocal{D}^{A}+1}^{\pazocal{D}}\\}.$
(13)

In the following, we shall specifically outline the steps to trace out modes
in systems of fermions and bosons. Creation and annihilation operators on the
subspaces $\pazocal{H}^{A}$ and $\pazocal{H}^{B}$ will allow us to construct a
product Hilbert space that is isomorphic to the total Hilbert space.

####  II.3.1 Fermionic Mode-Reduced Density Operators

Fermionic annihilation and creation operators satisfy canonical anti-
commutation rules [34]:

$\displaystyle f_{i}f_{j}^{{\dagger}}+f_{j}^{{\dagger}}f_{i}=$
$\displaystyle\delta_{ij}\mathds{1},$ (14) $\displaystyle
f_{i}f_{j}+f_{j}f_{i}=$ $\displaystyle
f_{i}^{{\dagger}}f_{j}^{{\dagger}}+f_{j}^{{\dagger}}f_{i}^{{\dagger}}=0.$

They can be used to define occupation number states, which form a basis of the
Fock space $\pazocal{F}$ associated with
$\pazocal{H}^{(1)}=sp\\{{\ket{\phi_{i}^{(1)}}}_{i=1}^{\pazocal{D}}\\}$,

$\displaystyle\ket{n_{1},\dots,n_{\pazocal{D}}}:=(f_{1}^{{\dagger}})^{\left(n_{1}\right)}\dots(f_{\pazocal{D}}^{{\dagger}})^{\left(n{D}\right)}\ket{\Omega}.$
(15)

Here, $\ket{\Omega}$ denotes the vacuum state and $n_{i}=\\{0,1\\}$. Fermionic
creation and annihilation operators associated with the modes of $H^{A}$ and
$H^{B}$ give rise to Fock spaces $\pazocal{F}_{A}$ and $\pazocal{F}_{B}$. A
natural isomorphism between the product of these Fock spaces and the total
Fock space

$\pazocal{F}\cong\pazocal{F}_{A}\otimes\pazocal{F}_{B}$ (16)

can be constructed by mapping number states as follows:

$\ket{n_{1},\dots,n_{\pazocal{D}}}\mapsto\underbrace{\ket{n_{1},\dots,n_{\pazocal{D}_{A}}}}_{\pazocal{F}_{A}}\otimes\underbrace{\ket{n_{\pazocal{D}_{A+1}},\dots,n_{\pazocal{D}}}}_{\pazocal{F}_{B}}.$
(17)

Note that for the isomorphism to be well-defined, a strict order of the mode
labels must be adhered to. This becomes obvious when swapping two creation
operators which introduces an additional overall minus sign, cf. Eq.(14).

Given the product structure of the isomorphic Hilbert space
$\pazocal{F}_{A}\otimes\pazocal{F}_{B}$ it is possible to define reduced
density operators for specific subsets of modes. In this paper, we consider
mode reduced and mode-mode reduced density operators. Using the canonical
anti-commutation relations (14), matrix representations of the reduced density
operators can be formed out of expectation values.

The mode-reduced density matrix of mode $\ket{\phi_{m}}$ in the local basis
$\\{\ket{\Omega},f^{\dagger}_{m}\ket{\Omega}\\}$ reads [33]:

$\rho_{\\{m\\}}=\begin{pmatrix}\left\langle
f_{m}f_{m}^{{\dagger}}\right\rangle&\left\langle
f_{m}^{{\dagger}}\right\rangle\\\ \left\langle f_{m}\right\rangle&\left\langle
f_{m}^{{\dagger}}f_{m}\right\rangle\end{pmatrix}.$ (18)

The mode-mode reduced density matrix of modes $\ket{\phi_{i}}$ and
$\ket{\phi_{j}}$ in the local basis
$\\{\ket{\Omega},f^{\dagger}_{i}\ket{\Omega},f^{\dagger}_{j}\ket{\Omega},f^{\dagger}_{i}f^{\dagger}_{j}\ket{\Omega}\\}$
is given by:

$\displaystyle\;\;\;\rho_{\\{i,j\\}}=$ (19)
$\displaystyle\begin{pmatrix}\left\langle
f_{i}f_{i}^{{\dagger}}f_{j}f_{j}^{{\dagger}}\right\rangle&\left\langle
f_{i}f_{i}^{{\dagger}}f_{j}^{{\dagger}}\right\rangle&\left\langle
f_{i}^{{\dagger}}f_{j}f_{j}^{{\dagger}}\right\rangle&\left\langle
f_{i}^{{\dagger}}f_{j}^{{\dagger}}\right\rangle\\\ \left\langle
f_{i}f_{i}^{{\dagger}}f_{j}\right\rangle&\left\langle
f_{i}f_{i}^{{\dagger}}f_{j}^{{\dagger}}f_{j}\right\rangle&\left\langle
f_{i}^{{\dagger}}f_{j}\right\rangle&\left\langle
f_{i}^{{\dagger}}f_{j}^{{\dagger}}f_{j}\right\rangle\\\ \left\langle
f_{i}f_{j}f_{j}^{{\dagger}}\right\rangle&\left\langle
f_{j}^{{\dagger}}f_{i}\right\rangle&\left\langle
f_{i}^{{\dagger}}f_{i}f_{j}f_{j}^{{\dagger}}\right\rangle&\left\langle
f_{i}^{{\dagger}}f_{j}^{{\dagger}}f_{i}\right\rangle\\\ \left\langle
f_{j}f_{i}\right\rangle&\left\langle
f_{i}f_{j}^{{\dagger}}f_{j}\right\rangle&\left\langle
f_{i}^{{\dagger}}f_{j}f_{i}\right\rangle&\left\langle
f_{i}^{{\dagger}}f_{i}f_{j}^{{\dagger}}f_{j}\right\rangle\end{pmatrix}.$

It should be noted that, in general, fermionic modes comprise spatial and spin
degrees of freedom. Therefore, for a spin-$1/2$ particle, the reduced density
operator of a spatial mode $\phi$ is given by Eq. (19) with
$i\equiv\\{\phi,\uparrow\\}$ and $j\equiv\\{\phi,\downarrow\\}$. The reduced
density operator of two spatial modes becomes a $16\times 16$ matrix. In this
paper, we explicitly compute both types of reduced density operators for the
ground state of Harmonium. Further details will be presented in Sec. III and
the evaluation of the matrix elements for the Harmonium ground state is
outlined in App. B.

####  II.3.2 Bosonic Mode-Reduced Density Operators

Reduced density operators for subsets of modes can also be defined in bosonic
systems. Bosons satisfy the well-known canonical commutation relations [34]:

$\displaystyle b_{i}b_{j}^{{\dagger}}-b_{j}^{{\dagger}}b_{i}=$
$\displaystyle\delta_{ij}\mathds{1},$ (20) $\displaystyle
b_{i}b_{j}-b_{j}b_{i}=$ $\displaystyle
b_{i}^{{\dagger}}b_{j}^{{\dagger}}-b_{j}^{{\dagger}}b_{i}^{{\dagger}}=0.$

These can be used to define an (orthonormal) number-state basis of the bosonic
Fock-space $\pazocal{B}$:

$\ket{n_{1},\dots,n_{\pazocal{D}}}_{b}=\frac{(b_{i}^{{\dagger}})^{\left(n_{1}\right)}}{\sqrt{n_{1}!}}\dots\frac{(b_{\pazocal{D}}^{{\dagger}})^{\left(n{D}\right)}}{\sqrt{n{D}!}}\ket{\Omega},$
(21)

where $n_{i}=\\{0,\cdots,N\\}$ .

Since the occupation number of bosonic modes is not constrained by the Pauli
exclusion principle, the general expression of mode-reduced density operators
in terms of expectation values of creation and annhilation operators is more
complex.

For instance, the matrix elements of the mode-reduced density operator of mode
$\ket{\phi_{m}}$ in the local basis
$\\{\ket{\Omega},b_{m}^{\dagger}\ket{\Omega},\frac{(b_{m}^{\dagger})^{2}}{\sqrt{2!}}\ket{\Omega},\ldots\\}$
read:

$\displaystyle\left(\rho_{\\{m\\}}\right)_{(k,l)}=\frac{1}{\sqrt{k!\;l!}}\Bigg{\langle}(b^{\dagger}_{m})^{k}\left(\prod_{i=1}^{N_{max}}(1-b^{\dagger}_{m}b_{m}/i)\right)(b_{m})^{l}\Bigg{\rangle}$
(22)

Here, $N_{max}$ represents the highest particle number sector in which the
support of the global bosonic state does not vanish. In particular, for pure
states of fixed particle number such as the $N$-Harmonium ground state,
$N_{max}=N$. Also, note that for finite $N$, $\rho_{\\{m\\}}$ is a
$(N+1)\times(N+1)$ matrix.

Similarly, bosonic mode-mode reduced density operators $\rho_{\\{i,j\\}}$ can
be expressed as expectation values of creation and annihilation operators.
However, the size of $\rho_{\\{i,j\\}}$ scales quadratically with the particle
number $N$.

Unlike particle-reduced density operators, mode-reduced density operators can
be mixed states of various particle numbers, as evident by Eqs. (18), (19) and
(22). When assessing entanglement and correlations of mode-reduced density
operators, it becomes necessary to take into account that some parts of the
density operators are not accessible by physical observations. More
specifically, super-selection rules prohibit certain (projective)
measurements, which in turn constrains the algebra of observables (2). This
shall be discussed in the following section.

###  II.4 Super-selection Rules

The study of entanglement in non-relativistic systems of indistiguishable
particles has motivated much discussion with regards to determining the part
of the physical state space, which is in principle experimentally accessible.

It is widely accepted that certain fermionic and bosonic states are un-
physical. More specifically, it appears impossible to form coherent
superpositions of states possessing different values of a conserved quantity
such as the (global) particle number. For instance, a coherent superposition
of a one-electron and two-electron state
$\ket{\psi}=\alpha\ket{e^{-}}+\beta\ket{e^{-},e^{-}}$ cannot be realised in
laboratories. This is an example of states that are prohibited by the global
particle number super-selection rule (g)N-SSR: for certain types of particles
(e.g. massive fermions), isolated systems have a fixed and conserved total
number of particles. Mathematically, this amounts to the particle number
operator commuting with the density operator of the system,
$[\hat{N},\rho]=0$. For systems which are partitioned into subsets $A$ and $B$
with respect to a particular mode, this further implies that states must be of
the form [15]:

$\ket{\Phi_{N}}=\sum_{k}\alpha_{k,N}\ket{k}_{A}\otimes\ket{N-k}_{B},$ (23)

where $k\in[0,1]$ for fermions, whereas $k\in[0,N]$ for bosons. In this paper,
we consider ground states of bosonic and fermionic Harmonium which have a
fixed particle number and are subject to the (g)N-SSR. In addition to the
(g)N-SSR, local parity and particle numbers with regards to a particular mode
sub-system, further constrain the physically realisable density operators.
This will be discussed subsequently for both fermions and bosons.

####  II.4.1 Fermions

Within the realm of experimental realisation, two constraints relevant for the
construction of mode reduced density operators are known for fermionic
systems: the local parity-super-selection rule (P-SSR) [35] and the local
particle number-super-selection rule (N-SSR) [36].

P-SSR precludes a superposition of even and odd-number fermion states. Under
local conservation of parity $q$ (P-SSR), observables $\hat{O}_{A/B}$ must
satisfy:

$\hat{O}_{A/B}^{P}=\sum_{q}P_{q}\hat{O}_{A/B}P_{q},$ (24)

where the $P_{q}$ represent projectors onto subspaces of definite parity
number. Therefore, the observables $\hat{O}_{A/B}^{P}$ are block-diagonal and
the algebra of observables decomposes into an even an odd parity sector
$\pazocal{A}_{ee}+\pazocal{A}_{oo}$. Note that a violation of local parity
super-selection would have wide-ranging consequences such as super-luminal
signalling [37, 5, 10] and is therefore considered impossible.

Local N-SSR is an even more stringent constraint and it requires the local
observables $\hat{O}_{A/B}^{N}$ to be block-diagonal in the particular
fermion-N-particle number sectors:

$\hat{O}_{A/B}^{N}=\sum_{N\geq 0}P_{N}\hat{O}_{A/B}P_{N},$ (25)

where the $P_{N}$ denote the projectors onto the local N-fermion subspace.

Local number super-selection has been the subject of many discussions. Though
it cannot be related to fundamental laws of nature such as no-super-luminal-
signalling, it applies to most practical scenarios. Recently, works have
suggested that the status of N-SSR is challenged by the existence of pure
states with indefinite particle numbers (e.g. Majorana fermions, cf. Ref.
[38]). One can however say, that the creation of useful coherent
superpositions of states with different numbers of fermions is incredibly
challenging. Within the realm of current experimentation, N-SSR tends to be
satisfied [10] as any results claiming to realise N-SSR violating states have
not been reproduced [39].

####  II.4.2 Bosons

Unlike fermions, bosonic systems do not exhibit a local parity super-selection
rule; however, a bosonic local N-SSR for the construction of physical density
operators has been scrutinuously discussed (cf. [15, 40]). In general, one
must discern systems of massive and massless bosons. For reasons similar to
the ones outlined in Sec. II.4.1, a local N-SSR should be applied to massive
bosons; in particular, the realisation of coherent superpositions of massive
bosons with different local particle numbers has proven to be extremely
challenging. Thus, the local observables, allowed under local N-SSR, are those
which are block-diagonal with respect to the bosonic N-particle number
sectors. This algebra is constructed analogously to its fermionic counterpart
defined in Eq. (25).

The introduction of a local N-SSR is motivated by the fact that there is no
known methodology to easily prepare massive bosonic states which violate local
N-SSR [15]. Note that the status of local N-SSR, particularly for massless
bosons, is still subject to much debate [40, 41, 42], hinging upon a choice of
phase reference. However, the status of N-SSR for massless bosons is not of
further relevance to the Harmonium model which describes massive particles.

###  II.5 Symmetries and Entanglement

In this section, we briefly state a theorem which becomes instrumental when
calculating entanglement. Genererally, due to the possibly high dimensionality
of the state space, it is very challenging to explicitly carry out the
minimisation in Eq. (6). However, for states which exhibit local unitary
symmetries, the computation simplifies significantly [43]:

###### Theorem II.1

If a density operator $\rho$ is invariant under local unitary operations U(g),
representing elements of a group $g\in G$, with generator $\mathscr{T}_{G}$,
s.t. ${U^{\dagger}(g)\rho U(g)=\rho}$, then the closest separable state with
respect to $\rho$, $\sigma^{*}$, possesses the same symmetries as $\rho$, such
that $\mathscr{T}_{G}(\sigma*)=\sigma*$.

Without proving the statement which is done elsewhere (cf. [44]), we merely
reiterate that the local unitaries are not entanglement generating. In this
paper, we focus on the symmetries of local parity and particle number, as
defined in Eqs. (24) and (25). In some cases, these symmetries even permit an
analytical evaluation of the relative entropy of entanglement.

##  III Results

We demonstrate the impact of particle number and parity constraints on
entanglement in fermionic and bosonic Harmonium. The ground state of Harmonium
is a pure and fixed particle number state; it therefore obeys the global
N-SSR. Due to local super-selection rules, not all parts of mode-reduced
density operators are accessible by physical measurements (cf. Sec. II.4). As
a consequence, care must be taken when determining the amount of physically
meaningful entanglement. Our results show a striking impact of super-selection
rules and reveal the difference between ‘apparent’ and physically accessible
entanglement in Harmonium. Since the Harmonium model of harmonically
interacting and trapped particles is considered as a reasonable first
approximation to systems such as quantum dots (cf. Refs. [23, 24]), our
results demonstrate the need to re-evaluate the notion of entanglement in
various other continuous systems.

Hermite functions constitute a natural single-mode basis of the Harmonium
model. In what follows, we shall use bases constructed out of Hermite
functions $\phi_{m}^{(\tilde{l})}$ of natural length scale $\tilde{l}$ (cf.
App. A). Note that by Eq. (11), selecting a different arbitrary length scale
$\bar{l}$ would merely amount to ‘stretching’ or ‘quenching’ the $\kappa$ axes
in the correlation and entanglement measure plots below. Using the basis of
Hermite functions, Slater determinants (and permanents) will be abbreviated by
their respective mode number, e.g.

$f^{\dagger}_{0}f^{\dagger}_{1}f^{\dagger}_{2}\ket{\Omega}\equiv\ket{\phi_{0},\phi_{1},\phi_{2}}\equiv\ket{m_{0},m_{1},m_{2}}$
(26)

for ease of notation.

###  III.1 Fermions

In the presence of an external magnetic field, an additional Zeeman term
$\pazocal{H}_{B}=\frac{c_{S}}{\hbar}\sum_{i}^{N}\hat{S_{i}}\cdot\vec{B_{ext}}$
must be added to the Harmonium Hamiltonian [32], where $\hat{S_{i}}$
represents the spin angular momentum operator of the $i$-th particle, $c_{S}$
is a constant and $B_{ext}$ represents the external magnetic field. The
additional Zeeman term commutes with the $N$-Harmonium Hamiltonian and leads
to a splitting of the single-particle states to
$\bar{\phi}_{i}(x;\uparrow\\!\\!/\\!\\!\downarrow)=\phi_{i}(x)\otimes\ket{\uparrow\\!\\!/\\!\\!\downarrow}$.
In spin-ful Harmonium, the ground state can be determined analogously to Eq.
((10)) by replacing $\phi$ with $\bar{\phi}$. Note that the relative strength
and orientation of the magnetic field determines the ground-state
configuration.

In the regime of strong magnetic fields, i.e.
${\absolutevalue{B_{ext}}\gg\frac{\hbar}{c_{s}}\omega\sqrt{1+\kappa}}$, the
system’s spin degree of freedom will be quenched and all spins aligned. The
resulting ground state equals the ground state of spin-less Harmonium.

For smaller field strengths, the individual spins will only partially be
aligned and the configuration of the ground state is modified accordingly. For
instance, if
$\frac{\hbar}{c_{s}}\omega\sqrt{1+\kappa}<\absolutevalue{B_{ext}}<\frac{2\hbar}{c_{s}}\omega\sqrt{1+\kappa}$,
the ground-state eigenfunction for 3-Harmonium in one spatial dimension
($n=1$) reads:

$\pazocal{N}_{0}^{s}\begin{vmatrix}\phi_{0}^{\tilde{l}}(x_{1})\otimes\ket{\downarrow}&\phi_{0}^{\tilde{l}}(x_{2})\otimes\ket{\downarrow}&\phi_{0}^{\tilde{l}}(x_{3})\otimes\ket{\downarrow}\\\
\phi_{1}^{\tilde{l}}(x_{1})\otimes\ket{\downarrow}&\phi_{1}^{\tilde{l}}(x_{2})\otimes\ket{\downarrow}&\phi_{1}^{\tilde{l}}(x_{3})\otimes\ket{\downarrow}\\\
\phi_{0}^{\tilde{l}}(x_{1})\otimes\ket{\uparrow}&\phi_{0}^{\tilde{l}}(x_{3})\otimes\ket{\uparrow}&\phi_{0}^{\tilde{l}}(x_{3})\otimes\ket{\uparrow}\\\
\end{vmatrix}\cdot e^{\vec{X}^{T}\bm{B}\vec{X}}.$ (27)

####  III.1.1 Mode Entanglement

Figure 2: Correlation and entanglement of spatial modes $\phi_{0}$ (top),
$\phi_{2}$ (bottom) in one-dimensional $3$-fermion Harmonium as a function of
the relative particle-particle interaction strength $\kappa$. Notably,
physically accessible correlation ($I^{N/P}$) and entanglement ($E^{N/P}$) are
significantly smaller by comparison to the measures calculated without
invoking SSR’s ($E$ and $I$). Figure 3: Relative amount of physically
accessible entanglement $E^{P/N}/E$ of spatial modes $\phi_{0}$ to $\phi_{6}$
as a function of the coupling strength $\kappa$ in one-dimensional 4-fermion
Harmonium.

Figure 4: Relative amount of physically accessible entanglement $E^{P/N}/E$ of
spatial modes $\phi_{0}$ (top) and $\phi_{2}$ (bottom) as a function of the
coupling strength $\kappa$ for one-dimensional Harmonium of different fermion
numbers $N=3,4,5$.

In the presence of a weak external field in spin-ful Harmonium, the mode
reduced density operator of spatial mode $m$ becomes diagonal in the local
basis
$\\{\ket{\Omega},\ket{m_{\uparrow}},\ket{m_{\downarrow}},\ket{m_{\uparrow},m_{\downarrow}}\\}$:

$\rho_{\\{m_{\uparrow},m_{\downarrow}\\}}=\text{diag}(q_{0},q_{1},q_{2},q_{3})$
(28)

The diagonal representation arises from the spin symmetries and fixed total
particle number of the Harmonium ground state (cf. Eq. (19)). The structure of
the Schmidt-decomposed $N$-particle state reads:

$\displaystyle\ket{\psi}=$
$\displaystyle\sqrt{q_{0}}\ket{\Omega}\otimes\ket{\psi_{N}}+\sqrt{q_{1}}\ket{m_{\uparrow}}\otimes\ket{\psi_{N-1}}$
(29)
$\displaystyle+\sqrt{q_{2}}\ket{m_{\downarrow})}\otimes\ket{\psi_{N-1}^{\prime}}+\sqrt{q_{3}}\ket{m_{\uparrow},m_{\downarrow})}\otimes\ket{\psi_{N-2}}.$

In our analysis, we have determined the amplitudes $\sqrt{q_{i}}$
analytically. This has become possible since the integral expressions for
expectation values of fermionic ladder operators (cf. App. B) amount to
integrals of a Gaussian factor multiplied by a polynomial. The analytic forms
of the $q_{i}$ have then been used to derive an analytic expression for the
spatial-mode entanglement and mutual information. Utilising the techniques
introduced in Ref. [5], super-selection rules are taken into account by
removing the unphysical parts of the reduced density operators. For the
reader’s reference, we give an account of the explicit expressions for the
evaluation of correlation and entanglement measures in App. C. Note that the
respective measures are labeled as $I,E$ (without super-selection rules), and
$I^{N/P},E^{N/P}$ (including the N-SSR, and P-SSR, respectively). Though the
analytic form of the $q_{i}$ contains a large number of terms (cf. App. D for
the special case of $N=3$ fermions), it allows us to derive an analytic
expression for the entanglement measures in an infinite-dimensional
interacting many-body system. It should be stressed that this itsself is a
remarkable result.

Figure 5: Relative amount of physically accessible entanglement $E^{P/N}/E$ in
two-dimensional $4$-fermion Harmonium as a function of trapping frequency
detuning $\omega_{y}/\omega_{x}$ and coupling strength $\kappa$. In two
spatial dimensions, the basis of spatial modes consists of products of two
Hermite functions, i.e. $\phi_{(i,j)}=\phi_{(i)}\phi_{(j)}$. Note that
$\omega_{y}/\omega_{x}$ is bounded from above by $3$. For relative frequency
detunings exceeding this value, the system has a different ground state which
is effectively one-dimensional. Colour online.

The broad hierarchy and variation of the spatial mode entanglement and
correlation measures with $\kappa$ (cf. Eq. (11)) is qualitatively similar for
various modes as shown in Fig. 2. As an example, this figure displays
correlation and entanglement of modes $\phi_{0}$ and $\phi_{2}$ in one-
dimensional 3-fermion Harmonium. Unsuspectedly, all correlation and
entanglement measures assume a local maximum at a finite positive value of
$\kappa$ before monotonically decreasing with further growing interaction
strengths. This is a peculiar feature of the Harmonium system, insofar as one
would expect mode entanglement to grow monotonically with the particle-
particle interaction strength $\kappa$, bounded, only by the theoretical
maximal entanglement of $E_{max}=ln(4)$. Furthermore, one observes that the
amount of physically accessible entanglement decreases significantly once
SSR’s are taken into account; (in Fig. 2 compare $E$ and $I$ with $E^{P/N}$
and $I^{P/N}$). Henceforth, we shall therefore evaluate and display the
relative amount of physically accessible entanglement $E^{P}/E$ and
$E^{N}/E$.222At $\kappa=0$ the entanglement $E$ vanishes and the fractions are
not defined. Thus, we shall provide the graphs of $E^{P}/E$ and $E^{N}/E$ with
analytic continuation at zero interaction strength.

Fig. 3 shows a comparison of the relative amount of physically accessible
entanglement between different spatial modes in one-dimensional 4-fermion
Harmonium. With the exception of modes $\phi_{1}$ and $\phi_{2}$, the relative
amount of entanglement accessible under P-SSR approaches zero for small
interaction strengths. For large $\kappa\gtrsim 10^{3}$, $E^{P}/E$ is
hierarchically ordered with the mode number, which is not the case in the
intermediate regime $\kappa\approx 1\ldots 100$. Interestingly, the relative
amount of entanglement of mode $\phi_{1}$ accessible under P-SSR remains
approximately constant, while for $\phi_{2}$ it assumes a maximum around small
interaction strengths. A slightly different behaviour is observed for the
relative amount of entanglement accessible under N-SSR. For all modes, a local
minimum is located around vanishing $\kappa$. Remarkably, the magnitudes of
the maxima at large $\kappa$ grow with increasing mode number.333One might
speculate that this result is a reminiscence of the impact of SSRs becoming
negligible in the many-copy limit, cf. Ref. [25]. For instance, $E^{N}$ of
mode $\phi_{6}$ at around $\kappa\approx 10^{4}$ accounts for approximately
$30\%$ of the entanglement $E$, the highest value of all modes under
consideration. It should be noted that the relative amount of entanglement
makes no statement on its absolute value, which decreases significantly for
higher mode numbers. For the sake of completeness, additional plots of the
entanglement measures for different particle numbers are shown in App. E.

The variation of spatial mode entanglement in one-dimensional Harmonium with
the particle number is shown in Fig. 4. For mode $\phi_{0}$, the relative
amount of entanglement under P-SSR and N-SSR assumes a minimum around
vanishing $\kappa$. Throughout, $E^{N}$ is smaller than $E^{P}$ and only
amounts to less than approximately $10\%$ of the total entanglement. For mode
$\phi_{2}$, the relative amount of entanglement assumes a maximum around
vanishing $\kappa$ for some particle numbers while being minimal for others.
Note that no clear hierarchy of the relative amount of accessible entanglement
with the particle number is observed.

In two-dimensional Harmonium, the form of the ground state is parametrised by
the detuning or quench of the trapping potential frequencies
$\omega_{y}/\omega_{x}$ and the relative interaction strength $\kappa$. In
Fig. 5 the relative amount of physically accessible entanglement is shown for
4-fermion Harmonium. Note that for detunings exceeding a certain threshold
value, e.g. $\omega_{y}/\omega_{x}=3$ for $N=4$ particles, the structure of
the two-dimensional Harmonium ground state changes and the system becomes
effectively one-dimensional. Generally, we observe that for fixed frequency
detuning the amount of accessible entanglement qualitatively behaves similarly
to the case of one-dimensional Harmonium and assumes a local maximum at finite
positive $\kappa$. For fixed interaction strength, increasing the quench of
the trap generally increases the relative amount of accessible entanglement.
As in the case of one-dimensional Harmonium, $E^{N}$ is always smaller than
$E^{P}$.

####  III.1.2 Mode-Mode Entanglement

In order to evaluate entanglement between two spatial fermionic modes, one has
to consider the two-mode reduced density operator. For spin-$1/2$ fermions,
this operator has a $16\times 16$ dimensional matrix representation, since
each spatial mode can be occupied by two anti-aligned spins. As in the case of
spatial-mode entanglement, we gauge the impact of local super-selection rules
by computing the entanglement measures with and without projecting the two-
mode reduced density operator onto the parity and number eigenstate basis. The
representation of such a two-mode-reduced density operator is well known in
quantum chemistry [45] and becomes block-diagonal in the number state basis
for pure states of fixed (global) particle number such as the Harmonium ground
state. Note that this implies that Harmonium cannot exhibit any mode-mode
correlation between number blocks with different particle numbers.

All coefficients of the two-mode reduced density operator can be computed
analytically. Using fermionic anti-commutation relations (14) the matrix
elements can be expressed in terms of particle-reduced density operators (cf.
App. F). Equipped with the analytic form of the two-mode reduced density
operator and its marginals, the evaluation of the entanglement measures
involves a parametrised optimisation to determine the closest separable state
(cf. Eq. (6)). Though generally challenging, this task simplifies
significantly once symmetries are taken into account (cf. Sec. II.5). The
optimisation is carried out using the PPT criterion [27, 46] as a necessary
condition for separability. When evaluating $E^{P/N}$, this criterion provides
an exact value of the entanglement since all matrix blocks are at most of size
$2\times 2$. When evaluating $E$, the criterion only provides a lower bound
[47]. However, since we are interested in gauging the relative amount of
physically accessible entanglement this lower bound is sufficient.
Furthermore, numerical computations reveal that this lower bound is almost
exact.

Mode-mode entanglement in $N$-fermion Harmonium is found to behave
qualitatively like in the specific case of $N=4$ fermions in one spatial
dimension, Fig. 6. We observe that the relative amount of physically
accessible entanglement $E^{P/N}$ is smaller by several orders of magnitude
than $E$. It has a minimum around vanishing coupling strength $\kappa$,
followed by an increase in the entanglement to a maximum around $\kappa\approx
10^{3}$ and a subsequent decrease for larger interaction strengths.
Interestingly, mode-mode entanglement between the two modes $\phi_{0}$ and
$\phi_{1}$ that form the ground state at $\kappa=0$ never exceeds $E\approx
10^{-5}$, even in the regime of strong couplings.

Figure 6: Mode-mode entanglement between mode $\phi_{0}$ and mode $\phi_{1}$
(top, black), and between mode $\phi_{0}$ and mode $\phi_{2}$ (bottom, grey)
in 4-fermion Harmonium as a function of $\kappa$. Note that the evaluation of
the entanglement measures involves a parametrised minimisation. Here, we have
used the PPT criterion [27] for separability which yields an exact value for
$E^{N/P}$ and provides a lower bound for $E$.

So far, we have assessed the amount of physically accessible entanglement in
systems of spin-ful fermions. The case of spin-less fermions can be considered
by applying an external magnetic field that aligns all spins. However, due to
the resulting symmetries no physically accessible mode entanglement is
observed while mode-mode entanglement is similar to the case of mode
entanglement for spin-ful fermions.

###  III.2 Bosons

We now consider the entanglement of modes in bosonic Harmonium with integer
spin $S=1$. As for fermions, the Zeeman term defined in Sec. III.1 is added to
the Hamiltonian which causes a splitting of the eigenstates with respect to
the magnetic quantum numbers $\sigma_{b}=\\{-1,0,1\\}$ (27). The relative
orientation of the magnetic field determines the ground-state spin
configuration of bosonic Harmonium. Thus, in the presence of an external
magnetic field, however weak it may be, it is energetically favourable to
occupy a particular spin mode. This scenario is physically realistic, because
all physical systems are subject to some external magnetic field, which will
perturb the energy levels of the system. For ease of notation, the spin degree
of freedom is omitted in what follows since all spins are assumed to be
aligned. This is equivalent to considering massive spin-less bosons.

####  III.2.1 Mode Entanglement

Figure 7: E of various spatial modes as a function of $\kappa$ in one-
dimensional 3-boson Harmonium. Figure 8: E of spatial modes $m_{0}$ (top) and
$m_{2}$ (bottom) as a function of $\kappa$ for one-dimensional Harmonium with
different boson numbers $N=3,4,5$.

Mode correlation and entanglement are encapsulated in the $(N+1)$ dimensional
mode reduced density operator (cf. Sec. II.3.2). In the basis
$\\{(b^{\dagger})^{i}\ket{\Omega}\\}^{N}_{i=0}$ this operator becomes diagonal
for the pure, fixed particle number ground state of Harmonium.

For example, for three bosons the mode-reduced density operator of mode $m$
can be parametrised as $\rho_{\\{m\\}}=diag(p_{0},p_{1},p_{2},p_{3})$.
Consequently, the Schmidt-decomposition of the three particle state can be
written as:

$\displaystyle\ket{\psi_{b}^{N=3}}=$
$\displaystyle\sqrt{p_{0}}\ket{\Omega}\otimes\ket{\psi_{3}^{\prime}}+\sqrt{p_{1}}\ket{m}\otimes\ket{\psi_{2}^{\prime}}$
(30)
$\displaystyle+\sqrt{p_{2}}\ket{m,m}\otimes\ket{\psi_{1}^{\prime}}+\sqrt{p_{3}}\ket{m,m,m}\otimes\ket{\Omega^{\prime}}.$

Similar to fermions, the matrix elements $p_{i}$ can be computed by evaluating
particle-reduced density operators. The related integrals contain a polynomial
multiplied by a Gaussian factor and therefore can be calculated analytically
(cf. App. B).

Due to the Schmidt-decomposition, the mutual information and entanglement
read:

$E(\rho_{b})=\frac{1}{2}I(\rho_{b})=-\sum_{k}p_{k}\ln(p_{k}).$ (31)

Once local N-SSR constraints are imposed, the Schmidt-decomposition of the
$N$-boson state also implies that the $N+1$ sectors with fixed local particle
number factorise individually. Therefore, we find that the amount of
entanglement accessible under N-SSR vanishes, $E^{N}=0$.

Generally, bosonic mode entanglement behaves qualitatively similar for
different particle numbers and modes. In Fig. 7, the entanglement of modes
$\phi_{0}$ to $\phi_{8}$ is displayed as a function of $\kappa$. While
vanishing for zero interaction strength, it attains a maximum at finite
$\kappa$ in the regime of strong couplings. We also observe a hierarchy of the
entanglement which strictly decreases with increasing mode number.

In Fig. 8, the influence of the boson number on mode entanglement is
demonstrated for modes $\phi_{0}$ and $\phi_{2}$ and $N=3,4,5$. No clear
hierarchy with the particle number is observed.

####  III.2.2 Mode-Mode Entanglement

Figure 9: Mode-mode entanglement between modes $\phi_{0}$ and $\phi_{1}$
(black, top), and between modes $\phi_{1}$ and $\phi_{2}$ (grey, bottom) in
3-boson Harmonium as a function of interaction strength $\kappa$.

To assess bosonic mode-mode entanglement, one has to determine the two-mode
reduced density operator. For pure fixed-particle number states, many of its
matrix elements in the particle number basis vanish. The example of two bosons
is outlined in Tab. 1. As above, one can compute the matrix elements of
$\rho_{\\{m,n\\}}$ by evaluating particle-reduced density operators in the
basis
$\\{(b_{m}^{\dagger})^{i}(b_{n}^{\dagger})^{j}\ket{\Omega}\\}^{2}_{i=0,j=0}$.

Note that once a local N-SSR is introduced, all local number sectors of
$\rho_{\\{m,n\\}}$ become separable in the fixed local-particle number basis.
In Tab. 1 this amounts to deleting the greyly marked matrix elements.
Therefore the entanglement accessible under local N-SSR vanishes, i.e.
$E^{N}=0$. Equipped with the analytic form of the two-mode reduced density
operator, one can calculate the entanglement between two modes by carrying out
a constrained minimisation in accordance with Eq. (6). Here, we again use the
PPT criterion which provides a lower bound for the entanglement [27, 46].

N | 0 | 1 | 2
---|---|---|---
Modes |  $\Omega$ | $m$ |  $n$ | $m,m$ | $m,n$ |  $n,n$
$\Omega$ |  |  |  |  |  |
$m$ |  |  |  |  |  |
$n$ |  |  |  |  |  |
$m,m$ |  |  |  |  |  |
$m,n$ |  |  |  |  |  |
$n,n$ |  |  |  |  |  |
Table 1: Elements of the boson 2-mode reduced density operator
$\rho_{\\{m,n\\}}$ with non-zero expectation values, without local N-SSR (grey
and black) and with local N-SSR (black).

Bosonic mode-mode entanglement qualitatively behaves similar for different
pairs of modes and particle numbers. Its key features resemble those of
fermion mode-mode entanglement. As an example, mode-mode entanglement in one-
dimensional $3$-boson Harmonium is shown in Fig. 9. Generally, the amount of
entanglement is small and decreases with increasing mode numbers.

##  IV Summary and Outlook

In this paper, we have investigated mode correlation and entanglement
contained in the ground state of bosonic and fermionic Harmonium. In
particular, we have evaluated the amount of entanglement that remains
accessible under physical constraints in form of super-selection rules. Our
main findings are that the amount of physically accessible entanglement is
considerably smaller than the ‘fluffy bunny’ entanglement which is calculated
without taking parity (fermions) and local number (fermions and bosons)
constraints into account. Remarkably, it was possible to derive analytic
expressions of the entanglement measures which is unusual given that Harmonium
is an interacting many-body system on a continuous state space. While
considering different modes, different numbers of particles and spatial
dimensions, we observed that mode and mode-mode entanglement in Harmonium
attain a maximum in the strong coupling regime while approaching zero again in
the limit of infinite interaction strength.

Our investigation of the physically accessible amount of mode-entanglement in
Harmonium provides a first analytic insight into mode-entanglement in
continuous variable systems. Furthermore, it is of interest to other
continuous variable systems whose Hamiltonians have Taylor expansions with
leading terms of quadratic order [23]. For such systems, our results provide
the basis for a first assessment of the physically meaningful entanglement.
This in particular, could be instrumental for gauging the suitability of such
systems for quantum information processing tasks, such as quantum
teleportation [9] or measurement based quantum computation [11]. Future
research should expand the above concepts to similar Hamiltonians by using
state perturbation methods. We remark that our findings indicate a general
necessity to re-evaluate the amount of physically meaningful entanglement in
systems of bosons and fermions.

###### Acknowledgements.

The authors gratefully thank Lexin Ding, Christian Schilling, Ben Yadin and
Axel Kuhn for useful comments and discussions. JOE acknowledges support from
the Studienstiftung des Deutschen Volkes. FT acknowledges support from the
Royal Society.

## References

  * Schrödinger [1935] E. Schrödinger, Discussion of probability relations between separated systems, Mathematical Proceedings of the Cambridge Philosophical Society 31, 555–563 (1935).
  * Amico  _et al._ [2008] L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Entanglement in many-body systems, Rev. Mod. Phys. 80, 517 (2008).
  * Deutsch [1995] D. Deutsch, Universality in quantum computation, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 449, 669 (1995).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Ding  _et al._ [2021] L. Ding, S. Mardazad, S. Das, S. Szalay, U. Schollwoeck, Z. Zimboras, and C. Schilling, Concept of orbital entanglement and correlation in quantum chemistry, Journal of Chemical Theory and Computation 17, 79 (2021), pMID: 33430597, https://doi.org/10.1021/acs.jctc.0c00559 .
  * Ding and Schilling [2020] L. Ding and C. Schilling, Correlation paradox of the dissociation limit: A quantum information perspective, Journal of Chemical Theory and Computation 16, 4159 (2020), pMID: 32433873, https://doi.org/10.1021/acs.jctc.0c00054 .
  * Hofer  _et al._ [2017] P. P. Hofer, D. Dasenbrook, and C. Flindt, On-demand entanglement generation using dynamic single-electron sources, physica status solidi (b) 254, 1600582 (2017).
  * Galler and Thunström [2021] A. Galler and P. Thunström, Orbital and electronic entanglement in quantum teleportation schemes, Phys. Rev. Research 3, 033120 (2021).
  * Olofsson  _et al._ [2020] E. Olofsson, P. Samuelsson, N. Brunner, and P. P. Potts, Quantum teleportation of single-electron states, Phys. Rev. B 101, 195403 (2020).
  * Debarba  _et al._ [2020] T. Debarba, F. Iemini, G. Giedke, and N. Friis, Teleporting quantum information encoded in fermionic modes, Phys. Rev. A 101, 052326 (2020).
  * Chiu  _et al._ [2013] Y.-J. Chiu, X. Chen, and I. L. Chuang, Fermionic measurement-based quantum computation, Phys. Rev. A 87, 012305 (2013).
  * Gross  _et al._ [2011] C. Gross, H. Strobel, E. Nicklas, T. Zibold, N. Bar-Gill, G. Kurizki, and M. K. Oberthaler, Atomic homodyne detection of continuous-variable entangled twin-atom states, Nature 480, 219 (2011).
  * Gottesman  _et al._ [2001] D. Gottesman, A. Kitaev, and J. Preskill, Encoding a qubit in an oscillator, Phys. Rev. A 64, 012310 (2001).
  * Benatti  _et al._ [2020] F. Benatti, R. Floreanini, F. Franchini, and U. Marzolino, Entanglement in indistinguishable particle systems, Physics Reports 878, 1 (2020).
  * Dalton  _et al._ [2017] B. J. Dalton, J. Goold, B. M. Garraway, and M. D. Reid, Quantum entanglement for systems of identical bosons: I. general features, Physica Scripta 92, 023004 (2017).
  * Morris  _et al._ [2020] B. Morris, B. Yadin, M. Fadel, T. Zibold, P. Treutlein, and G. Adesso, Entanglement between identical particles is a useful and consistent resource, Phys. Rev. X 10, 041012 (2020).
  * Cunha  _et al._ [2007] M. O. T. Cunha, J. A. Dunningham, and V. Vedral, Entanglement in single-particle systems, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 463, 2277 (2007).
  * Cioslowski and Pernal [2006] J. Cioslowski and K. Pernal, Wigner molecules: The strong-correlation limit of the three-electron harmonium, The Journal of Chemical Physics 125, 064106 (2006), https://doi.org/10.1063/1.2222361 .
  * Laguna and Sagar [2011] H. G. Laguna and R. P. Sagar, Statistical correlations in the moshinsky atom, Phys. Rev. A 84, 012502 (2011).
  * Nagy and Aldazabal [2011] I. Nagy and I. Aldazabal, Metric measures of interparticle interaction in an exactly solvable two-electron model atom, Phys. Rev. A 84, 032516 (2011).
  * Pipek and Nagy [2009] J. Pipek and I. Nagy, Measures of spatial entanglement in a two-electron model atom, Phys. Rev. A 79, 052501 (2009).
  * Bouvrie  _et al._ [2012] P. A. Bouvrie, A. P. Majtey, A. R. Plastino, P. Sánchez-Moreno, and J. S. Dehesa, Quantum entanglement in exactly soluble atomic models: the moshinsky model with three electrons, and with two electrons in a uniform magnetic field, The European Physical Journal D 66, 10.1140/epjd/e2011-20417-4 (2012).
  * Johnson and Payne [1992] N. Johnson and M. Payne, Exactly solvable model of interacting particles in a quantum dot, Superlattices and Microstructures 11, 309 (1992).
  * Teichmann  _et al._ [2013] K. Teichmann, M. Wenderoth, H. Prueser, K. Pierz, H. W. Schumacher, and R. G. Ulbrich, Harmonic oscillator wave functions of a self-assembled inas quantum dot measured by scanning tunneling microscopy, Nano Letters 13, 3571 (2013), pMID: 23777509, https://doi.org/10.1021/nl401217q .
  * Bartlett and Wiseman [2003] S. D. Bartlett and H. M. Wiseman, Entanglement constrained by superselection rules, Phys. Rev. Lett. 91, 097903 (2003).
  * Nielsen and Chuang [2000] M. Nielsen and I. Chuang,  _Cambridge University Press_ (2000).
  * Peres [1996] A. Peres, Separability criterion for density matrices, Phys. Rev. Lett. 77, 1413 (1996).
  * Henderson and Vedral [2001] L. Henderson and V. Vedral, Classical, quantum and total correlations, Journal of Physics A: Mathematical and General 34, 6899 (2001).
  * Piani [2009] M. Piani, Relative entropy of entanglement and restricted measurements, Physical Review Letters 103, 1 (2009), arXiv:0904.2705 .
  * Adesso  _et al._ [2016] G. Adesso, M. Cianciaruso, and T. R. Bromley, Quantum discord and nonclassical correlations beyond entanglement, in  _Quantum Information_ (John Wiley and Sons, Ltd, 2016) Chap. 9, pp. 175–194.
  * Schilling [2013] C. Schilling, Natural orbitals and occupation numbers for harmonium: Fermions versus bosons, Phys. Rev. A 88, 042105 (2013).
  * Tennie  _et al._ [2016] F. Tennie, V. Vedral, and C. Schilling, Pinning of fermionic occupation numbers: Higher spatial dimensions and spin, Phys. Rev. A 94, 012120 (2016).
  * Amosov and Filippov [2016] G. G. Amosov and S. N. Filippov, Spectral properties of reduced fermionic density operators and parity superselection rule, Quantum Information Processing 16, 10.1007/s11128-016-1467-9 (2016).
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Ueber das paulische equivalenzverbot, Zeitschrift fuer Physik 47, 631 (1928).
  * Wick  _et al._ [1952] G. C. Wick, A. S. Wightman, and E. P. Wigner, The intrinsic parity of elementary particles, Phys. Rev. 88, 101 (1952).
  * Wick  _et al._ [1970] G. C. Wick, A. S. Wightman, and E. P. Wigner, Superselection rule for charge, Phys. Rev. D 1, 3267 (1970).
  * Vidal  _et al._ [2021] N. T. Vidal, M. L. Bera, A. Riera, M. Lewenstein, and M. N. Bera, Quantum operations in an information theory for fermions, Phys. Rev. A 104, 032411 (2021).
  * Mourik  _et al._ [2012] V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Signatures of majorana fermions in hybrid superconductor-semiconductor nanowire devices, Science 336, 1003 (2012), https://www.science.org/doi/pdf/10.1126/science.1222360 .
  * Frolov [2021] S. Frolov, Quantum computing’s reproducibility crisis: Majorana fermions, Nature 592, 350 (2021).
  * Bartlett  _et al._ [2007a] S. D. Bartlett, T. Rudolph, and R. W. Spekkens, Reference frames, superselection rules, and quantum information, Reviews of Modern Physics 79, 555–609 (2007a).
  * Sanders  _et al._ [2003] B. C. Sanders, S. D. Bartlett, T. Rudolph, and P. L. Knight, Photon-number superselection and the entangled coherent-state representation, Physical Review A 68, 10.1103/physreva.68.042329 (2003).
  * Bartlett  _et al._ [2007b] S. D. Bartlett, T. Rudolph, and R. W. Spekkens, Reference frames, superselection rules, and quantum information, Rev. Mod. Phys. 79, 555 (2007b).
  * Vollbrecht and Werner [2001] K. G. H. Vollbrecht and R. F. Werner, Entanglement measures under symmetry, Phys. Rev. A 64, 062307 (2001).
  * Ding [2020] L. Ding,  _Fermionic Entanglement and Correlation_ , Master’s thesis, LMU München (2020).
  * Boguslawski and Tecmer [2015] K. Boguslawski and P. Tecmer, Orbital entanglement in quantum chemistry, International Journal of Quantum Chemistry 115, 1289 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24832 .
  * Girard  _et al._ [2015] M. W. Girard, Y. Zinchenko, S. Friedland, and G. Gour, Erratum: Numerical estimation of the relative entropy of entanglement [phys. rev. a 82, 052336 (2010)], Phys. Rev. A 91, 029901 (2015).
  * Horodecki [1997] P. Horodecki, Separability criterion and inseparable mixed states with positive partial transposition, Physics Letters A 232, 333 (1997).

##  Appendix A Hermite Polynomials

The Hermite polynomials are an orthogonal polynomial series. Different
conventions of indexing and prefactors can be found in the literature. For the
purposes of this paper the $\nu$-th Hermite polynomial is defined as follows:

$H_{\nu}(x)=(-1)^{\nu}e^{x^{2}}\frac{d^{\nu}}{dx^{\nu}}e^{-x^{2}}.$ (32)

The first few Hermite polynomials read:

$\displaystyle H_{0}(x)$ $\displaystyle=1,$ (33) $\displaystyle H_{1}(x)$
$\displaystyle=2x,$ $\displaystyle H_{2}(x)$ $\displaystyle=4x^{2}-2,$
$\displaystyle H_{3}(x)$ $\displaystyle=8x^{3}-12x,$ $\displaystyle H_{4}(x)$
$\displaystyle=16x^{4}-48x^{2}+12,$ $\displaystyle H_{5}(x)$
$\displaystyle=32x^{5}-160x^{3}+120x,$ $\displaystyle H_{6}(x)$
$\displaystyle=64x^{6}-480x^{4}+720x^{2}-120.$

The $\nu$-th Hermite function $\phi_{\nu}^{(l)}(x)$ is defined as follows:

$\phi_{\nu}^{(l)}(x)=(2^{\nu}\nu!\sqrt{\pi}l)^{-\frac{1}{2}}e^{\frac{x^{2}}{2l^{2}}}H_{\nu}\left(\frac{x}{l}\right).$
(34)

##  Appendix B Integrals for expectation values of mode reduced density
operators

The matrix elements of mode- and mode-mode reduced density operators (cf. Eqs.
(18) and (19)) can be readily computed using the integral kernel expressions
of the one- and two-particles reduced density operators (cf. definition of
$\rho^{(p)}$ in Eq. (12)). More specifically, the diagonal elements of the
mode-reduced density operator for mode $\varphi_{m}(\vec{x})$ read:

$\displaystyle\left\langle f_{m}^{{\dagger}}f_{m}\right\rangle$
$\displaystyle=\int_{-\infty}^{+\infty}d\vec{x}_{1}d\vec{x}^{\prime}_{1}\rho^{(1)}(\vec{x}_{1},\vec{x}^{\prime}_{1})\varphi_{m}(\vec{x}_{1})\varphi_{m}(\vec{x}^{\prime}_{1}),$
(35) $\displaystyle\left\langle f_{m}f_{m}^{{\dagger}}\right\rangle$
$\displaystyle=\sum_{i\neq
m}\int_{-\infty}^{+\infty}d\vec{x}_{1}d\vec{x}^{\prime}_{1}\rho^{(1)}(\vec{x}_{1},\vec{x}^{\prime}_{1})\varphi_{i}(\vec{x}_{1})\varphi_{i}(\vec{x}^{\prime}_{1}).$

Similarly, for the fermionic mode-mode reduced density operator (19) one
obtains an expression for the lower diagonal matrix element:

$\displaystyle\left\langle
f_{i}^{{\dagger}}f_{i}f_{i}^{{\dagger}}f_{i}\right\rangle=$ (36)
$\displaystyle\int_{-\infty}^{+\infty}d\vec{x}_{1}d\vec{x}^{\prime}_{1}d\vec{x}_{2}d\vec{x}^{\prime}_{2}\bra{\varphi_{i},\varphi_{j}}\rho^{(2)}(\vec{x}_{1},\vec{x}^{\prime}_{1},\vec{x}_{2},\vec{x}^{\prime}_{2})\ket{\varphi_{i},\varphi_{j}}$
with
$\displaystyle\braket{\vec{x}_{1},\vec{x_{2}}}{\varphi_{i},\varphi_{j}}=\frac{1}{\sqrt{2}}(\varphi_{i}(\vec{x}_{1}),\varphi_{j}(\vec{x}_{2})-\varphi_{i}(\vec{x}_{2}),\varphi_{j}(\vec{x}_{1})),$
and
$\displaystyle\braket{\varphi_{i},\varphi_{j}}{\vec{x}_{1}^{\prime},\vec{x}_{2}^{\prime}}=\frac{1}{\sqrt{2}}(\varphi_{i}(\vec{x}_{1}^{\prime}),\varphi_{j}(\vec{x}_{2}^{\prime})-\varphi_{i}(\vec{x}_{2}^{\prime}),\varphi_{j}(\vec{x}_{1}^{\prime})).$

All remaining matrix elements can now be computed using the canonical
commutation/anti-commutation relations, e.g.,

$\displaystyle\left\langle f_{k}f_{k}^{{\dagger}}\right\rangle=\left\langle
f_{k}f_{k}^{{\dagger}}f_{l}f_{l}^{{\dagger}}\right\rangle+\left\langle
f_{k}f_{k}^{{\dagger}}f_{l}^{{\dagger}}f_{l}\right\rangle,$ (37)
$\displaystyle\left\langle f_{k}^{{\dagger}}f_{k}\right\rangle=\left\langle
f_{k}^{{\dagger}}f_{k}f_{l}f_{l}^{{\dagger}}\right\rangle+\left\langle
f_{k}^{{\dagger}}f_{k}f_{l}^{{\dagger}}f_{l}\right\rangle.$

Other matrix elements, for instance those given by expectation values of an
odd number of ladder operators, evidently vanish for pure states of fixed
particle number, such as the Harmonium ground state.

For bosonic Harmonium one can calculate analogous integral expressions by
considering the bosonic p-particle reduced density operator and replacing the
Slater determinants in Eq. (36) with permanents of basis functions.

##  Appendix C Explicit expressions for fermionic mode entanglement

Without imposing SSR’s, the mutual information and entropy of the mode reduced
density operator 28 read:

$\displaystyle I(\rho_{\\{m_{\uparrow},m_{\downarrow}\\}})$
$\displaystyle=2S(\rho_{\\{m_{\uparrow},m_{\downarrow}\\}})=-2\sum_{i}q_{i}\ln(q_{i}),$
(38) $\displaystyle E(\rho_{\\{m_{\uparrow},m_{\downarrow}\\}})$
$\displaystyle=S(\rho_{\\{m_{\uparrow},m_{\downarrow}\\}})=-\sum_{i}q_{i}\ln(q_{i}).$

A detailed calculation presented in Ref. [5] establishes the density operators
accessible under consideration of P-SSR and N-SSR. Under local parity
constraints, the physical part of the density operator $\rho^{P}$ can be
expressed as:

$\displaystyle\rho^{P}=$
$\displaystyle(q_{1}+q_{2})\ket{\Psi_{odd}}\bra{\Psi_{odd}}$ (39)
$\displaystyle+(q_{0}+q_{3})\ket{\Psi_{even}}\bra{\Psi_{even}},$
$\displaystyle\ket{\Psi_{odd}}=$
$\displaystyle\sqrt{\frac{q_{1}}{q_{1}+q_{2}}}\ket{m_{\uparrow}}\otimes\ket{\psi_{N-1}}$
$\displaystyle+\sqrt{\frac{q_{2}}{q_{1}+q_{2}}}\ket{m_{\downarrow}}\otimes\ket{\psi_{N-1}^{\prime}},$
$\displaystyle\ket{\Psi_{even}}=$
$\displaystyle\sqrt{\frac{q_{0}}{q_{0}+q_{3}}}\ket{\Omega}\otimes\ket{\psi_{N}}$
$\displaystyle+\sqrt{\frac{q_{3}}{q_{0}+q_{3}}}\ket{m_{\uparrow},m_{\downarrow}}\otimes\ket{\psi_{N-2}}.$

Similarly, the introduction of local particle number constraints gives rise to
the physically accessible density operator $\rho^{N}$:

$\displaystyle\rho^{N}$
$\displaystyle=q_{0}\ket{\Psi_{0}}\bra{\Psi_{0}}+(q_{1}+q_{2})\ket{\Psi_{1}}\bra{\Psi_{1}}+q_{3}\ket{\Psi_{2}}\bra{\Psi_{2}},$
(40) $\displaystyle\ket{\Psi_{0}}$
$\displaystyle=\ket{\Omega}\otimes\ket{\psi_{N}},$
$\displaystyle\ket{\Psi_{1}}$ $\displaystyle=\ket{\Psi_{odd}},$
$\displaystyle\ket{\Psi_{2}}$
$\displaystyle=\ket{m_{\uparrow},m_{\downarrow}}\otimes\ket{\psi_{N-2}}.$

Consequently, the physically accessible mutual information can be written as:

$\displaystyle
I^{PSSR}(\rho_{\\{m_{\uparrow},m_{\downarrow}\\}})=-2\sum_{i}q_{i}\ln(q_{i})+(q_{1}+q_{2})\ln(q_{1}+q_{2})$
(41) $\displaystyle+(q_{0}+q_{3})\ln(q_{0}+q_{3}),$ $\displaystyle
I^{NSSR}(\rho_{\\{m_{\uparrow},m_{\downarrow}\\}})=-2\sum_{i}q_{i}\ln(q_{i})+(q_{1}+q_{2})\ln(q_{1}+q_{2})$
$\displaystyle+q_{0}\ln(q_{0})+q_{3}\ln(q_{3}).$

The physically accessible entanglement between a spin-ful mode subsystem and
the remaining system is given by:

$\displaystyle E^{PSSR}(\rho_{\\{m_{\uparrow},m_{\downarrow}\\}})=$
$\displaystyle-\sum_{i}q_{i}\ln(q_{i})+(q_{1}+q_{2})\ln(q_{1}+q_{2})$ (42)
$\displaystyle+(q_{0}+q_{3})\ln(q_{0}+q_{3}),$ $\displaystyle
E^{NSSR}(\rho_{\\{m_{\uparrow},m_{\downarrow}\\}})=$
$\displaystyle(q_{1}+q_{2})\ln(q_{1}+q_{2})-q_{2}\ln(q_{2})$ $\displaystyle-
q_{1}\ln(q_{1}).$

It shall be stressed, that only due to the spin-symmetries of this particular
configuration it becomes possible to determine analytic expressions of the
relative entropy of entanglement. In general, in order to find the closest
separable state, a minimisation over a possibly large parameter space is
necessary.

##  Appendix D Analytic expressions for the coefficients of the fermionic mode
reduced density operator

One can express the coefficients of the mode-reduced density operator in Eq.
(29) as a function of the dimensionless constant $\kappa$ (cf. Eq. (11)). It
is illustrative to consider the three-fermion Harmonium ground state in one
spatial dimension, as the expressions become increasingly complex for
increasing particle numbers and dimensions. For mode $0$ the coefficients of
the mode RDO of three-Harmonium read:

$\displaystyle q_{1}$
$\displaystyle=\frac{2\sqrt{a}\sqrt{\frac{a(a-3b)}{a-2b}}}{\sqrt{\frac{2a^{2}-5ab+b^{2}}{a-2b}}\sqrt{\frac{a\left(4a^{2}-12ab+5b^{2}\right)}{2a^{2}-5ab+b^{2}}}},$
(43) $\displaystyle q_{2}$
$\displaystyle=\frac{2a^{3/2}\sqrt{\frac{a(a-3b)}{a-2b}}\left(4a^{2}-12ab+7b^{2}\right)}{(a-2b)\left(\frac{2a^{2}-5ab+b^{2}}{a-2b}\right)^{3/2}\left(\frac{a\left(4a^{2}-12ab+5b^{2}\right)}{2a^{2}-5ab+b^{2}}\right)^{3/2}},$
$\displaystyle q_{3}$
$\displaystyle=\frac{\sqrt{a-3b}(2a-3b)^{2}\sqrt{\frac{4a^{2}-6ab+b^{2}}{a-b}}\sqrt{\frac{a\left(2a^{2}-4ab+b^{2}\right)}{4a^{2}-6ab+b^{2}}}}{4a(a-2b)^{2}(a-b)^{\frac{3}{2}}}$
$\displaystyle\times\sqrt{\frac{a\left(a^{2}-3ab+2b^{2}\right)}{4a^{2}-10ab+5b^{2}}}\sqrt{\frac{a\left(4a^{2}-10ab+5b^{2}\right)}{2a^{2}-4ab+b^{2}}},$
$\displaystyle q_{0}$ $\displaystyle=1-q_{1}-q_{2}-q_{3},$

where $a=\frac{1}{2\tilde{l}^{2}}$ and
$b=-\frac{1}{6}\left(\frac{1}{l^{2}}-\frac{1}{\tilde{l}^{2}}\right)$ and $N=3$
in this particular case.

##  Appendix E Fermionic spatial mode entanglement: additional plots

In this appendix, we display additional graphics on mode entanglement. These
are Figs. 10 and 11.

Figure 10: Relative amount of physically accessible entanglement $E^{P/N}/E$
of spatial modes $\phi_{0}$ to $\phi_{6}$ as a function of the coupling
strength $\kappa$ in one-dimensional 3-particle Harmonium. Figure 11: Relative
amount of physically accessible entanglement $E^{P/N}/E$ of spatial modes
$\phi_{0}$ to $\phi_{6}$ as a function of the coupling strength $\kappa$ in
one-dimensional 5-particle Harmonium.

##  Appendix F Fermionic mode-mode entanglement: deriving the two mode-reduced
density operator and introducing local SSR’s

The fermionic two-mode reduced density operator is $16\times 16$ dimensional
and can be divided into $5$ block diagonal sectors of fixed global particle
number $\\{N=0,1,2,3,4\\}$. The elements of the two-mode reduced density
operator can be written as expectation values of the creation and annihilation
operators which allows for an explicit matrix representation. For the local
basis
$\\{{\ket{i_{\uparrow}},\ket{i_{\downarrow}},\ket{j_{\uparrow}},\ket{j_{\downarrow}}\\}}$
with corresponding creation operators
$\\{{f_{1}^{\dagger},f_{2}^{\dagger},f_{3}^{\dagger},f_{4}^{\dagger}\\}}$ this
yields:

$\displaystyle($
$\displaystyle\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}})_{n_{1}^{\prime},n_{2}^{\prime},n_{3}^{\prime},n_{4}^{\prime}}^{n_{1},n_{2},n_{3},n_{4}}=$
(44) $\displaystyle\langle$
$\displaystyle(f_{1}^{\dagger})^{n_{1}^{\prime}}(f_{2}^{\dagger})^{n_{2}^{\prime}}(f_{3}^{\dagger})^{n_{3}^{\prime}}(f_{4}^{\dagger})^{n_{4}^{\prime}}f_{4}f_{3}f_{2}f_{1}\cdot$
$\displaystyle\cdot
f_{1}^{\dagger}f_{2}^{\dagger}f_{3}^{\dagger}f_{4}^{\dagger}(f_{4}^{\dagger})^{n_{4}}(f_{3}^{\dagger})^{n_{3}}(f_{2}^{\dagger})^{n_{2}}(f_{1}^{\dagger})^{n_{1}}\rangle$

The super-indices represent the occupation of various spin-ful modes. Non-
vanishing expectation values only exist for blocks of fixed global particle
number $N$ such that
$n_{1}+n_{2}+n_{3}+n_{4}=n_{1}^{\prime}+n_{2}^{\prime}+n_{3}^{\prime}+n_{4}^{\prime}=N$.
For example, the expectation value corresponding to the matrix element
$\ket{i_{\uparrow,\downarrow}}\otimes\ket{\Omega}\bra{i_{\uparrow,\downarrow}}\otimes\bra{\Omega}$
would correspond to the indices
$n_{1}=n_{2}=n_{1}^{\prime}=n_{2}^{\prime}=1,n_{3}=n_{4}=n_{3}^{\prime}=n_{4}^{\prime}=0$.
Practically speaking, expanding the expression in Eq. (44) gives rise to
integrals of the various particle reduced density operators (cf. Eq. (12)), in
analogy to the explicit expressions outlined in App. B.

The introduction of a local number or parity constraint severely constrains
the amount of entanglement between the modes. The physically accessible
entanglement arises from particular sub-sectors of the two-mode reduced
density operator and reads:

$\displaystyle E(\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}}^{N})$
$\displaystyle=S(\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}}^{N}\|\sigma_{N})$
(45) $\displaystyle
E(\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}}^{P})$
$\displaystyle=S(\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}}^{N}\|\sigma_{N})+S(\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}}^{even}\|\sigma_{even}).$

Here, the sector with fixed local particle number
$\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}}^{N}=\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}}^{odd}$
coincides with the odd particle number sector and is spanned by:
$\\{{\ket{i_{\uparrow}}\otimes\ket{j_{\uparrow}},\ket{i_{\uparrow}}\otimes\ket{j_{\downarrow}},\ket{i_{\downarrow}}\otimes\ket{j_{\uparrow}},\ket{i_{\downarrow}}\otimes\ket{j_{\downarrow}}\\}}$.
The sector with even local parity
$\rho_{m_{\uparrow\downarrow}m_{\uparrow\downarrow}}^{even}$ is spanned by
$\\{{\ket{\Omega}\otimes\ket{\Omega},\ket{i_{\uparrow\downarrow}}\otimes\ket{\Omega},\ket{\Omega}\otimes\ket{j_{\uparrow\downarrow}},\ket{i_{\uparrow,\downarrow}}\otimes\ket{j_{\uparrow\downarrow}}\\}}$.
The PPT-criterion here is a sufficient condition for separability [27] and
thus, we use it to determine the closest separable states $\sigma_{N}$ and
$\sigma_{even}$.