expressed as $\chi^{s}_{s}\propto T^{-\frac{1}{2}}(\Delta
B/T)^{-\frac{1}{2}}=(\Delta B)^{-\frac{1}{2}}$. Regarding the definition of
the critical exponent $\gamma$ with respect to the general form
$\chi\propto(g-g_{c})^{-\gamma}$, here $g$ is the driving parameter, we find
$\gamma=\frac{1}{2}$ for the second order derivatives of the free energy. On
the other hand, the correlation length can be expressed as $\xi\propto
T^{-\frac{1}{z}}$. In summary, at the critical point $\Delta B=0$, we also
find that thermodynamical properties $\frac{C_{v}}{T},\,\chi_{s},\,\xi$
satisfy the following scaling laws

$\displaystyle\frac{C_{v}}{T}\propto T^{\frac{d-z}{z}},\quad\chi_{s}\propto
T^{-\frac{\gamma}{vz}},\quad\xi\propto T^{-\frac{1}{z}},$ (195)

which signify the non-Fermi liquid behaviour at QCPcontinentino1989critical ;
continentino2005quantum ; continentino2017quantum .

###  III.6 III.6 Universal scaling functions at quantum criticality

In the previous subsection, we presented some analytic results for each phase
transition. We observe that the coefficients $a_{0}$ in (164) and $b_{0}$ in
(185) of free energies solely rely on the root densities (97) and (99) with
$a_{0}\approx\sigma(0),b_{0}\approx 2\pi\rho(\pi)$. Significantly, it is found
that the free energies are related to densities and dressed energies in
compact forms for different phase transitions

$\displaystyle\text{I-II: }\,f$ $\displaystyle=$ $\displaystyle
u+T^{\frac{3}{2}}\pi^{\frac{1}{2}}\rho(0)\left(\frac{\kappa^{{}^{\prime\prime}}(0)}{2}\right)^{-\frac{1}{2}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{-\frac{\kappa(0)}{T}}\right),$
(196) $\displaystyle\text{II-III: }\,f$ $\displaystyle=$ $\displaystyle
f_{0}+T^{\frac{3}{2}}\pi^{\frac{1}{2}}\rho(\pi)\left(\frac{-\kappa^{{}^{\prime\prime}}(\pi)}{2}\right)^{-\frac{1}{2}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{\kappa(\pi)}{T}}\right),$
(197) $\displaystyle\text{V-III: }\,f$ $\displaystyle=$ $\displaystyle
f_{0}+T^{\frac{3}{2}}\pi^{\frac{1}{2}}\sigma_{1}(0)\left(\frac{\varepsilon^{{}^{\prime\prime}}_{1}(0)}{2}\right)^{-\frac{1}{2}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{-\frac{\varepsilon_{1}(0)}{T}}\right),$
(198) $\displaystyle\text{II-IV: }\,f$ $\displaystyle=$ $\displaystyle
f_{0}-\frac{\pi
T^{2}}{6v_{c}}+T^{\frac{3}{2}}\pi^{\frac{1}{2}}\sigma_{1}(0)\left(\frac{\varepsilon^{{}^{\prime\prime}}_{1}(0)}{2}\right)^{-\frac{1}{2}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{-\frac{\varepsilon_{1}(0)}{T}}\right),$
(199) $\displaystyle\text{V-IV: }\,f$ $\displaystyle=$ $\displaystyle
f_{0}-\frac{\pi
T^{2}}{6v_{s}}+T^{\frac{3}{2}}\pi^{\frac{1}{2}}\rho(\pi)\left(\frac{-\kappa^{{}^{\prime\prime}}(\pi)}{2}\right)^{-\frac{1}{2}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{\kappa(\pi)}{T}}\right),$
(200)

where $f_{0}$ comes from the ground state, the terms with $T^{2}$ reflect the
contributions from the background parts. $\sigma_{1}(0)$ denotes the density
of length-$1$ spin strings at $\Lambda=0$, the second derivative
$\varepsilon^{{}^{\prime\prime}}_{1}(0)\equiv\left.\frac{d^{2}\varepsilon_{1}}{d\Lambda^{2}}\right|_{\Lambda=0}$,
and $\rho(0),\,\rho(\pi)$ denotes the charge density at $k=0,\,\pi$,
respectively. Similarly, for the charge dressed energy,
$\kappa^{{}^{\prime\prime}}(0)\equiv\left.\frac{d^{2}\kappa}{dk^{2}}\right|_{k=0},\kappa^{{}^{\prime\prime}}(\pi)\equiv\left.\frac{d^{2}\kappa}{dk^{2}}\right|_{k=\pi}$.
The polylog function $\mathrm{Li}_{\frac{3}{2}}$ represents the generating
function of free fermion criticality. The above scaling functions of the free
energy for different quantum phase transitions are valid for arbitrary
interaction strengths and fillings, revealing a microscopic origin of the
quantum phase transitions associated with the dressed energies. The functions
$-\varepsilon_{1}(0),\,-\kappa(0),\,\kappa(\pi)$ serve as criticality and
depend on the energy gaps away from the QCPs, i.e.

$\displaystyle-\varepsilon_{1}(0),\,-\kappa(0),\,\kappa(\pi)$
$\displaystyle\approx$ $\displaystyle\alpha_{B}\Delta
B+\alpha_{\mu}\Delta\mu+\alpha_{u}\Delta u.$ (201)

The factors $\alpha_{(B,\mu,u)}$ represent the different transition paths in
the vicinities of QCPs driven by external fields. These expressions
(196)-(200) display concise and elegant configurations independent of specific
details for arbitrary filling and interaction strength, and can apply to other
models with second order phase transitions associated with the dynamical
critical experiment $z=2$ and correlation length exponent $\nu=1/2$. The
derivations for these formulas are given in LPG:2022 .

##  IV IV. Spin incoherent liquid

Although spin incoherent Luttinger liquid has been studied in literature
cheianov2004nonunitary ; fiete2004green ; fiete2007colloquium ;
fiete2005theory ; fiete2005transport , almost all those works are based on the
framework of bosonization. There still lacks a study of such novel phenomenon
from the Bethe ansatz perspective. In Section II, we used the variations of
$\eta$-pair and spin magnetizations $(\Delta\eta^{z},\Delta S^{z})$ to
characterize the fractional charge and spinon excitations. Such fractionalized
quasi-particles reveal fermionic nature of quasiparticles, forming the
Luttinger liquid. Meanwhile we found the only possible fractional spin
excitations which can lead to the spin incoherent liquid at low temperatures.
Besides the fractional excitations, here we present rigorous results of the
SILL in terms of specific heat, criticality and correlation function.

###  IV.1 IV.1 Thermodynamics in SILL

In the previous analysis given in section III.5, we observe that a crossover
region fanning out from the critical point does show the existence of the SILL
above the phase boundary of the TLL and up to a critical temperature. This
phenomenon can be revealed through thermodynamic quantities of the model, such
as specific heat. We now identify different energy scaling of the SILL from
the TLL and quantum criticality. We first analyze the variables
$\varepsilon_{1}(0)$ emerged in the polylog functions of II-IV transition in
(199),

$\varepsilon_{1}(0)=2B-\int_{-\pi}^{\pi}\mathrm{d}k\cos ka_{1}(\sin
k-\Lambda)T\ln(1+e^{-\frac{\kappa(k)}{T}})|_{\Lambda=0},$ (202)

where $\varepsilon_{1}(0)=0,<0,>0$ at exact critical point, at phase IV and
phase II, representing the distance away from the phase boundary of II-IV. In
this case the Fermi point of the charge $\kappa(Q)=0$ gives $2\cos
Q=-\mu-2u-B$, leading to $\kappa(k)=-2\cos k+2\cos Q$, here $Q$ is the Fermi
point in charge sector. Thus $\varepsilon_{1}(0)$ at zero temperature can be
simplified by

$\varepsilon_{1}(0)=2B-\frac{4u}{\pi}\int_{0}^{Q}\mathrm{d}k\frac{\cos
k}{u^{2}+\sin^{2}k}(\cos k-\cos Q)+O(T^{2})$ (203)

with the temperature term omitted. By further performing Taylor expansion
around the critical point $B_{c}$ for fixed $\mu$ and $u$, thus we have

$\displaystyle Q\approx Q_{c}+\frac{1}{2\sin Q_{c}}\Delta B$ (204)

with $Q_{c}=\arccos\left(-\frac{1}{2}(\mu_{c}+2u_{c}+B_{c})\right)$. After
some algebra, the quantity $\varepsilon_{1}(0)$ in the vicinity of the
critical point of phase transition II-IV can be found to be given by

$\displaystyle\varepsilon_{1}(0)$ $\displaystyle\approx$ $\displaystyle
2B-\frac{4u}{\pi}\int_{0}^{Q_{c}}\mathrm{d}k\frac{\cos
k}{u^{2}+\sin^{2}k}(\cos k-\cos Q_{c}+\sin Q_{c}(Q-Q_{c}))$ (205)
$\displaystyle\approx$ $\displaystyle 2\Delta
B-\frac{2u}{\pi}\int_{0}^{Q_{c}}\mathrm{d}k\frac{\cos
k}{u^{2}+\sin^{2}k}\Delta B$ $\displaystyle\approx$ $\displaystyle 2\Delta
B\left[1-\frac{1}{\pi}\arctan\left(\frac{\sin Q_{c}}{u}\right)\right],$

where the second line brings in the condition
$\varepsilon_{1}(0)|_{B,Q_{c}}=2\Delta B$ with limitation
$\varepsilon_{1}(0)|_{B_{c},Q_{c}}=0$. Therefore
$\alpha_{B}=-2\left[1-\arctan\left(\sin Q_{c}/u\right)/\pi\right]$ can be
obtained in terms of the definition
$-\varepsilon_{1}(0)\equiv:\alpha_{B}\Delta B$.

Figure 13: (a) 3D plot of specific heat in the $T-B-C_{v}$ coordinate at
$\mu=-2,u=1$ for IV-II phase transition. The blue circled symbols present the
critical temperatures determined by the maximum values of the specific heat
(equation (206)), consistent with numerical simulation. The red lines mark the
TLL phase boundaries below which the specific heat show a linear temperature
dependent. Crossover regimes between the blue circled symbols and the red
lines denote the SILL phase (on the left) and the spin-gapped phase (on the
right). (b) and (c) show the specfic heat v.s. temperature near the phase
transition IV-II for the magnetic field is less than (greater than) critical
magnetic field $B_{c}$, respectively. The temperature-independent region and
the temperature-square-dependent region in the ratio $C_{v}/T$ can be visible,
charactering the thermodynamics of the TLL and SILL, respectively.

On the other hand, in the quantum critical region, heat capacity satisfies the
universal scaling form

$C_{v}/T=c_{0}+c_{1}T^{-1/2}\left[\frac{3}{4}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{x}\right)-x\mathrm{Li}_{\frac{1}{2}}\left(-\mathrm{e}^{x}\right)+x^{2}\mathrm{Li}_{-\frac{1}{2}}\left(-\mathrm{e}^{x}\right)\right]+O((\Delta
B/T)^{5/2}),$ (206)

where $c_{0}$ is the zero temperature background, $c_{1}$ a coefficient
depending on the transition point, explicitly see (169), and
$x=-\varepsilon_{1}(0)/T=\alpha_{B}\Delta B/T$. A brief discussion about
interaction-driven case is given in LPG:2022 . One characteristic of heat
capacity is that it displays bimodal structures around QCP. Thus it is
efficient to mark the QC boundaries in terms of the maxima points given by
heat capacity. The local maxima can be determined by $\partial C_{v}/\partial
B=0$, i.e.

$\frac{1}{4}\mathrm{Li}_{\frac{1}{2}}\left(-\mathrm{e}^{x}\right)-x\mathrm{Li}_{-\frac{1}{2}}\left(-\mathrm{e}^{x}\right)-x^{2}\mathrm{Li}_{-\frac{3}{2}}\left(-\mathrm{e}^{x}\right)=0,$
(207)

which gives two solutions $x_{1}=-1.5629,\,x_{2}=3.6205$. Figure 13 (a) shows
an overall behavior of quantum criticality in the vicinity of the IV-II phase
transition in the $T-B-C_{v}$ coordinate. The blue circled symbols denote the
maxima of specific heat from analytical results (206), showing a good
agreement between numerical TBA equations. In this critical region,
$T\gg\Delta B=B-B_{c}$, other thermodynamic properties also show universal
scaling behaviour given by equations (165)-(169).

As temperature decreases gradually from the QC part to a certain extent, TLL
regions appear, see the areas below the red lines in figure 13 (a). Around the
critical point $B_{c}=0.55$, there emerges two TLLs. In the region with
$B<B_{c}$, denoted as TLLSC, the system lies in phase IV with spin and charge
degrees of freedom coexisting. More than $B_{c}$, the system lies in phase II
with only charge degree, denoted as TLLC. In the TLL region, the specific heat
$C_{v}$ is linearly dependent on $T$, also see discussion in Section III. In
the crossover region between QC and the TLLC phase, the spin sector is gapped.
By utilizing the asymptotic behaviour of polylog function mohankumar2007two
and expanding (169), the specific heat is given by

$C_{v}\approx\frac{\pi
T}{3v_{c}}+\frac{3\pi^{\frac{1}{2}}a_{0}}{4}\left(\eta_{1}\right)^{-\frac{1}{2}}T^{\frac{1}{2}}\mathrm{e}^{\alpha_{B}\Delta
B/T}+O\left(\Delta B\mathrm{e}^{\alpha_{B}\Delta B/T}\right).$ (208)

