# Characteristic features of the strongly-correlated regime: Lessons from a
3-fermion one-dimensional harmonic trap

Victor Caliva  Universidad de Buenos Aires, Facultad de Ciencias Exactas y
Naturales, Departamento de Física. Buenos Aires, Argentina  CONICET -
Universidad de Buenos Aires, Instituto de Física de Buenos Aires (IFIBA).
Buenos Aires, Argentina  Johanna I. Fuks  Universidad de Buenos Aires,
Facultad de Ciencias Exactas y Naturales, Departamento de Física. Buenos
Aires, Argentina  CONICET - Universidad de Buenos Aires, Instituto de Física
de Buenos Aires (IFIBA). Buenos Aires, Argentina

###### Abstract

The transition into a strongly-correlated regime of 3 fermions trapped in a
one-dimensional harmonic potential is investigated. This interesting, but
little-studied system, allows us to identify characteristic features of the
regime, some of which are also present in strongly-correlated materials
relevant to the industry. Furthermore, our findings describe the behavior of
electrons in quantum dots, ions in Paul traps, and even fermionic atoms in
one-dimensional optical lattices. Near the ground state, all these platforms
can be described as fermions trapped in a harmonic potential. The correlation
regime can be controlled by varying the natural frequency of the trapping
potential, and to probe it, we propose to use twisted light. We identify 4
signatures of strong correlation in the one-dimensional 3-fermion trap, which
are likely to be present for any number $N$ of trapped fermions: i) the ground
state density is strongly localized with $N$ maximally separated peaks (Wigner
Crystal) ii) the symmetric and antisymmetric ground state wavefunctions become
degenerate (bosonization) iii) the von Neumann entropy grows, iv) the energy
spectrum is fully characterized by $N$ normal modes or less.

##  I Introduction

Forty years ago Feynman envisioned quantum emulators (also known as analogue
quantum computers) in his famous statement: “Nature isn’t classical, dammit,
and if you want to make a simulation of nature, you’d better make it quantum
mechanical, and by golly it’s a wonderful problem, because it doesn’t look so
easy” [1]. An analogue quantum computer is a controllable quantum system able
to simulate the behavior of another, more challenging (less controllable,
larger, more complex) quantum system [2, 3, 4]. Such a quantum emulator is not
universal, and unlike gate quantum computers, it can not be programmed and
there is no translation of an algorithm into a quantum circuit. Therefore, to
simulate a given problem, one needs to devise an experimentally feasible and
controllable quantum system able to emulate the behaviour of the original
problem. Several studies show the potential of noisy cold atoms and ion traps
as analogue quantum computers [8, 5, 6, 7]. Quantum advantage has been
demonstrated for the Gaussian boson sampling problem using a photon quantum
emulator[9].

Strong correlation is particularly challenging and has been widely studied,
both in-silico and experimentally [10, 11, 12, 13, 14, 15]. Stretched
molecular bonds involved in biochemical processes and transition-metal
catalysts for energy storage devices are some prominent examples of strongly-
correlated problems that are relevant to the industry. In the strongly-
correlated limit some materials exhibit novel phenomena not observed in weakly
correlated materials. Some examples are metal-insulator transitions [16],
magnetic order, unconventional superconductivity [17], superfluidity, quantum
phase transitions, and topological order [18, 19]. Designing and optimizing
materials with these properties is key to advancing technology and poses a
great challenge for the simulation community. Hartree Fock is a poor solution
for these types of problems and DFT, which is very successful in predicting
properties of weakly correlated materials, requires sophisticated exchange-
correlation functionals to achieve qualitatively correct results in this
regime [21, 20, 22]. Exchange-correlation functionals designed for strong
correlation, on the other hand, do not perform well for more weakly correlated
materials [23]. More expensive wavefunction methods such as CCSD also struggle
to simulate strong correlation [24] and Complete Active Space approaches
require large active spaces to converge [25]. The strongly-correlated regime
of quantum many-body systems is therefore an interesting application for
quantum computing.

In this work, we study one-dimensional harmonic traps of trapped fermions as a
testbed to learn about the characteristic features of the strongly-correlated
regime in fermionic many-body systems. The 2-fermion one-dimensional harmonic
trap was previously studied in Ref. [26]. We generalize the results to 3
fermions and infer the behaviour of $N$-trapped fermions. Analytical results
with cold fermionic atoms trapped in optical lattices have revealed that
phenomena previously investigated in the presence of many atoms may be studied
in the limit of a few particles as well [27, 28, 29, 30, 31]. Near the ground
state, electrons and holes in quantum wires and dots, ions in Paul traps, and
cold fermionic atoms in optical lattices are all realizations of the fermion
harmonic trap [32, 33, 34, 35]. Besides the importance of quantum dots in the
display industry, and of ion traps in the development of digital quantum
computers [36], these platforms can all serve as analogue quantum computers,
also known as quantum emulators. For atoms, the interaction is modeled as a
contact potential [37], whereas for electrons and ions the interaction is via
a Coulomb potential [38, 39] or soft-Coulomb [40, 41] for one-dimensional
traps. The confining potential is parabolic and its strength is characterized
by the natural frequency $\omega_{0}$. The correlation regime in harmonic
traps can be controlled by varying the particle-particle interaction strength
or the confinement strength of the parabolic potential $\omega_{0}$, and it
also depends on the number of trapped particles. By virtue of the high
tunability and the availability of analytic (weak and strong interaction
limits) [37, 42, 43] and very accurate one-dimensional numerical solutions
[44, 45], these systems provide a means to understand the physics and
interpret the experiments on more complex quantum many-body systems. The
problem of fermions trapped in a parabolic potential has been studied
theoretically for 3 electrons in a three-dimensional trap in the limits of
weak and strong confinement [42], in the context of electrons [27] [28], for 2
electrons in one and two dimensions [29], for ions [30] and for fermionic
atoms [31].

In this work, we derive expressions for the density and normal modes of 3
fermions in a one-dimensional harmonic trap in the strongly correlated limit
and use them to validate the numerical simulations. We then perturb the system
and analyze the transition into a strongly-correlated regime in the absorption
spectrum.

The paper is organized as follows: in section II we discuss the separability
of the harmonic trap hamiltonian into center of mass (CM) and internal degrees
of freedom and the antisymmetrization of the total wavefunction, in section
III we study the degeneration of the ground state and derive the normal modes
and an expression for the spatial wavefunction and the ground state density of
the 3-fermion harmonic trap in the strongly correlated limit $\omega_{0}\to
0$. In section IV we validate the finite $\omega_{0}$ numerical simulations
against the analytic results derived for the strongly correlated limit: we
confirmed the ground state density is localized as expected and the energy
spectrum contains the predicted frequencies in this limit. We then analyze the
quadrupole spectrum obtained from the numerical time-evolution of the density
after perturbation with a quadrupole field and discuss the selection rules. In
section V we discuss an experimental setup to probe correlation effects in the
quadrupole spectrum using twisted light and in section VI we conclude.

##  II Separability into center-of-mass (CM) and internal coordinates

The Hamiltonian describing a harmonic trap with $N$ particles is separable
into a center-of-mass (CM) term and a term dependent on the internal degrees
of freedom of the system. The CM system is equivalent to a particle with mass
$M=Nm$ and charge $Ne$,

${H_{0}}^{CM}=\frac{{P_{X}}^{2}}{2M}+\frac{1}{2}M{{\omega}_{0}}^{2}{X}^{2}$
(1)

with $X=\sum_{i}^{N}x_{i}$ y $P_{X}=\sum_{i}^{N}p_{i}$. For $N=3$ and defining
the Jacobi coordinates as

$q_{1}=\frac{x_{1}-x_{2}}{\sqrt{2}},q_{2}=\frac{x_{1}+x_{2}}{\sqrt{6}}-\sqrt{\frac{2}{3}}x_{3},$
(2)

the internal Hamiltonian reads,

${H_{0}}^{q}=\frac{{P_{1}}^{2}}{2m}+\frac{{P_{2}}^{2}}{2m}+\frac{1}{2}M{{\omega}_{0}}^{2}{q_{1}}^{2}+\frac{1}{2}M{{\omega}_{0}}^{2}{q_{2}}^{2}+V(q_{1},q_{2}).$
(3)

where $V(q_{1},q_{2})$ is the particle-particle interaction potential which is
the Coulomb potential [38] [39] $V(q_{1},q_{2})=-1/|q_{1}-q_{2}|$ in the case
of quantum dots and ion traps 111Rydberg atoms are also represented by this
Hamiltonian with $V(q_{1},q_{2})\propto 1/(q_{1}-q_{2})^{6}$ [37].. Atomic
units, $\hbar=e=m_{e}=a_{0}=1$, are used throughout the work.

Notice that due to the separability of the Hamiltonian the total wavefunction
$\Phi(x_{1},x_{2},x_{3})$ can be written as a product of a wavefunction
associated with the internal coordinates $\xi(q_{1},q_{2})$ and another
associated with the CM coordinates $\zeta(X)$:

