Superfluid instability in ultra-cold gas of fermionic atoms with attractive potential in a one-dimensional trap

In the context of ultracold fermionic atoms with effective spin S confined to an elongated trap we study the one-dimensional gas interacting via an attractive δ-function potential using the Bethe ansatz solution. There are N = 2S + 1 fundamental states: The particles can either be unpaired or clustered in bound states of 2, 3, ..., 2S and 2S + 1 atoms. In a magnetic field, the rich ground state phase diagram consists of these N states and various mixed phases in which combinations of the fundamental states coexist. The phase diagram simplifies considerably in zero-field, where only bound states of N atoms can exist. Due to the harmonic confinement and within the local density approximation, the density profile of bound states decreases along the tube from the center of the trap to its boundaries. In an array of tubes with weak Josephson tunneling superfluid order may arise. In zero-field the response functions determining the superfluid and density wave order are calculated using conformal field theory and the exact Bethe ansatz solution. The response function for superfluidity consists of a power law with distance, while the correlation function for density waves is a power law of distance times a sinusoidal factor oscillating with distance with a period given by two times the Fermi momentum. For S = 1/2 superfluidity is a possibility for all densities and density waves can be excluded. For S ≥ 3/2 superfluidity may occur at low densities but at high densities it gives way to density waves. We discuss the scenario of phase separation where for S ≥ 3/2 the system has superfluid long-range order toward the trap boundaries and density waves at its center.


Introduction
Spin-imbalanced ultracold gases of fermionic atoms confined to one-dimensional traps have been the subject of several recent studies [1,2,3]. Confinement to nearly one-dimensional tubes can be achieved if the ultracold cloud of atoms is subjected to a two-dimensional optical lattice, which defines a two-dimensional array of tubes [1]. The tubes can be regarded as isolated if the confinement by the laser beams is strong enough to suppress tunneling between tubes. The scattering between atoms under transverse harmonic confinement is subject to a confinementinduced resonance [4]. Fine-tuning this Feshbach-type resonance, the interaction between particles can be made attractive and its strength can be varied [5]. The interaction is local and can be approximated by a δ-function potential in space. The weak confinement along the tube is roughly harmonic and can be locally incorporated into the chemical potential. Consequently, these systems are only locally homogeneous and within the local density approximation display phase separation with the variation of the chemical potential along the tube [2,3].
For fermions with a spin S an attractive potential may lead to bound states of up to N = 2S+1 atoms, extending the concept of preformed Cooper pairs to larger clusters. The ground state phase diagram will have large number of possible pure (consisting of only one kind of clusters) and mixed phases (coexistence of two or more basic states) [6,7]. Interesting physics arises from the transitions between these phases as a function of the chemical potential and the external magnetic field. In this paper we consider an attractive δ-function potential between the atoms of spin S and assume the model has SU(N ) symmetry, i.e. the scattering is independent of the spin components.
One-dimensional spinless bosons with δ-function interaction were first studied by Lieb and Liniger [8] using Bethe's ansatz. These results where extended to spin-1/2 gases by M. Gaudin [9] and C.N. Yang [10]. For an attractive interaction the S = 1/2 ground state has two classes of solutions of the discrete Bethe ansatz equations, namely, real charge rapidities and paired complex conjugated rapidities [9,11,12], representing spin polarized particles and bound states of the Cooper type, respectively. The preformed Cooper pairs are gapped in a magnetic field smaller than a critical one (energy required to break up the pairs) and do not display longrange order. Similar results were obtained for the Hubbard model with attractive U [13,14]. Sutherland [15] generalized C.N. Yang's solution for spin 1/2 [10] to an arbitrary number of colors N = 2S + 1 [SU(N )-symmetry]. Takahashi [16] derived the integral equations for the ground state density functions for bound states up to N = 2S + 1 particles. The classification of states, the thermodynamics, the ground state equations and elementary excitations of the gas with arbitrary number of colors have been derived by Schlottmann [6,17,18] for both attractive and repulsive potentials. Several of these results have been recently rederived in the context of ultracold fermion gases [19,20,21,22].
For an isolated tube there is no long-range order of the bound states. A weakly interacting array of tubes, e.g. via tunneling of bound states between tubes, however, increases the effective dimension of the system so that long-range order can arise [23,24]. For a spin-imbalanced system this could lead to the realization of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) [25] phase for superfluidity in an ultracold gas of atoms [1,26]. The system remains strongly anisotropic and pure (there are almost no impurities) and is hence favorable for inhomogeneities like modulations of the order parameter of the FFLO type in the presence of an external magnetic field. As for unconventional superconductors, the instability to superconductivity is determined by the properties of the normal phase. In this case the normal phase is a Luttinger liquid. Evaluating the superfluid response functions for the Luttinger liquid, one can then determine which one is the first superfluid instability reached from the normal phase. There is a large number of possible superfluid order parameters. The correlation functions have been studied in detail for S = 3/2 and S = 5/2 using conformal field theory and the Bethe ansatz solution [27,28]. It was concluded that in a strong magnetic field the dominant phase is the pairing of atoms with spin components S z = S and S z = S − 1, one atom being a forward mover and the other one a backward mover. Since the Fermi gas is spin-imbalanced in a magnetic field, the Fermi momenta of the different spin components are different yielding sinusoidal oscillations in space, characteristic of an FFLO state. In the one-dimensional case the periods of oscillation can be extracted from the corresponding correlation functions. The correlation functions are the product of a power-law dependence of the distance and a cosine term with the desired periodicity.
In this paper we investigate the possibility of superfluidity in zero magnetic field. In zero-field only bound states of N atoms, each with a different spin-component, exist. Using the Bethe ansatz equations we show that the critical exponent of the superfluid response function increases with the particle density. The superfluid order parameter will be homogeneous and not of the FFLO type. On the other hand, the exponent of the particle density wave response function decreases with the particle density and for S ≥ 3/2 density waves are the dominating form order at high densities. As a consequence of the longitudinal confining harmonic potential, phase separation is expected, with the trap boundaries having a higher chance of being superfluid for S ≥ 3/2 than the center of the trap, which is more prone to density waves.
The remainder of the paper is organized as follows. In Sect. II we restate the model, the classification of states of the Bethe ansatz equations and the integral equations satisfied by the dressed energies and the densities for bound states of N particles in zero magnetic field. We also discuss the density profile due to the varying chemical potential along the trap in zero magnetic field. In Sect. III we study the response functions for superfluidity and density waves in zero magnetic field using the conformal invariance and the Bethe ansatz solution of the system. Conclusions are presented in Sect. IV.