By comparison, in the crossover region between QC and the TLLSC phase, the
SILL lies in the temperature range $E_{s}\sim k_{F}v_{s}\ll k_{B}T\ll
E_{c}\sim k_{F}v_{c}$. This area corresponds to SILL (spin-incoherent
Luttinger liquid) cheianov2004nonunitary ; fiete2004green ;
fiete2007colloquium ; fiete2005theory ; fiete2005transport  with specific heat

$C_{v}\approx\frac{\pi
T}{3v_{c}}+\frac{\pi^{2}a_{0}\left(\eta_{1}\right)^{-\frac{1}{2}}(-\varepsilon_{1}(0))^{-\frac{1}{2}}T}{3}\left[1+\frac{21\pi^{2}}{40}(-\varepsilon_{1}(0))^{-2}T^{2}\right]+O(T^{4}),$
(209)

which manifests a gas-liquid co-existence. We note that the coefficient before
the square bracket should be $\pi T/(3v_{s})$. The proof is given in
following. The spin dressed energy around QCP is written as
$\varepsilon_{1}(\Lambda)=D_{1}\Lambda^{2}+D_{2}$. Thus for the SILL region,
the spin Fermi point $A$ is related with
$\varepsilon_{1}(A)=D_{1}A^{2}+D_{2}=0$, resulting in
$A=(-D_{2}/D_{1})^{1/2}$. On the other hand, by the definition of $D_{1,2}$
given in (139) and (140), we have
$D_{1}=\varepsilon^{{}^{\prime\prime}}_{1}(0)/2=\eta_{1},D_{2}=\varepsilon_{1}(0)$.
Therefore the Fermi point $A$ is expressed as
$A=(-\varepsilon_{1}(0)/\eta_{1})^{1/2}$. Using the definition of $a_{0}$, see
(164), we have $a_{0}=\sigma_{1}(0)\approx\sigma_{1}(A)$ and
$\varepsilon_{1}^{{}^{\prime}}(A)=2\eta_{1}A$, which renders

$\displaystyle\frac{\pi^{2}a_{0}\left(\eta_{1}\right)^{-\frac{1}{2}}(-\varepsilon_{1}(0))^{-\frac{1}{2}}T}{3}$
(210) $\displaystyle=$
$\displaystyle\frac{\pi^{2}\sigma_{1}(A)T\left(\eta_{1}\right)^{\frac{1}{2}}}{3}\frac{1}{\eta_{1}(-\varepsilon_{1}(0))^{\frac{1}{2}}}$
$\displaystyle=$
$\displaystyle\frac{\pi^{2}\sigma_{1}(A)T}{3\eta_{1}A}=\frac{\pi\cdot
2\pi\sigma_{1}(A)T}{3\varepsilon_{1}^{{}^{\prime}}(A)}$ $\displaystyle=$
$\displaystyle\frac{\pi T}{3v_{s}}.$

This immediately gives an universal thermodynamic relation of the SILL

$\displaystyle C_{v}\approx\frac{\pi
T}{3}\left(\frac{1}{v_{c}}+\frac{1}{v_{s}}\right)+\frac{7\pi^{3}T^{3}}{40v_{s}(-\varepsilon_{1}(0))^{2}}+O(T^{4}).$
(211)

Furthermore, in figure 13 (b) and (c), we plot the specific heat below and
above the QCP for different values of magnetic fields. It is obvious to see
the region of the linear-temperature dependent specific heat, a crossover
region of the SILL with both the linear- and cubic- temperature-dependent
specific heat. The latter marks the crossover region $E_{s}\sim k_{F}v_{s}\ll
k_{B}T\ll E_{c}\sim k_{F}v_{c}$, showing the thermodynamic behaviour of the
SILL.

###  IV.2 IV.2 Correlation functions in SILL

The SILL regime with $k_{F}v_{s}<k_{B}T\ll k_{F}v_{c}$ lies between the
boundaries of TLL and QC. In this region, the spin degree of freedom behaves
like hot spins, whereas the charge still behaves as a collective motion of
bosons. As a result, the SILL largely behaves like a spin and charge decohered
liquid, i.e., possessing solely a propagating charge mode but not a spin mode,
see cheianov2004nonunitary ; fiete2004green ; fiete2007colloquium ;
fiete2005theory ; fiete2005transport . In the spin sector, the magnetic
exchange energy is lower than the Fermi energy, resulting in spin thermally
excited with equal probability fiete2005transport . In the strong coupling
regime, the spin degrees of freedom shows a spin dynamics of Heisenberg chain
with nondispersive spinon band due to small effective exchange coupling
$J=4t^{2}/U$. Whereas the charge acts as noninteracting fermions with
dispersive spectrum $\kappa(k)=-2t\cos(k)$. This SILL theory can also be
captured in the excitation spectrum in the low density regime in figure 3 (b).
The concept of SILL is helpful to explain the appearance of conductance
plateau in a quantum wirehew2008spin . When interaction increases or in the
vicinity of QCP, the spin velocity progressively dwindles to zero. This means
that the spin sector loses dynamics, indicating the emergence of the SILL.

In this regime $E_{s}\ll T\ll E_{c}$, the spin sector is totally thermal
averaged and equally excited, whereas the charge remains at relevant low-
energy performance, rendering the correlation functions independent of
temperature fiete2007colloquium . We note that the TLL theory has its own
applicable condition $T\ll E_{c},\,E_{s}$. Whereas for the region $E_{s}\ll
T\ll E_{c}$, the energy scales of charge and spin degrees of freedom can be
dealt with separately cheianov2004nonunitary . Under such a circumstance, the
finite-temperature correlation functions in terms of conformal field theory
schulz1990correlation ; frahm1990critical ; frahm1991correlation  still remain
valid for the charge and spin degrees of freedom operating with different
limits, i.e.

$\displaystyle|x\pm\mathrm{i}v_{c}t|\ll
v_{c}/T,\qquad|x\pm\mathrm{i}v_{s}t|\gg v_{s}/T.$ (212)

This is essential to capture the asymptotic behaviour of the SILL. Here we
would like to mention that the remaining temperature term in the spin sector
can be replaced by the typical energy scale of spin $T\sim E_{s}\sim
J\sim\left(k_{F\uparrow}+k_{F\downarrow}\right)/2\cdot v_{s}\equiv
k_{F}v_{s}$. With the help of the finite temperature asymptotics of
correlation functions under the condition (212), the two-point correlation
functions of prime fields can be obtained as

$\displaystyle\left\langle\phi(x,t)\phi(0,0)\right\rangle=\sum
A(D_{c},D_{s},N^{\pm}_{c},N^{\pm}_{s})\text{exp}(-2\mathrm{i}D_{c}k_{F,\uparrow}x)\text{exp}(-2\mathrm{i}(D_{c}+D_{s})k_{F,\downarrow}x)$
$\displaystyle\times\frac{1}{(x-\mathrm{i}v_{c}t)^{2\Delta^{+}_{c}}(x+\mathrm{i}v_{c}t)^{2\Delta^{-}_{c}}}\left(2\pi\alpha
k_{F}\right)^{2\Delta^{+}_{s}+2\Delta^{-}_{s}}\mathrm{e}^{-\pi\alpha\left(2\Delta^{+}_{s}+2\Delta^{-}_{s}\right)k_{F}x},$
(213)

where the conformal dimensions in gapless phases read

$\displaystyle 2\Delta^{\pm}_{c}(\Delta\bm{N},\bm{D})$ $\displaystyle=$
$\displaystyle\left(\mathrm{Z}_{cc}D_{c}+\mathrm{Z}_{sc}D_{s}\pm\frac{\mathrm{Z}_{ss}\Delta
N_{c}-\mathrm{Z}_{cs}\Delta
N_{s}}{2\mathrm{det}\mathrm{Z}}\right)^{2}+2N^{\pm}_{c},$ $\displaystyle
2\Delta^{\pm}_{s}(\Delta\bm{N},\bm{D})$ $\displaystyle=$
$\displaystyle\left(\mathrm{Z}_{cs}D_{c}+\mathrm{Z}_{ss}D_{s}\pm\frac{\mathrm{Z}_{cc}\Delta
N_{s}-\mathrm{Z}_{sc}\Delta
N_{c}}{2\mathrm{det}\mathrm{Z}}\right)^{2}+2N^{\pm}_{s}.$ (214)

Here $\mathrm{Z}$ is the dressed charge matrix given in (115), while
$N_{\alpha}^{\pm},\Delta\vec{N},\vec{D}$ are related to the three types of
excitations: adding particle from left and right Fermi points, changing total
particle number, moving of the particle center displacement, also see
essler2005one . From equation (213), the spin-spin correlation in the spin
sector displays exponential decay as a function of distance, while the
correlation in the charge sector shows a power-law decay. Explicitly, based on
the correlation functions obtained from the finite temperature CFT, we may
calculate various two-point correlations of field operators close to critical
field $B_{c}$ in SILL regime, including the single particle Green’s function

  1. (1)

$G^{\uparrow}$

$G^{\uparrow}_{B\rightarrow
B_{c}}\sim\frac{\text{exp}(-\mathrm{i}k_{F,\uparrow}x)}{(x-\mathrm{i}v_{c}t)^{1-\frac{2}{\pi}\sqrt{1-\frac{B}{B_{c}}}}}\left(2\pi\alpha
k_{F}\right)^{\frac{1}{2}-\frac{1}{\pi}\sqrt{1-\frac{B}{B_{c}}}}\mathrm{e}^{-\pi\alpha\left(\frac{1}{2}-\frac{1}{\pi}\sqrt{1-\frac{B}{B_{c}}}\right)k_{F}x}+h.c..$
(215)

This result was derived for the region close to the critical field $B_{c}$,
corresponding to the phase transition II-IV. Similarly, for other correlations

  2. (2)

$G^{\downarrow}$

$\displaystyle G^{\downarrow}_{B\rightarrow B_{c}}$ $\displaystyle\sim$
$\displaystyle\frac{\text{exp}(-\mathrm{i}k_{F,\downarrow}x)}{(x-\mathrm{i}v_{c}t)^{\frac{1}{4}+\frac{1}{\pi}\sqrt{1-\frac{B}{B_{c}}}}(x+\mathrm{i}v_{c}t)^{\frac{1}{4}-\frac{1}{\pi}\sqrt{1-\frac{B}{B_{c}}}}}$
(216) $\displaystyle\times\left(2\pi\alpha
k_{F}\right)^{1-\frac{2}{\pi}\sqrt{1-\frac{B}{B_{c}}}}\mathrm{e}^{-\pi\alpha\left(1-\frac{2}{\pi}\sqrt{1-\frac{B}{B_{c}}}\right)k_{F}x}+h.c..$

  3. (3)

$G^{n}$

$\displaystyle G^{n}_{B\rightarrow B_{c}}$ $\displaystyle\sim$ $\displaystyle
n^{2}+\left(\text{exp}(2\mathrm{i}k_{F,\uparrow}x)+\text{exp}(-2\mathrm{i}k_{F,\uparrow}x)\right)\frac{1}{|x-\mathrm{i}v_{c}t|^{2-\frac{8}{\pi}\sqrt{1-\frac{B}{B_{c}}}}}$
(217) $\displaystyle\times\left(2\pi\alpha
k_{F}\right)^{2-\frac{4}{\pi}\sqrt{1-\frac{B}{B_{c}}}}\mathrm{e}^{-\pi\alpha\left(2-\frac{4}{\pi}\sqrt{1-\frac{B}{B_{c}}}\right)k_{F}x}$
$\displaystyle+\left(\text{exp}(2\mathrm{i}k_{F,\downarrow}x)+\text{exp}(-2\mathrm{i}k_{F,\downarrow}x)\right)\left(2\pi\alpha
k_{F}\right)^{2-\frac{4}{\pi}\sqrt{1-\frac{B}{B_{c}}}}\mathrm{e}^{-\pi\alpha\left(2-\frac{4}{\pi}\sqrt{1-\frac{B}{B_{c}}}\right)k_{F}x}$
$\displaystyle+\left(\text{exp}(2\mathrm{i}\left(k_{F,\uparrow}+k_{F,\downarrow}\right)x)+\text{exp}(-2\mathrm{i}\left(k_{F,\uparrow}+k_{F,\downarrow}\right)x)\right)\frac{1}{|x-\mathrm{i}v_{c}t|^{2}}.$

  4. (4)

$G^{\perp}$

$\displaystyle G^{\perp}_{B\rightarrow B_{c}}$ $\displaystyle\sim$
$\displaystyle\left(\text{exp}(\mathrm{i}\left(k_{F,\uparrow}+k_{F,\downarrow}\right)x)+\text{exp}(-\mathrm{i}\left(k_{F,\uparrow}+k_{F,\downarrow}\right)x)\right)\frac{1}{|x-\mathrm{i}v_{c}t|^{\frac{1}{2}}}$
(218) $\displaystyle\times\left(2\pi\alpha
k_{F}\right)^{\frac{1}{2}+\frac{1}{\pi}\sqrt{1-\frac{B}{B_{c}}}}\mathrm{e}^{-\pi\alpha\left(\frac{1}{2}+\frac{1}{\pi}\sqrt{1-\frac{B}{B_{c}}}\right)k_{F}x}.$

  5. (5)

$G^{p}$

