Universal thermodynamics of the one-dimensional attractive Hubbard model
Song Cheng Yi-Song Yu Murray T Batchelor Xi-Wen Guan




State Key Laboratory of Magnetic Resonance and Atomic and Molecular Physics,
Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, China
University of Chinese Academy of Sciences, Beijing 100049, China.
Department of Theoretical Physics, Research School of Physics and Engineering,
Australian National University, Canberra ACT 0200, Australia


State Key Laboratory of Magnetic Resonance and Atomic and Molecular Physics,
Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, China
University of Chinese Academy of Sciences, Beijing 100049, China.



Centre for Modern Physics, Chongqing University, Chongqing 400044, China
Department of Theoretical Physics,
Research School of Physics and Engineering,
Australian National University, Canberra ACT 0200, Australia
Mathematical Sciences Institute,
Australian National University, Canberra ACT 0200, Australia




[]xiwen.guan@anu.edu.au
State Key Laboratory of Magnetic Resonance and Atomic and Molecular Physics,
Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, China
Department of Theoretical Physics,
Research School of Physics and Engineering,
Australian National University, Canberra ACT 0200, Australia

Center for Cold Atom Physics, Chinese Academy of Sciences, Wuhan 430071, China






71.10.Fd, 75.40.Cx,02.30.Ik




The one-dimensional (1D) Hubbard model, describing electrons on a lattice with an on-site repulsive interaction,
provides a paradigm for the physics of quantum many-body phenomena.




Here by solving the thermodynamic Bethe ansatz equations we study the universal thermodynamics, quantum criticality and  magnetism
of the 1D attractive Hubbard model.

We show  that the compressibility and the susceptibility  of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO)-like state obey simple additivity rules at low temperatures,
indicating an existence of two  free quantum fluids.

The magnetic properties, such as magnetization and susceptibility, reveal three physical regions: quantum fluids  at low temperatures,
a non-Fermi liquid at high temperatures and the quantum fluid to non-Fermi liquid crossover in between.

The lattice interaction is seen to significantly influence the nature of the FFLO-like state in 1D.

Furthermore, we show  that the dimensionless Wilson ratio provides an ideal parameter to map out the various phase boundaries and to
characterize the two  free  fluids   of the FLLO-like state.

The quantum scaling functions for the thermal and magnetic properties yield the same dynamic critical exponent z=2 and
correlation critical exponent ν=1/2 in the quantum critical region whenever a phase transition occurs.

Our results provide a rigorous understanding of quantum criticality and free fluids of many-body systems on a 1D lattice.


 

Universal thermodynamics of the one-dimensional attractive Hubbard model
    Xi-Wen Guan
    December 30, 2023
========================================================================







§ INTRODUCTION


How to capture the essential features of many-body physics through a simple model is always of great importance in condensed matter physics.

In this regard,  the Hubbard model <cit.> has long provided an active area of research since it was put forward
as an instance of a Mott insulator and later considered as a potential high-T_c superconductor.
 
 The Hubbard model has thus become a prototypical strongly correlated system which provides rich many-body phenomena,  such as a Mott transition,
 superconductivity, spin-charge separation and a Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state.
 
However, the Hubbard model, as a simplification of interacting fermions on realistic lattices, can be analytically resolved in neither two-dimensions (2D) nor three-dimensions (3D).

The one-dimensional (1D) case within a single band is integrable, firstly solved by Lieb and Wu in terms of  the
Yang-Baxter equation <cit.> and the nested Bethe ansatz <cit.> (see Ref. Ess05 for an extensive review).

More specifically, since Lieb and Wu's seminal work, the 1D repulsive Hubbard model has been investigated in various aspects,
including, but not restricted to, thermodynamic properties in the ground state <cit.>,
low-lying excitations <cit.>,
finite temperature thermodynamics <cit.>
and correlation functions <cit.>.


The thermodynamics of the 1D Hubbard model  is accessible through two alternative approaches – the thermodynamic Bethe ansatz (TBA) equations <cit.>
and the quantum transfer matrix method <cit.>.