$\Phi(x_{1},x_{2},x_{3})=\xi(q_{1},q_{2})\zeta(X).$ (4)

###  II.1 Fermionic wave function

In order for the total wave function Eq. (4) to represent fermions it must be
antisymmetric. The symmetry of the total wave function is obtained from the
product of a spin and a spatial function,

$\Psi(x,\sigma)=\Lambda(\sigma)\Phi(x).$ (5)

To characterize the symmetries of the spin wave function $\Lambda(\sigma)$ we
use the concept of Young diagrams, discussed in some depth in Appendix VII.
For the 3-fermion system, we obtain 8 spin states: 4 are pure symmetric, 2 are
symmetric mixed, and 2 are antisymmetric mixed. There are no pure
antisymmetric states for this problem [42].

The symmetries that define the spatial wave function $\Phi(x)$ are the Pauli
principle and the Parity [27],

$\Phi(-q_{1},q_{2},X)=-\Phi(q_{1},q_{2},X)\>\textup{( Pauli principle)}$ (6)
$\Phi(-q_{1},-q_{2},X)=\pm\Phi(q_{1},q_{2},X)\>\textup{(Parity)}$ (7)

The total wave functions $\Psi$ must be assembled [42]:

${\Psi}_{i}(1,2,3)=\hat{{\textit{A}}_{a}}{\Lambda}_{i}(1,2,3)\Phi(1,2,3)$ (8)

with the antisymmetrizer
$\hat{{\textit{A}}_{a}}=\frac{1}{\sqrt{N!}}\sum_{\textit{P}}(-1)^{p}\hat{\textit{P}}$,
and the permutation operator $\hat{\textit{P}}$ [42, 46]. ${\Lambda}_{i}$ is
the spin function, $i$ is equal to 1 for $S=3/2$ in pure symmetric case and
$S=1/2$ for mixed symmetric case, and equal to 2 and $S=1/2$ for mixed
antisymmetric case. The labels $(1,2,3)$ are spin or spatial coordinates
depending on the function in question. Spin functions are discussed in the
appendix.

##  III Analytic solution in the strong correlation limit

###  III.1 Degeneration of the ground state

In the strong correlation limit, when particle-particle repulsion dominates
over parabolic confinement, the two-particle problem has two degenerate ground
states with opposite symmetries in the spatial part of the wavefunction. A
solution differs from another solution only by a sign. There are odd functions
and even functions in the spatial part, with their corresponding change in the
symmetry of the spin part. The same procedure can be applied to larger
systems, $N>2$. For each pair of particles with opposite spins, it is possible
to change the symmetry of the spatial part of the wave function and build new
solutions in the limit of infinitely strong interaction. In this limit, there
are multiple degenerate ground states $D$, the value of $D$ is given by the
number of different spin configurations, counting the different possibilities
of dividing $N$ particles between the particles with spin-up $N_{\uparrow}$
and spin-down $N_{\downarrow}$ [47].

$D=\frac{N!}{N_{\uparrow}!N_{\downarrow}!}$ (9)

For a system with $N_{\uparrow}=1$ and $N_{\downarrow}=2$, or vice versa (the
account of $D$ contemplates the two configurations) we have $D=3$, i.e. three
wavefunctions with equal energy value.

###  III.2 Normal modes

In the work of Balzer [29] the normal modes of the one-dimensional 2-fermion
harmonic trap in the strong correlation limit are derived. The assumption is
that the wavefuntion is highly localized at $X_{0}$ and $-X_{0}$ in such limit
and the solution is found using a first order Taylor expansion around these
equilibrium positions. In this work we generalized Balzer’s approach for the
case of 3 fermions to find the charactersitic frequencies (normal modes) of
the system. The potential in the original coordinates reads,

$\begin{split}V(x_{1},x_{2},x_{3})=\frac{1}{2}{{\omega}_{0}}^{2}({x_{1}}^{2}+{x_{2}}^{2}+{x_{3}}^{2})\\\
+\frac{1}{|x_{1}-x_{2}|}+\frac{1}{|x_{2}-x_{3}|}+\frac{1}{|x_{1}-x_{3}|}\end{split}$
(10)

where $x_{1}$, $x_{2}$ y $x_{3}$ are the positions of each of the 3 electrons.
We set the gradient of the potential to zero to obtain the minimum equilibrium
positions $\overrightarrow{x_{0}}=(x_{1}=x_{0},x_{2}=0,x_{3}=-x_{0})$, with
$x_{0}=\frac{\sqrt[3]{5}}{2^{\frac{2}{3}}\sqrt[3]{{{\omega}_{0}}^{2}}}$, then
we calculate the Hessian matrix of the potential and evaluate it at
$\overrightarrow{x_{0}}$. The resulting matrix $A$ reads,

$A={\omega_{0}}^{2}\begin{bmatrix}14/5&-8/5&-1/5\\\ -8/5&21/5&-8/5\\\
-1/5&-8/5&14/5.\end{bmatrix}$ (11)

Diagonalizing matrix A, we get the normal modes

$\Omega_{0}=\omega_{0};\quad\Omega_{1}=\sqrt{3}\omega_{0};\quad\Omega_{2}=\sqrt{\frac{29}{5}}\omega_{0}=\sqrt{5.8}\omega_{0}.$
(12)

The values obtained for the normal modes in Eqs. (12) coincide with those
obtained by Daniel F. V. James in 1997 [30]. The first mode
$\Omega_{0}=\omega_{0}$ is associated with the center of mass. $\Omega_{1}$ is
associated with Coulombian interaction between particles and is called Thomas-
Fermi limit [49, 50] or Universal Breathing Mode (UBM), $\Omega_{1}$ is
present in all electronic systems and is independent of the number of
particles $N$ [49, 51, 48]. $\Omega_{2}$ on the other hand is specific to the
3-fermion one-dimensional harmonic trap [30].

With the normal modes obtained in Eqs. (12) we deduce the spectrum for 3
fermions,

$E_{\eta_{0}\eta_{1}\eta_{2}}=V(\overrightarrow{x_{0}})+{\omega}_{0}(\eta_{0}+\frac{1}{2})+\Omega_{1}(\eta_{1}+\frac{1}{2})+\Omega_{2}(\eta_{2}+\frac{1}{2}),$
(13)

The eigenenergies of the total system $E_{\eta_{0}\eta_{1}\eta_{2}}$ can be
expressed as the sum of the eigenenergies of the CM coordinates
${\epsilon_{\eta_{0}}}^{CM}$ plus the eigenenergies associated with the
internal coordinates $q_{1},q_{2}$, namely $\epsilon_{\eta_{1}\eta_{2}}$:

$E_{\eta_{0}\eta_{1}\eta_{2}}=\epsilon_{\eta_{1}\eta_{2}}+{\epsilon_{\eta_{0}}}^{CM}=\epsilon_{\eta_{1}\eta_{2}}+{\omega}_{0}(\eta_{0}+\frac{1}{2}),$
(14)

with $\eta_{0}\in\mathbb{N}$ being the quantum number associated with the CM,
and $\eta_{1},\eta_{2}\in\mathbb{N}$ being the quantum numbers associated to
the internal degrees of freedom. Applying Balzer’s method only to the internal
potential $V(q_{1},q_{2})$, we obtain a minimum position
$\overrightarrow{q_{0}}=(q_{1}=\pm q_{0},q_{2}=0)$, with $q_{0}=\sqrt{2}x_{0}$
and

$\epsilon_{\eta_{1}\eta_{2}}=V(\overrightarrow{q_{0}})+\Omega_{1}(\eta_{1}+\frac{1}{2})+\Omega_{2}(\eta_{2}+\frac{1}{2}),$
(15)

where $V(\overrightarrow{q_{0}})$ is the internal potential evaluated at the
minimum position $\overrightarrow{q_{0}}$, and $\eta_{1}$,$\eta_{2}$ represent
transitions between internal states.

Table 1: Analytical values of $\Omega_{1}$ and $\Omega_{2}$ for various values of $\omega_{0}$ (in atomic units) $\omega_{0}$ | 0.5 | 0.1 | 0.05 | 0.005
---|---|---|---|---
$\Omega_{1}$ | 0.8660 | 0.1732 | 0.0866 | 0.0087
$\Omega_{2}$ | 1.2042 | 0.2408 | 0.1204 | 0.0120

###  III.3 Spatial wavefunctions

According to Balzer [29], in the strongly correlated limit we can construct
the total spatial wave function of the system as a multiplication of the
elementary wave functions of the harmonic oscillator without Coulombian
interaction, decoupled, with each eigenfunction characterized by one of the
normal modes $\Omega_{i}$ in Eqs. (12),

$\begin{split}\phi_{\eta_{i}}(R_{i})=\prod_{i=0}^{Nd-1}(\frac{\Omega_{i}}{\pi\omega_{0}})^{\frac{1}{4}}\frac{1}{\sqrt{2^{\eta_{i}}\eta_{i}!}}e^{\frac{-\Omega_{i}{R_{i}}^{2}}{2}}H_{\eta_{i}}(\sqrt{\Omega_{i}}R_{i})\end{split}$
(16)