$\displaystyle G^{p}_{B\rightarrow B_{c}}$ $\displaystyle\sim$
$\displaystyle\text{exp}(-\mathrm{i}\left(k_{F,\uparrow}+k_{F,\downarrow}\right)x)\frac{1}{(x-\mathrm{i}v_{c}t)^{\frac{9}{4}}(x+\mathrm{i}v_{c}t)^{\frac{1}{4}}}$
(219) $\displaystyle\times\left(2\pi\alpha
k_{F}\right)^{\frac{1}{2}-\frac{3}{\pi}\sqrt{1-\frac{B}{B_{c}}}}\mathrm{e}^{-\pi\alpha\left(\frac{1}{2}-\frac{3}{\pi}\sqrt{1-\frac{B}{B_{c}}}\right)k_{F}x}+h.c..$

It should be noticed that although these results are obtained in a very rough
approximation, they capture the essential features of the SILL
fiete2007colloquium , i.e., spin-spin correlation decays exponentially while
the correlation in the charge degree of freedom behaves like spinless
noninteracting fermions. The field theory approach was given in
cheianov2004nonunitary ; fiete2004green ; fiete2007colloquium . The reason why
the crossover in the vicinity of phase transition between IV and V can not
develop such a similar concept as charge incoherent Luttinger liquid is that
near a half-filled lattice, the charge velocity monotonically and
exponentially tends to zero giamarchi2003quantum . This means subtlety
occurring near the Mott phase transition due to rapid vanishing charge
collective mode. This limitation is called as the holon confinement. In the
next section, we will develop a new concept, Contact susceptibility, to study
the Mott phase transition.

##  V V. Contact and Contact susceptibility

In this section, we focus on an experimentally measurable quantity double
occupancy campo2012thermal ; van2018competing , as well as the associated
Contact, revealing the competition between thermal fluctuations and quantum
fluctuations from different sources, i.e., external fields and interaction.
The double occupancy serves as an efficient instrument to demarcate phase
transition, especially the Mott phase transition in the extended Hubbard model
with long range interaction glocke2007half ; de2011thermodynamics . On the
other hand, the partial wave Contact, which was first proposed in ultracold
Fermi gas tan2008energetics ; ZhangSZ:2009 , has become an important theme in
the study of ultracold atoms Braaten:2009 ; Stewart:2010 ; Yu:2015 ;
Yoshida:2015 ; Cui:2016 ; chen2014critical . Here we will introduce the
concept of Contact susceptibilities with respect to the temperature, magnetic
field and chemical potential, and investigate applications of these Contact
susceptibilities. We will show that the susceptibilities build up a general
connection between interaction-driven quantum criticality and the phase
transitions induced by external fields. Using these relations, we will also
obtain caloric effect in interaction-driven quantum refrigeration and scaling
laws at quantum criticality.

Figure 14: (a): Double occupancy maps out the phase diagram of the 1D Hubbard
model plotted. The white dot-dashed lines show the zero temperature phase
boundaries. While the thick dashed lines show the contour lines of $d=0.04$
and $d=0.07$, respectively. At the critical points, the double occupancy
suddenly changes. Parameter setting for numerics is the same as figure (8)
(a), i.e. $T=0.005$ and $u=1$. (b): Double occupancy versus magnatization in
Mott phase for various coupling strengths. The numerics from the TBA equations
confirms the accuracy of the analytical result (220) with the Contact given by
(222).

###  V.1 V.1 Double occupancy and Mott phase

In the Hubbard model, the lattice filling parameter, interaction strength and
external fields can drive different phase transitions giamarchi1991umklapp ;
giamarchi1997mott . In contrast to the dilute limit case, the Mott insulator
phase induced by interaction is much less understood. When the interaction
strength increases up to a critical value $u_{c}$glocke2007half ;
de2011thermodynamics ; giamarchi2003quantum  in the canonical ensemble, the
system reaches a Mott insulator state in the extended Hubbard model. The half-
filled phase can be delineated by the Luther-Emery liquid with one gapless and
one gapped sectors luther1974backward , in which the excitations comply with
bosonic and fermionic statistics, respectively. In order to explore the
magnetic order and detect the Mott phase transition, one can introduce the
double occupancy $d=\frac{1}{N}\sum_{i}\left\langle
n_{i,\uparrow}n_{i,\downarrow}\right\rangle$, which depicts the probability of
two electrons with opposite spin occupying a single lattice site with
potential energy $E_{\mathrm{pot}}=ud$. The double occupancy has been invested
extensively in theoretical and experimental research on the Hubbard model at
half filling. It appears to show discontinuity as the interaction approaches
critical coupling in interaction-induced Mott transition (note that $u_{c}=0$
for this model in the canonical ensemble in which this phenomenon is
unconspicuous) and exhibits nonmonotonic behavior in the half-filled state as
temperature varies, similar to the Pomeranchuk effect (the melting pressure of
liquid Helium-3 shows a trend of first decreasing and then increasing with
temperature) glocke2007half ; de2011thermodynamics ; wietek2021mott . It also
reveals the competing physics from the charge and spin fluctuations. For
doping-induced Mott transition, $d$ is significantly suppressed in the area of
$\mu<u$ and generates a pronounced signal when density exceeds unity
de2011thermodynamics . It is useful to locate the Mott transition point
through the detection of double occupancy.

In the canonical ensemble, the double occupancy can be obtained from the free
energy $f$ by

$\displaystyle d$ $\displaystyle=$ $\displaystyle\frac{1}{4}\frac{\partial
f}{\partial u}-\frac{1}{4}+\frac{n_{c}}{2},$ (220)

in which the last two terms $-1/4+n_{c}/2$ stem from the extra terms $-2uN+uL$
in the Hamiltonian. We can define $C=\partial f/\partial u$ to be the lattice
version of Contact $C$,

$\displaystyle C=\partial f/\partial u=4d-2n_{c}+1,$ (221)

which is analogous to the Tan’s Contact tan2008energetics ; chen2014critical
in the continuous systems of ultracold atoms. Regarding the role that the
double occupancy can reflect the phase information, we expect that $d$
distinguishes different phases with and without internal degrees of freedom,
see figure 14 (a). $d$ always vanishes for phases I, II, and III and has no
demarcation lines between these phases. The boundaries of II and IV or III and
V are obvious, while for IV and V the inflection point of the contour line
marks the phase transition. By the fact that the double occupancy essentially
reflects charge fluctuations at Mott insulator accompanying with the existence
of the antiferromagnetic order, we demonstrate the roles of magnetic field $B$
and interaction strength $u$ on $d$ in figure (14) (b). It is clear to observe
that repulsive interaction always lowers the possibility of double occupancy
which is quite intuitive. When the interaction increases, the electrons with
different spins are repel each other more strongly and the system becomes more
incompressible. The infinity coupling $u\to\infty$ makes the particles fully
localized. The decrease caused by magnetization originates from the tendency
that the magnetic field tends to align the spins of electrons in the field
direction. Through the derivative of the background term in (185) with respect
to $u$, the exact analytic expression of the Contact $C_{0}$ at Mott phase V
of ground state can be obtained as

$C_{0}=-1-\frac{4}{3}\frac{\partial}{\partial
u}\left(\lambda_{1}\eta_{1}\right)A^{3}-4\lambda_{1}\eta_{1}A^{2}\frac{\partial
A}{\partial u}-\frac{2}{15}\frac{\partial}{\partial
u}\left(\lambda_{2}\eta_{1}\right)A^{5}-\frac{2}{3}\lambda_{2}\eta_{1}A^{4}\frac{\partial
A}{\partial u},$ (222)

where $\partial A/\partial u$ is given in (231) below. We plotted the double
occupancy through (222) (solid lines) in figure 14 (b), showing a good
agreement between this analytical result and the numerical calculation from
the TBA equations. Based on the relation from (221), it is more essential to
investigate the variations of the Contact with respect to the changes of the
external fields and temperature.

Figure 15: Universal Scaling behaviours of susceptibilities: (a) Contact
susceptibilities with respect to chemical potential, magnetic field and
temperature v.s. magnetic field at the quantum criticality of phase transition
IV-V. (b) Contact susceptibilities v.s. the argument $(B-B_{c})/T$ show
universal scaling forms at quantum criticality. Parameter setting reads $u=1,$
and $\mu=-0.8272$. Numerical results obtained from the TBA equations (symbols)
confirm the scaling functions of Contact susceptibilities (226)-(229) (solid
lines).

###  V.2 V.2 Contact susceptibilities and Mott phase transition

In comparison to the double occupancy $d$, the Contact $C$ is more essential
to capture many-body effects induced by interaction strength. Using the
thermal potential $f=e-\mu n_{c}-2Bm-Ts-uC$ and the Maxwell relations in its
derivatives with respect to the temperature, chemical potential and magnetic
fields, we may build up general relations between Contact susceptibilities and
interaction-driven variations of density, magnetization and entropy:

$\displaystyle\frac{\partial n_{c}}{\partial u}$ $\displaystyle=$
$\displaystyle-\frac{\partial C}{\partial\mu},$ (223)
$\displaystyle\frac{\partial m}{\partial u}$ $\displaystyle=$
$\displaystyle-\frac{\partial C}{\partial(2B)},$ (224)
$\displaystyle\frac{\partial s}{\partial u}$ $\displaystyle=$
$\displaystyle-\frac{\partial C}{\partial T}.$ (225)

We prove that these Contact susceptibilities will provide striking features of
interaction effect in the thermodynamics of the model.

Building on the Contact susceptibilities, we further show that the Contact
susceptibilities exhibit universal scaling behaviour in the quantum critical
region chen2014critical . The Mott insulator phase with average one electron
occupying one site is of particular interest in many-body physics. Here using
the free energy (185) for the transition from IV to V, we obtain the Contact
susceptibilities

$\displaystyle C$ $\displaystyle=$
$\displaystyle-1-\frac{4}{3}f_{11}A^{3}-4\lambda_{1}\eta_{1}A^{2}\frac{\partial
A}{\partial
u}-\frac{2}{15}f_{21}A^{5}-\frac{2}{3}\lambda_{2}\eta_{1}A^{4}\frac{\partial
A}{\partial
u}+\frac{T^{\frac{1}{2}}b_{0}b_{1}}{2\pi^{\frac{1}{2}}}f_{\frac{1}{2}},$ (226)
$\displaystyle\frac{\partial C}{\partial\mu}$ $\displaystyle=$
$\displaystyle-\frac{b_{0}b_{1}}{2\pi^{\frac{1}{2}}T^{\frac{1}{2}}}f_{-\frac{1}{2}},$
(227) $\displaystyle\frac{\partial C}{\partial(2B)}$ $\displaystyle=$
$\displaystyle\frac{\partial C_{0}}{\partial
2B}-\frac{b_{0}b_{1}}{4\pi^{\frac{1}{2}}T^{\frac{1}{2}}}\left(1-\frac{4A_{c}}{\pi
u+A_{c}}\right)f_{-\frac{1}{2}},$ (228) $\displaystyle\frac{\partial
C}{\partial T}$ $\displaystyle=$
$\displaystyle\frac{b_{0}b_{1}}{4\pi^{\frac{1}{2}}T^{\frac{1}{2}}}f_{\frac{1}{2}}+\frac{b_{0}b_{1}}{2\pi^{\frac{1}{2}}T^{\frac{3}{2}}}\left[\Delta\mu+\left(1-\frac{4A_{c}}{\pi
u+A_{c}}\right)\Delta B\right]f_{-\frac{1}{2}},$ (229)

where $C_{0}$ denotes the background terms of ground state, that is equation
(222). We have denoted the background of the magnetic susceptibility and other
functions as

$\displaystyle\frac{\partial C_{0}}{\partial 2B}$ $\displaystyle=$
$\displaystyle-4A^{2}\left(f_{11}+\frac{1}{6}f_{21}A^{2}\right)\frac{\partial
A}{\partial(2B)}-8\eta_{1}A\left(\lambda_{1}+\frac{1}{3}\lambda_{2}A^{2}\right)\frac{\partial
A}{\partial(2B)}\frac{\partial A}{\partial u},$ (230)
$\displaystyle-4\eta_{1}A^{2}\left(\lambda_{1}+\frac{1}{6}\lambda_{2}A^{2}\right)\frac{\partial^{2}A}{\partial(2B)\partial
u},$ $\displaystyle f_{11}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
u}\left(\lambda_{1}\eta_{1}\right),f_{21}=\frac{\partial}{\partial
u}\left(\lambda_{2}\eta_{1}\right),\qquad\frac{\partial A}{\partial
B}=-\frac{1}{\eta_{1}A\left(1+\frac{A}{\pi u}\right)},$ $\displaystyle b_{1}$
$\displaystyle=$ $\displaystyle\frac{\partial(2+C_{2})}{\partial
u}=-2-\frac{8}{\pi u\left(1+u^{2}\right)^{\frac{3}{2}}}A^{2}\frac{\partial
A}{\partial u}.$

Solving (179), we have

$\displaystyle\frac{\partial A}{\partial u}$ $\displaystyle=$
$\displaystyle\left(1+u^{2}\right)\frac{u-\left(1+u^{2}\right)^{\frac{1}{2}}}{A\left(1+\frac{A}{\pi
u}\right)}+\frac{3u}{2\left(1+u^{2}\right)}\frac{A\left[1+\frac{2}{9\pi
u}\left(4+\frac{1}{u^{2}}\right)A\right]}{1+\frac{A}{\pi u}},$ (231)
$\displaystyle\frac{\partial^{2}A}{\partial B\partial u}$ $\displaystyle=$
$\displaystyle\frac{1+\frac{2A}{\pi u}}{\eta_{1}A^{2}\left(1+\frac{A}{\pi
u}\right)^{2}}\frac{\partial A}{\partial u}-\frac{1}{\eta_{1}\pi
u^{2}\left(1+\frac{A}{\pi u}\right)^{2}}.$ (232)

