Exotic pairing in 1D spin-3/2 atomic gases with $SO(4)$ symmetry

Tuning interactions in the spin singlet and quintet channels of two colliding atoms could change the symmetry of the one-dimensional spin-3/2 fermionic systems of ultracold atoms while preserving the integrability. Here we find a novel $SO(4)$ symmetry integrable point in thespin-3/2 Fermi gas and derive the exact solution of the model using the Bethe ansatz. In contrast to the model with $SU(4)$ and $SO(5)$ symmetries, the present model with $SO(4)$ symmetry preserves spin singlet and quintet Cooper pairs in two sets of $SU(2)\otimes SU(2)$ spin subspaces. We obtain full phase diagrams, including the Fulde-Ferrel-Larkin-Ovchinnikov like pair correlations, spin excitations and quantum criticality through the generalized Yang-Yang thermodynamic equations. In particular, various correlation functions are calculated by using finite-size corrections in the frame work of conformal field theory. Moreover, within the local density approximation, we further find that spin singlet and quintet pairs form subtle multiple shell structures in density profiles of the trapped gas.


Introduction
Large-spin atomic fermions displaying rich pairing structures and diverse many-body phenomena could be realized through controlling interactions in spin scattering channels. Promising theoretical progress has been made on large-spin atomic fermions in the context of multi-particle clustering superfluidity [1, 2,3,4]. In particular, fermionic alkaline-earth atoms can exhibit an exact SU(κ) symmetry with κ = 2I +1 [5,6], where I is pure nuclear spin. Experimental explorations of these fermionic systems have been reported on the trimer state of 6 Li atoms [7,8,9], the SU(10) symmetry fermionic gas of 87 Sr atoms with I = 9/2 [10] and the two-orbital magnetism of SU(κ) symmetry in 87 Sr atoms [11], the SU(2) ⊗ SU(6) symmetry fermionic atoms of 173 Yb with its spin-1/2 isotope [12], the SU(6) Mott-insulator state of 173 Yb atoms [13] and two-orbital Hubbard model [14], etc. Despite much theoretical and experimental efforts on the large spin atomic systems, understanding spin pairs and large spin magnetism is still rather elusive [15], see a recent review [16].
In contrast to the conventional spin-1/2 electronic magnetism [17,18], one-dimensional (1D) large-spin ultracold atomic fermions exhibit richer high spin phenomena [19,20]. Preliminary study of multi-colour superfluidity and quartetting ordering [21,22] was also carried out through few simplest integrable large spin fermionic systems, such as the SU(4)-and SO(5)-invariant spin-3/2 fermions [23,24,25]. However, there is still little progress toward to the understanding of large spin pairing in non-SU(κ) symmetry models [26,27,28]. To this regard, the spin-3/2 fermionic system comprises an ideal model towards to the precise understanding of large spin non-SU(κ) symmetries. From experimental point of review, spin-3/2 systems can be realized with alkali atoms of 6 Li, 132 Cs and alkaline-earth atoms of 9 Be, 135 Ba and 137 Ba [19], see the recent experiment [29] for more details.
Advances in controlling ultracold atoms provide us promising opportunities of realizing quantum many-body systems with more tunable parameters [30,31]. For example, S-wave scattering in these spin-3/2 atomic fermionic systems acquires six effective interaction channels according to the Clebsch-Gordan coefficients. One is the spin singlet channel with total spin-0 and the others are the spin quintet channels with total spin-2, while the S-wave scattering in the channels of total spin-1 and -3 is forbidden. When the interactions in all scattering channels are the same, the system exhibits SU(4) symmetry and is integrable in 1D [23,24], see model (i) in Table I. However, if the interaction in quintet channels is different from that in the singlet channel, then the system can exhibit a SO(5) symmetry without fine tuning [19], see model (ii) in Table I. The model with SO (5) symmetry is integrable at a special point [27]. No. g 00 g 2,2 g 2,1 g 2,0 g 2,−1 g In fact, the scattering lengths in the quintet channels may be tuned to breakdown the SU(4) or SO(5) symmetries while preserving the integrability. In this paper, we propose an integrable spin-3/2 atomic Fermi gas with SO(4) symmetry, see model (iii) in Table I. Based on the exact solution, we find that the model has subtle spin pairing phases and exhibits Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) like pair correlations. For SU(4) symmetric model, there are no pairs when the interactions are repulsive, while for attractive interactions, four-body bound states appear in the ground state [23,24]. With the competition between chemical µ and external magnetic field h, three-body bound states and two-body pairs also could emerge in the ground state [25]. However, for the present SO(4) symmetry model with either repulsive interaction (c > 0) or attractive interaction (c < 0), bound pairs could emerge in the ground state regardless of the sign of c. In contrast to the SU(4) case, we also find that there are no four-and three-body bound states in the SO(4) system. We further show that the integrable SO(4) model presents universal quantum criticality with dynamic exponent z = 2 and correlation exponent ν = 1/2 in terms of different pairing states.
The paper is organized as following. In Section 2, we discuss the symmetry and conserved quantities of the integrable SO(4) model in the frame work of Yang-Baxter equation. The exact solution, thermodynamic limit and thermodynamic Bethe anatz equations are derived in Section 3. In Section 4, ground state properties and full phase diagrams are discussed in grand canonical ensemble. The elementary excitations of this model and Luttinger liquid behaviour are studied in detail in Section 5. In Section 6, various asymptotics of correlation functions for several quantum phases in strong coupling limit are calculated explicitly. The equation of states and the quantum criticality of the system are further studied in Section 7, followed by discussions in Section 8.

