mathematics

Article
Static and Dynamic Properties of a Few Spin 1/2
Interacting Fermions Trapped in a Harmonic Potential
Abel Rojo-Francàs, Artur Polls and Bruno Juliá-Díaz *
Departament de Física Quàntica i Astrofísica, Facultat de Física, and Institut de Ciències del Cosmos (ICCUB),
Universitat de Barcelona, E–08028 Barcelona, Spain; abel.9.1995@gmail.com (A.R.-F.); artur@fqa.ub.edu (A.P.)
* Correspondence: bruno@fqa.ub.edu




Received: 21 May 2020; Accepted: 15 July 2020; Published: 21 July 2020


Abstract: We provide a detailed study of the properties of a few interacting spin 1/2 fermions trapped
in a one-dimensional harmonic oscillator potential. The interaction is assumed to be well represented
by a contact delta potential. Numerical results obtained by means of direct diagonalization techniques
are combined with analytical expressions for both the non-interacting and strongly interacting regime.
The N = 2 case is used to benchmark our numerical techniques with the known exact solution
of the problem. After a detailed description of the numerical methods, in a tutorial-like manner,
we present the static properties of the system for N = 2, 3, 4 and 5 particles, e.g., low-energy spectrum,
one-body density matrix, ground-state densities. Then, we consider dynamical properties of the
system exploring first the excitation of the breathing mode, using the dynamical structure function
and corresponding sum-rules, and then a sudden quench of the interaction strength.

Keywords: few-body systems; one-dimensional trap; trapped atoms; direct diagonalization

1. Introduction
The theoretical study of one-dimensional systems has always attracted a lot of attention. Very early
in the development of quantum mechanics it was recognized that reducing the dimensionality of
the systems enhances the quantum effects and gives origin to a plethora of interesting phenomena,
see for instance Ref. [1] and references therein. After the pioneering work of Tonks, who investigated
the equation of state of hard-rods, hard disks and hard-spheres in one, two and three dimensions,
respectively [2], Girardeau established the map between impenetrable bosons and free fermions in
one dimension [3]. However, it took a long time until the experimental realization of these systems
became a reality [4–6]. This experimental breakthrough took place in the context of the ultracold gases
that have developed a frenetic experimental and theoretical activity during the last years achieving
unbelievable possibilities to control the geometry, the interactions and the number of particles for
many types of setups [7–10]. More recently, the groundbreaking experiments of Jochim’s group in
Heildelberg, have opened new theoretical challenges to study one-dimensional fermionic systems.
They have been able to precisely control the number of atoms and the strength of the interactions [11].
These experiments have provided deep insight in the expected phenomena of fermionization of two
distinguishable fermions [12], i.e., with different third spin component. In addition, they have been
able to change the number of particles, going from few to many particles, and being able to study the
appearance of correlations in the many-body wave function. In fact, by adding particles to the system
one by one, it was shown how the many-body correlated system is built [13]. These experiments used
6 Li atoms and took advantage of the Feshbach resonance mechanism to modulate the interactomic

interaction [14]. Consequently, many theoretical problems have been addressed and few- and
many-body techniques for ab-initio descriptions used in other fields have been specially tailored
for these many-body one-dimensional systems. The impurity problem, the presence of pairing

Mathematics 2020, 8, 1196; doi:10.3390/math8071196 www.mdpi.com/journal/mathematics

Mathematics 2020, 8, 1196 2 of 38

phenomena in few-fermion systems [15], the behaviour of two fermions in a double-well potential [16],
the realization of an antiferromagnetic spin chain of few cold atoms [17], all belong to the long list of
theoretical challenges that the experimental advances are offering to theoreticians. They are looking
at those systems as a very versatile laboratory where they apply ab-initio techniques, without the
complexity of the inter-particle interactions that appear in other fields, as for instance in nuclear
physics [18,19]. Reciprocally, the technological control is such that ultracold atomic systems can also
be used as quantum simulators to solve very intricate condensed-matter problems [9]. Although we
will not pay attention to them in this paper it is worth mentioning that trapping and control of the
interactions have been extended also to mixtures of different nature and statistics: mixtures of bosons,
and bosons and fermions [20–25].
In this work, we consider a system of a few spin-1/2 fermions trapped in a one-dimensional
harmonic trap. The interaction between the fermions is modeled by a contact potential. In the literature
there are some studies under similar conditions [25–33]. We pay attention both to the ground state and
energy spectrum of the system and also to the dynamical properties of few-fermion systems depending
on the interaction strength and the number of particles [34,35].
The structure of this work is the following. In Section 2, we present the theoretical description of
the fermionic system, composed by a fixed number of atoms with spin up and down. In particular,
we write down the full Hamiltonian both in first and second quantized form. Some limiting scenarios
are considered, e.g., both zero and strong interaction. In Section 3, we describe the numerical tools used
to study the ground state and low-excited states of the fermionic mixture. The main method employed
is direct diagonalization of the Hamiltonian matrix in a truncated Hilbert space. In Section 4, we study
the ground state of the system for different number of fermions and polarizations. In Section 5 we
turn our attention to the lower part of the spectrum, computing the spectra for the same few-fermion
systems considered. In Section 6, we compute the response of the system to two different external
perturbations. Finally in Section 7 we provide the main conclusions of our work.

2. Theoretical Approach
In this section we introduce the theoretical tools to describe a system composed of a few fermions
trapped in a one dimensional harmonic oscillator potential. The goal is to provide a self-contained
explanation of the methods employed, detailing the numerical subtleties encountered.

2.1. The Hamiltonian of the System

The Hamiltonian of the system can be split in two parts, the non-interacting single-particle part,
which describes particles trapped in a harmonic oscillator potential, and the two-body interaction
part, which describes the interactions between the fermions. The Hamiltonian, in first quantized form,
can be written as [33,35]
N
H = Hho + ∑ Vint ( xi − x j ) , (1)
i< j

where Hho = ∑i hho (i ) is the sum of single-particle harmonic oscillator Hamiltonians. The Hamiltonian,
hho , is a one-body operator, which includes the kinetic- and the harmonic-potential energy,

h̄2 ∂2 1
hho = − + mω 2 x2 , (2)
2m ∂x2 2
where m is the particle mass and ω is the harmonic oscillator frequency. The eigenvalues of this
Hamiltonian are well known,
 
1
en = n + h̄ω n = 0, 1, 2, ... , (3)
2
Mathematics 2020, 8, 1196 3 of 38

and the corresponding eigenfunctions are written as,

 mω 1/4 r 
1 − mωx
2 mω
Φn ( x ) = √ e 2h̄ Hn x , (4)
2n n! πh̄ h̄

where Hn ( x ) is the n-th Hermite polynomial. These single-particle wave functions will be used to
build the many-body basis. All these wave functions have a similar structure: a Gaussian function
multiplied by a Hermite polynomial with a normalization factor, that depends on n.
We consider identical spin 1/2 fermions, i.e., with two possible spin states: |1/2, mi, m =
1/2, −1/2, where m is the spin projection on the z-axis. Sometimes, these two spin states are denoted
by: |↑i and |↓i, respectively. The fermions are assumed to interact via a contact spin-independent
delta potential [36]. Therefore two fermions interact only if they are at the same position. However,
the many-body wave function of a system of identical fermions should be antisymmetric, preventing
two fermions with the same spin from being at the same position. Therefore, fermions with the same
spin projection do not interact, and thus the contact interaction is only acting between fermions with
different spin projection.
In any case, the total wave function including the spin degree of freedom should be antisymmetric.
This requirement allows the system to have two particles with different spin projection in the
same position.
The interaction term of Equation (1) is a contact interaction, which is expressed as [36]

Vint ( xi − x j ) = gδ( xi − x j ) , (5)

where g characterizes the strength of the interaction and δ is a Dirac delta function.
Taking into account the identity [34]

N
1
∑ xi2 = NX2 + N ∑(xi − x j )2 , (6)
i =1 i< j

1
where X = N ∑i xi is the center of mass coordinate, the Hamiltonian can be split in two pieces,

H = HCM + Hr , (7)

where HCM in harmonic oscillator units reads:

1 d2 N
HCM = − 2
+ X2 , (8)
2N dX 2
and governs the center-of-mass motion. Hr affects only the relative coordinates and is
translationally invariant.
One way to create a one-dimensional trap is using a cigar-shaped trapping potential, with the
transverse trap frequency ω⊥ (in the radial direction) much larger than the trap frequency ω in the
axial direction. In this situation, the coupling constant g is related to the one-dimensional scattering
length (a1d ) as [36,37],
h̄2
g=− . (9)
ma1d
From the relation between the one-dimensional scattering length and the three-dimensional
scattering length [38], we can write the coupling constant as

2h̄2 a3d 1
g= √ , (10)
ma2⊥ 1 − |ξ (1/2)| a3d /( 2a⊥ )
Mathematics 2020, 8, 1196 4 of 38

√
with a⊥ = h̄/mω⊥ and ξ is the Riemann zeta function.
In a trap under these conditions, the system can be treated as a one-dimensional system. As the
trapping in the transverse dimension is very strong, all particles occupy the lowest state of the
transverse harmonic oscillator and the physics takes place in the longitudinal direction [38].
For convenience we use the harmonic oscillator units, in which the energy is measured
√ p in h̄ω
units, the length in units of h̄/mω and the coupling constant g is expressed in units of ωh̄3 /m.
From this moment on, all magnitudes will be expressed in these units.

2.2. Non-Interacting and Infinite Interaction Limits

In general, it is not possible to solve analytically the Schrödinger equation for an arbitrary value
of the interaction strength. However, there are two limits, i.e., the non-interacting and the infinite
interaction cases, in which one has explicit analytical solutions [39,40].
For these two limits, one can easily determine the energy and the density profile of the ground
state of the system, and explicitly write the wave function.

2.2.1. Non-Interacting Case: Fermi Gas

In the non-interacting case, the system is in a state called Fermi gas. In these conditions,
the behavior of the system is given by the single-particle states of the harmonic oscillator with the
restriction of the Pauli principle that does not allow two fermions with the same spin projection in the
same harmonic oscillator state. More specifically, in the ground state of a system with Nd particles of
spin down and Nu particles of spin up, the particles occupy the lowest single-particle energy states, Nd
and Nu , respectively. Therefore, the energy of the ground state is

Nd −1  Nu −1  Nd2 + Nu2
 
1 1
Egs ( g = 0) = ∑ n+
2
+ ∑ n+
2
=
2
. (11)
n =0 n =0

Actually the wave function of the system is given by a Slater determinant built with the lowest Nu and
Nd single-particle wave functions with spin up and down, respectively. The total third component of
the spin of this wave function is M = Nu −2
Nd
and the total spin is the minimum S compatible with M,
S = M [41]. For a non-polarized system, Nu = Nd , then M = 0 and S = 0.
The density profile associated to this wave function is the sum of the probabilities of finding a
particle in each occupied state in the position x,

Nd −1 Nu −1
ρ( x ) = ∑ |Φn ( x )|2 + ∑ |Φn ( x )|2 , (12)
n =0 n =0

where Φn ( x ) are the 1D harmonic oscillator wave functions, Equation (4). The density profile is
normalized to the number of particles.

2.2.2. Infinite Interaction Case

In the infinite interaction limit, the system experiences a phenomenon reminiscent of the
fermionization process in Bose systems [2,3]. The ground state of the Hamiltonian in the infinite
interaction limit for a fully- polarized state is a Slater determinant of the N first harmonic oscillator
states, Ψ A . For a non fully-polarized case the ground state many-fold is degenerate, with degeneration
D = N!/( Nu !Nd !) [33] and total third spin component M = ( Nu − Nd )/2. The states can be written
as [40],
D
Ψ= ∑ ak θ (xPk (1) , ..., xPk ( N ) )Ψ A (x1 , ..., x N ) , (13)
k =1

where θ ( x1 , ..., x N ) = 1 when x1 < x2 < ... < x N and zero otherwise, Pk are the D permutations of the
N components and the ak are the coefficients for each permutation. In this expression we have fixed
Mathematics 2020, 8, 1196 5 of 38

the first Nu particles to be spin projection up and the remaining Nd particles to have spin projection
down. Their energy depends only on the single-particle states present in Ψ A ,

N −1 
( Nd + Nu )2

1
Egs ( g → ∞) = ∑ n+
2
=
2
. (14)
n =0

The density profile associated to these wave functions is the sum of the probabilities of finding a
particle in each occupied state in the position x [42],

N −1
ρ( x ) = ∑ |Φn ( x )|2 , (15)
n =0

where Φn ( x ) are the 1D harmonic-oscillator wave functions, Equation (4).

2.3. Second Quantization

When working with several particles, it is useful to use the second-quantization formalism. One of
the main ingredients of this formalism is the many-body Fock space. In order to define the Fock space,
we need to have a single-particle basis, such as the harmonic oscillator states. The basis of the Fock
space is constructed by indicating how many particles occupy each single-particle state. In addition,
the antisymmetry rules are implemented by the anticommutation rules fulfilled by the creation and
annihilation operators [43].
A particular state is the vacuum state denoted by |0i. This vector represents a state without
particles, i.e., all its occupation numbers are zero. This state has norm equal to 1.

2.3.1. Creation and Annihilation Operators

The main tool in second quantization are the creation ( ai† ) and annihilation ( ai ) operators,
which act on the Fock space. Any observable, which is represented by an operator, can be expressed in
terms of them [43]. The action of a creation operator on the vacuum state, denoted as ai† |0i, creates a
particle in the state i. In general, acting with a creation operator on a Fock space vector, creates a
particle in the i-th state, whenever possible. Or in other words, it increases the occupation number of
this single-particle state by one. On the other hand, the action of an annihilation operator (ai ) on a state
destroys a particle in the i-th state, whenever possible. In other words, it decreases the occupation
number of the state by one.
In addition, a Fock space vector with N particles can be created by the action of N creation
operators on the vacuum space, as follows,