In figure 15, we demonstrate scaling behaviour of Contact susceptibilities,
confirming the analytic expressions (226)-(229) (solid lines) with the
numerical data from TBA (symbols).

###  V.3 V.3 Interaction-driven quantum cooling

It is remarkable to observe that equation (225) essentially relates entropy,
temperature and interaction strength. This relation is of great importance for
interaction driven quantum cooling since the interaction can be tuned in cold
atom experiments via Feshbach resonance, and temperature plays a vital role in
controlling the entropy. Therefore an adiabatic cycle process can be used in
realization of quantum refrigeration in real physical systems. Inspired by
this, here we consider an isentropic process by ramping up or down the
interaction strength. Focusing on the $(T,u)$ coordinates for a fixed magnetic
field, we perform total derivatives on entropy $s$, i.e.,

$\mathrm{d}s=\frac{\partial s}{\partial u}\mathrm{d}u+\frac{\partial
s}{\partial T}\mathrm{d}T=0.$ (233)

Using equation (225) and the relation $\partial s/\partial T=C_{v}/T$, the
points on the isentropic line in the $(u,T)$ coordinates satisfy the relation:

$\frac{C_{v}}{T}\frac{\partial T}{\partial u}=\frac{\partial C}{\partial T}.$
(234)

Thus the interaction driven Grüneisen ratio $\Gamma_{\mathrm{int}}$ is related
to $\partial C/\partial T$ via the relation $\Gamma_{\mathrm{int}}=\partial
C/\partial T\cdot u/{C_{v}}$. Note that $\partial C/\partial T$ dramatically
changes near quantum phase transition (see equation (229) and figure 15),
where $C_{v}$ has a minimum (see figure 12 ) and the entropy has a maximum
(see figure 16). A good cooling effect is observed when the interaction drives
the system approaching a critical point. The lowest temperature point can be
estimated by $\partial C/\partial T=0$. To characterise the refrigerating
efficiency, we give a brief description on the lowest reachable temperature
during an adiabatic cooling cycle. From the results (229) of $\partial
C/\partial T$, the condition for determining the lowest temperature reads

$\frac{1}{2}\mathrm{Li}_{\frac{1}{2}}\left(-\mathrm{e}^{x}\right)-x\mathrm{Li}_{-\frac{1}{2}}\left(-\mathrm{e}^{x}\right)=0,$
(235)

where $x=\alpha_{u}\Delta u/T$. It is found that, to an good approximation,
$x\approx 1.3117$. Thus the extreme point can be determined. In figure 16, we
plot the isentropic lines near the phase boundaries of II-IV (a) and IV-V (b),
respectively. It is obvious that entropy shows a minimum in the quantum
critical regions. Away from QCP with temperature $T\ll|u-u_{c}|$, entropy
linearly depends on temperature, see the left part of figure 16 (a) and right
part of figure 16 (b). In figure 16 (b) we draw an Otto cycle for the
interaction-driven refrigeration process. The stages A, D lie around QCP,
whereas the stages B,C are located in the TLL area. There contains four steps
to cool the target material. For A $\rightarrow$ B, the working substance is
adiabatically ramped up from the target temperature $T_{\rm target}$ to the
nonthermal higher temperature stage B. Then through a hot isochore process B
$\rightarrow$ C, the working substance comes into contact with the ambient,
transferring heat to the high temperature source. While the temperature of
working substance reduces to the one at the thermal state C. Next, for the
isentrope process C $\rightarrow$ D, the working substance is adiabatically
ramped down to the low temperature stage D. This is an opposite process
contrast to the A $\rightarrow$ B. Finally, for the isochore process D
$\rightarrow$ A, the working substance contacts with the target object,
absorbing heat from the target material and reaching the thermal state A.
Consequently, the target object is cooled down by this cycle.

Figure 16: Plot of isentropic lines in $T-u$ plane for the interaction-driven
phase transitions II-IV (a) and IV-V (b) at $B=0.15$ and $\mu=-2.5$, of which
the temperature on the vertical axis is logarithmic. The purple dotted lines
in (b) denote the interaction driven Otto cycle, where the stages A and B (or
stages C, D) lie on the isentropic lines, see the main text.

Now let us determine the lowest temperature which can be reached through an
isentropic process indicated in the figure 16. From equations (196)-(200), the
phase II (V) contains one charge (spin) degrees of freedom. Consequently,
their entropy $s_{\mathrm{L1}}$ and $s_{\mathrm{L2}}$ are given by (61) and
(62), respectively, namely,

$\displaystyle s_{\mathrm{L1}}$ $\displaystyle\approx$ $\displaystyle\frac{\pi
T_{\mathrm{L1}}}{3v_{c}},$ (236) $\displaystyle s_{\mathrm{L2}}$
$\displaystyle\approx$ $\displaystyle\frac{\pi T_{\mathrm{L2}}}{3v_{s}}.$
(237)

Comparing isentropic lines with the same entropy for phase II, IV and V, for
example $s=0.00008$ in figure 16, the temperature of TLLC (phase II) is higher
than that of TLLS (phase V) since charge velocity $v_{c}$ changes faster than
spin velocity $v_{s}$ when the interaction is changed around the critical
point, i.e.,

$T_{\mathrm{L1}}>T_{\mathrm{L2}}.$ (238)

On the other hand, when it approaches the QCP, i.e., at the extreme low
temperature for each isentropic line, entropy have explicit expressions for
the transition II-IV and IV-V

$\displaystyle s_{\mathrm{s1}}$ $\displaystyle\approx$
$\displaystyle\lambda_{3}\pi^{1/2}\sigma_{1}(0)(\varepsilon^{{}^{\prime\prime}}_{1}(0)/2)^{-1/2}T_{c1}^{1/2},$
(239) $\displaystyle s_{\mathrm{s2}}$ $\displaystyle\approx$
$\displaystyle\lambda_{3}\pi^{1/2}\rho(\pi)(-\kappa^{{}^{\prime\prime}}(\pi)/2)^{-1/2}T_{c2}^{1/2},$
(240)

respectively, where
$\lambda_{3}=x\mathrm{Li}_{1/2}\left(-\mathrm{e}^{x}\right)-3/2\mathrm{Li}_{3/2}\left(-\mathrm{e}^{x}\right)\approx
1.3467$.

With the above analysis, we observe that the entropy shows a square root
dependence on the temperature at extreme point, it is proportional to
temperature in the Luttinger liquid. Therefore, considering an isentropic
cooling process through the ramping up or down in the $T-u$ plane around
critical phase transitions from II to IV or from V to IV, see figure 16, the
minimum temperatures can be reached

$\displaystyle\text{II-IV: }\frac{T_{c1}^{1/2}}{T_{L1}}$ $\displaystyle=$
$\displaystyle\frac{\pi^{1/2}(\varepsilon^{{}^{\prime\prime}}_{1}(0)/2)^{1/2}}{3\lambda_{3}v_{c}\sigma_{1}(0)},$
(241) $\displaystyle\text{V-IV: }\frac{T_{c2}^{1/2}}{T_{L2}}$ $\displaystyle=$
$\displaystyle\frac{\pi^{1/2}(-\kappa^{{}^{\prime\prime}}(\pi)/2)^{1/2}}{3\lambda_{3}v_{s}\rho(\pi)},$
(242)

respectively. A brief discussion about interaction-driven quantum cooling is
given in LPG:2022 . Based on previous results (96)- (99), the leading
contributions to minimum temperature near quadruple critical point are

$\displaystyle\text{II-IV: }\frac{T_{c1}^{1/2}}{T_{L1}}$
$\displaystyle\approx$
$\displaystyle\frac{\pi^{1/2}\eta_{1}^{1/2}}{6\lambda_{1}\lambda_{3}\delta},$
(243) $\displaystyle\text{V-IV: }\frac{T_{c2}^{1/2}}{T_{L2}}$
$\displaystyle\approx$
$\displaystyle\frac{2\lambda_{1}\pi^{5/2}}{3\lambda_{3}\eta_{1}A}.$ (244)

These suggest that the lowest temperature can be reached around QCPs, which
can be also reachable through adiabatic demagnetizaion cooling, see the study
of the Grüneisen parameters in yu2020gruneisen . Moreover, we further note
that the location of the minimum temperature for each contour entropy line in
the $T$-$u$ plane is governed by the relation $\alpha_{u}\Delta u/T\approx
1.3117$ in (235), where $\alpha_{u}$ needs to be determined.

###  V.4 V.4 Calculation of $\alpha_{u}$ for phase transition IV-V

Figure 17: (a) 3D plot of density $n$ in $u-\mu-n$ coordinate with fixed
$B=0.82714$: showing quantum phases I, II (or IV) and III (or V). The red
dash-dotted line in the $u-\mu$ plane represents the projection of 3D when the
density firstly reaches unity from phase II or IV. And the green circle on the
$n=1$ line denotes small areas around $u_{c}=1$. (b) Quantum scaling behaviour
of density near phase transition from IV to V driven by interaction. Numerical
result from the TBA equations and analytical result are in good agreement.

One remaining problem from the discussion in the previous subsection is that
the underlying coefficient $\alpha_{u}$ is not given. It not only determines
the scaling factor of the critical temperature at quantum criticality driven
by the interaction, but also determines the extreme point associated with the
lowest temperature in an isotropic process in the refrigeration cycle. We know
that $\alpha_{u}$ emerged in polylog function is related to the phase
transition line via the formula $\kappa(\pi)=\alpha_{B}\Delta
B+\alpha_{\mu}\Delta\mu+\alpha_{u}\Delta u$. For transition IV-V, the phase
boundary between phase IV and V is marked by $\kappa(\pi)=0$, where
$\kappa(\pi)$ is expressed via $\varepsilon_{1}(\Lambda)$ by

$\displaystyle\kappa(\pi)$ $\displaystyle=$ $\displaystyle
2-\mu-2u-B+\int_{-A}^{A}\mathrm{d}\Lambda
a_{1}(\Lambda)\varepsilon_{1}(\Lambda),$ (245)
$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$ $\displaystyle
2B-4\textmd{Re}\sqrt{1-(\Lambda-\mathrm{i}u)^{2}}+4u-\int_{-A}^{A}\mathrm{d}\Lambda^{\prime}a_{2}\left(\Lambda-\Lambda^{\prime}\right)\varepsilon_{1}(\Lambda^{\prime}).$
(246)

However for the transition IV-V, the corresponding equations (245) and (246)
are not in polynomial forms in terms of interaction $u$. In this part, we
present a numerical scheme with equation (223) based on the quality that phase
transition from IV to V occurs due to the emergence of Mott insulator with
constant density. In addition to the magnetic field $B$ and chemical potential
$\mu$ that can drive the phase transition, the interaction also drives the
system from one phase to others. In figure 17 (a), we plot the phase diagram
in the ($u-\mu-n_{c}$) coordinate at fixed magnetic field $B=0.82714$. The
vacuum phase I has $n_{c}=0$. The density is $0<n_{c}<1$ corresponding to
phase II or IV and $n_{c}=1$ for the phase III or V. The boundary line for IV
and V (or II and III) is marked with a constant density $n_{c}=1$. On this
transition line, we perform total derivatives of $n_{c}$, yielding

$\mathrm{d}n_{c}=\frac{\partial n_{c}}{\partial u}\mathrm{d}u+\frac{\partial
n_{c}}{\partial\mu}\mathrm{d}\mu=0.$ (247)

Substituting (223) for $\partial n_{c}/\partial u=-\partial C/\partial\mu$
with (227) and $\partial n_{c}/\partial\mu=\chi_{c}$ with (188) into (247),
near the boundary line $n_{c}=1$ in the $u-\mu$ plane, we have

$\alpha_{u}=-\alpha_{\mu}\frac{\partial\mu}{\partial u}.$ (248)

where we use $b_{1}=\alpha_{u},-1=\alpha_{\mu}$. Near the phase transition IV
to V, we see $\alpha_{\mu}=-1$ since the chemical potential $\mu$ only appears
in the charge leading term, no contribution to the driving terms in spin
string, see (245) and (246). At this point, we relate the value of
$\alpha_{u}$ to the slope at QCP in the $u_{c}-\mu_{c}$ plane. Therefore
equation (248) provides us with a way to obtain the unknown coefficient
$\alpha_{u}$. Moreover, $\partial\mu/\partial u$ can also be obtained directly
from (245) and (246)

$\displaystyle 2\frac{\partial\mu}{\partial u}$ $\displaystyle=$
$\displaystyle-2+\int_{-A}^{A}\mathrm{d}\Lambda\left[\frac{\partial
a_{1}(\Lambda)}{\partial
u}\varepsilon_{1}(\Lambda)+a_{1}(\Lambda)\frac{\partial\varepsilon_{1}(\Lambda)}{\partial
u}\right],$ (249)
$\displaystyle\frac{\partial\varepsilon_{1}(\Lambda)}{\partial u}$
$\displaystyle=$ $\displaystyle-4\frac{\partial}{\partial
u}[\textmd{Re}\sqrt{1-(\Lambda-\mathrm{i}u)^{2}}-u]-\int_{-A}^{A}\mathrm{d}\Lambda^{\prime}\frac{\partial}{\partial
u}[a_{2}\left(\Lambda-\Lambda^{\prime}\right)\varepsilon_{1}(\Lambda^{\prime})],$
(250)