with $H_{\eta_{i}}(\sqrt{\Omega_{i}}R_{i})$ an Hermite polynomial, $N$ the
number of particles, $d$ the number of dimensions, and
$\overrightarrow{\eta}=(\eta_{0},...,\eta_{N-1})$ a set of labels that
represent the transition frequencies of the total system. In our case we have
3 fermions in one dimension, therefore 3 characteristic frequencies, whose
total spatial wave function can be written as the product of three elementary
wave functions of the decoupled harmonic oscillator and
$\overrightarrow{\eta}=(\eta_{0},\eta_{1},\eta_{2})$. Notice that the
wavefunction Eq. (16) is not antisymmetrized, but since in the strongly-
correlated limit symmetric and antisymmetric wavefunctions are degenerate in
energy the lack of antisymmetrization does not affect the estimation of the
energy (see Table 3). In Table 2 we show the analytic spatial wavefunctions,
Eq. (16), derived for the lowest 5 spatial wavefunctions of our 3-fermion one-
dimensional system, along with the corresponding transition frequencies.

Table 2: Total spatial wavefunctions in the strongly-correlated limit for 3-fermions trapped in a one-dimensional harmonic trap, derived using Balzer’s method [29], Eq. 16. The corresponding transition frequencies $(\epsilon_{\eta_{0},\eta_{1},\eta_{2}}-\epsilon_{0,0,0})$ are shown in the second column. Eigenfunctions $\ket{\phi_{\eta_{0}\eta_{1}\eta_{2}}}$ (only spatial part) | frequency
---|---
$\ket{\phi_{000}}=(\frac{\sqrt{3}\sqrt{5.8}}{\pi^{3}})^{\frac{1}{4}}{2}^{-\frac{3}{2}}e^{-\frac{(\omega_{0}{X}^{2}+\Omega_{1}(q_{1}-\frac{2x_{0}}{\sqrt{2}})^{2}+\Omega_{2}{q_{2}}^{2})}{2}}$ | 0
$\ket{\phi_{100}}=(\frac{\sqrt{3}\sqrt{5.8}}{\pi^{3}})^{\frac{1}{4}}{2}^{-\frac{1}{2}}\sqrt{\omega_{0}}Xe^{-\frac{(\omega_{0}{X}^{2}+\Omega_{1}(q_{1}-\frac{2x_{0}}{\sqrt{2}})^{2}+\Omega_{2}{q_{2}}^{2})}{2}}\propto X\ket{\phi_{000}}$ | $\omega_{0}$
$\ket{\phi_{010}}=(\frac{\sqrt{3}\sqrt{5.8}}{\pi^{3}})^{\frac{1}{4}}{2}^{-\frac{1}{2}}\times{3}^{\frac{1}{4}}\sqrt{\omega_{0}}(q_{1}-\frac{2x_{0}}{\sqrt{2}})e^{-\frac{(\omega_{0}{X}^{2}+\Omega_{1}(q_{1}-\frac{2x_{0}}{\sqrt{2}})^{2}+\Omega_{2}{q_{2}}^{2})}{2}}\propto(q_{1}-\frac{2x_{0}}{\sqrt{2}})\ket{\phi_{000}}$ | $\Omega_{1}$
$\ket{\phi_{200}}=(\frac{\sqrt{3}\sqrt{5.8}}{\pi^{3}})^{\frac{1}{4}}{2}^{-\frac{5}{2}}(4\omega_{0}X^{2}-2)e^{-\frac{(\omega_{0}{X}^{2}+\Omega_{1}(q_{1}-\frac{2x_{0}}{\sqrt{2}})^{2}+\Omega_{2}{q_{2}}^{2})}{2}}\propto(4\omega_{0}X^{2}-2)\ket{\phi_{000}}$ | $2\omega_{0}$
$\ket{\phi_{001}}=(\frac{\sqrt{3}\sqrt{5.8}}{\pi^{3}})^{\frac{1}{4}}{2}^{-\frac{1}{2}}\times{5.8}^{\frac{1}{4}}\sqrt{\omega_{0}}q_{2}e^{-\frac{(\omega_{0}{X}^{2}+\Omega_{1}(q_{1}-\frac{2x_{0}}{\sqrt{2}})^{2}+\Omega_{2}{q_{2}}^{2})}{2}}\propto q_{2}\ket{\phi_{000}}$ | $\Omega_{2}$

###  III.4 Ground state density

As $\omega_{0}\to 0$ and repulsion dominates over confinement, the fermions
localize maximally far from each other. This regime is known as Wigner crystal
in the condensed matter community [28, 40, 52, 53]. There has been recent
experimental evidence of this phenomenon [54].

To compute the ground state density we adapted Balzer’s Equation (24) in Ref.
[29] to 3 fermions:

$\rho_{\overrightarrow{\eta}}(x)=\int{dx_{1}dx_{2}dx_{3}|\Phi_{\overrightarrow{\eta}}(x_{1},x_{2},x_{3})}|^{2}\sum_{i=1}^{3}\delta(x_{i}-x).$
(17)

After replacing the wave functions with Eq. (16) we obtain the density as a
sum of 3 Gaussian functions:

$\begin{split}\rho(x)=A_{1}(\omega_{0},\Omega_{1},\Omega_{2})e^{-f_{1}(\omega_{0},\Omega_{1},\Omega_{2})(x-x_{0})^{2}}\\\
+A_{2}(\omega_{0},\Omega_{2})e^{-f_{2}(\omega_{0},\Omega_{2})x^{2}}\\\
+A_{1}(\omega_{0},\Omega_{1},\Omega_{2})e^{-f_{1}(\omega_{0},\Omega_{1},\Omega_{2})(x+x_{0})^{2}},\end{split}$
(18)

with
$f_{1}(\omega_{0},\Omega_{1},\Omega_{2})=\frac{\Omega_{1}\Omega_{2}}{\frac{\Omega_{2}}{2}+\frac{\Omega_{1}}{6}+\frac{\Omega_{1}\Omega_{2}}{3\omega_{0}}}$
and
$f_{2}(\omega_{0},\Omega_{2})=\frac{3\Omega_{2}}{2+\frac{\Omega_{2}}{\omega_{0}}}$.According
to the definition of a Gaussian function:
$\sigma_{i}^{2}=\frac{1}{2f_{i}(\omega_{0},\Omega_{1},\Omega_{2})}$, and
$A_{i}(\omega_{0},\Omega_{1},\Omega_{2})=\frac{1}{\sigma_{i}\sqrt{2\pi}}$. The
first and third terms represent each of the end bell curves, and the second
term represents the middle bell curve. It can be noted that the number of
peaks appearing in the density is equal to the number of fermions in the
system.

Figure 1: Comparison between the ground state density of the 3-fermion one-
dimensional harmonic trap with $\omega_{0}=0.005$ in the purely symmetric
state, obtained by numerical simulation in octopus (black dotted), and the
analytic expression found in the limit of $\omega_{0}\to 0$ Eq. (18) (red
line). Figure 2: Ground state density of the 3-fermion one-dimensional
harmonic trap with $\omega_{0}=0.005$ in the purely symmetric (upper panel)
and mixed (lower panel) cases of lowest energy. The calculations were
performed in Octopus using a simulation box of size $acell=100$ a.u and
$spacing=0.8$ a.u. Figure 3: The numerically simulated ground state density
for various correlation regimes corresponding to increasingly weaker trapping
potentials (increasingly smaller $\omega_{0}$), purely symmetrical case.
Figure 4:  The internal spectrum of the system for 3 different values of
$\omega_{0}$ to observe the transition towards the high correlation limit. At
$\omega_{0}=0.005$ it is observed that the spacing between eigenenergies is
equal to the frequencies of the normal modes.

##  IV Numerical simulations

The total Hamiltonian as well as the CM and internal Hamiltonians,
$H_{CM}(X)$, $H_{q}(q_{1},q_{2})$ and
$H_{Total}(X,q_{1},q_{2})=H_{CM}(X)+H_{q}(q_{1},q_{2})$ (Eqs. 1-3), were
solved numerically using the open-source software Octopus [55], with the
interaction between particles modeled as Soft-Coulomb,

$\begin{split}V_{S-Cou}(q_{1},q_{2})=\frac{1}{\sqrt{(\sqrt{2}q_{1})^{2}+{a}^{2}}}\\\
+\frac{1}{\sqrt{(\frac{-q_{1}+\sqrt{3}q_{2}}{\sqrt{2}})^{2}+{a}^{2}}}\\\
+\frac{1}{\sqrt{(\frac{q_{1}+\sqrt{3}q_{2}}{\sqrt{2}})^{2}+{a}^{2}}},\end{split}$
(19)

to avoid the singularities in the one-dimensional Coulombian potential 222The
Soft-Coulomb parameter $a$ models interaction properties. In the case of
quantum wires, the parameter $a$ mimics the short-range cutoff distance caused
by the transverse dimension of the system (for example the diameter of a
carbon nanotube or a quantum wire semiconductor) [40] [40, 41]. The parameter
$a$ is the Soft-Coulomb parameter, which we set to $a=1$.

