
The lowest excited configuration of harmonium

C. L. Benavides-Riveros,1,2 J. M. Gracia-Bondı́a2,3 and J. C. Várilly4

1 Zentrum für Interdisziplinäre Forschung, Wellenberg 1, Bielefeld 33615, Germany

2 Departamento de Fı́sica Teórica, Universidad de Zaragoza, 50009 Zaragoza, Spain

3 Instituto de Fı́sica Teórica, CSIC–UAM, Madrid 28049, Spain

4 Escuela de Matemática, Universidad de Costa Rica, 11501 San José, Costa Rica

Phys. Rev. A 86 (2012), 022525

Abstract

The harmonium model has long been regarded as an exactly solvable laboratory bench for
quantum chemistry [1]. For studying correlation energy, only the ground state of the system
has received consideration heretofore. This is a spin singlet state. In this work we exhaustively
study the lowest excited (spin triplet) harmonium state, with the main purpose of revisiting the
relation between entanglement measures and correlation energy for this quite different species.
The task is made easier by working with Wigner quasiprobabilities on phase space.

1 Introduction
Replacing the wave function of electronic systems by the reduced 2-body density matrix 𝛾2 tremen-
dously saves computation without losing relevant physical information. Till very recently, the
solutions to the 𝑁-representability for that matrix [2,3] were impractical. This certainly did not im-
pede great advances in the use of 𝛾2 for many-electron quantum systems – see for instance [4]. Now
a constructive solution [5] to that representability problem, leading to a hierarchy of constraints [6]
on the variation space for 𝛾2, has been unveiled.

At any rate, the last fifteen years have witnessed a justifiable amount of work in trying to obtain
the 2-body matrix as a functional of the 1-body density matrix 𝛾1. Starting with the pioneer work
by Müller [7], several competing functionals have been designed, partly out of theoretical prejudice,
partly with the aim of improving predictions for particular systems. We shall discuss pure state
representability for 𝛾1 in the case of our interest in Section 6.

Two-electron systems are special in that 𝛾2 is known “almost exactly” in terms of 𝛾1. Let us
express 𝛾1 by means of the spectral theorem in terms of its natural orbitals and occupation numbers.
For instance, the ground state of the system admits a 1-density matrix:

𝛾1(𝒙, 𝒙′) =
(
↑1↑1′ + ↓1↓1′

)
𝛾1(𝒓, 𝒓′) =

(
↑1↑1′ + ↓1↓1′

) ∑︁
𝑖

𝑛𝑖 𝜙𝑖 (𝒓)𝜙∗𝑖 (𝒓′). (1)

1


Here
∑

𝑖 𝑛𝑖 = 1. Mathematically this a mixed state. The corresponding 2-density matrix is given by

𝛾2(𝒙1, 𝒙2; 𝒙′1, 𝒙
′
2) =

(
↑1↓2 − ↓1↑2

) (
↑1′↓2′ − ↓1′↑2′)

∑︁
𝑖 𝑗

𝑐𝑖𝑐 𝑗

2
𝜙𝑖 (𝒓1)𝜙𝑖 (𝒓2)𝜙∗𝑗 (𝒓′1)𝜙

∗
𝑗 (𝒓′2),

with coefficients 𝑐𝑖 = ±√𝑛𝑖 . (2)

The expression is exact, but the signs of the 𝑐𝑖 need to be determined to find the ground state [8,9].
Note that 𝛾2

2 = 𝛾2. The first excited state of the system admits a reduced 1-density matrix of the
kind:

𝛾1(𝒙; 𝒙′) = (spin factor) ×
∑︁
𝑖 𝑗

𝑛𝑖
(
𝜙2𝑖 (𝒓)𝜙∗2𝑖 (𝒓

′) + 𝜙2𝑖+1(𝒓)𝜙∗2𝑖+1(𝒓
′)
)

with
∑

𝑖 𝑛𝑖 = 1 and spin ∈ { ↑1↑1′ ,
1
2 (↑1↑1′ + ↓1↓1′), ↓1↓1′ }. The corresponding spinless 2-density

matrix 𝛾2(𝒓1, 𝒓2; 𝒓′1, 𝒓
′
2) is given by∑︁

𝑖 𝑗

𝑐𝑖𝑐 𝑗

2
[
𝜙2𝑖 (𝒓1)𝜙2𝑖+1(𝒓2)𝜙∗2 𝑗 (𝒓

′
1)𝜙

∗
2 𝑗+1(𝒓

′
2) + 𝜙2𝑖+1(𝒓1)𝜙2𝑖 (𝒓2)𝜙∗2 𝑗+1(𝒓

′
1)𝜙

∗
2 𝑗 (𝒓

′
2)

− 𝜙2𝑖 (𝒓1)𝜙2𝑖+1(𝒓2)𝜙∗2 𝑗+1(𝒓
′
1)𝜙

∗
2 𝑗 (𝒓

′
2) − 𝜙2𝑖+1(𝒓1)𝜙2𝑖 (𝒓2)𝜙∗2 𝑗 (𝒓

′
1)𝜙

∗
2 𝑗+1(𝒓

′
2)
]
,

with coefficients 𝑐𝑖 = +√𝑛𝑖 . Due to the antisymmetry of this state, there is no ambiguity in the
choice of sign.

A completely integrable analogue of a two-electron atom, here called harmonium, describes two
fermions interacting with an external harmonic potential and repelling each other by a Hooke-type
force; thus the harmonium Hamiltonian in Hartree-like units is

𝐻 =
𝑝2

1
2

+
𝑝2

2
2

+ 𝑘

2
(𝑟2

1 + 𝑟2
2) −

𝛿

4
𝑟2

12, (3)

where 𝑟12 := |𝒓1 − 𝒓2 |. This model is rooted in the history of quantum mechanics: Heisenberg first
invoked it to approach the spectrum of helium [1].

Several problems related with this model – although not quite the present one – are analytically
solved; and so it is tempting to employ it as a testing ground for methods used in other systems,
such as the helium series. Indeed, Moshinsky [10] reintroduced it with the purpose of calibrating
correlation energy. There is considerable interest nowadays on learning from harmonium, including
further study of correlation [11–13], approximation of functionals [14, 15], and beyond quantum
chemistry, questions of entanglement [16–19] and black hole entropy [20].

In the past, harmonium problems have been attacked with ordinary wave mechanics [21]. Now,
for the analysis of harmonium the phase space representation of quantum mechanics recommends
itself. The deep reason for this is the metaplectic invariance of that formalism [22], hidden in
the standard approach: this made it easy to solve the sign dilemma in the exact Löwdin–Shull–
Kutzelnigg formula [8,9] for 𝛾2 in terms of 𝛾1, for two-electron systems [23,24]. We come to this at
the end of the next section. Such a phase-space description was taken up first by Dahl [25], and then
developed, within the context of a phase-space density functional theory (WDFT), by Blanchard,
Ebrahimi-Fard and ourselves [23, 24, 26–28].

Our goal in this article is to understand, in WDFT terms, the first excited state of harmonium.
As for helium-like atoms, we expect it to be the lowest spin triplet state, to which we refer simply

2


as the triplet. Particularly we make clear the nonexistence of a phase dilemma in this situation, and
pinpoint the similarities and differences between the relative behavior of entropy and correlation
energy for the (spin singlet) ground state and for the triplet. Again, and essentially for the same
reason, WDFT shows its worth here – see Section 6.

The customary plan of the paper follows. In Section 2 we briefly recall for the benefit of the
reader our treatment for the singlet ground state; this helps to introduce the notation. Sections 3 and 4
deal with the general mathematical structure of triplet 1-body Wigner functions. Section 5 computes
the Wigner quasiprobabilities for the harmonium triplet. Section 6 deals with the corresponding
natural orbitals. In Section 7 the behaviour of the occupation numbers, obtained numerically, is
compared to that of the ground state. Section 8 continues this comparison in the setting of quantum
information theory. The relative correlation energy for the triplet is smaller than for the singlet, just
as is the purity parameter. The proportionality between entropy and correlation energy, observed
in the weak correlation limit for the singlet, fails for the triplet state. Section 9 is the conclusion.

2 Wigner natural orbitals for the harmonium ground state
Given any interference operator |Ψ⟩⟨Φ| acting on the Hilbert space of a two-electron system, we
denote

𝑃2ΨΦ(𝒓1, 𝒓2; 𝒑1, 𝒑2; 𝜍1, 𝜍2; 𝜍1′ , 𝜍2′) (4)

:=
∫

Ψ(𝒓1 − 𝒛1, 𝒓2 − 𝒛2; 𝜍1, 𝜍2)Φ∗(𝒓1 + 𝒛1, 𝒓2 + 𝒛2; 𝜍1′ , 𝜍2′) 𝑒2𝑖( 𝒑1·𝒛1+ 𝒑2·𝒛2) 𝑑𝒛1 𝑑𝒛2.

These are 4 × 4 matrices on spin space. When Ψ = Φ we speak of Wigner quasiprobabilities,
which are always real, and we write 𝑑2 for 𝑃2. The extension of this definition to mixed states is
immediate. The corresponding reduced 1-body functions are found by

𝑃1ΨΦ(𝒓1; 𝒑1; 𝜍1; 𝜍1′) = 2
∫

𝑃2ΨΦ(𝒓1, 𝒓2; 𝒑1, 𝒑2; 𝜍1, 𝜍2; 𝜍1′ , 𝜍2) 𝑑𝒓2 𝑑 𝒑2 𝑑𝜍2.

These are 2 × 2 matrices on spin space. When Ψ = Φ we write 𝑑1 for 𝑃1. The associated
spinless quantities 𝑑2(𝒓1, 𝒓2; 𝒑1, 𝒑2) and 𝑑1(𝒓; 𝒑) are obtained by tracing on the spin variables.
The marginals of 𝑑2 give the pairs densities 𝜌2(𝒓1, 𝒓2), 𝜋2( 𝒑1, 𝒑2). The marginals of 𝑑1 give
the electronic density, namely 𝜌(𝒓1) =