|n1 , n2 , n3 , ..., nm i = ( a1† )n1 ( a2† )n2 ( a3† )n3 ...( a†m )nm |0i (16)

with N = ∑i ni .
In the case of fermions, to satisfy the Pauli principle, i.e., the antisymmetry of the wave function,
the creation and annihilation operators should fulfill the following anti-commutation relations
n o
ai† , a j = δij
n o (17)
ai† , a†j

= ai , a j = 0 .

Therefore the occupation number of the Fock basis can take the values ni = 0 or ni = 1. Finally,
the action of the creation and annihilation operators can be written as

ai† |n1 , n2 , ...ni−1 , ni , ...i = (−1) Ni (1 − ni ) |n1 , n2 , ...ni−1 , (1 − ni ), ...i

(18)
ai |n1 , n2 , ...ni−1 , ni , ...i = (−1) Ni ni |n1 , n2 , ...ni−1 , (1 − ni ), ...i ,
Mathematics 2020, 8, 1196 6 of 38

where Ni is the number of occupied states with index lower than i,

i −1
Ni = ∑ nk . (19)
k =1

The phase (−1) Ni is due to the anti-commutation relations. Notice that the factor in front of the state
becomes zero if one tries to put two particles in the same single-particle state or one tries to annihilate
a particle in a single-particle state which is not occupied.

2.3.2. Fock Space

The harmonic oscillator eigenfunctions are used to build the many-body basis of the Fock space.
In addition, the spin degrees of freedom are also incorporated to the single particle wave functions
which are defined as: ϕn,m ( x, s) = Φn ( x )χm (s), where Φn ( x ) is the harmonic oscillator wave function
of the level n and χm is the spin wave function, where m =↑ (↓) is the spin projection.
With these single-particle states, we build the many-body basis of the Fock space:

n0↓ , n0↑ , n1↓ , n1↑ , ... , (20)

where ni,m indicates the number of particles in the ϕi,m single-particle state.

2.3.3. Operators in Second Quantization

In second quantization, operators can be expressed in terms of creation and annihilation operators.
A general one-body operator Ô is expressed as

Ô = ∑ hi | Ô | ji ai† a j , (21)
i,j

where |i i and | ji are single-particle states, hi | Ô | ji = ϕi∗ ( x )O( x ) ϕ j ( x )dx is the one-body matrix
R

element and ϕi ( x ) are the single-particle wave functions used to build the Fock space. Notice that the
sub-index i stands for all the quantum numbers necessary to specify the single-particle state.
A general two-body operator V̂ is expressed as

1 1
V̂ = ∑
2 i,j,k,l
hi, j| V̂ |k, l i ai† a†j al ak = ∑ vij,kl ai† a†j al ak ,
2 i,j,k,l
(22)

where vij,kl is the two-body matrix element defined as

Z
vij,kl = ϕi∗ ( x1 ) ϕ∗j ( x2 )V ( x1 , x2 ) ϕk ( x1 ) ϕl ( x2 )dx1 dx2 . (23)

2.3.4. The Hamiltonian in Second Quantization

As already mentioned, the Hamiltonian describing the harmonic oscillator is a one-body operator.
On the other hand, the interaction term Equation (1) is a two-body operator. In our case, choosing the
harmonic-oscillator eigenfunctions as the single-particle basis, the Hamiltonian part corresponding to
the harmonic oscillator is diagonal, with eigenvalues ei = i + 1/2. Therefore in second quantization
it reads,
ĥho = ∑ ei ai† ai = ∑ ei n̂i , (24)
i i
Mathematics 2020, 8, 1196 7 of 38

where n̂i = ai† ai is the number operator associated to the single-particle state |i i. For the interaction
term, the two-body matrix elements are expressed as
Z
∗
vij,kl = dx1 dx2 ϕi,m i
( x1 , s1 ) ϕ∗j,m j ( x2 , s2 ) gδ( x1 − x2 ) ϕk,mk ( x1 , sk ) ϕl,ml ( x2 , sl )
Z D E
= g dxΦi∗ ( x )Φ∗j ( x )Φk ( x )Φl ( x ) χmi χm j χmk χml (25)
Z
= gδmi ,mk δm j ,ml dxΦi∗ ( x )Φ∗j ( x )Φk ( x )Φl ( x ) .

Notice that the interaction

D does not affect
E the spin of the particles, and we have used the orthogonality
of the spin functions: χmi χm j χmk χml = δmi ,mk δm j ,ml . The calculation of the integral can be found in
Section 3.2. Notice that the indices labeling the single-particle states run over the harmonic oscillator
wave functions and the spin projections. Finally, the Hamiltonian reads

1
H= ∑ n̂i ei + 2 ∑ vij,kl ai† a†j al ak . (26)
i ijkl

3. Numerical Methods
In this chapter we describe the method named direct diagonalization to find an approximate
solution of the many-body Schrödinger equation. The technical aspects on how to construct the
matrix Hamiltonian, to choose a suitable many-body basis of the Fock space, and how to calculate the
two-body matrix elements are discussed in detail. Finally, we use the two-particle system, that has an
analytical solution [44], as a benchmark of our numerical procedure.

3.1. Direct Diagonalization

In order to study the ground state and the low-excited states of a system, one has to solve the
many-body Schrödinger equation. With this objective, we use a direct diagonalization in a truncated
basis space. This method needs to built the Hamiltonian matrix in an appropriate subspace and
diagonalize it. In order to obtain accurate results we need to use a large basis, generating large
matrices. To diagonalize these large matrices, we use a Lanczos method implemented by the ARPACK
package [45]. This method allows us to diagonalize large matrices, and obtain the lowest eigenvalues
with high accuracy. Hence, the direct diagonalization method allows one to obtain the ground state
and the low-energy spectrum and the corresponding states. Other strategies can also provide the
lowest part of the spectrum, such as the Multi Configuration Time Dependent Hartree (MCTDH) [46],
used in Ref. [27] for few-fermion systems. In contrast, methods such as Monte Carlo [47], or Density
Matrix Renormalization Group (DMRG) [48], are usually used to compute ground state properties.
We consider systems with 2 to 5 fermions. These particles have spin 1/2, therefore for N particles
we can have N + 1 possible total spin projections. As the Hamiltonian commutes with the total third
spin component, the Hamiltonian is built in boxes with well defined M. Thus, it is possible to treat
each total spin projection independently [33]. The simplest case is when all fermions have spin up,
then M = N/2, and the wave function factorizes in an antisymmetric spatial wave function and a
symmetric spin wave function χ(S = N/2, M = N/2) for N particles. The antisymmetric spatial
wave function is a Slater determinant built with the lowest single-particle states. The antisymmetry
of the wave function prevents two particles from being in the same position and the particles in
this wave function do not feel the contact interaction. Therefore, its energy is independent of the
interaction strength.
The other trivial cases are the negative M values because due to the spin symmetry of
the Hamiltonian, the properties of the states do not depend on the sign of the spin projection.
The diagonalization of a Hamiltonian box with a given M would provide the same eigenvalues
that the box with total spin projection − M. Therefore we concentrate in the study of cases with M ≥ 0.
Mathematics 2020, 8, 1196 8 of 38

Basis Truncation
To obtain the exact results, when using a diagonalization method, we would need to use a
complete basis. However this is not possible because of the infinite dimension of the Hilbert space.
Therefore we are forced to diagonalize the Hamiltonian in a finite subspace. Usually, this procedure
does not provide the exact eigenvalue. However, the lowest-energy obtained by diagonalization is still
an upper-bound to the ground-state energy. Note this is in contrast to exact diagonalization problems
in small and discrete spaces, modelling finite optical lattices [49], where the method provides the exact
values within machine precision.
A common way to construct a finite many-body basis is considering a finite number of
single-particle states, usually, those with the lowest energy [49]. In our case, we diagonalize the
Hamiltonian in a subspace with well-defined total third spin component. We start considering the
lowest NM single-particle states and take into account an energy constraint in the construction of
the many-body basis: the energy of the many-body basis states, which is given by the sum of the
single-particle energies, should be smaller or equal than a fixed energy Emax , this procedure is described
in [50,51].
This maximum energy depends on the number of particles and spin configuration. Notice that we
consider always M > 0 therefore, Nu > Nd . The maximum energy, Emax , corresponds to the energy of a
non-interacting many-body state built as follows: one spin-up particle in the maximum single-particle
energy state and the remaining (Nu − 1) spin-up particles in the lowest Nu − 1 single-particle states.
On the other hand the Nd spin-down particles are located in the lowest Nd single-particle states.
Therefore, the maximum energy considered in the construction of the many-body basis is

( Nd )2 +( Nu −1)2
       
N −1
Emax = ∑n=d 0 n+ 1
2 + ∑nN=u −
0
2
n + 12 + ( NM − 1) + 12 = 2 + NM − 12 , (27)

where NM is the number of single-particle states used.

One can see in Table 1 that the dimension of the many-body basis is strongly reduced if the energy
restriction discussed above is taken into account. Nevertheless, the reduction of the size of the space
considered does not affect the quality of the results in the low-energy regime of the spectrum that
we are exploring in this paper. The dimension of the many-body basis in the case without energy
restriction is
( NM ) ! ( NM ) !
, (28)
( NM − Nu )!Nu ! ( NM − Nd )!Nd !
which grows very rapidly.

Table 1. Number of single-particle states used in the construction of the many-body basis states in the
second column. The number of many-body basis states, with and without energy restriction, are shown
in the third and fourth columns, respectively.

Number of Number of Many-Body Basis States

Single-Particle States
with Energy Restriction without Restriction
2 particles, M = 0 100 5050 10,000
3 particles, M = 1/2 50 10,725 61,250
4 particles, M = 0 40 30,800 608,400
4 particles, M = 1 40 19,530 395,200
5 particles, M = 1/2 30 22,923 1,766,100
5 particles, M = 3/2 30 11,349 822,150

In the table we have omitted the configurations with a maximum value of M, which correspond
to trivial cases, with only one configuration built as the product of a symmetric spin function,
Mathematics 2020, 8, 1196 9 of 38

with all spins parallel and an antisymmetric spatial function built with the lowest single-particle
energy wave functions. As has been said previously, these configurations do not feel the effects of a
contact interaction.

3.2. The Two-Body Matrix Elements of the Interaction

The evaluation of the two-body interaction on the many-body Fock basis, requires the calculation
of two-body interaction matrix elements, Equation (25). These matrix elements contain an integral
with four wave functions,
Z
vij,kl = gδmi ,mk δm j ,ml dxΦi∗ ( x )Φ∗j ( x )Φk ( x )Φl ( x )
(29)
= gδmi ,mk δm j ,ml Iijkl .

The integral, that we label as Iabcd can be solved analytically, because the wave functions are the
harmonic oscillator wave functions Equation (4),
Z ∞
Iabcd = Φ a ( x )Φb ( x )Φc ( x )Φd ( x )dx
−∞
Z ∞ (30)
1 2
= √ e−2x Ha ( x ) Hb ( x ) Hc ( x ) Hd ( x )dx .
π 2a+b+c+d a!b!c!d! ∞

Let us consider, Z ∞
0 2
Iabcd = e−2x Ha ( x ) Hb ( x ) Hc ( x ) Hd ( x )dx . (31)
−∞
We notice that this integral has to be zero if (a + b + c + d) is an odd number. This is because the
Φi wave functions have well-defined parity. Therefore if the multiplication of the four functions is
odd, the integral value is zero. This is why we only calculate the integrals with (a + b + c + d) even.
Using the properties [52]
Z ∞
2
e−2x Ha ( x ) Hb ( x ) Hc ( x ) =
−∞
(32)
2(a+b+c−1)/2 a+b−c+1 a−b+c+1 −a + b + c + 1
     
= Γ Γ Γ ,
π 2 2 2

and
n
m! Hm−n+2r ( x )
Hm ( x ) Hn ( x ) = 2n n! ∑ (n − r )!(m − n + r )! 2r r!
,n ≤ m. (33)
r =0
0
We can express the integral Iabcd as

0
d
c! 1 2(a+b+c−d+2r−1)/2
Iabcd =2d d! ∑
(d − r )!(c − d + r )! 2r r! π
r =0
a+b−c+d+1 a−b+c−d+1
   
×Γ −r Γ +r (34)
2 2
−a + b + c − d + 1
 
×Γ +r .
2
Mathematics 2020, 8, 1196 10 of 38

Notice that the arguments of the gamma functions (Γ) are half-integers. Because (a + b + c + d) is even,
then when changing any sign, the result will remain even, and the gamma functions read,

√ (2n)! 
 
1
Γ n+ = π n 

2 2 n! 
,n ≥ 0. (35)
√ (−4)n n! 
 
1
Γ −n = π 

2 (2n)!

Finally, the integral Iabcd can be expressed as

r
d
a+b−c+d+1
 
1 c!d! 1
Iabcd =
π2 2
√ ∑ r!(d − r)!(c − d + r)!
a!b!
Γ
2
− r
r =0 (36)
a−b+c−d+1 −a + b + c − d + 1
   
×Γ +r Γ +r .
2 2

The presence of factorials in the expressions of the harmonic oscillator wave functions can cause
fake overflows in the calculation. One possible solution to this problem is by taking logarithms of the
expressions to be evaluated. As shown before, a suitable way to write the integrals entering in the
two-body matrix elements is
r
d
a+b−c+d+1
 