The model
The Hamiltonian for a gas of nonrelativistic particles with (2S + 1) colors (spin S) interacting via an attractive δ-function potential is where x i are the coordinates, N t is the total number of particles and c is the interaction strength. By fine-tuning the confinement-induced resonance [4] the interaction can become attractive and its strength can be varied. Hereh 2 /2m, where m is the mass of the particles, has been equated to 1, or alternatively it has been scaled into H and c.

Classification of Bethe states
The states of the coordinate Bethe ansatz are plane waves constructed from the two-particle scattering matrix. This scattering matrix satisfies the so-called Yang-Baxter triangular relation, which is a necessary condition for integrability. As a consequence of the triangular relation many-particle scattering processes can be factorized into two-particle processes and the order in which the individual scattering processes take place can be interchanged (the order becomes arbitrary). A solution can be constructed by iteratively applying (nesting) N = 2S + 1 Bethe ansätze, such that one color is eliminated at each step [15]. Each Bethe ansatz gives rise to a new set of rapidities, {k j }, j = 1, · · ·, N t for the charges (coordinate Bethe ansatz) and {Λ (l) α }, l = 1, · · ·, N − 1, with α = 1, · · ·, M (l) for the internal degrees of freedom (spin). Here M (l) is the number of rapidities in the l th set and α is the running index within each set. All rapidities within a given set have to be different to ensure linearly independent solutions.
For an attractive interaction and large length of the trap, L, the solutions of the discrete Bethe equations can be classified according to (i) real charge rapidities, belonging to the set {k j }, associated with unpaired propagating spin-polarized particles, (ii) complex spin and charge rapidities, which correspond to bound states of particles with different spin components, and (iii) strings of complex spin rapidities, which represent bound spin states [17,18]. States in class (iii) are not represented in the ground state; these states correspond to excited states and are not considered here. The classification of states is completely determined by the SU(N )-invariance of the model and the attractive effective potential. The classification of states is therefore completely analogous to that of the Anderson impurity of arbitrary spin in the U → ∞ limit [29,30,31,32,33] and the one-dimensional degenerate supersymmetric t − J model [34,35]. The interaction strength in the case of the Anderson impurity model is determined by the hybridization, while the supersymmetric t − J model has no energy scale and the interaction strength is equal to one. Since only particles with different spin components experience an effective attractive interaction, we may build bound states of up to (2S + 1) particles. A bound state of n particles (n ≤ N = 2S + 1) is characterized by one real ξ (n−1) rapidity and in general complex Λ (l) rapidities, l < n − 1, given by [16] These spin and charge strings form classes (i) and (ii), which are present in the ground state [16,17]. The real rapidities ξ (n−1) have all to be different and satisfy the Fermi-Dirac statistics, i.e. the states are either occupied or empty. In the ground state the rapidities are densely distributed in the interval [−B l , B l ] and we denote with ε (l) (ξ), l = 0, 1, · · ·, 2S, the dressed energy potentials (entering the Fermi-Dirac distribution).