∫
𝑑1(𝒓1, 𝒑1) 𝑑 𝒑1, and the momentum density 𝜋( 𝒑1) =∫

𝑑1(𝒓1, 𝒑1) 𝑑𝒓1.
It should be obvious how to extend the definitions to 𝑁-electron systems and their reduced

quantities; the combinatorial factor for 𝑑𝑁 ↦→ 𝑑𝑛 is
(𝑁
𝑛

)
.

Putting together (2) and (1) with (4), one arrives [23] at:

𝑑2(𝒓1, 𝒓2; 𝒑1, 𝒑2; 𝜍1, 𝜍2; 𝜍1′ , 𝜍2′) = (spin factor) ×
∑︁
𝑖 𝑗

𝑐𝑖 𝑐 𝑗

2
𝜒𝑖 𝑗 (𝒓1; 𝒑1)𝜒𝑖 𝑗 (𝒓2; 𝒑2), (5)

and

𝑑1(𝒓1; 𝒑1; 𝜍1, 𝜍1′) = 2
∫

𝑑2(𝒓1, 𝒓2; 𝒑1, 𝒑2; 𝜍1, 𝜍2; 𝜍1′ , 𝜍2) 𝑑𝜍2 𝑑𝒓2 𝑑 𝒑2

=
(
↑1↑1′ + ↓1↓1′

) ∑︁
𝑖

𝑛𝑖 𝜒𝑖 (𝒓1; 𝒑1).

3


Here 𝑛𝑖 are the occupation numbers with
∑

𝑖 𝑛𝑖 = 1, the 𝜒𝑖 𝑗 the natural Wigner interferences and
𝜒𝑖 := 𝜒𝑖𝑖 denote the natural Wigner orbitals; the spin factor is that of (2). Evidently

(
↑1↑1′ + ↓1↓1′

)
is a rotational scalar. We replace it by 2 in what follows.

The relation 𝑐𝑖 = ±√𝑛𝑖 holds. In principle there still remains the problem of determining the
signs of the infinite set of square roots, to find the ground state. To recover 𝑑2 from 𝑑1 is no mean
feat, since it involves going from a statistical mixture to a pure state – see below.

Bringing in extracule and intracule coordinates, respectively given by

𝑹 =
1
√

2
(𝒓1 + 𝒓2), 𝒓 =

1
√

2
(𝒓1 − 𝒓2),

𝑷 =
1
√

2
( 𝒑1 + 𝒑2), 𝒑 =

1
√

2
( 𝒑1 − 𝒑2),

the harmonium Hamiltonian is rewritten:

𝐻 = 𝐻𝑅 + 𝐻𝑟 :=
𝑃2

2
+ 𝜔2𝑅2

2
+ 𝑝2

2
+ 𝜇2𝑟2

2
.

We have introduced the frequencies 𝜔 :=
√
𝑘 and 𝜇 :=

√
𝑘 − 𝛿. Assume 𝛿 < 𝑘 , so both “electrons”

remain in the potential well. For the harmonium ground state the (spinless) Wigner 2-body
quasiprobability is readily found [25]:

𝑑2(𝒓1, 𝒓2; 𝒑1, 𝒑2) =
1
𝜋6 exp

(
−2𝐻𝑅

𝜔

)
exp

(
−2𝐻𝑟

𝜇

)
. (6)

The reduced 1-body phase space quasiprobability for the ground state is thus obtained:

𝑑1(𝒓1; 𝒑1) =
2
𝜋3

(
4𝜔𝜇

(𝜔 + 𝜇)2

)3/2
𝑒−2𝑟2

1𝜔𝜇/(𝜔+𝜇)𝑒−2𝑝2
1/(𝜔+𝜇) .

For its natural orbital expansion, with 𝑖 integer ⩾ 0 and 𝐿𝑖 the corresponding Laguerre polynomial,
one finds [23]

𝑐2
𝑖 = 𝑛𝑖 =

4√𝜔𝜇(√
𝜔 + √

𝜇
)2 (√𝜔 − √

𝜇
√
𝜔 + √

𝜇

)2𝑖
=: (1 − 𝑡2) 𝑡2𝑖 ; (7)

𝑓𝑖 (𝒓1; 𝒑1) = 𝑓𝑖 (𝑥1; 𝑝1𝑥) 𝑓𝑖 (𝑦1; 𝑝1𝑦) 𝑓𝑖 (𝑧1; 𝑝1𝑧), where

𝑓𝑖 (𝑥; 𝑝𝑥) =
1
𝜋
(−1)𝑖𝐿𝑖

(
2
√
𝜔𝜇 𝑥2 + 2𝑝2

𝑥/
√
𝜔𝜇

)
𝑒−

√
𝜔𝜇 𝑥2−𝑝2

𝑥/
√
𝜔𝜇 .

The functions 𝑓𝑖 determine up to a phase the interferences: for 𝑗 ⩾ 𝑘 ,

𝑓 𝑗 𝑘 (𝑥, 𝑝𝑥) =
1
𝜋
(−1)𝑘

√
𝑘!√︁
𝑗!

(
2
√
𝜔𝜇 𝑥2 + 2𝑝2

𝑥/
√
𝜔𝜇

) ( 𝑗−𝑘)/2
× 𝑒−𝑖( 𝑗−𝑘)𝜗𝐿 𝑗−𝑘

𝑘

(
2
√
𝜔𝜇 𝑥2 + 2𝑝2

𝑥/
√
𝜔𝜇

)
𝑒−

√
𝜔𝜇 𝑥2−𝑝2

𝑥/
√
𝜔𝜇,

where 𝜗 := arctan
(
𝑝𝑥/

√
𝜔𝜇 𝑥

)
.

4


The 𝐿
𝑗−𝑘
𝑘

are associated Laguerre polynomials. The 𝑓𝑘 𝑗 are complex conjugates of the 𝑓 𝑗 𝑘 .
Now, with the alternating choice (unique up to a global sign):

𝑐𝑖 = (−)𝑖 √𝑛𝑖 =
√︁

1 − 𝑡2 (−𝑡)𝑖,

and the above 𝑓 𝑗 𝑘 , formula (5) does reproduce (6). This was originally proved in [23], and verified
by minimization in [24]; we refer the reader to those papers. Trivially, the same sign rule holds for
natural orbitals of the garden variety (2).

3 Generalities on the triplet state
For a general two-electron system in a triplet spin state the reduced 1-density possesses three
different spin factors, say

↑1↑1′ and 1
2
(
↑1↑1′ + ↓1↓1′

)
and ↓1↓1′ .

While the spatial function for the ground state is symmetric, and consequently its spin part anti-
symmetric, for the first excited state the situation is exactly the opposite: the spatial function is
antisymmetric and its spin part is symmetric. This leads to important differences between both
cases for the natural orbital decomposition.

General triplet states are describable in the form [8, 21]:

Ψ𝑡1(𝒓1, 𝒓2; 𝜍1, 𝜍2) = ↑1↑2

∑︁
𝑖 𝑗

1
2
𝑐𝑖 𝑗 [𝜓𝑖 (𝒓1)𝜓 𝑗 (𝒓2) − 𝜓 𝑗 (𝒓1)𝜓𝑖 (𝒓2)],

Ψ𝑡0(𝒓1, 𝒓2; 𝜍1, 𝜍2) =
1
√

2
(
↑1↓2 + ↓1↑2

) ∑︁
𝑖 𝑗

1
2
𝑐𝑖 𝑗 [𝜓𝑖 (𝒓1)𝜓 𝑗 (𝒓2) − 𝜓 𝑗 (𝒓1)𝜓𝑖 (𝒓2)],

Ψ𝑡,−1(𝒓1, 𝒓2; 𝜍1, 𝜍2) = ↓1↓2

∑︁
𝑖 𝑗

1
2
𝑐𝑖 𝑗 [𝜓𝑖 (𝒓1)𝜓 𝑗 (𝒓2) − 𝜓 𝑗 (𝒓1)𝜓𝑖 (𝒓2)],

where 𝑐𝑖 𝑗 = −𝑐 𝑗𝑖. Here {𝜓𝑖} is a complete orthonormal set. In the absence of magnetic fields, the
wave functions can be taken real. We thus assume that the matrix 𝐶 = [𝑐𝑖 𝑗 ] is real, as well as the
functions 𝜓𝑖. Wave function normalization gives rise to Tr(𝐶𝑡 𝐶) = ∑

𝑖 𝑗 𝑐
2
𝑖 𝑗
= 1.

For the spin part, a less conventional and more cogent description is found in terms of polarization
vectors and the correlation tensor [29, App. F]; however, it is hardly worthwhile to introduce it here.
So we shall be content with presenting the Wigner 2-body quasiprobabilities for triplet states under
the matrix form

𝑃2Ψ𝑡1Ψ𝑡1 = ↑1↑2↑1′↑2′ 𝑑2 =

©­­­«
𝑑2 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0

ª®®®¬ , 𝑃2Ψ𝑡 ,−1Ψ𝑡 ,−1 = ↓1↓2↓1′↓2′ 𝑑2 =

©­­­«
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 𝑑2

ª®®®¬ ,
𝑃2Ψ𝑡0Ψ𝑡0 =

1
2
(
↑1↓2 + ↓1↑2

) (
↑1′↓2′ + ↓1′↑2′

)
𝑑2 =

1
2

©­­­«
0 0 0 0
0 𝑑2 𝑑2 0
0 𝑑2 𝑑2 0
0 0 0 0

ª®®®¬ ;

5


where 𝑑2 is the spinless Wigner 2-body quasiprobability, given by the expression

𝑑2(𝒓1, 𝒓2; 𝒑1, 𝒑2)

=
1
4

∑︁
𝑖 𝑗 ,𝑘𝑙

𝑐𝑖 𝑗 𝑐𝑘𝑙