1c!d! 1
Iabcd = ∑ 2 2 a!b! r!(d − r)!(c − d + r)!
√ Γ
2
− r
r =0 π (37)
a−b+c−d+1 −a + b + c − d + 1
   
×Γ +r Γ +r .
2 2

This expression can be written as I = ∑rd f ( a, b, c, d, r ), where f is a well-behaved function.

Furthermore, we can write this result as

d
Iabcd = ∑ exp{log [ f (a, b, c, d, r)]} , (38)
r =0

where f must be positive.

The logarithm log ( f ) is given by

1
log ( f ) = − log (2) − 2 log (π )
2
1
+ (log (c!) + log (d!) − log ( a!) − log (b!))
2
− log ((d − r )!) − log ((c − d + r )!) − log (r!)
a+b−c+d+1
  
(39)
+ log Γ −r
2
a−b+c−d+1
  
+ log Γ +r
2
−a + b + c − d + 1
  
+ log Γ +r .
2

Notice that due to the presence of the gamma functions, the previous expression could require the
evaluation of a logarithm of a negative quantity. For this reason, we calculate instead log (|Γ(n)|), and
compute the associated phase separately. This phase is introduced in the final expression after the
Mathematics 2020, 8, 1196 11 of 38

exponentiation. Thus, the final expression of the integral taking into account the possible negative
values of gamma functions reads

d
Iabcd = ∑ p exp{log [ f (a, b, c, d, r)]} , (40)
r =0

where p is the phase, given by p = ∏3n=1 pn , with pn being the phase generated by each
gamma function. The logarithm of the gamma functions and the associated phases are calculated as
  
1 1

log Γ n + = log (π ) + log ((2n)!) − n log (2) − log ((n)!) , p = 1 

2 2 
   ,n ≥ 0. (41)
1 1
log Γ = log (π ) + 2n log (2) + log ((n)!) − log ((2n)!) , p = (−1)n 

−n 
2 2

The logarithm of a factorial is calculated as:

N
log ( N!) = ∑ log (n) . (42)
n =1

3.3. A Benchmark for the Two-Particle Case

In general, it is not possible to solve exactly the energy spectrum. In our case, to determine
the energy spectrum we have to diagonalize the Hamiltonian in a large Hilbert space by using
sophisticated numerical techniques. However the case of two particles has been analytically solved in
the literature [44].

3.3.1. Theoretical Spectrum for Two Particles

For two fermions with opposite spin, M = 0, the energies of the relative motion are obtained by
solving the transcendental equation [44]

Γ(− Er /2 + 3/4) g
= − 3/2 , (43)
Γ(− Er /2 + 1/4) 2

where Γ are gamma functions, Er is the energy of the relative system and g is the interaction strength.
In addition, the center-of-mass motion is governed by an harmonic oscillator Hamiltonian
Equation (8), and its energy is given by Equation (3). Then, the energy of a two-body state is the
sum of its relative and center-of-mass energies. Therefore, each relative state has its corresponding
center-of-mass excitations.
These analytical results are used as a test for our numerical calculations and allow us to critically
analyze the dependence of the numerical results on the size of the Fock subspace used in the
diagonalization. As g increases one needs a larger subspace.

3.3.2. Comparison of Analytical and Numerical Results

In Figure 1 we report as a function of g the low-energy part of the two-particle energy spectrum
calculated by diagonalization in a truncated basis space (green dots) and the spectrum obtained by
solving Equation (43). Notice that this equation provides the energies of the relative motion to which
we add the possible energies of the center of mass: ECM . For the ground state, ECM = 1/2. On the
contrary, the diagonalization provides the total energies of the system. The diagonalization has been
performed using 100 single-particle modes that translates when the energy restriction is taken into
account into a dimension of 5050 of the matrix to be diagonalized.
Mathematics 2020, 8, 1196 12 of 38

3
Energy

1
Numerical diagonalization
Analytical
0
-5 0 5 10 15
g
Figure 1. Energy spectrum of the two-particle system as a function of the interaction strength.
The calculations are performed using 100 single-particle modes (green dots). The analytical spectrum
(dashed red line) is obtained with Equation (43) for the energy of the relative system and adding after
the possible energies of the center of mass. Energy given in harmonic oscillator units. The energies
obtained via Equation (43) correspond to even-parity states, the corresponding odd-parity states do not
depend on the interaction, and result in horizontal lines in close agreement with the numerical ones.

We explore both the attractive and the repulsive interaction regimes. In general, for small values
of | g|, both attractive and repulsive, the agreement between both types of calculations is very good.
Apparently, the quality of the results deteriorates faster for attractive interaction. The calculations
discussed in this paper will consider mainly repulsive interactions which will be explored up to the
interaction strength limit which is almost reached for the highest values of g considered, g = 15.
In the plot, we also include the horizontal lines which are energies of states not affected by the
interaction. The first one corresponds to the state with S = 1 and M = 1 whose wave function can
be decomposed as an antisymmetric wave function in coordinate space built with the single-particle
states: n = 0, 1 times the triplet spin state with M = 1. The energy of this state, E = 2, does not depend
on g. In fact, all the states described by horizontal lines in the figure, are eigenstates of the Hamiltonian
with M = 1. Notice also that in the limit g → ∞ states with M = 0 become degenerate with states with
M = 1.
To complete the study of the accuracy of our calculations we investigate the convergence of
the energy of the ground state and first excited state of the two particle system as a function of the
number of harmonic oscillator modes used (NM ). Note that for N = 2, the energy constraint for the
construction of the two-body basis is given by the maximum energy Emax = 1/2 + ( NM − 1/2) = NM .
To this end, in Figure 2 we report the difference between the ground-state energy obtained by
the diagonalization procedure and the analytical ones, for two values of the interaction strength.
The figure also reports the difference of the first excited state. For N = 2, it is possible to establish
the dependence of the dimension of the Hilbert space on NM , which is given by NM ( NM + 1)/2.
As expected, the difference between the calculated and the exact value decreases with the number
1/2
of modes. Following Ref. [53], we fit functions of the type C1 /NM + C2 /NM with excellent results.
1/2
In addition, for large values of NM , the convergence of the energy goes as 1/NM . This dependence
indicates a slow convergence of the numerical results by increasing the number of single-particle
modes. As expected, the differences are larger for the larger strength, but in both cases the differences
Mathematics 2020, 8, 1196 13 of 38

are very small. In all cases the difference is positive, indicating that the numerical results are upper
bounds to the exact ones.

g=1, g.s.
0.2 g=1, ex.
g=5, g.s.
g=5, ex.
0.1
Ediag.-Eanalytical

0.05

0.02

0.01

0.005
20 30 40 50 60 70 100
Number of modes
Figure 2. Differences between the numerical results and the exact ones as a function of the number of
single particle modes, for two values of the interaction strength, g = 1 and g = 5, for the ground (g.s)
1/2
and first excited (ex.) state. The lines are fits of the type C1 /NM + C2 /NM to the calculated points.
Both axes are in a logarithmic scale.

4. Ground State Properties

In this section we discuss and analyze different observables characterizing the ground state
of the system. This analysis is done for several numbers of particles and total spin projections,
devoting special attention to the two-particle system, as it is the only one with analytical solution.
As the interaction is spin independent, the Hamiltonian for a given number of particles is
calculated in a subspace with a well-defined total spin projection M. After diagonalization of the
Hamiltonian, the lowest eigenvalue and its corresponding eigenvector is identified as the ground state
of the system for this spin projection.
The ground state is written in the Fock basis as

|Ψi = ∑ Cn |ψn i , (44)

n

where |ψn i is a vector of the many-body basis of the Fock space, all with the same spin projection and
constructed with the single particle eigenstates of the harmonic oscillator.

4.1. Energy and Virial Theorem

The energy of the ground-state as a function of the interaction strength g for several number
of particles and spin configurations is shown in Figure 3. The first thing to observe is that for each
number of particles the states with maximum spin projection are not affected by the interaction. In fact,
the wave function of these states can be factorized as a symmetric spin function, with all spins up and
an antisymmetric spatial wave function. Due to the Pauli principle, this spatial wave function is built
with the first N single-particle harmonic oscillator states. As discussed in Section 2 its energy is given
by E = N 2 /2. In all other cases, the ground-state energy grows when g increases and tends to saturate
when g → ∞. More precisely, for the N = 2, M = 0, the energy evolves from E = 1 at g = 0, to E = 2
Mathematics 2020, 8, 1196 14 of 38

in the limit g → ∞. Later in this section, we will discuss the structure of both the spatial and the spin
part of the wave function in both limits. For N = 3, M = 1/2 the energy goes from E = 5/2, to the
E = 9/2 for g → ∞. In a similar way, for N = 4, M = 0 the energy goes from E = 4 at g = 0 to E = 8
at g → ∞ while for N = 4, M = 1 the energy goes from E = 5 to E = 8. Notice that in the limit g → ∞
the energy depends only on the number of particles and it is the same for the different spin projections.
The values of the energy at g → ∞ were predicted by Equation (14), E = ( Nd + Nu )2 /2 where Nd and
Nu are the number of particles with spin up and down respectively. On the other hand, in the case
of an attractive interaction, the binding energy increases as the interaction becomes more attractive.
In this case the energies belonging to different number of particles cross in the figure, as the systems
with more particles can accumulate more attraction.
In Figure 3 we also show the results for the cases N = 3, M = 1/2, N = 4, M = 0 and
N = 4, M = 1 reported in Ref. [54] for the repulsive regime. For small interaction strengths, g . 5,
both calculations agree fairly well, but for larger interactions we can appreciate that our calculations
provide larger energy values due to the use of a truncated basis. This effect is more pronounced in the
four particles case, which is performed with less harmonic oscillator single-particle states.

14

12

10

6
Energy

−2 2 particles, M=0 4 particles, M=1

2 particles, M=1 4 particles, M=2
−4 3 particles, M=1/2 5 particles, M=1/2
3 particles, M=3/2 5 particles, M=3/2
4 particles, M=0 5 particles, M=5/2
−6
−5 0 5 10 15
g

Figure 3. Ground-state energy as a function of the interaction strength for different number of
particles and spin configurations. The calculations have been performed by using 100, 50, 40 and
30 single-particle modes to construct the many-body basis for N = 2,3,4 and 5 particles respectively.
The black squares are the corresponding results reported in Ref. [54].

After this general view of the dependence of the ground-state energy on the strength of the
interaction for different systems, we are going to analyse in more detail the two-particle system with
M = 0, which energy as a function of g is plotted in Figure 3. In the non-interacting case, the wave
function can be factorized

Ψ g=0 (1, 2) = Φ0 ( x1 )Φ0 ( x2 )χ(S = 0, M = 0) , (45)

as a symmetric function in space with both particles in the same single-particle state and the two-body
spin function corresponding to a singlet state which is antisymmetric: then the total wave function is
Mathematics 2020, 8, 1196 15 of 38

antisymmetric and its energy is E = 1. In the limit of a very strong interaction, the wave function can
be factorized as [3]

1
Ψ g=∞ (1, 2) = √ |Φ0 ( x1 )Φ1 ( x2 ) − Φ1 ( x1 )Φ0 ( x2 )|χ(S = 0, M = 0) (46)
2

i.e., the absolute value of the determinant built with the ground and first excited single-particle states of
the confining harmonic oscillator times the singlet two-particle spin function. In this case, the presence
of the absolute value ensures that the spatial part of the wave function is symmetric. Notice also
that this wave function does not allow two particles to be at the same position and therefore the
two particles do not feel the interaction. The final energy of this wave function is the sum of the
single-particle energies, E = 1/2 + 3/2 = 2. However, even if the wave function of Equation (46) looks
rather simple, due to symmetry arguments, the diagonalization procedure requires a very large Fock
space to asymptotically approach the solution.
It is also interesting to study the decomposition of the energy in different pieces: kinetic,
harmonic-oscillator potential energy and interaction energy, which are constrained by the
virial relation:
2 h T i − 2 hVho i + hVint i = 0 . (47)

The fulfilment of this relation reinforces the accuracy of the calculations. The derivation of the virial
relation can be found in the Appendix A, and the evaluation of the kinetic energy and the harmonic
trapping potential energy in Appendix B.
In Figure 4 we report, as a function of the interaction strength, the different contributions to
the ground-state energy, for different number of particles and total third-spin component. The total
energy is also reported. In general, for all cases considered, the behavior of the different contributions
is rather similar. For a repulsive interaction, the total energy increases. In the repulsive interaction
range, both the kinetic energy and the oscillator-potential energy are similar. Actually, they are equal
at g = 0 and also tend to the same value in the limit of infinite interaction. The oscillator-potential
energy is greater than the kinetic energy in this range. The interaction energy starts from zero in the
non-interacting case, has a maximum, and goes to zero as the interaction strength tends to infinity.

N=2, M=0 Virial Eint N=3, M=1/2

5 Ekin E
4 Epot
3
2
1
0
8 N=4, M=0 N=4, M=1
6
Energy

4
2
0
-2
12 N=5, M=1/2 N=5, M=3/2
9
6
3
0
-3
0 2 4 6 8 10 0 2 4 6 8 10
g

Figure 4. Different energy contributions to the total ground-state energy, as a function of the interaction
strength for different number of particles and total third-spin component. The total energy and the
fulfilment of the virial theorem is also shown (Equation (47)). The legend is common to all graphics.
Mathematics 2020, 8, 1196 16 of 38