Integrable SO(4) model and conserved quantities
We consider the model Hamiltonian that describes dilute spin-3/2 atomic gases of N fermions with contact interaction constrained by periodic boundary conditions to a line of length L.
Here we let 2m = = 1. The last term in the Hamiltonian (2.1) is the Zeeman energy. The spin polarization (magnetization) is given byM = jf z j , wheref z is hyperfine spin moment in he z-direction. The projection opera-torP lm ij = |lm lm| projects the total spin-l state onto the spin-m state in the z-direction of two colliding atoms i and j. In the above equation, the summation lm is carried out for l = 0, 2 and m = −l, −l + 1, · · · , l. The interaction strength in the channel |lm lm| is given by g lm = −2 2 /(ma lm 1D ), where is the effective 1D scattering length depending on the 3D scattering length a lm in the channel |lm lm| [32]. The constant C ≈ 1.4603 and (l, m) stands for the interaction channel. Remarkably, the model (2.1) exhibits three classes of mathematical symmetries that preserve the integrability by appropriate choices of the interaction strength via c lm = −2/a lm 1D = mg lm / 2 , see Table  1.
The attractive potentials may lead to spin-J pairs with J = 0, 2 in the quasi-momentum space. Consequently, the generatorsφ † of these spin-J pairs comprise the following two sets of representations regarding to the attractive channelsV 2 andV 1φ Using the language of second quantization, the pseudo spin operators can be expressed asĵ d,α =ψ † j d,αψ , whereψ = (ψ 3/2 ,ψ 1/2 ,ψ −1/2 ,ψ −3/2 ) t and the superscript t denote the transposition. One can check that following commutation relations are valid The pseudo spin operators are commutative with the interaction potential. Therefore, the model has SO(4) symmetry. This symmetry does not preserve the number of atoms with each internal degree of freedom. But the total number of atoms are conserved. Meanwhile, the quantitieŝ are conserved. HereN i is the number of atoms with spin-i component. The spin polarization along the z-direction, which is given byM = 3 2Ĵ 3/2 + 1 2Ĵ 1/2 , is also conserved.