where the factor $2$ emerged at the leftmost end of (249) arises from the
derivative of $\kappa(\pi)$ and $\mu$. Then $\partial\mu/\partial u$ can be
obtained by iteratively solving the above two equations. Here, we use color
map to get $\alpha_{u}$ instead, which can apply to arbitrary systems without
knowing explicit analytic formulas. In figure 17 (a), we plot density phase
diagram containing QCP $(u_{c},\mu_{c})$ of interest. In order to evaluate
$\alpha_{u}(u_{c},\mu_{c})$ numerically for the interaction-driven phase
transitions, compared to figure 12 (c) (d) driven by external fields, we
choose two adjacent points $(u_{1},\mu_{1}),(u_{2},\mu_{2})$ near QCP
$u_{c}=1$, see the green areas in figure 17 (a). It is evaluated from (248)
that
$\alpha_{u}(u_{c}=1)\approx(\mu_{1}-\mu_{2})/(u_{1}-u_{2})\approx-1.9627$.
Using this value as the argument factor of the scaling function for the phase
transition IV-V with the other part already known in equations (185) and
(200), we plot the density scaling law in figure 17 (b) in terms of the
variation of interaction strength. One can see that the obtained $\alpha_{u}$
captures well the thermodynamic scaling law.

###  V.5 V.5 Calculation of $\alpha_{u}$ for phase transition II-IV

Similar to the analysis of the parameter $\alpha_{u}$ for the transition of IV
to V, the relevant coefficients $\alpha_{B},\alpha_{\mu},\alpha_{u}$ for
transition II to IV are related to each other. Contrast to the Mott-insulator
transition, this transition arises from the introduction of spin-down
electrons, i.e., the phase transition points are determined by zero spin-down
particle density $n_{\downarrow}=0$. Considering a magnetic field $B$ and
coupling strength $u$ driven phase transition, we perform total derivative of
$n_{\downarrow}$, i.e.,

$\mathrm{d}n_{\downarrow}=\frac{\partial n_{\downarrow}}{\partial
u}\mathrm{d}u+\frac{\partial n_{\downarrow}}{\partial(2B)}\mathrm{d}(2B)=0.$
(251)

Substituting $\partial m/\partial B=\chi_{s}$ with (168) and (225) into the
(251), with the help of $n_{\downarrow}=n_{c}/2-m$ we obtain

$\alpha_{u}=-\alpha_{B}\frac{\partial B}{\partial u}.$ (252)

Recall that (205) gives us the value of
$\alpha_{B}=-2\left[1-\arctan\left(\sin Q_{c}/u\right)/\pi\right]$. We also
note that the quantity $\partial B/\partial u$ can be obtained from (203),
namely, near the critical line, we have the following relation

$B=\frac{2u}{\pi}\int_{0}^{Q}\mathrm{d}k\frac{\cos k}{u^{2}+\sin^{2}k}(\cos
k-\cos Q).$ (253)

Therefore, $\partial B/\partial u$ is easily obtained through the derivative
on both sides of this equation with respect to $u$. From equation (252), we
finally get the following analytic expression

$\displaystyle\alpha_{u}$ $\displaystyle=$
$\displaystyle\frac{4}{\pi}\left[\int_{0}^{Q_{c}}\mathrm{d}k\frac{\cos^{2}k}{u_{c}^{2}+\sin^{2}k}-2u_{c}^{2}\int_{0}^{Q_{c}}\mathrm{d}k\frac{\cos^{2}k}{(u_{c}^{2}+\sin^{2}k)^{2}}\right.$
(254) $\displaystyle\left.+\arctan\left(\frac{\sin
Q_{c}}{u_{c}}\right)+\frac{\sin(2Q_{c})}{2(u_{c}^{2}+\sin Q_{c}^{2})}\right].$

In summary, there exist simple relations between any pair of these scaling
factors $\alpha_{B},\,\alpha_{\mu}$ and $\alpha_{u}$. In the previous two
subsections, via Contact susceptibilities expressions (223) and (225), we
build up relationships between $\alpha_{u}$ and $\alpha_{\mu}$ as well as
between $\alpha_{u}$ and $\alpha_{B}$, see (248) and (252). In fact, in terms
of thermodynamic potential, Maxwell relations can also build a connection
between the charge susceptibility $\partial n_{c}/\partial(2B)$ and
magnetization susceptibility $\partial m/\partial\mu$. Moreover, (248) and
(252) remain valid for other transitions, like I-II, II-III and III-V. Based
on this feature, we can conclude that the following relations

$\displaystyle\frac{\alpha_{u}}{\alpha_{\mu}}$ $\displaystyle=$
$\displaystyle-\frac{\partial\mu}{\partial u},$ (255)
$\displaystyle\frac{\alpha_{u}}{\alpha_{B}}$ $\displaystyle=$
$\displaystyle-\frac{\partial B}{\partial u},$ (256)
$\displaystyle\frac{\alpha_{B}}{\alpha_{\mu}}$ $\displaystyle=$
$\displaystyle-\frac{\partial\mu}{\partial B}$ (257)

hold true in general, i.e., for arbitrary interaction strength and density. In
the above, we have discussed the exact values of $\alpha_{B}$ for the phase
transition II-IV and $\alpha_{\mu}$ for the phase transition IV-V. The
remaining coefficients can be derived through equations (255), (256), (257) in
a straightforward way. Here, we give the values of $\alpha_{\mu}$ for II-IV
and $\alpha_{B}$ for IV-V

$\displaystyle\text{II-IV: }\alpha_{\mu}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\arctan\left(\frac{\sin Q_{c}}{u}\right),$ (258)
$\displaystyle\text{IV-V: }\alpha_{B}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}+\int_{0}^{A}\mathrm{d}\Lambda
a_{1}(\Lambda)\frac{\partial\varepsilon_{1}(\Lambda)}{\partial B}.$ (259)

Together with equations (196)-(200), the results obtained here offer more
general description of quantum criticality in terms of full internal and
external potentials in the 1D Hubbard model in arbitrary experimentally
controllable parameters. These hold true for the second order phase transition
in higher dimensional quantum systems too.

##  VI VI. Conclusion and remarks

We have presented analytical results of thermal and magnetic properties of 1D
Hubbard model, ranging from elementary spin and change excitations, to the
spin incoherent Luttinger liquid, universal thermodynamics, quantum
criticality and interaction-driven refrigeration. A summary of our new results
is the following:
1) We have studied elementary excitations involving fractional spin and charge
excitations, gapped charge $k-\Lambda$ strings and spinon bound state
excitations as well as some combinations of them in section II. These are
complementary to that studied in essler2005one . Based on the study of such
excitations, together with the dimensionless Wilson ratio, universal scalings
of thermodynamics of the model, we have a rigorous investigation of the spin
incoherent Luttinger liquid, which was previously studied only in the
framework of effective theory via bosonization. From the conformal field
theory point of view, we have given rigorously the characteristics of various
correlation functions near the phase transition from phase II to phase IV,
showing an existence of collective mode in charge degrees of freedom rather
than the spin degrees of freedom.
2) We have presented general analytical results of thermodynamics, independent
of microscopic details of the model. Explicitly, we have determined the
additivity rules of charge and spin susceptibilities in grand canonical
ensemble equations (103) (104) and canonical ensemble equations (113) (114).
Away from the critical point and in the low-energy regime, in general,
coherent spin and charge degrees of freedom give rise to the well known
phenomenon of spin-charge separation, indicating the nature of the TLL. The
crossover regime between the TLL and the quantum critical region belongs to
the SILL, indicating a coexistence of liquid and gas, see equation (211).
Besides the thermodynamics, the quantum criticality and universal scaling laws
induced by the variation of magnetic field and chemical potential have been
obtained analytically and confirmed numerically, see general results of
criticality equations (196)-(200). In addition, the interaction-driven quantum
critical behaviour has been studied too.
3) We have introduced the lattice version of the Contact and Contact
susceptibilities with respect to the external fields and temperature. In
particular, we have investigated applications of Contact susceptibilities,
that build up a general connection between interaction-driven quantum
criticality and the phase transitions induced by external fields. By virtue of
the Contact susceptibilities, we have discussed Mott transition, quantum
refrigeration and interaction-induced quantum transitions, see section V.

In view of the rapid advances in trapping and controlling ultracold atoms in
experiment, the results obtained here will provide direct guidance to explore
experimentally various many-body phenomena in the 1D Hubbard model, such as
quantum criticality, spin coherent and incoherent Luttinger liquids,
generalized hydrodynamics and transport properties. These relations reveal
deep insights into quantum criticality of the Hubbard model in higher
dimensions. Furthermore, applications of our method to quantum metrology and
other quantum technologies are also highly desirable.

## Acknowledgements

We thank the Integrable Theory Group at APM for their helpful discussions.
X.W.G. acknowledges kind hospitalities of Rice University and Institute for
Advanced Study, Tsinghua University during his visits in 2018 and 2021,
respectively.

J.J.L. performed the analytical and numerical study of the model. X.W.G.
supervised J.J.L. for conducting the results reported in this paper. Both
X.W.G. and H.P. initiated this study.

J.J.L. and X.W.G. is supported by the NSFC key grant No. 12134015, the NSFC
grant No. 11874393 and No. 12121004. X.W.G. is also partially supported by the
Innovation Program for Quantum Science and Technology 2021ZD0302000. H.P.
acknowledges support from the US NSF (PHY-2207283) and the Welch Foundation
(Grant No. C-1669).

##  Appendix A APPENDIX A: Universal thermodynamics

Universal low energy physics can be derived from the TBA equations under the
conditions of $B/T\gg 1$ and $|\mu|/T\gg 1$. In such low energy limits, the
TBA equations (14) and (15) reduce to two coupled equations in terms of charge
quasimomentum $k$ and spin length-1 $\Lambda$ string, see (130) and (131).
Analysing the structure of the simplified TBA equations (130) and (131), we
find that it is efficient to take Sommerfeld expansion of the kernel functions
in terms of the Fermi points. Let us denote the fermi point $Q,\,A$ for charge
and spin degrees of freedom satisfying $\kappa(Q)=0,\varepsilon_{1}(A)=0$. For
our convenience, we denote $\frac{\partial{\bar{a}_{1}(\sin
k,\Lambda)}}{\partial{\Lambda}}=a_{1}(\sin k-\Lambda)$,
$\frac{\partial{\bar{a}(\sin k,\Lambda)}}{\partial{\Lambda}}=a_{1}(\sin
k-\Lambda)+a_{1}(\sin k+\Lambda)$. Now we focus on the second term in the
(130)

$\displaystyle-\int_{-\infty}^{\infty}\mathrm{d}\Lambda
a_{1}(\sin{k}-\Lambda)T\ln(1+\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}})$
(A1) $\displaystyle=$ $\displaystyle-\bar{a}_{1}(\sin
k,\Lambda)T\ln(1+\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}})|_{-\infty}^{\infty}-\int_{0}^{\infty}\mathrm{d}\Lambda\bar{a}(\sin{k},\Lambda)\frac{1}{1+\mathrm{e}^{\frac{\varepsilon_{1}(\Lambda)}{T}}}\frac{\partial\varepsilon_{1}(\Lambda)}{\partial{\Lambda}}$
$\displaystyle=$
$\displaystyle-\int_{0}^{\infty}\mathrm{d}\Lambda\bar{a}(\sin{k},\Lambda)\frac{1}{1+\mathrm{e}^{\frac{\varepsilon_{1}(\Lambda)}{T}}}\frac{\partial\varepsilon_{1}(\Lambda)}{\partial{\Lambda}}$
$\displaystyle=$
$\displaystyle-\int_{\varepsilon_{1}(0)}^{\varepsilon_{1}(\infty)}\mathrm{d}\varepsilon_{1}\bar{a}(\sin{k},\Lambda(\varepsilon_{1}))\frac{1}{1+\mathrm{e}^{\frac{\varepsilon_{1}}{T}}}$
$\displaystyle=$
$\displaystyle-\int_{\frac{\varepsilon_{1}(0)}{T}}^{\frac{\varepsilon_{1}(A)}{T}}\mathrm{d}xT\bar{a}(\sin{k},\Lambda(Tx))\frac{1}{1+\mathrm{e}^{x}}-\int_{\frac{\varepsilon_{1}(A)}{T}}^{\frac{\varepsilon_{1}(\infty)}{T}}\mathrm{d}xT\bar{a}(\sin{k},\Lambda(Tx))\frac{1}{1+\mathrm{e}^{x}}$
$\displaystyle=$
$\displaystyle-\int_{\frac{\varepsilon_{1}(0)}{T}}^{\frac{\varepsilon_{1}(A)}{T}}\mathrm{d}xT\bar{a}(\sin{k},\Lambda(Tx))\left(1-\frac{1}{1+\mathrm{e}^{-x}}\right)-\int_{\frac{\varepsilon_{1}(A)}{T}}^{\frac{\varepsilon_{1}(\infty)}{T}}\mathrm{d}xT\bar{a}(\sin{k},\Lambda(Tx))\frac{1}{1+\mathrm{e}^{x}}$
$\displaystyle=$
$\displaystyle-\int_{\frac{\varepsilon_{1}(0)}{T}}^{0}\mathrm{d}xT\bar{a}(\sin{k},\Lambda(Tx))+\int_{0}^{-\frac{\varepsilon_{1}(0)}{T}}\mathrm{d}xT\bar{a}(\sin{k},\Lambda(-Tx))\frac{1}{1+\mathrm{e}^{x}}$
$\displaystyle-\int_{0}^{\frac{\varepsilon_{1}(\infty)}{T}}\mathrm{d}xT\bar{a}(\sin{k},\Lambda(Tx))\frac{1}{1+\mathrm{e}^{x}}$
$\displaystyle=$
$\displaystyle-\int_{\varepsilon_{1}(0)}^{0}\mathrm{d}\varepsilon_{1}\bar{a}(\sin{k},\Lambda(\varepsilon_{1}))+\int_{0}^{\infty}\mathrm{d}xT\frac{\bar{a}(\sin{k},\Lambda(-Tx))-\bar{a}(\sin{k},\Lambda(Tx))}{1+\mathrm{e}^{x}}$
$\displaystyle=$
$\displaystyle\int_{\varepsilon_{1}(0)}^{0}\mathrm{d}\varepsilon_{1}\frac{\partial{\Lambda}}{\partial{\varepsilon_{1}}}\varepsilon_{1}(a_{1}(\sin{k}-\Lambda)+a_{1}(\sin{k}+\Lambda))$
$\displaystyle-\int_{0}^{\infty}\mathrm{d}xT\frac{(a_{1}(\sin{k}-\Lambda(0))+a_{1}(\sin{k}+\Lambda(0)))\Lambda^{{}^{\prime}}(0)2Tx}{1+\mathrm{e}^{x}}$
$\displaystyle=$
$\displaystyle\int_{0}^{A}\mathrm{d}\Lambda\varepsilon_{1}(\Lambda)(a_{1}(\sin{k}-\Lambda)+a_{1}(\sin{k}+\Lambda))-2T^{2}(a_{1}(\sin{k}-A)+a_{1}(\sin{k}+A))\Lambda^{{}^{\prime}}(0)\frac{\pi^{2}}{12}$
$\displaystyle=$
$\displaystyle\int_{0}^{A}\mathrm{d}\Lambda\varepsilon_{1}(\Lambda)(a_{1}(\sin{k}-\Lambda)+a_{1}(\sin{k}+\Lambda))$
$\displaystyle-\frac{\pi^{2}T^{2}}{6}\frac{1}{\varepsilon_{1}^{{}^{\prime}}(A)}(a_{1}(\sin{k}-A)+a_{1}(\sin{k}+A)).$