∫
[𝜓𝑖 (𝒓1 − 𝒛1)𝜓 𝑗 (𝒓2 − 𝒛2) − 𝜓 𝑗 (𝒓1 − 𝒛1)𝜓𝑖 (𝒓2 − 𝒛2)]

× [𝜓∗
𝑘 (𝒓1 + 𝒛1)𝜓∗

𝑙 (𝒓2 + 𝒛2) − 𝜓∗
𝑙 (𝒓1 + 𝒛1)𝜓∗

𝑘 (𝒓2 + 𝒛2)] 𝑒2𝑖( 𝒑1·𝒛1+ 𝒑2·𝒛2) 𝑑𝒛1 𝑑𝒛2

=
1
4

∑︁
𝑖 𝑗 ,𝑘𝑙

𝑐𝑖 𝑗 𝑐𝑘𝑙 [𝑃𝑖𝑘 (𝒓1; 𝒑1)𝑃 𝑗 𝑙 (𝒓2; 𝒑2) − 𝑃𝑖𝑙 (𝒓1; 𝒑1)𝑃 𝑗 𝑘 (𝒓2; 𝒑2)

− 𝑃 𝑗 𝑘 (𝒓1; 𝒑1)𝑃𝑖𝑙 (𝒓2; 𝒑2) + 𝑃 𝑗 𝑙 (𝒓1; 𝒑1)𝑃𝑖𝑘 (𝒓2; 𝒑2)] . (8)

By integrating out one set of coordinates, we obtain the 1-body quasiprobabilities:

𝑃1Ψ𝑡1Ψ𝑡1 = ↑↑′ 𝑑1 =

(
𝑑1 0
0 0

)
, 𝑃1Ψ𝑡 ,−1Ψ𝑡 ,−1 = ↓↓′ 𝑑1 =

(
0 0
0 𝑑1

)
,

𝑃1Ψ𝑡0Ψ𝑡0 =
1
2
(
↑↑′ + ↓↓′

)
𝑑1 =

1
2

(
𝑑1 0
0 𝑑1

)
.

Here 𝑑1 is the spinless 1-body quasidensity corresponding to the triplet:

𝑑1(𝒓; 𝒑) = 2
∫

𝑑2(𝒓, 𝒓2; 𝒑, 𝒑2) 𝑑𝒓2 𝑑 𝒑2

=
1
2

∑︁
𝑖 𝑗 ,𝑘𝑙

𝑐𝑖 𝑗 𝑐𝑘𝑙

∫
[𝑃𝑖𝑘 (𝒓; 𝒑)𝑃 𝑗 𝑙 (𝒓2; 𝒑2) − 𝑃𝑖𝑙 (𝒓; 𝒑)𝑃 𝑗 𝑘 (𝒓2; 𝒑2)

− 𝑃 𝑗 𝑘 (𝒓; 𝒑)𝑃𝑖𝑙 (𝒓2; 𝒑2) + 𝑃 𝑗 𝑙 (𝒓; 𝒑)𝑃𝑖𝑘 (𝒓2; 𝒑2)] 𝑑𝒓2 𝑑 𝒑2

=
1
2

∑︁
𝑖 𝑗 ,𝑘𝑙

𝑐𝑖 𝑗 𝑐𝑘𝑙 [𝑃𝑖𝑘 (𝒓; 𝒑) 𝛿 𝑗

𝑙
− 𝑃𝑖𝑙 (𝒓; 𝒑) 𝛿 𝑗

𝑘
− 𝑃 𝑗 𝑘 (𝒓; 𝒑) 𝛿𝑖𝑙 + 𝑃 𝑗 𝑙 (𝒓; 𝒑) 𝛿𝑖𝑘 ]

= 2
∑︁
𝑖 𝑗 ,𝑘

𝑐𝑖𝑘 𝑐 𝑗 𝑘 𝑃𝑖 𝑗 (𝒓; 𝒑) = 2
∑︁
𝑖 𝑗

𝑑𝑖 𝑗 𝑃𝑖 𝑗 (𝒓; 𝒑),

where 𝐷 = 𝐶𝐶𝑡 = −𝐶2 is a positive definite matrix.

4 The Schmidt decomposition of the triplet
Let 𝐶 be any real antisymmetric square matrix. It is well known that there exists a real orthogonal
matrix 𝑄 such that 𝐴 = 𝑄𝑡𝐶𝑄, with 𝐴 a real block-diagonal matrix:

𝐴 = diag[𝐴0, 𝐴1, . . . ], 𝐴0 = 0, 𝐴𝑖 =

(
0 𝑎𝑖
−𝑎𝑖 0

)
.

6


By convention, here 𝑎𝑖 ⩾ 0. Therefore∑︁
𝑖 𝑗 ,𝑘𝑙

𝑐𝑖 𝑗𝑐𝑘𝑙𝑃𝑖𝑘 (𝒓1; 𝒑1)𝑃 𝑗 𝑙 (𝒓2; 𝒑2)

=
∑︁

𝑖 𝑗 ,𝑘𝑙,𝑣𝑤

𝑎𝑣𝑎𝑤
[
𝑞𝑖,2𝑣𝑞 𝑗 ,2𝑣+1 − 𝑞𝑖,2𝑣+1𝑞 𝑗 ,2𝑣

] [
𝑞𝑘,2𝑤𝑞𝑙,2𝑤+1 − 𝑞𝑘,2𝑤+1𝑞𝑙,2𝑤

]
𝑃𝑖𝑘 (𝒓1; 𝒑1)𝑃 𝑗 𝑙 (𝒓2; 𝒑2)

=
∑︁

𝑖 𝑗 ,𝑘𝑙,𝑣𝑤

𝑎𝑣𝑎𝑤
[
𝑞𝑖,2𝑣𝑃𝑖𝑘 (𝒓1; 𝒑1)𝑞𝑘,2𝑤𝑞 𝑗 ,2𝑣+1𝑃 𝑗 𝑙 (𝒓2; 𝒑2)𝑞𝑙,2𝑤+1

− 𝑞𝑖,2𝑣𝑃𝑖𝑘 (𝒓1; 𝒑1)𝑞𝑘,2𝑤+1𝑞 𝑗 ,2𝑣+1𝑃 𝑗 𝑙 (𝒓2; 𝒑2)𝑞𝑙,2𝑤
− 𝑞𝑖,2𝑣+1𝑃𝑖𝑘 (𝒓1; 𝒑1)𝑞𝑘,2𝑤𝑞 𝑗 ,2𝑣𝑃 𝑗 𝑙 (𝒓2; 𝒑2)𝑞𝑙,2𝑤+1

+ 𝑞𝑖,2𝑣+1𝑃𝑖𝑘 (𝒓1; 𝒑1)𝑞𝑘,2𝑤+1𝑞 𝑗 ,2𝑣𝑃 𝑗 𝑙 (𝒓2; 𝒑2)𝑞𝑙,2𝑤
]
.

Let us now make the definition 𝜒𝑟 𝑝 (𝒓; 𝒑) :=
∑

𝑚𝑘 𝑞𝑚𝑟 𝑃𝑚𝑘 (𝒓; 𝒑) 𝑞𝑘 𝑝, so that

𝑃𝑚𝑘 (𝒓; 𝒑) =
∑︁
𝑟 𝑝

𝑞𝑚𝑟 𝜒𝑟 𝑝 (𝒓; 𝒑) 𝑞𝑘 𝑝 .

This is the set of Wigner natural orbitals, and has the following nice property:∫
𝜒𝑟 𝑝 (𝒓; 𝒑) 𝑑𝒓 𝑑 𝒑 =

∫ ∑︁
𝑚𝑘

𝑞𝑚𝑟𝑃𝑚𝑘 (𝒓; 𝒑)𝑞𝑘 𝑝 𝑑𝒓 𝑑 𝒑 =
∑︁
𝑚𝑘

𝑞𝑚𝑟𝑞𝑘 𝑝 𝛿
𝑚
𝑘 = 𝛿𝑟𝑝 .

Hence,∑︁
𝑖 𝑗 ,𝑘𝑙

𝑐𝑖 𝑗𝑐𝑘𝑙𝑃𝑖𝑘 (𝒓1; 𝒑1)𝑃 𝑗 𝑙 (𝒓2; 𝒑2)

=
∑︁
𝑣𝑤

𝑎𝑣𝑎𝑤
[
𝜒2𝑣,2𝑤 (𝒓1; 𝒑1)𝜒2𝑣+1,2𝑤+1(𝒓2; 𝒑2) − 𝜒2𝑣,2𝑤+1(𝒓1; 𝒑1)𝜒2𝑣+1,2𝑤 (𝒓2; 𝒑2)

− 𝜒2𝑣+1,2𝑤 (𝒓1; 𝒑1)𝜒2𝑣,2𝑤+1(𝒓2; 𝒑2) + 𝜒2𝑣+1,2𝑤+1(𝒓1; 𝒑1)𝜒2𝑣,2𝑤 (𝒓2; 𝒑2)
]
.

The other three summands in (8) yield the same expression. For instance, the third is

−
∑︁
𝑖 𝑗 ,𝑘𝑙

𝑐𝑖 𝑗𝑐𝑘𝑙𝑃𝑖𝑙 (𝒓1; 𝒑1)𝑃 𝑗 𝑘 (𝒓2; 𝒑2)

= −
∑︁

𝑖 𝑗 ,𝑘𝑙,𝑣𝑤

𝑎𝑣𝑎𝑤
[
𝑞𝑖,2𝑣𝑞 𝑗 ,2𝑣+1 − 𝑞𝑖,2𝑣+1𝑞 𝑗 ,2𝑣

] [
𝑞𝑘,2𝑤𝑞𝑙,2𝑤+1 − 𝑞𝑘,2𝑤+1𝑞𝑙,2𝑤

]
𝑃𝑖𝑙 (𝒓1; 𝒑1)𝑃 𝑗 𝑘 (𝒓2; 𝒑2)

= −
∑︁

𝑖 𝑗 ,𝑘𝑙,𝑣𝑤