In Table 3 we compare the internal eigenenergies derived in the strong
correlation limit for Coulomb and Soft-Coulomb interaction with the ones
obtained numerically with Soft-Coulomb and a finite but small $\omega_{0}$ in
Octopus. The analytic Coulomb versus Soft-Coulomb eigenergies differ by less
than $10^{-4}$ Hartree. These results confirm the overlap between the three
wavefunctions is very small (Wigner Crystal). Because of the small overlap the
lack of antisymmetrization of Eq. (16) does not affect the ground state
energy. The results also show the approximation introduced by the Soft-Coulomb
potential is very small. The discrepancy between the analytic and numerical
soft-Coulomb solutions, of the order of $10^{-3}$ Hartree, is due to the
finite value of $\omega_{0}$ used in the simulations, as can be seen from the
fact that the discrepancy is smaller for $\omega_{0}=0.001$ than it is for the
less correlated $\omega_{0}=0.005$ system.

Table 3: Internal eigenenergies derived for the strongly-correlated limit for Coulomb interaction (4th column) and Soft-Coulomb interaction (5th column) versus numerical solution for $\omega_{0}=0.005$ and $\omega_{0}=0.001$ and Soft-Coulomb interaction (6th column) Corr. regime | $\eta_{1}$ | $\eta_{2}$ | $\epsilon_{\eta_{1}\eta_{2}}$ |  $\epsilon_{\eta_{1}\eta_{2}}$ S-Cou | Numerical
---|---|---|---|---|---
$\omega_{0}=0.005$ | 0 | 0 | 0.112142 | 0.112120 | 0.112465
| 1 | 0 | 0.120802 | 0.120781 | 0.121274
| 0 | 1 | 0.124184 | 0.124162 | 0.124927
| 2 | 0 | 0.129463 | 0.129441 | 0.130133
| 1 | 1 | 0.132844 | 0.132822 | 0.133879
| 0 | 2 | 0.136225 | 0.136204 | 0.137636
$\omega_{0}=0.001$ | 0 | 0 | 0.034561 | 0.036881 | 0.036921
| 1 | 0 | 0.036293 | 0.038613 | 0.038671
| 0 | 1 | 0.036969 | 0.039290 | 0.039380
| 2 | 0 | 0.038025 | 0.040345 | 0.040427

Since the overlap between the three wavefunctions is very small (Wigner
Crystal), the lack of antisymmetrization of Eq. (16) does not affect the
ground state energy, which is well captured by this model, as can be confirmed
by the similar results obtained for analytic and numerical solutions in Table
3.

###  IV.1 Ground state density

In Fig. 1 we compare the density derived for the strong correlation limit Eq.
(18) with a numerically simulated density for a small $\omega_{0}$. We find
two distinct solutions for the numerical ground state density, corresponding
to the purely symmetric and mixed fermionic wavefunctions, depicted for the
strongly-correlated $\omega_{0}=0.005$ a.u. case in Fig. 2. The relationship
between the heights of the three peaks in the density follows the pattern
described in [40, 56]. In the strongly correlated regime, the pure symmetric
states shows a subtle difference in the peak heights, with all three peaks
being of similar height (Fig.1 in [56]). For mixed states, peaks with more
differentiated heights are observed in the density, whith the middle peak
being much shorter than the external peaks. We notice that the analytic
wavefunctions Eq. (16) are not antisymmetrized, therefore there is a unique
solution instead of two. Comparing the analytic and numerical densities we
observe that the Balzer wavefunctions best represent the pure antisymmetric
case. For the aim of comparison, we present the equilibrium positions
$\overrightarrow{x_{0}}$ and the width $\sigma$ of the bell curves calculated
from Eq. (18) in Tables 4 and 5 respectively. In Fig. 3 we show the purely
symmetrical ground state density for different correlation regimes
corresponding to increasing correlated systems $\omega_{0}=0.5;0.2;0.05;0.005$
a.u. We observe that as we approach the strongly correlation limit the number
of peaks in the density grows from one to two to three.

Table 4: $\omega_{0}=0.005$ a.u. (strong correlation regime): Absolute value $|x_{0}|$ (a.u.) of the equilibrium position $\overrightarrow{x_{0}}=(-x_{0},0,x_{0})$, from the generalization of Balzer’s solution to 3 fermions (analytic, not anisymetrized) versus numerical solution for symmetric and mixed cases. $\omega_{0}$ |  $\overrightarrow{x_{0}}$ analytic |  $\overrightarrow{x_{0}}$ symmetric |  $\overrightarrow{x_{0}}$ mixed
---|---|---|---
0.05 | 7.937 | 8.739 | 8.739
0.005 | 36.840 | 38.907 | 38.926
Table 5: Width Comparison ($\sigma$) of the density bell curves obtained analytically and numerically for a value of $\omega_{0}=0.005$ (high correlation). $\omega_{0}$ | Bell curves | $\sigma_{analitico}$ | $\sigma_{sim}$ | $\sigma_{mix}$
---|---|---|---|---
0.05 | Middle | 2.470 | 2.099 | 2.099
0.05 | End | 2.629 | 2.707 | 2.707
0.005 | Middle | 7.8113 | 7.4129 | 6.9870
0.005 | End | 8.2705 | 8.2676 | 8.2956

###  IV.2 von Neumann entropy

In the case of the 2-fermion trap it was observed that, as we approach the
strongly-correlated regime, the von Neumann entropy grows [26]. We found that
this is also the case for the 3-fermion trap and is likely to be the case for
any $N$. The von Neumann entropy is a measure of the entanglement present in
the system [57]. In this problem, strong-correlation implies large
entanglement, two deeply related concepts that are sometimes elusive to
differentiate. In Table 6 we show the von Neuman entropy $s^{\omega_{0}}$ for
various values of $\omega_{0}$,

$s^{\omega_{0}}=\frac{1}{N}\sum_{i=1}^{N}\eta^{\omega_{0}}_{i}ln(\eta^{\omega_{0}}_{i})$
(20)

where $\eta^{\omega_{0}}_{i}$ is the Natural Orbital Occupation Number of
Natural Orbital $i$ , and $N$ is the number of particles in the system. To
obtain $\eta^{\omega_{0}}_{i}$ we diagonalized the 1-Body Reduced Density
Matrix $\Gamma^{(1)}(x;x^{\prime})$ for various $\omega_{0}$,

$\begin{split}\Gamma_{\omega_{0}}^{(1)}(x;x^{\prime})&=\int
dx_{2}dx_{3}\Psi_{\omega_{0}}^{*}(x_{1}^{\prime},x_{2}^{\prime},x_{3}^{\prime})\Psi_{\omega_{0}}(x_{1},x_{2},x_{3})\\\
&=\sum_{j=1}^{\infty}\eta_{j}^{\omega_{0}}\varphi_{j}^{\omega_{0}}(x)\varphi_{j}^{\omega_{0},*}(x^{\prime})\end{split}$
(21)

Table 6: Von Neumann entropy for various values of $\omega_{0}$ in a.u. $\omega_{0}$ | 0.5 | 0.1 | 0.05 | 0.005
---|---|---|---|---
$s^{\omega_{0}}$ | 0.169 | 0.329 | 0.367 | 0.386

Notice that in the Hartree-Fock and DFT methods the von Neumann entropy is
always zero, because all $\eta_{i}$ are either 0 or 1 [58]. Thus for these
efficient ab-initio methods simulating the strongly-correlated regime is a
challenge. In Ref. [26] it was shown that DFT using LDA approximation fails to
capture the localization of the density in the strongly correlated limit in
the 2-fermion harmonic one-dimensional harmonic trap case (see inset Fig.3 in
Ref. [26]). Even if other exchange-correlation approximations may capture the
strong-correlation limit [23] the simulation of the transition from a weakly
to a strongly-correlated regime (for example dissociation curves of molecules)
remains a challenge.

###  IV.3 Energy spectrum and selection rules

The internal energy spectrum (ground state plus 18 unoccupied states) was
calculated numerically in Octopus and is shown in Fig. 4. We observe that as
we make $\omega_{0}$ smaller and approach the strongly-correlated regime, the
spacing between internal eigenenergies converges to the normal modes. for
computational efficiency we calculated the spectrum of the internal system
instead of the spectrum of the total system. Since the Hamiltonian is
separable the total response of the system is the sum of the CM and internal
responses. The CM spectrum is that of a harmonic oscillator of natural
frequency $\omega_{0}$.

