Mode entanglement in fermionic and bosonic Harmonium

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.


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].
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⟩ ± |10⟩), 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,25]).Previous work has considered particle entanglement of the Harmonium ground state [26].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 [27].
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.

Correlation and Entanglement Measures
The state of any quantum system can be described by a positive semi-definite density operator ρ = i p i |ψ i ⟩ ⟨ψ i | with unit trace i p i = 1 [28].Physical observables are associated with operators Ô which form an algebra A(H).The expectation value of any observable Ô with respect to ρ reads: Provided that the total Hilbert space can be expressed as a product of two Hilbert spaces H = H A ⊗ H B , and that the set of physical observables form an algebra which correspondingly decomposes as , the reduced density operator ρ A = tr B (ρ) of a particular subsystem A is uniquely determined by the following condition on the observables ÔA [28]: 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 2.1 (Uncorrelated states) A state represented by density operator ρ is uncorrelated iff: The set of uncorrelated states is denoted by S unc .A state is correlated iff its density operator ρ ̸ ∈ S unc .
The notion of separable states allows one to distinguish classical correlations and quantum entanglement.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. ‡ Classical correlations are quantified by the Shannon entropy.Its quantum analogue is the von Neumann entropy: where λ i are the eigenvalues of ρ [28].The quantum relative entropy is a measure which captures our ability to distinguish two quantum states ρ and σ: The quantum mutual information I(ρ) is a widely accepted measure of the total correlation [30,5] between two subsystems A and B: (5) Geometrically, the quantum mutual information describes the distance to the closest uncorrelated state [5] I(ρ) = min σ∈Sunc S(ρ∥σ).We consider the relative entropy of entanglement, E [31] as our measure of ‡ The separability of a state can be assessed by a number of criteria.In this paper, we shall use the PPT criterion [29] which is a necessary and sufficient condition for separability in twodimensional subsystems and provides a lower bound for systems in larger state spaces.entanglement.E is defined as the quantum relative entropy of the state ρ and the closest separable state σ * ∈ S sep : E is strictly positive if and only if ρ is non-separable, i.e. entangled [31].Geometrically, the relative entropy of entanglement describes the distance to the closest separable state σ * , cf. the dotted red line in Fig. 1.
The relative entropy of entanglement for a bipartite system in a pure state simplifies to E(ρ) = S(ρ A/B ).
The classical correlation of a state ρ is defined as the quantum relative entropy of the separable state σ * and the closest uncorrelated state ρ A ⊗ ρ B [30]: The classical correlation can be visualised by the geometrical distance between σ * and the closest uncorrelated state.
In general, the amount of entanglement and classical correlation does not sum up to the total correlation (I ̸ = E + C) as mixed quantum states give rise to correlation described by the quantum discord [32].

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: 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 {|ϕ 1 i ⟩ D i=1 } form a basis of modes of the single-particle Hilbert space H (1) .Observations confined to a subset of modes i=1 } are related to a decomposition of the singleparticle Hilbert space into an orthogonal sum: 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 H A and H B will allow us to construct a product Hilbert space that is isomorphic to the total Hilbert space.

Fermionic Mode-Reduced Density Operators
Fermionic annihilation and creation operators satisfy canonical anti-commutation rules [34]: They can be used to define occupation number states, which form a basis of the Fock space F associated with Here, |Ω⟩ 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 F A and F B .A natural isomorphism between the product of these Fock spaces and the total Fock space can be constructed by mapping number states as follows: 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. (10).
Given the product structure of the isomorphic Hilbert space F A ⊗ 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 anticommutation relations (10), matrix representations of the reduced density operators can be formed out of expectation values.
The mode-reduced density matrix of mode |ϕ m ⟩ in the local basis {|Ω⟩ , f † m |Ω⟩} reads [33]: The mode-mode reduced density matrix of modes 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 ϕ is given by Eq. ( 15) with i ≡ {ϕ, ↑} and j ≡ {ϕ, ↓}.The reduced density operator of two spatial modes becomes a 16 × 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. 3 and the evaluation of the matrix elements for the Harmonium ground state is outlined in Appendix C.

2.2.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]: These can be used to define an (orthonormal) numberstate basis of the bosonic Fock-space B: where 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 modereduced density operator of mode 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 , ρ {m} is a (N + 1) × (N + 1) matrix.
Similarly, bosonic mode-mode reduced density operators ρ {i,j} can be expressed as expectation values of creation and annihilation operators.However, the size of ρ {i,j} scales quadratically with the particle number N .
Unlike particle-reduced density operators, modereduced density operators can be mixed states of various particle numbers, as evident by Eqs. ( 14), ( 15) and (18).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 Sec.2.4.