𝑎𝑣𝑎𝑤
[
𝑞𝑖,2𝑣𝑃𝑖𝑙 (𝒓1; 𝒑1)𝑞𝑙,2𝑤+1𝑞 𝑗 ,2𝑣+1𝑃 𝑗 𝑘 (𝒓2; 𝒑2)𝑞𝑘,2𝑤

− 𝑞𝑖,2𝑣𝑃𝑖𝑙 (𝒓1; 𝒑1)𝑞𝑙,2𝑤𝑞 𝑗 ,2𝑣+1𝑃 𝑗 𝑘 (𝒓2; 𝒑2)𝑞𝑘,2𝑤+1

− 𝑞𝑖,2𝑣+1𝑃𝑖𝑙 (𝒓1; 𝒑1)𝑞𝑙,2𝑤+1𝑞 𝑗 ,2𝑣𝑃 𝑗 𝑘 (𝒓2; 𝒑2)𝑞𝑘,2𝑤
+ 𝑞𝑖,2𝑣+1𝑃𝑖𝑙 (𝒓1; 𝒑1)𝑞𝑙,2𝑤𝑞 𝑗 ,2𝑣𝑃 𝑗 𝑘 (𝒓2; 𝒑2)𝑞𝑘,2𝑤+1

]
.

This leads to the same contribution as the first summand. Then use symmetry under the interchange
of the two particles. In summary,

𝑑2(𝒓1, 𝒓2; 𝒑1, 𝒑2)
=
∑︁
𝑣𝑤

𝑎𝑣𝑎𝑤
[
𝜒2𝑣,2𝑤 (𝒓1; 𝒑1)𝜒2𝑣+1,2𝑤+1(𝒓2; 𝒑2) − 𝜒2𝑣,2𝑤+1(𝒓1; 𝒑1)𝜒2𝑣+1,2𝑤 (𝒓2; 𝒑2)

− 𝜒2𝑣+1,2𝑤 (𝒓1; 𝒑1)𝜒2𝑣,2𝑤+1(𝒓2; 𝒑2) + 𝜒2𝑣+1,2𝑤+1(𝒓1; 𝒑1)𝜒2𝑣,2𝑤 (𝒓2; 𝒑2)
]
. (9)

7


The reduced 1-body phase space (spinless) quasidensity for the triplet is obtained, as before,

𝑑1(𝒓; 𝒑) = 2
∫

𝑑2(𝒓, 𝒓2; 𝒑, 𝒑2) 𝑑𝒓2 𝑑 𝒑2 = 2
∑︁
𝑤

𝑎2
𝑤 [𝜒(2𝑤,2𝑤) (𝒓; 𝒑) + 𝜒(2𝑤+1,2𝑤+1) (𝒓; 𝒑)] . (10)

Notice that in the previous equation each occupation number 𝑛𝑖 := 2𝑎2
𝑖

appears twice. This is a
consequence of the Pauli exclusion principle.

Unlike the singlet case, there is no sign rule to be deciphered here. Instead there are the
ambiguities:

𝜒2𝑤,2𝑤 = 𝜒′2𝑤,2𝑤 cos2 𝜃𝑤 − (𝜒′2𝑤,2𝑤+1 + 𝜒′2𝑤+1,2𝑤) sin 𝜃𝑤 cos 𝜃𝑤 + 𝜒′2𝑤+1,2𝑤+1 sin2 𝜃𝑤,

𝜒2𝑤+1,2𝑤+1 = 𝜒′2𝑤,2𝑤 sin2 𝜃𝑤 + (𝜒′2𝑤,2𝑤+1 + 𝜒′2𝑤+1,2𝑤) sin 𝜃𝑤 cos 𝜃𝑤 + 𝜒′2𝑤+1,2𝑤+1 cos2 𝜃𝑤 .

They clearly leave the form (10) untouched. We see here the action of 𝑆𝑂 (2) on each invariant
block. One may choose the angles as to maximize their overlap with the leading natural orbitals for
the ground state, as done in the seminal paper by Löwdin and Shull [8]. We omit that. Let us define

𝐴𝑤 :=
(
cos 𝜃𝑤 − sin 𝜃𝑤
sin 𝜃𝑤 cos 𝜃𝑤

)
.

The above transformation can be construed as

𝜒 = (𝐴𝑣 ⊗ 𝐴𝑤) 𝜒′ =
©­­­«
cos 𝜃𝑣 cos 𝜃𝑤 − cos 𝜃𝑣 sin 𝜃𝑤 − sin 𝜃𝑣 cos 𝜃𝑤 sin 𝜃𝑣 sin 𝜃𝑤
cos 𝜃𝑣 sin 𝜃𝑤 cos 𝜃𝑣 cos 𝜃𝑤 − sin 𝜃𝑣 sin 𝜃𝑤 − sin 𝜃𝑣 cos 𝜃𝑤
sin 𝜃𝑣 cos 𝜃𝑤 − sin 𝜃𝑣 sin 𝜃𝑤 cos 𝜃𝑣 cos 𝜃𝑤 − cos 𝜃𝑣 sin 𝜃𝑤
sin 𝜃𝑣 sin 𝜃𝑤 sin 𝜃𝑣 cos 𝜃𝑤 cos 𝜃𝑣 sin 𝜃𝑤 cos 𝜃𝑣 cos 𝜃𝑤

ª®®®¬ 𝜒
′,

with

𝜒 :=
©­­­«

𝜒2𝑣,2𝑤
𝜒2𝑣,2𝑤+1
𝜒2𝑣+1,2𝑤
𝜒2𝑣+1,2𝑤+1

ª®®®¬ , and similarly for 𝜒′,

in the case 𝑣 = 𝑤.
To similarly examine the symmetry of expression (9), again one does not have to contend with

the whole tensor product matrix, since most contributions vanish. As regards the sum in (9), one
can write in compressed form:

𝜒𝜒 =

©­­­«
cos2 𝜃𝑣 cos2 𝜃𝑤 − cos2 𝜃𝑣 sin2 𝜃𝑤 − sin2 𝜃𝑣 cos2 𝜃𝑤 sin2 𝜃𝑣 sin2 𝜃𝑤
− cos2 𝜃𝑣 sin2 𝜃𝑤 cos2 𝜃𝑣 cos2 𝜃𝑤 sin2 𝜃𝑣 sin2 𝜃𝑤 − sin2 𝜃𝑣 cos2 𝜃𝑤
− sin2 𝜃𝑣 cos2 𝜃𝑤 sin2 𝜃𝑣 sin2 𝜃𝑤 cos2 𝜃𝑣 cos2 𝜃𝑤 − cos2 𝜃𝑣 sin2 𝜃𝑤
sin2 𝜃𝑣 sin2 𝜃𝑤 − sin2 𝜃𝑣 cos2 𝜃𝑤 − cos2 𝜃𝑣 sin2 𝜃𝑤 cos2 𝜃𝑣 cos2 𝜃𝑤

ª®®®¬ 𝜒
′𝜒′,

with

𝜒𝜒 :=
©­­­«
𝜒2𝑣,2𝑤 (𝒓1; 𝒑1) 𝜒2𝑣+1,2𝑤+1(𝒓2; 𝒑2)
𝜒2𝑣,2𝑤+1(𝒓1; 𝒑1) 𝜒2𝑣+1,2𝑤 (𝒓2; 𝒑2)
𝜒2𝑣+1,2𝑤 (𝒓1; 𝒑1) 𝜒2𝑣,2𝑤+1(𝒓2; 𝒑2)
𝜒2𝑣+1,2𝑤+1(𝒓1; 𝒑1) 𝜒2𝑣,2𝑤 (𝒓2; 𝒑2)

ª®®®¬ ; and similarly for 𝜒′𝜒′.

One verifies that (9) is invariant under this set of transformations.

8


5 Lowest triplet state of harmonium
The energy spectrum for harmonium is obviously (ℕ + 3

2 )𝜔 + (ℕ + 3
2 )𝜇. Since 𝜇 < 𝜔, the energy

of the first excited states is 𝐸fs = (3𝜔 + 5𝜇)/2. For our present purposes, it is enough to choose an
intracule excitation state along the 𝑥-axis (say). The corresponding 2-quasidensity is given by:

2
𝜋6 exp

(
−2𝐻𝑅

𝜔

)
exp

(
−2𝐻𝑟

𝜇

) ( (𝑝1𝑥 − 𝑝2𝑥)2 + 𝜇2(𝑥2
1 − 𝑥2

2)
2

𝜇
− 1

2

)
. (11)

Henceforth we work in the chosen nontrivial mode, since the problem factorizes completely. By
integrating one set of variables, the reduced one-body spinless quasidensity is obtained, after some
work:

𝑑1(𝑟; 𝑝) = 2
∫

𝑑2(𝑟, 𝑟2; 𝑝, 𝑝2) 𝑑𝑟2 𝑑𝑝2 =
2
𝜋

(2√𝜔𝜇

𝜔 + 𝜇

)3
𝑒
− 2𝜔𝜇

𝜔+𝜇 𝑟
2− 2

𝜔+𝜇 𝑝
2
(
𝜔𝑟2 + 1

𝜔
𝑝2
)
. (12)

The marginals of 𝑑1 give the electronic density and momentum density:

𝜌(𝑟) =
∫

𝑑1(𝑟; 𝑝) 𝑑𝑝 =
2
𝜋

(2√𝜔𝜇

𝜔 + 𝜇

)3
𝑒
− 2𝜔𝜇

𝜔+𝜇 𝑟
2
∫

𝑒
− 2

𝜔+𝜇 𝑝
2
(
𝜔𝑟2 + 1

𝜔
𝑝2
)
𝑑𝑝

=
2
𝜋

(2√𝜔𝜇

𝜔 + 𝜇

)3
𝑒
− 2𝜔𝜇

𝜔+𝜇 𝑟
2
(
𝜋(𝜔 + 𝜇)

2

)1/2 (
𝜔𝑟2 + 𝜔 + 𝜇

4𝜔

)
,

𝜋(𝑝) =
∫

