Exact density profiles and symmetry classification for strongly interacting multi-component Fermi gases in tight waveguides

Ultracold atomic gases made with rare-Earth elements cooled to quantum degeneracy and subjected to two-dimensional optical lattices provide a beautiful realization of the model of one-dimensional multicomponent Fermi gases with strong and equal intercomponent repulsion among the species [1].

In the absence of external harmonic confinement in the longitudinal direction, the system of multi-component Fermions with intercomponent delta interactions is a generalization of the Yang–Gaudin Hamiltonian [2–4] and can be solved by nested Bethe Ansatz [5]. The homogeneous system may also be described at low energy as a multi-component Luttinger liquid. This model has been extensively studied in the context of electronic multichannel systems [6, 7] and exotic condensed matter materials [7–11].

The experiments with ultracold atoms are characterized by the presence of an external longitudinal confinement, which can be well approximated as harmonic. In this case, the Bethe Ansatz solutions do not apply, however, a special integrable case is provided by the limit of infinitely strong repulsions among the species. In this regime, corresponding to the Tonks-Girardeau regime for ultracold bosonic atoms, Fermions belonging to different components further Fermionize, as has been experimentally shown in [12]: the wavefunction vanishes at contact and they can be mapped onto a noninteracting Fermi gas in the same external confinement with particle number corresponding to the total number of Fermions in the mixture, following the same idea as the original Girardeau solution for bosons [13]. An additional difficulty in the multicomponent case is that the manifold displays a large degeneracy [14]. This is associated to the arbitrariness in fixing the relative phase once two Fermions belonging to different components get in contact and then exchange their mutual position. The degeneracy is however broken at finite interactions, where one expects a unique ground state [15, 16]. The ground state branch can be obtained by performing a strong-coupling $1/g$ expansion, g being the interaction strength among the Fermions. This provides a unique way to identify in the degenerate manifold what would be the ground and excited states at finite interactions [16, 17].

In this work, we explore this type of solution for multicomponent Fermi mixtures with a number of components that ranges from two to six as in the ¹⁷³Yb experiment of [1]. We study in particular the density profiles of the imbalanced mixtures, thus generalizing the works of [18, 19]. As main result, we find a complex, inhomogeneous spatial structure. We then explore the symmetry property of both the ground and excited state branches of the solution, introducing and evaluating suitable operators for the mixture, i.e. the transposition class-sum and the three-cycle class sum operators [20, 21]. This allows us to test (and generalize) the Lieb–Mattis theorem [22] for the trapped multicomponent mixtures, showing that the ground state carries the most symmetric configuration allowed by the imbalance among the components.

We consider a system of N Fermions of equal mass m, divided in r species with population ${N}_{1},{N}_{2},\ldots ,{N}_{r}$. We assume that all components are subjected to the same harmonic potential $V(x)=m{\omega }^{2}{x}^{2}/2$, as is the case of Fermions in optical traps. The Fermions belonging to different species interact with each other via the contact potential $v(x-x^{\prime} )=g\delta (x-x^{\prime} )$, where g is the interaction strength, and $\delta (x)$ is the Dirac delta function. The total Hamiltonian reads

$$\begin{eqnarray}&&{ \mathcal H }=\displaystyle \sum _{i=1}^{N}\left[-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}^{2}}+\displaystyle \frac{1}{2}\;m{\omega }^{2}{x}_{i}^{2}\right]+g\displaystyle \sum _{i\lt j}\delta ({x}_{i}-{x}_{j}),\end{eqnarray} \tag{ 1 }$$

where ${x}_{1},\ldots ,{x}_{N}$ are the coordinates of the Fermions. The effect of contact interactions can be replaced by a cusp condition on the many-body wavefunction:

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}}{2\;{mg}}\left[{\left(\displaystyle \frac{\partial {\rm{\Psi }}}{\partial {x}_{i}}-\displaystyle \frac{\partial {\rm{\Psi }}}{\partial {x}_{j}}\right)}_{{x}_{i}-{x}_{j}\to {0}^{+}}-{\left(\displaystyle \frac{\partial {\rm{\Psi }}}{\partial {x}_{i}}-\displaystyle \frac{\partial {\rm{\Psi }}}{\partial {x}_{j}}\right)}_{{x}_{i}-{x}_{j}\to {0}^{-}}\right]={\rm{\Psi }}({x}_{i}={x}_{j}).\end{eqnarray} \tag{ 2 }$$