Dressed energy potential in zero magnetic field
In zero magnetic field and in the thermodynamic limit the number of particles of each spincomponent (color) is the same, i.e. there is no spin imbalance. It is then energetically favorable to bind all the particles into clusters of N atoms. Each cluster has one atom of each color. These clusters are the generalization of preformed Cooper pairs to N colors. Hence, B l = 0 for l < 2S and only B 2S ≥ 0. Only one energy potential ε (2S) (ξ) needs to be considered, which satisfies the integral equation [6,18] where the D 2S (ξ) is the driving term given by [6,18] and K 2S,2S (ξ) is the integration kernel, which can be written in a compact form [6] Here μ is the chemical potential, playing the role of the Lagrange parameter for the conservation of the total number of particles. μ determines the integration limit B 2S through the condition that ε (2S) (±B 2S ) = 0. Here B 2S determines the two Fermi points of the system, since occupied states correspond to ε (2S) (ξ) < 0 and for empty states ε (2S) (ξ) is positive. Since bound states of N particles have one particle with each spin component, these clusters do not have magnetic response and Eq. (3) does not depend on the magnetic field (assumed to be zero or less than the binding energy of the bound states). If ε (2S) is rescaled to ε (2S) /c 2 , μ to μ/c 2 , B 2S to B 2S /|c| and ξ to ξ/|c|, the equation (3) is universal, i.e., independent of the magnitude of |c|. Hence, within the framework of the grand canonical ensemble, without loss of generality, it is sufficient to present the results for |c| = 1.

Distribution density of rapidities
The density function of the rapidities is obtained from the dressed energy ε (2S) (ξ) by differentiation with respect to μ, i.e., [6,18] where ρ (2S) (ξ) is the particle density and ρ (2S) h (ξ) the corresponding hole density for bound states involving 2S + 1 particles. The integral equation satisfied by the density function is where the integration kernel and the integration limits are the same ones as for the dressed energy potential, but the driving term is different [16]. The number of bound states per unit length is given by If ξ and B 2S are rescaled to ξ/|c| and B 2S /|c|, respectively, as discussed above, then ρ (2S) remains invariant and n 2S rescales to |c|n 2S . The density of the total number of particles is given by N t /L = (2S + 1)n 2S and the ground state energy of the system by where ∞ is the energy density. If the band is partially filled, the group velocity for bound states of N particles is defined as [6,18] The low lying excitation states form a simple Dirac sea with two Fermi points at ξ = ±B 2S . The Fermi momentum is p F,2S = πn 2S , i.e., it is completely determined by the occupation of the band.