The Bethe ansatz solution and thermodynamic limit
The model (2.1) with SO(4) symmetry can be solved by means of the coordinate Bethe ansatz (BA) [33,17]. The BA wave function of the model reads where x j and σ j are the coordinate and the spin in z-direction of the j-th atom respectively. Here k j are pseudo-momenta of the particles, j = 1, 2, · · · , N, and Q and P are the permutations of {1, 2, · · · , N}. Θ(Q) is a series of production of step functions, i.e.
. As usual, θ(x) = 1 when x ≥ 1 and θ(x) = 0 otherwise. In particular, the key ingredient of BA is the superposition coefficients of these plane waves A(Q, P), which can be consequently determined by the two-body scattering relation below.
The wave function (3.1) of the many-body interacting fermions (2.1) should satisfy the fermionic statistics and time independent Schrödinger equationĤΨ = EΨ . These restrictions lead to the scattering relation between two superposition coefficients A(Q (ba) , P (ba) ) =Ŝ Qa,Q b (k Pa − k P b ) A(Q (ab) , P (ab) ), where the two-body scattering matrix satisfies the Yang-Baxter equation S ab (λ)S ac (λ + ν)S bc (ν) = S bc (ν)S ac (λ + ν)S ab (λ), (3.3) which guarantees the integrability of the model. In the above equations we denoted Q (ab) = {· · · , Q a−1 , Q a Q b , Q b+1 , · · · }, Q (ba) = {· · · , Q a−1 , Q b Q a , Q b+1 , · · · }, and P (ab) and P (ba) have a similar definition. Eq. (3.2) indicates that there is no interaction in the total spin-1 and -3 channels. This is the consequence of the symmetry of the wave function. The periodic boundary conditions give rise to the following eigenvalue problem which is used to determine the quasi-momenta ks. By using the nested algebraic BA [34,35] and after complicated calculations, we obtain the energy of the system as E = N j=1 k 2 j − hM, where the quasi-momenta {k i } should satisfy following BA equations, Here λ and ν are the spin rapidities, and we denoted the function e a (x) = x−iac x+iac . The BA equations coincide with the two-band model of electrons in 1D [36]. From these BA equations we can obtain the full phase diagrams and thermodynamics of the spin-3/2 interacting fermions with SO(4) symmetry. In Eq. (3.5), the quasi-momenta {k i } are nested with the spin rapidities λ and ν. For c > 0, the spin rapidity λ provides the spin flip with a total spin change of 1 whereas ν gives the spin flip with a total spin change of 2, see Finding root patterns of the BA equations (3.5) presents a big theoretical challenge towards to understanding the physics of the model. Here we find that the quasi-momenta {k i } can be either real or complex conjugated pairs at the thermodynamic limit, L → ∞ and n = N/L is a constant. The real roots k u,z ∈ R with z = 1, 2, · · · , N u characterize the quasi-momenta for unpaired fermions. The complex conjugated roots k p,z ∈ C with z = 1, 2, · · · , N p denote the quasi-momenta of bound states, i.e., spin-J pairs with total spin J = 0, 2. We observe that the BA equations (3.5) present the spin pair bound states solely depending on the attractive interaction channels ofV 1 andV 2 (or say c < 0 or c > 0), see Eqs. (2.8) and (2.9). Explicitly, for c > 0, the momenta for a bound pair of two-atom with different spins have a complex conjugate pattern k ± p,z = k p,z ± ic/2 with a binding energy ǫ p = c 2 /2 in the thermodynamic limit. Where a real spin rapidity λ p,z = k p,z is coupled together with k ± p,z . So do the spin rapidities ν p,z in the case of c < 0. For the ground state, we observe that the BA roots only involve real k, complex conjugated k, 1-string ν and 2-string ν for c > 0, see Section 4. Here we denote the corresponding density distribution functions as ρ (u) , ρ (p) , ρ (1) and ρ (2) . Whereas, for c < 0, the 1-string and 2-string λ populate the ground state instead of the ν. At finite temperatures, the spin rapidities λ's and ν's evolve different lengths of spin strings similar to the spin-1/2 Heisenberg integrable model [37]: λ n,z,j = λ n,z + 1 2 (n + 1 − 2j)i|c|, z = 1, 2, · · · , M 1,n , ν n,z,j = ν n,z + 1 2 (n + 1 − 2j)i|c|, z = 1, 2, · · · , M 2,n , λ n,z , ν n,z ∈ R, j = 1, 2, · · · n, (3 .6) where M 1,n and M 2,n are the numbers of n-string of λ and ν, respectively. Thus we have n M l,n = M l (l = 1, 2). The BA equations in thermodynamic limit are obtained by substituting the string hypothesis into eqs. (3.5). By taking the logarithm of the BA equations and using the relation when c > 0, we get the logarithmic form of the BA equations Here the quantum number I (u), (p), (1), (2) takes integers or half odd integers as the following The BA wave function acquires that any two rapidities of each branch are not equal. Otherwise the wave function is zero. According to this restriction, the quantum numbers Is of each branch take different numbers in parameter space [34]. The solutions of the BA equations (3.5) are classified by a set of quantum number I (u), (p), (1), (2) . Thus the momentum of the BA eigenstate is given by Accordingly, the energy E and the magnetization M of the model are given by respectively. At finite temperature, some quantum numbers Is (called vacancies) are occupied in parameter spaces whereas some of those quantum numbers are not occupied. The unoccupied BA roots are called unoccupied vacancies. In spin sector, the unoccupied vacancies are regarded as holes of the spin strings. The spin strings are nothing but spin wave bound states. For thermodynamic limit, i.e. L, N → ∞ and N/L is finite, the strings and holes can be treated as density distribution functions ρ(k) and ρ h (k), respectively. It follows that d dk . Consequently, we obtain the integral form of the BA equations in the thermodynamic limit, Here for c > 0, ρ 1,n and ρ 2,n are the string densities of n-strings of λ and nstrings of ν, respectively. Whereas for c < 0, ρ 1,n and ρ 2,n are the densities of n-strings of ν and n-string λ, respectively. The * denote the convolutionf * and a n (k) = 1 π n|c| (nc) 2 +k 2 . The ρ l n and ρ l n h are the densities of particle and hole for the length n-strings in spin sector. The energy, momentum and total spin of the Hamiltonian are given by where α 1 = 1, α 2 = 2 for c > 0, whereas α 1 = 2, α 2 = 1 for c < 0. We also denote n u,p = dkρ u,p (k) as the density of atoms in the scattering state and the density of the pairs, respectively. Moreover, m l,n = dkρ l,n (k) with l = 1, 2 are the densities of n-strings in the spin sector. At finite temperature T , following Yang and Yang's grand canonical description [38], the thermodynamic Bethe ansatz (TBA) equations can be obtained from minimization of the Gibbs free energy Ω = E − T S − µN with S the entropy of the system. Here we define the dressed energies of the charge rapidities as ε α (k) = T ln ρ α h (k)/ρ α (k) with α = u, p for unpaired fermions and pairs, respectively, thus the TBA equations are given by where ε − = − ln(1 + ρ/ρ h ). The effective chemical potentials of unpaired fermions and pairs are defined by µ u = µ + (α 1 + α 2 )h/2 and µ p = c 2 /2 + 2µ + α 2 h. For c > 0, ε 1,n and ε 2,n are the dressed energies for n-strings of λ and n-strings of ν, respectively. Whereas for c < 0, ε 1,n and ε 2,n are the dressed energies for n-strings of ν and n-strings of λ, respectively.
The TBA equations (3.15) involve infinite branches of ε l,n coupled together. This imposes a big theoretical challenge to find the exact solutions of the TBA equations (3.15). We observe that the dressed energy of the spin string tends to be a constant when the spin string length increases. The larger strings contribute the smaller energy as the temperature decreases. In order to discuss the solutions of the TBA equations, we present the TBA equations (3.15) in the recursive form Here the convolution kernel is given by G(k) = (1/2|c|)sech(πk/c) and the function ε + is defined by ε + (k) = T ln(1 + ρ h /ρ). The thermal potential per unit length is p = p u + p p and the effective pressures are given by for the unpaired fermions (the bound pairs). The TBA equations (3.16)-(3.23) provide not only ground state properties at zero temperatures but also full thermodynamics at finite temperatures. In the following Sections, we will solve the complicated TBA equations and discuss the physics of the spin-3/2 fermionic system with SO(4) symmetry.