Obviously, the virial theorem is trivially respected for the systems with maximum spin projection,
as in these cases the interaction energy is zero and therefore the kinetic and harmonic potential
energy coincide. In general the virial theorem is well fulfilled, although it deteriorates a little when
the interaction strength is increased. The behavior of the different contributions to the total energy
can be explained by taking into account the virial theorem and knowing the non-interacting and
infinite-interacting limits. In the non-interacting case, the kinetic energy is equal to the oscillator
potential energy. In the infinite interaction limit, the wave function is given by Equation (46),
that prevents two particles from having the same position, therefore, the interaction energy is zero.
As a consequence, the kinetic and the harmonic oscillator energies are also equal and fulfill the virial
relation. For intermediate interactions, the interaction energy is positive, and the kinetic energy is
lower than the oscillator-potential energy in agreement with the virial theorem. The continuity of
the interaction energy between zero and infinite interaction allows us to predict the existence of a
maximum of the interaction energy as a function of g. On the other hand, when the interaction is
attractive, the interaction energy takes negative values, increasing the binding energy as the interaction
becomes more attractive. At the same time, the kinetic energy grows and the harmonic potential
energy decreases, as a consequence of the decrease of the size of the system.

4.2. One Body Density Matrix

The one-body density matrix (OBDM) provides a non-trivial insight on the many-body structure
of the system. The OBDM allows one to obtain the natural orbits, i.e., the eigenvectors of the OBDM,
the momentum distribution and the density profile of the system. In second quantization, the matrix
elements of the OBDM associated to |Ψi are defined as,

ρij = hΨ| a†j ai |Ψi = ∑ Cn∗ Cm hψn | a†j ai |ψm i . (48)

n,m

4.2.1. Density Profile

The first observable to be calculated is the density profile of the ground state of the system,
that provides information on how the particles are spatially distributed in the trap.
The density profile in terms of the matrix elements of the OBDM is given by,

ρ( x ) = ∑ Φi∗ (x)ρij Φ j (x) , (49)

i,j

where Φn ( x ) are the single-particle wave functions used to construct the many-body basis.
This expression corresponds to the diagonal elements of the one-body density matrix in
spatial representation.
In Figure 5 we report the density profile of the ground state of the system for different number
of particles and spin configurations and for several values of the interaction strength. As expected,
the density profiles at g = 0 for the different systems fully agree with the analytical expressions
given in Equation (12). For large values of g, g = 14 in the figure, the density profiles are very close
to the ones predicted in Equation (15). The small differences are due to the finite value of g and
also to the finite size of the Hilbert space where we are diagonalizing the Hamiltonian. In any case,
both for the energy and the spatial distribution of the particles, g = 14 can be considered a strongly
interacting regime. Notice, that for the same number of particles and different spin configurations
the density profiles are rather different, when g is small. However, for g → ∞ the profile for a given
number of particles does not depend on the spin projection. In the strongly interacting limit, the
density profile shows as many peaks as the number of particles, i.e., the particles try to be distributed
equidistantly to minimize the repulsion. The size of the system also increases with the interaction
strength. For attractive interactions, g = −2 is shown in the figure, the density profile is narrower and
reaches higher values than in the repulsive case.
Mathematics 2020, 8, 1196 17 of 38

1.4 N=2, M=0 g=-2 N=3, M=1/2

1.1 g=0
g=2
0.8 g=6
g=14
0.5
0.2
1.6 N=4, M=0 N=4, M=1
1.3
1.0
ρ

0.7
0.4
0.1
1.6 N=5, M=1/2 N=5, M=3/2
1.3
1.0
0.7
0.4
0.1
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
x

Figure 5. Density profiles of the ground state for different number of particles and spin projections.
The density profiles are shown for various values of the interaction strength. The legend is common to
all panels.

In Figure 6 we report the densities of the spin up and spin down particles. The cases N = 2,
M = 0 and N = 4, M = 0 are not reported, because these cases have the same particles with spin
up and spin down, and for symmetry reasons, the densities of both types of particles are equal and
these cases do not provide new information. For the remaining cases, we can see that in every density
profile there are the same number of peaks as particles with the correspondent spin projection. As the
interaction increases, the peaks separate to mimimize the interaction energy.

4.2.2. Natural Orbits

The diagonalization of the density matrix provides its eigenvectors, i.e., natural orbits, and its
eigenvalues, which can be interpreted as the occupation number of the natural orbits. The sum of the
eigenvalues of the natural orbits is normalized to the total number of particles.
In the non-interacting case, the ground state of an N-particle system, which corresponds to a
Slater determinant built with N single-particle harmonic oscillator functions with their spin projection,
produces an OBDM with N-natural orbits with eigenvalue 1 and the all others with eigenvalue
zero. In these cases the natural orbits can be identified with the single-particle wave functions used
in the construction of the Slater determinant. When turning on the interaction, one gets a set of
N-natural orbits with eigenvalues smaller than 1 and additional natural orbits with eigenvalues
significantly smaller. The natural orbits are expressed as linear combinations of the harmonic oscillator
single-particle basis and the Slater determinant built with the N natural orbits with the highest
eigenvalues define the wave function of this type with largest overlap with the ground state. However,
this condition does not imply that the energy corresponding to this wave function is the minimum
energy for a wave function built with a single Slater determinant. The fact that the OBDM has
eigenvalues smaller than 1 points out the existence of correlations beyond the mean-field in the
ground-state of the system which translate into the impossibility to express the ground-state wave
function as a single Slater determinant of single-particle wave functions.
Mathematics 2020, 8, 1196 18 of 38

1.0 1.0
N=3, M=1/2 N=3, M=1/2 g=−2
0.7 0.7 g=0
g=2
g=6
0.4 0.4 g=12
0.1 0.1

1.0 N=4, M=1 1.0 N=4, M=1

0.7 0.7
0.4 0.4
0.1 0.1
ρ↑

ρ↓
1.0 N=5, M=1/2 1.0 N=5, M=1/2

0.7 0.7
0.4 0.4
0.1 0.1
1.0 N=5, M=3/2 1.0 N=5, M=3/2

0.7 0.7
0.4 0.4
0.1 0.1
−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3

x x

Figure 6. Densities of the spin up (ρ↑ , left panels) and spin down (ρ↓ , right panels) particles in
the ground state. The density profiles are shown for various values of the interaction strength
and for different number of particles and spin projection. The legend is common to all panels.
The corresponding total densities are reported in Figure 5.

The eigenvalues associated to the natural orbits can be interpreted as the occupation numbers of
the single-particle states defined by the natural orbits. Notice also that a given natural orbit does not
mix single particle states with different parity or spin projection.
In Figure 7 we report the dependence of the ten largest eigenvalues of the one-body density
matrix as a function of the interaction strength for the ground-state of systems with different number
of particles and spin projections. As explained above, at g = 0 we have N, where N is the number
of particles, eigenvalues equal to unity, and the rest equal to zero. When the interaction increases,
the highest eigenvalues decrease indicating the presence of correlations in the system. Then as the
eigenvalues are normalized to N, the smaller eigenvalues start to increase, indicating the impossibility
to describe the wave function with only one Slater determinant. The fact that the eigenvalues of the
OBDM could be the same, does not necessarily imply that the natural orbits are the same. However,
for N = 2 and M = 0 the eigenvalues appear doubly degenerate independently of the value of g.
One corresponds to spin up and the other to spin down that have identical natural orbits and also the
same spatial distribution. The same is true for the N = 4 and M = 0 case. For g = 0 we have four
natural orbits, two corresponding to the single-particle state n = 0, both with spin up and spin down
and the other to n = 1 also with both spin projections, all with eigenvalue equal to unity. Then when
the interaction increases, the natural orbits mix higher excited harmonic oscillator states but do not
mix the spin projection. The natural orbits are independent of the spin projection and the eigenvalues
separate in two groups each one with degeneracy two. They have not only the same eigenvalue
but also the same spatial structure. Instead, for N = 3 and M = 1/2 the natural orbits of each spin
projection are different.
Mathematics 2020, 8, 1196 19 of 38

N=2, M=0 N=3, M=1/2

1.0

0.7
1 5 9
2 6 10
0.4 3 7
4 8
0.1
N=4, M=0 N=4, M=1
1.0
Eigenvalues

0.7

0.4

0.1
N=5, M=1/2 N=5, M=3/2
1.0

0.7

0.4

0.1
0 2 4 6 8 10 0 2 4 6 8 10
g

Figure 7. The largest 10 eigenvalues of the one-body density matrix (OBDM) of the ground state of
systems with different number of particles and spin configurations as a function of the interaction
strength.

The reason behind the previous observations is that the OBDM is constructed in boxes of well
defined third spin component,

a†↑ a↑ a†↓ a↑
! !
ρ↑↑ ρ↓↑
ρ↑↓ ρ↓↓
= ∑ Cn∗ Cm hψn |
a†↓ a↑ a†↓ a↓
|ψm i . (50)
n,m

However, the crossed terms are zero and only the diagonal terms contribute. Therefore the matrix
separates in two pieces, one with spin up and the other with spin down. If the system is symmetric in
the third spin component, both subspaces are equal and the eigenvalues and spatial eigenfunctions of
each subspace are also equal.

5. Low-Energy Excited States

To complete the study of the structure of the few-body fermion systems, in this section we discuss
the low energy excitation spectrum which will be useful to understand the dynamics and the response
of the system to an external perturbation. In this section we present an algorithmic procedure to assign
quantum labels to the different states, a more formal study of the symmetries in few-atom systems
with examples for three to five atoms is given in Ref. [32].

5.1. Energy Spectrum

Once the Hamiltonian for a given number of particles (N) and total third spin component (M)
has been diagonalized, we have access to the spectrum and the structure of the eigenstates. We pay
particular attention to the dependence on the interaction strength of the low-energy part of the
spectrum, which is shown in Figure 8 for several particles and spin configurations.
The first observation is the existence of states which are not affected by the interaction and
therefore appear as horizontal lines in the plots. These states correspond to antisymmetric wave
functions in space and give zero probability of having two particles in the same position. Actually,
these wave functions can be factorized in an antisymmetric wave function in space and a spin
Mathematics 2020, 8, 1196 20 of 38

symmetric wave function. Later on, in Section 5.2 we discuss how to assign the total spin quantum
number to the different eigenstates.
Another general observation is the existence of excitations with a constant energy shift
independent of the interaction strength which is associated to excitations of the center of mass of the
system. Notice that the interaction which is invariant under space translations does not affect the
center-of-mass motion.
An interesting feature is that the energy of all states saturates when the interaction strength tends
to infinity. Actually, some states merge to the same energy when g → ∞ increasing the degeneracy of
the energy levels in this limit.
Once more, let us analyze more carefully the case N = 2 and M = 0. As explained in the previous
section, the non-interacting ground state corresponds to

Ψ(1, 2) = Φ0 ( x1 )Φ0 ( x2 )χ(S = 0, M = 0) , (51)

with energy E = 1. The two first excited levels have the same energy, E = 2 at g = 0, one corresponds
to a state
1
Ψ(1, 2) = √ (Φ0 ( x1 )Φ1 ( x2 ) − Φ1 ( x1 )Φ0 ( x2 ))χ(S = 1, M = 0) , (52)
2
which is a product of an antisymmetric function in coordinate space and a symmetric one in spin.
This state is not affected by the contact interaction and its energy is independent of g. The energy of
this state merges with the ground state energy when g → ∞. As the Hamiltonian commutes with the
spin, this state is degenerate with states having the same spatial wave function and the spin functions
corresponding to S = 1 but M = 0, ±1.

5
4
3
2
1 N=2, M=0 N=3, M=1/2

9
8
Energy

7
6
5
N=4, M=0 N=4, M=1
4
13
12
11
10
9
N=5, M=1/2 N=5, M=3/2
8
0 2 4 6 8 10 0 2 4 6 8 10
g

Figure 8. Low-energy spectrum for several number of particles and spin configurations as a function
of the interaction strength. Notice that the excitation energies are measured with respect to the
ground-state energy of the non-interacting system.

The other excited state at g = 0 corresponds to

1
Ψ(1, 2) = √ (Φ0 ( x1 )Φ1 ( x2 ) + Φ1 ( x1 )Φ0 ( x2 ))χ(S = 0, M = 0) , (53)
2
Mathematics 2020, 8, 1196 21 of 38