𝑑1(𝑟; 𝑝) 𝑑𝑟 =
2
𝜋

(2√𝜔𝜇

𝜔 + 𝜇

)3
𝑒
− 2

𝜔+𝜇 𝑝
2
∫

𝑒
− 2𝜔𝜇

𝜔+𝜇 𝑟
2
(
𝜔𝑟2 + 1

𝜔
𝑝2
)
𝑑𝑟

=
2
𝜋

(2√𝜔𝜇

𝜔 + 𝜇

)3
𝑒
− 2

𝜔+𝜇 𝑝
2
(
𝜋(𝜔 + 𝜇)

2𝜔𝜇

)1/2 (
𝜔 + 𝜇

4𝜇
+ 1
𝜔
𝑝2
)
.

Finally, as expected, we get∫
𝜋(𝑝) 𝑑𝑝 =

∫
𝜌(𝑟) 𝑑𝑟 = 2

𝜋

(2√𝜔𝜇

𝜔 + 𝜇

)3 (
𝜋(𝜔 + 𝜇)

2

)1/2 (
𝜋(𝜔 + 𝜇)

2𝜔𝜇

)1/2 (
𝜔 + 𝜇

4𝜇
+ 𝜔 + 𝜇

4𝜔

)
= 2.

From the viewpoint of WDFT, the most interesting part of the energy corresponds to the
interelectronic repulsion of this first excited state 𝐸2fs. The 1-body Hamiltonian is given by
ℎ(𝑟, 𝑝) = 𝑝2/2 + 𝜔2𝑟2/2. It is a simple exercise to obtain the 1-body energy 𝐸1fs by integrating
expression (12) with this observable:

𝐸1fs =
𝜔

2
+ 3(𝜇2 + 𝜔2)

4𝜇
.

The interelectronic potential in (3) is (𝜇2 − 𝜔2)𝑟2
12/4, so to obtain the repulsion energy 𝐸2fs, one

has just to integrate expression (11) with this observable:

𝐸2fs =

∫
2
𝜋2 exp

(
−2𝐻𝑅

𝜔

)
exp

(
−2𝐻𝑟

𝜇

) [
2𝐻𝑟

𝜇
− 1

2

]
𝜇2 − 𝜔2

4
𝑟2

12 𝑑𝑅 𝑑𝑟 𝑑𝑃 𝑑𝑝

=
1
𝜋
(𝜇2 − 𝜔2)

∫
exp

(
−2𝐻𝑟

𝜇

) [
𝑟2𝑝2

𝜇
+ 𝜇𝑟4 − 𝑟2

2

]
𝑑𝑟 𝑑𝑝 =

3
4
𝜇2 − 𝜔2

𝜇
,

9


which is 3 times the interelectronic repulsion energy for the corresponding mode of the singlet [23].
This is not surprising, since, in the triplet configuration the electrons tend to be mutually farther
apart than in the singlet.1

6 Spectral analysis of the 1-body triplet state
In order to determine the occupation numbers of this system, first we have to find the good
coordinates. Let us perform the transformation

(𝑄, 𝑃) :=
(
(𝜔𝜇)1/4𝑟, (𝜔𝜇)−1/4𝑝

)
; or, in shorthand, 𝑈 = 𝑆𝑢,

where 𝑆 is symplectic and 𝑢 = (𝑟, 𝑝). We may also write 𝜗 := arctan(𝑃/𝑄), so that

𝑃 = 𝑈 sin 𝜗 and 𝑄 = 𝑈 cos 𝜗.

Recalling 2√𝜔𝜇/(𝜔+ 𝜇) = (1− 𝑡2)/(1+ 𝑡2) from (7), the 1-quasidensity (12) takes the simple form:

𝑑1(𝑈, 𝜗) := 𝑑1(𝑢(𝑈, 𝜗)) = 2(1 − 𝑡2)3

𝜋(1 + 𝑡2)3 𝑒−(1−𝑡
2)𝑈2/(1+𝑡2)𝑈2

(
1 + 𝑡

1 − 𝑡
cos2 𝜗 + 1 − 𝑡

1 + 𝑡
sin2 𝜗

)
=

2(1 − 𝑡2)3

𝜋(1 + 𝑡2)3 𝑒−(1−𝑡
2)𝑈2/(1+𝑡2)𝑈2

(
1 + 𝑡2

1 − 𝑡2
+ 2𝑡

1 − 𝑡2
cos 2𝜗

)
.

The one-body quasidensity may be expanded as follows:

𝑑1(𝑈, 𝜗) =
∑︁
𝑟𝑠

𝑓𝑟𝑠 (𝑈, 𝜗) 𝑑𝑟𝑠 where 𝑑𝑟𝑠 = 2𝜋
∫

𝑑1(𝑈, 𝜗) 𝑓 ∗𝑟𝑠 (𝑈, 𝜗)𝑈 𝑑𝑈 𝑑𝜗.

Then, for 𝑟 ⩾ 𝑠,

2𝜋
∫

𝑓 ∗𝑟𝑠 (𝑈, 𝜗) 𝑑1(𝑈, 𝜗)𝑈 𝑑𝑈 𝑑𝜗

=
4(1 − 𝑡2)3

𝜋(1 + 𝑡2)3 (−1)𝑠
√
𝑠!

√
𝑟!

∫ ∞

0
𝑒−(1−𝑡

2)𝑈2/(1+𝑡2)𝑒−𝑈
2 (2𝑈2) (𝑟−𝑠)/2𝐿𝑟−𝑠

𝑠 (2𝑈2)𝑈3 𝑑𝑈

×
∫ 𝜋

−𝜋
𝑒𝑖(𝑟−𝑠)𝜗

[
1 + 𝑡2

1 − 𝑡2
+ 2𝑡

1 − 𝑡2
cos 2𝜗

]
𝑑𝜗

=
4(1 − 𝑡2)3

𝜋(1 + 𝑡2)3 (−1)𝑠
√
𝑠!

√
𝑟!

∫ ∞

0
𝑒−(1−𝑡

2)𝑈2/(1+𝑡2)𝑒−𝑈
2 (2𝑈2) (𝑟−𝑠)/2𝐿𝑟−𝑠

𝑠 (2𝑈2)𝑈3 𝑑𝑈

× 𝜋

[
2(1 + 𝑡2)

1 − 𝑡2
𝛿𝑠𝑟 +

2𝑡
1 − 𝑡2

(𝛿𝑠+2
𝑟 + 𝛿𝑠−2

𝑟 )
]
,

so that

𝑑1(𝑈, 𝜗) =
∑︁
𝑠

𝑑𝑠𝑠 (𝑡) 𝑓𝑠𝑠 (𝑈, 𝜗) + 𝑑𝑠+2,𝑠 (𝑡) 𝑓𝑠+2,𝑠 (𝑈, 𝜗) + 𝑑𝑠,𝑠+2(𝑡) 𝑓𝑠,𝑠+2(𝑈, 𝜗),

where actually 𝑑𝑠+2,𝑠 = 𝑑𝑠,𝑠+2.

1Interestingly, (12) is a non-Gaussian Wigner function taking only positive values. This prompts two remarks. First,
in consonance with common wisdom [30, 31], it is confirmed that as of itself 𝑑1 is a nearly classical state. Second,
there are mathematical recipes that produce such positive-valued Wigner functions representing mixed states [32]. It
would be good to know whether or not (12) can be obtained as such an output.

10


Using the standard Mellin transform [33, 34]:∫ ∞

0
𝑥𝛼−1 𝑒−𝑝𝑥𝐿𝜆

𝑛 (𝑐𝑥) 𝑑𝑥 =
Γ(𝛼)
𝑝𝛼

𝑃
(𝜆,𝛼−𝜆−𝑛−1)
𝑛

(
1 − 2𝑐

𝑝

)
=
Γ(𝛼)
𝑝𝛼

(𝜆 + 1)𝑛
𝑛! 2𝐹1

(
−𝑛, 𝛼
𝜆 + 1 ;

𝑐

𝑝

)
,

we obtain by fairly easy manipulations,

𝑑𝑠𝑠 (𝑡) = (1 − 𝑡2)2 (𝑠 𝑡2𝑠−2 + (1 + 𝑠) 𝑡2𝑠
)
;

𝑑𝑠,𝑠+2(𝑡) = (1 − 𝑡2)2
√︁
(𝑠 + 1) (𝑠 + 2) 𝑡2𝑠+1.

This means that, to find the occupation numbers, one has to diagonalize a symmetric pentadiagonal
matrix:

𝐷 = (1 − 𝑡2)2

©­­­­­­­­­­«

1 0 𝛼0𝑡 0 0 0 · · ·
0 1 + 2𝑡2 0 𝛼1𝑡

3 0 0 · · ·
𝛼0𝑡 0 2𝑡2 + 3𝑡4 0 𝛼2𝑡

5 0 · · ·
0 𝛼1𝑡

3 0 3𝑡4 + 4𝑡6 0 𝛼3𝑡
7 · · ·

0 0 𝛼2𝑡
5 0 4𝑡6 + 5𝑡8 0 · · ·

0 0 0 𝛼3𝑡
7 0 5𝑡8 + 6𝑡10 · · ·

...
...

...
...

...
...

. . .

ª®®®®®®®®®®¬
, (13)

where 𝛼𝑠 :=
√︁
(𝑠 + 1) (𝑠 + 2) .

It is readily checked that the trace of this matrix is 2, as it should be. Its eigenspaces split into
two parts: ℓ2 = 𝑉1 ⊕ 𝑉2, where 𝑉1 = { 𝒙 : all 𝑥2𝑛 = 0 } and 𝑉2 = { 𝒙 : all 𝑥2𝑛+1 = 0 }. They
correspond respectively to the matrices

𝐷even = (1 − 𝑡2)2

©­­­­­­­­«

1 𝛼0𝑡 0 0 0 · · ·
𝛼0𝑡 2𝑡2 + 3𝑡4 𝛼2𝑡