Quantum phase diagram and pairing signature
In this section, we will study the phase diagram and the pairing signature of the model. A rigorous way of finding the ground state in grand canonical ensemble is to take the zero temperature limit in the recursive TBA eqs. (3.16)-(3.23). In the limit of T → 0, we have ε + (k) ≥ 0 and ε − (k) ≤ 0. For a nonzero magnetic field h and c > 0, from Eqs. (3.23), we observe that the large strings of λ and ν always have positive dressed energies. Therefore there are no such strings in the ground state. From Eqs. (3.21)- (3.22), it is obvious that ε 1,n (n ≥ 2) and ε 2,n (n ≥ 3) are positive so that the corresponding strings do not populate in the ground state. Even the magnetic field tends to zero, the dressed energies of these strings approach to 0 + . Therefore there are still no such strings in the ground state. The Eq. (3.18) also gives a nonnegative dressed energy for the real λ (or ν) rapidity for c > 0 (or c < 0). Thus the dressed energy ε 1,1 is also excluded from the ground state.
In the T → 0 limit, the TBA equations (3.21)-(3.22) reduce to a new set of dressed energy equations that characterize the Fermi seas of the paired states and single atoms in terms of chemical potential µ and magnetic field h. The band filling of those Fermi seas with respect to µ and h provide an analytical way to determine the full phase diagram in the µ-h plane. For c > 0, the dressed energy equations involve quasi-momenta k u and k p , and spin rapidities ν 1 and ν 2 . The corresponding density distributions of them are denoted as ρ u , ρ p , ρ (1) and ρ (2) . For c < 0, the dressed energy equations involve k u , k p , ν 1 and ν 2 . We can define the ground state densities in a similar way as the ones for c > 0. The densities of these rapidities satisfy the following integral BA equations h ) t , ρ = (ρ u , ρ p , ρ (1) , ρ (2) ) t , ρ 0 = (1/2π, 1/π, 0, 0) t . Here the superscript t means the transposition and integral kernelK is given byK The TBA eqs. (3.16)-(3.23) in the T → 0 limit reduce to the dressed energy equations which can be written in the following form − ) t and ε 0 = (k 2 − µ 1 , 2k 2 − µ 2 , α 1 h, α 2 h) t . These equations are very convenient to analyze quantum phase diagram and quantum phase transitions in terms of chemical potential and magnetic field. We observe that only the atomic pairs and the unpaired atoms populate the ground state. The absence of three-or four-body charge bound states is the consequence of the repulsive interacting channels due to the SO(4) symmetry. For c > 0, the existing pairs are φ † 2,±2 and φ † 2,0 , whereas for c < 0, the paired states comprise φ † 2,±1 and φ † 0,0 . The competition among the pairing, chemical potential µ and external magnetic field h gives rise to a rich phase diagram of this model. Following the method presented in [39], full phase diagrams of the model can be both analytically and numerically worked out from the density equations (4.1) or from the dressed energy equations (4.3). The phase diagram of the system with c > 0 is shown in the left panel in Fig.  2. For c < 0, different BA root patterns lead to a different phase diagram which is presented in the right panel in Fig. 2.
In the following, we only consider the regime for c > 0. In Fig. 2, V stands for the vacuum. The phase FF denotes a fully-polarized phase of single |3/2 atoms in the region of strong magnetic field and small chemical potential. The FFLO-like phase composes of the largest spin-2 component pairsφ 2,2 and excess fermions of |3/2 atoms that form two mismatched Fermi surfaces in quasi-momentum space. In this phase, atoms can be in either the scattering state |3/2 or bounded pair stateφ 2,2 . We show that spatial oscillation of the pair correlation function solely depends on the mismatch between the Fermi surfaces of |3/2 and |1/2 atoms, i.e. ∆k F = π(n 3/2 − n 1/2 ). In the strong coupling region, the slowest decaying term of the pair correlation function G p (x, t) = G|φ † 2,2 (x, t)φ 2,2 (0, 0)|G is obtained, G p ≈ A 0 cos(π∆k F x)/(|x + iv u t| θu |x + iv p t| θp ), where θ u = 1/2, θ p = 1/2 + n p /c and n u,p are the densities of unpaired |3/2 atoms and atomic pairsφ 2,2 , respectively. More detailed analysis will be further discussed in Section 6. Thus the momentum pair distribution has peaks at the mismatch of the Fermi surfaces ∆k F , which presents a characteristic of the FFLO phase.
Moreover, the FP phase is the fully-paired stateφ 2,2 , where the pair correlation function decays as a power of distance. The leading term of pair correlation function G p is given by G p = A/(|x + iv p t| θp ) with θ p = (1/2)(1+ 1/γ). This indicates that the pair correlation in the fully-paired phase has a longer correlation length than the one in the phase FFLO. According to the Luttinger liquid theory, the long distance behavior of the pair correlation of this type gives rise to the algebraic form G p ≈ |x| −1/(2K) , where the parameter K ≈ 1 − 1/γ. However, the correlation function for the one particle Green's function decays exponentially. When the interaction is very strong, i.e. near the phase transition line from vacuum into the FP phase, the system behaves like a Fermionic super Tonks-Girardeau gas [40,41].
The MP phase denotes a fully paired phase, but not a fully-polarizedφ 2,2 state. For c > 0, three kinds of spin pairsφ 2,2 ,φ 2,0 andφ 2,−2 are involved in the MP phase. By increasing the chemical potential µ and magnetic field h, the system may enter into three other mixed phases denoted as V, FF and FFLO Phases, see Fig. 2. The phase transition from FFLO phase into V phase occurs as the 2-string ν emerges in the ground state. While the one from MP phase into V phase is driven by the appearance of real k, i.e. involvement of pair breaking. The one from V phase into FFLO phase is driven by spin flipping processes involving the real ν and the one from FF phase into FFLO phase is driven by spin flipping process involving the 2-string ν.
In a harmonic trap, quantum phase segments in density profiles can be used to identify different quantum phases in experiments with cold atoms. We can extract the threshold values of the phase boundaries in Fig. 2 through the density profiles of the trapped gas. Within the local density approximation, the local chemical potential is replaced by µ(x) = µ 0 − 1 2 ω 2 0 x 2 for the trapped gas. The total density of fermions is given by n = n u + 2n p . For fixed particle numbers, we plot the density profiles of the trapped gas in Fig.  3. In contrast to spin-1/2 Fermi gas [42,43,44,45], spin quintet pairs lead to multiple shell structures in the density profiles. The trapping center can be a mixture of three types of spin pairs as well as unpaired single atoms accompanied by the different shells involving the fully-polarized spin-2 pairs, the FFLO state, and the FF phase, see Fig. 3.