In the following, we will focus on the impenetrable limit $g\to \infty$. In this case the cusp condition imposes that the many-body wavefunction vanishes when two particles meet, i.e. ${\rm{\Psi }}({x}_{j}={x}_{{\ell }})=0$ for each pair $\{j,{\ell }\}$. This condition is exactly satistisfied by the fully antisymmetric solution ${{\rm{\Psi }}}_{A}({x}_{1},\ldots ,{x}_{N})$ of $N={N}_{1}+{N}_{2}\quad +\cdots +\quad {N}_{r}$ noninteracting Fermions in the same confining potential, corresponding to a Slater determinant constructed with single-particle wavefunctions ${\phi }_{0},\ldots ,{\phi }_{N-1}$. To construct the exact solution, we further note that the behavior of the many-body wavefunction under exchange between two Fermions belonging to different components is not fixed by symmetry, and requires to be fixed by additional conditions. Hence, we consider a general solution of the form [17, 23]

$$\begin{eqnarray}&&{\rm{\Psi }}({x}_{1},\ldots ,{x}_{N})=\displaystyle \sum _{P\in {S}_{N}}{a}_{P}\chi ({x}_{P(1)}\lt \cdots \lt \quad {x}_{P(N)}){{\rm{\Psi }}}_{A}({x}_{1},\ldots ,{x}_{N}),\end{eqnarray} \tag{ 3 }$$

where $\chi ({x}_{1}\lt \cdots \lt \quad {x}_{N})$ is the indicator function of the sector $\{{x}_{1}\lt \cdots \lt \quad {x}_{N}\}\subset {{\mathbb{R}}}^{N}$, i.e. it is 1 within the sector and 0 everywhere else and S_N the permutation group of N elements. Here, both Ψ and ${{\rm{\Psi }}}_{A}$ are assumed to have unit normalization, and the choice of the coefficients a_P will be detailed below. Note that we need only to determine $S=\tfrac{N!}{{N}_{1}!\ldots {N}_{r}!}$ coefficients a_P: the wavefunction is antisymmetric under exchange of Fermions belonging to the same component, and this property is already encoded in ${{\rm{\Psi }}}_{A}$. This observation allows us to restrict ourself to the so-called snippet basis [23, 24], i.e. consider only the global permutations modulo the permutations of particles belonging to the same species.

General solution

In order to determine the a_P coefficients for the ground state manifold, denoted ${a}_{1},\ldots ,{a}_{S}$ in the snippet basis, we use the same method as in [17]. It consists in a perturbative expansion of the energy to first order in $1/g\to 0$, i.e. we write $E={E}_{A}-({{\hslash }}^{4}/{m}^{2})K/g,$ where E_A is the energy associated with ${{\rm{\Psi }}}_{A}$ and $K=-({m}^{2}/{{\hslash }}^{4})(\partial E/\partial {g}^{-1})$ is proportional to the interaction energy, related to the Tan's contact coefficient in the two-component case [25]. We then use the Hellmann–Feynman theorem and the cusp condition to write

$$\begin{eqnarray}&&K=\displaystyle \frac{1}{4\langle {\rm{\Psi }}| {\rm{\Psi }}\rangle }\displaystyle \sum _{i\lt j}\int {\rm{d}}{x}_{1}\ldots {\rm{d}}{x}_{N}\delta ({x}_{i}-{x}_{j}){\left[{\left(\displaystyle \frac{\partial {\rm{\Psi }}}{\partial {x}_{i}}-\displaystyle \frac{\partial {\rm{\Psi }}}{\partial {x}_{j}}\right)}_{{x}_{i}-{x}_{j}\to {0}^{+}}-{\left(\displaystyle \frac{\partial {\rm{\Psi }}}{\partial {x}_{i}}-\displaystyle \frac{\partial {\rm{\Psi }}}{\partial {x}_{j}}\right)}_{{x}_{i}-{x}_{j}\to {0}^{-}}\right]}^{2},\end{eqnarray} \tag{ 4 }$$