N -Harmonium
The Hamiltonian of the N -Harmonium model system in n spatial dimensions reads: where Ω = diag(ω 1 , • • • , ω n ) parametrises the harmonic trap strength and K parametrises the particleparticle interactions; (⃗ x 1 , . . ., ⃗ x N ) and (⃗ p 1 , . . ., ⃗ p N ) denote particle positions and their respective momenta.Coupling can be attractive K < 0 or repulsive K > 0. One requires −Ω 2 mN < K for the existence of bound states.The Hamiltonian can be separated and diagonalised by a suitable coordinate transformation (cf.Ref. [35] for further details).
The bosonic ground state has the profile of a Gaussian and reads: with ⃗ X ≡ 1 N (⃗ x 1 + . . .⃗ x N ) representing the centre-ofmass vector, and constants A ≡ diag(A (1) , . . ., A (n) ), where where , and 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. [36]).More explicitly, the N -fermion ground state Ψ (f ) of the Harmonium model is given by: Here, the ϕ µ k represent products of Hermite functions of the relevant spatial dimensions: ) are Hermite functions of degree µ α and natural length scale lα .
The quantum number vectors µ 1 , • • • , µ 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 µ i can be identical.Note that l (α) is an intrinsic length scale and one can set l (1) = 1 without loss of generality.The lengths l(α) are then determined by the dimensionless coupling constants: κ (α) is bounded, as one requires κ (α) > −1 for bound state solutions.We can thus express the bosonic and fermionic ground states as a function of the dimensionless constant κ (α ).

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 |ψ⟩ = α |e − ⟩ + β |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, [ N , ρ] = 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]: where k ∈ [0, 1] for fermions, whereas k ∈ [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.

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) [37] and the local particle number-superselection rule (N-SSR) [38].P-SSR precludes a superposition of even and oddnumber fermion states.Under local conservation of parity q (P-SSR), observables ÔA/B must satisfy: where the P q represent projectors onto subspaces of definite parity number.Therefore, the observables ÔP A/B are block-diagonal and the algebra of observables decomposes into an even an odd parity sector A ee + A oo .Note that a violation of local parity superselection would have wide-ranging consequences such as super-luminal signalling [39,5,10] and is therefore considered impossible.
Local N-SSR is an even more stringent constraint and it requires the local observables ÔN A/B to be blockdiagonal in the particular fermion-N-particle number sectors: where the P N denote the projectors onto the local Nfermion 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-luminalsignalling, 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. [40]).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 [41].

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,42]).In general, one must discern systems of massive and massless bosons.For reasons similar to the ones outlined in Sec.2.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 blockdiagonal 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 [42,43,44], 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.

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 [45]: Theorem 2.1 If a density operator ρ is invariant under local unitary operations U(g), representing elements of a group g ∈ G, with generator T G , s.t.U † (g)ρU (g) = ρ, then the closest separable state with respect to ρ, σ * , possesses the same symmetries as ρ, such that T G (σ * ) = σ * .
Without proving the statement which is done elsewhere (cf.[46]), 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.

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 modereduced density operators are accessible by physical measurements (cf.Sec.2.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.[24,25]), our results demonstrate the need to re-evaluate the notion of entanglement in various other continuous systems.
Hermite functions constitute a natural singlemode basis of the Harmonium model.
In what follows, we shall use bases constructed out of Hermite functions ϕ ( l) m of natural length scale l (cf.Appendix B).Note that by Eq. ( 22), selecting a different arbitrary length scale l would merely amount to 'stretching' or 'quenching' the κ 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.
for ease of notation.