Elementary spin and charge excitations
We consider the elementary excitations in the fully paired MP phase. Suppose that the external magnetic field is zero. The excitations in the charge sector of the system are similar to the ones in the SU(4) interacting fermions. However, the spin excitations in the spin sector are quite different from those of SU(4) Fermi gas. In the fully paired MP phase, the elementary charge excitation has an energy gap. In contrast to the SU(2) attractive Fermi gas [46], there exist gapless low-lying spin excitations triggered by changing the quantum number of 2-string ν in this paired MP phase. The change of the number of 2-string ν results in necessary spin configuration changes in the paired states ofφ 2,2 ,φ 2,0 andφ 2,−2 . In particularly, breaking down a 2-string ν leads to the total spin change of 4, i.e. creation of 4 spinons, see Fig. 5. Each spinon carries spin-1, and we denote this type of excitations as large spinon excitations. There is no such kind of spin excitation in the attractive SU(2) Fermi gas. These spinons can be directly calculated from the BA equation (4.1). One 2-string of ν is excited that leads to four spin-1 holes of ν: δM 2 = δM h 2 /4. If one adds d holes of 2-string ν into the ground state, the excited energy and momentum are given by respectively. Where ν h j are the positions of the added holes. For h = 0 and γ = c/n = 1, the 4-hole exciting spectra is shown in Fig. 4(b). However, for c < 0, one 2-string λ excitation splits into four spin-1/2 holes which give an excitation with total spin-2.
In order for our convenience to discuss correlation function, we will focus on three phases FF, FFLO and FP for the strong coupling regime and a finite external magnetic field h. These phases near the vacuum phase boundary in  It is important to note that in these three phases both spin rapidities of ν and λ are absent in the ground state. Thus the dressed energy equations (4.3) reduce to the similar structure as the one for the attractive spin-1/2 Fermi gas [39,47,25]. However, spin excitations and finite temperature behaviour of these states in the SO(4) symmetry model are quite different from that of the attractive spin-1/2 Fermi gas. In general, at zero temperature, the Fermi points Q u,p are determined by the condition ε u (±Q u ) = 0 and ε p (±Q p ) = 0, where Q u,p > 0. At the ground state, all the quantum numbers are filled up to the Fermi points, i.e. I u (k u,z ) − I u (k u,z−1 ) = 1 and I p (k p,z ) − I p (k p,z−1 ) = 1 with k u,z > k u,z−1 and k p,z > k u,z−1 . All excitations can be classified as three types of excitations, i.e. particle-hole excitations, backwards scatting process and adding particles at different Fermi points. We denote these excitations as Type 1, Type II and Type III, respectively, see Fig. 5.
The Type I elementary excitations are characterized by moving atoms (pairs) close the left or right pseudo Fermi point outside the Fermi sea, see Fig. 5 (1). For example, we consider the excitations close to the right Fermi point Q u in real quasi-momentum space k. The process is to create N + u holes at positions I u (k u,Nu ), I u (k u,Nu ) − 1, · · · , I u (k u,Nu ) − N + u + 1 and to add N + u real ks outside the Fermi sea where the corresponding quantum numbers are I u (k u,Nu ) + 1, I u (k u,Nu ) + 2, · · · , I u (k u,Nu ) + N + u . Based on this setting, one can perform standard calculation of finite-size corrections to the energy and momentum [48,47,25]. This type of excitations can take place close to the left pseudo Fermi point −Q u with adding N − u holes below the Fermi point. In the branch of atomic pairs, particle-hole excitations also occur close to the left and right pseudo Fermi points with adding the quantum numbers N − p and N + p , respectively. The Type II excitations arise from the changes of total number of unpaired fermions or bound pairs. This type of excitations are characterized by adding (or removing) atoms/pairs close to the Fermi points, see Fig. 5 (2). We denote the change of particle number as ∆N u (or ∆N p ).
The Type III excitations are created by moving particles from the left Fermi point to the right Fermi point and vice versa. This process is also known as backscattering, see Fig. 5 (3). These two types of excitations change the pseudo Fermi points. The quasi-momenta of single atoms/pairs with positive real parts are viewed as the right-going atoms/pairs whereas the ones with negative real parts are denoted as the left-going atoms/pairs. Both the Type II and the Type III excitations lead to the particle number difference between the right-and left-going atoms/pairs. We denote these number differences as 2∆D u and 2∆D p for the unpaired atoms and the pairs, respectively.
Based on the above descriptions, all the three types of excitations can be unified in the following form of the finite-size corrections for the total momentum and energy of elementary excitations where the indices α, β, γ = u, p and v β is the Fermi velocities of the corresponding states. In the above equations Z = Z (w) is the dressed charge which is determined by For the pure paired phase, w=FP and we have For the phase, w=FFLO and we have , The value of dressed charges can be obtained by solving the dressed charge equation (5.4) numerically. The asymptotic forms of the dressed charges in the strong coupling regime is given in Section 6. On the other hand, the excited energy and excited momentum can also be expressed by the conformal dimensions ∆ ± according to the conformal field theory [49,50] where k F,α is the Fermi momenta, k F,α = πN α /L. It follows that These relations provide us a analytical way to calculate the conformal dimensions via finite-size corrections