To measure the internal response we performed a time evolution of the internal
density after perturbation with an instantaneous quadrupole external field
(quadrupole kick). For a sufficiently small applied field, the response can be
expanded in a series of powers in the intensity $E_{0}$ of the field. The
linear ($\propto E_{0}$) response of an observable $\hat{O}$ is proportional
to the matrix element
$\braket{{\hat{O}}^{(1)}}\propto\sum_{k=1}^{N}\bra{\phi_{0}}\hat{O}\ket{\phi_{k}}\bra{\phi_{k}}\hat{H_{I}}\ket{\phi_{0}}$
[26][59]. By performing the Fourier transform of the evolved density, we
obtain the frequency spectrum, where the peaks correspond to transition
frequencies between the ground and excited states that are not forbidden by
symmetry. The allowed transitions depend on the symmetry of the system and
also of the excitation field (selection rules). Notice that to excite internal
degrees of freedom and study correlation effects we need a quadrupole or
higher field because a dipole field (plane waves) can only excite CM
transitions (Kohn’s theorem [32, 38]). To study correlation effects in the
quadrupole spectrum, a light field with spatial structure, for example,
optical vortices (also known as twisted light), is proposed, the details of
the setup are discussed in the next section.

In Fig. 5 we show the internal quadrupole spectrum obtained from Fourier
transformation of the time-resolved internal density for $\omega_{0}=0.02$ and
$\omega_{0}=0.005$ a.u.. The computational details are included in the
Appendix. We plot the absolute value of the quadrupole spectrum against a
dimensionless frequency (expressed in units of $\omega_{0}$). The peaks are
close but not on top for the two values of $\omega_{0}$. We notice that for
$\omega_{0}=0.02$ a.u. there is an additional peak, this is because for this
less correlated system, the spectrum is not yet the simple equidistant
spectrum of the strongly-correlated limit (see Fig. 4). We observe that the
position of the peaks (corresponding to the allowed transition frequencies)
approximately obey:
$n.(\Omega_{1}+\Omega_{2})/(3!\omega_{0})=n.(\sqrt{3}+\sqrt{5.8})/6=n.\Omega_{eff}$,
with $n\in\mathbb{Z}$ and $3!=N!$. This would mean that there is only one
independent internal quantum number $n_{eff}$ instead of two. In Table 7 we
show the numerical values of the peak positions obtained by Octopus and
compare them with the analytical values. We observe that for $n_{eff}=3;4$ the
behavior of $n_{eff}$ is reassuring, namely the more correlated
$\omega_{0}=0.005$ system shows very good agreement with the strong-correlated
limit analytic solution, better than the less correlated $\omega_{0}=0.02$
system. But the behavior of $n_{eff}=1;2$ doesn’t fit so well into the
hypothesis of an equidistant spectrum. The explanation of these observations
is a work in progress and will be the topic of a future publication.

Table 7: Peak positions for $\omega_{0}=0.02$ and $\omega_{0}=0.005$ a.u. (in units of $\omega_{0}$). The first column is the internal quantum number $n$, such that the peaks are at $n\Omega_{eff}$ (see discussion in the text). n | $\omega_{0}$ | Numerical | Analytic ($\omega_{0}\to 0)$)
---|---|---|---
1 | 0.02 | 0.6284 | 0.6901
| 0.005 | 0.7540 | 0.6901
2 | 0.02 | 1.3198 | 1.3801
| 0.005 | 1.2567 | 1.3801
3 | 0.02 | 1.8845 | 2.0702
| 0.005 | 2.0101 | 2.0702
4 | 0.02 | 2.6389 | 2.7602
| 0.005 | 2.7647 | 2.7602

Figure 5: Quadrupole spectra (in arbitrary units). The frequency is divided by
$\omega_{0}$. Upper panel $\omega_{0}=0.02$, middle panel $\omega_{0}=0.005$
a.u. for various values of the spin angular momentum $m$: $Q_{m=0}$ (red),
$Q_{m=2}$ (blue) and $Q_{m=-2}$ (purple). Lower panel: comparison of
$Q(l=m=2)$ for $\omega_{0}=0.005$ (blue) versus $\omega_{0}=0.02$ a.u. (cyan).

##  V Measurement of the quantum observables

The transition into a strongly-correlated regime can be observed in the
absorption spectrum, i.e. the response of the system to a perturbation with a
light field. A homogenous light field, however, can only excite the CM degrees
of freedom in harmonic traps, because it couples to the dipole moment of the
system, which is a CM variable (see Kohn’s theorem [32]). We therefore need a
quadrupole field or higher to study the correlation between particles, such
that internal degrees of freedom do get excited. We propose to use the spatial
structure of a twisted light field to excite dipole-forbidden transitions and
study the correlation effects [26, 60, 63, 64, 61, 62]. Such a setup could be
used to probe the correlation regime but also to selectively excite the system
into a desired dipole-forbidden state, i.e. preparing the system in a
particular initial state. The latter may have applications in digital ion trap
quantum computers, whose working regime is the strongly correlated limit we
have investigated in this paper [65, 8]

Twisted light, also known as optical vortices, designates a family of highly
non-homogeneous optical beams with single or multiple phase singularities
carrying orbital angular momentum (OAM), among other interesting features [61,
62, 69, 70, 68, 66, 67]. We propose to use the spatial inhomogeneity
associated with the topological charge $l$ of the optical vortices to study
the internal energy spectrum of the trap. This idea was first introduced in
Ref. [26] to study 2 harmonically trapped fermions. We model the interaction
between the system and an optical vortex of topological charge $l=1$ and $l=2$
as,

${H_{I}}^{l=1}=\frac{eE_{0}q_{r}}{2\sqrt{2}}(X^{2}+{q_{1}}^{2}+{q_{2}}^{2})sen(\omega
t)$ (22)
${H_{I}}^{l=2}=\frac{eE_{0}{q_{r}}^{2}}{12\sqrt{2}}X(X^{2}+{q_{1}}^{2}+{q_{2}}^{2})sen(\omega
t)$ (23)

where $E_{0}$ and $q_{r}$ are the amplitude and the waist of the beam
respectively. Details on the derivation of Eqs. (22)(23) and the experimental
setup for the twisted light probe can be found in Ref. [26]. The
specifications for the twisted light probe depend on the size of the system
(which can be estimated from the density) rather than on the number of
fermions in the trap.

In Fig. 5 we present the internal quadrupole spectrum of a 3-fermion harmonic
trap after perturbation with a twisted light field of $l=1$, whose interaction
with the system is described by Eq. (22). The CM quadrupole spectrum has peaks
at even multiples of the natural frequency, i.e. $2\eta_{0}\omega_{0}$ where
$\eta_{0}$ is the CM quantum number (not shown).

##  VI Conclusions

We were able to elucidate 4 features that characterize the strongly-correlated
regime in one-dimensional fermion harmonic traps. Our findings could be used
to probe the interactions between the electrons in quantum dots or between the
ions in Paul traps. Any of these platforms could be used as an analogue
quantum computer to simulate the transition into a strongly correlated regime
in more complex quantum many-body systems, something hard to explore with
classical processors. Wigner crystallization [52, 53, 54], bosonization [72,
71, 73, 74, 75, 76, 77, 78], and metal-insulator transitions [16] are some
examples of phenomena that can emerge as we approach the strongly-correlated
regime.

In addition, identification of the characteristic signatures of strong
correlation allows us to identify the quantum observables we need to track and
is helpful for the interpretation of experimental results on strongly-
correlated materials. Some of the features we identified as signatures of the
strongly-correlated regime are common to strongly-correlated materials
relevant to the industry, such as transition metals (used as catalysts in
energy storage devices) and stretched molecules involved in chemical
reactions. In these systems, Hartree Fock is a bad approximation (or what is
the same, the von Neumann entropy is non-zero) and the energy levels become
degenerate, which makes convergence for ab initio methods such as DFT,
challenging 333Degenerate states make the convergence of the self-consistent
field calculation in HF and DFT hard, requiring tricks such as level shift or
smearing..

In this work we generalize the findings of [26] from 2 to 3-fermion one-
dimensional traps. We expect the same signatures of strong correlation to be
present in the more general case of $N$-fermion traps. Generalization to
higher dimensional traps is the scope of future work.

We have derived an analytical expression for the spatial part of the 3-fermion
wavefunction in the highly correlated limit and we show that, in this limit,
the energy spectrum of the 3-fermion one-dimensional harmonic trap is fully
determined by 3 normal modes: two coincide with the modes of the 2-fermion
trap [26, 29, 43], namely $\omega_{0}$ associated with the center of mass and
the universal breathing mode $\Omega_{1}=\sqrt{3}\omega_{0}$, and the third
normal mode $\Omega_{2}=\sqrt{5.8}\omega_{0}$ is specific to the 3-fermion
one-dimensional trap [30]. The many-body wavefunction can be approximated by a
product of Gaussian functions[45] (independent oscillators with characteristic
frequencies equal to the normal modes) [79, 29]. We used these analytic
results to validate the numerical simulations.

The ground state wavefunction was simulated numerically using soft-Coulomb
interaction and the open source Octopus code [55] for different values of the
trap frequency $\omega_{0}$. As a check, we ensured the simulated internal
wavefunction respects parity (see Fig. 6 in the Appendix). We also checked
that the use of Soft-Coulomb interaction, instead of full Coulomb, does not
significantly affect the density. We found good agreement between the
analytical density in the strongly-correlated limit and the simulated
symmetric ground state density for small but finite $\omega_{0}$ (see Fig.1).