5 0 0 · · ·
0 𝛼2𝑡

5 4𝑡6 + 5𝑡8 𝛼4𝑡
9 0 · · ·

0 0 𝛼4𝑡
9 6𝑡10 + 7𝑡12 𝛼6𝑡

13 · · ·
0 0 0 𝛼6𝑡

13 8𝑡14 + 9𝑡16 · · ·
...

...
...

...
...

. . .

ª®®®®®®®®¬
and

𝐷odd = (1 − 𝑡2)2

©­­­­­­­­«

1 + 2𝑡2 𝛼1𝑡
3 0 0 0 · · ·

𝛼1𝑡
3 3𝑡4 + 4𝑡6 𝛼3𝑡

7 0 0 · · ·
0 𝛼3𝑡

7 5𝑡8 + 6𝑡10 𝛼5𝑡
11 0 · · ·

0 0 𝛼5𝑡
11 7𝑡12 + 8𝑡14 𝛼7𝑡

15 · · ·
0 0 0 𝛼7𝑡

15 9𝑡16 + 10𝑡18 · · ·
...

...
...

...
...

. . .

ª®®®®®®®®¬
.

It is easily checked that these matrices have the same set of eigenvalues, as they should, since the
occupation numbers must appear twice.

11


As was shown in Section 3, there is a skewsymmetric matrix 𝐶 such that 𝐷 = 𝐶𝑡𝐶. This matrix
is tridiagonal, and is the sum of two skew-symmetric matrices whose diagonalization is trivial:

𝐶 = (1 − 𝑡2)

©­­­­­­­­«

0 −1 0 0 0 · · ·
1 0 0 0 0 · · ·
0 0 0 −

√
3𝑡2 0 · · ·

0 0
√

3𝑡2 0 0 · · ·
0 0 0 0 0 · · ·
...

...
...

...
...
. . .

ª®®®®®®®®¬
+ (1 − 𝑡2)

©­­­­­­­­­«

0 0 0 0 0 · · ·
0 0

√
2𝑡 0 0 · · ·

0 −
√

2𝑡 0 0 0 · · ·
0 0 0 0

√
4𝑡3 · · ·

0 0 0 −
√

4𝑡3 0 · · ·
...

...
...

...
...

. . .

ª®®®®®®®®®¬
=: 𝐴 + 𝐵.

Also, 𝐷 is the sum of two Hermitian matrices, namely 𝐴𝑡𝐴+𝐵𝑡𝐵, which is diagonal, and 𝐴𝑡𝐵+𝐵𝑡𝐴.
One is reminded here of the Weyl problem: given two 𝑛 × 𝑛 Hermitian matrices 𝐴, 𝐵 whose

spectra are known, what could the spectrum of their sum 𝐶 := 𝐴 + 𝐵 be? Some facts are clear: with
an obvious notation for the eigenvalues, these must satisfy

𝑐1 + · · · + 𝑐𝑛 = 𝑎1 + · · · + 𝑎𝑛 + 𝑏1 + · · · + 𝑏𝑛; 𝑐1 ⩽ 𝑎1 + 𝑏1;

less clear, but also true, are
𝑐2 ⩽ 𝑎1 + 𝑏2; 𝑐2 ⩽ 𝑎2 + 𝑏1;

and so on. The conditions written above are already optimal for 𝑛 = 2. The necessary constraints
are all linear homogeneous inequalities, bounding convex polyhedra. Horn made a conjecture for
the general form of such inequalities, which was eventually proved [35].

The pure-state 𝑁-representability problem in quantum chemistry (or “quantum marginal prob-
lem”, in the jargon of information theory) should be considered as solved, after the work by
Klyachko [36, 37]. It is of the same type and answered by similar inequalities. Both questions
reduce to finding moment polyhedra for coadjoint orbits of unitary groups (associated to pertinent
Hilbert spaces), which are computed by Duistermaat–Heckman measures [38]. A very readable and
up-to-date account of all this is [39]. The Hilbert spaces considered are finite-dimensional. How-
ever, the results are valid for finite-rank approximations in the chemical context, and the patterns of
the inequalities extend in a rather obvious way. Thus it is scarcely surprising that the Weyl problem
surfaces in this simple instance. We leave for the future consideration of the moment polytopes
for the occupation numbers,2 and choose in this paper a direct approach to the eigenpair problem,
completed by numerical analysis.

The matrices 𝐷even and 𝐷odd are tridiagonal symmetric real matrices. The general eigenvalue
problem for a matrix 𝑇 of this kind reduces to solving the following set of recurrence equations:

©­­­­­­­­«

𝑑0 𝑡1 0 0 · · ·
𝑡1 𝑑1 𝑡2 0 · · ·
0 𝑡2 𝑑2 𝑡3 · · ·
0 0 𝑡3 𝑑3 · · ·
0 0 0 𝑡4 · · ·
...

...
...

...
. . .

ª®®®®®®®®¬

©­­­­­­­­«

𝜙0(𝑛𝑟)
𝜙1(𝑛𝑟)
𝜙2(𝑛𝑟)
𝜙3(𝑛𝑟)
𝜙4(𝑛𝑟)

...

ª®®®®®®®®¬
= 𝑛𝑟

©­­­­­­­­«

𝜙0(𝑛𝑟)
𝜙1(𝑛𝑟)
𝜙2(𝑛𝑟)
𝜙3(𝑛𝑟)
𝜙4(𝑛𝑟)

...

ª®®®®®®®®¬
,

2The number of their extremal edges grows very quickly with 𝑁 and the rank; this makes for precision, but also for
strenuous work.

12


where 𝑛𝑟 is an eigenvalue. The general solution is completely given in terms of the occupation
numbers, by the following formula [40, Sect. 5.48]:

𝜙𝑚 (𝜆) =
𝜙0(𝜆)

𝑡1𝑡2 . . . 𝑡𝑚
det[𝜆𝐼 − 𝑇]𝑚𝑚, for each 𝑚 ⩾ 1,

where [𝜆𝐼 − 𝑇]𝑚𝑚 is the upper left 𝑚 × 𝑚 submatrix of (𝜆𝐼 − 𝑇), and 𝜙0(𝜆) ≠ 0 is chosen so as to
normalize the eigenvector.

This result implies that 𝑇 = 𝑄𝐷𝑄𝑡 , where 𝑑𝑖 𝑗 = 𝑛𝑖𝛿
𝑖
𝑗

is the diagonal matrix whose entries are
the eigenvalues and 𝑞𝑖 𝑗 = 𝜙𝑖 (𝑛 𝑗 ). Since 𝑄𝑄𝑡 = 𝑄𝑡𝑄 = 1, the following orthogonality relations
hold:

∞∑︁
𝑟=0

𝜙𝑚 (𝑛𝑟)𝜙𝑙 (𝑛𝑟) = 𝛿𝑚𝑙 ,

∞∑︁
𝑚=0

𝜙𝑚 (𝑛𝑟)𝜙𝑚 (𝑛𝑠) = 𝛿𝑟𝑠,

∞∑︁
𝑟=0

𝑛𝑟 𝜙𝑚 (𝑛𝑟)𝜙𝑙 (𝑛𝑟) = 𝑑𝑚𝛿
𝑚
𝑙 + 𝑡𝑚𝛿

𝑚−1
𝑙 .

In summary, for 𝑑1 we obtain

𝑑1(·) =
∑︁
𝑟

𝑛𝑟

[ ∞∑︁
𝑖=0

𝑓2𝑖,2𝑖 (·) 𝜙2
even,𝑖 (𝑛𝑟) +

∞∑︁
𝑖=0

( 𝑓2𝑖,2𝑖+2 + 𝑓2𝑖+2,2𝑖) (·) 𝜙even,𝑖 (𝑛𝑟)𝜙even,𝑖+1(𝑛𝑟)

+
∞∑︁
𝑖=0

𝑓2𝑖+1,2𝑖+1(·) 𝜙2
odd,𝑖 (𝑛𝑟) +

∞∑︁
𝑖=0

( 𝑓2𝑖+1,2𝑖+3 + 𝑓2𝑖+3,2𝑖+1) (·) 𝜙odd,𝑖 (𝑛𝑟)𝜙odd,𝑖+1(𝑛𝑟)
]
.

Here 𝑛𝑟 depends solely on the parameter 𝑡 of (7).

7 Numerical analysis of the occupation numbers
As advertised, to find the 𝑛𝑟 we fall back on numerical computation. Figure 1 shows the behavior
of the rank-eight approximation of the eigenvalues, as 𝑡 is varied. Note that the first eigenvalue is
very close to 1 in the neighborhood of 𝑡 = 0, while the others are very small. As the value of 𝑡 rises,
the first eigenvalue begins to decrease and the others rise for a while. In the neighborhood of 𝑡 = 1
all eigenvalues approach zero.

Note that 𝑡 is a very nonlinear parameter: although 𝑡 ∼ 𝛿/8𝑘 for small 𝛿, the value 𝑡 = 1/2
means 𝜇/𝜔 = 1/9 or 𝛿/𝑘 = 80/81. This shows that, unless 𝛿 is pretty close to the dissociation
value, the harmonium triplet is not badly described by a Hartree–Fock state. Whenever 𝑡 ≲ 0.6, that
is, 𝛿/𝑘 ≲ 255/256, the first two occupation numbers contain almost all the physical information for
the system.

Also, one we can show that whenever 𝑡 ≲ 0.5, a good approximation to the five first occupation
numbers is

𝜆1 ≈ 1 − 3𝑡4 + 8𝑡6, 𝜆2 ≈ 3𝑡4 − 8𝑡6, 𝜆3 ≈ 5𝑡8, 𝜆4 ≈ 7𝑡12 and 𝜆5 ≈ 9𝑡16.

Figure 2 compares the behavior of the first two eigenvalues for the singlet and triplet states of
harmonium. In this sense, the Hartree–Fock approximation works better in the triplet case than for

13


Figure 1: First six eigenvalues of the matrix 𝐷even.