which is symmetric in spatial coordinates and antisymmetric in spin. It can be shown that this state
corresponds to the first excitation of the center of mass times the intrinsic wave function of the ground
state described above. In the ground state, the center of mass is always in the ground state of HCM .
However, in the state we are considering the center of mass occupies the first excited state of HCM .
This is true for all values of g and the energy of this state can be written as E = Eg.s. + 1 for all values
of g. This state has degeneracy one. Actually, the wave function of this state can be written as the
ground-state wave function at a given g times the wave function of the first excited state of HCM
This type of arguments explain the full spectrum for N = 2. It is worth reminding that
our numerical procedure is constructed in the Fock space using single-particle wave functions,
without decomposing the Hamiltonian in the center of mass and intrinsic parts. Therefore the
proper localization of the center of mass excitations provides a test on the numerical accuracy of
the calculations.
We continue our analysis by considering the case for N = 3 and M = 1/2. The ground state
at g = 0 has an energy E = 1/2 + 1/2 + 3/2 = 5/2, which corresponds to a state that populates
the single-particle harmonic-oscillator states: {0 ↑, 0 ↓, 1 ↑}, that results in M = 1/2. Then when g
increases, the state mixes more complicated configurations, the energy increases monotonously and
tends to Eg→∞ = 1/2 + 3/2 + 5/2 = 9/2. The spin of the ground state is S = 1/2, which is the
minimum compatible with the value of M = 1/2 and its parity is negative all along the values of
g. At g → ∞ the ground state merges with the state with S = 3/2, M = 1/2 which is a state built
with a Slater determinant with the single-particle harmonic-oscillator states: n = 0, 1, 2 times the
symmetric wave function of three spins S = 3/2, M = 1/2. The particles in this state do not feel a
contact interaction and the energy of this state is independent of the interaction strength. Notice also
the existence of two states at g = 0 with M = 1/2 built with the single-particle harmonic-oscillator
states {0 ↑, 0 ↓, 2 ↑} and {0 ↑, 1 ↓, 1 ↑} which give origin at g = 0 to two states with total energy
E = 7/2. One describes an excitation of the center of mass built on top of the ground state which
evolves with g keeping always the same energy shift respect to the ground state, and the other state
merges with the ground state when g → ∞.
Next, we discuss the case N = 4 and M = 0. Of course, the number of levels increases with the
number of particles and the number of levels is larger when the total spin projection, M is smaller.
One immediately detects states not affected by the interaction. In particular, the lowest one that
corresponds to a state with S = 2, the maximum spin for N = 4 and M = 0. This state is described by a
wave function that is the product of an antisymmetric wave function given by a Slater determinant built
with the single-particle harmonic-oscillator states: n = 0, 1, 2, 3 and a symmetric spin function of four
particles with S = 2 and M = 0. The total energy of this state, E = 1/2 + 3/2 + 5/2 + 7/2 = 8, is given
by the sum of the first four single-particle energies of the harmonic-oscillator trapping potential.
The ground state is a state with S = 0, i.e., the minimum total spin compatible with the value
of M. The energy of the state at g = 0, E = 4, corresponds to the occupation of the single particle
levels: {0 ↑, 0 ↓, 1 ↑, 1 ↓}. The energy gets more repulsive as g increases and for g → ∞ the energy
merges with the energy E = 8 of the previously discussed state. The parity of the ground-state is
positive. One can also identify the excitations of the center of mass characterized by a constant shift of
the energy respect to the ground-state or on top of a given state. For instance, the first one corresponds
to a center of mass excitation of one unit of energy of the Hamiltonian associated to the center of mass,
HCM . There are other levels whose quantum numbers can be identified and many of them merge
together when g → ∞.
Similar comments apply to the case N = 4, M = 1. One can also identify states not affected by the
interaction. In particular, the lowest one corresponds to S = 2, M = 1, that has the same energy as the
state S = 2, M = 0 of the previous panel. Actually, both states S = 2, M = 1 and S = 2, M = 0 can
be built from the state S = 2, M = 2 constructed with an antisymmetric space wave function times
a symmetric spin wave function with all spins up. These states are obtained by applying the ladder
operator S− to the state S = 2, M = 2. When the Hamiltonian commutes with S− this operation does
Mathematics 2020, 8, 1196 22 of 38

not change the energy of the state. These arguments are discussed at length in Section 5.2. In this
M −box, the ground state starts with E = 5 at g = 0 which corresponds to the Slater determinant built
with the single-particle states: {0 ↑, 0 ↓, 1 ↑, 2 ↑} and when g → ∞ merges with several states that tend
to E = 1/2 + 3/2 + 5/2 + 7/2 = 8.
Finally let us discuss very briefly the case N = 5, S = 1/2. In this case, there are many levels
below the first level which is not sensitive to the interaction strength, with an energy E = 1/2 + 3/2 +
5/2 + 7/2 + 9/2 = 25/2 and S = 5/2, M = 1/2. The ground state, has S = 1/2, and positive parity.
At g = 0 it has an energy E = 13/2 and when g → ∞ it merges with the state that is non-sensitive to
the interaction discussed above, with an energy E = 25/2.
For the case with M = 3/2 there are considerably fewer levels. The ground state at g = 0 has
more energy than the ground state for M = 1/2, basically due to the Pauli principle. In fact, at g = 0
the ground-state is built with the single particle states {0 ↑, 0 ↓, 1 ↑, 2 ↑, 3 ↑}, which gives an energy
E = 17/2. However, when g is increased the state merges with the state S = 5/2, M = 3/2 which is
also degenerate with the sate S = 5/2, M = 1/2 discussed in the previous panel.
Whenever possible, we have compared our results with the ones provided in Ref. [33] obtaining a
good agreement in all cases considered. In fact our results, which provide upper bounds, are in general
slightly below the results of Ref. [33].
Finally, let us point out that in the limit of large interaction several groups of states merge to the
same energy giving rise to an increase of the energy degeneracy. In particular, for the ground state this
degeneracy is expressed as D = NuN! !N ! .
d

5.2. Spin Determination

As the Hamiltonian does not depend on the spin, it commutes with the spin operators and should
be possible to assign a total spin to the eigenstates of the Hamiltonian. In particular, it commutes
with the ladder lowering S− and raising S+ spin operators. Therefore, if we have a state with a total
spin S and spin projection M with energy E, due to the commutation relations of S+(−) with the
Hamiltonian, it turns out that by applying the operators S+(−) to this state one obtains a state with the
total spin projection increased (reduced) by 1, i.e., M ± 1, with the same energy E. Therefore, for each
eigenstate with a well-defined spin, there will be 2S + 1 states with the same energy. We can use this
fact to determine the total spin of the states. The argument is as follows: Given N fermions, we start
considering the maximum possible spin projection, that will be the state with all spins pointing up,
with M = N/2. In this case, the wave function factorizes in a spatial antisymmetric wave function
times a symmetric spin wave function. This wave function has a well defined spin which coincides
with the maximum value of the spin projection, S = N/2. Then if we apply 2S times to this state the
operator S− , we will get the remaining states belonging to this spin multiplet, all having the same
energy and all of them factorized as an antisymmetric wave function in space times a spin symmetric
wave function. Obviously, these states are found in different M-box subspaces. In fact, if we consider
the Hamiltonian in the box with M = S − 1, among the energies obtained in this box, we find the
previous energy attributed to S = N/2, corresponding in this case to the state with M = N/2 − 1.
Besides, if the remaining states in this M-box are not degenerate, one can assume that the total
spin of each one will be the minimum S compatible with the value of M, i.e., S = N/2 − 1. One can
have more than one state with this spin and one would need another quantum number to distinguish
between them. Next, we take the box of the Hamiltonian associated to the value of M, M = N/2 − 2.
Here we immediately identify the energies obtained for the previous spin values and again the rest of
states have a spin that is the minimum compatible with the value of M. One can continue this process
until the total spin projection has the positive lowest half-integer value for odd N or zero value in case
of even N.
To illustrate this procedure, we consider the lowest-energy states of the case N = 4 shown in
Figure 9 as a function of the interaction strength. The maximum possible spin projection is N/2 = 2.
The wave functions with S = 2 and M = 2 can be factorized in an antisymmetric spatial wave function
Mathematics 2020, 8, 1196 23 of 38

times a spin wave function with all spins pointing up. In particular, the lowest energy for this states
(S = 2 and M = 2), E = 8, is the sum of the first four single-particle harmonic oscillator energies,
E = 1/2 + 3/2 + 5/2 + 7/2. In addition, this state is not affected by the interaction and therefore
appears as the lowest horizontal line in figure. This energy level corresponds also to the energy of the
different possible 2S + 1 states of the S = 2 multiplet. In the next step we consider the box with M = 1.
First we identify the previous energy and find out several non-degenerate eigenstates which can all be
associated to S = 1 (dot-dashed red lines in the figure). Notice that all of them have different energies,
but S = 1. For instance, at g = 4 we find three S = 1 levels with energies below E = 8. In the next step
we diagonalize the subspace with M = 0 and identify the M = 0 states corresponding to S = 2 and to
the S = 1 states detected before. The rest of levels have S = 0. At g = 5, we find three S = 0 levels
below E = 8, all with different energy, and non degenerate. Notice that one of them corresponds to a
center of mass excitation and that merges with other levels at E = 9 when g → ∞. Using these types of
arguments one can assign the spin to the eigenstates of the Hamiltonian even if the total wave function
could be complicated and most of the times non-factorizable in a spin and a space part.

10

8
Energy

5 S=0
S=1
S=2
4
0 2 4 6 8 10 12 14
g
Figure 9. Low-energy spectrum of the fermionic system with four particles as a function of the
interaction strength. The states are labelled according to their total spin, as explained in the text.

6. Dynamical Excitation
In Sections 4 and 5 we have concentrated on the static properties of the fermionic system.
In particular, we have analyzed the ground state and also the main features of the low-energy spectrum.
In the present section we turn our attention to the dynamics of the fermionic system in the harmonic
trap. Simulating the dynamics of quantum many-body systems is important for current quantum
technological applications where one could for instance design quantum systems to produce desired
many-body states after a certain evolution [55–57]. A second relevant aspect is that the dynamics of the
system reflects in many ways its internal structure. In this sense, one can devise a dynamical evolution
in order to unveil the spectral structure of the system [58]. This is the main goal of this section.
To this aim we consider two different perturbations to the ground state of the system, which are
sensitive to the low-energy spectra presented before. The first is a sudden change in the trap frequency.
This perturbation excites both center-of-mass and relative modes of the system and can be analyzed by
computing the dynamic structure function. In the second one, we numerically obtain the time evolution
of the system after a sudden quench of the interaction strength. In this case, the center-of-mass modes
Mathematics 2020, 8, 1196 24 of 38

are not excited and we look for traces of the internal structure on the time evolution of the central
density of the cloud.

6.1. Sudden Change in the Trap Frequency: Breathing Mode

A well-known way to study the internal structure of a quantum many-body system trapped in
a harmonic oscillator is by exciting the so called breathing mode. In contrast with the Kohn (dipole)
mode, where the cloud is initially displaced from the minimum of the harmonic potential [59], in our
case we study the response of the system to a change in the trapping frequency. For the two particle
case there are studies such as [60,61], and using different confining potential such as [62]. The main
tool we consider is the dynamic structure function associated to the mono-polar excitation operator

N
F̂ = ∑ x2j . (54)
j

This perturbation preserves spin and parity, therefore, it only connects the ground state of the system
with other excited states with the same total spin and parity.
This operator can be separated in two pieces: a center-of-mass and an intrinsic one, i.e.,
2 + 1 2
F̂ = NXCM N ∑ j<i ( xi − x j ) . The center of mass of the system is described by an harmonic oscillator
Hamiltonian, and the wave functions associated to this part are the harmonic oscillator wave functions.
Note that this perturbation can excite the center of mass or the intrinsic part, but not both at the
same time.

6.1.1. Dynamic Structure Function

The dynamic structure function, normalized to the number of particles, is defined as

1 2
SF ( E) =
N ∑ hΨi | F̂ |Ψ0 i δ ( E − ( Ei − E0 )) , (55)
i >0

where F̂ is the excitation operator associated to the external perturbation, see Equation (54), |Ψ0 i is the
ground state of the unperturbed system and |Ψi i are its excited states.
In the second quantization formalism, the structure function reads,

2
1
SF ( E) =
N ∑ ∑ Cn,0 Cm,i
∗
∑ hk| x2 |l i hψm | a†k al |ψn i δ (E − (Ei − E0 )) , (56)
i >0 n,m k,l

where |ki and |l i are harmonic-oscillator single-particle states and |ψn i are the Fock states of our
2 |j
basis. For the center of mass, the hiCM | XCM CM i are different from zero only when iCM = jCM or
iCM = jCM ± 2. This implies that the center of mass can be excited at most by two energy quanta.
In some limiting cases we know the analytic value of the dynamic structure function. These cases
can be used to benchmark our numerical calculations [34]. For the non-interacting case, we expect a
single peak in the dynamic structure function with an energy E − E0 = 2. The intensity of this peak
depends on the number of particles and the total spin projection. Likewise, in the infinite interaction
limit, we also expect a single peak in the dynamical structure function with energy E − E0 = 2.
In addition, a peak with energy E − E0 = 2 is expected for any value of the interaction strength, due to
the center-of-mass excitation, with a constant intensity as a function of the interaction strength. In the
range of repulsive interaction, we expect more peaks associated to states with the same spin and parity
as the ground state. These correspond to intrinsic excitations.
In Figure 10 we present the dynamic structure function for the two fermions, for the
non-interacting case and for values of the interaction strength g = 1, g = 2, and g = 15. We are
considering the total spin projection S = 0, i.e., the ground state and all the states excited by the
Mathematics 2020, 8, 1196 25 of 38

monopolar operator to have total S = 0 and even parity. As expected, for the non-interacting case
we only have one peak at excitation energy E − E0 = 2. For any interaction strength there is a peak
at E − E0 = 2 with the excitation energy associated to the center-of-mass excitation with an intensity
independent of the interaction strength. The rest of the peaks observed are associated to intrinsic
excitations. In addition, for large interaction strengths, we recover only one peak at E − E0 = 2, and the
rest of the peaks vanish. The intensity of the excitations decreases for the states with more energy,
and there is a dominant peak near the energy of the center-of-mass excitation. In the non-interacting
case, the single peak has two contributions: the center of mass, and the dominant intrinsic peak.
When we turn on the interaction, the dominant intrinsic peak moves to a lower energy, but for large
interaction, returns to the same energy than the center of mass excitation. This re-entrant behaviour,
reported for bosons in Ref. [63], is also reflected in the energy spectrum for the two particles case.
For the remaining excitations, we can identify the correspondent states via the energy spectrum: all
of them are excitations of the relative motion, with the same spin and parity than the ground state.
Notice that the good location of the center-of-mass excitation is a good test on the numerical evaluation
of S F ( E), which in Equation (56) is expressed in terms of single-particle matrix elements and the Fock
basis. It is also worth mentioning that in this case S F ( E) for two particles with M = 0 is the same
as if the particles were two bosons Ref. [34]. In the case of N = 2 with M = 0, it has been possible
to identify the center-of-mass and intrinsic excitations, because we know the analytic spectrum [44].
This is not the case for the rest of configurations. The N = 4 with M = 0 configuration is a good
example to explore the behaviour of the dynamic structure function of the breathing mode in a more
complex case.