One of the 4 characteristic features of the strongly-correlated limit is that
the fermions localize maximally far from each other (Wigner Crystal, see Fig.
3). This localization renders the fermions distinguishable, because it becomes
possible to label them. We show that for 3 fermions the total wavefunction can
be written as a product of 3 Gaussian functions with little overlap among them
in the Wigner Crystal regime. As a consequence, the antisymmetrization of the
wavefunction does not affect the energy, and thus symmetric (bosons) and
antisymmetric (fermions) solutions become degenerate, which is the second
feature we identify as characteristic of the strongly-correlated limit. There
has been experimental evidence of this symmetry [72, 80], also known as
bosonization [72, 71, 73, 74, 75, 76, 77, 78]. The third characteristic
feature we have identified is the growth of the van Neumann entropy [57] (see
table 6), which is a measure of the entanglement present in the system. The
forth characteristic feature is the simplification of the energy spectrum. As
we approach $\omega_{0}$ the spectrum is fully characterized by $N$
characteristic frequencies corresponding to the $N$ normal modes. For
intermediate $\omega_{0}$ the energy spectrum has many more distinct energy
levels (see Fig. 4).

To study the fingerprints of particle-particle correlation in the spectrum we
perturbed the system with a quadrupole field. A twisted light probe could be
used for this purpose as proposed in Ref. [26, 81]. We observe that as we
approach $\omega_{0}\to 0$, the peaks in the simulated spectrum approach
combinations of the normal modes $\omega_{0}$, $\Omega_{1}$, $\Omega_{2}$, as
expected. The positions of the peaks in the simulated quadrupole spectrum (see
Fig. 5) are observed to be approximately proportional to the sum of
$\Omega_{1}$ and $\Omega_{2}$ divided by the factorial of the number of
particles $N!=3!$, which suggests that we can define this sum as an effective
internal natural frequency $\Omega_{eff}=(\Omega_{1}+\Omega_{2})/3!$. The
explanation of these selection rules is the scope of future work.

##  VII Appendix: Young tableaux of the wavefunctions of the system of 3
electrons.

A Young diagram is a mathematical object made up of rows of boxes, the
arrangement of boxes is such that each row has a quantity less than or equal
to the top row. This object is related to Group Theory. A Young Tableau is
defined through the placement of symbols, particularly numbers, into the boxes
of the corresponding Young diagram[82].

centertableaux 1 & 2 3 4 5
6 7 8
9 10

Different diagrams have several standard tableaux, which are the result of the
different ways of arranging the numbers in Young’s tableaux [46, 82] but
arranged in such a way that an ascending sequence can be read both by rows
(from left to right) and columns (from top to bottom). For example, for the
previous example, there are 10 standard tableaux.

For 2 electrons, two diagrams can be associated with the singlet and triplet
configurations, respectively. When there are 3 electrons involved, there are
two diagrams in which 2 electrons occupy one spin channel while the remaining
electron occupies the other channel. Additionally, there is a single diagram
depicting all electrons with the same spin orientation. Our research focuses
on a system consisting of 3 electrons, where the electron is affiliated with
the second-order special unitary group SU(2). An electron is represented by a
box $\square$. Initially, by multiplying two boxes $\square\otimes\square$, we
obtain the diagrams for the case of two electrons. These diagrams result in
the tableaux illustrating the symmetrical case (triplet), depicted as three
rows with two boxes each, and for the antisymmetrical case (singlet), the
tableau associated is a column with two boxes, we obtain

$2\otimes 2=3\oplus 1$

. For 3 electrons:

$2\otimes 2\otimes 2=4\oplus 2\oplus 2$

where $\otimes$ is the direct product and $\oplus$ the direct sum. We get 4
symmetrical states, 2 symmetric mixed states, and 2 antisymmetric mixed states
[42].

For the 4 symmetrical states, we have a Young diagram represented by a single
row (for pure antisymmetric states there would be a single column):
centertableaux 1 & 2 3

and using the notation $\ket{S,M_{S}}$, $S$ spin and $M_{S}$ spin projection,
the states are:

$\ket{\frac{3}{2},\frac{3}{2}}=\uparrow\uparrow\uparrow$
$\ket{\frac{3}{2},-\frac{3}{2}}=\downarrow\downarrow\downarrow$
$\ket{\frac{3}{2},\frac{1}{2}}=\frac{1}{\sqrt{3}}(\uparrow\downarrow\uparrow+\downarrow\uparrow\uparrow+\uparrow\uparrow\downarrow)$
$\ket{\frac{3}{2},-\frac{1}{2}}=\frac{1}{\sqrt{3}}(\uparrow\downarrow\downarrow+\downarrow\uparrow\downarrow+\downarrow\downarrow\uparrow)$

To define symmetric and antisymmetric mixed states, since the spin space is
degenerate, there is more than one orthogonal spin eigenstate $\Lambda_{i}$
for $S$ y $M_{S}$ set. We use the notation $\ket{S,M_{S},i}$ for those states.

For the 2 symmetric mixed states, we have the tableau: centertableaux 1 & 2
3
with states:

$\ket{\frac{1}{2},\frac{1}{2},1}=\frac{1}{\sqrt{6}}(2\uparrow\uparrow\downarrow-\uparrow\downarrow\uparrow-\downarrow\uparrow\uparrow)$
$\ket{\frac{1}{2},-\frac{1}{2},1}=\frac{1}{\sqrt{6}}(2\downarrow\downarrow\uparrow-\uparrow\downarrow\downarrow-\downarrow\uparrow\downarrow)$

For the 2 antisymmetric mixed states, we have the tableau: centertableaux 1 &
3
2

with states:

$\ket{\frac{1}{2},\frac{1}{2},2}=\frac{1}{\sqrt{2}}(\uparrow\downarrow\uparrow-\downarrow\uparrow\uparrow)$
$\ket{\frac{1}{2},-\frac{1}{2},2}=\frac{1}{\sqrt{2}}(\uparrow\downarrow\downarrow-\downarrow\uparrow\downarrow)$

As in Taut’s work [42], $i=1$ labels symmetric mixed states and $i=2$ labels
the antisymmetric states.

##  VIII Computational Details

The ground state calculations were computed in Octopus version 10 [55] using
the following simulation boxes: i) $\omega_{0}=0.5$ a.u.: $20$ a.u. simulation
box with $0.5$ a.u. spacing. ii) $\omega_{0}=0.1$ a.u.: $20$ a.u. simulation
box with $0.5$ a.u. spacing. iii) $\omega_{0}=0.05$ a.u.: $100$ a.u.
simulation box with $0.8$ a.u. spacing. iv) $\omega_{0}=0.005$ a.u. : $100$
a.u. simulation box with $0.8$ a.u. spacing.

The exact calculations were done using the modelmb functionality in Octopus
[83] and the code was modified to output the quadrupole moment in one
dimension.

The interaction with the TL field is modeled as in Ref. [26],

$H_{I}^{l}=E_{0}f(t)\sum_{i}^{N}x_{i}^{(l+1)}$ (24)

assuming dimension and position of the harmonic trap w.r.t. to the twisted
light beam fulfills the requirements stated in the same Ref. When $l=0$ the
field couples to the dipole and only the center of mass degrees of freedom get
excited (Kohn’s theorem [32]). For the time dependence, we choose a step
function $f(t)=\text{rect}(t-\tau/2)$, where $\tau=0.01a.u.$ is the duration
of the pulse. The first-order quadrupole response to the TL probe can be
computed as

$Q^{(1)}(\omega)=FT[Q^{(1)}(t)]/FT[{\cal E}(t)]$ (25)

where ${\cal E}(t)=E_{0}f(t)$ is the external perturbation due to the
interaction between the harmonic trap and the TL field. The time-dependent
quadrupole can be computed from the evolution of the density $\rho(x,t)$ as

$Q(t)=-e\langle\sum_{i}^{N=2}\hat{x}_{i}^{2}\rangle=-e\int\rho(x,t){\bf
x}^{2}dx.$ (26)

$Q^{(1)}(t)$ in Eq. (25) is the variation of the quadrupole moment due to the
action of a TL probe of odd topological charge $l$, i.e. $Q^{(1)}(t)=\delta
Q(t)=Q(t)-Q^{0}$, with $Q^{0}=-e\int\rho^{0}(x,t){\bf x}^{2}dx$ being the
unperturbed quadrupole moment. The contribution of $FT[\text{rect}(t-\tau/2)]$
is only significant around $\omega=0$ and will be ignored here. The simulation
was set with an value of $E_{0}=0.001$, $dt=0.01$ the timestep, $T=6000$ the
total propagation time, all these values in a.u. The propagator used was
Approximated Enforced Time-Reversal Symmetry (AETRS).

###  VIII.1 Parity check on simulated wave function

Figure 6: First 7 internal spatial eigenstates and their respective parities
for $\omega_{0}=0.005$ a.u. (strong-correlation regime). In the strongly-
correlated limit each state has degeneracy 2.