Local density approximation
An ultra-cold atom system is inherently inhomogeneous since the diameter of the tube gradually changes with position from the center of the trap to its boundaries. As a consequence of the changing diameter of the tube, the quantization in the plane transversal to the tube gradually changes the zero of energy. This change can be represented by a harmonic potential of frequency ω ho , so that the actual local chemical potential μ is a function of x given by Within the local density approximation, it is μ(x) that enters the Bethe equations, (3)-(5). The Bethe ansatz solution is then exact for the one-dimensional system, but approximate for the trap. This approximation [1,2,3] is expected to be good if the variation of μ with x is slow, i.e. it is the largest length scale in the system other than the length of the trap. The approximation neglects the quantization of the harmonic confinement, which is treated classically and is locally incorporated into the chemical potential. Assuming the trap is symmetric and has length L, then given μ(0) and μ(L/2), i.e. the chemical potential at the center and boundary of the trap, the position along the trap can be written as (from Eq. (11)) This parametrization eliminates the constants m and ω ho from Eq. (11), and the local density profile can be expressed as a function of the position along the trap, x. The occupation of the band of bound states of N particles grows monotonically with the chemical potential. For μ ≤ −S(S + 1)c 2 /3 the band is empty, i.e. there are no particles in the system. The integration limit B 2S grows with μ, being zero for μ 0 = −S(S + 1)c 2 /3. The variation of μ and (2S + 1)n 2S as a function of B 2S is shown in Fig. 1 for various S. The quantities have been scaled with |c| so that the results are universal for all interaction strengths. Note that μ − μ 0 grows approximately quadratically with B 2S as a consequence of the parabolic dispersion of the Hamiltonian. For small B 2S the dependence of n 2S on B 2S is roughly linear.
In Fig. 2 we show the density profile of bound states of N atoms, n 2S , as a function of the position along the trap x, S = 1/2, · · ·, 7/2, μ(L/2) = −S(S + 1)c 2 /3 and μ(0) = c 2 . Note that for the parameters chosen the density of particles decreases rapidly when the boundary of the trap is approached.

Conformal towers and correlation functions
In an array of tubes (optical lattice) interactions between particles in different tubes and/or Josephson tunneling between tubes may give rise to a dimensional crossover from one dimension to three dimensions and opens the possibility for superfluid long-range order [23,26]. As for high temperature superconductors, where the instability to superconductivity and properties of the superconducting phase are determined by the properties of the normal phase, in the present case superfluidity is determined by the Luttinger liquid. To investigate the possibility of superfluidity we need to calculate the superfluid response function in the disordered phase, i.e., for the one-dimensional gas. The asymptotic behavior of the equal time correlation function at large distances can be determined using conformal field theory. Zero-field density profile N t /L within the local particle density approximation as a function of the position along the trap x for S = 1/2 (black), S = 3/2 (red), S = 5/2 (green) and S = 7/2 (blue). N t has been scaled with |c| so that the plot is universal. We have chosen μ(L/2) = −S(S + 1)c 2 /3 and μ(0) = 1.0c 2 .

Conformal towers
The low-lying excitations of the system are described by the finite size corrections to the ground state energy [36,37,38]. While the ground state energy for the bulk, E GS given by Eq. (9), is an extensive quantity, the mesoscopic corrections are of order 1/L, where L is the length of the system. In zero magnetic field, we need to consider only the contribution of the populated Dirac sea of bound states of N atoms to the finite size corrections. For periodic boundary conditions, four quantum numbers determine the excitations of this branch, namely, ΔN 2S corresponds to the added or removed number of bound states, D 2S is the parity variable, i.e. 2D 2S is the difference between forward and backward movers, and the n ± 2S count the number of particle and hole excitations about each Fermi point (+ for forward movers and − for backward movers). The cases S = 1/2 [27,39,40], S = 3/2 [27] and S = 5/2 [28] in a magnetic field have been extensively investigated previously.
The finite size corrections to the ground state energy are given by [27,36,37,38,39,40] (13) and the corresponding change in momentum for the excitations is In Eq. (13) v 2S is the group velocity and z 2S,2S is the generalized dressed charge for the band.
Note that D 2S is only determined modulo 1.
The dressed generalized charge is obtained as z 2S,2S = Ξ 2S,2S (B 2S ), where the Ξ 2S,2S is the solution of the integral equation Here the integration kernel is the same as for the integral equation for ε (2S) (ξ), i.e. Eq. (5). Comparing Eqs. (7) and (16) it follows that Ξ 2S,  Fig. 3(a) for various spin values. Note that the dependence of z 2S,2S with the integration limit increases dramatically with increasing S. For very large B 2S the Fredholm integral equation (16) can be reduced to a Wiener-Hopf integral equation, by shifting the variable ξ by the amount B 2S . The Wiener-Hopf equation can be solved analytically with standard methods (see e.g. Appendix B of [33]). This way we obtain asymptotically z 2S,2S (B 2S → ∞) = 1/ √ 2S + 1.