The former is established on the so-called `string  hypothesis' and Yang-Yang grand canonical ensemble approach <cit.>,
whereas  the latter  stems from the lattice path integral formulations for the partition function <cit.>.

 In principle, the low-lying excitations can be constructed with the help of the TBA equations in the zero temperature limit and by the logarithm
 of the Lieb-Wu equations <cit.>.
 
Despite these systematic approaches and other methods  employed for the study of the
ground state properties <cit.>
and low-lying excitations <cit.>,
a complete understanding of the universal thermodynamics and quantum criticality of the 1D Hubbard model has not yet been achieved.
 
The key reason for preventing the solution of this problem is the  difficulty of finding  a suitable generating function for the equation of state at low temperatures.


On the other hand, the correlation functions are also extremely difficult to calculate  directly using the Bethe wave function.

For  a 1D conformally invariant  system, the critical exponents determining the power law decay of correlation functions
are connected with finite-size corrections to the ground state energy <cit.>.

The 1D repulsive Hubbard model is conformally invariant only in the vicinity of Fermi points.

 The conformal field theory (CFT) approach provides one method to obtain the asymptotics of correlation functions <cit.>.

The low-lying excitations provide a practicable opportunity for  an investigation of long distance asymptotics of correlation functions <cit.>,
where the finite-size corrections are  accessible through the Bethe ansatz method <cit.>.

However, difficulties involved in the actual calculations of correlation functions usually  prevent full
access to the many-body correlations <cit.>.



The mechanism of Cooper pairing in the 1D attractive Hubbard model has attracted attention <cit.> due to the discovery of high-temperature superconductors.

 In particular, the FFLO-like pair correlation and spin correlations are consequently  investigated by various methods,
 such as density-matrix renormalization group <cit.>, quantum Monte Carlo <cit.> and CFT <cit.>.

 The nature of the FFLO-like pair correlation was predicted in expansion dynamics of the attractive Hubbard model  trapped in 1D <cit.>.
 
Very recently,  trapping cold atoms on  optical lattices becomes a promising  method to simulate the many-body physics of  the Hubbard model <cit.>.

In particular,  ultracold  atoms offer an ideal platform for testing results predicted from  1D exactly solvable  models <cit.>.


It is understood that the macroscopic behaviour of 1D materials, such as  the spin compound Cu(C_4H_4N_2)(NO_3)_2 <cit.>
and the heavy fermion material YbNi_4P_2 <cit.>, demonstrates a type of 3D Fermi liquid behaviour <cit.>.

The motivation of the present work is to provide understanding of free fluid  nature and quantum criticality in the context of the 1D attractive Hubbard model.

Firstly, the 1D attractive Hubbard model plays an important role in understanding many-body phenomena such as superconductivity,
BEC-BCS crossover and FFLO-like  correlation <cit.>,
with several publications touching upon it <cit.>.

 Secondly, one expects to find universal behaviour for this model, including thermodynamics, quantum criticality and Luttinger liquid properties.
 
 Thirdly, regarding  the complicated FFLO state, it is highly desirable to obtain simple rules to describe the nature of quantum liquids  in the attractive Hubbard model.
 
 Last, but not least, the interplay of this work with experiments with ultracold atoms <cit.> may broaden our knowledge of many-body physics through 1D exactly solvable models.


This paper is organized as follows. In section II, we present  a derivation of the TBA equations for the 1D attractive Hubbard model and
determine the ground state phase diagram.

In section III, we derive the equation of state in the strong coupling regime.

In section IV, using  the equation of state, we obtain various analytical results for the thermodynamics and magnetism which are relevant to experimental study.

We also investigate quantum criticality and obtain the universal scaling forms of thermodynamic quantities.

In section V, we demonstrate  the free fluid nature of the FFLO phase through  the simple additivity rules of  the thermodynamic quantities.

We find that  the compressibility Wilson ratio is  very powerful in identifying the Fermi liquid/Tomonaga-Luttinger liquid phases in the low temperature phase diagram.

The last section VI is reserved for  a summary  and conclusion.

We conclude this section by noting that this article provides a fuller and more detailed account of our key results presented elsewhere <cit.>.





§ THERMODYNAMICS: THE YANG-YANG APPROACH




 §.§ Thermodynamic Bethe Ansatz equations



The 1D Hubbard model is described by the Hamiltonian

    H=   -∑_j=1^L ∑_a=↑,↓( c_j,a^† c_j+1,a + c_j+1,a^† c_j,a) 
       + u ∑_j=1^L ( 1-2n_j,↑) ( 1-2n_j,↓),

where c_j,a^† and c_j,a are the creation and annihilation operators of fermions  with spin a (a=↑ or a=↓) at site j in a
1D  periodic lattice of length L, n_j,a=c_j,a^† c_j,a is the corresponding particle number operator,
and u represents an  on-site interaction between particles (u>0 for repulsion and u<0 for attraction).

By means of the  Bethe ansatz, the eigenenergies of the Hamiltonian are given by E=-2∑_j=1^N cos k_j+u(L-2N), where
the  quasimomenta  {k_j}  satisfy the Lieb-Wu equations <cit.>

    exp(i  k_j  L)   =∏_α=1^M sin k_j-Λ_α+i  u/sin k_j-Λ_α - i  u, 
    ∏_j=1^N sin k_j-Λ_β + i  u/sin k_j-Λ_β - i  u   =-∏_α=1^M Λ_α-Λ_β+2  i   u/Λ_α-Λ_β-2  i  u,

where {Λ_β} denote spin rapidities, j=1,2,…,N, β=1,…,M, with N and M the total particle number and spin down particle number, respectively.





Similar to the analysis<cit.> used for  the repulsive case u>0,
one finds that the roots to the Bethe ansatz equations (<ref>) and (<ref>)
for  the attractive Hubbard model  can be  divided into three categories: single real k, k-Λ string and Λ-Λ string,
which constitute the  string hypothesis.

They are given by <cit.>


  
  * single real k's.
 
  * the α-th k-Λ string of length-m, for which there are  2m k's,

    k_α^1   =arcsin(Λ_α^'^m+ i  m   |u|),
    
    k_α^2   =arcsin(Λ_α^'^m+ i   (m-2)    |u|),
    
    k_α^3   =π-k_α^2,
    ⋮
    
    k_α^2m-2   =arcsin(Λ_α^'^m- i (m-2)  |u|),
    
    k_α^2m   =arcsin(Λ_α^'^m-i  m   |u|),

accompanied by m spin-rapidities

    Λ_α^'^m,j   =   Λ_α^'^m+ i  (m+1-2j)    |u|,

 in  Λ space, where  j=1,2,3,…, m and Λ_α^'^m is the real center of the k-Λ string, see Fig. <ref>.
 
  
  * the β-th Λ-Λ string of length-m,
  
    Λ_β^m,j=Λ_β^m+ i  (m+1-2j)  |u|,

where j=1,2,3, …, m, and Λ_α^m is the real center of the Λ string.

The Λ strings represent the spin wave bound states in the spin sector.



In the above equations we  denoted M_m, M_m^', and ℳ_e as the number of Λ strings of length m,
of k-Λ strings of length-m, and of single real k's, respectively.

It is easy to see that M=∑_m=1^∞ m(M_m+M_m^') and N=ℳ_e+∑_m=1^∞ 2m M_m^'.


Substituting the string hypothesis into the Lieb-Wu equations and taking logarithms leads to the discrete nested BA equations


    k_j  L =2π  I_j+∑_m=1^∞∑_α=1^M_m^'θ( sin k_j-Λ_α^'^m/m |u|)
    +∑_m=1^∞∑_α=1^M_mθ( sin k_j-Λ_α^m/m |u|),
          ∑_j=1^N-2M^'θ( Λ_α^n-sin k_j/n |u|)
    =2π  J_α^n+∑_m=1^∞∑_β=1^M_mΘ_n m( Λ_α^n-Λ_β^m/n |u|),
          
    2L [ arcsin(Λ_α^'^n+i  n   |u|) ]
    =2π J_α^'^n+∑_j=1^N-2M^'θ( Λ_α^'^n-sin k_j/n |u|)
    +∑_m=1^∞∑_β=1^M_m^'Θ_nm( Λ_α^'^n-Λ_β^'^m/|u|),


where M^'=∑_m=1^∞ m M_m^' is the total number of Λ's involved in the k-Λ strings, θ(x)=2 arctan(x), and


    Θ_nm(x)=θ(x/|n-m|) + 2θ( x/|n-m|+2)+⋯+  2θ( x/n+m-2) + θ( x/n+m)    if n≠ m
    
    2θ( x/2)+2θ( x/4)+⋯+  2θ( x/2n-2)+θ(x/2n)    if n=m.



The quantum numbers I_j, J_α^n and J_α^'^n are either integers or half-odd integers, stemming from the multivaluedness of
the log functions. They are determined by the  relations

    I_j=   integers   if ∑_m=1^∞ (M_m^'+M_m) is even
    half-odd integers   if ∑_m=1^∞ (M_m^'+M_m) is odd,
    
    J_α^n=   integers   if N-M_n is odd
    half-odd integers   if N-M_n is even,
    J_α^'^n=   integers   if L-N+M_n^' is odd
    half-odd integers   if L-N+M_n^' is even.



With the help of the string hypothesis, the eigenenergies are

    E   =   -2∑_j=1^N-2M^'cos k_j-2uN+uL 
          -4∑_n=1^∞∑_α=1^M_n^'[ √(1-(Λ_α^'^n+i   n   |u|)^2)].


We now introduce counting functions for the quantum numbers,
y(k_j)=2π  I_j/L, z_n(Λ_α^n)=2π  J_α^n/L and z_n^'(Λ_α^'^n)=2π J_α^'^n/L.

Considering the thermodynamic limit, N, M, L→∞ with N/L,  M/L finite, we further define the distributions

    y(k)/k   =   
    2π[ρ^p(k)+ρ^h(k) ],
    z_n(Λ)/Λ   =   
    2π[σ_n^p(Λ)+σ_n^h(Λ)],
    z_n^'(Λ)/Λ   =   
    2π[σ_n^'^p(Λ)+σ_n^'^h(Λ)],

where ρ^p, σ_n^p, σ_n^'^p (ρ^h, σ_n^h, σ_n^'^h) are root densities of particles (holes) in quasimomenta of excess fermions,
Λ-string parameter space  and k-Λ string space, respectively.

Then one can derive the densities of excess fermions, Λ-spin strings and k-Λ strings, with

    ρ^p(k)+ρ^h(k)=1/2π
          -cos k∑_n=1^∞∫_-∞^∞Λ  a_n(sin k-Λ) [ σ_n^p(Λ)+σ_n^'^p(Λ) ], 
          σ_n^h(Λ)= ∫_-π^πk   a_n(sin k-Λ) ρ^p(k)
          -∑_m=1^∞ A_nm∗σ_m^p(Λ),
          σ_n^'^h(Λ) =1/π[1/√(1-(Λ+i  n  |u|)^2)] 
          -∑_m=1^∞ A_nm∗σ_m^'^p (Λ) -∫_-π^πk   a_n(sin k-Λ)ρ^p(k),

where the function

    a_n(x)=1/2π2n|u|/(n|u|)^2+x^2.

 As usual, ∗ stands for the convolution (f∗ g)(Λ)=∫_-∞^∞f(Λ-Λ')g(Λ')dΛ', namely,

    A_nm∗ f (x)=δ_n, m  f(x)+∫_-∞^∞y/2π/xΘ_nm( x-y/|u|)f(y).

Here we denoted the derivative of the function  Θ_nm as

    1/2π/xΘ_nm(x)
    =
         a_|n-m|(x)+2a_|n-m|+2(x)
    
         +…+a_n+m(x)     if n≠ m
    
         2a_2(x)+2a_4(x)+…
    
         +2a_2n-2(x)+a_2n(x)    if n=m.


The root  distribution functions (<ref>)-(<ref>) determine spin and charge excitations, spin dynamics and full energy spectra.
In  the grand canonical ensemble, the Gibbs free energy per site  can be expressed in terms of these root densities in different sectors

    f   =     e-μ  n_c-2B   m-T  s 
       =   ∫_-π^πk   (-2cos k-μ-2u-B)ρ^p(k)
           -∑_n=1^∞∫_-∞^∞Λσ_n^'^p(Λ)
    [ 4√(1-(Λ_α^'^n+ i  n   |u|)^2). 
          . +   n(2μ+4u) ] 
          +∑_n=1^∞∫_-∞^∞Λ  2n  B  σ_n^p(Λ)-T  s+u,

where μ is the chemical potential, B the magnetic field and T the temperature.

In the above equations n_c is the particle density, m=N-2M/2L the magnetization and s the entropy per site.


Following the Yang-Yang grand can­onical descrip­tion <cit.>, the entropy per site is explicitly given by


    s   =   ∫_-π^πk  {( ρ^p(k)+ρ^h(k))ln( ρ^p(k)+ρ^h(k)) -ρ^p(k)lnρ^p(k)-ρ^h(k)lnρ^h(k) }
           +∑_n=1^∞∫_-∞^∞Λ {(σ_n^'^p(Λ)+σ_n^'^h(Λ)) ln(σ_n^'^p(Λ)+σ_n^'^h(Λ)) -σ_n^'^p(Λ)lnσ_n^'^p(Λ) -σ_n^'^h(Λ)lnσ_n^'^h(Λ)}
           +∑_n=1^∞∫_-∞^∞Λ {( σ_n^p(Λ)+σ_n^h(Λ)) ln( σ_n^p(Λ)+σ_n^h(Λ)) -σ_n^p(Λ)lnσ_n^p(Λ) -σ_n^h(Λ) lnσ_n^h(Λ) }.



In the following, we only consider the physics with B≥ 0 and μ≤ 0.



In the  thermodynamic equilibrium, the true equilibrium state can be determined by the minimization of the free energy with respect to  the densities.

Carrying out a  variation of (<ref>) under the restriction  of  (<ref>)-(<ref>), we obtain  the TBA equations for the attractive Hubbard model
in the form













    ε^u(k)   =   -2cos k-μ-2u-B  
          +∑_n=1^∞∫_-∞^∞Λ  a_n(sin k-Λ)ε_n^' -(Λ)
          -∑_n=1^∞∫_-∞^∞Λ  a_n(sin k-Λ)ε_n^ -(Λ), 
    ε_n(Λ)   =   ∫_-π^πk  cos k   a_n(sin k-Λ)ε^u -(k)
          + 2nB+∑_m=1^∞ T_nm∗ε_m^ -(Λ),
    ε_n^'(Λ)   =   -4√(1-(Λ+ i   n  |u|)^2)-n(2μ+4u)
          +∫_-π^πk  cos k   a_n(sin k-Λ)ε^u -(k)
          +∑_m=1^∞ T_nm∗ε_m^' -(Λ),


where we have denoted

    ε^u -(x)    =    Tln(1+e^-ε^u(x)/T),
    ε_n^' -(x)    =    Tln(1+e^-ε_n^'(x)/T),
    ε_n^ -(x)    =   Tln(1+e^-ε_n(x)/T).


In the above equations,  we  defined the dressed energies

    ε^u(k)   =   Tlnζ(k)=Tlnρ^h(k)/ρ^p(k),
    ε_n(Λ)   =   Tlnη_n(Λ)=Tlnσ_n^h(Λ)/σ_n^p(Λ),
    ε_n^'(Λ)   =   Tlnη_n^'(Λ)=Tlnσ_n^'^h(Λ)/σ_n^'^p(Λ).


The  convolution T_nm∗ f(x)=A_nm∗ f(x)-δ_n,m f(x) is defined by convention.



The TBA equations (<ref>)-(<ref>) indicate that the dressed energies
ε^u(k), ε_n(Λ), ε_n^'(Λ)
describe the excitation energies which are subject to interactions among the bound states of electrons, spin wave fluctuations, magnetic field and chemical potential.

 They contain   full  thermal and magnetic fluctuations in both spin and charge degrees of freedom.
 
 Therefore from these equations we can determine the thermal and magnetic properties of the model in full temperature regimes.
 
 After some algebra, the Gibbs free energy per site is consequently given by


    f=u-∫_-π^πk/2πε^u -(k)
          
    -∑_n=1^∞∫_-∞^∞Λ/π[1/√(1-(Λ+ i   n  |u|)^2)]ε_n^' -(Λ).


This result builds up analytical access to the full thermodynamics of the model.




 §.§ Zero Temperature Phase Diagram 


In  the zero temperature limit, most dressed energies are nonnegative and thus make no significant contributions to the free energy (<ref>).

We observe that  in the ground state, there exist only unpaired fermions and bound pairs of fermions.

The spin Λ-Λ strings Eq.(<ref>)  are suppressed due to the fact that in the FFLO-like phase IV, the spin wave bound states ferromagnetically couple to the Fermi sea of  the unpaired fermions.

The driving term in the TBA equation (<ref>) is positive due to this ferromagnetic ordering. At T→ 0, the  Λ-Λ strings  are gapped.

The driving term in the TBA equation (<ref>) can be positive when n≥ 2 due to the negative chemical potential.

Taking the  limit T→ 0, the corresponding TBA equations (<ref>) and (<ref>) thus reduce to coupled linear integral equations,
called the dressed energy equations,

    ε^u (k)=   -2cos k-μ-2u-B 
       -∫_-A^A Λ  a_1(sin k-Λ) ε_1^' (Λ),
    ε_1^' (Λ)=   -2μ-2∫_-π^πk  cos^2 k  a_1(sin k-Λ)
       -∫_-Q^Q k  cos k   a_1(sin k-Λ) ε^u (k) 
       - ∫_-A^A Λ^'  a_2(Λ-Λ^')ε_1^'(Λ^'),

where the integration boundaries Q and A represent the Fermi points of these two kinds of states (pairs and single fermions).

In Eq. (<ref>) we used the expression

    4√(1-(Λ-i  n |u|)^2)-4n|u|  
           =∫_-π^πk/πcos^2 k   2 n |u|/(nu)^2+(sin k-
    Λ)^2.


 The integration boundaries  are  determined by  ε^u (± Q)=0 and ε^'_1(± A)=0.
 
 Within the intervals [-Q,Q] and [-A,A], the dressed energies  are negative, i.e., ε^u(k)≤ 0 and ε^'_1(Λ)≤ 0.
 
 This means that   particle  states occupy all vacancies in the two Fermi seas.


With the help of (<ref>)-(<ref>), the root densities for quasimomentum  k and
spin rapidity Λ in the k-Λ string of length-1 at  zero temperature are expressed as

    ρ(k)=   1/2π-cos k ∫_-A^A Λ  a_1(sin k-Λ) σ_1^'(Λ),
    σ_1^'(Λ)=   1/π1/√(1-(Λ+ i   |u|)^2)-∫_-Q^Q k   a_1(sin k-Λ)ρ(k)
       -∫_-A^A Λ^'  a_2(Λ-Λ^')σ_1^'(Λ^').


In the grand canonical ensemble, we explicitly write down the above root densities,
which  satisfy the two conditions ∫_-Q^Q k  ρ(k) + 2∫_-A^A Λ σ_1^'(Λ)=N/L and
∫_-A^A Λ σ_1^' (Λ)=M/L=N_↓/L.

Thus the total particle density is given by n_c = N/L =∫_-Q^Q k  ρ(k) + 2∫_-A^A Λ σ_1^'(Λ)
 and the magnetization per site by m=(N-2M)/(2L)=1/2∫_-Q^Q k  ρ(k).

By varying the integration boundaries Q and A, the system possesses different fillings and quantum phases.

 A phase transition occurs when the dressed energies exactly satisfy ε^u(0)=0, ε^u(π)=0 or ε_1^'(0)=0.

Consequently we can determine  five phases, () vacuum, () fully polarized state, () half-filling state,
() partially polarized state, i.e., FFLO-like state, () fully paired state.

The phase boundary between () and () is determined by ε_1^'(0)=0  together with the condition  Q=0.

Then the TBA equation (<ref>) leads to the critical field value μ_c=2|u|-2√(1+u^2).

The phase boundary  between () and () and between () and ()
are determined by the conditions A=0, ε^u(0)=0 and by A=0, ε^u(π)=0, respectively.

With regard to the boundaries for the FFLO-like phase, the situation is much more  subtle.

The phase boundary between () and () is determined by ε_1^'(0)=0 and ε^u(Q)=0,
while the phase boundary between () and () is determined by ε^u(0)=0 and ε_1^'(A)=0.


The phase boundaries in the ground state phase diagram Fig. <ref> are summarized as follows


  * (-)

    μ_c1=2|u|-2√(1+u^2).


  * (-)

    μ_c2=-B-2u-2.


  * (-)

    μ_c3 = 2-B-2u.


  * (-)

    μ_c4=   2|u|-2√(1+u^2)
       -∫_-Q^Q k  cos k   a_1(sin k)[cos Q-cos k],
    
    B_c4=   2√(1+u^2)-2cos Q 
       -∫_-Q^Q k  cos k   a_1(sin k)[cos Q-cos k],

with Q∈ [0,π].


  * (IV-V) This phase transition occurs if the critical magnetic field is sufficient to break the bound state of fermions, whose boundary in principle is fixed by

    ε_1^' (Λ)=   -2μ-2∫_-π^πk  cos^2 k  a_1(sin k-Λ) 
       -∫_-A^A Λ^'  a_2(Λ-Λ^')ε_1^'(Λ^'), 
    ε_1^'(A)=     0, 
    μ=    -2-2u-B-∫_-A^A Λ  a_1(Λ) ε_1^' (Λ).


 When  A ≪ 1,  the density for pairs of fermions is low, the phase boundary could be obtained  by iteration,
 i.e., by applying Taylor expansion to (<ref>) with respect to Λ, it can be approximately resolved by iteration.

The  solution of (<ref>) gives  A in terms of μ and B, then we derive the phase boundary by
substituting the above results for A into the Eqs. (<ref>) and  (<ref>). By iteration, we finally obtain

    μ_c5≈     2|u|-B-2 
       +4√(2)/π |u|α_1[μ_c5+2(√(1+u^2)-|u|)]^3/2.

Here at low energy physics only length-1 k-Λ strings  are involved.

From the  TBA equations (<ref>)-(<ref>),  we may introduce the  parameters α_n  and β_n to indicate
the interacting effect of the length-n k-Λ bound states on a lattice in the low density regiem.

They are given by

    α_n    =   ∫_-π^π k    2n |u| cos^2 k (n^2u^2-3 sin^2 k)/π(n^2u^2+sin^2 k)^3,
    β_n   =   ∫_-π^π k   a_n(sin k).

In general, α_n represents the lattice effect in the length-n k-Λ strings.



Meanwhile, if A ≫ 1, the phase boundary is given  by (<ref>) and (<ref>)
in Appendix (<ref>), where we have used the Wiener-Hopf method to solve the TBA integral equations.



From the dressed energy equations (<ref>) and (<ref>) the complete phase diagram at zero temperature is shown in Fig. <ref>.

This phase diagram was also obtained by the Shiba transformation, which  builds  up a mapping between repulsive and attractive regions in the ground state
of the Hubbard model <cit.>.

However, once we are concerned with the low temperature thermodynamics, correlation functions and quantum criticality,
the Shiba transformation does not work in actual calculations,
see the analysis of the ground state properties of the attractive Hubbard model <cit.>.

This is mainly because   the different  spin-spin strings, k-Λ  strings and excess fermions  have different  cut-off  processes (the cut-off strings, see Appendix B) at low temperature physics.

For example, the spin fluctuation term (the third term)  in the unpaired dressed energy can be safely  ignored in the strongly attractive Hubbard model at low temperatures.

However, the counterpart of such a spin  fluctuation term in  the repulsive Hubbard model essentially determines the  antiferromagnetic ordering.

Even in the repulsive regime, such spin string dynamics, quantum criticality and scaling functions still lack an analytical calculation.

The ground state  properties of the attractive Hubbard model were initially studied by Woynarovich <cit.>.

In this paper, using the TBA equations (<ref>)-(<ref>), we obtain exact results for the FFLO pairing correlation,
universal thermodynamics and quantum criticality of the 1D attractive Hubbard model.

Our study provides a precise  understanding of the universal  low energy physics of interacting fermions with pairing and depairing on  a 1D lattice.





§ EQUATION OF STATE


The TBA equations describe the full  thermodynamics of the model.

At low temperatures quantum liquid behavior and critical scaling in the thermodynamics should be obtained from the TBA equations (<ref>)-(<ref>).

However, the analysis of such coupled nonlinear integral equations provides a formidable challenge.

In particular, it is challenging to solve infinitely many coupled nonlinear integral TBA equations, i.e., the desired analytical or numerical solution is not achievable by solving the whole set of TBA equations.
 
 This obstacle prevents us to understand the microscopic Cooper pairing mechanism and many-body phenomena for this model.
 
 On the other hand,  in the FFLO-like phase IV, the spin wave bound states ferromagnetically couple to the Fermi sea of  the unpaired fermions.
 In this phase, except two gapless excitations in the sectors of bound pairs and excess fermions, there exist a spin wave  ferromagnetic fluctuation, which is no longer a linear dispersion.
 
 Bosonization or Tomonaga-Luttinger liquid (TLL) theory <cit.> are not available once such a  ferromagnetic ordering is involved in the low temperature physics.
 
 A similar situation was studied in the 1D two-component Bose gas  <cit.>.


 Moreover, the TLL  is not  applicable to the quantum critical region near a phase transition.
 
Here we proceed with an analytical  investigation of the low energy physics of the 1D attractive Hubbard model
beyond the scope of the TLL approaches.
 
In order to obtain the universal thermodynamics and quantum criticality of the 1D attractive Hubbard model,
we first solve the TBA equations  (<ref>)-(<ref>) analytically in the strong coupling regime.

We will  derive the equation of state which is crucial for the investigation of the quantum criticality of the model.

These results can be helpful to understand current experimental developments in ultra-cold atoms <cit.>.


In the following discussion we  mainly concentrate on the low density regime.

In general, it is very difficult to find universal  characteristics of  quantum liquids in quantum many-body systems,
for example, for  the Gaudin-Yang Fermi gas <cit.>.

Under the assumption that the density of  pairs and the bound states of multiple fermions are low and the interaction is strong,
the TBA equations (<ref>)-(<ref>)  can be rewritten as

    ε^u(k)=    -2cos k + 2a̅cos^2 k - μ -2u -B +∑_n=1^∞ p_n^b 
       + a̅ - T e^-2B/T e^-K̅ I_0(K̅)+o(1/|u|^4),
    ε_n^'(Λ)=   -2nμ+η_n-d_1/π n |u| - d_2/π (n|u|)^3
       +Λ^2 [ d_1/π(n|u|)^3-φ_n ] +o(1/|u|^4),





where K̅=∫_-π^πk/2π cos k  ln(1+e^-ε^u(k)/T)
and I_0(x) is the zeroth order modified Bessel function, which stems  from the spin-wave contributions.

In the above equations, we denoted

    d_1   =   2π-∫_-π^πk  cos k  ε^u -(k),
    
     d_2   =   -π/2-∫_-π^πk  cos k  sin^2 k  ε^ u-(k),
    η_n   =   ∑_m=1^∞∫_-∞^∞Λ  T_nm(Λ)ε_m^' -(Λ),
    a̅   =   1/2∑_n=1^∞∫_-∞^∞Λ[b_n(Λ)-4 b_n(Λ)Λ^2/(n u)^2+Λ^2]
      ε_n^'  -(Λ),
    φ_n   =   ∑_m=1^∞∫_-∞^∞Λ  Q_nm(Λ)ε_m^'  -(Λ),

 with  b_n(Λ)=a_n(Λ)/(n u)^2+Λ^2 and

    Q_nm(x)=
    
    b_|n-m|(x)+2b_|n-m|+2(x)+…
    
    +2b_n+m-2(x)+b_n+m(x)   if n ≠ m
    
    2b_2(x)+2b_4(x)+…
    
    +2b_2n-2(x)+b_2n(x)     if n=m.

The results (<ref>) and (<ref>) are valid for the low density limit and strong interaction regime.


Substituting <ref> into (<ref>), the pressure per unit length is given by  p=p^u+∑_n=1^∞ p_n^b+|u| with


    p^u=T∫_-π^πk/2πln(1+ e^-ε^u(k)/T),   
    p_n^b=T∫_-∞^∞Λ/π[1/√(1-(Λ+ i  n  |u|)^2)]ln(1+e^-ε_n^'(Λ)/T).



Using the results (<ref>) and (<ref>) and taking integration by parts within the above expressions for the effective pressures (<ref>),
we then obtain the set of coupled equations

    p^u   =     T ln( 1+e^(μ+2u+B-∑_n=1^∞ p_n^b + a̅ -2 )/T)  
           -2a̅/π∫_-1^1 x  x^2/√(1-x^2)/1+e^2x/T/z +2a̅/1+e^4/T/z^2
          +2/π∫_-1^1 x  arccos(-x)/1+e^2x/T/z +o(1/u^4), 
    
    p_n^b   =     T[1-1/4(nu)^2]ln(1+e^2nμ/T) 
          +2 D_n/π[1-1/4(nu)^2] ∫_0^∞x  arctan√(x)/1+e^D_n x/T/ζ_n
          +o(1/u^4),


 which serve as the equations of state.
 
In the above equations,

    D_n   =   d_1/π n|u|-(n u)^2φ_n,
    
    z   =   e^( μ+2u+B-∑_n=1^∞ p_n^b +a̅)/T,
    ζ_n   =   e^( 2nμ-η_n+d_1/π n |u|+d_2/π(n|u|)^3) /T.


We also  defined  the auxiliary functions

    d_1   =    2π-4 ∫_-1^1 x  √(1-x^2)/1+e^2x/T/z
          - 4a̅∫_-1^1 x  x^3/√(1-x^2)/1+e^2x/T/z+o(1/u^4),
    
    d_2   =     -π/2 - 4/3∫_-1^1 x  (1-x^2)^3/2/1+e^2x/T/z+o(1/u^4),
    a̅   =    ∑_n=1^∞D_n/π (n u)^2∫_0^∞x  √(x)/(1+x)^2/1+e^D_n x/T/ζ_n+o(1/u^4),
    η_n   =   ∑_m=1^∞𝔗_nm^ξ (m)+o(1/u^4), 
    φ_n   =   ∑_m=1^∞𝔗_nm^ϕ (m)+o(1/u^6).

In these equations, we define 𝔗_nm^x (m)= x_|n-m|^m+2x_|n-m|+2^m + ⋯ + 2x_n+m-2^m + x_n+m^m, with x_0^m=0 (x=η,ϕ) and  auxiliary functions

    ξ_p^m=     Tln(1+e^2mμ/T) 
       + 2 D_m/π∫_0^∞x  arctan( m/p√(x))/1+e^D_m x/T/ζ_m,
    ϕ_p^m=    T/2(p u)^2(1+e^2mμ/T) 
       + m/pD_m/π u^2∫_0^∞x  √(x)/(p^2+m^2 x)/1+e^D_m x/T/ζ_m
       + D_m/π (p u)^2∫_0^∞x  arctan( m/p√(x))/1+e^D_m x/T/ζ_m.


These functions are indicative of the  sophisticated many-body effects induced by k-Λ strings of different lengths.

A more detailed derivation of the above result is presented in Appendix <ref>.


In order to conceive the universal behavior of the system, we need to further simplify the equations of state (<ref>) and  (<ref>).

 To this end, we utilize the conditions |μ/T| ≫ 1 and strong interaction |u|≫ 1, which suppress the large length  k-Λ strings in  this physical regime.
 
 We observe  that  no larger length-n  k-Λ bound states than n=1 exist in the FFLO  phase IV at low temperatures.
 
 Then the pressure per unit length  simplifies to p= p^u + p^b + |u|, where  p^u and p^b are given by

    p^u=     Tln(1+e^(μ+2u+B-p^b-2)/T) 
       +2/π∫_-1^1 x  arccos(-x)/1+e^2x/T/z_1+o(1/u^2), 
    
    p^b=    2 D_1/π∫_0^∞x  arctan√(x)/1+e^D_1 x/T/ζ+o(1/u^2),

where z_1=e^( μ+2u+B-p^b )/T, ζ=e^( 2μ-η+d_1/π |u|)/T
and the above  auxiliary functions with n=1 read

    D_1   =   d_1/π |u|-u^2φ,
    
    d_1   =    2π-4 ∫_-1^1 x  √(1-x^2)/1+e^2x/T/z_1+o(1/u^2),
    η   =    2 D_1/π∫_0^∞x  arctan( 1/2√(x))/1+e^D_1 x/T/ζ+o(1/u^2), 
    φ   =   D_1/2π u^2∫_0^∞x  √(x)/(4+x)/1+e^D_1 x/T/ζ
          + D_1/4π u^2∫_0^∞x  arctan( 1/2√(x))/1+e^D_1 x/T/ζ+o(1/u^4).


Here we only consider the corrections up to order 1/|u| in the strong coupling regime |u|≫ 1.

The equations of state (<ref>) and (<ref>)  give a very good approximation of the low energy physics.

In Fig. <ref>, we demonstrate the  accuracy of these  equations  compared to the numerical results obtained from the TBA equations (<ref>)-(<ref>).

The peaks in the susceptibility and the discontinuities of the first derivative  of  the density reveal important behavior of the model near quantum phase transitions.








The pressures (<ref>) and (<ref>) could be further approximately resolved by appropriate iteration.

For the low density regime n_c≪ 1,  we  expand the numerators  in the pressure p^b
and the auxiliary functions η and φ  with respect to a small value of   x in these integrals.
Then we can represent η and φ in terms of p^b.

After iteration we  thus obtain  p^b≈ -T^3/2  f_3/2/√(d_1/|u|-π p^b/8),
where we have defined f_s=_s [-exp( 1/T( 2μ-p^b/2+d_1/π |u|) )]
in terms of the polylog function _s(x).

Using this expression for p^b and after some lengthy algebra,  we finally obtain the closed form expressions

    p^b=   -1/√(π D_0)  T^3/2f̃_3/2+o(1/u^2,T^2 ),
    
    p^u=    Tln( 1+e^(μ+2u+B- p^b-2)/T)
       +2/π∫_-1^1 x  arccos(-x)/1+e^2x/T/(z_0  e^-p^b/T)+o(1/u^2,T^2 )

for the two pressures, with the auxiliary functions

    d_0=   2π-4 ∫_-1^1 x  √(1-x^2)/1+e^2x/T/z_0+o(1/u^2,T^2 ),
    
    D_0=   d_0/|u|π+1/8√(|u|/d_0)  T^3/2g_3/2+o(1/u^2,T^2 ).


In results (<ref>) and (<ref>),
f̃_s=g_s-1/2√(|u|/d_0)  T^1/2 g_s  g_s-1, z_0=e^(μ+2u+B)/T
and g_s=Li_s[ -e^(2μ+d_0/|u|π)/T].

The pressures (<ref>) and (<ref>) give deep insight into quantum scaling in the critical regimes.




§ QUANTUM CRITICALITY


 Quantum phase transitions occur in the attractive Hubbard model at zero temperature as the external magnetic field and chemical potential
 are varied across any phase boundary in Fig. <ref>.
 
In general, near a quantum critical point, the model  is expected to show universal scaling behaviour in the thermodynamic quantities due to the collective nature of many-body effects <cit.>.

We see that the 1D attractive Hubbard model is an ideal model to explore such a universal  scale-invariant description on a 1D lattice,
which  can be determined by the power-law scaling of the various thermodynamic properties.

The behavior of the thermodynamic quantities  is governed by scaling functions with critical exponents
in the V-shaped region fanning out to  finite temperatures from the quantum critical point.
 
 In order to calculate the thermodynamic quantities which contain enough thermal and quantum fluctuations to describe quantum criticality,
 we here use the form of the equation of state with the results given in (<ref>), (<ref>), and (<ref>)-(<ref>) for the pressure terms.
 
 We observe that  first-order derivatives of these pressures  with respect to μ or B form a set of linear equations.
 
 Solution to this set of linear equations directly leads  to the particle density n_c=(∂ p/∂μ)_B and magnetization m=1/2(∂ p/∂ B)_μ.
 
 Similarly, one can  derive the second-order derivatives of the pressures, the compressibility κ=(∂ n/∂μ)_B and  the susceptibility χ=(∂ m/∂ B)_μ.
 
The corresponding scaling laws can be obtained in the different physical regimes.


At very low temperatures, spin fluctuation in the FFLO-like phase is suppressed,
as are the bound states of higher k-Λ strings for |μ/T| ≫ 1.

In this regime, the thermodynamics of the model  is governed by a two-component TLL
or say two-component Fermi liquid consisting   of excess fermions and of hard-core bosonic  charge bound states.

 The leading low-temperature correction to the free energy is given by

    f ≈ f_0 - π T^2/6( 1/v_1 + 1/v_2),

where f_0 is the ground state free energy and v_1 (v_2) is the sound velocity of excess fermions (bound pairs).
This result is valid for arbitrary interaction strength.

When  the particle density is very low, i.e., n_1,2≪ 1,  we   explicitly obtain the two velocities

    v_1≈     2π n_1 [1+4 n_2/|u|+12(n_2/|u|)^2  ], 
    
    v_2≈    π  n_2 √(2α_1)/β_1[1+1/β_1(2 n_1/|u|+n_2/|u|) . 
       . + 3/β_1^2(2 n_1/|u|+n_2/|u|)^2 ],

where the lattice  parameters

    α_1   =   ∫_-π^πk  2|u|cos^2 k   (u^2-3sin^2 k)/π(sin^2 k+ u^2)^3,
    β_1   =   ∫_-π^πk  a_1(sin k),

 are  functions of |u| representing the lattice effect <cit.>.
 
 In the above equations, n_1,2 stands for the densities of excess fermions and the bound pairs, respectively.

We plot the two lattice parameters against interaction strength in Fig. <ref>.

We shall see that the  critical exponents and thermodynamics of the model are subject to these two   parameters.

The susceptibility is independent of temperature so that the dimensionless Wilson ratio reaches a constant (we will study this nature of the Fermi liquid in the next section).

The TLL validates only in the region below the crossover temperatures, where the entropy or specific heat retains a linear temperature-dependence, see the dashed lines in Fig. <ref>.

The entropy in the temperature-magnetic field plane displays the visible areas of the critical regions (QC) near different critical points.

In what follows, we will derive the scaling functions for the critical regions.







Using the equation of state with the pressures (<ref>), (<ref>), and (<ref>)-(<ref>),
we can further derive the scaling forms of the thermodynamic quantities in the  critical regimes.
 
Analytic results for the scaling functions help to understand the  microscopic origin of quantum criticality of  the 1D attractive Hubbard model.
 
For convenience, we first simplify the auxiliary functions

    δ=   1/π∫_-1^1 x  1/1+e^2x/T/z̅_0x/√(1-x^2),
    γ=   1/π∫_-1^1 x  1/1+e^2x/T/z̅_01/√(1-x^2),
    γ'=   1/π T∫_-1^1 x  e^2x/T/z̅_0/(1+e^2x/T/z̅_0)^21/√(1-x^2),
    δ'=   1/π T∫_-1^1 x  e^2x/T/z̅_0/(1+e^2x/T/z̅_0)^2x/√(1-x^2),

where z̅_0=exp( μ+2u+B+ T^3/2 g_3/2 / √(π D_0)) with D_0 given in (<ref>).

By virtue of results (<ref>), (<ref>), and (<ref>)-(<ref>), the closed form expressions for thermodynamic quantities
can be derived.

For strong attraction, we have the relations

    n_c=    γ+(1-γ)∂ p^b/∂μ,
    m=    1/2[ γ+(1-γ)∂ p^b/∂ B],
    κ=    (1- ∂ p^b/∂μ)^2 γ' + (1-γ) ∂^2 p^b/∂μ^2,
    χ=    1/2[ (1- ∂ p^b/∂ B)^2 γ' + (1-γ) ∂^2 p^b/∂ B^2],

for thermodynamic quantities. Here we calculated the derivatives of the pressures









    ∂ p^b/∂μ   =    -τ^1/2(2 f_1/2-τ f_3/2)/Δ_t ,
    ∂ p^b/∂ B   =    -2 δτ^1/2(f_1/2 - τ f_3/2)/|u|Δ_t ,
    ∂^2 p^b/∂μ^2   =    - f_-1/2(16 π - √(π)τ^3/2 f_3/2)/4D_0 τ^1/2Δ_t^3,
    ∂^2 p^b/∂ B^2   =    -2δ' τ^1/2(f_1/2-τ f_3/2)/|u|Δ_t
    -πδ^2 f_-1/2/u^2 D_1 τ^1/2(4 f_1/2-τ f_3/2) Δ _t^4
          ×( 16 √(π)f_1/2 -8 τ^1/2 f_1/2^2 -5  τ ^3/2 f_1/2 f_3/2 -4 √(π)τ  f_3/2)

with
τ=T/D_0 and
Δ_t=√(π)- 1/2τ^1/2 f_1/2+1/8τ^3/2 f_3/2.
These results constitute very accurate results for the thermodynamics.

The asymptotic results for the thermodynamic properties (<ref>)-(<ref>) have been demonstrated in Fig. <ref>.









The universality class of quantum criticality is determined by the critical exponents.

As we have seen in Fig. <ref>, the 1D attractive Hubbard model has a rich phase diagram.

At least one branch of the density of states shows sudden change when the driving parameters vary across the phase boundary in the phase diagram.

The singular behavior of thermodynamic properties is uniquely determined by the critical exponents,  which are independent of the microscopic details of the system.

Indeed, quantum criticality of quantum many-body systems  depends  solely  on the  dimensionality and the symmetry of the Hamiltonian.

Here we expand the above equations of state for the thermodynamic quantities in the limit |μ-μ_c| ≪ T.

We  derive the scaling forms of the thermodynamics at quantum criticality and thus read off the critical exponents.


 
We find that the suddenly changed density of state usually results in a quantum phase transition,  so that the thermodynamical properties can be cast into  the forms of universal quantum scaling functions in the critical region.
  
For example,  for the phase transition from the fully-paired phase V to the FFLO-like state IV,
 thermodynamic quantities of  excess fermions  display   the singular parts in the scaling functions,
 whereas the thermodynamic properties of the bound pairs present  the  regular parts.
  
In contrast to the attractive SU(2) Fermi gas, the half-filling  phase in the attractive Hubbard model
contributes a constant regular part to the thermodynamic quantities due to its unique band-filling.


Our results for the scaling functions of particle density, magnetization, compressibility and susceptibility are summarized as follows:

  
∙ phase transition  (I-V),

    n_c=   -√(2|u|/π)  T^1/2 Li_1/2( -exp( 2μ-2μ_c1/T) ),
    
    m≈     0 ,
    κ=    -2 √(2|u|/π)  T^-1/2 Li_-1/2( -exp( 2μ-2μ_c1/T) ),
    χ≈     0.

∙ phase transition (I-II),

    n_c=    -1/2√(π) T^1/2Li_1/2(-exp( μ-μ_c2/T)  ), 
    
    m=   -1/4√(π) T^1/2Li_1/2(-exp( μ-μ_c2/T)  ), 
    κ=   -1/2√(π)  T^-1/2Li_-1/2( -exp( μ-μ_c2/T) ), 
    χ=    -1/4√(π)  T^-1/2Li_-1/2( -exp( μ-μ_c2/T) ).

∙ phase transition (II-III),

    n_c=     1+1/2√(π)  T^1/2 Li_1/2( - exp( -μ-μ_c3/T) ),
    
    m=   1/2+1/4√(π)  T^1/2 Li_1/2( - exp( -μ-μ_c3/T) ),
    κ=    -1/2√(π)  T^-1/2Li_-1/2( - exp( -μ-μ_c3/T) ),
    χ=    -1/4√(π)  T^-1/2Li_-1/2( - exp( -μ-μ_c3/T) ).

∙ phase transition (II-IV),

    n_c=   n_b4+λ_1 T^1/2Li_1/2( -exp( 2(μ-μ_c4)/T) ),
    
    m=   m_b4+λ_2 T^1/2Li_1/2( -exp( 2(μ-μ_c4)/T) ),
    κ=   κ_b4+λ_3 T^-1/2Li_-1/2(-exp( 2(μ-μ_c4)/T) ),
    χ=   χ_b4+λ_4   T^-1/2Li_-1/2( -exp(2(μ-μ_c4)/T) ).

∙ phase transition (V-IV),

    n_c=   n_b5+λ_5   T^1/2Li_1/2( -exp( μ-μ_c5/T) ),
    
    m=    -1/4√(π)  T^1/2 Li_1/2( -exp( μ-μ_c5/T) ),
    κ=   κ_b5+ λ_6   T^-1/2 Li_-1/2( -exp( μ-μ_c5/T)  ),
    χ=     -1/4√(π)  T^-1/2Li_-1/2( -exp( μ-μ_c5/T) ).


In the above scaling forms some constants are given in Appendix <ref>.

These  scaling forms can be cast into the form of well known universal scaling laws.

For example, the universal scaling laws for  the density and compressibility read <cit.>

    n(μ,B,T)=     n_0(μ,B,T)+T^d/z+1-(1/ν z)𝒢( μ-μ_c/T^1/ν z), 
    κ(μ,B,T)=    κ_0(μ,B,T)+T^d/z+1-(2/ν z)ℱ( μ-μ_c/T^1/ν z),


where n_0 and κ_0 are the regular parts,  i.e., the background values before the phase transition.

Meanwhile 𝒢(x)=Li_1/2(x), ℱ(x)=Li_-1/2(x) give the scaling functions in the singular parts.

From the above scaling forms, we  read off the dynamical exponent z=2 and correlation critical exponent ν=1/2.

This scaling theory is valid for all  phase transitions across the phase boundaries in the phase diagram <ref>.

Such universal scaling laws are demonstrated in Fig. <ref> for various phase transitions.


The above scaling forms are observed to give the same critical exponents which characterize the universality class of free-fermion criticality.

An intuitive explanation for this result is that the phase transitions occurred in the 1D Hubbard model  have a common feature:  at least one branch of Fermi sea vanishing, namely ε^u,b(0) =0.

This naturally leads to a change in dispersion, i.e., a linear dispersion vanishes while a quadratic dispersion is created when the phase transition occurs.

This change in dispersion underlies a universality class of quantum criticality, see also the recent studies of the 1D interacting Bose gas <cit.>
and the 1D Heisenberg spin chain <cit.>.





Moreover, the phase V in the phase diagram Fig. <ref> shows a gapped phase (fully paired phase),
where the susceptibility reveals a particular exponential decay at low temperatures.

Using the equation of state with the pressures (<ref>), (<ref>), and (<ref>)-(<ref>),
we further  show  that the susceptibility decays exponentially with the energy gap induced by the ferromagnetic ordering, namely


    χ≈T^-1/2/4√(π)e^-Δ/T,


where the energy gap is given by Δ=ε^u(0)=-2-μ-2u-B+p^b with  p^b=4(2π-q^3/3)/3|u|π( 1+2|u|πμ/2π-q^3/3)^3/2.

This result can also be obtained by applying Sommerfeld expansion in  <ref>.

We  approximately obtain   the susceptibility
χ≈1/2γ' ≈ -T^-1/2/2√(π)Li_-1/2(-e^-Δ/T) from <ref>.

In the next section, we further demonstrate the macroscopic nature of the susceptibility in the FFLO phase.





§ FREE FLUIDS AND ADDITIVITY RULES


Fermi liquid theory is believed to break down in 1D strongly correlated systems due to the absence of well defined quasi-particles <cit.>.

Consequently the TLL theory is generally believed to describe the collective low-lying excitations in 1D many-body systems.

Despite such a big difference in the microscopic origins of the two low-energy theories,
both the Fermi liquid and the TLL share a common feature – a small  distortion of the Fermi surface or Fermi points results in
the universal low-energy physics of many-body systems.

From the results of the last section, we observed that at very low temperatures the low-energy physics of the FFLO-like  state
is governed by the universality class of a two-component TLL.

However, in view of  the macroscopic properties of the 1D attractive Hubbard model,
we argue that such a universality class of two-component TLL reveals an important nature of  free fluids.

In order to show this elegant nature, we will introduce two  effective chemical potentials for the excess fermions and bound pairs on a 1D lattice.

Then we will show  that the thermodynamic properties in the FFLO-like phase   behave like  two independent free fluids.

In particular, we find simple  additivity rules for the compressibility and susceptibility which represent a universal characteristic of quantum liquids at the renormalization fixed point.


Prior to a discussion of the free fluids, we first make an approximation for the zero temperature TBA equations <ref> in the low density regime,

    ε^1(k)=    k^2-μ_1 -a_1 ⋆ε^2(k), 
    ε^2(Λ)=    α_1 Λ^2-α_1 μ_2 - a_1 ⋆ε^1(Λ) - a_2 ⋆ε^2(Λ),


where a_m⋆ε^n(x)=∫_-y_c^y_cy   a_m(x-y)ε^n(y) with y_c being the Fermi point of ε^n(y), i.e., ε^n(y_c)=0.

In the above equations, we introduced the effective chemical potentials for excess fermions and bound pairs as

    μ_1=   μ+B+2u+2,
    μ_2=   2/α_1(μ+2√(u^2+1)-2|u| ).


The effective chemical potential of the bound pairs  reveals a deep physical insight into the crossover from Bose-Einstein condensate (BEC) to Bardeen-Cooper-Schrieffer (BCS) superconductor.

Later we shall see these effective chemical potentials reveal an important free quantum liquid nature.

We will show  that for the balanced case the effective chemical potential μ_2 varies from the kinetic energy of bound pairs to the free Fermi energy when the interaction changes from negative infinity to zero.

This reveals a 1D analogue of the BEC-BCS crossover.

This form of  the TBA equations is useful to access the ground state properties, such as sound velocities, stiffness and effective chemical potentials.

By virtue of Eqs. (<ref>) and  (<ref>) we rewrite the free energy per site (<ref>) as

    f=u+∫_-k_c^k_ck/2πε^1 (k) + ∫_-Λ_c^Λ_cΛ/2πβ_1 ε^2(Λ),

where ε^u=ε^1 and ε_1^'=ε^2.


We now proceed to  calculate  the TLL parameters of the model and compare them with those of the 1D  attractive SU(2) Fermi gas <cit.>.

The basic idea is to express the effective chemical potentials in terms of the Fermi points by employing iteration
of <ref>.

Using  the fact  that the two  dressed energies vanish at their corresponding Fermi points, we  express those Fermi points in terms of  the densities of excess fermions and bound pairs.

We shall  see that this process leads to a separation of two free fluids in the ground state energy per site.


In order to simplify the lengthy iterations, we firstly rescale the TBA equations (<ref>) and (<ref>) by defining ε̃^n=ε^n/u^2, μ̃_n=μ_n/u^2, ỹ_c=y_c/|u|, and ã_n(x)=n/π1/n^2+x^2. Then we introduce a vector presentation of the rescaled TBA equations.

In view of the properties of even functions, we utilize the base {k^2n} and {Λ^2n} (n=0,1,2,…)
to expand these scalar equations, thus we  have,  respectively

    ε⃗^ 1=    V⃗^1-𝐀^1(Λ̃_c) ε⃗^ 2,
    ε⃗^ 2=    V⃗^2-𝐀^1(k̃_c) ε⃗^ 1-𝐀^2(Λ̃_c) ε⃗^ 2.

The vectors  V⃗^1=[-μ̃_1,1,0,…]^t and V⃗^2=[-α_1μ̃_2,α_1,0,…]^t
are the driving terms and the  superscript t represents transpose operation.

The matrix 𝐀^n(ỹ_c) ε⃗ corresponds to the integral ∫_-ỹ_c^ỹ_cy  ã_n(x-y) ε̃(y).




Furthermore, as we only retain the first few leading terms, ε⃗^ n and 𝐀^n(ỹ_c)
can be expanded as sums of a few leading orders with  respect to y_c, i.e.,
ε⃗^ n=ε⃗_(0)^ n+ε⃗_(1)^ n+ε⃗_(2)^ n+… and
𝐀^n(ỹ_c)=𝐀_(1)^n(ỹ_c)+𝐀_(3)^n(ỹ_c)+𝐀_(5)^n(ỹ_c)+….

More details of the latter expansion are presented  in <ref>.

Substituting  these expansions into the TBA equations of vectorial form <ref>,
and sorting  terms order by order,  leads to   the set of  equations



    ε⃗_(0)^ 1=    V⃗^1,     ε⃗_(0)^ 2= V⃗^2, 
    ε⃗_(1)^ 1=   -𝐀^1_(1)(Λ̃_c) ε⃗_(0)^ 2, 
    ε⃗_(1)^ 2=   -𝐀^1_(1)(k̃_c) ε⃗_(0)^ 1 - 𝐀^2_(1)(Λ̃_c) ε⃗_(0)^ 2, 
    ε⃗_(2)^ 1=   -𝐀^1_(1)(Λ̃_c) ε⃗_(1)^ 2, 
    ε⃗_(2)^ 2=   -𝐀^1_(1)(k̃_c) ε⃗_(1)^ 1 - 𝐀^2_(1)(Λ̃_c) ε⃗_(1)^ 2, 
    ε⃗_(3)^ 1=   -𝐀^1_(1)(Λ̃_c) ε⃗_(2)^ 2 - 𝐀^1_(3)(Λ̃_c) ε⃗_(0)^ 2,  
    ε⃗_(3)^ 2=   -𝐀^1_(1)(k̃_c) ε⃗_(2)^ 1 - 𝐀^2_(1)(Λ̃_c) ε⃗_(2)^ 2 - 𝐀^1_(3)(k̃_c) ε⃗_(0)^ 1
       - 𝐀^2_(3)(Λ̃_c) ε⃗_(0)^ 2.



It is easy to solve the above vectorial forms ε⃗_(r)^ 1 and ε⃗_(r)^ 2 with r=1,2,3.
We then substitute these results into the scalar expression of the rescaled TBA equations.

Together with  ε̃^n(ỹ_c)=0 and the expansion μ̃_n=μ̃_n^(2)+μ̃_n^(3)+μ̃_n^(4)+⋯,
we then obtain a set of recurrence equations for μ̃_n^(2), μ̃_n^(3) and μ̃_n^(4).

Here we observe that the expansions for chemical potentials begin from n=2 due to the fact that ε̃^1(k̃_c)=-μ̃_1+k̃_c^2+o( k̃_c^3)=0, i.e.,  μ̃_1=k̃_c^2+o( k̃_c^3). Similarly for μ̃_2.

We solve these equations and then express the solution as the vectorial equation

    [
    [     μ̃_1; α_1 μ̃_2;          ]]
    =(𝐈+2/3𝐓)
    [
    [     k̃_c^2; α_1 Λ̃_c^2;            ]],

where the matrix 𝐓 is given by

    𝐓 =1/π[
    [     0 2Λ̃_c; 2k̃_c  Λ̃_c;       ]].

Finally, the free energy per site <ref>  is expressed in terms of Fermi points k_c and Λ_c, with result

    f=-2/3π(k_c^3+α_1β_1 Λ_c^3)+u.


We now proceed  to obtain the particle densities  in terms of  the Fermi points.

To this end,  we turn to the total particle density n_c and magnetization m̅ per site based on <ref>,

    n_c=   -∂ f/∂μ= -∫_-k_c^k_ck/2π∂ε^1/∂μ - β_1 ∫_-Λ_c^Λ_cΛ/2π∂ε^2/∂μ, 
    m̅=   -∂ f/∂ B= -∫_-k_c^k_ck/2π∂ε^1/∂ B - β_1 ∫_-Λ_c^Λ_cΛ/2π∂ε^2/∂ B.


In order to get  closed forms for these two properties,
we  first take partial derivatives of <ref> with respect to μ and B, respectively.
Then we   rewrite these integral equations in terms of the vectorial forms similar to  <ref>.

Finally, by lengthy  iteration and after some manipulations, we obtain

    [
      [ ñ_1; ñ_2;     ]]
    =1/π( 𝐈-𝐓+𝐓^2)^t
    [
      [     k̃_c; β_1 Λ̃_c;          ]],

where ñ_r=n_r/|u| (r=1,2) with n_1=m̅ and n_2=(n_c-n_1)/2 being respectively the densities for the excess fermions and bound pairs.

Here the redefined m̅=2m is introduced according to the original TBA equations.

 An inverse of <ref> gives the cut-off momenta in terms of the densities n_1,2

    k_c   ≈   π n_1 ∑_n=0^3 ( 2n_2/|u|)^n, 
    Λ_c    ≈   π n_2/β_1∑_n=0^3 [ 2n_1+n_2/β_1|u|]^n.



For the next step, substituting Eq. (<ref>) into (<ref>),
leads to separating the ground state energy per site into the energies of excess fermions and bound pairs, with result

    e=e_1+e_2+e_b.

Here e_b is the binding energy and the subscripts 1 and  2 denote the  excess fermions and  bound pairs, respectively.

The terms are given explicitly by

    e_1=   π^2/3  n_1^3 [1+2(2n_2/|u|)+3(2n_2/|u|)^2], 
    
    e_2=   π^2/3 α_1 n_2^3/β_1^2[1+2(2n_1+n_2/β_1|u|)+3(2n_1+n_2/β_1|u|)^2], 
    
    e_b=   -(2u+2)n_1-4(u+√(u^2+1))n_2.

As usual,  we define a dimensionless interaction strength  γ_s=2|u|/n_s (s=1,2) <cit.>. Using the relation

    K_s=π/√(3e(γ_s)-2γ_s e(γ_s)/γ_s+1/2γ_s^2^2 e(γ_s)/γ_s^2),


the Luttinger parameters for the excess fermions and bound pairs can be directly worked out to be

    K_1=1,    K_2=2√(2) β_1/√(α_1)[1-2/β_1 γ_2+1/(β_1 γ_2)^2].



We note that the Luttinger parameter K_2 in the fully paired phase depends explicitly on the lattice parameters α_1 and β_1.

This behavior is different from the constant value K_2=4 for the bound pairs phase of the strongly attractive SU(2) Fermi gas <cit.>.

In the limits  u→ 0 and n_s/|u|  small,  the lattice parameters α_1→ 2, β_1→ 2.

Thus  we have  K_2=4 which is the same as for the SU(2) Fermi gas.

The two limits u → 0 and n_s/|u|≪ 1  represent the lattice-gas mapping  between 1D attractive Hubbard model and SU(2) Fermi gas <cit.>.


Beside this framework of the TLL theory, we also find that for  the  low density case,
the chemical potentials for the unpaired fermions and pairs are given explicitly by

    μ_1   =   π n_1^2A_1^2+4π^2α_1/3β_1^3|u|n_2^3A_2^3, 
    μ_2   =   π^2n_2^2/β_1^2A_2^2+4π^2/3α_1|u|n_1^3A_1^3+2π^2/3β_1^3|u|n_2^3A_2^3 ,

where    A_1=1+2n_2/|u|+(2n_2/|u|)^2 and A_2=1+2n_1+n_2/β_1 |u|+(2n_1+n_2/β_1 |u|)^2,
which  indicate  interacting effects among pairs and unpaired fermions.

We observe that the chemical potential μ_2 tends to  the kinetic energy of bound pairs  in the BEC limit |u| →∞.

Whereas in the weak coupling limit, |u| → 0, μ_2 tends to  the Fermi energy of the free fermions on a 1D lattice.



The effective chemical potentials (<ref>) and (<ref>) reveal  that the thermodynamic quantities could be separable,
i.e., the total is equal to a sum of the effective thermodynamic quantities of  two individual constituents.


Here we further derive the additivity rules for the compressibility and susceptibility.

For the compressibility, using the standard thermodynamic relation κ=(∂ n_c/∂μ)_B,
the derivatives of the density and effective chemical potentials for fixed magnetic field could be further expressed as
n_c=n_1+2 n_2 and μ_1=α_1/2μ_2=μ, respectively.

Inserting these relations into the definition of compressibility, κ=∂ n_c/∂μ|_B=n_1+2n_2/μ,
we thus obtain

    κ=κ_1+2/α_1κ_2.

Here the effective compressibilities of excess fermions and bound pairs are defined as
κ_1=(∂ n_1/∂μ_1)_B and κ_2=2(∂ n_2/∂μ_2)_B.

Details are given in see in Appendix  <ref>.

The additivity rule (<ref>) for the compressibility can be confirmed numerically, as shown  in Fig. <ref>(a).





For the susceptibility in the canonical ensemble, defined as χ̅=(∂m̅/∂ B)_n_c,
it is straightforward to see n_c=n_1+2n_2=0 and B=μ_1-α_1/2μ_2,
and thus the additivity rule

    1/χ̅=1/χ̅_1+α_1/21/χ̅_2.


Here χ̅_1=(∂ n_1/∂μ_1)_n_c and χ̅_2=2(∂ n_2/∂μ_2)_n_c
are the effective susceptibilities for excess fermions and bound pairs, respectively.

These explicit expressions for the effective thermodynamic quantities can be found in Appendix <ref>.

The additivity rule (<ref>) for the susceptibility can  also be confirmed  numerically, as shown in Fig. <ref>(b).

Similar to the observation concerning TLL parameters,
the additivity rules for the 1D attractive  Hubbard model also reduce to those for the SU(2) Fermi gas through the lattice-gas mapping.

In Appendix  <ref> we calculate the individual compressibility and susceptibility explicitly.



The simple additivity nature of the thermodynamics at low temperatures
characterizes the universal low energy physics of the FFLO-like state of the 1D attractive Hubbard model.

In this sense, the additivity rules reflect a universal nature of the multicomponent  TLL in 1D.

The simple additivity rule thus reveals the significant two free fluid nature of the FFLO phase,
as predicted in expansion dynamics of the FFLO state in 1D <cit.>.

The  macroscopic magnetic properties in the FFLO-like phase show the properties of the ordinary higher-dimensional Fermi liquid,
see Fig. <ref>.
 
This figure shows that in  the free fluids  region the magnetization is nearly temperature independent.
 
In the non-Fermi liquid region thermal fluctuations gradually overwhelm quantum fluctuations.
 
 Thus the magnetization has a uniform temperature dependence for different magnetic fields, indicating paramagnetism.

The non-Fermi liquid crossover region reveals a scaling invariance, which was studied in Section IV.

Such Fermi liquid-like features have been found in the spin compound Cu(C_4H_4N_2)(NO_3)_2 <cit.> and the heavy fermion material YbNi_4P_2 <cit.>.

The study of Fermi and non-Fermi liquids in 1D has received significant recent interest <cit.>.





Using the explicit expressions for the compressibility (<ref>) and susceptibility (<ref>),
we may calculate the Wilson ratio, which is a dimensionless ratio defined as the susceptibility or compressibility over the specific heat divided by the temperature.

The Wilson ratio is the ratio describing quantum fluctuations and energy thermal fluctuations.

Both the Fermi liquid and TLL give a constant Wilson ratio <cit.>, i.e., two types of fluctuations are on equal-footing in temperature scaling.

However, near a critical point, the dimensionless Wilson ratios exhibits a sudden enhancement indicating a sudden change in the density of state.
  
Therefore the Wilson ratios serves as a powerful tool for distinguishing the phases of a quantum liquid  and for determining the finite temperature phase diagram as well.


The compressibility Wilson ratio R_W is determined by

    R_W^κ   =   π^2 k_B^2/3 κ/C_v/T
       =   π(κ_1+2/α_1κ_2)/( 1/v_1+1/v_2),

where we have used <ref> to calculate the specific heat and set the Boltzmann constant to k_B=1.

This Wilson ratio vanishes in both phases I (vacuum) and III (half-filling phase).

In the limit n_c/|u| → 0 the compressibility Wilson ratio for phases II and IV are respectively,
R_W^κ=1 and R_W^κ=2√(2)β_1/√(α_1).

These results turn out to be the same as for the  strongly attractive SU(2) Fermi gas <cit.>
when  the limit u → 0 is applied.
 
On the other hand, the susceptibility Wilson ratio is defined by
R_W^χ=4/3(π k_B/μ_B  g_L)^2 χ/C_v/T
with Bohr magneton μ_B and Lande factor g_L.

Fig. <ref> shows a contour plot of each type of Wilson ratio which demonstrates the macroscopic feature of the Fermi liquid nature.

This figure also presents the low-temperature phase diagram in the B-μ plane.






§ CONCLUSION

In summary we have presented a framework to determine the nature of quantum criticality and quantum liquids in the 1D attractive Hubbard model.

We have obtained the universal thermodynamics of the  model by solving  the TBA equations.

In particular, we have analytically derived the equation of state at low temperatures,
from which we have obtained effective chemical potentials of excess fermions and bound pairs, along with
the density, compressibility, susceptibility and specific heat in terms of the chemical potential μ, magnetic field B,  temperature T and interaction strength constant.

At quantum criticality  the  scaling forms of these thermal and magnetic properties have been obtained.

The dynamical exponent z=2 and correlation critical exponent ν=1/2, indicating the universality class of criticality of free fermion theory.


Our results provide strong evidence for the existence of two free fluids of bound pairs and of unpaired fermions,
which were noticed in the expansion dynamics of the FFLO state in 1D <cit.>.

Regarding the nature of the two fluids in the attractive Hubbard model,  we have  shown  that   in the low-density regime the interaction
effect resulting from the paired and unpaired fermions   can be absorbed into  effective chemical potentials
of two non-interacting ideal gases.

Consequently, the  additivity rules in the compressibility and susceptibility of the 1D attractive Hubbard model
hold as long as  the dimensionless Wilson ratio remains a constant.

This behavior significantly reflects the free fluids  nature in  thermodynamic properties of  the model.

In this phase, the FFLO pair correlation function

    G_p(x,t)   =   ⟨Ψ_↑^†(x,t) Ψ_↓^†(x,t) Ψ_↑(0,0) Ψ_↓(0,0) ⟩
       ≈    A_p,1cos( π(n_↑-n_↓)x )/|x+i  v_1   t|^2θ_1 |x+i  v_2   t|^2θ_2
          + A_p,2cos( π(n_↑-3n_↓)x )/|x+i  v_1   t|^2θ_3 |x+i  v_2   t|^2θ_4,


shows a typical spatial oscillation which  is a characteristic  of the FFLO state.

In the above equation, the exponents θ_1 ≈ 1/2, θ_2 ≈ 1/2+n_2/|u|β_1,  θ_3 ≈1/2-4  n_2/|u|β_1 and θ_4 ≈5/2-4 n_1/|u|-3  n_2/|u|β_1  depend essentially on the lattice parameter β_1.

Here n_2,1=N_2,1/L are  the dimensionless  densities  of pairs and  unpaired fermions, with the sound velocities v_1,2 given in (<ref>).

The study of the FFLO pair correlation is presented elsewhere <cit.>.

To conclude, we note that our work provides benchmark physics of the 1D attractive Hubbard model of relevance to experiments with ultracold  fermionic atoms on  lattices.



Acknowledgments.  The authors SC and YCY contributed equally to the calculations in this paper.
The authors  thank R. Hulet  for helpful discussion. This work is supported by
Key NNSFC grant number 11534014, MOST grant number 2017YFA0304500,
NNSFC grant numbers 11374331, 11174375 and ARC Discovery Projects DP130102839, DP170104934.






Afigure
Atable
Aequation
Alabel












§ APPENDICES




 §.§ Wiener-Hopf method



The phase boundary between phases IV and V is determined by the conditions ε^u(0)=0 and ε^b(0)< 0,
which imply that Q=0 and A is finite.

Thus at zero temperature the TBA equations are simplified to

    ε^u(k)=   -2cos k-μ-2u-B 
       -∫_-A^A Λ  a_1(sin k-Λ)ε_1^'(Λ), 
    ε_1^'(Λ)=   -2μ-2∫_-π^πk  cos^2 k    a_1(sin k -Λ) 
       -∫_-A^A Λ^'  a_2(Λ-Λ^') ε_1^'(Λ^').

Particularly, if A=∞, it follows that the chemical potential  μ=0.
The intersection of the phase boundary with  the B-axis  could be calculated exactly by the  Fourier transformation

    B_c1=2|u|-2+2∫_0^∞ω J_1(ω)exp(-|u|ω)/wcosh(u ω).


Now we  consider the more general case A ≫ 1, for  which  the phase boundary can be resolved using the Wiener-Hopf method.

By applying Fourier transformation on <ref> and after some algebraic manipulations, we have


    ε_1^'(Λ) = -μ-∫_-∞^∞ω J_1(ω)/ωcosh(u ω)exp(i ω Λ)
    +∫_0^∞Λ^' ε_1^'(Λ^'+A)  [R(Λ-Λ^'-A)+R(Λ+Λ^'+A)],


where we have introduced the function

    R(x)=∫_-∞^∞w/2π exp(i w x)/1+exp(2|u w|).


Substituting  y(Λ)=ε_1^'(Λ+A) and expanding  y(Λ)=∑_n=0^∞ y_n(Λ)
in terms of powers of Λ in <ref>, the result can be  separated
 into a series of Wiener-Hopf integral equations in terms of the functions  y_n(Λ), namely

    y_n(Λ)=g_n(Λ)+∫_0^∞Λ^'  R(Λ-Λ^')   y_n(Λ^').

Here we denote the driving terms

    g_0(Λ)   =   -μ-∫_-∞^∞ω J_1(ω) e^iω(Λ+A)/ωcosh(u ω), 
    
    g_n(Λ)   =   ∫_0^∞Λ^'  R(Λ+Λ^'+2A)   y_n-1(Λ^').


To solve these integral equations for y_n(Λ), we begin by defining

    ỹ_n^±(ω)=   ∫_-∞^∞Λ θ_H(±Λ)  y_n(Λ) e^iωΛ,

where ỹ_n^+(ω) (ỹ_n(ω)) is an analytic function in  the upper (lower) half-plane.

It is  obvious that  the Fourier transformation of y_n(x) satisfies  the relation ỹ_n(ω)=ỹ_n^+(ω)+ỹ_n^-(ω).

From <ref> it follows that

    ỹ_n^+(ω)  1/1+exp(-2|u| |ω|) + ỹ_n^-(ω)=g̃_n(ω),

by applying Fourier transformation. We further decompose the denominator  1+exp(-2|u| |ω|) into a product of two pieces,

    1+exp(-2|u|  |ω|)=G^+(ω)G^-(ω),

where G^+(ω) (G^-(ω)) is an analytic function in the upper (lower) half-plane.

Then substituting this last equation into <ref> results in the  form

    ỹ_n^+(ω)/G^+(ω)+G^-(ω)ỹ_n^-(ω)=G^-(ω)g̃_n(ω).


Furthermore, we decompose  G^-(ω)g̃_n(ω) into a sum of two pieces,

    G^-(ω)g̃_n(ω)=Q_n^+(ω)+Q_n^-(ω),

where similarly Q_n^+(ω) (Q_n^-(ω)) is an analytic function in the upper (lower) half-plane.

Then substitution of this last equation into <ref> gives

    ỹ_n^+(ω)=    G^+(ω)Q_n^+(ω),
    ỹ_n^-(ω)=    Q_n^-(ω)/G^-(ω).

In this way we can  work out the Fourier transformation of y_0(Λ) and  y_n(Λ) itself.

To this end, recalling (<ref>), we firstly decompose 1+exp(-2|u|  |ω|) as

    G^+(ω)=G^-(-ω) 
    
    =   √(2π)/Γ(1/2-  i  |u|  ω/π)( -i  |u|ω/π)^-i |u|ω/πexp(i |u|ω/π),

where we should note that lim_ω→∞ G^±(ω)=1, along with the special values
G^±(0)=√(2) and G^±(±iπ/2|u|)=√(π/e) of these functions.

The decomposition for G^-(ω)g̃_n(ω) in general is subtle,
however the leading case G^-(ω)g̃_0(ω) is accessible.
We start analysis from the Fourier transformation of g_0(Λ),

    g̃_0(ω)=-μ  2πδ_D(ω)-2π J_1(ω)exp(-i ω A)/ωcosh(uω),

where on the rhs the δ_D function could be decomposed as

    2πδ_D(ω)=i( 1/ω+i ϵ-1/ω-i ϵ)     (ϵ→ +0 ).

The second term on the rhs is a meromorphic function of ω with poles located at

    ω_n=iπ/2|u|(2n+1)     (n ∈𝐙)

originating from the term 1/cosh(uω), implying the decomposition

    1/cosh(uω)=    χ^+(ω)+χ^-(ω), 
    
    χ^+(ω)=   i/|u|∑_n=0^∞ (-1)^n 1/ω+ω_n, 
    
    χ^-(ω)=   1/cosh(uω)-i/|u|∑_n=0^∞ (-1)^n 1/ω+ω_n,

where χ^+(ω) and χ^-(ω) are analytic functions in the upper and lower half-planes, respectively.

With the help of <ref>, as for any analytic and bounded function f^-(ω) in the lower half-plane, the decomposition of f^-(ω)/cosh(u ω) is

    f^-(ω)/cosh(u ω)=    F^+(ω)+F^-(ω),
    
    
    F^+(ω)=    i/|u|∑_n=0^∞ (-1)^n f^-(-ω_n)/ω+ω_n,
    
    
    F^-(ω)=    f^-(ω)/cosh(u ω)-F^+(ω).


By virtue of <ref>, we make the following decomposition for G^-(ω)g̃_0(ω),

    Q_0^+(ω)=   -i μ  G^-(0)/ω+i ϵ-q(ω)
    
    Q_0^-(ω)=   i μ  G^-(0)/ω+i ϵ
    - 2π J_1(ω)exp(-i ω A)G^-(ω)/ωcosh(u ω) + q(ω)

where q(ω)=4i∑_n=1^∞ (-1)^n G^-(-i h_n)I_1(h_n)exp(-h_n A)/(2n+1)(ω+i h_n), I_1(z)
is the first order modified Bessel function, h_n=π/2|u|(2n+1), with the series converging only if A>1.

If A ≫ 1, using <ref>, we have

    y_0^+(ω)=G^+(ω)[ -i μ G^-(0)/ω+i ϵ- q(ω) ].

Obviously, we know y(0)=ε_1^'(A)=0, which implies

    0=y(0)=lim_ω→∞ -i ω ỹ^+(ω).

Hereafter we replace y(ω) with y_0(ω), which is a reasonable approximation if A ≫ 1. Therefore <ref> give rise to

    μ=-4∑_n=0^∞G^-(-i h_n)I_1(h_n)exp(-h_n A)/(2n+1)G^-(0).


Since  we have obtained a parametric expression for the critical chemical potential, we turn  to the expression for the magnetic field.
Due to the fact that  the phase boundary is determined by ε^u (0)=0, we thus use  <ref> to determine the magnetic field.

For simplicity, we rewrite <ref> as

    ε^u(k)=   -2cos k-μ-2u-B 
       + ∫_A^∞Λ [ a_1(sin k - Λ) + a_1(sin k + Λ) ] ε^'_1(Λ) 
       -∫_-∞^∞Λ  a_1(sin k - Λ) ε_1^'(Λ), 
    ε_1^'(Λ)=    ε_1^' (0)(Λ) - ∫_-A^A Λ^'  a_2(Λ-Λ^') ε_1^'(Λ^'),

where we have denoted

    ε_1^'(0)(Λ)=-2μ-2∫_-π^πk  cos^2 k    a_1(sin k - Λ).


Substituting  <ref> into the last term on the rhs of <ref> gives

    ε^u(k) 
    
    =   -2cos k-μ-2u-B 
       +∫_0^∞Λ [s(Λ+A-sin k)+ s(Λ+A+sin k) ]   y(Λ) 
       -∫_-∞^∞Λ  s(Λ-sin k) ε_1^' (0)(Λ),

where we have introduced the function s(x)=1/4|u|cosh(π x/2|u|) and made use of the two identities

    1/4|u|cosh(π x/2|u|)   =   ∑_n=0^∞ (-1)^n a_2n+1(x), 
    ∫_-∞^∞y   a_n(x-y) a_m(y-z)   =    a_m+n(x-z).


Substituting the expansion s(x)=1/2|u|∑_n=0^∞ (-1)^n exp(-h_n x),
where | π x/u| < 1 and <ref> into <ref>,
and after some algebraic manipulations, we arrive at the result

    ε^u(k)=   -2cos k-2u-B 
       +∑_n=0^∞(-1)^n/|u| ỹ^+(i h_n) cosh(h_n sin k) exp(-h_n A) 
       + 2∫_0^∞ω J_1(ω) cos(ωsin k) exp(-|u|ω) /ωcosh(uω),


Using <ref> and ε^u (0)=0 we derive the  expression

    B=   -2+2|u|+∑_n=0^∞(-1)^n/|u| ỹ^+(i h_n) exp(-h_n A) 
       +2∫_0^∞ω J_1(ω) exp(-|u|ω) /ωcosh(uω)

for determining the critical magnetic field.

Here we denoted  h_n= π/2|u|(2n+1).

The equation (<ref>) sets up a relation between the magnetic field and the chemical potential.

In summary, the phase boundary between phase IV and V is determined  by <ref>  for  A ≫ 1.



Bequation



 §.§  Derivation of the Equation of State 



The derivation of the equation of state is rather involved.

Here we sketch the calculations for the terms p^u and p^b_n.

Prior to substituting the dressed energies into the definitions of p^u and p^b_n integrating by parts,
we first need to find a suitable form of the TBA equations for this procedure, i.e., (<ref>) and (<ref>).

For simplicity in later discussion, we approximate the definition of p^b_n as

    p^b_n=    ∫_-∞^∞Λ/π1/√(Λ+i  n   |u|)ε_n^' -(Λ)  
    
    =    T∫_-∞^∞Λ/2π∫_-π^πk   a_n(Λ-sin k) 
    
    =    ∫_-∞^∞Λ ε_n^' -(Λ)  Δ_n(Λ) + o( 1/|u|^4),


where Δ_n(Λ)=a_n(Λ)-1/2 b_n(Λ)+ 2Λ^2 b_n(Λ).

In the above equations, we used the abbreviations

    ε_n^' -(x)    =   Tln(1+e^-ε_n^'(x)/T),
    ε_n^ -(x)    =   Tln(1+e^-ε_n(x)/T).



To obtain the result (<ref>), regarding the first series of integral terms on the rhs of (<ref>),
we expand them in the strong coupling regime as

    ∑_n=1^∞∫_-∞^∞Λ  a_n(sin k-Λ)ε_n^' -(Λ) 
    
    =   ∑_n=1^∞∫_-∞^∞Λ Δ_n(Λ) ε_n^' -(Λ)  ×
       a_n(Λ){ 1+ 2Λsin k - sin^2 k/(nu)^2+Λ^2 + [2Λsin k - sin^2 k/(nu)^2+Λ^2]^2 }/Δ_n(Λ) + o( 1/|u|^4) 
    
    =   ∑_n=1^∞ p_n^b + a̅ + 2a̅cos^2 k + o( 1/|u|^4),

where we have inserted (<ref>) and a̅ has been defined in <ref>.

While for the second series of integral terms, it is easy to see that under the assumption B/T ≫ 1,
these spin-wave contributions are no more than -T  e^-2B/Te^-K̅ I_0(K̅),
which is accessible through simple iteration of (<ref>).
In fact, the spin degree of freedom is frozen here, and thus this term could be neglected in later discussion.

We can rewrite (<ref>) as

    ε^u(k)=ε^u_0(k)-A^u,

where ε^u_0(k)=-2cos k + 2 a̅cos^2 k and A^u=μ+2u+B-∑_n=1^∞ p^b_n + a̅.

Integrating by parts in p^u, we obtain

    p^u=   T ln( 1+e^(μ+2u+B-∑_n=1^∞ p_n^b + a̅ -2)/T)
       +1/π∫_2a̅-2^2a̅+2ε^u_0 k(ε_0^u)/1+e^ε_0^u/T/z,

where z=e^A^u/T and k(ε_0^u)=arccos( 1-√(1+2a̅ε_0^u)/2a̅) represents the inverse function of ε_0^u (k).

By taking account of a̅∼∑_n=1^∞p^b_n/(nu)^2 for the  strong coupling regime, the integral in the above equation can be further simplified,

    1/π∫_2a̅-2^2a̅+2ε^u_0 k(ε_0^u)/1+e^ε_0^u/T/z
    =2/π∫_a̅-1^a̅+1x k(2x)/1+e^2x/T/z
          =   2/π∫_-1^1x arccos(-x)/1+e^2x/T/z-2a̅/π∫_-1^1x  x^2/√(1-x^2)/1+e^2x/T/z
             +2a̅/1+e^2ε^u(π)/T+o(1/|u|^4),

where we have changed the integration variable ε^u_0=2x
and then applied Taylor expansion with respect to a̅.
See ε^u(π) in <ref>. The result (<ref>) is therefor achieved.

We then turn to the transformation of ε_n^'(Λ). Similar to the treatment for ε^u(k),
we employ Taylor expansion to expand (<ref>) in  the strong coupling region, with result

    ε^'_n(Λ)   =   -2nμ- a_n(Λ) [ 2π-∫_-π^πk cos k ε^u -(k)   ] 
           - b_n(Λ) [ -π/2-∫_-π^πk cos k  sin^2 k  ε^u -(k) ] 
          + ∑_m=1^∞ T_nm∗ε_m^' -(Λ)   + o( 1/|u|^4),


where the integral terms in the brackets are denoted as d_1 and d_2, respectively.

Here d_1 and d_2 can be calculated via integration by parts, similar to that done for p^u above,
see the explicit expressions in <ref>.

With respect to the convolution term, due to the condition of low density,
the cut-off of the dressed energy ε^'_n(Λ) is small, thus in general we can make the approximations

    ∫_-∞^∞Λ^'  a_p(Λ-Λ^')ε_q^' -(Λ^')  
       =   ∫_-∞^∞Λ^'  a_p(Λ^') ε_q^' -(Λ^') 
          -Λ^2 ∫_-∞^∞Λ^'  b_p(Λ^') ε_q^' -(Λ^') +o ( 1/|u|^4),

which results in <ref> in the main text.


We next rewrite (<ref>) as

    ε^'_n(Λ)=D_n (Λ/n|u|)^2 - A^b_n,

where A^b_n=2nμ-η_n+d_1/π n |u|+d_2/2π(n|u|)^3 and D_n was defined in <ref>.

<ref> is arrived at by substituting the above equation into the definition of p^b_n and integrating by parts.

Lastly, a̅,
ξ_p^m=T∫_-∞^∞Λ^'  a_p(Λ^') ε_m^' -(Λ^')
and ϕ_p^m=T∫_-∞^∞Λ^'  b_p(Λ^') ε_m^' -(Λ^') are calculated in a similar way.



Cequation



 §.§ Some constants in the scaling functions



The constants used  in the  scaling forms for the phase transition (II-IV)  are given explicitly by

    λ_1   =   -2√(|u|)(1-q_c4/π)/√(2π-q_c4^3/3), 
    λ_2   =   (1-q_c4/π)q_c4/π/√(|u|)√(2π-q_c4^3/3),
    
    n_b4   =   γ,    m_b4=1/2γ,
    κ_b4   =   γ' (1+ 4/√(π)τ^1/2f̃_1/2 + 6/πτf̃_1/2^2 + 5/π^3/2τ^3/2f̃_1/2^3 - 2/√(π)τ^3/2f̃_3/2), 
    λ_3   =   -4√(|u|)(1-q_c4/π)/√(2π-q_c4^3/3), 
    χ_b4   =   1/2γ^ ' + 2δ γ^ '-(1-γ)δ'/|u|√(π)( τ^1/2f̃_1/2 + 1/2√(π)τf̃_1/2^2 + 1/4πτ^3/2f̃_1/2^3 . 
          .- τ^3/2f̃_3/2), 
    λ_4   =   -2(q_c4/π)^2(1-q_c4/π)/√(|u|^3)√(2π-q_c4^3/3).

Here the parameter q_c4=√(B+2-2√(1+u^2)) +1/3π|u|( B+2-2√(1+u^2)).

The constants used  in the  scaling forms for the phase transition (V-IV)  are given explicitly by

    n_b5   =   -2/√(π)τ^1/2f̃_1/2 -1/πτf̃_1/2^2 - 1/2π^3/2τ^3/2f̃_1/2^3 + 1/√(π)τ^3/2f̃_3/2, 
    λ_5   =    -1/2√(π)( 1- 4/π√(1+2π|u|μ̃_c5/2π-q_c5^3/3)),
    λ_6   =   -1/2√(π)( 1 - 8/π√(1+2π|u|μ̃_c5/2π-q_c5^3/3))
    κ_b5   =    -(1-γ) f̃_-1/2/D_0 τ^1/2( 4/√(π) + 6/πτ^1/2f̃_1/2 + 6/π^3/2τf̃_1/2^2  .
          .- 7/4πτ^3/2f̃_3/2+ 5/π^2τ^3/2f̃_1/2^3 ),

where  μ̃_c5≈ 2|u|-B-2+8√(2)/3π|u|α_1(2√(1+u^2)-B-2)^3/2
and q_c5=√(μ̃_c5+2u+B+2).


Dequation


 §.§  Explicit forms of the additivity rules 



The vectorial forms of <ref> are accessible by expanding the rescaled TBA equations in terms of
{k^2n} and {Λ^2n} (n=0,1,2…).

We here give the explicit expression for the matrix 𝐀^n(ỹ_c) (n=1,2) with the elements

    {𝐀^n (ỹ_c) }_jl=2/π∑_0≤ j ≤ i < ∞(-1)^i   C_2i^2j ỹ_c^2i-2j+2l+1/n^2i+1(2i-2j+2l+1),

with j,l=0,1,2,….
Thus the first two orders of 𝐀_(q)^n (y_c) are written as


    𝐀^1_(1)(ỹ_c)=   2/π[
    [  1  0  ⋯; -1  0  ⋯;  ⋮         ;    ]]
     ỹ_c,
       𝐀^2_(1)(ỹ_c)=    1/π[
    [    1    0    ⋯; -1/4    0    ⋯;    ⋮               ;      ]]
     ỹ_c,
    
    𝐀^1_(3)(ỹ_c)=   2/π[
    [ -1/3  1/3    ⋯;    2  1/3    ⋯;    ⋮               ;      ]]
     ỹ_c^3,
       𝐀^2_(3)(ỹ_c)=   1/π[
    [ -1/12   1/3     ⋯;   1/8 -1/12     ⋯;     ⋮                  ;       ]]
     ỹ_c^3.



Next, the partial derivatives of <ref> read


    ∂ε⃗^ 1/∂μ̃
    =   ∂V⃗^1/∂μ̃
    -𝐀^1(Λ̃_c) ∂ε⃗^ 2/∂μ̃,
    ∂ε⃗^ 2/∂μ̃
    =   V⃗^2
    -𝐀^1(k̃_c) ∂ε⃗^ 1/∂μ̃
    -𝐀^2(Λ̃_c) ∂ε⃗^ 2/∂μ̃.


With the help of the explicit forms of  𝐀^n (ỹ_c), i.e., (<ref>) and (<ref>), we can
obtain <ref>, which relates the densities and the cutoffs.

Together with <ref>, we then obtain the relation between the densities and the effective chemical potentials (<ref>) and (<ref>).


On this basis, we now proceed to  derive the explicit expressions for the effective compressibility and susceptibility
in terms of densities of bound pairs and excess fermions.

Apparently, the densities of bound pairs and excess fermions rely on the chemical potential and the magnetic field, and vice versa,
which in fact indicates under fixed magnetic field one could obtain the following results through the total derivatives,

    κ_1=( ∂ n_1/∂μ_1)_B=n_1/μ,   κ_2=2( ∂ n_2/∂μ_2)_B=α_1n_2/μ,

where we keep d B=∂ B/∂ n_1d n_1+∂ B/∂ n_2d n_2=0. Thus we have

    n_1/μ=1/J( ∂ B/∂ n_2)_n_1,   n_2/μ=-1/J( ∂ B/∂ n_1)_n_2.

Here the Jacobian determinant

    J=   ( ∂μ/∂ n_1)_n_2( ∂ B/∂ n_2)_n_1 - ( ∂ B/∂ n_1)_n_2( ∂μ/∂ n_2)_n_1
    
    =   -α_1/2[ ( ∂μ_1/∂ n_1)_n_2( ∂μ_2/∂ n_2)_n_1 - ( ∂μ_2/∂ n_1)_n_2( ∂μ_1/∂ n_2)_n_1],

where we have used <ref>.

Similarly, the magnetic field is dependent on the effective chemical potentials while the latter is dependent on densities of bound pairs and excess fermions.

Therefore by application of chain rule we have

    (∂ B/∂ n_1)_n_2=   (∂μ_1/∂ n_1)_n_2
    -α_1/2(∂μ_2/∂ n_1)_n_2, 
    (∂ B/∂ n_2)_n_1=   (∂μ_1/∂ n_2)_n_1
    -α_1/2(∂μ_2/∂ n_2)_n_1.


It is obvious that once the explicit expression of μ_r in terms of n_s (r,s=1,2) is known,
our goal of the effective compressibilities is easy to achieve.
We use <ref> to derive

    μ_1=    π^2 n_1^2[1+2(2n_2/|u|)+3(2n_2/|u|)^2] 
       +4π^2α_1/3β_1^3|u|n_2^3[1+32n_1+n_2/β_1|u|],
    μ_2=    π^2 n_2^2/β_1^2[1+22n_1+n_2/β_1|u|+3(2n_1+n_2/β_1|u|)^2] 
       +4π^2/3α_1 |u|n_1^3[1+3(2n_2/|u|)] 
       +2π^2/3β_1^3 |u|n_2^3[1+32n_1+n_2/β_1|u|],

and thus

    κ_1=   π^2/J[ -α_1 n_2/β_1^2-4 α_1 n_1 n_2/|u|β_1^3+4 n_1^2/|u|-4 n_1^3/u^2+24 n_1^2 n_2/u^2. 
       .-12 α_1 n_1^2 n_2/u^2 β_1^4+6 α_1 n_2^3/u^2 β_1^4], 
    κ_2=   -2α_1π^2/J[ n_1-n_1^2/|u|+4 n_1 n_2/|u|-6 n_1^2 n_2/u^2-α_1 n_2^2/u β_1^3. 
       .+12 n_1 n_2^2/u^2-6 α_1 n_1 n_2^2/u^2 β_1^4],
    
    J=   -2π^4α_1/β_1^2n_1 n_2
    [
    1+4 n_1/|u|β_1+12 n_1^2/u^2 β_1^2+4 n_2/|u|+4 n_2/|u| β_1. 
       .+24 n_1 n_2/u^2 β_1^2+8 n_1 n_2/u^2 β_1+12 n_2^2/u^2+10 n_2^2/u^2 β_1^2+16 n_2^2/u^2 β_1].


The situation for effective susceptibilities is rather simple. With fixed total particle density, one confirms that n_1+2n_2=0,
and thus the total derivative of the effective chemical potentials with respect to n_r (r=1,2) is

    μ_1
    =   [ ( ∂μ_1/∂ n_1)_n_2 -1/2( ∂μ_1/∂ n_2)_n_1]  n_1, 
    μ_2
    =   -2[(∂μ_2/∂ n_1)_n_2
    -1/2(∂μ_2/∂ n_2)_n_1] n_2.


After some algebraic manipulations we then obtain

    χ̅_1=   1 / (∂μ_1/∂ n_1
    -1/2∂μ_1/∂ n_2), 
    χ̅_2=   -1/ (∂μ_2/∂ n_1
    -1/2∂μ_2/∂ n_2),

which together with <ref> result in


    χ̅_1=   1/(2π^2)/n_1-n_1^2/|u|+4 n_1 n_2/|u|-6 n_1^2 n_2/u^2-α_1 n_2^2/|u| β_1^3+12 n_1 n_2^2/u^2-6 α_1 n_1 n_2^2/u^2 β_1^4,
    χ̅_2=   1/π^2/n_2/β_1^2-4 n_1^2/|u| α_1+4 n_1^3/u^2 α_1 +4 n_1 n_2/|u| β_1^3-24 n_1^2 n_2/u^2 α_1+12 n_1^2 n_2/u^2 β_1^4-6 n_2^3/u^2 β_1^4.




99

Hubbard1963  J. Hubbard, J. Proc. R. Soc. A 276, 237 (1963); ibid
281 401 (1964).

Yang:1967C. N. Yang,  Phys. Rev. Lett. 19, 1312 (1967).

Baxter:1972R. J. Baxter,  Ann. Phys. (N.Y.) 70, 193 (1972).


Lieb1968 E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).


Ess05F. H. L. Essler, H. Frahm, F. Göhmann, A. Klümper  and V. E. Korepin,
The One-Dimensional Hubbard Model (Cambridge University Press, Cambridge, 2005).

MT1969s M. Takahashi, Prog. Theo. Phys. 42, 1098 (1969); ibid 43, 860 (1970);
ibid 43, 1619 (1970).

MT1971M. Takahashi, Prog. Theo. Phys. 45, 756 (1971).

Shiba1972 H. Shiba, Phys. Rev. B 6, 930 (1972).


Krivnov1975 V. Ya. Krivnov and A. A. Ovchinnikov, Zh. Eksp. Teor. Fiz. 67, 1568 (1974).

Usuki1989 T. Usuki, N. Kawakami and A. Okiji, Phys. Lett. A 135, 476 (1989).

Woyna1991 F. Woynarovich and K. Penc, Z. Phys. B 85, 269 (1991). 




Coll1974 C. F. Coll, Phys. Rev. B 9, 2150 (1974).

Ovchi1970 A. A. Ovchinnikov, Sov. Phys. JETP 30, 1160 (1970).

Woyna1982s F. Woynarovich, J. Phys. C 15, 85 (1982); ibid 15, 97 (1982); ibid 16, 5293 (1983).

Klumper1990 A. Klümper, A. Schadschneider and J. Zittartz, Z. Phys. B 78, 99 (1990).

Degu2000 T.Deguchi, F. H. L. Essler, F.Göhmann, A. Klümper, V. E. Korepin and K. Kusakabe, Phys. Rep. 331, 197 (2000). 


Woyna1983 F. Woynarovich, J. Phys. C 16, 6593 (1983).


Ess1994a F. H. L. Essler and V. E. Korepin, Phys. Rev. Lett. 72, 908 (1994).

Ess1994b F. H. L. Essler and V. E. Korepin, Nucl. Phys. B 426, 505 (1994). 


MT1972 M. Takahashi, Prog. Theo. Phys. 47, 69 (1972). 

MT1974 M. Takahashi, Prog. Theo. Phys. 52, 103 (1974).


Lee1988 K.-J.-B. Lee and P. Schlottmann, Phys. Rev. B 38, 11566 (1988).

Okiji1989 N. Kawakami, T. Usuki and A. Okiji, Phys. Lett. A 137, 287 (1989);
T. Usuki, N. Kawakami and A. Okiji, J. Phys. Soc. Japan 59, 1357 (1990). 

Sacra1995 P. D. Sacramento, J. Phys. C 7, 143 (1995). 



Frahm1990 H. Frahm and V. E. Korepin, Phys. Rev. B 42, 10553 (1990); ibid 43, 5653 (1991).


Woyna1987  F. Woynarovich and H. P. Eckle, J. Phys. A 20, L443 (1987); F. Woynarovich, J. Phys. A 22 4243 (1989).


Ogata1990 M. Ogata and H. Shiba, Phys. Rev. B 41, 2326 (1990); M. Ogata, T. Sugiyama and H. Shiba, Phys. Rev. B
43, 8401 (1991).

Parola1990 A. Parola and S. Sorella, Phys. Rev. Lett. 64, 1831 (1990); Phys. Rev. Lett. 45, 13156 (1992); Phys. Rev. B 57, 6444 (1997).

Penc1995 K. Penc, F. Mila and H. Shiba, Phys. Rev. Lett. 75, 894 (1995);
K. Penc, K. Hallberg, F. Mila and H. Shiba, Phys. Rev. Lett. 77, 1390 (1996);
Phys. Rev. B 55, 15475 (1997).

Gohmann1998 F. Göhmann, A. R. Its and V. E. Korepin, Phys. Lett. A 249, 117 (1998);
F. Göhmann and V. E. Korepin, Phys. Lett. A 260, 516 (1999).


Ess2000 E. Jeckelmann, F. Gebhard and F. H. L. Essler, Phys. Rev. Lett. 85, 3910 (2000);
F. H. L. Essler and A. M. Tsvelik, Phys. Rev. B 65, 115117 (2002);
D. Controzzi and F. H. L. Essler, Phys. Rev. B 66, 165112 (2002);
B. Doyon and S. Lukyanov, Nucl. Phys. B 644, 451 (2002);
F. H. L. Essler and A. M. Tsvelik, Phys. Rev. Lett. 90, 126401 (2003).


Bogoliubov1988 N. M. Bogoliubov and V. E. Korepin, Mod. Phys. Lett. B 1, 349 (1988);
Teor. i Mat. Fiz. 82, No. 3, 331 (1990);
Int. J. Mod. Phys. B 3, 427 (1994).


DMRG A. E. Feiguin, and F. Heidrich-Meisner, Phys. Rev. B 76, 220508(R) (2007);
A. Lüscher, R. M. Noack and A. M. Läuchli, Phys. Rev. A 78, 013637 (2008);
M. Rizzi, M. Polini, M. A. Cazalilla, M. P. Tosi and R. Fazio, Phys. Rev. B 77, 245105 (2008);
M. Tezuka and M. Ueda, Phys. Rev. Lett. 100, 110403 (2008);
M. Tezuka and M. Ueda, New J. Phys. 12, 055029 (2010).


QMC G. G. Batrouni, M. H. Huntley, V.G. Rousseau and R. T. Scalettar, Phys. Rev. Lett. 100, 116405 (2008);
S. K. Baur, J. Shumway and E. J. Mueller, Phys. Rev. A 81, 033628 (2010);
M. J. Wolak, V. G. Rousseau, C. Miniatura, B. Gémaud, R. T. Scalettar and G. G. Batroun, Phys. Rev. A 82, 013614 (2010).

SC S. Cheng, Y.-Z. Jiang, Y.-C. Yu, M. T. Batchelor and X.-W. Guan,  
Nucl. Phys. B 929, 353 (2018).

Klumper1996 A. Klümper and R.Z. Bariev, Nucl. Phys. B 458, 623 (1996);
G. JüKttner, A. Klümper and J. Suzuki, Nucl. Phys. B 522, 471 (1998).

Yang1969 C. N. Yang and C. P. Yang, J. Math. Phys. 10, 1115 (1969).

Koma1987 T. Koma, Prog. Theor. Phys. 78, 1213 (1987); ibid 81, 783 (1989).






Cardy1984 J. L. Cardy, J. Phys. A 17, L385 (1984); Nucl. Phys. B 270 186, (1986).
Blote1986 H. W. Blöte, J. L. Cardy and M. P. Nightingale, Phys. Rev. Lett. 56, 742 (1986).
Aff1986 I. Affleck, Phys. Rev. Lett. 56, 746 (1986).
Belavin1984 A. A. Belavin, A. M. Polyakov and A. B. Zamolodchikov, Nucl. Phys. B 241, 333 (1984).
ZhaoE:2008E. Zhao and W. Vincent Liu, Phys. Rev. A 78, 063605 (2008).

Kajala:2011J. Kajala, F. Massel and P. Törmä, Phys. Rev. A 84, 041601 (R) (2011).


Singha2011 A. Singha, M. Gibertini, B. Karmakar, S. Yuan, M. Polini, G. Vignale, M. I. Katsnelson, A. Pinczuk, L. N. Pfeiffer, K. W. West and V. Pellegrini, Science  332, 1176 (2011).

Hart2015 R. A. Hart, P. M. Duarte, T.-L. Yang, X. Liu, T. Paiva, E. Khatami, R. T. Scalettar, N. Trivedi, D. A. Huse and R. G. Hulet, Nature 519, 211 (2015).

Greif2016 D. Greif, M. F. Parsons, A. Mazurenko, C. S. Chiu, S. Blatt, F. Huber, G. Ji and M. Greiner, Science 351, 953 (2016).

Parsons2016 M. F. Parsons, A. Mazurenko, C. S. Chiu, G. Ji, D. Greif and M. Greiner, Science 353, 1253 (2016).

Lawrence2016 L. W. Cheuk, M. A. Nichols, K. R. Lawrence, M. Okan, H. Zhang, E. Khatami, N. Trivedi, T. Paiva, M. Rigol and M. W. Zwierlein, Science 353, 1260 (2016).

Boll2016 M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo, L. Pollet, I. Bloch and C. Gross, Science 353, 1257 (2016).



Zhang:2015 R. Zhang, Y. Cheng, H. Zhai and P. Zhang, Phys. Rev. Lett. 115, 135301 (2015).

Mazurenko:2017 A. Mazurenko,  C.  S. Chiu, G. Ji, M. F. Parsons, M.  Kanász-Nagy, R.  Schmidt, F.  Grusdt, E.  Demler, D.  Greif and  M.  Greiner,
Nature 545, 462 (2017).

Kono:2015 Y. Kono, T. Sakakibara, C. P. Aoyama, C. Hotta, M. M. Turnbull,
C. P. Landee and Y. Takano,

  Phys. Rev. Lett. 114, 037202 (2015).

Krellner:2016C. Krellner, S. Lausberg, A. Steppke, M. Brando, L. Pe- drero, H. Pfau, S. Tencé, H. Rosner, F. Steglich and C Geibel, New J. Phys. 13, 103014 (2011).

Carmelo:1992J. M. P. Carmelo, P. Horsch and A. A. Ovchinnikov, Phys. Rev. B 45, 7899 (1992).

Shaginyan:2016V. R. Shaginyan, V. A. Stephanovich, K. G. Popov, E. V. Kirichenko and S. A. Artamonov, Ann. Phys. (Berlin) 528, 483 (2016).

short S. Cheng, Y.-Z. Jiang, M. T. Batchelor and X.-W. Guan, “FFLO correlation and free fluids in the one-dimensional attractive Hubbard model", arXiv:1708.07776.

Giamarchi:2004T. Giamarchi, Quantum Physics in one  dimension (Oxford University Press, Oxford, 2004).

Matveev:2008K. A. Matveev and A. Furusaki, Phys. Rev. Lett. 101, 170403 (2008).

YCY2016 Y.-C. Yu, Y.-Y. Chen, H.-Q. Lin, R. A. Roemer and X.-W. Guan, Phys. Rev. B 94, 195129 (2016).



Fisher1989 M. P. A. Fisher, P. B. Weichman, G. Grinstein and D. S. Fisher, Phys. Rev. B 40, 546 (1989).

note1 The sound velocity here is accessible by using the formula v_r=√(L/m_p  n_r   r∂^2 E_r/∂ L^2),
where the particle mass has been rescaled as m_p=1/2 and E_r=e_r L. See e_r in <ref>.

Sachdev S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, 1999).

Ho2010 Q. Zhou and T.-L. Ho, Phys. Rev. Lett. 105, 245702 (2010).

Hazzard2011 K. R. A. Hazzard and E. J. Mueller, Phys. Rev. A 84, 013604 (2011).


Yang:2017B. Yang, Y.-Y. Chen, Y.-G. Zheng, H. Sun, H.-N. Dai, X.-W. Guan, Z.-S. Yuan and J.-W. Pan,  Phys. Rev. Lett. 119, 165701 (2017)

He:2017 F. He, Y.-Z. Jiang, Y.-C. Yu, H.-Q. Lin and X.-W. Guan,  arXiv:1702.05903.

Rozhkov:2014A. V. Rozhkov, Phys. Rev. Lett. 112, 106403 (2014).

Lebed:2015 A. G. Lebed, Phys. Rev. Lett. 115, 157001 (2015).