Figure 2: First and second occupation numbers of the ground state and of the first excited state.

the singlet. Around 𝑡 = 0.4 the second approximated occupation number for the latter is above 0.13,
and for the former is below 0.052. The same behaviour was also observed in the toy model studied
in [42]. This does not mean, however, that correlation is always weaker in the triplet state – see the
next section.

8 Spatial entropy and correlation energies
We move towards the comparison of the triplet system with the singlet system in regard to disorder
(suppressing the spin variables). To measure this, a useful quantity is the linear entropy 𝑠 associated
to the 1-body function:

𝑠 = 1 − Π(𝑑1),
whereΠ(𝑑1) is the purity of the system – see below.3 Mathematically, the quantity 𝑠 is a lower bound
for the Jaynes entropy, which has been used to quantify the entanglement between one particle and
the other 𝑁−1 particles of the system [42], and proposed as a handle on the correlation energies [43].

3Truth to be told, the notion of entropy native to the Wigner quasiprobability approach is the one discussed in [41].
We put aside the question of its eventual usefulness here.

14


In this paper the singlet has been modelled in such a way that, for each one-dimensional mode:

Πgs,1(𝑑1) =
∫

𝑑2
1 (𝑟; 𝑝) 𝑑𝑟 𝑑𝑝 =

∑︁
𝑖

𝑛2
𝑖 .

Instead, for the triplet one should take for the excited mode:

Πfs,𝑥 (𝑑1) =
1
2

∫
𝑑2

1 (𝑟𝑥; 𝑝𝑥) 𝑑𝑟𝑥 𝑑𝑝𝑥 =
∑︁
𝑖

𝑛2
𝑖 .

This second definition is natural in that correlations due solely to the antisymmetric character of the
wave function do not contribute to the entanglement of the system [18, 44, 45]. This ensures that
the entropy for a 1-body function of the Hartree–Fock type is zero.

In the singlet case, the occupation numbers are equal to (1 − 𝑡2) 𝑡2𝑖. Thus, the purity of this
system is easily computable, to wit, Πgs,1(𝑑1) = (1 − 𝑡2)/(1 + 𝑡2) for each mode. This quantity
coincides with the quotient of the geometric and arithmetic means of the frequencies, that is,
Πgs,1 = 2√𝜔𝜇/(𝜔 + 𝜇). For 𝑛 modes one just takes the 𝑛th power [17]. Moreover, for small values
of the coupling 𝛿, we obtain

𝑠gs,1 ∼ 1
32

𝛿2

𝜔4 , (14)

which for this approximation is exactly the absolute value of the (dimensionless) correlation en-
ergy [26]. This appears to vindicate the contention of [43]. (Actually, for the singlet it is not difficult
to compute the Jaynes entropy, given by

−
∑︁
𝑖

𝑛𝑖 log 𝑛𝑖 = − log(1 − 𝑡2) − 𝑡2 log 𝑡2

1 − 𝑡2
.

This was done by Srednicki [20] some time ago.)

For the triplet state, we have to compute Tr(𝑑2
1) for the matrix given in (13). Since

𝑑2
1 = (1 − 𝑡2)4

©­­­­«
1 + 𝛼2

0𝑡
2 0 𝛼0𝑡 (1 + 2𝑡2 + 3𝑡4) · · ·

0 (1 + 2𝑡2)2 + 𝛼2
1𝑡

6 0 · · ·
𝛼0𝑡 (1 + 2𝑡2 + 3𝑡4) 0 𝛼2

0𝑡
2 + (2𝑡2 + 3𝑡4)2 + 𝛼2

2𝑡
10 · · ·

...
...

...
. . .

ª®®®®¬
,

we get

Tr(𝑑2
1) = (1 − 𝑡2)4

[
4

∞∑︁
𝑖=0

𝛼2
𝑖 𝑡

2(2𝑖+1) + 2
∞∑︁
𝑖=1

𝑖2 𝑡4(𝑖−1)
]

= 2(1 − 𝑡2)4
∞∑︁
𝑖=1

[
2𝑖(𝑖 + 1) 𝑡2(2𝑖−1) + 𝑖2 𝑡4(𝑖−1)

]
=

2(1 − 𝑡2)
1 + 𝑡2

[
1 + 2𝑡2

(1 + 𝑡2)2

]
,

after some calculation. So the purity of the first excited mode is

Πfs,𝑥 =
1 − 𝑡2

1 + 𝑡2

[
1 + 2𝑡2

(1 + 𝑡2)2

]
= Πgs,1

[
1 + 2𝑡2

(1 + 𝑡2)2

]
=

2√𝜔𝜇

𝜔 + 𝜇

(
1 + 1

2

(𝜔 − 𝜇

𝜔 + 𝜇

)2
)
.

15


Since the other two modes contribute with two ground state factors, the total purity can be written
as Πfs = Πfs,𝑥Πgs,𝑦Πgs,𝑧. For the purity parameter, one obtains finally

𝑠gs = 1 −
(
1 − 𝑡2

1 + 𝑡2

)3
and 𝑠fs = 1 −

(
1 − 𝑡2

1 + 𝑡2

)3 [
1 + 2𝑡2

(1 + 𝑡2)2

]
= 𝑠gs −

2𝑡2(1 − 𝑡2)3

(1 + 𝑡2)5 .

In conclusion, 𝑠fs ⩽ 𝑠gs.
At long last, we may go back to Moshinsky’s starting point, the assessment of electron corre-

lation, only now for the excited state. The Hartree–Fock approximation for the relevant mode, in
view of (8), is of the form

𝑊HF(𝑟1, 𝑟2; 𝑝1, 𝑝2) =
1
2
[
𝑊00(𝑟1; 𝑝1)𝑊11(𝑟2; 𝑝2) −𝑊01(𝑟1; 𝑝1)𝑊10(𝑟2; 𝑝2)

−𝑊10(𝑟1; 𝑝1)𝑊01(𝑟2; 𝑝2) +𝑊11(𝑟1; 𝑝1)𝑊00(𝑟2; 𝑝2)
]
,

where

𝑊00(𝑟; 𝑝) =
1
𝜋
𝑒−𝜂𝑟

2−𝑝2/𝜂, 𝑊11(𝑟; 𝑝) =
2
𝜋
𝑒−𝜂𝑟

2−𝑝2/𝜂 (𝜂𝑟2 + 𝑝2/𝜂 − 1
2 ),

with their corresponding interferences. Remember that
∫
𝑊𝑖 𝑗 𝑑𝑟 𝑑𝑝 = 𝛿𝑖 𝑗 . In intracule-extracule

coordinates:
𝑊HF(𝑅, 𝑟; 𝑃, 𝑝) =

2
𝜋2

(
𝜂𝑟2 + 𝑝2/𝜂 − 1

2
)
𝑒−𝜂𝑅

2−𝑃2/𝜂−𝜂𝑟2−𝑝2/𝜂 .

The parameter 𝜂 is determined by minimization. The mean value of the energy predicted by
this function is:

𝐸HF =
1
2

∫
(𝑝2 + 𝜔2𝑟2) [𝑊00(𝑟; 𝑝) +𝑊11(𝑟; 𝑝)] 𝑑𝑟 𝑑𝑝 − 𝛿

4

∫
(𝑟1 − 𝑟2)2 𝑊HF(1, 2) 𝑑1 𝑑2

=

(
𝜂 + 𝜔2

𝜂

)
− 3𝛿

4𝜂
= 𝜂 + 𝜔2 + 3𝜇2

4𝜂
.

The minimum 𝑑𝐸/𝑑𝜂 = 0 occurs when 𝜂 = 1
2

√︁
𝜔2 + 3𝜇2. Therefore, the energy predicted by

Hartree–Fock is
√︁
𝜔2 + 3𝜇2. Thus, the “correlation energy” for the lowest excited state of harmo-

nium is
𝐸c,fs = 𝐸fs − 𝐸HF =

3𝜔 + 5𝜇
2

−
√︃
𝜔2 + 3𝜇2 − 2

√︃
(𝜔2 + 𝜇2)/2 ∼ − 7

64
𝛿2

𝜔3 .

Thus, the relative correlation energies are

Efs :=
|𝐸c,fs |
𝐸fs

∼ 7
256

𝛿2

𝜔4 and Egs :=
|𝐸c,gs |
𝐸gs

∼ 1
32

𝛿2

𝜔4 .

Both quantities are related by a factor of 7/8. For this approximation, as one would have expected,
Efs ⩽ Egs.

Figure 3 shows the exact dependence of the relative correlation energy for both systems as a
function of 𝛿/𝜔2. The relative correlation energy for the singlet is greater than for the triplet, just as
the purity parameter for the singlet is greater than the one for the triplet. At 𝛿/𝜔2 = 0.67 the relation
between these two quantities changes and the relative correlation energy for the triplet is greater than

16


Figure 3: Relative correlation energy of the singlet and of the triplet excited mode. As expected,
the relative correlation energy for the singlet is greater than for the triplet for small values of the
coupling. At 𝛿/𝜔2 ∼ 0.67 the order is inverted.

Figure 4: Moshinsky’s hole for the triplet: (𝜌(𝑟) − 𝜌HF(𝑟))/𝜔1/2 as a function of 𝜔1/2𝑟.

for the singlet. Note however that the entropy depends only the behavior of the occupation numbers,
while the correlation energy has to do with the natural orbitals as well. Such a nice proportionality
as (14) fails for the triplet state.

Finally, Figure 4 shows the difference between the exact profile 1-density and the Hartree–
Fock profile 1-density for the harmonium triplet, 𝜌HF(𝑟) :=

∫
𝑊HF(𝑟, 𝑟2; 𝑝1, 𝑝2) 𝑑𝑝1 𝑑𝑟2 𝑑𝑝2. This

description goes back to the Coulson–Neilson classic paper [46] on the helium Coulomb system.
The “Moshinsky’s hole” observed in the neighborhood of 𝑟 = 0 graphically shows the Hartree–
Fock underestimation of the mean distance between the fermions, for the excited configuration of
harmonium as well.