In the above calculations, the following equations are used:

$\displaystyle\frac{\bar{a}(\sin{k},\Lambda(Tx))-\bar{a}(\sin{k},\Lambda(-Tx))}{2Tx}=\frac{\partial{\bar{a}}}{\partial{\Lambda}}|_{\Lambda=\Lambda(0)}\Lambda^{{}^{\prime}}(0),$
$\displaystyle\Lambda^{{}^{\prime}}(0)=\frac{1}{\frac{d\varepsilon_{1}}{d\Lambda}|_{\varepsilon_{1}=0}}=\frac{1}{\varepsilon^{{}^{\prime}}_{1}(A)},\qquad\int_{0}^{\infty}\mathrm{d}x\frac{x}{1+\mathrm{e}^{x}}=\frac{\pi^{2}}{12}.$
(A2)

Conducting similar calculation in the integrals of $\varepsilon_{1}(\Lambda)$,
the dressed energies can be converted to

$\displaystyle\kappa(k)$ $\displaystyle=$
$\displaystyle-2\cos{k}-\mu-2u-B+\int_{-A}^{A}\mathrm{d}\Lambda\varepsilon_{1}(\Lambda)a_{1}(\sin{k}-\Lambda)$
(A3)
$\displaystyle-\frac{\pi^{2}T^{2}}{6}\frac{1}{\varepsilon_{1}^{{}^{\prime}}(A)}(a_{1}(\sin{k}-A)+a_{1}(\sin{k}+A)),$
$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$ $\displaystyle
2B-\frac{\pi^{2}T^{2}}{6\kappa^{{}^{\prime}}(Q)}\cos
Q(a_{1}(\sin{Q}-\Lambda)+a_{1}(\sin{Q}+\Lambda))$ (A4)
$\displaystyle+\frac{\pi^{2}T^{2}}{6\varepsilon_{1}^{{}^{\prime}}(A)}(a_{2}(\Lambda-A)+a_{2}(\Lambda+A))$
$\displaystyle+\int_{-Q}^{Q}\mathrm{d}k\kappa(k)\cos ka_{1}(\sin{k}-\Lambda)$
$\displaystyle-\int_{0}^{A}\mathrm{d}\Lambda^{{}^{\prime}}\varepsilon_{1}(\Lambda^{{}^{\prime}})(a_{2}(\Lambda-\Lambda^{{}^{\prime}})+a_{2}(\Lambda+\Lambda^{{}^{\prime}})).$

Compare above equations with root densities (20) and (21), we gain insight
that there is an implicit connection between the two sets of equations which
are highly symmetric. In order to get a compact and closed form of the free
energy, combing equations (20) and (21) with (A3) and (A4) and using
Sommerfeld expansion, we obtain the free energy

$\displaystyle f$ $\displaystyle=$
$\displaystyle-\int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi}T\ln(1+\mathrm{e}^{-\frac{\kappa(k)}{T}})+u$
(A5) $\displaystyle=$
$\displaystyle-\frac{\pi^{2}T^{2}}{6\kappa^{{}^{\prime}}(Q)}+\int_{0}^{Q}\frac{\mathrm{d}k}{\pi}\kappa(k)+u$
$\displaystyle=$
$\displaystyle-\frac{\pi^{2}T^{2}}{6\kappa^{{}^{\prime}}(Q)}+u$
$\displaystyle+\int_{-Q}^{Q}\mathrm{d}k\left[-2\cos{k}-\mu-2u-B-\frac{\pi^{2}T^{2}}{6}\frac{1}{\varepsilon_{1}^{{}^{\prime}}(A)}(a_{1}(\sin{k}-A)+a_{1}(\sin{k}+A))\right]\rho(k)$
$\displaystyle+\int_{-A}^{A}\mathrm{d}\Lambda\left[2B-\frac{\pi^{2}T^{2}}{6\kappa^{{}^{\prime}}(Q)}\cos
Q(a_{1}(\sin{Q}-\Lambda)+a_{1}(\sin{Q}+\Lambda))\right.$
$\displaystyle\left.+\frac{\pi^{2}T^{2}}{6\varepsilon_{1}^{{}^{\prime}}(A)}(a_{2}(\Lambda-A)+a_{2}(\Lambda+A))\right]\sigma_{1}(\Lambda)$
$\displaystyle=$
$\displaystyle-\frac{\pi^{2}T^{2}}{6\kappa^{{}^{\prime}}(Q)}+u+\left[\int_{-Q}^{Q}\mathrm{d}k\left(-2\cos{k}-\mu-2u-B\right)\rho(k)+2B\int_{-A}^{A}\mathrm{d}\Lambda\sigma_{1}(\Lambda)\right]$
$\displaystyle-\frac{\pi^{2}T^{2}}{6\varepsilon_{1}^{{}^{\prime}}(A)}\left[\int_{-Q}^{Q}\mathrm{d}k(a_{1}(\sin{k}-A)+a_{1}(\sin{k}+A))\rho(k)-(a_{2}(\Lambda-A)+a_{2}(\Lambda+A))\sigma_{1}(\Lambda)\right]$
$\displaystyle-\frac{\pi^{2}T^{2}}{6\kappa^{{}^{\prime}}(Q)}\left[\int_{-A}^{A}\mathrm{d}\Lambda\cos
Q(a_{1}(\sin{Q}-\Lambda)+a_{1}(\sin{Q}+\Lambda))\right]$ $\displaystyle=$
$\displaystyle-\frac{\pi^{2}T^{2}}{6\kappa^{{}^{\prime}}(Q)}+f_{0}-\frac{\pi^{2}T^{2}}{6\varepsilon_{1}^{{}^{\prime}}(A)}\left[\sigma_{1}(A)+\sigma_{1}(-A)+\int_{-A}^{A}\mathrm{d}\Lambda(a_{2}(\Lambda-A)+a_{2}(\Lambda+A))\sigma_{1}(\Lambda)\right]$
$\displaystyle-\frac{\pi^{2}T^{2}}{6\kappa^{{}^{\prime}}(Q)}\left[\rho(Q)+\rho(-Q)-\frac{1}{\pi}\right]+\frac{\pi^{2}T^{2}}{6\varepsilon_{1}^{{}^{\prime}}(A)}\int_{-A}^{A}\mathrm{d}\Lambda(a_{2}(\Lambda-A)+a_{2}(\Lambda+A))\sigma_{1}(\Lambda)$
$\displaystyle=$ $\displaystyle
f_{0}-\frac{\pi^{2}T^{2}\rho(Q)}{3\kappa^{{}^{\prime}}(Q)}-\frac{\pi^{2}T^{2}\sigma_{1}(A)}{3\varepsilon_{1}^{{}^{\prime}}(A)}$
$\displaystyle=$ $\displaystyle f_{0}-\frac{\pi
T^{2}}{6}\left(\frac{1}{v_{c}}+\frac{1}{v_{s}}\right),$

where
$f_{0}=\int_{-Q}^{Q}\mathrm{d}k\left(-2\cos{k}-\mu-2u-B\right)\rho(k)+2B\int_{-A}^{A}\mathrm{d}\Lambda\sigma_{1}(\Lambda)$
is the background contribution from zero temperature. The form of the free
energy serves as universal thermodynamics of the TLL in terms of spin and
charge degrees of freedom. Notice that in the above calculation process,
following integral identity was used

$\left\\{\begin{array}[]{l}f_{1}=f^{(0)}_{1}+K_{12}*f_{2}\\\
f_{2}=f^{(0)}_{2}+K_{21}*f_{1}+K_{22}*f_{2}\\\
g_{1}=g^{(0)}_{1}+K_{21}*g_{2}\\\
g_{2}=g^{(0)}_{2}+K_{12}*g_{1}+K_{22}*g_{2}\end{array}\rightarrow\int
f_{1}g^{(0)}_{1}+\int f_{2}g^{(0)}_{2}=\int g_{1}f^{(0)}_{1}+\int
g_{2}f^{(0)}_{2}.\right.$ (A6)

In spin polarized band II, spin degree vanishes. The free energy in low
temperature is simplified as $f=f_{0}-\frac{\pi T^{2}}{6}\frac{1}{v_{c}}$,
specific heat $\frac{C_{v}}{T}=\frac{\pi}{3}\frac{1}{v_{c}}$ accordingly.
While in phase V with charge half filled, the free energy is denoted as
$f=f_{0}-\frac{\pi T^{2}}{6}\frac{1}{v_{s}}$ and specific heat
$\frac{C_{v}}{T}=\frac{\pi}{3}\frac{1}{v_{s}}$.

##  Appendix B APPENDIX B: Free energy of phase III

In previous consideration, the effects from $k-\Lambda$ strings are ignored at
low-temperature limit $T\ll B,\,\mu$. Meanwhile, one can infer from equations
(15) and (16) that the leading orders of $\varepsilon_{n}^{\prime}$ and high
$\Lambda$ string with length $n$ contributing to free energy are
$\mathrm{e}^{\frac{2n\mu}{T}}$ and $\mathrm{e}^{-\frac{2nB}{T}}$,
respectively, rendering these gapped terms omitted reasonably. In order to
estimate those energy scales, in this appendix, we rigorously calculate the
contributions to the thermodynamics from the lowest gapped $k-\Lambda$ strings
$\varepsilon_{1}^{\prime}$ and from the real spin rapidity $\varepsilon_{1}$.
To make calculations less difficult, we choose area that $\varepsilon_{1}$ is
strictly positive and make the two driving terms
$\mathrm{e}^{\frac{2\mu}{T}},\,\mathrm{e}^{-\frac{2B}{T}}$ commensurable. Here
we treat the problem in the weak chemical potential region and high particle
density regime, i.e. close to the quadruple critical point. Consequently, we
are intend to discuss the case in phase III. On the basis of above
limitations, for the phase III, the TBA equations in terms of length-1
$\varepsilon_{1}$ and $\varepsilon_{1}^{\prime}$ are rewritten as

$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$ $\displaystyle
2B-4\operatorname{Re}\sqrt{1-(\Lambda-
iu)^{2}}+4u+\sum_{m=1}^{\infty}A_{1m}*T\ln\left(1+e^{-\frac{\varepsilon_{m}}{T}}\right)$
(B1)
$\displaystyle-2J_{1}a_{1}(\Lambda)+2J_{2}\left[\frac{\pi}{u}a^{2}_{1}(\Lambda)-\frac{4\pi^{2}}{u^{2}}\Lambda^{2}a^{3}_{1}(\Lambda)\right],$
$\displaystyle\varepsilon^{{}^{\prime}}_{1}(\Lambda)$ $\displaystyle=$
$\displaystyle-2\mu-2J_{1}a_{1}(\Lambda)+2J_{2}\left[\frac{\pi}{u}a^{2}_{1}(\Lambda)-\frac{4\pi^{2}}{u^{2}}\Lambda^{2}a^{3}_{1}(\Lambda)\right],$
(B2) $\displaystyle\kappa(k)$ $\displaystyle=$
$\displaystyle-2\cos{k}-C_{1}2\sin^{2}k+C_{2},$ (B3)

in which $C_{1},C_{2}$ are the integral functions related to the gapped
dressed energies $\varepsilon_{1}$ and $\varepsilon_{1}^{\prime}$

$\displaystyle C_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{d}\Lambda\left[\frac{\pi}{u}a^{2}_{1}(\Lambda)-\frac{4\pi^{2}}{u^{2}}\Lambda^{2}a^{3}_{1}(\Lambda)\right]\left[T\ln(1+\mathrm{e}^{-\frac{\varepsilon^{{}^{\prime}}_{1}(\Lambda)}{T}})-T\ln(1+\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}})\right],$
(B4) $\displaystyle C_{2}$ $\displaystyle=$
$\displaystyle-\mu-2u-B+\int_{0}^{\infty}\mathrm{d}\Lambda
2a_{1}(\Lambda)\left[T\ln(1+\mathrm{e}^{-\frac{\varepsilon^{{}^{\prime}}_{1}(\Lambda)}{T}})-T\ln(1+\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}})\right],$
(B5)