where we have used the natural units, namely the harmonic oscillator length ${a}_{\mathrm{ho}}=\sqrt{{\hslash }/m\omega }$ as unit length and the harmonic oscillator energy ${\hslash }\omega$ as unit energy. Separating the integral over the different sectors $\{{x}_{P(1)}\lt \cdots \lt \quad {x}_{P(n)}\}\subset {{\mathbb{R}}}^{N}$ and recalling that ${{\rm{\Psi }}}_{A}$ is normalized to one, we finally obtain

$$\begin{eqnarray}&&K=\displaystyle \frac{{\displaystyle \sum }_{P,Q\in {S}_{N}}{({a}_{P}-{a}_{Q})}^{2}{\alpha }_{P,Q}}{\tfrac{1}{N!}{\displaystyle \sum }_{P\in {S}_{N}}{a}_{P}^{2}},\end{eqnarray} \tag{ 5 }$$

where the ${\alpha }_{P,Q}$ coefficients are non-zero if the sectors P and Q differ only by transposing two adjacent coordinates, and in this case we have, using permutational symmetry

$$\begin{eqnarray}&&{\alpha }_{P,Q}={\alpha }_{k}={\int }_{{x}_{1}\lt \cdots \lt {x}_{N}}{\rm{d}}{x}_{1}\ldots {\rm{d}}{x}_{N}\delta ({x}_{k}-{x}_{k+1}){\left[\displaystyle \frac{\partial {{\rm{\Psi }}}_{A}}{\partial {x}_{k}}\right]}^{2},\end{eqnarray} \tag{ 6 }$$

for $k\in \{1,\ldots ,N-1\}$. Intuitively, ${\alpha }_{k}$ can be seen as the energy cost of an exchange between particles of different species at positions k and $k+1$. Note that, thanks to the parity invariance of this problem, we also have ${\alpha }_{k}={\alpha }_{N-k}$, so that we have only $\lfloor N/2\rfloor$ coefficients to compute.

In order to find the wavefunction which corresponds to the ground state at finite, large interactions, the next step is to find the solutions that minimize the energy, i.e. that maximize K. To do so, we impose that $(\partial K/\partial {a}_{i})=0$ for all a_i. This turns out to be equivalent to the diagonalization problem

$$\begin{eqnarray}&&{\bf{V}}\vec{A}=K\vec{A},\end{eqnarray} \tag{ 7 }$$

with $\vec{A}={({a}_{1},\ldots ,{a}_{S})}^{{\rm{T}}}$ and ${\bf{V}}$ is a S × S matrix depending only on the ${\alpha }_{k}$ coefficients. More precisely, the ${\bf{V}}$ matrix is defined in the snippet basis by