Asymptotics of correlation functions
In this section, we will calculate various correlation functions of the system with strong interactions in three phases FF, FFLO and FP. At zero temperature, all correlation functions reveal a power law decay in the phases FF and FFLO. In the pure paired FP phase the pair correlation decay as a power law of distance whereas the single particle Green's function decays exponentially. However, at finite temperatures, all correlation functions exponentially decay in long distance.
From the conformal field theory, the two-point correlation function for primary fields with the conformal dimensions ∆ ± is given by [49] Here the conformal dimensions ∆ ± are determined by (5.11) in terms of N ± , ∆N, ∆D and the dressed charge Z.
For strong interaction regime, the leading order of the dressed charge is given by Z αβ = δ α,β , where δ α,β is the Kronecker delta function. From this leading order of the dressed charge, we obtain the conformal dimensions as 2∆ ± α = 2N ± α +(∆D α ±∆N α /2) 2 +O(c −1 ). For the strong coupling regime, the order of 1/c corrections to the correlations are presented via the conformal dimensions where δ α with α = FF, FP and FFLO are the first order corrections, see Table 2.
In the fully-polarized phase FF, fermions are unpaired and there does not exist contact interaction. Thus the dressed charge reads Z u u = 1. In the phase FFLO, the dressed charges are given by For the fully-paired phase FP we have the dressed charge of pairs Z p p = 1 − 2Q p /πc + O(c −3 ).
Based on the above settings, we may calculate the asymptotics of the correlation functions G|Ô † (x, t)Ô(0, 0)|G . For a given correlation function, the operator O(x, t) act on the ground state that gives the quantum numbers ∆N u,p . Expanding the operator O(x, t) with respect to the primary fields with conformal dimensions ∆ ± and their descendent fields, we may calculate the asymptotics of the correlation functions by the formula (6.1). We present the conformal dimensions of various correlators in Table 2. In general, the real part of a correlation can be expressed as where λ s = ( α N α ∆D α /L) −1 is the wave length of the oscillation and the exponent θ α = 2∆ + α +2∆ − α . The Fourier transform of the correlation function near the Fermi velocity k 0 is given by [51] where ν s = 2∆ + + 2∆ − − 1, s = α ∆ + α − ∆ − α is the conformal spin and 2s is always integer.
The various correlation functions, for example, single particle Green's function, pair-pair correlation, density-density correlation are listed in Table 2. The fully polarized case FF is regarded as a single component free fermionic gas. The correlation function of the field operator is given by This excitation associates with the quantum number ∆N u = 1 and ∆D u ∈ Z+1/2. A few leading orders of the asymptotics of the single particle Green's function are obtained by choices of the values of the ∆D, see the rows 1-3 in the Table 2. The leading term in the correlation function G u is given by with the setting ∆D u = ±1/2 and N + u = N − u = 0, see the row-1 and the row-2 in Table 2. Here n u is the density of single atoms and A 0 is some constant. The row-1 shows a left-going wave with ∆D u = −1/2, or say that adding an unpaired atom near the left Fermi point. The row-2 indicates a right-going excitation wave with ∆D u = 1/2. We also consider the charge density correlation function G n (x, t) = n(x, t)n(0, 0) in this phase, namely where A n1 and A n2 are constants, see the rows 4-8 in the Table 2.
In phase FFLO, both the spin-2 pairs and single spin-3/2 atoms coexist. The single particle Green's function G u acquires ∆N u = 1, ∆N p = 0 and ∆D u,p ∈ Z + 1/2. The few leading orders of the asymptotics of G u are indicated in the rows 9-12 in Table 2. The first two leading terms read where A 1 and A 2 are constants and the exponents read In the correlation function (6.9), the first term presents the result listed in the rows 9 and 10 in Table 2, and the second term is presented in the rows 11 and 12. For finitely strong interaction, the single particle correlation function decays slower than that of the free Fermions due to the correlations between pairs and between pairs and unpaired atoms.
In the FFLO phase, the asymptotic pair-pair correlation is given by with the quantum numbers ∆N u = 0, ∆N p = 1, ∆D u ∈ Z + 1/2 and ∆D p ∈ Z. The leading order of the pair-pair correlation function is given explicitly by where A 3 is a constant and θ u3 = 1/2, θ p3 = 1/2 + n p /c, (6.13) see the rows 13 and 14 in Table 2. Here n 3/2 and n 1/2 are the densities of spin-3/2 and -1/2 atoms, respectively. We see that the long distance asymptotics for the pair correlation function oscillates with the wave number π(n 3/2 − n 1/2 ), a mismatch between the two Fermi surfaces of the two species. This is a characteristic of the FFLO-like correlation. We would also like to address that the FFLO-like oscillation terms arise from Type III excitations, i.e. backscattering for bound pairs and unpaired fermions. The charge density correlation function G n (x, t) is presented in the rows 17-19 in Table 2.
In the phase FP, the pair-pair correlation function G p acquires the quantum number ∆N p = 1. We present the result of the long distance asymptotics of pair-pair correlations in the rows 24-27 in Table 2. The leading term of the pair-pair correction reads where A 4 is a constant. This shows that the exponent θ p4 is smaller than that in the phase FFLO. In this phase the pair-pair correlation dominates the quantum correlation. Finally, we present the charge density-density correlations with the correlation g(x, t) = cos(2πn p x) |x + iv p t| θ 1 , θ 1 = 2 − 4n p /c, (6.16) where A n2 is constant. The related quantum numbers are listed in the rows 28-31 in Table 2.