and $J_{1},J_{2}$ come from charge degree

$\displaystyle J_{1}$ $\displaystyle=$
$\displaystyle\frac{T^{\frac{3}{2}}}{\sqrt{1+2C_{1}}}\Gamma(\frac{3}{2})\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{2+C_{2}}{T}}\right)-\frac{T^{\frac{5}{2}}}{8(1+2C_{1})^{\frac{5}{2}}}\Gamma(\frac{5}{2})\mathrm{Li}_{\frac{5}{2}}\left(-\mathrm{e}^{\frac{2+C_{2}}{T}}\right),$
(B6) $\displaystyle J_{2}$ $\displaystyle=$
$\displaystyle\frac{T^{\frac{5}{2}}}{3(1+2C_{1})^{\frac{3}{2}}}\Gamma(\frac{5}{2})\mathrm{Li}_{\frac{5}{2}}\left(-\mathrm{e}^{\frac{2+C_{2}}{T}}\right).$
(B7)

Above expressions are similar to previous equations (135)-(144), except that
the contribution from $k-\Lambda$ string is taken into account.

In terms of the above formations of TBA equations within high density region,
the Gibbs free energy per site is given by

$\displaystyle f$ $\displaystyle=$ $\displaystyle-\mu-u-
B-T\int_{-\infty}^{\infty}\frac{\mathrm{d}\Lambda}{\pi}\operatorname{Re}\frac{1}{\sqrt{1-(\Lambda-\mathrm{i}u)^{2}}}\ln\left(1+\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}}\right)$
(B8)
$\displaystyle-T\ln\left(1+\mathrm{e}^{\frac{\kappa(0)}{T}}\right)+\frac{1}{\pi}\frac{T^{\frac{3}{2}}}{\sqrt{1+2C_{1}}}\Gamma(\frac{3}{2})\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{2+C_{2}}{T}}\right)$
$\displaystyle+\frac{1}{24\pi}\frac{T^{\frac{5}{2}}(1+8C_{1})}{(1+2C_{1})^{\frac{5}{2}}}\Gamma(\frac{5}{2})\mathrm{Li}_{\frac{5}{2}}\left(-\mathrm{e}^{\frac{2+C_{2}}{T}}\right).$

In order to get $C_{1,2}$ or $J_{1,2}$, we should solve the coupled equations
(B1)- (B3). We may ignore the convolution terms whose leading order is
$1+\mathrm{e}^{-\frac{2B}{T}}$. After carefully analysing the terms related to
the parameter $\Lambda$, the integral term
$\ln(1+\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}})$ is approximately
equal to $\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}}$, and with
expansions around $\Lambda=0$ we arrive at,

$\displaystyle\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}}$
$\displaystyle\approx$
$\displaystyle\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}\mathrm{e}^{-\frac{\eta_{1}\Lambda^{2}}{T}}\mathrm{e}^{\frac{g_{\Lambda}}{T}},$
(B9)

where we denote
$g_{\Lambda}=2J_{1}a_{1}(\Lambda)-\frac{2\pi}{u}J_{2}a^{2}_{1}(\Lambda)+\frac{8\pi^{2}}{u^{2}}J_{2}\Lambda^{2}a^{3}_{1}(\Lambda)$.
Conducting similar manipulation to charge bound states and using (B4) and (B5)
related with the integral $\varepsilon_{1}$ and $\varepsilon_{1}^{\prime}$, we
get the coefficients $C_{1},\,C_{2}$

$\displaystyle C_{1}$ $\displaystyle=$ $\displaystyle
T\mathrm{e}^{\frac{2\mu}{T}}\left[\frac{\pi}{u}t_{2}-\frac{4\pi^{2}}{u^{2}}t_{3}\right]-T\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}\left[\frac{\pi}{u}s_{2}-\frac{4\pi^{2}}{u^{2}}s_{3}\right],$
(B10) $\displaystyle C_{2}$ $\displaystyle=$
$\displaystyle-\mu-2u-B+2Tt_{1}\mathrm{e}^{\frac{2\mu}{T}}-2Ts_{1}\mathrm{e}^{-\frac{2B+\eta_{0}}{T}},$
(B11)

in which

$\displaystyle s_{1}(s_{2},s_{3})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{d}\Lambda
a_{1}(\Lambda)(a^{2}_{1}(\Lambda),\Lambda^{2}a^{3}_{1}(\Lambda))\mathrm{e}^{-\frac{\eta_{1}\Lambda^{2}}{T}}\mathrm{e}^{\frac{g_{\Lambda}}{T}},$
(B12) $\displaystyle t_{1}(t_{2},t_{3})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{d}\Lambda
a_{1}(\Lambda)(a^{2}_{1}(\Lambda),\Lambda^{2}a^{3}_{1}(\Lambda))\mathrm{e}^{\frac{g_{\Lambda}}{T}},$
(B13)

resulting from the subleading terms of
$\varepsilon_{1},\,\varepsilon_{1}^{\prime}$. To this point, we note that (B1)
and (B2) depend on the value of $J_{1},\,J_{2}$, while (B3) relies on
$C_{1},C_{2}$. Meanwhile, we observe that $C_{1},\,C_{2}$ are given by
$J_{1},\,J_{2}$, in turn $J_{1}$ and $J_{2}$ also contain $C_{1}$ and $C_{2}$.
Therefore, the remaining problem is to get the expressions for $C_{1},\,C_{2}$
and the entire story can be solved. By complicated and lengthy iterative
procedures, the $J_{1}$ and $J_{2}$ can be obtained explicitly

$\displaystyle J_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\pi^{\frac{1}{2}}T^{\frac{3}{2}}\tau_{1}+\frac{1}{2}\pi^{\frac{1}{2}}T^{\frac{5}{2}}f_{\frac{3}{2}}\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}\left(\frac{\pi}{u}s^{(0)}_{2}-\frac{4\pi^{2}}{u^{2}}s^{(0)}_{3}\right)-\frac{3}{32}\pi^{\frac{1}{2}}T^{\frac{5}{2}}\tau_{2},$
(B14) $\displaystyle J_{2}$ $\displaystyle=$
$\displaystyle\frac{\pi^{\frac{1}{2}}T^{\frac{5}{2}}}{4}\tau_{2},$ (B15)

where
$f_{n}=\mathrm{Li}_{n}\left(-\mathrm{e}^{\frac{2-\mu-2u-B}{T}}\right),\,s^{(0)}_{2}=\frac{1}{2\pi^{3/2}u}\left(a^{-\frac{1}{2}}-a^{-\frac{3}{2}}\right),\,s^{(0)}_{3}=\frac{\pi^{-5/2}a^{-3/2}}{4},a=\frac{\eta_{1}u^{2}}{T}$,
erfc is complementary error function, and $\tau_{1},\,\tau_{2}$ are defined by

$\displaystyle\tau_{1}$ $\displaystyle=$ $\displaystyle
f_{\frac{3}{2}}+\left(1+\frac{T^{\frac{1}{2}}f_{\frac{3}{2}}}{2\pi^{\frac{1}{2}}u}+\frac{3Tf^{2}_{\frac{3}{2}}}{16\pi
u^{2}}\right)\mathrm{e}^{\frac{2\mu}{T}}f_{\frac{1}{2}}-\left[\mathrm{e}^{a}\operatorname{erfc}(\sqrt{a})+\right.$
$\displaystyle\left.\left(\frac{T^{\frac{1}{2}}f_{\frac{3}{2}}}{\pi^{\frac{3}{2}}u}+\frac{Tf^{2}_{\frac{3}{2}}}{2\pi^{2}u^{2}}\right)\frac{\pi^{\frac{1}{2}}}{a^{\frac{1}{2}}}-\left(\frac{2T^{\frac{1}{2}}f_{\frac{3}{2}}}{\pi^{\frac{3}{2}}u}+\frac{3Tf^{2}_{\frac{3}{2}}}{2\pi^{2}u^{2}}\right)\frac{\pi^{\frac{1}{2}}}{2a^{\frac{3}{2}}}\right]\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}f_{\frac{1}{2}},$
$\displaystyle\tau_{2}$ $\displaystyle=$ $\displaystyle
f_{\frac{5}{2}}+\mathrm{e}^{\frac{2\mu}{T}}f_{\frac{3}{2}}-\mathrm{e}^{a}\operatorname{erfc}(\sqrt{a})\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}f_{\frac{3}{2}}.$
(B16)

Substituting these results obtained here into the expression of free energy
(B8), then the final solution (63) is finally obtained.