Asymptotic behavior of correlation functions
The asymptotic long-time and large distance behavior of the correlation function of a given conformal field operator O can be obtained using conformal field theory. The asymptote is determined by the low energy excitations of the system as a function of the set of quantum numbers ΔN 2S , D 2S , n + 2S and n − 2S . At T = 0 the space and time dependent correlation function for the operator O is given by [37,38,39,41,42] where p F,2S = πn 2S is the Fermi momentum defined previously. Here Δ ± 2S are the conformal dimensions that are obtained from the finite size corrections to the ground state energy, Eq. (13), and the momentum, Eq. (14), for the excitations [37,38,39,41,42] There are two factors in the denominator of Eq. (17), one for forward moving excitations and one for backward movers. The smallest conformal dimensions are obtained for n ± 2S = 0, which correspond to the smallest exponents in the correlation functions. Hence, the correlations with the longest range correspond to n ± 2S = 0. The above results are easily extended to finite temperature T by replacing each factor (x ± ivt) −2Δ ± in Eq. (17) by but this extension is not needed here. The leading term of the equal time correlation function for the operator O is then of the form where the amplitude A cannot be determined from conformal field theory. The exponent θ 2S is given by The factor cos(2D 2S p F,2S x) reflects the space modulation of the order parameter. For superfluidity in a finite magnetic field it arises from the spin-imbalance in the system and is reminiscent of the FFLO phase, where λ = 1/(2πD 2S n 2S ) is the distance between nodes. For a particle density wave response function, on the other hand, it reflects the modulation of the density wave.

Response function for superfluidity
Weak Josephson tunneling between tubes and interactions between particles in different tubes in an array of tubes [23,26], may give rise to a dimensional crossover from one-dimension to a higher dimension, yielding superfluid long-range order. The superfluid order parameter in zero magnetic field is a generalization of Cooper pairs for S = 1/2 to higher spin, in this case the bound states of N = 2S + 1 atoms. In the presence of a magnetic field there are many possible order parameters, which have been discussed in [27,28]. The correlation function for the bound states determines if long-range order is a possibility. A small exponent θ 2S favors superfluidity, in particular, the exponent θ 2S should be smaller than 2, which is the critical scaling dimension in one-dimension. This criterion is simple for the zero-field case, but has to be complemented by other criteria if the system is spin-imbalanced. Correlation functions for pure phases have been studied previously for other models [42,43]. The operator O corresponds then to adding or removing one cluster, i.e. ΔN 2S = ±1. Consequently, from Eq. (15) for a half-integer spin D 2S is an integer. The smallest exponent corresponds to D 2S = 0; it is easily seen that D 2S = ±1 yields a larger θ 2S . From Eqs. (18) and (21) we obtain that θ 2S = 1/(2z 2 2S,2S ). The exponent θ 2S as a function of B 2S /|c| for various spin values is shown in Fig. 3(b). and S = 7/2 (blue). The integration limit has been scaled with |c| so that the plot is universal. The dashed line at θ = 2 corresponds to the critical scaling dimension in 1D. The dotted line, θ = 1, represents equal exponents for superfluidity and density waves and is discussed in the text.
In the limit B 2S → 0 the exponent θ 2S tends to the value 0.5 for all values of S. With increasing atom density in the tube (increasing B 2S ) the exponent increases. In the limit B 2S → ∞ we obtain θ 2S = (2S + 1)/2. θ 2S increases monotonically with B 2S from 0.5 to (2S + 1)/2. Hence, for S = 1/2 and S = 3/2 the exponent θ 2S is always less than the critical scaling dimension and consequently the situation is favorable for superfluidity for all densities. The critical scaling dimension in 1D is 2, so that for θ > 2 the correlation function falls off too fast to yield long-range order. For S > 3/2 the exponent θ 2S is larger than 2 for larger densities. Hence, if S > 3/2 only low atom densities (or strong interactions) are favorable for superfluid long-range order in an array of tubes with Josephson tunneling. θ = 2 is indicated by the faint dashed line in Fig. 3(b).

