Strongly correlated one-dimensional Bose-Fermi quantum mixtures: symmetry and correlations

We consider multi-component quantum mixtures (bosonic, fermionic, or mixed) with strongly repulsive contact interactions in a one-dimensional harmonic trap. In the limit of infinitely strong repulsion and zero temperature, using the class-sum method, we study the symmetries of the spatial wave function of the mixture. We find that the ground state of the system has the most symmetric spatial wave function allowed by the type of mixture. This provides an example of the generalized Lieb-Mattis theorem. Furthermore, we show that the symmetry properties of the mixture are embedded in the large-momentum tails of the momentum distribution, which we evaluate both at infinite repulsion by an exact solution and at finite interactions using a numerical DMRG approach. This implies that an experimental measurement of the Tan's contact would allow to unambiguously determine the symmetry of any kind of multi-component mixture.


Introduction
Ultracold atom experiments allow to engineer and probe, with an incredible and alwaysimproving precision, a yet inaccessible variety of many-body quantum systems [1]. One-dimensional (1D) models, which display several unique features associated to the reduced dimensionality [2], are the object of intense theoretical and experimental interest. Quantum gases in one dimension can be realized in actual experiments by trapping atoms in tight optical waveguides [3][4][5][6]. Moreover, these experiments offer the possibility to tune the interactions between the atoms [1,7], and hence to access the strongly correlated regime for various kinds of quantum mixtures. For instance, a strongly repulsive two-component 1D Fermi gas was realized in [8], and a mixture with up to six fermionic components was realized using Ytterbium atoms [6]. Other experiments with Ytterbium isotopes offer promising perspectives for the realization of strongly interacting 1D Bose-Fermi mixtures [9].
In multi-component Bose-Fermi quantum mixtures, the exchange symmetry between identical particles is fixed by their bosonic or fermionic nature, but not between distinguishable particles. Thus, a natural question is to find how to characterize, both theoretically and experimentally, the global exchange symmetry of the many-body wavefunction. This fundamental property is directly related to the magnetic properties of the system [10][11][12]. Therefore, strongly interacting 1D atomic mixtures appear to be a perfect experimental and theoretical playground for studying quantum magnetism [13]. They also offer the opportunity to study itinerant magnetism phases, i.e. magnetism without a lattice [14].
A typical feature in strongly interacting gases is the so-called fermionization: at increasing repulsive interactions, particles which are not subjected to the Pauli principle tend to avoid each other. In the limit of infinite repulsion, contact interactions mimic the Pauli exclusion principle, e.g. they induce zeros in the many-body wave function when two particles meet. This phenomenon allows the construction of an exact solution for any type of quantum mixture through a mapping onto a non-interacting spinless Fermi gas [11,[15][16][17][18][19].
In this work, we focus on 1D multi-component Bose-Fermi quantum mixtures with strongly repulsive contact interactions, taking the experimentally relevant case of a harmonic external confinement. In the limit of infinite interactions, we provide an exact solution of the multi-component mixture with arbitrary number of components using the method developed in [11,19]. In order to address the symmetry of the many-body wave function of the ground and excited states we use the class-sum method. We show that the mixtures obey a generalized version of the Lieb-Mattis theorem [20]. Namely, as obtained for multi-component fermions in [21], the spatial ground state wave function of each system is the most symmetric one. Our analysis of the symmetries also confirms the ordering of the excited states recently obtained in [22] for fermionic mixtures and extends this concept to the case of Bose-Fermi mixtures. Furthermore, we explore how symmetry affects the properties of the system in coordinate and momentum space by studying density profiles and momentum distributions. For this purpose we combine the exact solution at infinite interactions with numerical Density Matrix Renormalization Group (DMRG) calculations at finite interactions. This allow us to analyze the symmetry structure of the mixture in the fermionization process. In particular, we show that a precise measurement of the asymptotic momentum distribution would allow to probe the symmetry of the mixture, thus generalizing our previous result on multi-component fermions [23].
The paper is organized as follows: In Section 2, we give a description of our model and of its solution. Then, in Section 3, we provide a detailed description of the classsum method that we use to extract the exchange symmetry of the many-body wave function. In Section 4, we provide our results for the density profiles and momentum distributions. Finally, Sec. 5 gives our concluding remarks and perspectives.

Model and general solution
We consider a mixture of and f flavors (e.g. spin components), all having the same mass m. For a Bose-Fermi mixture this is always an approximation, which is appropriate e.g. for mixtures of isotopes. This type of mixtures is experimentally accessible with Ytterbium atoms [9] for instance. The particles are confined in a tight atomic waveguide so that their motion can be considered one-dimensional. Their positions are given by the coordinates For the sake of simplicity of notations, the exponents specifying the type and flavor of the particle will be from now on omitted. All the particles are also subjected to the same one-dimensional confinement V (x) = 1 2 mω 2 x 2 , which could be realized e.g. in cigar-shaped optical traps. The interactions between the particles are modelled via the pseudo-potential V int (x, y) = gδ(x − y). Here g = −2 2 /ma 1D , and a 1D is the 1D effective scattering length [7], which accurately describes the s-wave scattering dominated collisions in a restricted geometry at low enough energies and densities [1]. We assume that the bosonboson, boson-fermion and fermion-fermion interactions are characterized by the same interaction strength g. In the fermionic case, no s-wave collisions occur among fermions belonging to the same flavor. The Hamiltonian describing the system is then given by The interaction term in H is equivalent to imposing the cusp condition on the many-body wave function Ψ = Ψ(x 1 , . . . , x N ),

Exact many-body wave function at infinite repulsion
As a consequence of the cusp condition, in the strongly correlated limit g → ∞, the many-body wave function vanishes whenever x i = x j , leading to fermionization features. In this limit one can build an exact many-body wave function satisfying the cusp condition by taking [11,16] Ψ(x 1 , . . . , where S N is the permutation group of N elements, θ P (x 1 , . . . , x N ) is equal to 1 if x P (1) < · · · < x P (N ) (coordinate sector, later indicated as (P (1), P (2), ...P (N ))) and 0 otherwise, and Ψ A is the fully antisymmetric fermionic wave function Ψ A (x 1 , . . . , ..,N where φ 0 , . . . , φ N −1 are the eigenfunctions of the single particle Hamiltonian The coefficients a P in Eq.(3) corresponding to exchange of particles belonging to the same bosonic (fermionic) component of the mixture are such that a P = 1 (a P = −1) respectively. Hence, the number of independent coefficients is reduced to D N,b+f = N ! N 1 !...N b+f ! . This allows to extremely reduce the cost of calculations, by identifying sectors that are equal modulo permutations of identical particles. This smaller basis is often called the snippet basis [16,18]. In this paper, we will switch from one formalism to the other depending on whether a distinction between identical bosons is necessary or not.
In order to find the coefficients a P we use the variational approach developed in [11]. This is based on a strong-coupling expansion of the energy to order 1/g, i.e. by setting E = E A − K/g where E A is the fermionic energy associated with Ψ A . Using the Hellmann-Feynman theorem together with (2) we obtain where α P,Q = dx 1 . . . dx N θ Id (x 1 , . . . , x N )δ(x k − x k+1 ) [∂Ψ A /∂x k ] 2 ≡ α k if P and Q are equal up to a transposition τ k of two consecutive particles that are either distinguishable particles or indistinguishable bosons at positions k and k + 1, and 0 otherwise (see [10,21] for computational methods). Note that the spatial parity of the trap implies that α k = α N −k for all k. Then, the variational condition ∂K/∂a P = 0 is shown to be equivalent to a simple diagonalization problem V a snip = Ka snip , where a snip is the vector of the D N,b+f independent a i coefficients and V is a D N,b+f × D N,b+f matrix defined in the snippet basis by where the index d means that the sum has to be taken over snippets k that transpose distinguishable particles as compared to snippet i, while b means that the sum is taken over sectors that transpose identical bosons. The highest eigenvalue of V yields the ground state of the system, whereas the other eigenvalues correspond to all the excited states that belong to the same D N,b+f -degenerate manifold at g → ∞.
Exchange symmetry Young diagram(s) symmetric antisymmetric mixed Mixture Young tableaux Table 1: Young diagrams and standard Young tableaux for different two-component mixtures with N = 4. The particles of the species are labeled a resp. b, and can be either bosonic or fermionic.

Symmetry of the mixture
In this section, we study the spatial symmetry of our solution under exchange of particles. To this purpose we use the class-sum method, inspired from nuclear physics [24,25], and previously used in [18,21,23].

General method
We consider a N -particle wave function Ψ(x 1 , . . . , x N ), identified by a vector a sect = (a 1 , . . . , a N ! ) as in Eq. (3). We want to determine its permutational symmetry, i.e. find to which irreducible representation(s) of the permutation group S N it belongs. Each irreducible representation can be labeled by a Young diagram, where, according to the standard notation, boxes in the same column (resp. row) correspond to an antisymmetric (resp. symmetric) exchange of particles [26]. Examples of Young diagrams and their corresponding symmetries for N = 4 are given in Table 1. Using the Young diagrams, it is possible to know a priori, given a configuration N B 1 + . . .
which irreducible representations are possible for Ψ, compatible with the antisymmetrization of fermions belonging to the same component, and symmetrization for same-component bosons. More precisely, in order to build the standard Young tableaux for a given mixture, we index by a letter the particles of each given species. We therefore order the species by decreasing order and label them alphabetically. We then label the boxes of the Young diagrams imposing that the entries of each row and column are always increasing (in weak sense), and with the symmetry constraint that no more than one same-component bosonic/fermionic label can be present in one column/row. Examples are provided for 4-particle mixtures in Table 1. We also recall that the number of standard Young tableaux corresponding to a given Young diagram is equal to the dimension of the associated irreducible representation [27]. The number d λ of possible standard Young tableaux associated with a Young diagram λ is given by the Hook-length formula where h λ (i, j) is equal to the number of cells below the box (i, j) + the number of cells at the right of the box (i, j) + 1. Interestingly, the dimension d λ counts the number of energy eigenstates having the symmetry λ [28,29]. In order to obtain the symmetry associated with a given wave function Ψ(x 1 , . . . , x N ) belonging to the degenerate manifold, we define a set of N ! × N ! matrices whose eigenvalues are directly connected to the irreducible representations of S N , namely the conjugacy class-sums [27,30]. Considering cycles σ of length |σ| = p, which permute p elements of {1, . . . , N } in a cyclic way, the p-cycle class-sum Γ (p) is represented in the coordinate-sector basis of N non-interacting fermions as M σ is an N ! × N ! matrix whose elements are (M σ ) P Q = (−1) |σ|−1 δ P,σ•Q , and the factor (−1) |σ|−1 is due to the inherent antisymmetry of the basis. As an example, for N = 3, taking as basis ((1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1)), one has and The key point of our method is that the class-sum eigenvalues γ (p) are directly related to the irreducible representations of S N : given a Young diagram λ = [λ 1 , . . . , λ n ], where λ i is the number of boxes in row i, the following relation holds [31]: where µ i = λ i − i + n. ‡ The fundamental reason behind the above connection (10) between class sums and Young tableaux is that there is a one-to-one correspondence of the eigenvalues of a class-sum and the characters of the irreducible representations of S N [25,30]. Once one considers the relation between the irreducible representations of S N ‡ Note that there is a notation misprint in Ref. [31] where n refers to the total number of atoms instead of the total number of rows. and the ones of SU (κ) (for the case of κ-species mixtures) [26], it is possible to connect the eigenvalues of the class-sum operators and those of the SU (κ) Casimir operators. The latter are the ones that, roughly speaking, describe the character of SU (κ) irreps. and thus generalise the SU (2) concept of spin length.
Using the above properties, in order to determine the symmetry of the wave function Ψ of a given mixture we proceed as follows: first, we determine the matrix Γ (2) on the coordinate-sector basis and compute its eigenvalues γ (2) . Then, we list all the Young diagrams allowed within the given type of mixture, and using Eq. (10), we associate each eigenvalue γ (2) to a given Young diagram. If different diagrams correspond to the same γ (2) , we repeat this operation to higher order in the size p of the cycles until there are no longer degeneracies. Finally, we expand the wave function Ψ, represented in the same basis as the vector a sect = (a 1 , . . . , a N ! ), on the basis of the eigenvectors of Γ (p) in order to determine the weights of the different irreducible representations, and hence the symmetry content of the wave function.
Finally, we would like to comment that the class-sum method has been here adapted to the ansatz (3), but it can be used as well for other ansatz, e.g. the Bethe ansatz [12,32], as long as the M σ matrices in Eq. (7) are defined correspondingly.

Results and discussion
The case of fermionic mixtures has already been studied in [21,23], and the case of simple Bose-Fermi mixtures N B 1 + N F 1 with N B 1 = N F 1 was studied in [18]. In this work, we extend these results to the more complex and general case of Bose-Fermi mixtures with more than one bosonic/fermionic components. In particular, we study the cases 2 B + 2 F , 3 B + 3 F , 2 B + 2 B + 2 F and 2 B + 2 F + 2 F . Note that we only need p max = 3 in order to discriminate between and or between and , which are degenerate in Γ (2) . We focus on the ground states of a given symmetry configuration, in the sense that, given a symmetry S, they have the highest eigenvalue of V denoted K(S), corresponding to the the lowest energy E(S). Our results are shown in Table 2. First, we observe that, as in the case of Fermi mixtures [21], these states have a pure symmetry, i.e. they belong to a unique irreducible representation. Moreover, a direct comparison of Table 2(a) and Table 2(b) shows that Bose-Fermi mixtures with more than two components have much complex symmetric structures than the two-component ones, the latter having always only two possible symmetries [18]. Another interesting observation is that a given symmetry S corresponds to the same ground state energy slope K(S), independently of the mixture. In that sense, the symmetry of the mixture appears to be a deeper physical characterization than the actual choice of the mixture.
Furthermore, it is clear in Table 2 that diagrams with higher K(S) or, equivalently, lower E(S), are more "horizontal", or more "symmetric". In the context of a two-spin component fermionic mixture, this observation is known as the Lieb-Mattis theorem and has very strong implications in the theory of ferromagnetism [20]. It states that Mixture Symmetry γ (2) K(S)/( 2 ω 2 a ho )  the mixture with the lowest total spin, and therefore the most symmetric spatial wavefunction, will display the lowest ground-state energy, which implies that the ground state is unmagnetized. On the other hand, it was proven in [33] that the ground state of an interacting bosonic mixture is always fully polarized, which again is equivalent to the property that the spatial wave function is the most symmetric. Our results show that the ground state corresponds to the most symmetric spatial configuration allowed by the Bose-Fermi mixture, generalizing the results obtained in the case of purely fermionic mixtures in [21,23], in agreement with the second Lieb-Mattis Theorem (LMT II) [20].
The notion that a state S is "more symmetric" than S is usually understood in the sense of the so-called pouring principle, which in terms of Young diagrams means that the diagram corresponding to S can be poured into the diagram corresponding to S by moving boxes up and right [20]: e.g.
can be poured into . The LMT II claims that E(S) > E(S ) if S can be poured in S . Nevertheless the pouring principle has its limitations: clearly and can not be compared according to this principle. Here, we show that these symmetries display different interaction energies and are therefore comparable (e.g. K < K ), confirming the ordering observed in [22] in the case of a Fermi gas and extending the analysis to the case of Bose-Fermi mixtures.

Density profiles and momentum distributions
This section is dedicated to the study of the effect of symmetries on the one-body density matrix ρ ν (x, y) in the strongly interacting regime. For any component ν ∈ 1, . . . , b + f, ρ ν (x, y) is defined by where the first coordinate of Ψ is chosen to belong to species ν. The one-body density matrix measures the first-order spatial coherence among two particles of species ν in positions x and y respectively. From the knowledge of the one-body density matrix one can obtain the density profile, according to This quantity measures the probability of finding a particle of species ν in position x. It is accessible in cold-atom experiments via in situ imaging [1]. The momentum distribution n ν (k) of species ν is also directly extracted from the one-density matrix by Fourier transform, according to The momentum distribution is one of the most common observables in cold-atom experiments. All the components of a mixture can be separately measured using spinresolved time-of-flight techniques [1,34]. The large-momentum tails of the momentum distribution for each component of a mixture follow a universal k −4 law, as first proven for the bosonic Tonks-Girardeau gas in [35,36]. The weights C ν = lim k→∞ k 4 n ν (k), known as Tan's contacts, are of primary importance because of their relations to many other physical quantities, as the two-point correlators or the interaction energy, in various quantum mixtures and trapping potential [23,[37][38][39][40]. One of these relations allows to link C ν to the interaction energy, and more specifically to the part of the interaction energy to which species ν contributes. More precisely, if we denote H int,νν the interaction energy between species ν and ν , the following relation holds

Exact results at infinite interactions
In order to compute the one-body density matrix in the infinitely repulsive case g → ∞, we start from the exact solution for the many-body wave function (3). Few derivation steps allow to reduce the many-body integration in (11) to the calculation of singleparticle integrals: we label the N ! permutations P ∈ S N by an index i k , where i ∈ {1, . . . , N } denotes the position of the first particle, and k ∈ {1, . . . , (N − 1)!} labels the permutations of the N − 1 other particles. Then, choosing x ≤ y, Eq. (11) becomes where with Noticing that dx 2 . . . dx N θ P i k (x, . . . , x N )θ P i l (y, . . . , x N ) ∝ δ kl , we obtain after some algebra [21,41]: (17) where a i k a j k (i−1)!(j−i)!(N −j)! , and where the integration limits (L ij , U ij ) are given by The density profile is readily obtained as n ν (x) = N ν i R (ii) ν (x, x), and the momentum distribution is obtained by numerically evaluating the Fourier transform (13) for the specific case of Eq. (15).
The exact solution allows also for an independent evaluation of the Tan's contact using the Tan relation on the interaction energy. Using Eq. (4), we find that the interaction energy H int,νν is given by where σ N (µ, ν, k) is the set of all permutations such that particles in k and k+1 positions are from the species µ and ν (or are particles i and j for σ N (i, j, k)), and τ k is the transposition of these two particles. Note that H int,νν = 0 if species ν is fermionic. The calculation of the interaction energy allows to obtain the Tan's contact using Eq.(14).

Numerical solution at arbitrary interactions
In realistic experimental conditions, the interactions are tunable, allowing to explore all the interactions regime from small to large values. In particular, a relevant question is to what extent the density profiles and the momentum distributions evaluated from the exact solution at g → ∞ describe accurately the ground states at large but finite interaction. In order to tackle the arbitrary interaction regime, we have performed numerical calculations using a two-site optimization scheme, according to the Density-Matrix Renormalization-Group (DMRG) method, based on matrix product states (MPS). The code allows the exploitation of Abelian symmetries (here particle-number conservation for each species). To address the continuous problem, we discretized the relevant region of the trap and calculated properties of the resulting lattice model. The continuum results were then obtained by increasing the number of lattice points in consecutive simulations (thereby decreasing the lattice spacing) and finite-size scaling the quantities of interest. Density distributions were measured as site-occupations in the discretized model; the Tan contacts were determined through the interaction part of the energy. The momentum distributions were obtained from Fourier-transformed two-point correlators n ν (k) = 1 L L j,l=1 e i(j−l)k c † ν,j c ν,l , where c ( †) ν,j are the annihilation (creation) operators for the species ν in the lattice model. The maximal number of lattice sites considered was 216, covering a region of 12 harmonical trap lengths. The virtual bond-dimension of the MPS was m ≈ 150 which appears of be sufficient for the very low density of particles we consider here.

Density profiles
The density profiles of a multi-component mixture display a wealth of information concerning the intercomponent interactions as well as the magnetic structure.
At mean-field level, by increasing the strength of repulsive intercomponent interactions [42] the system is predicted to undergo spatial separation. However, for strong interactions mean-field theory loses its validity due to the change of equation of state and the increased effect of quantum fluctuations. Luttinger liquid theory for a homogeneous system also predicts an instability towards spatial separation [43]. For a binary mixture in a harmonic trap, an analysis based on the local-density approximation of the exact equation of state was performed in [44,45], yielding a partial spatial separation of the two components at strong interactions. For small particle numbers, and finite interactions, spatial separation was found within some parameter regimes using mapping to a spin model and exact diagonalization in [46], and for finite interactions an exact solution showing spatial separation was obtained in [19].
We present in Figs. 1-2 our results for both binary and ternary multi-component mixtures made of six particles. The results for the binary mixtures agree with [19,46] Figure 1: (Color online.) Exact solution for the bosonic (brown) and fermionic (orange) ground state density profiles (in unit of a −1 ho ) in the limit g → ∞, as functions of the spatial coordinate x/a ho . From left to right, the mixtures are 3 B + 3 F , 2 B + 2 F + 2 F and 2 B + 2 B + 2 F . For the case of ternary mixtures, the density profiles of the two components with the same statistics and same particle number coincide.   configuration allowed by the type of mixture, then it is built in such a way to minimize the number of anti-symmetric exchanges, which is the case if the fermions are on the edges. For ternary mixtures, at fixed numbers of particles in each component, we find that the details of the density profiles depend on the type of mixture. This shows that the spatial separation is not only an effect of the interactions but also of the symmetry. Also, we remark that the results at finite interactions display the same features as the strongly-correlated regime already at intermediate interactions. This suggests that the spatial symmetry is already present for finite interactions.

Momentum distributions and Tan's contact
We now turn to the discussion of our results for the momentum distribution. Fig. 3 shows our predictions for the regime of infinite interaction strength. The momentum distribution of the bosonic components display one central peak, whereas the fermionic components display as many peaks as the number of the fermions in the corresponding species.
This feature, typical of non-interacting spinless fermions in harmonic confinement [47], has also been reported for multi-component interacting fermionic mixtures in harmonic traps [48]. The occurrence of different peaks in the fermionic momentum distribution can be seen as a consequence of the Pauli principle for fermions belonging to the same species. Fig. 4 shows the dependence of the momentum distributions on interaction strength. At increasing interactions we notice that the number of peaks in each distribution is conserved, and the shape of both fermionic and bosonic momentum distributions looks qualitatively the same. However, a detailed analysis shows that both bosonic and fermionic momentum distributions shrink in their central part. This is complementary to the coordinate space description, where repulsive interactions tend to broaden the profiles. Furthermore, the tails of the momentum distribution become more and more important at increasing interactions, as it is visible on Fig. 4. These tails are shown in Fig. 5 for a finite value of the interaction parameter g = 10.0 ωa ho and in Fig. 6 for the limiting case g → ∞. All the mixtures analysed exhibit the k −4 behavior with the expected Tan's contact coefficients. The asymptotic behaviour of the tails (dashed lines) has been obtained from Eq. (14), where the interaction energy for the mixtures at g = 10.0 ωa ho has been numerically calculated (Fig. 5), while for the mixtures at g → ∞ has been evaluated from the exact solution via Eq. (19) (Fig. 6). For balanced mixtures with the same number of particles in each component we notice that Tan's contact are not equal for bosonic and fermionic species, at difference from the case of multi-component fermionic mixtures [23]. This properties readily follows from Eq. (19), showing that in the bosonic case there is an additional term corresponding to intracomponent boson-boson interactions.
Finally, we notice that the knowledge of the Tan's contacts C ν for each component of the mixture allows to deduce the symmetry of the mixture, since the total contact C tot = b+f ν=1 C ν is linked to the energy slope K via Eq. (4) and (14) through and each K(S) corresponds to a unique symmetry S (see Table 2 in Sec. 3). Thus, an experimental measurement of Tan's contact, e.g. through spin-resolved time-of-flight measurements, or through the energy, would allow a symmetry spectroscopy of the Bose-Fermi mixture [49]. This generalizes our previous results on fermionic mixtures [23].  (14). The bending of the momentum curves for very large momenta is an artifact originating from the discretization procedure. Figure 6: (Color online) Bosonic (brown) and fermionic (orange) ground state k 4 n ν (k) functions in unit of a −3 ho and as functions of ka ho . From left to right, the mixtures are 3 B + 3 F , 2 B + 2 F + 2 F and 2 B + 2 B + 2 F . The horizontal lines are obtained from the exact solution via Eq. (19) in Eq. (14).

Conclusions
In this work we have studied the spatial and the momentum distributions of a zerotemperature multi-component boson-fermion mixture, trapped in a one-dimensional harmonic potential. We have considered the ideal case of equal-mass particles and equal contact-interaction strength, that could be obtained in a good approximation in the experiments by trapping and cooling large atomic-number isotopes (e.g. Ytterbium). At infinitely large repulsive interactions we have determined the ground-state properties by using an exact many-body wavefunction. At finite interactions, we have calculated the spatial density profiles and momentum distributions by numerical DMRG methods, obtaining an excellent agreement between the exact calculations and the numerical ones. We have found that at large interactions bosons and fermions are spatially phase separated in the trap and that the momentum distribution tails -the contacts -at fixed number of particles, increase by increasing the number of bosons and/or the number of fermionic components. Both effects are related to the wavefunction symmetry, that we have shown to be, by exploiting a sum-class method, the most possible symmetric one for each type of mixture. Our results are in agreement with a generalized version of the Lieb-Mattis theorem. Importantly, we have found that the exchange symmetry of the many-body wave-function is preserved for large but finite interactions, as in typical experimental conditions, and we have shown that this symmetry is measurable through the tails of the momentum distribution.