100 g=0 g=1

10-1
10-2
10-3
10-4
10-5
10-6
SF

100 g=2 g=15

10-1
10-2
10-3
10-4
10-5
10-6

1 2 4 6 10 1 2 4 6 10
E-E0

Figure 10. Dynamic structure function of a mono-polar excitation for the case of two fermions.
The total spin of the ground state and of all excited states connected to it through the mono-polar
excitation, are zero. The different panels correspond to four different values of interaction
strength. From the non-interacting case, g = 0 to a strongly interacting one g = 15.

In Figure 11 we report the dynamic structure function associated to the breathing mode for N = 4
with M = 0, for different values of the interaction strength, concretely, for the non-interacting case
(g = 0), for g = 1, g = 2 and g = 15. As in the case of two particles, in the non-interacting case, there is
only one peak with energy E − E0 = 2. In the interacting case, there is a peak at energy E − E0 = 2
Mathematics 2020, 8, 1196 26 of 38

with a constant intensity, which we identify with a center-of-mass excitation. All the peaks present
can be associated to excitations of the ground state to states with the same spin and parity as the
ground state. Using this information, we verify that the spin determinations in Section 5.2 are in
agreement with the results obtained here. Note that not all the states that verify these conditions (have
the same spin and parity as the ground state) are excited. This is because the state has to be an intrinsic
excitation or a center-of-mass excitation. In addition, in this case we also observe a dominant peak
with a re-entrant behaviour. This is common to all the cases with a different number of particles and
total spin projection we have considered. Finally, we can see that in the limit of strong interaction, the
intensity of all the peaks, except the dominant and the center of mass, decrease, tending to zero for the
infinite interaction limit.

100 g=0 g=1

-1
10
10-2
10-3
10-4
10-5
10-6
SF

100 g=2 g=15

10-1
-2
10
-3
10
10-4
10-5
10-6

1 2 4 6 10 1 2 4 6 10
E-E0

Figure 11. Dynamic structure function of a mono-polar excitation for the case of four particles. The total
spin of the ground state and of all excited states connected to it through the mono-polar excitation,
are zero. The different panels correspond to four different values of interaction strength. From the
non-interacting case, g = 0 to a strongly interacting one g = 15.

6.1.2. Sum Rules

The energy momenta of the dynamic structure function are a useful tool to analyze the response
of the system to an excitation. They are defined as
Z
Mn = En S F ( E)dE , (57)

where E is the excitation energy. They can be explicitly computed from the definition of S F ( E),
Equation (56), by using the second quantization formalism as:

2
1
Mn =
N ∑ ∑ Cn,0 Cm,i
∗
∑ hk| x2 |l i hψm | a†k al |ψn i (Ei − E0 )n . (58)
i >0 n,m k,l

If the dynamic structure function is known, the energy momenta can be computed directly using
Equation (57). In general however, the dynamic structure function is difficult to calculate and is not
exactly known. In these situations there are useful theorems that allow us to compute Mn using only
ground state properties [64]. These relations are named sum rules.
Mathematics 2020, 8, 1196 27 of 38

For the monopolar excitation operator, the lowest sum rules can be simply expressed as the
following expectation values on the ground state [34]

1 1 ∂2 E0 (λ)
M−1 = −
2 N ∂λ2 λ =0
4 (59)
M1 = hV i
N ho
4
M3 = (h T i + 3 hVho i) .
N
Using perturbation theory, we have an alternative expression for M−1 ,

1 1 ∂
M−1 = − h0̄| F̂ |0̄i , (60)
2 N ∂λ λ =0

where |0̄i and E0 (λ) are the ground state and its corresponding energy of the perturbed Hamiltonian
0
Ĥ = Ĥ + λ F̂.
In Figure 12 we report the energy momenta M−1 , M1 and M3 as a function of the interaction
strength. The sum rule M−1 has been calculated only using the explicit values of the dynamic structure
function. For M1 and M3 , the two methods, direct integration of the dynamic structure function and
sum rules, are in perfect agreement for the case of two particles. For N = 4 the agreement between
both methods is also very good. In all cases, the three momenta considered grow monotonically as the
interaction strength is increased.

0.5 1.0

0.4 0.8
M-1

0.3 0.6
N=2, M=0 N=4, M=0
0.2 0.4
4.0
2.0
1.8 3.5
1.6 3.0
M1

1.4 2.5
1.2 N=2, M=0 N=4, M=0
2.0
1.0
16.0
7.0
14.0
M3

6.0 12.0
5.0 10.0
N=2, M=0 N=4, M=0
4.0 8.0
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
g g

Figure 12. Values of the three energy momenta M−1 , M1 and M3 for the cases of two and four particles
as a function of the interaction strength. The calculations done using the explicit value of the dynamic
structure function, Equation (57) are represented by dots and the line is the value computed using the
sum rules, Equation (59). For M−1 the results from Equations (59) and (60) are indistinguishable on
the figure.

The sum rules can also be used to obtain approximate values of the excitation energies. This is a
tool which is for instance of interest for quantum Monte Carlo methods, where the low-energy spectrum
of the system cannot be computed, see for instance the bosonic case in [63]. Indeed, if only one peak
Mathematics 2020, 8, 1196 28 of 38

p
appears in the dynamic structure function, we can compute its excitation energy by Mn /Mn−2 [64].
In a more general case this ratio can at least be used as an estimation of the energy of the main peak.
p √
In Figure 13 we report the ratios M1 /M−1 , M3 /M1 and excitation energy of the main
contribution, which is associated to an intrinsic excitation, as a function of the interaction strength.
p √
We can observe that in the case of two particles, both estimations, M1 /M−1 and M3 /M1 , differ and
neither have the same energy of the intrinsic one, indicating that there is more than one peak in the
dynamic structure function. At g = 0 , there is only one peak in which two states coexist, i.e.,
an intrinsic and a center-of-mass excitation, both estimates coincide and provide the exact excitation
energy. Then the differences grow, indicating the existence of more than one peak in S F ( E). In fact,
there are higher intrinsic excitations whose strength decrease quickly besides the center-of-mass
excitation which is fixed in the excitation energy 2. The behavior of this difference is a consequence
of the re-entrance phenomena mentioned above. For large g the differences decrease again and for
g → ∞ both estimations coincide and provide the exact excitation energy which again is the same for
the center of mass and the only excited intrinsic state, similarly to the non-interacting case Figure 10.
When the number of particles increases, the strength of the center-of-mass peak decreases and the two
ratios provide closer estimates which in turn are close to the main intrinsic excitation energy. This can
be seen in Figure 11, where the intrinsic excitation is clearly dominant over the rest of excitations.
Similar to the two particle case, the interaction strength region where the estimates are worse is where
the intrinsic excitation energy takes its lower value.

N=2, M=0 N=3, M=1/2

2.0

1.9 Intrinsic
√(M1/M-1)
√(M3/M1)
1.8
N=4, M=0 N=4, M=1

2.0
Eexc

1.9

1.8
N=5, M=1/2 N=5, M=3/2

2.0

1.9

1.8
0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
g
p √
Figure 13. Ratios of the sum rules M1 /M−1 and M3 /M1 as a function of the interaction strength
for several number of particles and spin configurations. The purple + signs correspond to the excitation
energy of the main peak of the dynamic structure function.

6.2. Interaction Quench

An alternative way to explore the internal excitations of an interacting many-body system is by
performing an instantaneous quench of the interaction strength, which for the two particle case has
been studied in [65,66], and for different confining potential in [67]. In current experiments in ultracold
atomic gases this is feasible. Indeed, a great control on the interaction strength is achieved by means of
Feshbach resonances [14].
A possible protocol is to prepare the system in its ground state for a certain value of the interaction
strength, then the latter is suddenly changed and as a consequence the original state is non-stationary
Mathematics 2020, 8, 1196 29 of 38

anymore and evolves in time, populating the new set of eigenstates corresponding to the new value of
the interaction strength.
This way to probe the system has similarities with the breathing mode. It also preserves
spin and parity. However, due to the invariance of the interaction under translations, there are
no center-of-mass excitations.

6.2.1. Time Evolution of the Perturbed System

We want to study the time evolution of the ground state corresponding to a given value of g after
performing a sudden quench of the interaction strength to a new value gnew . The initial ground state
is no longer the ground state of the new Hamiltonian. However we can express the old ground state in
terms of the eigenvectors of the new Hamiltonian and calculate its time evolution:
0
|Ψ(t = 0)i = ∑ cn |Ψn i , (61)
n

0
where |Ψ(t = 0)i is the ground state of the Hamiltonian before the quenching, and the states |Ψn i are
the eigenfunctions of the Hamiltonian after the quenching
0 0 0 0
H |Ψn i = En |Ψn i , (62)
0
with eigenenergies En . The coefficients cn can be computed as:
0
cn = hΨn |Ψ(t = 0)i . (63)

The time evolution of the state after the quenching is calculated by taking into account that the states
0
|Ψn i are the eigenstates of the new Hamiltonian,
0 0
|Ψ(t)i = ∑ cn e−iEn t |Ψn i . (64)
n

0
In turn, the eigenstates |ψn i can be expressed in terms of the many-body basis used to describe
our system:
0
|Ψn i = ∑ Cn,i |ψi i , (65)
i

where |ψi i is a state of the many-body basis.

Knowing the time evolution of the state we can study the time dependence of any observable.
In particular we are going to consider the time evolution of the central density of our system:
0 0
ρ( x = 0, t) = ∑ ∑ ∑ c∗n cm e−i(Em −En )t Cn,k
∗
Cm,l hψk | a†j ai |ψl i Φi∗ (0)Φ j (0) , (66)
i,j k,l n,m

where Φi ( x ) are the single-particle harmonic oscillator wave functions. Notice that the time evolution
0 0
of the system is governed by Em − En , involving all energy differences of the states of the new
Hamiltonian that have non-zero overlap with the ground-state of the old Hamiltonian. In addition,
it is important to realize that the relevance of a given frequency 2πνn,m = En − Em depends on the
product of the projection coefficients c∗n cm .
For instance, it is plausible to think that, for a small change of g, the original ground state should
have a large overlap with the new ground state. In this situation the dominant frequencies will be
0 0
associated to the differences En − E0 . In this sense, the frequencies would be similar to the frequencies
0
involved in the calculation of the dynamical structure function associated to the breathing mode of H ,
with the exception that the excitations of the center of mass are not present.
Mathematics 2020, 8, 1196 30 of 38

In Figure 14, we report, for the two-particle system, the projections of the ground state of the
0
non-interacting system to the eigenstates of H ( g = 1). We also report the projection for the quench
g = 0 → 5. In the case g = 1 (small g) there is a dominant projection into the ground state of
0
H and then the value of the projections decreases rather quickly. For the case g = 5, the largest
projection is smaller than in the previous case but the projections to higher states decrease more slowly.
Therefore for g = 5 one also has relevant projections into higher excited states.

1 g=1
g=5
|cn|=|<Ψ‘n|Ψ(t=0)>|

0.1

0.01

0 5 10 15 20 25 30 35
‘
En

Figure 14. Projection of the initial ground state for g = 0 into the eigenstates of g = 1 and g = 5 as a
function of its energy, for a two particles systems with zero total spin. The y axis is in logarithmic scale.

6.2.2. Central Density Oscillations

In this section we propose a method to determine the low energy spectrum of a trapped few-body
system, by applying a sudden quench of the strength of the interatomic interactions. In principle,
the procedure can be implemented experimentally.
One can prepare a system with N particles trapped in a harmonic potential without interactions.
Then apply to the system a sudden quench of the interaction, let the system evolve and measure the
central density as a function of time. One should observe an oscillatory behavior.
In order to obtain information about the energy spectrum of the system at g = 1 (the final
value of the quench), we perform a Fourier analysis of the central density as a function of time,
see Figure 15. In this way one obtains the characteristic oscillation frequencies which in turn are
related to differences of the energies of the system for g = 1. Therefore, from the frequencies of the
Fourier analysis one can reconstruct the excitation energies ∆E = 2πν, mainly with respect to the
ground state. To extract the frequencies by Fourier analysis one needs to know the time evolution of
the central density in an interval of time much longer than the one shown in the inset of Figure 15.
Notice also that it is convenient to subtract the average value of the central density in order to avoid a
peak with zero frequency in the Fourier analysis. Actually, this peak at ν = 0 does not provide any
physical information about the excitation energy spectrum and could overlap with other peaks at
low frequencies. We observe a dominant peak at low energy, and more peaks with smaller intensity
at higher energies. Although the time interval used to perform the Fourier analysis is very large,
the Fourier frequencies still come out with an uncertainty width.
Mathematics 2020, 8, 1196 31 of 38

0.09
N=2, M=0 1.4
0.08 N=4, M=0 N=2, M=0
1.3
N=4, M=0

ρ(x=0,t)
0.07 1.2
1.1
0.06
Fourier transform

1
0.9
0.05
0 2 4 6 8 10 12
0.04 t
0.03

0.02

0.01

0
0 1 2 3 4 5
Excitation energy

Figure 15. Frequency analysis of the time dependence of the central density after an interaction quench
from g = 0 to g = 1. The average density has been subtracted. We consider the cases N = 2 (solid lines)
and N = 4 (dot-dashed lines). In the inset we depict the first oscillations of the signal used to compute
the Fourier analysis shown in the main panels. For the Fourier analysis we used a time interval, t = 500.
All magnitudes in the figure are in harmonic oscillator units.