$$\begin{eqnarray}{V}_{{ij}}=\left\{\begin{array}{ll}{\rm{-}}{\alpha }_{i,j} & {\rm{if}}\;{\rm{}}i\ne j,\\ \displaystyle \sum _{k\ne i}{\alpha }_{i,k} & {\rm{if}}\;{\rm{}}i=j,\end{array}\right.\end{eqnarray} \tag{ 8 }$$

where the α coefficents are defined as in equation (6) (see also [26, 27]).

In order to compute the ${\alpha }_{k}$ coefficients, we use the following expression for ${{\rm{\Psi }}}_{A}$, based on a Vandermonde determinant result [28, 29]

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{A}({x}_{1},\ldots ,{x}_{N})=\displaystyle \frac{1}{\sqrt{N!}}{C}_{N}\displaystyle \prod _{k=1}^{N}{{\rm{e}}}^{-{x}_{k}^{2}/2}\displaystyle \prod _{1\leqslant i\lt j\leqslant N}({x}_{i}-{x}_{j}),\end{eqnarray} \tag{ 9 }$$

where ${C}_{N}=\tfrac{1}{{\pi }^{N/4}{\prod }_{k=0}^{N-1}\sqrt{{2}^{k}k!}}$. We then have after some algebra

$$\begin{eqnarray}&&{\left[\displaystyle \frac{\partial {{\rm{\Psi }}}_{A}}{\partial {x}_{k}}\right]}_{{x}_{k}={x}_{k+1}}=\displaystyle \frac{1}{\sqrt{N!}}{C}_{N}{{\rm{e}}}^{-{x}_{k}^{2}}\displaystyle \prod _{\underset{i\ne k,k+1}{i=1}}^{N}{{\rm{e}}}^{-{x}_{i}^{2}/2}{({x}_{i}-{x}_{k})}^{2}\displaystyle \prod _{\begin{array}{c}1\leqslant j\lt {\ell }\leqslant N\\ j,l\ne k,k+1\end{array}}({x}_{j}-{x}_{{\ell }})\end{eqnarray} \tag{ 10 }$$

and thus, using again the Vandermonde formula

$$\begin{eqnarray}&&{\left[\displaystyle \frac{\partial {{\rm{\Psi }}}_{A}}{\partial {x}_{k}}\right]}_{{x}_{k}={x}_{k+1}}^{2}=\displaystyle \frac{{2}^{2N-3}{{\rm{e}}}^{-2{x}_{k}^{2}}}{\pi N!(N-1)!(N-2)!}{[\mathrm{det}[{\phi }_{i-1}({x}_{j})]]}_{N-2\times N-2}^{2}\displaystyle \prod _{\begin{array}{c}i=1\\ i\ne k,k+1\end{array}}^{N}{({x}_{i}-{x}_{k})}^{4}.\end{eqnarray} \tag{ 11 }$$

Finally, permutation and parity invariances yield

$$\begin{eqnarray}{\alpha }_{k} & = & \displaystyle \frac{{2}^{2N-3}}{\pi N!(N-1)!(N-2)!}\displaystyle \frac{1}{(k-1)!(N-k-1)!}{\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}{x}_{k}{{\rm{e}}}^{-2{x}_{k}^{2}}\displaystyle \sum _{P,Q\in {S}_{N-2}}\epsilon (P)\epsilon (Q)\\ & & \times \displaystyle \prod _{\underset{i\ne k,k+1}{i=1}}^{N}{\displaystyle \int }_{{L}_{k}}^{{U}_{k}}{\rm{d}}{x}_{i}{({x}_{i}-{x}_{k})}^{4}{\phi }_{P(i)-1}({x}_{i}){\phi }_{Q(i)-1}({x}_{i}),\end{eqnarray} \tag{ 12 }$$

where the integration limits are $({L}_{k},{U}_{k})=(-\infty ,{x}_{k})$ if $i\lt k$ and $({x}_{k},+\infty )$ otherwise. An alternative derivation of the coefficients ${\alpha }_{k}$ can be found in [26].

As the first application, using equation (3) together with the solution of equation (7) for the coefficients a_P and equation (12) for the weigths ${\alpha }_{k}$ we determine the exact density profile of each component of the mixture, according to the definition

$$\begin{eqnarray}&&{n}_{\nu }(x)={N}_{\nu }\int \displaystyle \prod _{j=1}^{N}{\rm{d}}{x}_{j}\delta (x-{x}_{\nu }){{\rm{\Psi }}}^{\ast }({x}_{1},\ldots {x}_{N}){\rm{\Psi }}({x}_{1},\ldots {x}_{N}),\end{eqnarray} \tag{ 13 }$$

where we have indicated by x_ν one of the coordinates corresponding to a Fermion belonging to the νth component of the mixture. In the following, for each given mixture, we focus on the ground state and first many-body excited state with a symmetry different from the ground state one.

Quite generally, for all the ground state density profiles we find that the total density $n(x)={\sum }_{\nu =1}^{r}{n}_{\nu }(x)$ coincides with the one of a noninteracting N-particle Fermi gas under external confinement, as previously reported [30]. For equal populations in the various species , the ground state density profile is the same for all the species and coincides apart to a normalization factor with the one of a noninteracting Fermi gas with N particles in harmonic confinement, which is characterized by N peaks [31], as shown in figure 1. The excited states display a variety of profiles, where in particular the two-component mixture does not change and the six-component one shows small deviations with respect to the corresponding ground-state density profiles, while the three-component mixture displays a different, two-peak structure. These differences may be accounted for by considering the different symmetry of the excited states in the three cases, see section 3.

Zoom In Zoom Out Reset image size

In the case of an imbalanced mixture, the partial density profiles display a rich structure, as presented in figure 2. We first consider the case of an imbalanced two-component mixture. In the polaron case, where ${N}_{1}=N-1$ and ${N}_{2}=1$ (top panels in figure 2), for the ground-state density we observe a spatial separation of the polaron density profile at the center of the trap and the majority component at the wings of the trap. By comparing with the results for the corresponding noninteracting gas, one clearly sees the effect of repulsive interactions: the mutual repulsion among the two components push the majority component to a larger region of the trap and a hole is created around the polaron. For the excited state, we find that the density profiles are proportional to the ones of a noninteracting gas. We understand this as being related to the symmetry properties of this particular excited state—as it will be discussed in detail in section 3 below, we find that it has the same symmetry as the noninteracting gas.

Zoom In Zoom Out Reset image size

In the case of more than one particle in the minority component (central panels in figure 2) for the ground state density profiles, we observe a partial demixing through a more complex structure, with the majority component occupying both the inner core and the external wings of the density profiles. This may be viewed as a mesoscopic realization of an antiferromagnetic configuration [32]. The excited state density profiles have instead a two-peak (majority component) or one-peak (minority component) structure which recalls the ground state of the polaron. Also in this case, the analysis of the symmetry of the profiles brings an explanation, since we find that this excited state has the same symmetry as the ground state of the polaron.

The three component imbalanced mixture (bottom panels in figure 2) displays an even more complex spatial-separated shell structure in the ground state profiles, generalizing the two component case: the minority component occupies the inner of the trap, and the majority component the outer shells, with the intermediate component placed spatially between the other two. For the excited state, we find again a mixed state whose density profiles are proportional to the ones of a noninteracting gas. As it will be discussed in detail in section 3 below, this is in agreement with the fact that this state has the same symmetry as the ground-state of the two-component balanced mixture.

We now analyze the symmetry properties of the quantum many body states under exchange of particles [33–36]. Obviously, single component free Fermions are described by a totally anti-symmetric wavefunction. In the presence of several components interacting among each other, the wavefunction must be totally antisymmetric under exchange of particles belonging to each component but the total wavefunction has a more involved symmetry under exchange of an arbitrary pair of particles. The actual symmetry of the ground state and the excited states can be deduced from the properties of the permutation group of N elements, S_N [37]. More precisely, one can expand the eigenstates of the quantum system over the different irreducible representations of this group. We follow this route below, and demonstrate that the ground state (and any excited state) has a well defined symmetry, i.e. can be described by a single irreducible representation, identified by a Young tableau.

The previous diagonalization process allows us to obtain a set of S values of $K,{K}_{1}\leqslant \cdots \leqslant \quad {K}_{S}$, with ${K}_{1}=0$ and ${K}_{S}={K}_{\mathrm{max}}$, and a set of associated eigenvectors ${\vec{A}}_{1},\ldots ,{\vec{A}}_{S}$ that correspond to decreasing-energy solutions at finite interactions. In order to completely characterize the symmetry of the various states $({K}_{{\ell }},{\vec{A}}_{{\ell }})$ of the degenerate manifold at $g=\infty$, we determine to which irreducible representation of S_N, i.e. to which Young tableau, the solution ${\vec{A}}_{{\ell }}$ corresponds. The Young tableaux are defined in the standard fashion [37]: elements belonging to the same line (column) are symmetric (antisymmetric) under exchange. Thus, for example, for six particles divided in six components, a completely antisymmetric state ${\vec{A}}_{{\ell }}$ will correspond to whereas a completely symmetric state will correspond to . Note that these two states are associated with a single-component non-interacting Fermi gas and a single component Tonks gas, respectively.

In order to classify those states according to their symmetry, we use the k-cycle class sums operators ${{\rm{\Gamma }}}^{(k)}$ [20, 21], defined by ${{\rm{\Gamma }}}^{(k)}={\sum }_{{i}_{1}\lt \cdots \lt {i}_{k}}({i}_{1}...{i}_{k})$, where $({i}_{1}...{i}_{k})$ is the cyclic permutation of particles ${i}_{1},\ldots ,{i}_{k}$. In our system, on the basis of the coordinate sectors, ${{\rm{\Gamma }}}^{(k)}$ is a $N!\times N!$ operator whose elements are ${{\rm{\Gamma }}}_{{ij}}^{(k)}={(-1)}^{(k-1)}$ if going from sector i to sector j only exchange k Fermions in a cyclic fashion, and ${{\rm{\Gamma }}}_{{ij}}^{(k)}=0$ otherwise. Note that, also in this case, we can reduce the dimension of such operator to an S × S matrix, by summing over the contributions of the sectors associated with one snippet.

The spectral decomposition of ${{\rm{\Gamma }}}^{(k)}$ allows to associate each Young tableau Y with a given eigenvalue ${\gamma }_{k}$ [20, 21, 24]. In the following, we shall use in particular the class-sum operator ${{\rm{\Gamma }}}^{(2)}$. Our method consists first in computing and diagonalizing the transposition class sum ${{\rm{\Gamma }}}^{(2)}$ for a given system. Its eigenvalues ${\gamma }_{2}$ can be linked to the Young tableaux according to the expression

$$\begin{eqnarray}&&{\gamma }_{2}=\displaystyle \frac{1}{2}\displaystyle \sum _{i}{\lambda }_{i}({\lambda }_{i}-2i+1),\end{eqnarray} \tag{ 14 }$$

where i and ${\lambda }_{i}$ refer respectively to the line and number of boxes in this line of Young tableau. Thus, projecting a given solution ${\vec{A}}_{{\ell }}$ over the eigenbasis of ${{\rm{\Gamma }}}^{(2)}$ allows to characterize its symmetry and to analyze it in terms of Young tableaux. In table 1 we summarize some of our results for the ground states $({\vec{A}}_{{\rm{max}}},{K}_{{\rm{max}}})$ and the first excited states with a different symmetry $({\vec{A}}_{{\rm{excited}}},{K}_{{\rm{excited}}})$ of different systems for N = 6 particles.

Table 1.  Symmetries and corresponding K eigenvalue (see equation (7)) for the multi-component Fermionic mixtures of figures 1 and 2.

System	${A}_{{\rm{max}}}$	${K}_{{\rm{max}}}$	${A}_{{\rm{excited}}}$	${K}_{{\rm{excited}}}$
$r=2,{N}_{1}={N}_{2}=3$	${Y}_{-3}$	24.97	${Y}_{-5}$	18.91
$r=3,{N}_{1}={N}_{2}={N}_{3}=2$	Y3	30.63	Y0	28.96
r = 6, ${N}_{1}\quad =\cdots =\quad {N}_{6}=1$	Y15	34.33	Y9	33.35
$r=2,{N}_{1}=5,{N}_{2}=1$	${Y}_{-9}$	14.60	${Y}_{-15}$	0
$r=2,{N}_{1}=4,{N}_{2}=2$	${Y}_{-5}$	22.70	${Y}_{-9}$	14.60
$r=3,{N}_{1}=3,{N}_{2}=2,{N}_{3}=1$	Y0	28.96	${Y}_{-3}$	24.97

We set [37]

and the ${Y}_{-\gamma }$ are obtained by taking the symmetric of ${Y}_{\gamma }$ with respect to the main diagonal.

Table 1 shows that the ground and excited states constructed with equation (3) have a definite symmetry which can be readily extracted from the associated Young tableau. An exception is provided by the ground state of the case ${N}_{1}={N}_{2}=3$ as well as the excited state of ${N}_{1}=3,{N}_{2}=2,{N}_{3}=1$, where the transposition class-sum operator ${{\rm{\Gamma }}}^{(2)}$ does not allow to uniquely associate a Young tableau to the wavefunction, since equation (14) gives two Young tableaux corresponding to the eigenvalue −3, and , and also to the eigenvalue 3, and . In this case, the ambiguity is lift off with the help of the 3-cycle class sum ${{\rm{\Gamma }}}^{(3)}$, since these two Young tableaux correspond to two different eigenvalues of ${{\rm{\Gamma }}}^{(3)}$ [20].

Furthermore, for all the mixture considered we note that the ground state, corresponding to $K={K}_{\mathrm{max}}$, is the most symmetric one [38, 39], and the excited state is obtained by decreasing the symmetry of the state by taking out one cell from the top row and putting it in the first available lower row of the Young tableau, thus making it more antisymmetric. The opposite case of K = 0 is associated with the most antisymmetric Young tableau. The above result for the ground state supports the observation of [40], where an ansatz for the ground state wavefunction was suggested, and is also in agreement with a general demonstration provided in [41]. We notice also that when the number of particles coincides with the number of components, i.e. r = N, the ground state is fully symmetric, and has the same symmetry as the bosonic Tonks–Girardeau gas.

This analysis provides a verification and an example of the generalization of the Lieb–Mattis theorem [22] for the case of multicomponent Fermionic mixtures. The theorem states that for N electrons in one-dimension, interacting by a symmetric potential, the energy of a state with total spin $S(S^{\prime} )$ is such that $E(S)\leqslant E(S^{\prime} )$ if $S\lt S^{\prime}$. It follows that the ground state has the smallest possible value for the spin S, which is realized by an antisymmetric spinor. Hence, the spatial wavefunction has the most symmetric configuration. In the case of more than two spin components, as also discussed in [38], the same feature occurs, which is displayed by our results. Since K tends to be maximized when the wavefunction is more symmetric, we can see K as an energetic indicator of the symmetry.

In this work we have considered a multicomponent strongly correlated Fermionic mixture with up to six components. Using a generalization of the pioneering solution due to Girardeau for the Bose gas and of the recently developed solution for the two-component Fermi gas we have determined the exact many-body wavefunctions for the degenerate manifold at infinite interactions. We have identified the one which corresponds to the ground state at finite interactions as the one which has the maximum value for the parameter K, related to the interaction energy. We have then obtained the density profiles for the mixture under harmonic confinement. For an imbalanced mixture we have found a partial phase separation among the components, which is an effect of the strong repulsive interactions. Furthermore, we have characterized the symmetry properties of the ground and of some excited state wavefunctions of the manifold by introducing suitable class-sum operators, and we have shown that the ground-state wavefunctions have a definite symmetry, corresponding to the most symmetric (or less antisymmetric) configuration compatible with the imbalance among the components. Our exact solution for the inhomogeneous multicomponent mixture provides an important benchmark for numerical simulations of strongly correlated multicomponent Fermi gases, in a regime where the presence of the quasi-degenerate manifold challenges the convergence of the calculations, as well as for quantum simulators.

This paper is dedicated to the memory of Marvin Girardeau, of whom AM recalls his enthusiastic guidance. We thank M Gattobigio, J Levinsen, P Massignan, M Parish and T Roscilde for fruitful discussions. AM acknowledges the hospitality of the Centre de Recherches Mathématiques de Montréal where part of this work was performed and financial support from ANR projects Mathostaq (ANR-13-JS01-0005-01) and SuperRing (ANR-15-CE30-0012-02).