Static and dynamic properties of a few spin $1/2$ interacting fermions trapped in an harmonic potential

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 exact 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.


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][5][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][8][9][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 phenomena in few-fermion systems [15], the behaviour of two fermions in a double-well potential [16], the realization of 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 very versatile laboratory where 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 to mention 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][21][22][23].
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 [20,[24][25][26]. 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 [27,28].
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 exact 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

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.

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 [26,28] where H ho = ∑ i h ho (i) is the sum of single-particle harmonic oscillator hamiltonians. The Hamiltonian, h ho , is a one-body operator, which includes the kinetic-and the harmonic-potential energy, where m is the particle mass and ω is the harmonic oscillator frequency. The eigenvalues of this Hamiltonian are well known, e n = n + 1 2 hω n = 0, 1, 2, ... , (3) and the corresponding eigenfunctions are written as, where H n (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 an 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, m , m = 1/2, −1/2, where m is the spin projection on the z-axis. Sometimes, these two spin states are denoted by: |↑ and |↓ , respectively. The fermions are assumed to interact via a contact spin-independent delta potential [29]. 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 to be at the same position. Therefore the 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 Eq. (1) is a contact interaction, which is expressed as [30] where g characterizes the strength of the interaction and δ is a Dirac delta function. 6 of 44 Taking into account the identity [27] where X = 1 N ∑ i x i is the center of mass coordinate, the Hamiltonian can be split in two pieces, where H CM in harmonic oscillator units reads: and governs the center-of-mass motion. H r 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 (a 1d ) as [29,31], From the relation between the one-dimensional scattering length and the three-dimensional scattering length [32], we can write the coupling constant as 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 [32].
For convenience we use the harmonic oscillator units, in which the energy is measured inhω 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.

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 [33].
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.

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 to have two fermions with the same spin projection in the same harmonic oscillator state. More specifically, in the ground state of a system with N d particles of spin down and N u particles of spin up, the particles occupy the lowest single-particle energy states, N d and N u , respectively. Therefore, the energy of the ground state is Actually the wave function of the system is given by a Slater determinant built with the lowest N u and N d single-particle wave functions with spin up and down, respectively. The total third component of the spin of this wave function is M = N u −N d 2 and the total spin is the minimum S compatible with M, S = M [34]. For a non-polarized system, N u = N d , 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, where Φ n (x) are the 1D harmonic oscillator wave functions, Eq. (4). The density profile is normalized to the number of particles.

Infinite interaction case: Tonks-Girardeau gas
In the infinite interaction limit, the system experiences a phenomenon similar to the fermionization process in Bose systems [2,3]. In fact, a system of fermions under this situation, are described as a system of non-interacting identical fermions [26]. This system is known as a Tonks-Girardeau gas. In our case, the wave function can be written as [33], where Ψ A is an antisymmetric wave function on the spatial coordinates and Ψ S is a symmetric wave function which depends on both spatial and spin coordinates. The antisymmetric part, Ψ A , is a Slater determinant built with the N = N d + N u harmonic oscillator wave functions. And the symmetric part, Ψ S , is a linear combination of products of sign functions of differences of the particle coordinates and their spin functions. The energy of the ground state depends only on the single-particle states present in The different Ψ S functions associated with the same Ψ A define the degeneracy of the ground state: [33]. The total third spin component of this wave function is However, different total spin functions coexist in the degenerated ground state.
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, where Φ n (x) are the 1D harmonic-oscillator wave functions, Eq. (4).

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 [35].
A particular state is the vacuum state denoted by |0 . This vector represents a state without particles, i.e., all its occupation numbers are zero. This state has norm equal to 1.

Creation and annihilation operators
The main tool in second quantization are the creation (a † i ) and annihilation (a i ) operators, which act on the Fock space. Any observable, which is represented by an operator, can be expressed in terms of them [35]. The action of a creation operator on the vacuum state, denoted as a † i |0 , 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 (a i ) 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.
Also, a Fock space vector with N particles can be created by the action of N creation operators on the vacuum space, as follows, with N = ∑ i n i . 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 Therefore the occupation number of the Fock basis can take the values n i = 0 or n i = 1. Finally, the action of the creation and annihilation operators can be written as where N i is the number of occupied states with index lower than i, The phase (−1) N i 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.

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: where n i,m indicates the number of particles in the ϕ i,m single-particle state.

Operators in second quantization
In second quantization, the operators can be expressed in terms of creation and annihilation operators. A general one-body operatorÔ is expressed aŝ where |i and |j are single-particle states, i|Ô |j = ϕ * i (x)O(x)ϕ j (x)dx is the one-body matrix 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 operatorV is expressed aŝ where v ij,kl is the two-body matrix element defined as

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 Eq. (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 ǫ i = i + 1/2. Therefore in second quantization it reads,ĥ wheren i = a † i a i is the number operator associated to the single-particle state |i . For the interaction term, the two-body matrix elements are expressed as Notice that the interaction does not affect the spin of the particles, and we have used the orthogonality of the spin functions: χ m i χ m j χ m k χ m l = δ m i ,m k δ m j ,m l . 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

Numerical methods
In this chapter we describe the method named exact diagonalization used to solve 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 [36], as a benchmark of our numerical procedure.

Exact 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 an exact diagonalization method. 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 [37]. This method allows us to diagonalize large matrices, and obtain the lowest eigenvalues with high accuracy.
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 [26]. 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.

Basis truncation
To obtain the exact results, when using a diagonalization method, we would need to use the 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, like finite optical lattices [38], 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 [38]. In our case, we diagonalize the Hamiltonian in a subspace with well-defined total third spin component. We start considering the lowest N M single-particle states and take into account an energy constrain 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 E max , this procedure is described in [39].
This maximum energy depends on the number of particles and spin configuration. Notice that we consider always M > 0 therefore, N u > N d . The maximum energy, E max , 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 (N u − 1) spin-up particles in the lowest N u − 1 single-particle states. On the other hand the N d spin-down particles are located in the lowest N d single-particle states. Therefore, the maximum energy considered in the construction of the many-body basis is where N M 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. Note that the dimension of the many-body basis in the case without energy restriction is which grows very rapidly. Notice also that 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, with all spins parallel and an antisymmetric spatial function built with the lowest single-particle energy wave functions. As it has been said previously, these configurations do not feel the effects of a contact interaction.

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, Eq. (25). These matrix elements contain an integral with four wave functions, The integral, that we label as I abcd can be solved analytically, because the wave functions are the harmonic oscillator wave functions Eq. (4), Let us consider, 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. and We can express the integral I ′ abcd as 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, Finally, the integral I abcd can be expressed as The presence of factorials in the analytical expression of the matrix elements Eq. (36), can originate numerical problems. In this section we propose a procedure to overcome this problem.
In general the matrix elements of the interaction should be well behaved. However, the presence of factorials in the expressions of the harmonic oscillator wave functions can originate 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 This expression can be written as where f is a well-behaved function. Furthermore, we can write this result as where f must be positive.
The logarithm log ( f ) is given by 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 exponentiation. Thus, the final expression of the integral taking into account the possible negative values of gamma functions reads where p is the phase, given by p = ∏ 3 n=1 p n , with p n being the phase generated by each gamma function. The logarithm of the gamma functions and the associated phases are calculated as The logarithm of a factorial is calculated as:

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 [36].

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 [36] where Γ are gamma functions, E r 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 Eq. (8), and its energy is given by Eq. (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.

Comparison of analytical and numerical results
In Fig. 1 we report as a function of g the low-energy part of the two-particle energy spectrum calculated by exact diagonalization (solid green line) and the spectrum obtained by solving Eq. (43). Notice that this equation provides the energies of the relative motion to which we add the possible energies of the center of mass: E CM . For the ground state, E CM = 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 in a dimension of 5050 of the matrix to be diagonalized.
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 strength interaction 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 state with S = 1 and M = 1 which 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 degenerated 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 (N M ). Note that for N = 2, the energy constraint for the construction of the two-body basis is given by the maximum To this end, in Fig. 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 N M , which is given by N M (N M + 1)/2. As expected, the difference between the calculated and the expected value decreases with the number of modes. Following Ref. [41], we fit functions of the type C 1 /N 1/2 M + C 2 /N M with excellent results. In addition, for large values of N M , the convergence of the energy goes as 1/N 1/2 M . 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 are very small. In all cases the difference is positive, indicating that the numerical results are upperbounds to the exact ones.

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 number 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 where |ψ n 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.

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 Fig. 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 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 the Eq. (14), E = (N d + N u ) 2 /2 where N d and N u 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.
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 Fig. 3. In the non-interacting case, the wave function can be factorized 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 antisymmetric and its energy is E = 1. In the limit of a very strong interaction, the wave function can be factorized as [3]  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 Eq. (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: The fulfilment of this relation reinforces the accuracy of the calculations. The derivation of the virial relation can be found in the Appendix A.
In Fig. 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. 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 Eq. (46), that prevents to have two particles at 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.

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 |Ψ are defined as,

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, 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 Fig. 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 Eq. (12). For large values of g, g = 14 in the figure, the density profiles are very close to the ones predicted in Eq. (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.

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 a 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 mean-field wave function with largest overlap with the ground state. However, this condition does not imply that the energy corresponding to this mean-field wave function is the minimum energy for a mean-field wave function. 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.
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 Fig. 6 we report the dependence of the ten larger 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, being N 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 degenerated 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. The reason behing the previous observations is that the OBDM is constructed in boxes of well defined third spin component, 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.

Low-energy excited states
To complete the study of the structure of the few-body fermions 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.

Energy spectrum
Once the Hamiltoninan 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 Fig. 7 for several number of 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 symmetric wave function. Later on, in Sect. 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 with energy E = 1. The two first excited levels have the same energy, E = 2 at g = 0, one corresponds to a state which is a product of antisymmetric 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 degenerated with states having the same spatial wave function and the spin functions corresponding to S = 1 but M = 0, ±1. The other excited state at g = 0 corresponds to 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 H CM . However, in the state we are considering the center of mass occupies the first excited state of H CM . This is true for all values of g and the energy of this state can be written as E = E g.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 H CM 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 E g→∞ = 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, H CM . There are other levels which 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. Then as the Hamiltonian commutes with S − this operation does not change the energy of the state. These arguments are discussed at length in Sect. 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 has an energy E = 13/2 and when g → ∞ merges with the state non-sensitive to the interaction discussed above, which an energy E = 25/2.
For the case with M = 3/2 there are much less 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 degenerated 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. [26] 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. [26].
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 degeneration. In particular, for the ground state this degeneration is expressed as D = N! N u !N d ! .

Spin determination
As the Hamiltonian does not depend on the spin, it commutes with the spin operators. In particular 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 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. One can then assign the spin quantum number to the remaining states, which is the minimum S compatible with the value of M, i.e. S = N/2 − 1. Notice that, 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 an 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 Fig. 8 as a function of the interaction strength. The maximum possible spin projection is N/2 = 2. The wave function of this state can be factorized as a Slater determinant with the first four single particle states of the harmonic oscillator times a spin wave function with all spins pointing up. The energy of this wave function, 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 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 type 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.

Dynamical excitation
In Sects. 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 [42][43][44]. 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 [45]. 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 are not excited and we look for traces of the internal structure on the time evolution of the central density of the cloud.

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 [46], in our case we study the response of the system to a change in the trapping frequency. The main tool we consider is the dynamic structure function associated to the mono-polar excitation operator 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.
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.

Dynamic structure function
The dynamic structure function, normalized to the number of particles, is defined as whereF is the excitation operator associated to the external perturbation, see Eq. (54), |Ψ 0 is the ground state of the unperturbed system and |Ψ i are its excited states.
In the second quantization formalism, the structure function reads, where |k and |l are harmonic-oscillator single-particle states and |ψ n are the Fock states of our basis. For the center of mass, the i CM | X 2 CM |j CM are different from zero only when i CM = j CM or i CM = j CM ± 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 [27]. For the non-interacting case, we expect a single peak in the dynamic structure function with an energy E − E 0 = 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 − E 0 = 2. In addition, a peak with energy E − E 0 = 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 Fig. 9 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 monopolar operator have total S = 0 and even parity. As expected, for the non-interacting case we only have one peak at excitation energy E − E 0 = 2. For any interaction strength there is a peak at E − E 0 = 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 − E 0 = 2, and the rest of 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 reentrant behaviour, reported for bosons in Ref. [47], 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), where in Eq. (56) is expressed in terms of single-particle matrix elements and the Fock basis. It is also worth to mention that in this case S F (E) with for two particles with M = 0 is the same as if the particles were two bosons Ref. [27]. 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 [36]. 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. In Fig. 10 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 − E 0 = 2. In the interacting case, there is a peak at energy E − E 0 = 2 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 Sect. 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 reentrant behaviour. This is common to all the cases with 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.

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 where E is the excitation energy. They can be explicitly computed from the definition of S F (E), Eq. (56), by using the second quantization formalism as: If the dynamic structure function is known, the energy momenta can be computed directly using Eq. (57). But in general, 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 M n using only ground state properties [48]. These relations are named sum rules. For the monopolar excitation operator, the lowest sum rules can be simply expressed as the following expectation values on the ground state [27] Using perturbation theory, we have an alternative expression for M −1 , where |0 and E 0 (λ) are the ground state and its corresponding energy of the perturbed Hamiltonian H ′ =Ĥ + λF.
In Fig. 11 we report the energy momenta M −1 , M 1 and M 3 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 M 1 and M 3 , 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.
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 [47]. Indeed, if only one peak appears in the dynamic structure function, we can compute its excitation energy by M n /M n−2 [48]. In a more general case this ratio can at least be used as an estimation of the energy of the main peak.
In Fig. 12 we report the ratios M 1 /M −1 , √ M 3 /M 1 and excitation energy of the main contribution, which is associated to an intrinsic excitation, as a function of the interaction strength. We can observe that in the case of two particles, both estimations, M 1 /M −1 and √ M 3 /M 1 , 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 which 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 reentrance 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 Fig. 9. 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 Fig. 10, 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.

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. 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 protocols 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 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.

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 g new . The initial ground state is not anymore 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: where |Ψ(t = 0) is the ground state of the Hamiltonian before the quenching, and the states |Ψ ′ n are the eigenfunctions of the Hamiltonian after the quenching with eigenenergies E ′ n . The coefficients c n can be computed as: The time evolution of the state after the quenching is calculated by taking into account that the states |Ψ ′ n are the eigenstates of the new Hamiltonian, In turn, the eigenstates |ψ ′ n can be expressed in terms of the many-body basis used to describe our system: |Ψ where |ψ 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: where Φ i (x) are the single-particle harmonic oscillator wave functions. Notice that the time evolution of the system is governed by E ′ m − E ′ n , involving all energy differences of the states of the new Hamiltonian that have non-zero overlap with the ground-state of the old Hamiltonian. Also important is to realize that the relevance of a given frequency 2πν n,m = E n − E m depends on the product of the projection coefficients c * n c m . 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 one. In this situation the dominant frequencies will be associated to the differences E ′ n − E ′ 0 . In this sense, the frequencies would be similar to the frequencies 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. In Fig. 13, we report, for the two-particle system, the projections of the ground state of the non-interacting system to the eigenstates of H ′ (g = 1). In this quench, the interaction strength suffers a sudden change from g = 0 to g = 1. We also report the projection for the quench g = 0 → 5. Notice that not all the eigenstates of H ′ (g = 1) are populated, one should consider only projections to states of H ′ with the same parity and total spin and third spin component. Besides, as the interaction is translationally invariant, the quench of the interaction does not affect the center-of-mass motion and the center of mass remains in the lowest state. In the case g = 1 (small g) there is a dominant projection into the ground state of 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.  Figure 13. 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.

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 Fig. 14. 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 Fig. 14. 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. 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 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. All magnitudes in the figure are in harmonic oscillator units.
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 reentrance phenomena, 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 to mention 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 Fig. 7). No center-of-mass excitations affect this analysis.
In the same figure, Fig. 14, 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 Fig. 8, 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. Finally, it is worth to mention that 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.

Summary and conclusions
In this paper we have presented a quantum microscopic description of few 1/2 spin 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 in repulsive interactions. This type 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 the 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 identify the excited states. The spin assignation was useful to identify the excited states by the breathing mode. We have also calculated the energy momenta (M n ) associated to the dynamic structure function and used the ratios M n /M n−2 to estimate the average excitation energies. The differences between the ratios M n /M n−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 estimulate 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.