9 Conclusion
From the very beginning of quantum mechanics, the fundamental state of harmonium has provided
a useful playground for learning about such questions as correlation energy, entanglement or hole
entropy (including black hole entropy). Here, for the first time, we rather exhaustively analyze the
(spin triplet) first excited configuration of harmonium, particularly the behaviour of its occupation

17


numbers and natural orbitals. This is a different chemical species altogether, due to the antisymmetric
character of the orbital wave function. When exactly reconstructing à la Löwdin–Shull–Kutzelnigg
the two-body density as a functional of the one-body density, instead of the sign dilemma (already
solved by two of us) for the lowest-energy state, we find, as expected on general grounds, an
ambiguity in the choice of natural orbitals.

Also as expected, in the triplet case the first occupation number plays a more dominant role than
for the singlet, up to fairly high values of the coupling parameter, 𝑡 ≲ 0.4. Thus, within this range,
modeling the excited configuration as a Hartree–Fock state introduces a lower error than doing so
for the ground state. In parallel, the linear entropy of the first excited configuration is lower than
that of the ground state, and the relative correlation energy for the excited state stays below that of
the ground state for such values of the coupling. The order reverses at higher values of 𝑡.

Acknowledgments

JMGB thanks the Zentrum für interdisziplinäre Forschung (ZiF) at Bielefeld, in whose welcoming
atmosphere this paper received its finishing touches. CMBR and JMGB are grateful to Andrés F.
Reyes-Lega for an illuminating discussion. JCV thanks the Departamento de Fı́sica Teórica of the
Universidad de Zaragoza for warm hospitality.

CLBR and JMGB have been supported by grant FPA2009–09638 of Spain’s central government.
CLBR thanks Banco Santander for support. JMGB owes to ZiF for support, as well. JCV
acknowledges support from the Dirección General de Investigación e Innovación of Aragón’s
regional government, and from the Vicerrectorı́a de Investigación of the University of Costa Rica.

Last, but not least, we thank the referee for very helpful criticism, questions and suggestions,
leading to an improved presentation.

References
[1] W. Heisenberg, “Mehrkörperproblem und Resonanz in der Quantenmechanik”, Z. Physik 38 (1926), 411.

[2] C. Garrod and J. K. Percus, “Reduction of the 𝑁-particle variational problem”, J. Math. Phys. 5 (1964), 1756–
1776.

[3] P. W. Ayers, S. Golden and M. Levy, “Generalizations of the Hohenberg–Kohn theorem: I. Legendre transform
constructions of variational principles for density matrices and electron distribution functions”, J. Chem. Phys.
124 (2006), 054101.

[4] D. A. Mazziotti, “Two-electron reduced density matrix as the basic variable in many-electron quantum chemistry
and physics”, Chem. Rev. 112 (2012), 244.

[5] D. A. Mazziotti, “Structure of fermionic density matrices: complete 𝑁-representability conditions”, Phys. Rev.
Lett. 108 (2012), 263002.

[6] D. A. Mazziotti, “Significant conditions for the two-electron reduced density matrix from the constructive solution
of 𝑁-representability”, Phys. Rev. A 85 (2012), 062507.

[7] A. M. K. Müller, “Explicit approximate relation between reduced two- and one-particle density matrices”, Phys.
Lett. A 105 (1984), 446.

[8] P.-O. Löwdin and H. Shull, “Natural orbitals in the quantum theory of two-electron systems”, Phys. Rev. 101
(1956), 1730–1739.

18


[9] W. Kutzelnigg, “Die Lösung des quantenmechanischen Zwei-Elektronenproblems durch unmittelbare Bestim-
mung der natürlichen Einelektronenfunktionen”, Theor. Chem. Acta 1 (1963), 327.

[10] M. Moshinsky, “How good is the Hartree–Fock approximation?”, Amer. J. Phys. 36 (1968), 52–53.

[11] N. H. March, A. Cabo, F. Claro and G. G. N. Angilella, “Proposed definitions of the correlation energy from a
Hartree–Fock starting point: the two-electron Moshinsky model as an exactly solvable model”, Phys. Rev. A 77
(2008), 042504.

[12] P.-F. Loos, “Hooke’s law correlation in two-electron systems”, Phys. Rev. A 81 (2010), 032510.

[13] I. Nagy and J. Pipek, “Hartree–Fock partitioning of the two-matrix for an exactly solvable two-electron model
atom”, Phys. Rev. A 83 (2011), 034502.

[14] C. Amovilli and N. H. March, “Exact density matrix for a two-electron model atom and approximate proposals
for realistic two-electron systems”, Phys. Rev. A 67 (2003), 022509.

[15] I. Nagy and J. Pipek, “Approximations for the interparticle interaction energy in an exactly solvable two-electron
model atom”, Phys. Rev. A 81 (2010), 014501.

[16] C. Amovilli and N. H. March, “Quantum information: Jaynes and Shannon entropies in a two-electron entangled
artificial atom”, Phys. Rev. A 69 (2004), 054302.

[17] J. Pipek and I. Nagy, “Measures of spatial entanglement in a two-electron atom”, Phys. Rev. A 79 (2009), 052501.

[18] R. J. Yáñez, A. R. Plastino and J. S. Dehesa, “Quantum entanglement in a soluble two-electron model atom”,
Eur. Phys. J. D 56 (2010), 141–150.

[19] P. A. Bouvrie, A. P. Majtey, A. R. Plastino, P. Sá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”, Eur. Phys. J. D 66 (2012), 15.

[20] M. Srednicki, “Entropy and area”, Phys. Rev. Lett. 71 (1993), 666–671.

[21] E. R. Davidson, Reduced Density Matrices in Quantum Chemistry, Academic Press, New York, 1976.

[22] J. M. Gracia-Bondı́a, “Generalized Moyal quantization on homogeneous symplectic spaces”, in Deformation
Theory and Quantum Groups with Applications to Mathematical Physics, J. Stasheff and M. Gerstenhaber, eds.,
American Mathematical Society, Providence, RI, 1992; pp. 93–114.

[23] Ph. Blanchard, J. M. Gracia-Bondı́a and J. C. Várilly, “Density functional theory on phase space”, Int. J. Quant.
Chem. 112 (2012), 1134–1164.

[24] J. M. Gracia-Bondı́a and J. C. Várilly, “Exact phase space functional for two-body systems”, arXiv:1011.4742,
Bielefeld, 2010.

[25] J. P. Dahl, “Moshinsky atom and density functional theory. A phase space view”, Can. J. Chem. 87 (2009),
784–789.

[26] K. Ebrahimi-Fard and J. M. Gracia-Bondı́a, “Harmonium as a laboratory for mathematical chemistry”, J. Math.
Chem. 50 (2012), 440–454.

[27] C. L. Benavides-Riveros and J. M. Gracia-Bondı́a, “Physical Wigner functions”, Phys. Rev. A 87 (2013), 022118.

[28] C. L. Benavides-Riveros and J. C. Várilly, “Testing one-body density functionals on a solvable model”, Eur.
Phys. J. D 66 (2012), 274.

[29] K. Blum, Density Matrix Theory and Applications, Springer, Berlin, 2012.

19


[30] A. Kenfack and K. Życzkowski, “Negativity of the Wigner function as an indicator of non-classicality”, J. Opt.
B 6 (2004), 396.

[31] J. P. Dahl, H. Mack, A. Wolf and W. P. Schleich, “Entanglement versus negative domains of Wigner functions”,
Phys. Rev. A 74 (2006), 042323.

[32] J. M. Gracia-Bondı́a and J. C. Várilly, “Nonnegative mixed states in Weyl–Wigner–Moyal theory”, Phys. Lett. A
128 (1988), 20–24.

[33] A. P. Prudnikov, Yu. A. Brychkov and O. I. Marichev, Integrals and Series, Bell and Bain, Glasgow, 1983.

[34] G. E. Andrews, R. Askey and R. Roy, Special Functions, Encyclopedia of Mathematics and its Applications 71,
Cambridge University Press, Cambridge, 1999.

[35] A. Knutson and T. Tao, “Honeycombs and sums of Hermitian matrices”, Notices Amer. Math. Soc. 48 (2001),
175–186.

[36] A. A. Klyachko, “Quantum marginal problem and 𝑁-representability”, J. Phys. Conf. Ser. 36 (2006), 72–86.

[37] A. A. Klyachko, “The Pauli exclusion principle and beyond”, Ankara, 2009, arXiv:0904.2009.

[38] J. J. Duistermaat and G. J. Heckman, “On the variation in the cohomology of the symplectic form of the reduced
phase space”, Invent. Math. 69 (1982), 259–269.

[39] M. Christandl, B. Doran, S. Kousidis and M. Walter, “Eigenvalue distributions of the reduced density matrices”,
Commun. Math. Phys. 332 (2014), 1–52.

[40] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, 1965.

[41] E. Lieb, “Integral bounds for radar ambiguity functions and Wigner distributions”, J. Math. Phys. 31 (1990),
594–599.

[42] N. Helbig, I. V. Tokatly and A. Rubio, “Physical meaning of the natural orbitals: Analysis of exactly solvable
models”, Phys. Rev. A 81 (2010), 022504.

[43] G. T. Smith, H. L. Schmider and V. H. Smith, “Electron correlation and the eigenvalues of the one-matrix”, Phys.
Rev. A 65 (2002), 032508.

[44] J. Naudts and T. Verhulst, “Ensemble averaged entanglement of two-particle states in Fock space”, Phys. Rev. A
75 (2007), 062104.

[45] A. P. Balachandran, T. R. Govindarajan, A. R. de Queiroz and A. F. Reyes-Lega, “Entanglement, particle identity
and the GNS construction: a unifying approach”, Phys. Rev. Lett. 110 (2013), 080503.

[46] C. A. Coulson and A. H. Neilson, “Electron correlation in the ground state of helium”, Proc. Phys. Soc. 78
(1961), 831–837.

20