The location of the peaks allows us to recover the excitation energies of the low-excited states
with respect to the ground state corresponding to g = 1. However, to find the absolute energies of
these states would be necessary to have access also to the ground-state energy of the system.
The figure reveals the existence of a dominant peak at an excitation energy a little smaller than
2, and other two main visible peaks at energies a little smaller than 4 and 6, with much smaller
strength. All these peaks correspond to energy differences of the first excited states, which have
non-zero overlap with the non-interacting ground state (g = 0), and the ground state of the system
corresponding to g = 1, respectively. Notice that the re-entrance phenomenon, discussed in the
previous section, is reflected in the fact that the excitation energies are a little smaller than 2, 4, and 6
respectively. It is also worth mentioning the existence of a smaller but distinguishable peak, very close
to the excitation energy 2, that corresponds to the difference between the second and first intrinsic
excitations (see Figure 8). No center-of-mass excitations affect this analysis.
In the same figure, Figure 15, we present also the case of four particles and spin projection M = 0.
The shape of ρ( x = 0, t) shows also an oscillatory behavior, which should be determined by several
frequencies. There are two dominant peaks which can be identified, by looking at the excitation
spectrum of this system for g = 1 in Figure 9, with the lowest excitation energies with respect to the
ground state. Actually, the small peak that appears at very low energy, corresponds to the energy
difference of these two dominant peaks. Notice also that as the ground state has S = 0, and the
interaction is spin independent, the quench cannot connect states with different spins and therefore
they do not participate in the time dependence of ρ( x = 0, t). As explained for the two-particle case,
due to the translation invariance of the interaction, the center-of-mass excitations do not play any role
in the analysis. A similar procedure for the time dependence of the mean square radius of the system
leads to a determination of the excitation energies consistent with the ones obtained by the analyzing
the time dependence of the maximum of the density.
Mathematics 2020, 8, 1196 32 of 38

7. Summary and Conclusions

In this paper we have presented a quantum microscopic description of few spin 1/2 fermions
trapped in a 1D harmonic oscillator potential. The fermions interact through a contact interaction
that, due to the Pauli principle, is active only between particles with opposite spin. Along the paper,
we have mainly concentrated on repulsive interactions. These types of systems have already been
realized experimentally and new experiments concerning the structure and mainly the dynamics are
expected for the near future. Our main objective has been to study both static and dynamic properties
of the system as a function of the interaction strength. This interaction, which provides a good
modelization of real interactions, can be experimentally controlled by means of Feshbach resonances.
After a brief introduction, in the second section, we have described the Hamiltonian of the system,
that includes both the harmonic oscillator confining potential and the interaction between the fermions.
Furthermore, two analytical limits have been discussed in detail: the non-interacting and the infinite
interacting limits. In both cases, the many-body wave function, the energy and the density profile
of the ground state have been derived. In addition, the second quantization formalism has been
introduced as the framework for the calculation of the Hamiltonian.
In the third section we have described in detail the numerical procedures used in the paper,
mainly the exact diagonalization techniques. Special attention has been devoted to the convergence
problems associated with the dimension of the subspaces of the Hilbert space used to diagonalize the
Hamiltonian. The benchmark of the analytical results for the two-particle system has been used to test
the accuracy of the numerical methods.
In the next section, we have studied several static properties of the ground state of the system,
mainly the total energy and their contributions. The energy contributions have been used to check
the fulfillment of the virial relation, obtaining very good results for repulsive interaction, specially in
the low-interacting regime. The fulfilment of the virial relation reinforced the consistency of the
calculations. We also have discussed the density profile of the ground state. In the non-interacting
and the infinite interacting cases our calculations are in perfect agreement with the analytical results
presented in the first section. In the next step, we have computed the natural orbits and its eigenvalues.
Their analysis reveals the presence of important correlations beyond mean field in the system when
the interaction is turned on.
In the last section, we have studied several dynamical properties. First we have studied the
response of the system to a mono-polar excitation over the ground state. We have computed the
dynamic structure function associated to this perturbation, and using the knowledge of the energy
spectra obtained in the previous section, we have identified the excited states. The spin assignment
was useful to identify the excited states by the breathing mode. We have also calculated the energy
p
momenta (Mn ) associated to the dynamic structure function and used the ratios Mn /Mn−2 to
p
estimate the average excitation energies. The differences between the ratios Mn /Mn−2 for different
n indicates the presence of more than one excited state in the dynamic structure function. In fact,
for the two particle case there is more than one dominant peak in the response of the system. However,
for larger number of particles, even if there are more excited states, the strength of the response is
mainly concentrated in one dominant peak.
Finally, we have studied the effect of a quench of the interaction strength on the system. We have
discussed the similarities and the differences between the quench and the breathing mode. In addition,
we have proposed a possible experimental procedure to measure the excitation energies based on the
measurement of the time evolution of the central density of the system after a quench of the interaction.
Analyzing the frequencies of the Fourier transform of the time dependence of the oscillation of the
central density after the quench one can obtain information about the excitation energies. We expect
that the methods and results presented in the paper stimulate the use of exact diagonalization methods
to study the dynamics and time evolution of systems with a few fermions in more complicated
geometries and different interactions.
Mathematics 2020, 8, 1196 33 of 38

Author Contributions: Conceptualization, A.R.-F., A.P. and B.J.-D.; software, A.R.-F.; writing, A.R.-F., A.P. and
B.J.-D. All authors have read and agreed to the published version of the manuscript.
Funding: This work is partially funded by MINECO (Spain) Grants No. FIS2017-87534-P and from Secretaria
d’Universitats i Recerca del Departament d’Empresa i Coneixement de la Generalitat de Catalunya, co-funded by
the European Union Regional Development Fund within the ERDF Operational Program of Catalunya (project
QuantumCat, ref. 001-P-001644).
Acknowledgments: We thank very useful correspondence from M. Simeon, T. Sowinski and N. Zinner.
Conflicts of Interest: The authors declare no conflict of interest.

Appendix A. Derivation of the Virial Theorem

The energy of the system is given by the expectation values h Ei = h T i + hVho i + hVint i, where the
operators are defined as

1 ∂2 xi2
T= ∑− 2 ∂xi2
, Vho = ∑ 2
, Vint = ∑ gδ(xi − x j ) . (A1)
i i j <i

We use the virial theorem to obtain a relation between the different energy contributions. The virial
theorem is based on a scaling transformation of the many-body wave function and the transformation
of the different energy contributions under the scaling transformations. We start by defining a
wave function
Ψλ ( x1 , x2 , ..., x N ) = λ N/2 Ψ(λx1 , λx2 , ..., λx N ) . (A2)

This wave function has the same normalization than Ψ.

First, we compute the transformation of the kinetic energy,

1 ∂2
h T (λ)i = hΨλ | ∑ − |Ψ i = λ2 h T (λ = 1)i . (A3)
i
2 ∂xi2 λ

The harmonic oscillator energy under this transformation is

xi2
hVho (λ)i = hΨλ | ∑ |Ψλ i = λ−2 hVho (λ = 1)i . (A4)
i
2

Finally, the interaction energy transforms as

hVint (λ)i = hΨλ | ∑ gδ( xi − x j ) |Ψλ i

j <i
N ( N − 1) N
Z
= gλ dx1 ...dx N Ψ(λx1 , λx2 , ..., λx N )δ ( x1 − x2 ) Ψ(λx1 , λx2 , ..., λx N )
2 (A5)
N ( N − 1) N 1
Z
= gλ dy1 ...dy N Ψ(y1 , y2 , ..., y N )λδ (y1 − y2 ) Ψ(y1 , y2 , ..., y N )
2 λN
= λ hVint (λ = 1)i ,

where we have implemented the change of variable λxi = yi for all i. The expected value of the
Hamiltonian with the scaled wave function can be written as

h Eλ i = hΨλ | H |Ψλ i = λ2 h T (λ = 1)i + λ−2 hVho i + λ hVint i , (A6)

and taking into account that at λ = 1 the energy has a stationary point, i.e.,

∂ h Eλ i
= 0, (A7)
∂λ λ =1
Mathematics 2020, 8, 1196 34 of 38

we derive the virial relation