Fermions
In the presence of an external magnetic field, an additional Zeeman term H B = c S ℏ N i Ŝi • ⃗ B ext must be added to the Harmonium Hamiltonian [36], where Ŝ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 φi (x; ↑ / ↓) = ϕ i (x) ⊗ |↑/↓⟩.In spin-ful Harmonium, the ground state can be determined analogously to Eq. (( 21)) by replacing ϕ with φ.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. |B ext | ≫ ℏ cs ω √ 1 + κ, the system's spin degree of freedom will be quenched and all spins aligned.The resulting ground state equals the ground state of spinless 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 ℏ cs ω √ 1 + κ < |B ext | < 2ℏ cs ω √ 1 + κ, the groundstate eigenfunction for 3-Harmonium in one spatial dimension (n = 1) reads: The diagonal representation arises from the spin symmetries and fixed total particle number of the Harmonium ground state (cf.Eq. ( 15)).The structure of the Schmidt-decomposed N -particle state reads: (29) In our analysis, we have determined the amplitudes √ q i analytically.This has become possible since the integral expressions for expectation values of fermionic ladder operators (cf.Appendix C) 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 spatialmode 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 Appendix D, cf.Eqs.(D.1), (D.4) and (D.5).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.Appendix E 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.
The broad hierarchy and variation of the spatial mode entanglement and correlation measures with κ (cf.Eq. ( 22)) is qualitatively similar for various modes as shown in Fig. 2. As an example, this figure displays correlation and entanglement of modes ϕ 0 and ϕ 2 in one-dimensional 3-fermion Harmonium.Unsuspectedly, all correlation and entanglement measures assume a local maximum at a finite positive value of κ 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 particleparticle interaction strength κ, 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.§ 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 ϕ 1 and ϕ 2 , the relative amount of entanglement accessible under P-SSR approaches zero for small interaction strengths.For large κ ≳ 10 3 , E P /E is hierarchically ordered with the mode number, which is not the case in the intermediate regime κ ≈ 1 . . .100.Interestingly, the relative amount of entanglement of mode ϕ 1 accessible under P-SSR remains approximately constant, while for ϕ 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 κ.Remarkably, the magnitudes of the maxima at large κ grow with increasing mode number.∥For instance, E N of mode ϕ 6 at § At κ = 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.∥ One might speculate that this result is a reminiscence of the around κ ≈ 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 Appendix F. The variation of spatial mode entanglement in one-dimensional Harmonium with the particle number is shown in Fig. 4. For mode ϕ 0 , the relative amount of entanglement under P-SSR and N-SSR assumes impact of SSRs becoming negligible in the many-copy limit, cf.Ref. [27].Harmonium as a function of trapping frequency detuning ω y /ω x and coupling strength κ.In two spatial dimensions, the basis of spatial modes consists of products of two Hermite functions, i.e. ϕ (i,j) = ϕ (i) ϕ (j) .Note that ω y /ω 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.a minimum around vanishing κ.Throughout, E N is smaller than E P and only amounts to less than approximately 10% of the total entanglement.For mode ϕ 2 , the relative amount of entanglement assumes a maximum around vanishing κ 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 contrast to the particle entanglement results of Ref. [26].
In two-dimensional Harmonium, the form of the ground state is parametrised by the detuning or quench of the trapping potential frequencies ω y /ω x and the relative interaction strength κ.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.ω y /ω 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 κ.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 .

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 × 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 [47] 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 modemode 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 (10) the matrix elements can be expressed in terms of particlereduced density operators (cf.Appendix G).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.2.5).The optimisation is carried out using the PPT criterion [29,48] 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 × 2. When evaluating E, the criterion only provides a lower bound [49].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, cf.Eq. (G.2) for an explicit definition of E P/N .It has a minimum around vanishing coupling strength κ, followed by an increase in the entanglement to a maximum around κ ≈ 10 3 and a subsequent decrease for larger interaction strengths.Interestingly, mode-mode entanglement between the two modes ϕ 0 and ϕ 1 that form the ground state at κ = 0 never exceeds E ≈ 10 −5 , even in the regime of strong couplings.Coupling strength / κ Figure 6: Mode-mode entanglement between mode ϕ 0 and mode ϕ 1 (top, black), and between mode ϕ 0 and mode ϕ 2 (bottom, grey) in 4-fermion Harmonium as a function of κ.Note that the evaluation of the entanglement measures involves a parametrised minimisation.Here, we have used the PPT criterion [29] 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.

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.
3.1 is added to the Hamiltonian which causes a splitting of the eigenstates with respect to the magnetic quantum numbers σ 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.

2.2.2).
In the basis {(b † ) i |Ω⟩} 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 ρ {m} = diag(p 0 , p 1 , p 2 , p 3 ).Consequently, the Schmidt-decomposition of the three particle state can be written as: |ψ 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.Appendix C).
Due to the Schmidt-decomposition, the mutual information and entanglement read: 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 ϕ 0 to ϕ 8 is displayed as a function of κ.While vanishing for zero interaction strength, it attains a maximum at finite κ 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 ϕ 0 and ϕ 2 and N = 3, 4, 5.No clear hierarchy with the particle number is observed.

Mode-Mode Entanglement
To assess bosonic mode-mode entanglement, one has to determine the two-mode reduced density operator.For pure fixedparticle 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 ρ {m,n} by evaluating particle-reduced density operators in the basis . Note that once a local N-SSR is introduced, all local number sectors of ρ {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 [29,48].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.

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 [24].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.
(D.5) It shall be stressed, that only due to the spinsymmetries 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.

Figure 1 :
Figure 1: Visual representation of sets of quantum states and correlation measures [5].

2 10 3 10 4 10 5 Figure 2 :
Figure2: Correlation and entanglement of spatial modes ϕ 0 (top), ϕ 2 (bottom) in one-dimensional 3fermion Harmonium as a function of the relative particle-particle interaction strength κ.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 :
Figure 3: Relative amount of physically accessible entanglement E P/N /E of spatial modes ϕ 0 to ϕ 6 as a function of the coupling strength κ in one-dimensional 4-fermion Harmonium.

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

Figure 5 :
Figure5: Relative amount of physically accessible entanglement E P/N /E in two-dimensional 4-fermion Harmonium as a function of trapping frequency detuning ω y /ω x and coupling strength κ.In two spatial dimensions, the basis of spatial modes consists of products of two Hermite functions, i.e. ϕ (i,j) = ϕ (i) ϕ (j) .Note that ω y /ω 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.

Figure 8 :
Figure 8: E of spatial modes m 0 (top) and m 2 (bottom) as a function of κ for one-dimensional Harmonium with different boson numbers N = 3, 4, 5.

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