Equation of state and quantum criticality
The 1D fermionic systems usually exhibit a rich resource of Luttinger liquids and show a novel magnetism in term of spin-change separation scenario, see recent review [46]. In order to understand large spin cold atoms in 1D, it is very important to investigate the low temperature behavior of Luttinger liquids. We will prove that the Luttinger liquids of different pairing states comprise a universal low energy physics of the model. In general, the specific heat can be written in terms of sound velocities of the liquid phases, i.e. c L = (πT L/3) α 1/v α , where the summation carries out over all the Fermi velocities in the charge and spin sectors. This result is valid for the temperature below a crossover temperature T * , i.e. T ≪ T * . Here the crossover temperature can be determined from the equation of states [52]. Near a critical point and for the temperature T ≫ T * , the system lies in the critical regime, where the crossover temperature characterizes the energy gap T * ∼ µ − µ c and/or T * ∼ H − H c . In the critical regime, universal scaling behaviour of thermodynamical properties is expected to distinguish quantum criticality with the critical dynamic exponent z = 2 from the Luttinger liquid criticality of z = 1 [39].
For strong coupling regimes, only charge rapidities are left at the ground states for the phases FF, FFLO and FP. For low temperatures, the spin wave contributions to the dressed energies of these phases can be analytically derived. In the same fashion as the SU(2) Fermi gas [39], we can express the TBA equations of the dressed energies as Here r α is viewed as the effective mass of the corresponding charges. For the unpaired atoms, r u = 1 and for the pairs, r p = 2. In the above equations, f α denotes the contributions of spin wave bound states to the dressed energies of the charge degree of freedoms. However, for h ≫ T , the spin contribution term f α will exponentially decay as e −h/T and it is negligible in low energy physics. Follow the method proposed in [39], for the strong interaction regime, we find the following equation of states Here the effective chemical potentials A α is given by and D (a) are constants resulted from the strong coupling expansion of the integral kernels K. F j (A α /T ) are Fermi-Dirac integrals which can also be written as polylogarithm functions [39]. In the FP phase, these constants are given by D (1) = −D (3) = 1. By iteration with p α and the effective chemical potentials A u,p , we obtain a close form of the equation of sates for the FFLO phase as where F α j = F j (A α /T ) and renormalized chemical potentials A u,p are given by This spin-3/2 Fermi gas exhibiting rich quantum phase transitions provides an ideal model to investigate quantum criticality of large spin interacting fermions. The phase diagram Fig. 2 has multifold critical points. In particular, the chemical potential and external field could drive the system from one phase to another as these parameter pass across the phase boundaries in the phase diagram Fig. 2 at T = 0. All phase transitions in Fig. 2 are of the second order. Near the critical point, thermodynamical properties evolve into certain universal scaling forms as temperature tends to zero. Below a crossover temperature (see the upper panel of Fig. 8 below) the quantum criticality of the Luttinger liquid gives the dynamical exponent z = 1. However, beyond the crossover temperature, a non-relativistic quadratic dispersion leads to the free fermion criticality with z = 2, where new band excitations are involved [24,25,39,53].