∂ h Eλ i
= 0 = 2λ h T (λ = 1)i + −2λ−3 hVho (λ = 1i + hVint (λ = 1)i
∂λ λ =1 λ =1 (A8)
=2 h T i − 2 hVho i + hVint i = 0 .

Appendix B. Evaluation of the One-Body Matrix Elements

The Hamiltonian in second quantization has been written in the harmonic oscillator basis and
requires the knowledge of the single-particle matrix elements of both the harmonic trapping potential
and the kinetic energy. Just for completeness, we include here a brief derivation.
The states of the harmonic oscillator single-particle basis read

1 −x 2 /2 2 /2
|ni = p √ Hn ( x )e |χmn i = An Hn ( x )e− x | χmn i , (A9)
2n n! π

where An is a normalization constant, which depends on n. On the other hand |χm i is a single-particle
spin-1/2 wave function, with spin projection m.
In order to operate with the Hermite polynomials, we use the following properties:

(a) The m-th derivative of a Hermite polynomial,

(m) m!
Hn = 2m Hn−m ( x ) , (A10)
(n − m)!

(b) The recurrence relation,

Hn+1 ( x ) = 2xHn ( x ) − 2nHn−1 ( x ) , (A11)

(c) The orthogonality of the Hermite polynomials

Z ∞ √
2
Hm ( x ) Hn ( x )e− x = δm,n 2n n! π , (A12)
−∞

(d) The n-th power of x expressed in terms of Hermite polynomials

n/2
n! 1
xn =
2n ∑ m! ( n − 2m)!
Hn−2m , (A13)
m =0

(e) The product of two Hermite polynomials as a function of the sum of Hermite polynomials
n
m! Hm−n+2r ( x )
Hm ( x ) Hn ( x ) = 2n n! ∑ ( n − r ) ! ( m − n + r ) ! 2r r!
,n ≤ m. (A14)
r =0

First, we derive the expression of the one-body harmonic oscillator matrix elements, that reads

x2
h j| Vho |i i = h j| |i i . (A15)
2

Using Equation (A13) we can write x2 as

 
2 1 1
x = H2 + H0 . (A16)
2 2
Mathematics 2020, 8, 1196 35 of 38

Then, writing explicitly the Equation (A15), and taking into account Equation (A14) and the
orthogonality Equation (A12), we get
 
1 1 2
h j| Vho |i i = h j| H2 + H0 Ai Hi e− x /2 |χmi i
4 2
 
A 1 2
= h j| i 2n(i − 1) Hi−2 + (2i + 1) Hi + Hi+2 e− x /2 |χmi i (A17)
4 2
Ai A j Z
 
1 2
D E
= dxHj ( x ) 2n(i − 1) Hi−2 ( x ) + (2i + 1) Hi ( x ) + Hi+2 ( x ) e− x χmi χm j .
4 2

Finally, the harmonic oscillator matrix element is given by

q 
1
q
h j| Vho |i i = i (i − 1)δj,i−2 + (2i + 1)δj,i + (i + 2)(i + 1)δj,i+2 δmi ,m j . (A18)
4

Next, we derive the expression of the one-body kinetic matrix elements, that read

1 ∂2
 
h j | T |i i = h j | − |i i . (A19)
2 ∂x2

Taking the first and second derivatives of the harmonic oscillator wave function and using
Equation (A10)
∂
Hn ( x ) = 2nHn−1 ( x ) , (A20)
∂x
∂2
Hn ( x ) = 4n(n − 1) Hn−2 ( x ) , (A21)
∂x2
we can express the kinetic energy as

Ai − x2 /2  
h j | T |i i = h j | e −4i (i − 1) Hi−2 + 4x i Hi−1 + Hi − x2 Hi |χmi i . (A22)
2
Using Equation (A11) we obtain:

1 x2
 
− x2 /2
h j | T |i i = h j | Ai e i+ − Hi ( x ) |χmi i , (A23)
2 2

and taking the results of the harmonic oscillator potential energy Equation (A18), the kinetic energy
matrix elements are expressed as
 q 
1
q
h j | T |i i = − i (i − 1)δj,i−2 + (2i + 1)δj,i − (i + 2)(i + 1)δj,i+2 δmi ,m j . (A24)
4

Obviously, both the kinetic energy and the harmonic oscillator potential are diagonal in the
spin projection.

References
1. Giamarchi, T. Quantum Physics in One Dimension; Clarendon: Oxford, UK, 2004.
2. Tonks, L. The complete equation of state of one, two and three-dimensional gases of hard elastic spheres.
Phys. Rev. 1936, 50, 995. [CrossRef]
3. Girardeau, M. Relationship between Systems of Impenetrable Bosons and Fermions in One Dimension.
J. Math. Phys. 1960, 1, 516. [CrossRef]
4. Kinoshita, T.; Wenger, T.; Weiss, D.S. Observation of a one-dimensional Tonks-Girardeau gas. Science 2004,
305, 1125–1128. [CrossRef] [PubMed]
Mathematics 2020, 8, 1196 36 of 38

5. Paredes, B.; Widera, A.; Murg, V.; Mandel, O.; Fölling, S.; Cirac, I.; Shlyapnikov, G.V.; Hänsch, T.W.; Bloch,
I. Tonks—Girardeau gas of ultracold atoms in an optical lattice. Nature 2004, 429, 277–281. [CrossRef]
[PubMed]
6. Kinoshita, T.; Wenger, T.; Weiss, D.S. Local pair correlations in one-dimensional Bose gases. Phys. Rev. Lett.
2005, 95, 190406. [CrossRef] [PubMed]
7. Dalfovo, F.; Giorgini, S.; Pitaevskii, L.P.; Stringari, S. Theory of Bose-Einstein condensation in trapped gases.
Rev. Mod. Phys. 1999, 71, 463. [CrossRef]
8. Giorgini, S.; Pitaevskii, L.P.; Stringari, S. Theory of ultracold atomic Fermi gases. Rev. Mod. Phys. 2008,
80, 1215. [CrossRef]
9. Bloch, I.; Dalibard, J.; Zwerger, W. Many-body physics with ultracold gases. Rev. Mod. Phys. 2008, 80, 885.
[CrossRef]
10. Lewenstein, M.; Sanpera, A.; Ahufinger, V. Ultracold Atoms in Optical Lattices: Simulating Quantum Many Body
Physics; Oxford University Press: Oxford, UK, 2012.
11. Serwane, F.; Zürn, G.; Murmann, S.; Brouzos, I.; Lompe, T.; Jochim, S. Deterministic preparation of a tunable
few-fermion system. Science 2011, 332, 336–338. [CrossRef]
12. Zürn, G.; Serwane, F.; Lompe, T.; Wenz, A.N.; Ries, M.G.; Bohn, J.E.; Jochim, S. Fermionization of two
distinguishable fermions. Phys. Rev. Lett. 2012, 108, 075303. [CrossRef] [PubMed]
13. Wenz, A.; Zürn, G.; Murmann, S.; Brouzos, I.; Lompe, T.; Jochim, S. From few to many: Observing the
formation of a Fermi sea one atom at a time. Science 2013, 342, 457–460. [CrossRef] [PubMed]
14. Chin, C.; Grimm, R.; Julienne, P.; Tiesinga, E. Feshbach Resonances in Ultracold Gases. Rev. Mod. Phys. 2010,
82, 1225. [CrossRef]
15. Zürn, G.; Wenz, A.N.; Murmann, S.; Bergschneider, A.; Lompe, T.; Jochim, S. Pairing in few-fermion systems
with attractive interactions. Phys. Rev. Lett. 2013, 111, 175302. [CrossRef]
16. Murmann, S.; Bergschneider, A.; Klinkhamer, V.M.; Zürn, G.; Lompe, T.; Jochim, S. Two fermions in a double
well: Exploring a fundamental building block of the Hubbard model. Phys. Rev. Lett. 2015, 114, 080402.
[CrossRef]
17. Murmann, S.; Deuretzbacher, F.; Zürn, G.; Bjerlin, J.; Reimann, S.M.; Santos, L.; Lompe, T.; Jochim, S.
Antiferromagnetic Heisenberg spin chain of a few cold atoms in a one-dimensional trap. Phys. Rev. Lett.
2015, 115, 215301. [CrossRef] [PubMed]
18. Hammer, H.-W.; Nogga, A.; Schwenk, A. Three-body forces: From cold atoms to nuclei. Rev. Mod. Phys.
2013, 85, 197. [CrossRef]
19. Ring, P.; Schuck, P. The Nuclear Many-Body Problem; Springer: Berlin/Heidelberg, Germany, 1980.
20. Truscott, A.G.; Strecker, K.E.; McAlexander, W.I.; Partridge, G.B.; Hulet, R.G. Observation of Fermi Pressure
in a Gas of Trapped Atoms. Science 2001, 291, 2570. [CrossRef]
21. Blume, D. Few-body physics with ultracold atomic and molecular systems in traps. Rep. Prog. Phys. 2012,
75, 046401. [CrossRef]
22. Mistakidis, S.I.; Katsimiga, G.C.; Koutentakis, G.M.; Schmelcher, P. Repulsive Fermi polarons and their
induced interactions in binary mixtures of ultracold atoms. New J. Phys. 2019, 21, 043032. [CrossRef]
23. Zöllner, S.; Meyer, H.-D.; Schmelcher, P. Correlations in ultracold trapped few-boson systems: Transition
from condensation to fermionization. Phys. Rev. A 2006, 74, 063611. [CrossRef]
24. Cheiney, P.; Cabrera, C.R.; Sanz, J.; Naylor, B.; Tanzi, L.; Tarruell, L. Bright Soliton to Quantum Droplet
Transition in a Mixture of Bose-Einstein Condensates. Phys. Rev. Lett. 2018, 120, 135301. [CrossRef] [PubMed]
25. Sowinski, T.; Garcia-March, M.A. One-dimensional mixtures of several ultracold atoms: A review.
Rep. Prog. Phys. 2019, 82, 104401. [CrossRef]
26. Laird, E.K.; Shi, Z.-Y.; Parish, M.M.; Levinsen, J. SU(N) fermions in a one-dimensional harmonic trap.
Phys. Rev. A 2017, 96, 032701. [CrossRef]
27. Brouzos, I.; Schmelcher, P. Two-component few-fermion mixtures in a one-dimensional trap: Numerical
versus analytical approach. Phys. Rev. A 2013, 87, 023605. [CrossRef]
28. Lindgren, E.J.; Rotureau, J.; Forssén, C.; Volosniev, A.G.; Zinner, N.T. Fermionization of two-component
few-fermion systems in a one-dimensional harmonic trap. New J. Phys. 2014, 16 063003. [CrossRef]
29. Andersen, M.E.S.; Dehkharghani, A.S.; Volosniev, A.G.; Lindgren, E.J.; Zinner, N.T. An interpolatory ansatz
captures the physics of one-dimensional confined Fermi systems. Sci. Rep. 2016, 6, 28362. [CrossRef]
[PubMed]
Mathematics 2020, 8, 1196 37 of 38

30. P˛ecak, D.; Dehkharghani, A.S.; Zinner, N.T.; Sowiński, T. Four fermions in a one-dimensional harmonic trap:
Accuracy of a variational-ansatz approach. Phys. Rev. A 2017, 95, 053632. [CrossRef]
31. Gordillo, M.C. One-dimensional harmonically confined SU(N) fermions. Phys. Rev. A 2019, 100, 023603.
[CrossRef]
32. Harshman, N.L. Spectroscopy for a few atoms harmonically trapped in one dimension. Phys. Rev. A 2014,
89, 033633. [CrossRef]
33. Sowinski, T.; Grass, T.; Dutta, O.; Lewenstein, M. Few interacting fermions in one-dimensional harmonic
trap. Phys. Rev. A 2013, 88, 033607. [CrossRef]
34. Ledesma, D.; Romero-Ros, A.; Polls, A.; Juliá-Díaz, B. Dynamic structure function of two interacting atoms
in 1D. EPL 2019, 127, 56001. [CrossRef]
35. Pyzh, M.; Kronke, S.; Weitenberg, C.; Schmelcher, P. Spectral properties and breathing dynamics of a
few-body Bose-Bose mixture in a 1D harmonic trap. New J. Phys. 2018, 20, 015006. [CrossRef]
36. Olshanii, M. Atomic scattering in the presence of an external confinement and a gas of impenetrable bosons.
Phys. Rev. Lett. 1998, 81, 938. [CrossRef]
37. Haller, E.; Gustavsson, M.; Mark, M.J.; Danzl, J.G.; Hart, R.; Pupillo, G.; Nagerl, H.-C. Realization of an
Excited, Strongly Correlated Quantum Gas Phase. Science 2009, 325, 1224. [CrossRef] [PubMed]
38. Yin, X.Y.; Yan, Y.; Smith, D.H. Dynamics of small trapped one-dimensional Fermi gas under oscillating
magnetic fields. Phys. Rev. A 2016, 94, 043639. [CrossRef]
39. Guan, L.; Chen, S.; Wang, Y.; Ma, Z. Exact solution for infinitely strongly interacting Fermi gases in tight
waveguides. Phys. Rev. Lett. 2009, 102, 160402. [CrossRef]
40. Volosniev, A.G.; Fedorov, D.V.; Jensen, A.S.; Valiente, M.; Zinner, N.T. Strongly interacting confined quantum
systems in one dimension. Nat. Commun. 2014, 5, 5300. [CrossRef]
41. Lieb, E.H.; Mattis, D. Theory of ferromagnetism and the ordering of electronic energy levels. Phys. Rev. 1962,
125, 164. [CrossRef]
42. Decamp, J.; Armagnat, P.; Fang, B.; Albert, M.; Minguzzi, A.; Vignolo, P. Exact density profiles and symmetry
classification for strongly interacting multi-component Fermi gases in tight waveguides. New J. Phys. 2016,
18, 055011. [CrossRef]
43. Dickhoff, W.H.; van Neck, D. Many-Body Theory Exposed! World Scientific: Singapore, 2008.
44. Busch, T.; Englert, B.G.; Rzazewski, K.; Wilkens, M. Two cold atoms in an harmonic trap. Found. Phys. 1998,
28, 549. [CrossRef]
45. Wu, K.; Simon, H. Thick-restart Lanczos method for large symmetric eigenvalue problems. SIAM J. Matrix
Anal. Appl. 2000, 22, 602. [CrossRef]
46. Meyer, H.D.; Manthe, U.; Cederbaum, L.S. The multi-configurational time-dependent Hartree approach.
Chem. Phys. Lett. 1990, 165, 73. [CrossRef]
47. Rammelmüller, L.; Porter, W.J.; Braun, J.; Drut, J. Evolution from few- to many-body physics in
one-dimensional Fermi systems: One- and two-body density matrices and particle-partition entanglement.
Phys. Rev. A 2017, 96, 033635. [CrossRef]
48. Bellotti, F.F.; Dehkharghani, A.S.; Zinner, N.T. Comparing numerical and analytical approaches to strongly
interacting two-component mixtures in one dimensional traps. Eur. Phys. J. D 2017, 71, 37. [CrossRef]
49. Raventos, D.; Grass, T.; Lewenstein, M.; Juliá-Díaz, B. Cold bosons in optical lattices: A tutorial for Exact
Diagonalization. J. Phys. B At. Mol. Opt. Phys. 2017, 50, 113001. [CrossRef]
50. Plodzien, M.; Wiater, D.; Chrostowski, A.; Sowinski, T. Numerically exact approach to few-body problems
far from a perturbative regime. arXiv 2018, arXiv:1803.08387.
51. Chrostowski, A.; Sowiński, T. Efficient construction of many-body Fock states having the lowest energies.
Acta Phys. Pol. A 2020, 136, 566. [CrossRef]
52. Titchmarsh, E.C. Some integral Involving Hermite Polynomials. J. Lond. Math. Soc. 1948, 23, 15. [CrossRef]
53. Grining, T.; Tomza, M.; Lesiuk, M.; Przybytek, M.; Musial, M.; Massignan, P.; Lewenstein, M.; Moszynski,
R. Many interacting fermions in a one-dimensional harmonic trap: A quantum-chemical treatment. New J.
Phys. 2015, 17, 115001. [CrossRef]
54. Gharashi, S.E.; Blume, D. Correlations of the upper branch of 1d harmonically trapped two-component
Fermi gases. Phys. Rev. Lett. 2013, 111, 045302. [CrossRef]
55. Li, X.; Pecak, D.; Sowinski, T.; Sherson, J.; Nielsen, A.E.B. Global optimization for quantum dynamics of
few-fermion systems. Phys. Rev. A 2018, 97, 033602. [CrossRef]
Mathematics 2020, 8, 1196 38 of 38

56. Fang, B.; Carleo, G.; Bouchoule, I. Quench-induced breathing mode of one-dimensional Bose gases.
Phys. Rev. Lett. 2014, 113, 035301. [CrossRef] [PubMed]
57. Menotti, C.; Stringari, S. Collective oscillations of a 1D trapped Bose gas. Phys. Rev. A 2002, 66, 043610.
[CrossRef]
58. Moritz, H.; Stoferle, T.; Kohl, M.; Esslinger, T. Exciting Collective Oscillations in a Trapped 1D Gas.
Phys. Rev. Lett. 2003, 91, 250402. [CrossRef] [PubMed]
59. Kohn, W. Cyclotron Resonance and de Haas-van Alphen Oscillations of an Interacting Electron Gas.
Phys. Rev. 1961, 123, 1242. [CrossRef]
60. Ebert, M.; Volosniev, A.; Hammer, H.-W. Two Cold Atoms in a Time-Dependent Harmonic Trap in One
Dimension. Ann. Phys. 2016, 528, 698. [CrossRef]
61. Gharashi, S.E.; Blume, D. Broken scale-invariance in time-dependent trapping potentials. Phys. Rev. A 2016,
94, 063639. [CrossRef]
62. Kwasniok, J.; Mistakidis, S.I.; Schmelcher, P. Correlated dynamics of fermionic impurities induced by the
counterflow of an ensemble of fermions. Phys. Rev. A 2020, 101, 053619. [CrossRef]
63. Gudyma, A.I.; Astrakharchik, G.E.; Zvonarev, M.B. Reentrant behavior of the breathing-mode-oscillation
frequency in a one-dimensional. Phys. Rev. A 2015, 92, 021601. [CrossRef]
64. Bohigas, O.; Lane, A.M.; Martorell, J. Sum rules for nuclear collective excitations. Phys. Rep. 1979, 51, 267.
[CrossRef]
65. Sowinski, T.; Brewczyk, M.; Gajda, M.; Rzazewski, K. Exact dynamics and decoherence of two cold bosons
in a 1D harmonic trap. Phys. Rev. A 2010, 82, 053631. [CrossRef]
66. Budewig, L.; Mistakidis, S.I.; Schmelcher, P. Quench Dynamics of Two One-Dimensional Harmonically
Trapped Bosons Bridging Attraction and Repulsion. Mol. Phys. 2019, 117, 2043. [CrossRef]
67. Erdmann, J.; Mistakidis, S.I.; Schmelcher, P. Phase-separation dynamics induced by an interaction quench of
a correlated Fermi-Fermi mixture in a double well. Phys. Rev. A 2019, 99, 013605. [CrossRef]

c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).