Response function for density waves
Particle density waves oscillate in space with a wave number 2p F,2S . The density wave operator does not change the number of particles, so that ΔN 2S = 0 and consequently D 2S is an integer. We choose D 2S = ±1, which corresponds to oscillations with 2p F,2S , and n ± 2S = 0. The equal time correlation function is then proportional to i.e. the exponent is just the inverse of the one of superfluidity. Hence, θ DW decreases from the value 2 as B 2S increases. Density waves are then, in principle, favorable for all spins and all band fillings. We now have to compare θ 2S with θ DW . The simplest criterion is that the smallest exponent yields the dominant long-range order. The exponents are equal if θ 2S = θ DW = 1 or z 2S,2S = 1/ √ 2. The corresponding line is shown as dotted in Fig. 3(b). We may conclude that for S = 1/2 superfluidity is more favorable than density waves for all particle densities and interaction strengths. For S ≥ 3/2 there is a crossover from superfluidity at low densities to density waves at large B 2S /|c|. The transition is where the exponents are equal one.
The focus of our previous papers [7,27,28] was the possibility of finding inhomogeneous phases of two types in the gas in a magnetic field: (a) We considered the scenario of phase separation along the tube, and (b) modulations of the order parameters of the FFLO type.
In case (a), the confining harmonic potential varies with the position x along the tube. Within the local density approximation, which absorbs this variation into the chemical potential, μ is a function of x and hence different phases for a given magnetic field appear along the trap giving rise to phase separation [1,2,3,7,28]. In zero-field there is only one phase, so that in the present study, such a phase separation cannot occur. The density profile, however, changes dramatically along the tube, being larger at the center of the tube than at its boundaries.
In case (b), inhomogeneities like modulations of the order parameter of the FFLO type, may arise due to the spin-imbalance of the atoms in a magnetic field. In general this requires an array of tubes with interactions between particles in different tubes or Josephson tunneling between tubes [26,23,27] giving rise to a dimensional crossover from one-dimension to a three dimensions. The dimensional crossover opens the door to the possibility of superfluid long-range order at finite temperature. There is a very large number of order parameters for superfluidity for the spin-imbalanced gas, but only one in zero magnetic field.
The response functions in the disordered phase are determined by the Luttinger properties of the one-dimensional gas. They determine the instability toward superfluidity from the normal phase. A combination of three criteria determine the dominating order in the case of a spin imbalanced gas: (i) the smallest exponent θ corresponds to the longest range and is hence favorable to order if it is less than 2, the critical scaling dimension, (ii) a large distance between nodes, λ, favors order because it is energetically less unfavorable, and (iii) the preformed bound states should carry a small momentum, since a large momentum of the bound states is energetically unfavorable to a condensate. Cooper-pairing of atoms with spin-component S and S − 1 (one particle being a forward mover and the other a backward mover) was found to be the first superfluid phase accessed from the normal state as the temperature is lowered in a sufficiently large field. Note that the two main conditions for realization of the FFLO phase are satisfied in cold atom tubes: (1) the system is very pure (no impurities) and (2) it has a low effective dimension (extreme anisotropy). In a sufficiently strong magnetic field there are in general several bands of bound states represented in the ground state. Each band contributes additively to the critical exponent of a correlation function. Hence, critical exponents tend to be larger in a magnetic field and superfluidity is less favorable than in zero-field.
In zero-field there is only one operator for superfluidity, namely the creation (annihilation) operator for one preformed bound state of (2S + 1) atoms. The superfluid phase is homogeneous and not of the FFLO type. For S = 1/2 superfluidity is favorable because the critical exponent is smaller than one and smaller than θ DW for all densities and coupling strengths. For S ≥ 3/2, on the other hand, superfluidity is only possible at sufficiently low atom densities and/or strong interactions |c|, while for large B 2S /|c| density waves are the likely candidate for long-range order. The critical exponent grows quite rapidly with increasing spin size and particle density, and eventually θ 2S crosses the value one (dotted line in Fig. 3(b)) and beyond that point θ DW = 1/θ 2S will be less than one and density waves are the dominating feature.
It should be mentioned that a related SU(4) invariant Hubbard model (attractive U ) has been studied on a lattice [44]. The methods employed are the density-matrix renormalization group technique and a phenomenological Luttinger liquid, where the charge sector is represented by a free massless bosonic field. Their phase diagram (S = 3/2) displays two regions, a superfluid quartet phase at low density, which gives way to a charge density wave phase at larger densities. This crossover and the trend that superfluidity and DW exclude each other is confirmed by our exact calculation for N = 4. A quantitative comparison is more difficult because the models are different.
A possible application of a SU(4) system with attractive interactions is the alpha-particle condensation in nuclear matter. The four degrees of freedom are made up of the spin and the pseudospin. The alpha particles are strongly bound and exist in a Bose-Einstein condensation phase, which implies a low density. Of course a three-dimensional system has to be considered here. An excellent review on the status of the subject can be found in [45].