Conclusion
In conclusion, we have proposed a 1D integrable spin-3/2 fermionic gas of cold atoms with spin SO(4) symmetry. The symmetry, conserved quantities and integrability of this model has been studied by means of the BA We have shown that the integrable SO(4) symmetry spin-3/2 Fermi gas exhibits the spin singlet and quintet Cooper pairs in the two sets of SU(2) ⊗ SU(2) spin subspaces. The Bethe ansatz equations and finite temperature thermodynamical equations have been derived for the model in an analytical way. In particular, using the Bethe ansatz exact solutions we have thoroughly investigated spin pairing, full phase diagrams, equation of states, elementary excitations, correlation functions, magnetism and quantum criticality of the model. We have also shown that the SO(4) symmetry Fermi gase possesses various phases of Luttinger liquids with novel a magnetism. It is particular interesting that in the pure paired phase, breaking a 2-string of ν leads to four spin-1 spinons in spin excitations. The Luttinger liquid physics and and universal scaling behaviours of thermodynamical properties in regard of different spin states provide insights into understanding the large spin phenomena in the 1D interacting fermions.
Moreover, long distance asymptotics of various correlation functions have been calculated by using conformal field theory. In particular, the FFLOlike pair-pair correlations exist in the mixed phase of spin-2 pairs and single spin-3/2 atoms. Furthermore, the density profiles of the trapped gas has been also discussed in the context of a experimental setting with ultracold atoms. It turns out that the trapping center can be a mixture of different types of spin pairs and unpaired single atoms accompanied by the multiple shells involving the fully-polarized spin-2 pairs, the FFLO state, and the FF phase. Our results open to further study of the non-SU(κ) symmetry interacting fermions of ultracold atoms in theory and experiment.