## References

  * (1) Pruschke T, Jarrell M and Freericks J K 1995 Adv. Phys. 44 187–210
  * (2) Dagotto E 2005 Science 309 257–62
  * (3) Yanase Y, Jujo T, Nomura T, Ikeda H, Hotta T and Yamada K 2003 Phys. Rep. 387 1–149
  * (4) Miyake K 2007 J. Phys.: Condens. Matter 19 125201
  * (5) Mott N F 1968 Rev. Mod. Phys. 40 677
  * (6) Schulz H J 1990 Phys. Rev. Lett. 64 2831
  * (7) Stafford C A and Millis A J 1993 Phys. Rev. B 48 1409
  * (8) Jördens R, Strohmaier N, Günter K, Moritz H and Esslinger T 2008 Nature 455 204–7
  * (9) Kokalj J and McKenzie R H 2013 Phys. Rev. Lett. 110 206402
  * (10) Ramirez A P 1997 J. Phys.: Condens. Matter 9 8171
  * (11) Hart R A, Duarte P M, Yang T-L, Liu X, Paiva T, Khatami E, Scalettar R T, Trivedi N, Huse D A and Hulet R G 2015 Nature 519 211–4
  * (12) Chiu C S, Ji G, Mazurenko A, Greif D and Greiner M 2018 Phys. Rev. Lett. 120 243201
  * (13) Giamarchi T 2017 Nature 545 414–5
  * (14) Boll M, Hilker T A, Salomon G, Omran A, Nespolo J, Pollet L, Bloch I and Gross C 2016 Science 353 1257–60
  * (15) Parsons M F, Mazurenko A, Chiu C S, Ji G, Greif D and Greiner M 2016 Science 353 1253–6
  * (16) Cheuk L W, Nichols M A, Lawrence K R, Okan M, Zhang H, Khatami E, Trivedi N, Paiva T, Rigol M and Zwierlein M W 2016 Science 353 1260–4
  * (17) Schneider U et al 2012 Nat. Phys. 8 213–18
  * (18) Salomon G, Koepsell J, Vijayan J, Hilker T A, Nespolo J, Pollet L, Bloch I and Gross C 2019 Nature 565 56–60
  * (19) Mitra D, Brown P T, Guardado-Sanchez E, Kondov S S, Devakul T, Huse D A, Schauß P and Bakr W S 2018 Nat. Phys. 14 173–7
  * (20) Lieb E H and Wu F Y 1968 Phys. Rev. Lett. 20 1445
  * (21) Takahashi M 1972 Prog. Theor. Phys. 47 69–82
  * (22) Essler F H L, Frahm H, Göhmann F, Klümper A and Korepin V E 2005 The One-Dimensional Hubbard Model (Cambridge: Cambridge University Press)
  * (23) Koepsell J, Vijayan J, Sompet P, Grusdt F, Hilker T A, Demler E, Salomon G, Bloch I and Gross C 2019 Nature 572 358–62
  * (24) Scherg S, Kohlert T, Sala P, Pollmann F, Madhusudhana B H, Bloch I and Aidelsburger M 2021 Nat. Commun. 12 4490
  * (25) Gall M, Wurz N, Samland J, Chan C F and Köhl M 2021 Nature 589 40–3
  * (26) Vijayan J, Sompet P, Salomon G, Koepsell J, Hirthe S, Bohrdt A, Grusdt F, Bloch I and Gross C 2020 Science 367 186–9
  * (27) Hirsch J E 1985 Phys. Rev. B 31 4403
  * (28) Jarrell M 1992 Phys. Rev. Lett. 69 168
  * (29) Georges M and Kotliar G 1992 Phys. Rev. B 45 6479
  * (30) Glocke S, Klümper A and Sirker J 2007 Phys. Rev. B 76 155121
  * (31) De Leo L, Bernier J-S, Kollath C, Georges A and Scarola V W 2011 Phys. Rev. A 83 023606
  * (32) Ibarra-García-Padilla E, Mukherjee R, Hulet R G, Hazzard K R A, Paiva T and Scalettar R T 2020 Phys. Rev. A 102 033340
  * (33) Wietek A, Rossi R, Šimkovic IV F, Klett M, Hansmann P, Ferrero M, Stoudenmire E M, Schäfer T and Georges A 2021 Phys. Rev. X 11 041013
  * (34) Frahm H and Korepin V E 1990 Phys. Rev. B 42 10553
  * (35) Frahm H and Korepin V E 1991 Phys. Rev. B 43 5653
  * (36) Qin M, Schäfer T, Andergassen S, Corboz P and Gull E 2022 Annu. Rev. Condens. Matter Phys. 13 275–302
  * (37) Kim A J, Šimkovic IV F and Kozik E 2020 Phys. Rev. Lett. 124 117602
  * (38) Ilievski E and De Nardis J 2017 Phys. Rev. B 96 081118(R)
  * (39) Ilievski E, De Nardis J, Medenjak M and Prosen T 2018 Phys. Rev. Lett. 121 230602
  * (40) Guardado-Sanchez E, Morningstar A, Spar B M, Brown P T, Huse D A, and Bakr W S 2020 Phys. Rev. X 10 011042
  * (41) Bertini B, Heidrich-Meisner F, Karrasch C, Prosen T, Steinigeweg R and Žnidarič M 2021 Rev. Mod. Phys. 93 025003
  * (42) Moeckel M and Kehrein S 2008 Phys. Rev. Lett. 100 175702
  * (43) Eckstein M, Kollar M and Werner P 2009 Phys. Rev. Lett. 103 056403
  * (44) Schlünzen N, Joost J-P, Heidrich-Meisner F and Bonitz M 2017 Phys. Rev. B 95 165139
  * (45) Landau L 1959 Sov. Phys. JETP 8 70
  * (46) Baym G and Pethick C 2008 Landau Fermi-Liquid Theory: Concepts and Applications (New York: Wiley)
  * (47) Tomonaga S-i 1950 Prog. Theor. Phys. 5 544–69
  * (48) Luttinger J M 1963 J. Math. Phys. 4 1154–62
  * (49) Haldane F D M 1981 J. Phys. C: Solid State Phys. 14 2585
  * (50) Scopa S, Calabrese P and Piroli L 2022 Phys. Rev. B 106 134314
  * (51) Lorenz T, Hofmann M, Grüninger M, Freimuth A, Uhrig G S, Dumm M and Dressel M 2002 Nature 418 614–7
  * (52) Demler E, Nayak C, Kee H-Y, Kim Y B and Senthil T 2002 Phys. Rev. B 65 155103
  * (53) Kollath C, Schollwöck U and Zwerger W 2005 Phys. Rev. Lett. 95 176401
  * (54) Jompol Y, Ford C J B, Griffiths J P, Farrer I, Jones G A C, Anderson D, Ritchie D A, Silk T W and Schofield A J 2009 Science 325 597–601
  * (55) Schmidt T L, Imambekov A and Glazman L I 2010 Phys. Rev. B 82 245104
  * (56) Ma Y et al 2017 Nat. Commun. 8 14231
  * (57) He F, Jiang Y-Z, Lin H-Q, Hulet R G, Pu H and Guan X-W 2020 Phys. Rev. Lett. 125 190401
  * (58) Senaratne R, Cavazos-Cavazos D, Wang S, He F, Chang Y-T, Kafle A, Pu H, Guan X-W and Hulet R G 2022 Science 376 1305–8
  * (59) Recati A, Fedichev P O, Zwerger W and Zoller P 2003 Phys. Rev. Lett. 90 020401
  * (60) Voit J 1993 J. Phys.: Condens. Matter 5 8305
  * (61) Kim C, Matsuura A Y, Shen Z-X, Motoyama N, Eisaki H, Uchida S, Tohyama T and Maekawa S 1996 Phys. Rev. Lett. 77 4054
  * (62) Auslaender O M, Steinberg H, Yacoby A, Tserkovnyak Y, Halperin B I, Baldwin K W, Pfeiffer L N and West K W 2005 Science 308 88–92
  * (63) Kim B Jet al 2006 Nat. Phys. 2 397–401
  * (64) Veness T and Essler F H L 2016 Phys. Rev. B 93 205101
  * (65) Cheianov V V and Zvonarev M B 2004 Phys. Rev. Lett. 92 176401
  * (66) Fiete G A and Balents L 2004 Phys. Rev. Lett. 93 226401
  * (67) Fiete G A 2007 Rev. Mod. Phys. 79 801
  * (68) Fiete G A, Qian J, Tserkovnyak Y and Halperin B I 2005 Phys. Rev. B 72 045315
  * (69) Fiete G A, Le Hur K and Balents L 2005 Phys. Rev. B 72 125416
  * (70) Hew W K, Thomas K J, Pepper M, Farrer I, Anderson D, Jones G A C and Ritchie D A 2008 Phys. Rev. Lett. 101 036801
  * (71) Meden V 1999 Phys. Rev. B 60 4571
  * (72) Giamarchi T 2003 Quantum Physics in One Dimension (Oxford: Oxford University Press)
  * (73) Miranda E 2003 Braz. J. Phys. 33 3–35
  * (74) Cavazos-Cavazos D, Senaratne R, Kafle A and Hulet R G 2022 arXiv:2210.06306
  * (75) Giamarchi T 1991 Phys. Rev. B 44 2905
  * (76) Hubbard J 1963 Proc. R. Soc. Lond. A 276 238–57
  * (77) Brinkman W F and Rice T M 1970 Phys. Rev. B 2 4302
  * (78) Han X J, Liu Y, Liu Z Y, Li X, Chen J, Liao H J, Xie Z Y, Normand B and Xiang T 2016 New J. Phys. 18 103004
  * (79) Luther A and Emery V J 1974 Phys. Rev. Lett. 33 589
  * (80) Voit J 1998 Eur. Phys. J. B 5 505–19
  * (81) Orgad D 2001 Philos. Mag. B 81, 377–98
  * (82) Sciolla B, Tokuno A, Uchino S, Barmettler P, Giamarchi T and Kollath C 2013 Phys. Rev. A 88 063629
  * (83) Paiva T, Scalettar R, Randeria M and Trivedi N 2010 Phys. Rev. Lett. 104 066406
  * (84) Gorelik E V, Rost D, Paiva T, Scalettar R, Klümper A and Blümer N 2012 Phys. Rev. A 85 061602(R)
  * (85) Werner F, Parcollet O, Georges A and Hassan S R 2005 Phys. Rev. Lett. 95 056401
  * (86) Daré A-M, Raymond L, Albinet G and Tremblay A-M S 2007 Phys. Rev. B 76 064402
  * (87) Deguchi T, Essler F H L, Göhmann F, Klümper A, Korepin V E and Kusakabe K 2000 Phys. Rep. 331 197–281
  * (88) Yang C N and Yang C P 1969 J. Math. Phys. 10 1115
  * (89) Klümper A and Bariev R Z 1996 Nucl. Phys. B 458 623–39
  * (90) Batchelor M T, Guan X-W, Oelkers N and Tsuboi Z 2007 Adv. Phys. 56 465–543
  * (91) Carmelo J, Horsch P, Bares P A and Ovchinnikov A A 1991 Phys. Rev. B 44 9967
  * (92) Essler F H L and Korepin V E 1994 Phys. Rev. Lett. 72 908
  * (93) Feiguin A E and Heidrich-Meisner F 2007 Phys. Rev. B 76, 220508(R)
Lüscher A, Noack R M and Läuchli A M 2008 Phys. Rev. A 78 013637
Rizzi M, Polini M, Cazalilla M A, Bakhtiari M R, Tosi M P and Fazio R 2008
Phys. Rev. B 77 245105
Tezuka M and Ueda M 2008 Phys. Rev. Lett. 100 110403
Tezuka M and Ueda M 2010 New J. Phys. 12 055029

  * (94) Penc K, Mila F and Shiba H 1995 Phys. Rev. Lett. 75 894
  * (95) Krivnov V Ya and Ovchinnikov A A 1975 Zh. Eksp. Teor. Fiz. 67 1568–81
  * (96) Woynarovich F 1983 J. Phys. C: Solid State Phys. 16 6593
  * (97) Bogoliubov N M and Korepin V E 1988 Mod. Phys. Lett. B 1 349
Bogoliubov N M and Korepin V E 1990 Theor. Math. Phys. 82 231-43

  * (98) Woynarovich F and Penc K 1991 Z. Physik B - Condensed Matter 85 269-80
  * (99) Essler F H L and Korepin V E 1994 Nucl. Phys. B 426 505–33
  * (100) Sacramento P D 1995 J. Phys.: Condens. Matter 7 143
  * (101) Cheng S, Yu Y-C, Batchelor M T and Guan X-W 2018 Phys. Rev. B 97 121111(R)
  * (102) Cheng S, Jiang Y-Z, Yu Y-C, Batchelor M T and Guan X-W 2018 Nucl. Phys. B 929 353–76
  * (103) Cheng S, Yu Y-C, Batchelor M T and Guan X-W 2018 Phys. Rev. B 97 125145
  * (104) White S R, Scalapino D J, Sugar R L, Loh E Y, Gubernatis J E and Scalettar R T 1989 Phys. Rev. B 40 506
  * (105) Caffarel M and Krauth W 1994 Phys. Rev. Lett. 72 1545
  * (106) Wang Y, Esterlis I, Shi T, Cirac J I and Demler E 2020 Phys. Rev. Research 2 043258
  * (107) Timrov I, Marzari N and Cococcioni M 2021 Phys. Rev. B 103 045141
  * (108) Sénéchal D, Perez D and Pioro-Ladrière M 2000 Phys. Rev. Lett. 84 522
  * (109) Ovchinnikov A A 1969 Sov. Phys. JETP 29
  * (110) Nocera A, Patel N D, Fernandez-Baca J, Dagotto E and Alvarez G 2016 Phys. Rev. B 94 205145
  * (111) Mölter J, Barthel T, Schollwöck U and Alba V 2014 J. Stat. Mech. 2014 P10029
  * (112) Wang Z et al 2018 Nature 554 219–23
  * (113) Kawakami N, Usuki T and Okiji A 1989 Phys. Lett. A 137 287–90
  * (114) Usuki T, Kawakami N and Okiji A 1990 J. Phys. Soc. Jpn. 59 1357–65
  * (115) Takahashi M 1974 Prog. Theor. Phys. 52 103–14
  * (116) Wilson K G 1975 Rev. Mod. Phys. 47 773
  * (117) Guan X-W, Yin X-G, Foerster A, Batchelor M T, Lee C-H and Lin H-Q 2013 Phys. Rev. Lett. 111 130401
  * (118) Yu Y-C, Chen Y-Y, Lin H-Q, Römer R A and Guan X-W 2016 Phys. Rev. B 94 195129
  * (119) Boehler R and Ramakrishnan J 1980 J. Geophys. Res.: Solid Earth 85 6996–7002
  * (120) Küchler R et al 2003 Phys. Rev. Lett. 91 066405
  * (121) Grüneisen E 1908 Ann. Phys. 331 211–216
  * (122) Yu Y-C, Zhang S and Guan X-W 2020 Phys. Rev. Research 2 043066
  * (123) Guan X-W, Batchelor M T and Lee C 2013 Rev. Mod. Phys. 85 1633
  * (124) Sachdev S 2011 Quantum Phase Transitions (Cambridge: Cambridge University Press)
  * (125) Hazzard K R A and Mueller E J 2011 Phys. Rev. A 84 013604
  * (126) Fukui T and Kawakami N 1995 Phys. Rev. B 51 5239
  * (127) Vitoriano C and Coutinho-Filho M D 2009 Phys. Rev. Lett. 102 146404
  * (128) Zhang X B, Chen Y Y, Liu L X, Deng Y J and Guan X-W 2022 Natl. Sci. Rev. 9
  * (129) Continentino M A 1994 Phys. Rep. 239 179–213
  * (130) Continentino M A, Japiassu G M and Troper A 1989 Phys. Rev. B 39 9734
  * (131) Continentino M A 2005 Braz. J. Phys. 35 197
  * (132) Luo J-J, Pu H and Guan X-W 2023 Phys. Rev. B 107 L201103
  * (133) Mohankumar N 2007 Comput. Phys. Commun. 176 665–9
  * (134) Campo Jr V L, Capelle K, Hooley C, Quintanilla J and Scarola V W 2012 Phys. Rev. A 85 033644
  * (135) van Loon E G C P, Rösner M, Schönhoff G, Katsnelson M I and Wehling T O 2018 npj Quant. Mater. 3 32
  * (136) Tan S 2008 Ann. Phys. 323 2952–70
Tan S 2008 Ann. Phys. 323 2987–90
Tan S 2008 Ann. Phys. 323 2971–86

  * (137) Zhang S and Leggett A J 2009 Phys. Rev. A 79 023601
  * (138) Braaten E and Platter L 2009 Laser Phys. 19 550–3
  * (139) Stewart J T, Gaebler J P, Drake T E and Jin D S 2010 Phys. Rev. Lett. 104 235301
  * (140) Yu Z, Thywissen J H and Zhang S 2015 Phys. Rev. Lett. 115 135304
  * (141) Yoshida S M and Ueda M 2015 Phys. Rev. Lett. 115 135303
  * (142) Cui X 2016 Phys. Rev. A 94 043636
  * (143) Chen Y-Y, Jiang Y-Z, Guan X-W and Zhou Q 2014 Nat. Commun. 5 5140
  * (144) Giamarchi T 1997 Physica B: Condens. Matter 230 975–80