To ensure the simulated wavefuction has the proper symmetry we performed a
parity check. Since the Hamiltonian is separable and the center of mass
wavefuction $\zeta(X)$ commutes with parity, we checked parity for the more
interesting internal wavefunction $\xi(q_{1},q_{2})$, which needs to preserve
parity as well, $\xi(-q_{1},-q_{2})=\pm\xi(q_{1},q_{2})$ [27]. We perform the
check for the most correlated case we computed, namely $\omega_{0}=0.005$
a.u., which is also the hardest to converge. In Fig. 6 we plot the wave
functions $\xi(q_{1},q_{2})$ and $\xi(-q_{1},-q_{2})$ of the internal system
for a harmonic trap of $\omega_{0}=0.005$ a .u. For each set of
$\\{\eta_{1},\eta_{2}\\}$, there are two wavefunctions: one even and one odd.

## References

  * [1] R. P. Feynman, ”Simulating physics with computers”, InterNatl. J. (Wash.) of Theoretical Physics 21 (1982) 467.
  * [2] Kendon, Vivien M. and Nemoto, Kae and Munro, William J., ”Quantum analogue computing”, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 368.1924, 2010:3609-3620,
  * [3] Thomas Ayral, Pauline Besserve, Denis Lacroix, Edgar Andres Ruiz Guzman, ”Quantum computing with and for many-body physics”, arxiv.org/abs/2303.04850 (2023)
  * [4] Antoine Michel, Loïc Henriet, Christophe Domain, Antoine Browaeys, Thomas Ayral, ”Hubbard physics with Rydberg atoms: using a quantum spin simulator to simulate strong fermionic correlations”, https://arxiv.org/abs/2312.08065 (2023)
  * [5] Fraxanet, Joana, Tymoteusz Salamon, and Maciej Lewenstein. ”The coming decades of quantum simulation.” Sketches of Physics: The Celebration Collection. Cham: Springer International Publishing, 2023. 85-125.
  * [6] Tangpanitanon, Jirawat, et al. ”Signatures of a sampling quantum advantage in driven quantum many-body systems.” Quantum Sci. Technol. 8.arXiv: 2002.11946 (2023).
  * [7] Rahul Trivedi, Adrian Franco Rubio, J. Ignacio Cirac, ”Quantum advantage and stability to errors in analogue quantum simulators”, https://arxiv.org/abs/2212.04924 (2022)
  * [8] Hempel, Cornelius, et al. ”Quantum chemistry calculations on a trapped-ion quantum simulator.” Physical Review X 8.3 (2018): 031022.
  * [9] Zhong, Han-Sen, et al. ”Quantum computational advantage using photons.” Science 370.6523 (2020): 1460-1463.
  * [10] Langer, Judith, et al. ”Present and future of surface-enhanced Raman scattering.” ACS nano 14.1 (2019): 28-117.
  * [11] Cioslowski, Jerzy, et al. ”Contactium: A strongly correlated model system.” The Journal of Chemical Physics 158.18 (2023).
  * [12] Grzeszczyk, M., et al. ”Strongly Correlated Exciton‐Magnetization System for Optical Spin Pumping in CrBr3 and CrI3.” Advanced Materials (2023): 2209513.
  * [13] Chatzieleftheriou, Maria, et al. ”Mott Quantum Critical Points at finite doping.” Physical Review Letters 130.6 (2023): 066401.
  * [14] Pavosevic, Fabijan, Ivano Tavernelli, and Angel Rubio. ”Spin-flip unitary coupled cluster method: Toward accurate description of strong electron correlation on quantum computers.” The Journal of Physical Chemistry Letters 14.35 (2023): 7876-7882.
  * [15] F.Serwane, G. Zürn, T. Lompe, T. B. Ottenstein, A. N. Wenz, S. Jochim, ”Deterministic preparation of a tunable few-fermion system”. Science 332.6027 (2011): 336-338.
  * [16] Feng, Yejun and Wang, Yishu and Silevitch, D. M. and Cooper, S. E. and Mandrus, D. and Lee, Patrick A. and Rosenbaum, T. F., ” A continuous metal-insulator transition driven by spin correlations”, Nature Communications 2779.12 (2021)
  * [17] Adler, Ran, et al. ”Correlated materials design: prospects and challenges.” Reports on Progress in Physics 82.1 (2018): 012504.
  * [18] Georgescu, Iulia M., Sahel Ashhab, and Franco Nori. ”Quantum simulation.” Reviews of Modern Physics 86.1 (2014): 153.
  * [19] Tarruell, Leticia, and Laurent Sanchez-Palencia. ”Quantum simulation of the Hubbard model with ultracold fermions in optical lattices.” Comptes Rendus Physique 19.6 (2018): 365-393.
  * [20] Gräfenstein, Jürgen and Kraka, Elfi and Cremer, Dieter, ”Effect of the Self-Interaction Error for Three-Electron Bonds: On the Development of New Exchange-Correlation Functionals”, Physical Chemistry Chemical Physics 6 (2024)
  * [21] Cohen, Aron J. and Mori-Sánchez, Paula and Yang, Weitao, ”Challenges for Density Functional Theory”, Chemical Reviews, 112.1 (2012):289-320
  * [22] Prokopiou, Georgia and Hartstein, Michal and Govind, Niranjan and Kronik, Leeor, ”Optimal Tuning Perspective of Range-Separated Double Hybrid Functionals”,Journal of Chemical Theory and Computation 18.4 (2022):2331-2340
  * [23] Malet, Francesc and Gori-Giorgi, Paola, ”Strong Correlation in Kohn-Sham Density Functional Theory”, Phys. Rev. Lett., 109.24, (2012):246402
  * [24] Bulik, Ireneusz W. and Henderson, Thomas M. and Scuseria, Gustavo E., ”Can Single-Reference Coupled Cluster Theory Describe Static Correlation?”, Journal of Chemical Theory and Computation 11.7 (2015):3171-3179
  * [25] Marti, Konrad H. and Ondík, Irina Malkin and Moritz, Gerrit and Reiher, Markus,”Density matrix renormalization group calculations on relative energies of transition metal complexes and clustes”, The Journal of Chemical Physics 128.1 (2008):014104
  * [26] Johanna I. Fuks, Guillermo F. Quinteiro Rosen, Heiko Appel, Pablo I. Tamborenea, ”Probing many-body effects in harmonic traps with twisted light”, Phys. Rev. B, 107.8 (2023):L081111
  * [27] N.J.S.Loft, A.S.Dehkharghani, N.P.Mehta, A.G.Volosniev, and N.T.Zinner, ”A variational approach to repulsively interacting three-fermion systems in one-dimensional harmonic trap”, The European Physical Journal D 69.3 (2015): 1-17.
  * [28] Jerzy Cioslowski, Katarzyna Pernal, ”Wigner molecules: The strong-correlation limit of the three-electron harmonium”, The Journal of chemical physics 125.6 (2006): 064106.
  * [29] K.Balzer, ”Energy spectrum of strongly correlated particles in quantum dots”, Journal of Physics: Conference Series Vol. 35, No. 1. IOP Publishing. 2006: p. 019.
  * [30] James, Daniel FV. ”Quantum dynamics of cold trapped ions with application to quantum computation”. No. quant-ph/9702053. 1997.
  * [31] D’Amico, Pino, and Massimo Rontani. ”Three interacting atoms in a one-dimensional trap: a benchmark system for computational approaches.” Journal of Physics B: Atomic, Molecular and Optical Physics 47.6 (2014): 065303.
  * [32] S.K.Yip, ”Magneto-optical absorption by electrons in the presence of parabolic confinement potentials”, Physical Review B 43.2 (1991): 1707.
  * [33] Jacak, Lucjan, Pawel Hawrylak, and Arkadiusz Wojs. _Quantum dots_. Springer Science & Business Media, 2013.
  * [34] Sarkisyan, Hayk A., et al. ”Realization of the Kohn’s theorem in Ge/Si quantum dots with hole gas: Theory and experiment.” Nanomaterials 9.1 (2019): 56.
  * [35] Blatt, Rainer, and David Wineland. ”Entangled states of trapped atomic ions.” Nature 453.7198 (2008): 1008-1015.
  * [36] Brown, K. R., Chiaverini, J., Sage, J. M., Häffner, H., ”Materials challenges for trapped-ion quantum computers”, Nature Reviews Materials, 6.10 (2021): 892-905
  * [37] Busch, Thomas, et al. ”Two cold atoms in a harmonic trap.” Foundations of Physics 28.4 (1998): 549-559.
  * [38] Brey, L., N. F. Johnson, and B. I. Halperin. ”Optical and magneto-optical absorption in parabolic quantum wells.” Physical Review B 40.15 (1989): 10647.
  * [39] Deng, X-L., Diego Porras, and J. Ignacio Cirac. ”Quantum phases of interacting phonons in ion traps.” Physical Review A 77.3 (2008): 033403.
  * [40] DinhDuy Vu, S.Das Sarma, ”One-dimensional few-electron effective Wigner crystal in quantum and classical regimes”, Physical Review B 101.12 (2020): 125113.
  * [41] Loudon, Rodney. ”One-dimensional hydrogen atom.” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472.2185 (2016): 20150534.
  * [42] M.Taut, K.Pernal, J.Cioslowski, V.Staemmier, ”Three electrons in a harmonic oscillator potential: Pairs versus single particles”, The Journal of chemical physics 118.11 (2003): 4861-4871.
  * [43] M.Taut, ”Two electrons in an external oscillator potential: Particular analytic solutions of a Coulomb correlation problem”,Physical Review A 48.5 (1993): 3561.
  * [44] Gharashi, Seyed Ebrahim, and Doerte Blume. ”Correlations of the upper branch of 1D harmonically trapped two-component Fermi gases.” Physical review letters 111.4 (2013): 045302.
  * [45] Tomasz Grining, Michał Tomza, Michał Lesiuk, Michał Przybytek, Monika Musiał, Pietro Massignan, Maciej Lewenstein, and Robert Moszynsk, ”Many interacting fermions in a one-dimensional harmonic trap: a quantum-chemical treatment”. New Journal of Physics 17.11 (2015): 115001.
  * [46] Zhong-Qi Ma,_Group Theory for physicists_ , first edition, chapter 6, 2007.
  * [47] Tomasz Sowinski, Tobias Grass, Omjyoti Dutta, and Maciej Lewenstein, ”Few interacting fermions in a one-dimensional harmonic trap”, Physical Review A 88.3 (2013): 033607.
  * [48] Sebastian Bauch, David Hochstuhl, Karsten Balzer, and Michael Bonitz, ”Quantum Breathing mode of interacting particles in harmonic traps”, Journal of Physics: Conference Series Vol. 220, No. 1. IOP Publishing, 2010: p. 012013.
  * [49] C.Henning, K.Fujioka, P.Ludwig, A.Piel, A. Melzer, and M.Bonitz, ”Existence and Vanishing of the Breathing Mode in Strongly Correlated Finite Systems”, Physical review letters 101.4 (2008): 045002.
  * [50] Abraham, Jan Willem, et al. ”Quantum breathing mode of interacting particles in a one-dimensional harmonic trap.” Physical Review B 86.12 (2012): 125112.
  * [51] Hashemi, Hossein, Servaas Kokkelmans, and Thomas Secker, ”N interacting electrons in a harmonic trap”, Eindhoven University of Technology (2017).
  * [52] Kashuba, Oleksiy, et al. ”Counting interacting electrons in one dimension.” arXiv preprint arXiv:2305.15906 (2023).
  * [53] Räsänen, Esa, et al. ”Wigner molecules in polygonal quantum dots: A density-functional study.” Physical Review B 67.3 (2003): 035326.
  * [54] I. Shapir and A. Hamo and S. Pecker and C. P. Moca and Ö. Legeza and G. Zarand and S. Ilani , ”Imaging the electronic Wigner crystal in one dimension”, Science 364.6443 (2019): 870-875.
  * [55] N. Tancogne-Dejean, M. J. T. Oliveira, X. Andrade, H. Appel, C. H. Borca, G. Le Breton, F. Buchholz, A. Castro, S. Corni, A. A. Correa, U. De Giovannini, A. Delgado, F. G. Eich, J. Flick, G. Gil, A. Gomez, N. Helbig, H. Hübener, R. Jestädt, J. Jornet-Somoza, A. H. Larsen, I. V. Lebedeva, M. Lüders, M. A. L. Marques, S. T. Ohlmann, S. Pipolo, M. Rampp, C. A. Rozzi, D. A. Strubbe, S. A. Sato, C. Schäfer, I. Theophilou, A. Welden, A. Rubio, ”Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems”, The Journal of Chemical Physics 152.12 (2020): 124119
  * [56] Liming Guan, Shu Chen, Yupeng Wang, Zhong-Qi Ma, ”Exact Solution for Infinitely Strongly Interacting Fermi Gases in Tight Waveguides”, Physical review letters 102.16 (2009): 160402.
  * [57] Ziesche, Paul. ”Correlation strength and information entropy.” International Journal of Quantum Chemistry 56.4 (1995): 363-369.
  * [58] Lathiotakis, Nektarios N. and Helbig, Nicole and Rubio, Angel and Gidopoulos, Nikitas I. ”Local reduced-density-matrix-functional theory: Incorporating static correlation effects in Kohn-Sham equations”, Phys. Rev. A 90.3 (2014):032511.
  * [59] Iwata, J-I., K. Yabana, and G. F. Bertsch. ”Real-space computation of dynamic hyperpolarizabilities.” The Journal of Chemical Physics 115.19 (2001): 8773-8783.
  * [60] Quinteiro, Guillermo Federico, D. E. Reiter, and T. Kuhn. ”Formulation of the twisted-light–matter interaction at the phase singularity: The twisted-light gauge.” Physical Review A 91.3 (2015): 033808.
  * [61] Rosen, Guillermo F. Quinteiro, Pablo I. Tamborenea, and Tilmann Kuhn. ”Interplay between optical vortices and condensed matter.” Reviews of Modern Physics 94.3 (2022): 035003.
  * [62] Prinz, Eva, et al. ”Orbital angular momentum in nanoplasmonic vortices.” ACS Photonics 10.2 (2023): 340-367.
  * [63] Rivera, Nicholas, et al. ”Shrinking light to allow forbidden transitions on the atomic scale.” Science 353.6296 (2016): 263-269.
  * [64] Schmiegelow, Christian T., et al. ”Transfer of optical orbital angular momentum to a bound electron.” Nature communications 7.1 (2016): 12998.
  * [65] Wineland DJ, Monroe C, Itano WM, Leibfried D, King BE, Meekhof DM., ”Experimental Issues in Coherent Quantum-State Manipulation of Trapped Atomic Ions”, J Res Natl Inst Stand Technol, 103.3 (1998):259-328
  * [66] Quinteiro, Guillermo Federico. ”Below-bandgap excitation of bulk semiconductors by twisted light.” Europhysics Letters 91.2 (2010): 27002.
  * [67] Quinteiro, Guillermo Federico, et al. ”Reexamination of Bessel beams: A generalized scheme to derive optical vortices.” Physical Review A 99.2 (2019): 023845.
  * [68] Konstantin Y. Bliokh and Franco Nori, ”Transverse and longitudinal angular momenta of light”, Physics Reports, 592 (2015):0370-1573,
  * [69] Allen, L. and Barnett, S.M. and Padgett, M.J., _Optical Angular Momentum_ , Taylor & Francis, 2003
  * [70] Andrews, David L. _Structured light and its applications: An introduction to phase-structured beams and nanoscale optical forces_. Academic press, 2011.
  * [71] Haldane, F. D. M. ”’Luttinger liquid theory’of one-dimensional quantum fluids. I. Properties of the Luttinger model and their extension to the general 1D interacting spinless Fermi gas.” Journal of Physics C: Solid State Physics 14.19 (1981): 2585.
  * [72] Cazalilla, M. A. ”Bosonizing one-dimensional cold atomic gases.” Journal of Physics B: Atomic, Molecular and Optical Physics 37.7 (2004): S1.
  * [73] Jin, Tony, Paola Ruggiero, and Thierry Giamarchi. ”Bosonization of the interacting Su-Schrieffer-Heeger model.” Physical Review B 107.20 (2023): L201111.
  * [74] Giamarchi, Thierry. _Quantum physics in one dimension_. Vol. 121. Clarendon press, 2003.
  * [75] Fajardo Montenegro, José Luis. ”Bosonization method in the study of the fulde-ferrell-larkin-ovchinnikov phase.”
  * [76] L. Sonderhouse, C. Sanner, R. B. Hutson, A. Goban, T. Bilitewski, L. Yan, W. R. Milner, A. M. Rey, and J. Ye, ”Thermodynamics of a deeply degenerate SU (N)-symmetric Fermi gas”. Nature Physics 16.12 (2020): 1216-1221.
  * [77] Song, Bo, et al. ”Evidence for Bosonization in a three-dimensional gas of SU (N) fermions.” Physical Review X 10.4 (2020): 041053.
  * [78] S. Weinberg, _The Quantum Theory of Fields_ (University of Texas, Austin, 2000).
  * [79] Ernst, Jan Ole, and Felix Tennie. ”Mode Entanglement in Fermionic and Bosonic Harmonium.” arXiv preprint arXiv:2211.09647 (2022).
  * [80] F. Serwane and G. Zürn and T. Lompe and T. B. Ottenstein and A. N. Wenz and S. Jochim , ”Deterministic Preparation of a Tunable Few-Fermion System”, Science 332.6027 (2011):336-338
  * [81] Wagner, Mathias, A. V. Chaplik, and U. Merkt. ”Quadrupole excitations of quantum dots.” Physical Review B 51.19 (1995): 13817.
  * [82] R. MCWeeny, _Methods of molecular quantum mechanics_ , second edition, chapter 4, page 100, 1992.
  * [83] N.Helbig, J.I.Fuks, M.Casula, M.J.Verstraete, M.A.L.Marques, I.V.Tokatly, A.Rubio, ”Density functional theory beyond the linear regime: Validating an adiabatic local density aproximation”, Physical Review A 83.3 (2011): 032503.