# Exact results of one-dimensional repulsive Hubbard model

Jia-Jia Luo  Innovation Academy for Precision Measurement Science and
Technology, Chinese Academy of Sciences, Wuhan 430071, China  University of
Chinese Academy of Sciences, Beijing 100049, China  Han Pu  Department of
Physics and Astronomy, and Rice Center for Quantum Materials, Rice University,
Houston, Texas 77251-1892, USA  Xi-Wen Guan<EMAIL_ADDRESS><EMAIL_ADDRESS>Innovation Academy for Precision Measurement Science and
Technology, Chinese Academy of Sciences, Wuhan 430071, China  Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China
Peng Huanwu Center for Fundamental Theory, Xi’an 710069, China  Department of
Fundamental and Theoretical Physics, Research School of Physics, Australian
National University, Canberra ACT 0200, Australia

###### Abstract

We present analytical results of fundamental properties of one-dimensional
(1D) Hubbard model with a repulsive interaction, ranging from fractional
excitations to universal thermodynamics, interaction-driven criticality,
correlation functions, Contact susceptibilities and quantum cooling. Using the
exact solutions of the Bethe Ansatz equations of the Hubbard model, we first
rigorously calculate the gapless spin and charge excitations, exhibiting
exotic features of fractionalized spinons and holons. While we investigate the
gapped excitations in terms of the spin string and the $k-\Lambda$ string
bound states at arbitrary driving fields, showing subtle differences of spin
magnons and charge $\eta$-pair excitations. Based on the analysis on the
fractional charge and spin excitations, the spin-incoherent Luttinger liquid
with only the charge propagation mode is elucidated by the asymptotic of the
two-point correlation functions with the help of the conformal field theory.
For a high density and high spin magnetization region, i.e. near quadruple
critical point, we then further analytically obtain the thermodynamical
properties, dimensionless ratios and scaling functions near quantum phase
transitions. These results of thermodynamics provide rich universal low energy
physics and critical scaling laws for the phase transitions near the quadruple
point. In particular, we determine additivity rules of spin and charge
susceptibilities, and thermodynamics of spin-incoherent Luttinger liquid.
Finally, in order to capture deeper insight into the Mott insulator and
interaction driven criticality, we further study the double occupancy and its
associated Contact and Contact susceptibilities through which an adiabatic
cooling scheme upon the quantum criticality is proposed Our methods provide
rich perspectives of quantum integrability in 1D strongly correlated systems
and offer promising guidance to future experiments with interacting electrons
and ultracold atoms with and without a lattice.

Keywords: repulsive Hubbard model, fractional excitations, thermodynamics,
spin-incoherent Luttinger liquid, quantum cooling, interaction driven
criticality

##  I I. Introduction

The strongly correlated electronic systems provide us with a powerful platform
to investigate novel and intriguing many-body phenomena, such as high-$T_{c}$
superconductor pruschke1995anomalous ; dagotto2005complexity ;
yanase2003theory ; miyake2007new , Mott metal-insualtor
transitionmott1968metal ; schulz1990correlation ; stafford1993scaling ;
jordens2008mott ; kokalj2013thermodynamics , colossal
magnetoresistancedagotto2005complexity ; ramirez1997colossal ,
antiferromagnetic correlations hart2015observation ; chiu2018quantum ;
giamarchi2017firmer ; boll2016spin ; parsons2016site ; schneider2012fermionic
; cheuk2016observation ; salomon2019direct ; mitra2018quantum  etc., which
cannot be merely explained by compounding of single particle motions. Among
the continuous and lattice models of interacting electrons and ultracold
atoms, one of the prototypical integrable models is the one-band Hubbard model
Lieb:1968 ; takahashi1972one ; essler2005one . Recently, particular attention
has been paid to this model because of its operability and accessibility in
trapping ultracold atoms on lattices jordens2008mott ; koepsell2019imaging ;
vijayan2020time ; scherg2021observing ; gall2021competing . There have been
numerous publications focusing on the 1D Hubbard model with on-site
interaction as well as its variants in one or higher dimensions hirsch1985two
; jarrell1992hubbard ; georges1992hubbard ; glocke2007half ;
de2011thermodynamics ; ibarra2020thermodynamics ; wietek2021mott , involving
the studies of spectral weights and critical exponents for correlation
functionsschulz1990correlation ; frahm1990critical ; frahm1991correlation ;
stafford1993scaling ; kim2020spin ; qin2022hubbard , dynamical
transportqin2022hubbard ; ilievski2017ballistic ; ilievski2018superdiffusion ;
guardado2020subdiffusion ; bertini2021finite , interaction
quenchmoeckel2008interaction ; eckstein2009thermalization ;
schlunzen2017nonequilibrium  etc. In this scenario, the well-celebrated Landau
theory of Fermi liquid landau1959theory ; baym2008landau  remarkably describes
universal low energy physics of interacting fermions and leads to a wide range
of applications in higher dimensional quantum many-body systems.

On the other hand, the low energy excitations of 1D many-body systems usually
lack individual motions of quasiparticles. Instead, they form collective
motions of bosons, i.e., Tomonaga-Luttinger liquid (TLL) tomonaga1950remarks ;
luttinger1963exactly ; haldane1981luttinger . A hallmark of interacting
electrons in 1D is the splitting of low-lying spin and charge excitations into
two separate collective motions lorenz2002evidence ;
demler2002fractionalization ; kollath2005spin ; jompol2009probing ;
schmidt2010spin ; ma2017angle ; he2020emergence ; senaratne2022spin ;
vijayan2020time ; Scopa2022 , i.e., one solely carrying spin and the other
carrying charge. This is one of the most significant features predicted by
both the integrable theory and the TLL theory, namely, the propagators
generated by excitations separate into two massless wave packets with
distinctive velocities $v_{c},\,v_{s}$ in charge and spin degrees of freedom,
respectively voit1993charge ; kollath2005spin ; Recati:PhysRevLett.90.020401 ;
lorenz2002evidence . The spin-charge separation phenomenon was experimentally
evidenced in literature Kim:1996 ; senaratne2022spin ; Kim:2006 ;
auslaender2005spin ; jompol2009probing . This phenomenon was further evidenced
by detecting spin and charge propagating velocities jompol2009probing ;
vijayan2020time  from quenching an initial Mott state into few hole doped
states in the 1D Hubbard model. While some high-energy or gapped excitations
may contribute to the optical conductivity veness2016mobile . The spin-charge
separation has been recently verified and confirmed by the virtue of spin and
charge dynamical structure factors in 1D Fermi gas of ultracold atoms
senaratne2022spin . In this research, the team used Bragg beams to excite spin
and charge density waves separately and measured the corresponding spin and
charge Bragg spectra at different interacting strengths.

By taking into account the curvature correction to the linear charge
excitation spectrum and the nonlinear effect of the spin backward scattering
in the spin sector, theoretical and experimental results are in good
agreement. Nevertheless, when the temperature $T$ exceeds a characteristic
energy scale of spin, but still less than typical energy of charge Fermi
energy, i.e. $E_{s}\ll T\ll E_{c}$, the spin degree of freedom is fully
disordered, while the charge degree of freedom still remains in propagating
mode. Such a state is referred to as spin incoherent Luttinger liquid
(SILL)cheianov2004nonunitary ; fiete2004green ; fiete2007colloquium ;
fiete2005theory ; fiete2005transport ; hew2008spin . The correlation function
of the SILL shows an exponential decay of distance in the spin sector, while
it remains a power-law decay in the charge sector, giving a characteristic of
the emergent SILL cheianov2004nonunitary ; fiete2004green ;
fiete2007colloquium ; fiete2005theory , in contrast to the case of the usual
TLL that correlation functions in both spin and charge sectors exhibit power-
law decaying in distance schulz1990correlation ; frahm1990critical ;
meden1999nonuniversality ; giamarchi2003quantum ; miranda2003introduction . An
evidence for the spin-incoherent Luttinger liquid was recently reported by R
Hulet’s group at Rice University Cavazos-Cavazos:2022 . Although there have
been extensive study on this model, a comprehensive study of the spin-coherent
and incoherent Luttinger liquids, universal thermodynamics, fractional
quasiparticles and Mott phase in the 1D Hubbard model still remains elusive.

Moreover, comparing with the continuum gas systems, one has a consensus that
the Mott insulator phase dominated by antiferromagnetic ordering is of great
importance on account of its uniqueness stemming from the umklapp scattering
in terms of bosonization giamarchi1991umklapp . There are ample thoughts and
debates on the mechanism for the onset of the Mott insulator for a long time
stafford1993scaling ; hubbard1963electron ; brinkman1970application ;
han2016charge . But it is still a big theoretical challenge to capture the
essential nature of the insulating phase. This massive charge sector with
massless spin degree belongs to the category of Luther-Emery model
luther1974backward ; voit1998dynamical ; orgad2001spectral  whose dynamical
spectral function for Mott phase displays one or two singularities without
anomalous dimension voit1998dynamical ; orgad2001spectral . A quantity which
is experimentally accessible is the double occupancy, embodying charge
fluctuations. It gives insights into the emergence of the Mott
phaseglocke2007half ; de2011thermodynamics ; sciolla2013competition  and
conveys plentiful information apart from Mott phenomenon, for instance, the
Pomeranchuk effect glocke2007half ; de2011thermodynamics ; wietek2021mott ,
antiferromagnetic correlations wietek2021mott ; paiva2010fermions ;
gorelik2012universal , adiabatic cooling werner2005interaction ;
dare2007interaction ; de2011thermodynamics  over a wide energy scales and
temperature regimes. One of the key contents in this paper is to define the
Contact and Contact susceptibilities through which we further study adiabatic
cooling scheme at quantum criticality.

The Bethe Ansatz provides us with a feasible way to solve the 1D Hubbard
Hamiltonian Lieb:1968 ; Degu2000 . The framework of building up the finite-
temperature thermodynamics of 1D integrable systems was attributed to C N Yang
and C P Yang’s seminal paper YangCN:1969 . Later M Takahashi made an important
step forward by treating the thermodynamics of some more complicated
integrable systems takahashi1972one . In this scenario, quantum transfer
matrix was successfully developed to the study of the thermodynamics of
integrable models in terms of path integral klumper1996exact ; Batchelor:2007
. These two approaches endow us with possibility to acquire the excitation
spectra carmelo1991renormalized ; Ess1994a , universal free Luttinger liquids
information he2020emergence  and correlation functions frahm1990critical ;
qin2022hubbard  for the repulsive Hubbard model and the attractive Hubbard
model DMRG ; Penc1995 ; Krivnov1975 ; Woyna1983 ; Bogoliubov1988 ; Woyna1991 ;
Ess1994a ; Ess1994b ; Sacra1995 ; Cheng:2018A ; Cheng:2018B ; Cheng:2018C
despite big challenges in actual calculations of physical properties due to
the complexity of the TBA equations at arbitrary temperature. For higher
dimensionality and various lattice symmetries hirsch1985two ;
georges1992hubbard ; jarrell1992hubbard ; ibarra2020thermodynamics , there
does not seem to have completely applicable theoretical schemes and commonly
exploited numerical simulation methods including quantum Monte
Carlojarrell1992hubbard ; dare2007interaction ; kim2020spin ; wietek2021mott ,
matrix-decomposition algorithms white1989numerical , dynamical mean-field
theory hirsch1985two ; ibarra2020thermodynamics ; wietek2021mott , exact
diagonalization caffarel1994exact ; wang2020zero , perturbation theory
timrov2021self ; senechal2000spectral , density-Matrix renormalization
groupkollath2005spin ; schlunzen2017nonequilibrium  and other generalization
georges1992hubbard ; senechal2000spectral .

The elementary excitations in spin and charge degree and universal
thermodynamics of the 1D repulsive Hubbard model are briefly discussed in
LPG:2022 . In this paper, we develop analytical methods to obtain fundamental
properties of the 1D Hubbard model with repulsive interaction and external
fields. In particular, using the exact solutions of the Bethe Ansatz equations
of the Hubbard model, we rigorously calculate fractional excitation spectra,
universal thermodynamics, scaling functions, incoherent correlation functions,
Contact, Contact susceptibilities, quantum cooling and interaction-driven
quantum phase transitions in a pedagodical way. Exact results of various
excitation spectra reveal exotic features of fractionalized spinons and holons
and show subtle differences of spin magnons and charge $\eta$-pair
excitations. Universal thermodynamical properties of spin and charge Luttigner
liquids, dimensionless ratios and scaling functions near quantum phase
transitions in terms of chemical potential, magnetic field and interaction are
obtained explicitly near quadruple critical point. The asymptotic two-point
correlation functions for the spin-incoherent Luttinger liquid are derived
with the help of the conformal field theory. Finally, we study the double
occupancy, the Contact, Contact susceptibilities, interaction-driven adiabatic
cooling scheme and criticality and Mott insulator behaviour of the model as
well.

The outline of this paper is the following. In section II, we introduce the
basic knowledge of exact Bethe ansatz solution of the repulsive Hubbard model
and calculate the excitation spectra of fractional holons and spinons as well
as the gapped excitations above the ground state. In section III, we initiate
our study on the thermal and magnetic properties of TLL, quantum criticality,
non-Fermi liquid behaviour and the universal laws in the 1D Hubbard model. In
section IV, we study the SILL from the Bethe ansatz perspective. The
asymptotic single particle Green’s function and pairing correlation function
are presented. A depiction of double occupancy and analytical calculation of
Contact and its susceptibilities with respect to the chemical potential,
magnetic field and interaction near Mott phase are given in section V. In this
section, we also carry out a brief analysis on quantum cooling and interaction
driven phase transitions. The last section VI remains for a summary and
outlooks.

##  II II. Fractional holons and spinons

Fundamental low energy physics is intimately related to the elementary
excitations in 1D many-body systems. In particular, the low-lying excited
states involve many-body correlations and determine the unique dynamics of the
1D Hubbard model, showing novel features of fractional holons and spinons. In
the Mott phase, the excitations show a charge gapped phase with
antiferromagnetic correlations in spin sector. Spin and charge particle-hole
excitations reveal the origins of the spin-coherent TLL and spin-incoherent
liquid in a crossover temperature regime. The purpose of this section is aimed
to calculate such typical excitation spectra of the repulsive Hubbard model,
which complement the study in literature, see a review essler2005one .

###  II.1 II.1 Bethe ansatz equations and root patterns

In the presence of external potentials, the 1D Hamiltonian of Hubbard model is
given by

$H=-\sum_{j=1}^{L}\sum_{a=\uparrow\downarrow}\left(c_{j,a}^{+}c_{j+1,a}+c_{j+1,a}^{+}c_{j,a}\right)+u\sum_{j=1}^{L}(1-2n_{j\uparrow})(1-2n_{j\downarrow})-\mu\hat{N}-2B\hat{S}_{z},$
(1)

where $c_{j,a}^{\dagger}$ and $c_{j,a}$ are the creation and annihilation
operators of fermions with spin $a$ ($a=\uparrow$ or $a=\downarrow$) at site
$j$ in a periodic lattice of length $L$, $n_{j\sigma}$ is particle density
operator of spin $\sigma$ at site $j$,
$\hat{N}=\sum_{j=1}^{L}(n_{j\uparrow}+n_{j\downarrow})$ is the total number
particle operator, and
$\hat{S}_{z}=\frac{1}{2}\sum_{j=1}^{L}(n_{j\uparrow}-n_{j\downarrow})$ is the
magnetization operator. The hopping amplitude is set to unity, and $u$
represents an on-site interaction between particles. The chemical potential
$\mu$ and magnetic field $B$ are rescaled with respect to the hopping term.
$u>0$ for repulsion and $u<0$ for attraction, they can be converted into each
other by the Shiba transformation. We choose a unit system by setting
$k_{B},\,\hbar$ to $1$ throughout the paper.

At an absence of the external potential ($\mu=0,B=0$) and even lattice sites,
the Hamiltonian 1 possesses SO(4)$\cong$SU(2)$\times$ SU(2)$/\mathbb{Z}_{2}$
full symmetry composed of spin rotational and $\eta$-pairing
invarianceEss1994a :

${[H,S^{\alpha}]}=0,\qquad[H,\eta^{\alpha}]=0,$ (2)

where the spin and $\eta$-pair operators are given by

$\displaystyle S^{\alpha}=$
$\displaystyle\frac{1}{2}\sum_{j=1}^{L}\sum_{a,b=\uparrow,\downarrow}c_{j,a}^{+}(\sigma^{\alpha})_{b}^{a}c_{j,b},$
(3) $\displaystyle\eta^{x}=$
$\displaystyle-\frac{1}{2}\sum_{j=1}^{L}(-1)^{j}\left(c_{j,\uparrow}^{+}c_{j,\downarrow}^{+}+c_{j,\uparrow}c_{j,\downarrow}\right),$
$\displaystyle\eta^{y}=$
$\displaystyle\frac{i}{2}\sum_{j=1}^{L}(-1)^{j}\left(c_{j,\uparrow}^{+}c_{j,\downarrow}^{+}-c_{j,\uparrow}c_{j,\downarrow}\right),$
$\displaystyle\eta^{z}=$ $\displaystyle\frac{1}{2}\left(N-L\right),$

which are useful to analyze elementary excitations at zero magnetic field and
half-filled band Ess1994a . In the above equation, $\sigma^{\alpha}$ with
$\alpha=x,y,z$ denote the Pauli matrices.

In this paper, we will consider the case where the number of particles $n_{c}$
is less than one with the magnetization less than $n_{c}/2$ for repulsive
interaction, i.e., it corresponds to a negative chemical potential and
positive magnetic field. The other argument spaces can be involved in terms of
discrete symmetries essler2005one .

The eigenvalues of {$k$}, {$\Lambda$}, respectively the solutions of charge
quasimomentum and spin rapidity of this model, can be obtained by solving the
Lieb-Wu equations Lieb:1968 . While Takahashi takahashi1972one  assumed that
the patterns of quasimomenta {$k$}, {$\Lambda$} exhibit different patterns,
which determine full states of the model. The root patterns of repulsive
Hubbard model can be categorized into three types of strings: single $k$,
length-$n$ $\Lambda$ strings composed of $n$ spin-down electrons, length-$m$
$k-\Lambda$ string containing $m$ spin-down and $m$ spin-up particles. Let
$M_{e},M_{n},M^{\prime}_{n}$ denote the number of $k$, $\Lambda$ strings and
$k-\Lambda$ strings of length $n$, resepctively. Therefore, the total particle
number $N$ and spin down electron number $M$ read:

$\displaystyle M=\sum_{n=1}^{\infty}n\left(M_{n}+M_{n}^{\prime}\right),$ (4)
$\displaystyle N=\mathcal{M}_{e}+\sum_{n=1}^{\infty}2nM_{n}^{\prime}.$

In order to carry easily out our study, we prefer to present the following
Bethe ansatz equations and the TBA equations of the 1D Hubbard model with a
repulsive interaction. Based on these, we will derive analytically physical
properties of the model. In terms of the string solutions, the real centers of
these roots satisfy the so-called Takahashi’s forms of Bethe ansatz equations
takahashi1972one :

$\displaystyle k_{j}L=2\pi
I_{j}-\sum_{n=1}^{\infty}\sum_{\alpha=1}^{M_{n}}\theta\left(\frac{\sin
k_{j}-\Lambda_{\alpha}^{n}}{nu}\right)-\sum_{n=1}^{\infty}\sum_{\alpha=1}^{M_{n}^{\prime}}\theta\left(\frac{\sin
k_{j}-\Lambda_{\alpha}^{\prime n}}{nu}\right),$ (5)
$\displaystyle\sum_{j=1}^{N-2M^{\prime}}\theta\left(\frac{\Lambda_{\alpha}^{n}-\sin
k_{j}}{nu}\right)=2\pi
J_{\alpha}^{n}+\sum_{m=1}^{\infty}\sum_{\beta=1}^{M_{m}}\Theta_{nm}\left(\frac{\Lambda_{\alpha}^{n}-\Lambda_{\beta}^{m}}{u}\right),$
$\displaystyle 2L\operatorname{Re}\left[\arcsin\left(\Lambda_{\alpha}^{\prime
n}+n\mathrm{i}u\right)\right]=2\pi J_{\alpha}^{\prime
n}+\sum_{j=1}^{N-2M^{\prime}}\theta\left(\frac{\Lambda_{\alpha}^{\prime
n}-\sin
k_{j}}{nu}\right)+\sum_{m=1}^{\infty}\sum_{\beta=1}^{M_{m}^{\prime}}\Theta_{nm}\left(\frac{\Lambda_{\alpha}^{n}-\Lambda_{\beta}^{\prime
m}}{u}\right),$

where $\theta(x)=2\arctan(x)$ and $\Theta_{nm}$ is defined as

$\Theta_{nm}(x)=\left\\{\begin{array}[]{l}\theta\left(\frac{x}{|n-m|}\right)+2\theta\left(\frac{x}{|n-m|+2}\right)+\cdots+2\theta\left(\frac{x}{n+m-2}\right)+\theta\left(\frac{x}{n+m}\right),\text{
if }n\neq m\\\
2\theta\left(\frac{x}{2}\right)+2\theta\left(\frac{x}{4}\right)+\cdots+2\theta\left(\frac{x}{2n-2}\right)+\theta\left(\frac{x}{2n}\right),\text{
if }n=m\end{array}\right..$ (6)

The counting numbers $I_{j},J_{\alpha}^{n},J_{\alpha}^{\prime n}$ are integer
or half-odd integers, which rely on the odevity of string number,

$\displaystyle I_{j}\text{ is }\begin{cases}\text{ integer }&\text{ if
}\sum_{m}\left(M_{m}+M_{m}^{\prime}\right)\text{ is even }\\\ \text{ half-odd
integer }&\text{ if }\sum_{m}\left(M_{m}+M_{m}^{\prime}\right)\text{ is odd
}\end{cases},$ (7) $\displaystyle J_{\alpha}^{n}\text{ is }\begin{cases}\text{
integer }&\text{ if }N-M_{n}\text{ is odd }\\\ \text{ half-odd integer
}&\text{ if }N-M_{n}\text{ is even, }\end{cases},$ (8) $\displaystyle
J_{\alpha}^{\prime n}\text{ is }\begin{cases}\text{ integer }&\text{ if
}L-N+M_{n}^{\prime}\text{ is odd }\\\ \text{ half-odd integer }&\text{ if
}L-N+M_{n}^{\prime}\text{ is even }\end{cases}.$ (9)

The classification of these quantum numbers will be needed for characterizing
the excitations which will be presented later.

Every selected set of quantum numbers $I_{j},J_{\alpha}^{n},J_{\alpha}^{\prime
n}$ uniquely determine the value of quasimomentum
$k_{j},\,\Lambda_{\alpha}^{n},\,\Lambda_{\alpha}^{\prime n}$, and thus give a
specific state, either the ground state or an excited state of various filling
of the model. These numbers are in the ranges of

$\displaystyle-\frac{L}{2}<I_{j}\leq\frac{L}{2},$ (10)
$\displaystyle\left|J_{\alpha}^{n}\right|\leq\frac{1}{2}\left(N-2M^{\prime}-\sum_{m=1}^{\infty}t_{nm}M_{m}-1\right),$
$\displaystyle\left|J_{\alpha}^{\prime
n}\right|\leq\frac{1}{2}\left(L-N+2M^{\prime}-\sum_{m=1}^{\infty}t_{nm}M_{m}^{\prime}-1\right),$

where $M^{\prime}=\sum_{n=1}^{\infty}nM_{n}^{\prime}$ is the entire number of
spin-down particles included in $k-\Lambda$ strings, and $\
t_{mn}=2\min(m,n)-\delta_{mn}$.

Denote $\rho^{p},\sigma_{n}^{p},\sigma_{n}^{\prime p}$
($\rho^{h},\sigma_{n}^{h},\sigma_{n}^{\prime h}$) as the densities of real
quasimomenta $k$, the real part of the length-$n$ spin strings and the real
part of length-$n$ $k-\Lambda$ strings for particles (holes), the root
distributions of these types of root patterns are given by takahashi1972one

$\displaystyle\rho^{p}(k)+\rho^{h}(k)=$ $\displaystyle\frac{1}{2\pi}+\cos
k\sum_{n=1}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}\Lambda
a_{n}(\Lambda-\sin k)\left(\sigma_{n}^{\prime
p}(\Lambda)+\sigma_{n}^{p}(\Lambda)\right),$ (11)
$\displaystyle\sigma_{n}^{p}(\Lambda)+\sigma_{n}^{h}(\Lambda)=$
$\displaystyle-\sum_{m=1}^{\infty}A_{nm}*\sigma_{m}^{p}(\Lambda)+\int_{-\pi}^{\pi}\mathrm{d}ka_{n}(\sin
k-\Lambda)\rho^{p}(k),$ $\displaystyle\sigma_{n}^{\prime
p}(\Lambda)+\sigma_{n}^{\prime h}(\Lambda)=$
$\displaystyle\frac{1}{\pi}\operatorname{Re}\frac{1}{\sqrt{1-(\Lambda-\mathrm{i}nu)^{2}}}-\sum_{m=1}^{\infty}A_{nm}*\sigma_{m}^{\prime
p}(\Lambda)-\int_{-\pi}^{\pi}\mathrm{d}ka_{n}(\sin k-\Lambda)\rho^{p}(k),$

respectively. In the above equations
$a_{n}(x)=\frac{1}{2\pi}\frac{2nu}{(nu)^{2}+x^{2}}$, and the convolution term
denoted by $*$ denotes a convolution

$A_{nm}*f(x)=\int_{-\infty}^{\infty}\frac{\mathrm{d}y}{2\pi}\frac{\mathrm{d}}{\mathrm{d}x}\Theta_{nm}\left(\frac{x-y}{u}\right)f(y).$
(12)

The function $A_{nm}$ denotes a derivative of the function $\Theta_{nm}$,
namely

$\displaystyle
A_{nm}\left(\frac{x-y}{u}\right)=\frac{1}{2\pi}\frac{\mathrm{d}}{\mathrm{d}x}\Theta_{nm}\left(\frac{x-y}{u}\right)=$
(13)
$\displaystyle\left\\{\begin{array}[]{l}a_{|n-m|}\left(x-y\right)+2a_{|n-m|+2}\left(x-y\right)+\cdots+2a_{n+m-2}\left(x-y\right)+a_{n+m}\left(x-y\right),\text{
if }n\neq m\\\
2a_{2}\left(x-y\right)+2a_{4}\left(x-y\right)+\cdots+2a_{2n-2}\left(x-y\right)+a_{2n}\left(x-y\right),\text{
if }n=m\end{array}\right..$

Following the Yang-Yang method YangCN:1969 , the true physical equilibrium
states at finite temperature can be determined by the minimization of the free
energy with respect to the densities. The TBA equations for dressed energies
$\kappa(k),\varepsilon_{n}(\Lambda),\varepsilon_{n}^{\prime}(\Lambda)$,
associated with the product of the logarithm of the ratio of hole density to
particle density and temperature
$\kappa=T\ln\left(\frac{\rho^{h}}{\rho^{p}}\right),\varepsilon_{n}=T\ln\left(\frac{\sigma_{n}^{h}}{\sigma_{n}^{p}}\right),\varepsilon_{n}^{\prime}=T\ln\left(\frac{\sigma_{n}^{\prime
h}}{\sigma_{n}^{\prime p}}\right)$, are given by the following form
takahashi1972one

$\displaystyle\kappa(k)$ $\displaystyle=$ $\displaystyle-2\cos
k-\mu-2u-B+\sum_{n=1}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}\Lambda
a_{n}(\sin
k-\Lambda)T\ln\left(1+\mathrm{e}^{-\frac{\varepsilon_{n}^{\prime}(\Lambda)}{T}}\right)$
(14)
$\displaystyle-\sum_{n=1}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}\Lambda
a_{n}(\sin
k-\Lambda)T\ln\left(1+\mathrm{e}^{-\frac{\varepsilon_{n}(\Lambda)}{T}}\right),$
$\displaystyle\varepsilon_{n}(\Lambda)$ $\displaystyle=$ $\displaystyle
2nB-\int_{-\pi}^{\pi}\mathrm{d}k\cos ka_{n}(\sin
k-\Lambda)T\ln\left(1+\mathrm{e}^{-\frac{\kappa(k)}{T}}\right)+\sum_{m=1}^{\infty}A_{nm}*T\ln\left(1+\mathrm{e}^{-\frac{\varepsilon_{m}(\Lambda)}{T}}\right),$
(15) $\displaystyle\varepsilon_{n}^{\prime}(\Lambda)$ $\displaystyle=$
$\displaystyle
4\textmd{Re}\sqrt{1-(\Lambda-\mathrm{i}nu)^{2}}-2n\mu-4nu-\int_{-\pi}^{\pi}\mathrm{d}k\cos
ka_{n}(\sin k-\Lambda)T\ln\left(1+\mathrm{e}^{-\frac{\kappa(k)}{T}}\right)$
(16)
$\displaystyle+\sum_{m=1}^{\infty}A_{nm}*T\ln\left(1+\mathrm{e}^{-\frac{\varepsilon_{m}^{\prime}(\Lambda)}{T}}\right).$

It is observed that the above dressed energies are even functions with respect
to their variables and admit monotonically increasing behavior in the positive
argument space. The Gibbs free energy per site is given by

$\displaystyle f$ $\displaystyle=$
$\displaystyle-T\int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi}\ln\left(1+\mathrm{e}^{-\frac{\kappa(k)}{T}}\right)+u$
(17)
$\displaystyle-T\sum_{n=1}^{\infty}\int_{-\infty}^{\infty}\frac{\mathrm{d}\Lambda}{\pi}\operatorname{Re}\frac{1}{\sqrt{1-(\Lambda-\mathrm{i}nu)^{2}}}\ln\left(1+\mathrm{e}^{-\frac{\varepsilon_{n}^{\prime}(\Lambda)}{T}}\right).$

(17) serves as the equation of state from which one can directly obtain the
profiles of ground state and thermodynamic properties of the Hubbard model at
zero and finite temperatures by solving (14), (15) and (16). However, finding
the solution to the TBA equations (14), (15) and (16) imposes a big
theoretical challenge. Analytical results of universal thermodynamics,
Luttinger liquid properties, magnetic properties, scaling functions in the
vicinities of the critical points, as well as behaviour of the Mott phase,
quantum cooling are still lacking of a comprehensive study in literature.

###  II.2 II.2 Ground state

In the limit of zero temperature, one can see that the dressed energies of the
length-$n$ $\Lambda$ strings with $n>1$ and all $k-\Lambda$ are always
positive in the parameter space. Consequently, they do not contribute to zero
temperature thermal and magnetic properties. The real $k$ and length-$1$
$\Lambda$ string should be cut off at the points where dressed energies switch
sign. Thereby, in the ground state, the coupled equations for dressed energies
and densities and hole distribution functions Lieb:1968  are greatly
simplified

$\displaystyle\kappa(k)$ $\displaystyle=$
$\displaystyle-2\cos{k}-\mu-2u-B+\int_{-A}^{A}\mathrm{d}\Lambda
a_{1}(\sin{k}-\Lambda)\varepsilon_{1}(\Lambda),$ (18)
$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$ $\displaystyle
2B+\int_{-Q}^{Q}\mathrm{d}k\cos
ka_{1}(\sin{k}-\Lambda)\kappa(k)-\int_{-A}^{A}\mathrm{d}\Lambda^{\prime}a_{2}\left(\Lambda-\Lambda^{\prime}\right)\varepsilon_{1}(\Lambda^{\prime}),$
(19) $\displaystyle\rho^{p}(k)$ $\displaystyle=$
$\displaystyle\theta_{H}(Q-|k|)\left[\frac{1}{2\pi}+\cos
k\int_{-A}^{A}\mathrm{d}\Lambda a_{1}(\sin
k-\Lambda)\sigma_{1}^{p}(\Lambda)\right],$
$\displaystyle\sigma_{1}^{p}(\Lambda)$ $\displaystyle=$
$\displaystyle\theta_{H}(A-|\Lambda|)\left[-\int_{-A}^{A}\mathrm{d}\Lambda^{\prime}a_{2}\left(\Lambda-\Lambda^{\prime}\right)\sigma_{1}^{p}\left(\Lambda^{\prime}\right)+\int_{-Q}^{Q}\mathrm{d}ka_{1}(\sin
k-\Lambda)\rho^{p}(k)\right],$ $\displaystyle\rho^{h}(k)$ $\displaystyle=$
$\displaystyle\theta_{H}(|k|-Q)\left[\frac{1}{2\pi}+\cos
k\int_{-A}^{A}\mathrm{d}\Lambda a_{1}(\sin
k-\Lambda)\sigma_{1}^{p}(\Lambda)\right],$
$\displaystyle\sigma_{1}^{h}(\Lambda)$ $\displaystyle=$
$\displaystyle\theta_{H}(|\Lambda|-A)\left[-\int_{-A}^{A}\mathrm{d}\Lambda^{\prime}a_{2}\left(\Lambda-\Lambda^{\prime}\right)\sigma_{1}^{p}\left(\Lambda^{\prime}\right)+\int_{-Q}^{Q}\mathrm{d}ka_{1}(\sin
k-\Lambda)\rho^{p}(k)\right],$ $\displaystyle\sigma_{n}^{h}(\Lambda)$
$\displaystyle=$ $\displaystyle-
A_{n1}*\sigma_{1}^{p}(\Lambda)+\int_{-Q}^{Q}\mathrm{d}ka_{n}(\sin
k-\Lambda)\rho^{p}(k),\,\,{\rm for},\,n\geq 2,$
$\displaystyle\sigma_{n}^{\prime h}(\Lambda)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\operatorname{Re}\frac{1}{\sqrt{1-(\Lambda-
inu)^{2}}}-\int_{-Q}^{Q}\mathrm{d}ka_{n}(\sin k-\Lambda)\rho^{p}(k),$

where $Q,\,A$ are the Fermi points of charge and spin, respectively, i.e.
$\kappa(Q)=0,\varepsilon_{1}(A)=0$. $\theta_{H}$ is the Heaviside step
function. We observe that particles only exist within Fermi points, while
holons locate outside in the ground state. In the above equations, we also
present the hole density distribution functions for $\Lambda$ strings and
$k-\Lambda$ strings for our later calculation. Defining the total quantities
$\rho=\rho^{p}+\rho^{k},\,\sigma_{1}=\sigma_{1}^{p}+\sigma_{1}^{h}$, thus the
density Bethe ansatz equations can be recast in a more succinct form

$\displaystyle\rho(k)$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi}+\cos
k\int_{-A}^{A}\mathrm{d}\Lambda a_{1}(\Lambda-\sin k)\sigma_{1}(\Lambda),$
(20) $\displaystyle\sigma_{1}(\Lambda)$ $\displaystyle=$
$\displaystyle-\int_{-A}^{A}\mathrm{d}\Lambda^{\prime}a_{2}(\Lambda-\Lambda^{\prime})\sigma_{1}(\Lambda^{\prime})+\int_{-Q}^{Q}\mathrm{d}ka_{1}(\sin
k-\Lambda)\rho(k).$ (21)

Thus, the particle number and spin-down number can be expressed as in terms of
(20) and (21)

$\displaystyle n_{c}$ $\displaystyle=$
$\displaystyle\int_{-Q}^{Q}\mathrm{d}k\rho(k),$ (22) $\displaystyle
n_{\downarrow}$ $\displaystyle=$
$\displaystyle\int_{-A}^{A}\mathrm{d}\Lambda\sigma_{1}(\Lambda).$ (23)

The momenta of the charge and spin quasiparticles for $|k|\leq
Q,\,|\Lambda|\leq A$, and holons for $|k|>Q,|\Lambda|>A$ in terms of the
quasi-momentum parameters
$k_{j},\Lambda_{\alpha}^{n},\,\Lambda_{\alpha}^{\prime n}$ can be calculated
via the following expressions

$\displaystyle p(k)$ $\displaystyle=$ $\displaystyle\frac{2\pi
I_{k}}{L}=2\pi\int_{0}^{k}\mathrm{d}k^{\prime}\rho\left(k^{\prime}\right),$
(24) $\displaystyle p_{1}(\Lambda)$ $\displaystyle=$ $\displaystyle\frac{2\pi
J_{1}}{L}=2\pi\int_{0}^{\Lambda}\mathrm{d}\Lambda^{\prime}\sigma_{1}\left(\Lambda^{\prime}\right),$
(25) $\displaystyle p_{n}(\Lambda)$ $\displaystyle=$ $\displaystyle\frac{2\pi
J_{n}}{L}=2\pi\int_{0}^{\Lambda}\mathrm{d}\Lambda^{\prime}\sigma_{n}^{h}\left(\Lambda^{\prime}\right),\
n\geq 2,$ (26) $\displaystyle p_{n}^{\prime}(\Lambda)$ $\displaystyle=$
$\displaystyle\frac{2\pi
J_{n}^{\prime}}{L}=-2\pi\int_{0}^{\Lambda}\mathrm{d}\Lambda^{\prime}\sigma_{n}^{\prime
h}\left(\Lambda^{\prime}\right)+\pi(n+1),\ n\geq 1.$ (27)

Let us denote $N_{GS}$ be the total particle number, and $M_{GS}$ the spin-
down number in the ground state. Consider the ground state structure in the
presence of magnetic field and $0<n_{c}<1$ as follows (assume lattice length
$L$ is even): $M_{e}=N_{GS}=2*M_{GS}<L$, $M_{1}=M_{GS}$ is an odd number and
less than $N_{GS}/2,~{}M_{n}=0$ for $n\geq 2,\,M^{\prime}_{n}=0$ for $n\geq
1$, that is, both charge and spin sectors are partially filled and all
particles form charge and spin Fermi seas, inside of which are referred as
antiholon states for charge and spinon states for spin degree of freedom. At
zero temperature, the system lies in its lowest energy state and particles are
symmetrically distributed around zero momentum. Using (7) (8) and (24)-(27),
we can evaluate that $\left\\{I_{j}\right\\}$ are half-odd integers and fill
in the interval $\left[-\frac{N_{GS}-1}{2},\frac{N_{GS}-1}{2}\right]$, whereas
$\left\\{J_{\alpha}^{1}\right\\}$ are integer within
$\left[-\frac{M_{GS}-1}{2},\frac{M_{GS}-1}{2}\right]$. The momentum
distribution $p(k),\,p_{1}(\Lambda)$ respectively covers the intervals of
$\left[-\pi n_{c},\pi n_{c}\right]$ and $\left[-\pi n_{\downarrow},\pi
n_{\downarrow}\right]$. It is inferred from (10) that the values for quantum
numbers $\left\\{J_{\alpha}^{1}\right\\}$ satisfiy
$\left|J_{\alpha}^{1}\right|\leq\frac{1}{2}\left(N_{GS}-M_{GS}-1\right)$. This
implies that there are $N_{GS}-M_{GS}$ vacancies for $J_{\alpha}^{1}$ and thus
$N_{GS}-2M_{GS}$ holes are left. Building on the quantum numbers of the ground
state, elementary excitations can be classified into two types: gapless and
gapped excitations. Particle-hole excitations as well as their combinations
belong to the gapless ones, which always occur at sufficient low energy
sector. Whereas the length-$n$ $k-\Lambda$ and $\Lambda$ strings with $n>1$
require an additional energy to be excited. The key points for excitation
spectrum calculation are stated below:

  * •

dressed energy is an exact excitation energy of the particle with the
corresponding momentum, and holon energy corresponds to the negative dressed
energy. In other words, energies of particle excitations for different types
of quasiparticles are $\kappa(k)$, $\varepsilon_{n}(\Lambda)$ and
$\varepsilon^{\prime}_{n}(\Lambda)$, while their corresponding holon
excitations are $-\kappa(k)$, $-\varepsilon_{n}(\Lambda)$ and
$-\varepsilon^{\prime}_{n}(\Lambda)$, respectively.

  * •

by changing quantum numbers $M_{e}$, $M_{n}$ and $M^{\prime}_{n}$, one can
analyze the excited structure with respect to ground state for capturing the
holon or particle excitations of each type strings.

  * •

if there is no extra hole or particle appearing over the ground state, the
particle-hole excitation can exist. This type of excitation does not alter the
total particle number or parity.

  * •

individual particle and holon excitations can be transformed into each other,
e.g., in figure 1, a holon excitation can be created by conducting two
particle-hole excitations over a particle excitation.

Figure 1: Illustration of two equivalent excitations: a transformation between
particle excitation and holon excitation.

  * •

when $I_{j}$ converts from half-odd (interger) to interger (half-odd), extra
constant momentum $\pm\pi n_{c}$ needs to be added to total momentum with a
change of parity (In the later study, we only consider the case of $\pi n_{c}$
which can find common in continuous model. The case for $-\pi n_{c}$ is
symmetric with the former.), see, for example, figure 2.

Figure 2: Illustration of the variation of parity in excitation.

  * •

a multi-parametric excitation can be decomposed into few-parametric elementary
excitations. Once few particle excitations are determined, analytical
expressions of multiple particle excitations are created by few elementary
excitations. We elucidate this approach in later, e.g., $M_{e}=N_{GS}+2$ is
compose of double excitations with $M_{e}=N_{GS}+1$.

The techniques mentioned above enable us to classify the processes of
complicated excitations, based on which we can conceive unique collective
nature of fractional excitations. In what follows we will present various
spectrum patterns of elementary excitations, which are helpful for latter
analysis.

###  II.3 II.3 Fractional charge and spin excitations

Figure 3: Elementary excitations in spin and charge. The shaded area in orange
represents a particle-hole excitation of charge, where one particle inside the
Fermi sea is excited to outside at different fillings. Whereas the shaded area
in green presents a two-holon excitation in spin sector. All figures are drawn
in the first Brillouin zone. Interaction $u$ is fixed to unity for all figures
except (b): (a) For the case away from half-filled band at zero magnetic field
with $B=0,\mu=-1.5765$, we have charge and spin sound velocities
$v_{c}=2.0907,\,v_{s}=1.2259$, respectivly; (b) For strong coupling at zero
magnetic field with $B=0,\mu=-19.22,u=10$, we have charge and spin sound
velocities $v_{c}=1.8533,v_{s}=0.1824$, respectively. (c) presents a general
case with a presence of magnetization, $B=0.555,\mu=-1.32$; (d) presents
charge and spin excitations near half filled phase with a finite
magnetization, $B=0.555,\mu=-0.722$; (e) shows excitations near a spin-
polarized band, $B=0.555,\mu=-1.77$. At a low energy scale, holon excition
energy (empty-circles) is higher than particle excitation energy (filled-
circles); (f) shows excitations at $B=0.3,\mu=-2.5$ for a comparison with the
figure (e).

Elementary excitations of the Hubbard model were studied in literature, see
review essler2005one . Here we present a comprehensive understanding of
fractional excitations and various gapped excitations which provide deep
insight into understanding Luttinger liquid, spin incoherent liquid and low
energy physics.

####  II.3.1 Two holon excitations

We first study the elementary excitation of flipping one spin-down particle,
for which we have a change $M_{e}=L,\,M_{1}=M_{GS}-1$. Here we take $M_{GS}$
an odd number. From (7), (8), we can determine that the quantum number
$\left\\{I_{j}\right\\}$ are integers and $J_{\alpha}^{1}$ are half-odd
integers within ranges $-\frac{L}{2}<I_{j}\leq\frac{L}{2}$ and
$\left|J_{\alpha}^{1}\right|\leq\frac{1}{2}\left(N_{GS}-M_{GS}\right)$,
respectively. Therefore the vacancies for spin are $N_{GS}-M_{GS}+1$. In
contrast to the ground state, one more vacancy arises with a diminished
particle in the spin sector, and thus two holes $\Lambda_{h1},\,\Lambda_{h2}$
carrying separately spin $\frac{1}{2}$ are created in the spinon band for
$\Lambda$. For this spin flipping excitation, $\eta$ pairing and spin
magnetization $\left(\Delta\eta^{z},\Delta S^{z}\right)$ are determined
through $\Delta\eta^{z}=\frac{\Delta N}{2},\,\Delta S^{z}=\frac{\Delta
N-2\Delta M}{2}$ with $N,\,M$ total particle number and spin-down number in
the excited state, i.e., $\left(\Delta\eta^{z},\Delta S^{z}\right)=(0,1)$ in
this sample.

On the other hand, the particle-hole excitation can be constructed in the
charge sector with a particle taken out from $|k|<Q$ to $|k|>Q$, leaving a
holon $k_{h}$ inside the antiholon band and a particle $k_{p}$ in the holon
states. The energies and momenta of the particle-hole excitation in the charge
sector and the two-spinon excitation in the spin sector are

$\displaystyle E_{ph}$ $\displaystyle=$
$\displaystyle\kappa(k_{p})-\kappa(k_{h}),\qquad P_{ph}=p(k_{p})-p(k_{h}),$
(28) $\displaystyle E_{ss}$ $\displaystyle=$
$\displaystyle-\varepsilon_{1}(\Lambda_{h1})-\varepsilon_{1}(\Lambda_{h2}),\qquad
P_{ss}=-p_{1}(\Lambda_{h1})-p_{1}(\Lambda_{h2})+\pi n_{c}.$ (29)

In figure 3, we demonstrate the results of the particle-hole and the two-
spinon excitations for various parameters. Obviously, one universal nature of
low-energy excitation spectra is the existence of spin-charge separation
behavior voit1993charge ; kollath2005spin ; lorenz2002evidence , visualizing
the theory of the spin-charge separated TLL. In the absence of the magnetic
field, the system has spin rotation symmetry. Figure 3 (a) shows the particle-
hole excitation in the charge sector when away from half-filled lattice and
the two-spinon excitation in the spin sector which fully accounts for the
first Brillouin zone.

In contrast, Figure 3 (b) shows the particle-hole excitation and the two-
spinon excitation for strong coupling limit, reminiscent of the spin-charge
separated spectra in the 1D Yang-Gaudin model he2020emergence . The slopes of
the dispersion are characterised by velocities $v_{c},\,v_{s}$, displaying
$v_{s}\ll v_{c}$, which were recently measured experimentally in a quasi-1D
trapped repulsive Fermi gas senaratne2022spin . This regime offers a
remarkable possibility to study the spin-incoherent Luttinger liquid (SILL),
which occurs when the temperature exceeds the spin characteristic energy, but
is still less than the charge energy. In the SILL, the spin degree of freedom
is frozen. However, the charge degree of freedom propagates as a sound mode,
behaving like spinless fermions cheianov2004nonunitary ; fiete2004green ;
fiete2007colloquium .

Figure 3 (c), (d) presents the particle-hole excitation and the two-spinon
excitation for an arbitrary magnetic field, away from half-filling and near
half-filling, respectively. In these two cases, the spin degree is highly
suppressed in the presence of the magnetic field. The spin and charge
velocities are significantly different, and the spectra of spin and charge
display a large separation. Exspecially for figure 3 (d), the charge band is
very narrow, behaving like a single particle, which was experimentally studied
in vijayan2020time . The spin velocity exceeds the sound velocity of charge.
Figure 3 (e), (f) shows that in the charge sector the energy of holon
excitation at small momentum can exceed that of particle excitation in the
Hubbard model due to the cosine term presented in the charge dispersion, which
is prohibited in continuous modelgiamarchi2003quantum ; he2020emergence .
Nevertheless, when the system approaches the Mott phase, charge excitation
gradually shrinks to the cases of the figure 3 (d), indicating an emergence of
single particle excitation in the charge sector vijayan2020time . In general,
the spin sector exhibits similar behaviour as that of Heisenberg
chainovchinnikov1969excitation ; nocera2016magnetic . When magnetic field is
applied, the two-spinon excitation is no longer the low-energy one in the
sense that it does not fill the entire Brillouin area. The role played by the
low-energy excitation in spin degree of freedom is the spin particle-hole
spectrum as depicted in figure 4.

Figure 4: The shaded area in orange represents a particle-hole excitation of
charge. The blue shaded area denotes spin particle-hole excitation. The green
color denotes the two-spinon excitation, same as figure 3 (c). Parameter
setting: $B=0.555,\mu=-1.32,u=1$.

####  II.3.2 Antiholon- and holon-spinon excitations

Let’s further consider adding or removing a particle with spin up or down over
the ground state in the Hubbard model. There are four configurations with the
following quantum numbers: $(a)M_{e}=N_{GS}+1,M_{1}=M_{GS}+1$; (b)
$M_{e}=N_{GS}-1,M_{1}=M_{GS}$; (c) $M_{e}=N_{GS}+1,M_{1}=M_{GS}$; (d)
$M_{e}=N_{GS}-1,M_{1}=M_{GS}-1$. These four excitations can be classified by
the charge and spin magnetization values $\left(\Delta\eta^{z},\Delta
S^{z}\right)$: correspondingly, (a) antiholon-spinon particle
$(\frac{1}{2},-\frac{1}{2})$, (b) holon-spinon particle
$(-\frac{1}{2},-\frac{1}{2})$, (c) antiholon-spinon hole
$(\frac{1}{2},\frac{1}{2})$, and (d) holon-spinon hole
$(-\frac{1}{2},\frac{1}{2})$. Consequently, (a) has one antiholon particle in
the charge and one spinon particle in the spin sector; (b) has one holon in
the charge and one spinon particle in the spin sector; (c) has an antiholon
particle in the charge and one spinon hole in the spin sector; (d) has one
holon in the charge and one spinon hole in the spin sector. Analogous to
previous discussion, we obtain the following results of excitation spectra

$\displaystyle\begin{array}[]{ll}(a)&\left\\{\begin{array}[]{l}E_{\overline{h}\overline{s}}=\kappa(k_{p})+\varepsilon_{1}(\Lambda_{p})\\\
P_{\overline{h}\overline{s}}=p(k_{p})+p_{1}(\Lambda_{p})+\pi
n_{c}\end{array},\right.\\\
(b)&\left\\{\begin{array}[]{l}E_{h\overline{s}}=-\kappa(k_{h})+\varepsilon_{1}(\Lambda_{p})\\\
P_{h\overline{s}}=-p(k_{h})+p_{1}(\Lambda_{p})\end{array},\right.\\\
(c)&\left\\{\begin{array}[]{l}E_{\overline{h}s}=\kappa(k_{p})-\varepsilon_{1}(\Lambda_{h})\\\
P_{\overline{h}s}=p(k_{p})-p_{1}(\Lambda_{h})\end{array},\right.\\\
(d)&\left\\{\begin{array}[]{l}E_{hs}=-\kappa(k_{h})-\varepsilon_{1}(\Lambda_{h})\\\
P_{hs}=-p(k_{h})-p_{1}(\Lambda_{h})+\pi n_{c}\end{array}.\right.\end{array}$
(42)

In the light of the above discussion, these four excitation patterns can be
converted to each other through certain numbers of particle-hole excitations.
Their continuum spectra are depicted in figure 5. Comparing upper (low spin-
down density) and lower (near half-filling) panels of the figure 5, we observe
that the cases (a)(c)(d), which are related to the antiholon excitations in
charge sector or spinon hole excitations in spin part, are severely affected
by the filling and magnetization. When spin sector tends to vanish or charge
sector becomes gapped, the profiles of excitations behave like single particle
behaviours.

Figure 5: From left to right: (a) antiholon-spinon particle excitation
$(\frac{1}{2},-\frac{1}{2})$, (b) holon-spinon particle
excitation$(-\frac{1}{2},-\frac{1}{2})$, (c) antiholon-spinon hole excitation
$(\frac{1}{2},\frac{1}{2})$ and (d) holon-spinon hole excitation
$(-\frac{1}{2},\frac{1}{2})$. Upper panel: excitation spectra in the low spin-
down density area with the potentials $B=0.555,\mu=-1.85,u=1$ and density
$n_{c}=0.6609,n_{\downarrow}=0.0962$; Lower panel: excitation spectra near
half-filling with the potentials $B=0.555,\mu=-0.722,u=1$ and density
$n_{c}=0.9775,n_{\downarrow}=0.298$. In each panel, the excitation bands in
(a) and (d) are centered around the $\pi n_{c}$, whereas the excitation bands
in (b) and (c) are centered around $\pi$.

In general, we observe that the gapless excitations occur only at a
sufficiently small energy scale. The charge and spin excitations are well
decoupled only in certain configurations demonstrated in figure 3 and figure
4. This spin-charge separation phenomenon is unique in 1D strongly correlated
electron systems and can be described by the spin-coherent and spin-incoherent
TLL theory. We will further discuss the spin-coherent and spin-incoherent TLLs
in the 1D Hubbard model below.

###  II.4 II.4 $k-\Lambda$ string excitations

In comparison with the continuous systems, one peculiar excitation of 1D
Hubbard model is the gapped bound states called $k-\Lambda$ string, with
length-1 formed by two electrons with opposite spins due to lattice effect.
Such states make important contributions to optical conductivity close to
half-filling veness2016mobile . The continuum spectra are shaped like a
parabola with a convex opening, see the bottom of the excitation spectra in
figure 6. The length-$1$ $k-\Lambda$ string shows to be the lowest charge
gapped excitation among the other $k-\Lambda$ strings, leading to the
appearance of quasiparticle $\Lambda_{1}$ in the distribution of quantum
number $J_{\alpha}^{\prime 1}$. Such a gapped excitation also leads to the
change of the parity of $I_{j}$. It can be seen that the excited energy of
this length-$1$ $k-\Lambda$ has an energy and momentum
$E=\varepsilon^{\prime}_{1}(\Lambda_{1}),\,P=p^{\prime}_{1}(\Lambda_{1})+\pi
n_{c}$. For a more general consideration of $k-\Lambda$ excitations, one
$k-\Lambda$ string is added to the occupation numbers as discussed in section
Antiholon- and holon-spinon excitations, i.e., for antiholon-spinon
$M_{e}=N_{GS}+1,M_{1}=M_{GS}+1,M^{\prime}_{1}=1$, of which we decompose the
multi-parameter excitation into elementary ones. Similarly, we obtain the
information of such gapped excited states. The energy and momentum of this
compound case read

$\displaystyle(a)\,\left\\{\begin{array}[]{l}E_{\overline{h}\overline{s}}=\kappa(k_{p})+\varepsilon_{1}(\Lambda_{p})+\varepsilon^{\prime}_{1}(\Lambda_{1})\\\
P_{\overline{h}\overline{s}}=p(k_{p})+p_{1}(\Lambda_{p})+p^{\prime}_{1}(\Lambda_{1})\end{array},\right.$
(45)
$\displaystyle(b)\,\left\\{\begin{array}[]{l}E_{h\overline{s}}=-\kappa(k_{h})+\varepsilon_{1}(\Lambda_{p})+\varepsilon^{\prime}_{1}(\Lambda_{1})\\\
P_{h\overline{s}}=-p(k_{h})+p_{1}(\Lambda_{p})+p^{\prime}_{1}(\Lambda_{1})+\pi
n_{c}\end{array},\right.$ (48)
$\displaystyle(c)\,\left\\{\begin{array}[]{l}E_{\overline{h}s}=\kappa(k_{p})-\varepsilon_{1}(\Lambda_{h})+\varepsilon^{\prime}_{1}(\Lambda_{1})\\\
P_{\overline{h}s}=p(k_{p})-p_{1}(\Lambda_{h})+p^{\prime}_{1}(\Lambda_{1})+\pi
n_{c}\end{array},\right.$ (51)
$\displaystyle(d)\,\left\\{\begin{array}[]{l}E_{hs}=-\kappa(k_{h})-\varepsilon_{1}(\Lambda_{h})+\varepsilon^{\prime}_{1}(\Lambda_{1})\\\
P_{hs}=-p(k_{h})-p_{1}(\Lambda_{h})+p^{\prime}_{1}(\Lambda_{1})\end{array}.\right.$
(54)

The corresponding excitation spectra are plotted in figure 6. The length-$1$
$k-\Lambda$ string leads to gapped continuum bands over the states (a)
antiholon-spinon excitation, (b) holon-spinon excitation, (c) antiholon-spinon
excitation and (d) holon-spinon excitation.

Figure 6: From left to right: length-1 $k-\Lambda$ string is added to the
excitations given in case ’Antiholon- and holon-spinon excitations’: (a)
antiholon-spinon particle excitation, (b) holon-spinon particle excitation,
(c) antiholon-spinon hole excitation and (d) holon-spinon hole excitation. The
parameters are $B=0.555,\mu=-1.32,u=1$.

###  II.5 II.5 High spin string excitations

Now we discuss the excitations of high spin string. Due to the
antiferromagnetic ordering, the spin $\Lambda$ string bound states in the 1D
Heisenberg chain have received a great deal of interest essler2005one ;
WangZ:2018 ; molter2014bound . Like the spin root patterns in 1D repulsive
Fermi gas, the $\Lambda$ strings with length greater than $1$ are gapped too.
We consider the case of the quantum numbers $M_{1}=M_{GS}-1,\,M_{n}=1$ with
$n>1$, leading to an extra particle $\Lambda_{n}$ created in the parameter
space of $J_{\alpha}^{n}$. This configuration corresponds to the process of
flipping $n-1$ spin-up and simultaneously removing a spin-down from the
channel of $J_{\alpha}^{1}$ to constitute a length-$n$ $\Lambda$ string, and
thus one has the energy $\varepsilon_{n}(\Lambda_{n})$ and momentum
$p_{n}(\Lambda_{n})$ for the newly generated quasiparticle. This situation is
different from a multi-spinon-magnon excitation, where the quantum number of
down-spin $M$ is fixed essler2005one . Meanwhile we can apply a particle-hole
excitation in spin sector together with the length-$n$ $\Lambda$ string. In
general, its energy and momentum are given by

$\displaystyle E$ $\displaystyle=$
$\displaystyle-\varepsilon_{1}(\Lambda_{h})+\varepsilon_{1}(\Lambda_{p})+\varepsilon_{n}(\Lambda_{n}),$
(55) $\displaystyle P$ $\displaystyle=$ $\displaystyle-
p_{1}(\Lambda_{h})+p_{1}(\Lambda_{p})+p_{n}(\Lambda_{n}).$ (56)

In figure 7, we present the continuum excitation spectra for high spin strings
of length-$2$ and -$3$. The magnitudes of the energy gaps depend on the length
of the spinon bound states.

Figure 7: Excitation spectra of length-$n$ spinon strings. From left to right:
one-parametric excitation of length-2 $\Lambda$ string (orange line) and
length-3 $\Lambda$ string (green line); the length-$1$ $\Lambda$ particle-hole
excitation with a length-$2$ $\Lambda$ particle excitation; length-$1$
$\Lambda$ particle-hole excitation with a length-$3$ $\Lambda$ particle
excitation. The parameters are $B=0.555,\mu=-1.32,u=1$.

##  III III. Universal thermodynamics

The 1D Hubbard model lays out profound many-body physics at zero and finite
temperatures. It exhibits rich quantum phases and universal thermodynamics
when the temperature is much less than the Fermi energy. A comprehensive
understanding of the universal low temperature behaviour still remains
challenging due to complexity of the Bethe ansatz equations and its $N!$-many
terms of the wave function. Most study on thermal and magnetic properties of
the 1D Hubbard model has been carried out at zero temperature and half-filling
or in dilute limit due to the complexity and intricacy of spin and charge root
string patterns in the Bethe ansatz equations and the TBA equations. Even in
the low-temperature regime, there is no well established understanding of
universal behaviour, such as TLL and quantum scaling functions etc. Likely,
this can be resorted to numerical simulation by iterative means
kawakami1989thermodynamic ; usuki1990thermodynamic ; mott1968metal ;
qin2022hubbard ; ilievski2017ballistic . However, analytical results of the
TLL and quantum criticality are still elusive. In this section, we will
present our analytical and numerical study of universal thermodynamics and
quantum criticality of the 1D Hubbard model by means of the TBA equations near
quadruple critical point. In this regime, the model is rarely studied.

###  III.1 III.1 Dimensionless ratios and phase diagram

In this subsection, we study the dimensionless ratios and their application in
the 1D Hubbard model. Before doing so, let us first give a brief view of the
ground state phase diagram determined by equations (18) and (19). In the limit
of zero temperature, the dressed energies $\kappa(k)$ and
$\varepsilon_{1}(\Lambda)$ can be used to analyze the phase transitions and
thermodynamic quantities. In the grand canonical ensemble, the magnetic field
and chemical potential determine the integration boundaries $Q,\,A$ through
the conditions $\kappa(Q)=0,\varepsilon_{1}(A)=0$. The particle density
$n_{c}=\int_{-Q}^{Q}\mathrm{d}k\rho(k)$, the spin-down electrons per site
$n_{\downarrow}=\int_{-A}^{A}\mathrm{d}\Lambda\sigma_{1}(\Lambda)$ and the
magnetization $m=n_{c}/2-n_{\downarrow}$ can be obtained from the root
densities equations (20) and (21). The magnetic field and chemical potential
drive the system from one phase to another at zero temperature. Phase
transitions usually occur with the conditions $\kappa(0)=0$ or
$\varepsilon_{1}(0)=0$, or with the presence of an energy gap $\kappa(\pi)=0$
in part. Using these conditions we may give five phases of states
essler2005one ; takahashi1974low : (I) vacuum with $Q=0,\,A=0$ or $n_{c}=m=0$;
(II) partially filled and spin fully polarized state with $0<Q<\pi,\,A=0$ or
$0<n_{c}<1,\,m=n_{c}/2$; the boundary condition is given by $Q=0$ and $Q=\pi$;
(III) half-filled and spin fully polarized state with $Q=\pi,\,A=0$ or
$n_{c}=1,m=1/2$; (IV) partially filled and magnetized band with
$0<Q<\pi,\,0<A\leq\infty$ or $0<n_{c}<1,0\leq m<n_{c}/2$; the phase boundary
between phase (II) and (IV) is given by
$\kappa(0)<0<\kappa(\pi),\,\varepsilon_{1}(0)=0$; (V) half-filled, magnetized
band, also referred to as Mott insulator with $Q=\pi,\,0<A\leq\infty$ or
$n_{c}=1,0\leq m<n_{c}/2$, the phase boundary between (V) and (IV) is given by
$\kappa(\pi)=0,\varepsilon_{1}(0)<0$. Substituting each phase boundary
conditions into the dressed energy equations (18) and (19), we can
analytically determine all critical fields in the $B-\mu$ plane, also see
essler2005one . As a matter of fact, thermal properties feature dramatic
quantum fluctuations around a quantum critical point (QCP), whereas the
dimensionless ratios, the ratios between the fluctuations of two types of
sources, can become comparable to temperature, for instance the susceptibility
Wilson ratio $R^{\chi_{s}}_{w}$ Wilson:1975  and the compressibility Wilson
ratio $R^{\chi_{c}}_{w}$ guan2013wilson ; yu2016dimensionless

$R^{\chi_{s}}_{w}=\frac{4}{3}\left(\frac{\pi
k_{B}}{\mu_{B}g}\right)^{2}\frac{\chi_{s}}{C_{v}/T},\qquad
R^{\chi_{c}}_{w}=\frac{\pi^{2}k^{2}_{B}}{3}\frac{\chi_{c}}{C_{v}/T},$ (57)

describing the competition between magnetic fluctuation or particle number
fluctuation and thermal fluctuations, respectively. In the above expressions,
$\chi_{s}$ ($\chi_{c}$) is the magnetic susceptibility (compressibility),
respectively.

These ratios can be constants independent of temperature in the TLL phase.
Therefore the ground state phase diagram can be characterized by the
dimensionless ratios.

Moreover, the Grüneisen ratio kuchler2003divergence ; boehler1980experimental
, which was introduced by Eduard Grüneisen Gruneisen-AdP-1908  in the
beginning of 20th Century in the study of the effect of volume change of a
crystal lattice on its vibrational frequencies,

$\Gamma=\frac{V\frac{\mathrm{d}p}{\mathrm{d}T}|_{V,N}}{\frac{\mathrm{d}E}{\mathrm{d}T}|_{V,N}}=\frac{1}{T}\frac{\frac{\partial^{2}p}{\partial\mu^{2}}\frac{\partial
p}{\partial T}-\frac{\partial^{2}p}{\partial\mu\partial T}\frac{\partial
p}{\partial\mu}}{\frac{\partial^{2}p}{\partial\mu^{2}}\frac{\partial^{2}p}{\partial
T^{2}}-\left(\frac{\partial^{2}p}{\partial\mu\partial T}\right)^{2}},$ (58)

has been widely used to study the caloric effect of solids and phase
transitions associated with the changes of volume, chemical potential,
interaction and magnetic field. Similar to the magnetic Grüneisen ratio, which
quantifies the magnetocaloric effect in the refrigeration with the variation
of magnetic field, the interaction driven Grüneisen ratio will be studied
later, quantifying the caloric effect in the refrigeration driven by the
variation of the interaction strength.

A significant aspect of these ratios is their characterization of the TLL and
their scaling behaviour near critical points guan2013wilson ;
yu2016dimensionless ; yu2020gruneisen ; Cheng:2018A . In figure 8 (a), we
present the three-dimensional plot of the magnetic Wilson ratio in the
$B-\mu-R^{\chi_{s}}_{w}$ parameter space, while in figure 8 (b) we show the
Wilson ratio as a function of $\mu$ for fixed $B$. It is remarkable that the
magnetic Wilson ratio is suddenly enhanced near the QCP, distinguishing
different phases. Moreover, in terms of bosonization results of the magnetic
susceptibility, specific heat, Luttinger parameter and velocity, the Wilson
ratios are given by

$\displaystyle\text{II:}\quad R^{\chi_{s}}_{w}$ $\displaystyle\approx$
$\displaystyle 2,$ $\displaystyle\text{IV:}\quad R^{\chi_{s}}_{w}$
$\displaystyle\approx$ $\displaystyle 4v_{c}K_{s}/(v_{s}+v_{c}),$
$\displaystyle\text{V:}\quad R^{\chi_{s}}_{w}$ $\displaystyle\approx$
$\displaystyle 8K_{s},$ $\displaystyle\text{I and III:}\quad R^{\chi_{s}}_{w}$
$\displaystyle=$ $\displaystyle 0,$ (59)

where $K_{s}$ is the Luttinger parameter for spin, and $v_{c,s}$ are sound
velocities for charge and spin, respectively. These relations (59) are in
excellent agreement with figure 8.

Figure 8: (a) The 3D plot of the Wilson ratio $R^{\chi_{s}}_{w}$ maps out
phase diagram of the 1D Hubbard model at low temperatures, which significantly
marks different phases of states. The parameter setting reads: $T=0.005$ and
$u=1$. (b) Plot of Wilson ratio $R^{\chi_{s}}_{w}$ at fixed $B=0.4,u=1$. As
chemical potantial varies, four phases I II IV V can be reached.

The ground state of the 1D repulsive Hubbard model consists of five phases:
empty lattice I, partially filled and fully-polarized phase II, fully-filled
and -polarized phase III, partially-filled and -polarized phase IV, and fully-
filled and partially-polarized phase V (Mott insulator) in figure 8. Among
these, the phase IV shows the most abundant physics In this phase, charge and
spin degrees of freedom coexist and dramatically make up the phase of spin-
charge separated TLLs tomonaga1950remarks ; luttinger1963exactly ;
haldane1981luttinger , showing non-Landau Fermi liquid behaviour in one
dimension. In the TLL phase, all quasiparticles form collective motions and
thus decouple into two propagating modes with different different velocities
$v_{c},\,v_{s}$voit1993charge ; kollath2005spin ; lorenz2002evidence . which
rely on the value of the coupling and can be effectively estimated from TBA
equations. The separated phenomenon can also be observed from their excitation
spectra. In APPENDIX A, we show that the phase IV can be seen as two-component
free fluids characterised by additivity of free energy at low-energy level,
while phase II and V are representatives of single TLL:

$\displaystyle f$ $\displaystyle=$ $\displaystyle f_{0}-\frac{\pi
T^{2}}{6}\left(\frac{1}{v_{c}}+\frac{1}{v_{s}}\right)\quad\text{phase IV},$
(60) $\displaystyle f$ $\displaystyle=$ $\displaystyle f_{0}-\frac{\pi
T^{2}}{6}\frac{1}{v_{c}}\quad\text{phase II},$ (61) $\displaystyle f$
$\displaystyle=$ $\displaystyle f_{0}-\frac{\pi
T^{2}}{6}\frac{1}{v_{s}}\quad\text{phase V},$ (62)

where $v_{c}$ and $v_{s}$ are defined through
$v_{c}=\kappa^{{}^{\prime}}(Q)/(2\pi\rho(Q))$,
$v_{s}=\varepsilon_{1}^{{}^{\prime}}(A)/(2\pi\sigma_{1}(A))$ with cut-off
points $Q,A$ respectively. The relation between free energy and velocity is a
universal common nature enjoyed by a large family of systemshe2020emergence ;
guan2013fermi ; Cheng:2018A . Taking derivative of free energy with respect to
temperature, specific heat is directly derived $C_{v}=\pi
T/3\left(1/v_{c}+1/v_{s}\right),C_{v}=\pi T/(3v_{c}),C_{v}=\pi T/(3v_{s})$ for
IV,II and V, respectively. The above results (60)-(62) provide the universal
leading order correction of the temperature at low energy of interacting many-
body systems, i.e., characteristic of thermodynamics in the TLLs.

###  III.2 III.2 Contributions from $k-\Lambda$ string in phase III

In the low-temperature limit $\mu,B>T$, the low energy behaviours of the
system are immune to the gapped string excitations when the states are away
from the QCPs in figure 8. By analysing the TBA equations (14), (15) and (16),
we observe that it is reasonable to ignore the gapped string states and retain
only the gapless string states in both $\Lambda$ and $k-\Lambda$ strings at
low energy. The driving terms in these equations significantly determine the
contributions to the low energy states. The length-$n$ $\Lambda$ strings and
$k-\Lambda$ string have lesser contributions to the low energy states when $n$
is larger. Here we just calculate the contributions from the length-$1$
$\Lambda$ and $k-\Lambda$ strings in the phase III in which both are gapped.
The TBA equations suggest that the greater the value of the chemical potential
$\mu$ and magnetic field $B$ relative to the temperature $T$, the smaller the
contributions of the $k-\Lambda$ and $\Lambda$ strings. In the phase III, the
absolute value of chemical potential is small, while the magnetic field $B$ is
large. Therefore, without losing generality, we here consider the length-$1$
$k-\Lambda$ and $\Lambda$ strings at low temperature $|\mu|\sim T$ and $B>T$
in the phase III. After tedious iterations on free energy shown in APPENDIX B,
we obtain a close form of the free energy

$f=-\mu-
u-B+g_{\frac{3}{2}}T^{\frac{3}{2}}+g_{2}T^{2}+g_{\frac{5}{2}}T^{\frac{5}{2}}+O(T^{3},\mathrm{e}^{-\frac{4B}{T}},\mathrm{e}^{\frac{4\mu}{T}}),$
(63)

where the coefficients are given explicitly by

$\displaystyle g_{\frac{3}{2}}$ $\displaystyle=$
$\displaystyle-\frac{\lambda_{1}\pi^{\frac{1}{2}}}{\eta^{\frac{1}{2}}_{1}}\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}+\frac{f_{\frac{1}{2}}}{2\pi^{\frac{1}{2}}}\mathrm{e}^{\frac{2\mu}{T}}+\frac{f_{\frac{3}{2}}}{2\pi^{\frac{1}{2}}},$
(64) $\displaystyle g_{2}$ $\displaystyle=$
$\displaystyle-\left(\frac{\lambda_{1}f_{\frac{3}{2}}}{\eta^{\frac{1}{2}}_{1}u}+\frac{f_{\frac{1}{2}}}{2\eta^{\frac{1}{2}}_{1}\pi
u}\right)\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}+\frac{f_{\frac{1}{2}}f_{\frac{3}{2}}}{4\pi
u}\mathrm{e}^{\frac{2\mu}{T}},$ $\displaystyle g_{\frac{5}{2}}$
$\displaystyle=$
$\displaystyle-\left(\frac{\lambda_{1}f^{2}_{\frac{3}{2}}}{2\eta^{\frac{1}{2}}_{1}\pi^{\frac{1}{2}}u^{2}}+\frac{\lambda_{2}\pi^{\frac{1}{2}}}{4\eta^{\frac{3}{2}}_{1}}+\frac{f_{\frac{1}{2}}f_{\frac{3}{2}}}{2\eta^{\frac{1}{2}}_{1}\pi^{\frac{3}{2}}u^{2}}\right)\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}+\left(\frac{3f_{\frac{1}{2}}f^{2}_{\frac{3}{2}}}{32\pi^{\frac{3}{2}}u^{2}}+\frac{f_{\frac{3}{2}}}{32\pi^{\frac{1}{2}}}\right)\mathrm{e}^{\frac{2\mu}{T}}+\frac{f_{\frac{5}{2}}}{32\pi^{\frac{1}{2}}},$

where we denote the function
$f_{n}=\mathrm{Li}_{n}\left(-\mathrm{e}^{\frac{2-\mu-2u-B}{T}}\right)$, and
the parameters $\lambda_{1}=\frac{1}{\pi\sqrt{1+u^{2}}}$,
$\lambda_{2}=\frac{1-2u^{2}}{\pi(1+u^{2})^{5/2}}$,
$\eta_{0}=4(u-\sqrt{1+u^{2}})$, and $\eta_{1}=\frac{2}{(1+u^{2})^{3/2}}$. By
analysing (64), we observe that in the functions
$g_{\frac{3}{2},2,\frac{5}{2}}$, the terms involving exponents
$\mathrm{e}^{-\frac{2B+\eta_{0}}{T}}$ and $\mathrm{e}^{\frac{2\mu}{T}}$ come
from gapped length-$1$ $\Lambda$ and $k-\Lambda$ strings, respectively, while
the terms that do not contain any exponential terms stem from the
contributions of gapless charge $k$. In the above results, we have omitted the
contributions from the strings with the lengths more than length-$2$.

Figure 9: Free energy for gapped length-$1$ $\Lambda$ and $k-\Lambda$ strings.
The black and purple solid curves denote the numerical result and analytical
results (63) of the free energy, respectively. The red empty-circle-solid line
denotes the only contribution from the charge $k$ degree of freedom. The green
filled-circle-solid line denotes the free energy without $\Lambda-\Lambda$
strings, while the blue triangular-solid line denotes the free energy without
$k-\Lambda$ strings. It is obvious that the $\Lambda$ string contribution is
negligible.

From figure 9, we observe that our analytic expression of the free energy is
consistent with the numerical result obtained from the TBA equations (14),
(15) and (16). The contributions from the $k-\Lambda$ string are greater than
that from the $\Lambda$ strings due to the choice of parameters $B>T$ and
$|\mu|\sim T$. Moreover, from the absolute value of free energy, we observe
that the order of the $k-\Lambda$ string contributions is $O(10^{-4})$ as
temperature tends to $T=1$. In general, it can be safely ignored in the low
temperature behaviour. In the phase III, the contributions from the gapped
$\Lambda$ strings are very small and negligible too. But the length-$1$
$\Lambda$ string plays an important role in the phase IV and V. More detailed
calculations of the free energy for the phase III are given in APPENDIX B.

###  III.3 III.3 Magnetic properties at zero temperature near quadruple
critical point

The magnetic properties of the 1D Hubbard model at zero temperature provide a
benchmark for understanding universal low energy physics. In order to get
insightful analytical results of magnetic properties of the model, here we try
to analytically calculate of the magnetic properties under the conditions of
high density (Fermi point $Q$ tends to $\pi$) and high magnetization limits
($A$ tends to $0$), namely, the region in figure 8 around the quadruple point.
The physics in this region has not been studied thoroughly in the literature.
The main idea is to expand the TBA equations in terms of small quantities
$\delta=\pi-Q$ and $A$. To this end, we need to calculate the integral
interval of charge from the range $[0,Q]$ to $[0,\pi-Q]$ and thus the second
term of (19) can be re-expressed by substituting $k$ by $\pi-k$, namely

$\displaystyle\int_{-Q}^{Q}\mathrm{d}k\cos ka_{1}(\sin{k}-\Lambda)\kappa(k)$
$\displaystyle=$ $\displaystyle\int_{-Q}^{0}\mathrm{d}k\cos
k(a_{1}(\sin{k}+\Lambda)+a_{1}(\sin{k}-\Lambda))\kappa(k)$ $\displaystyle=$
$\displaystyle-4\operatorname{Re}\sqrt{1-(\Lambda-\mathrm{i}u)^{2}}+4u+\int_{0}^{\pi-Q}\mathrm{d}k\cos
k(a_{1}(\sin{k}+\Lambda)+a_{1}(\sin{k}-\Lambda))\kappa(\pi-k).$

Using the above expression, the spin dressed energy is rewritten in terms of
small integral boundaries

$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$
$\displaystyle\varepsilon^{0}(\Lambda)+\int_{0}^{\pi-Q}\mathrm{d}k\cos
k(a_{1}(\sin{k}+\Lambda)+a_{1}(\sin{k}-\Lambda))\kappa(\pi-k)$ (66)
$\displaystyle-\int_{-A}^{A}\mathrm{d}\Lambda^{\prime}a_{2}\left(\Lambda-\Lambda^{\prime}\right)\varepsilon_{1}(\Lambda^{\prime}),$

where
$\varepsilon^{0}(\Lambda)=2B-4\operatorname{Re}\sqrt{1-(\Lambda-\mathrm{i}u)^{2}}+4u$
is the leading term. We further write the dressed energies
$\kappa(\pi-k)=\sum_{i=0}^{\infty}\kappa^{i}(\pi-k)\delta^{i},\varepsilon_{1}(\Lambda)=\sum_{i=0}^{\infty}\varepsilon^{i}(\Lambda)A^{i}$,
where the coefficients $\kappa^{i},\,\varepsilon^{i}$ can be determined by
iterations of the TBA equations. Here, we present $\kappa(k)$ in terms of
$\delta=\pi-Q$

$\displaystyle\int_{0}^{\pi-Q}\mathrm{d}k\cos
k(a_{1}(\sin{k}+\Lambda)+a_{1}(\sin{k}-\Lambda))\kappa(\pi-k)$ (67)
$\displaystyle=$
$\displaystyle\delta\,2a_{1}(\Lambda)\kappa^{0}(\pi)+\delta^{2}\,(2a_{1}(\Lambda)\kappa^{1}(\pi)-a_{1}(\Lambda)\kappa^{0^{\prime}}(\pi))$
$\displaystyle+\delta^{3}\,(2a_{1}(\Lambda)\kappa^{2}(\pi)-a_{1}(\Lambda)\kappa^{1^{\prime}}(\pi)+\frac{1}{3}(a_{1}^{{}^{\prime\prime}}(\Lambda)-a_{1}(\Lambda))\kappa^{0}(\pi)+\frac{1}{3}a_{1}(\Lambda)\kappa^{0^{\prime\prime}}(\pi)),$

where the prime ′ (′′) denotes first (second) derivative with respect to $k$.
In order to derive above expansion, the following derivative formula are used

$\displaystyle F(n)$ $\displaystyle=$
$\displaystyle\int_{0}^{a(n)}\mathrm{d}xf(x,n),$ $\displaystyle
F^{{}^{\prime}}(n)$ $\displaystyle=$
$\displaystyle\int_{0}^{a(n)}\mathrm{d}x\frac{\partial{f(x,n)}}{\partial{n}}+f(a(n),n)\frac{\mathrm{d}a(n)}{\mathrm{d}n},$
(68)

here $n$ stands for a small quantity. Using the above expressions, we may
calculate (66)

$\varepsilon_{1}(\Lambda)=\varepsilon^{0}(\Lambda)-\int_{-A}^{A}\mathrm{d}\Lambda^{\prime}a_{2}\left(\Lambda-\Lambda^{\prime}\right)\varepsilon_{1}(\Lambda^{\prime}),$
(69)

where the leading term $\varepsilon^{0}(\Lambda)$ incorporating the role of
charge is modified as:

$\displaystyle\varepsilon^{0}(\Lambda)$ $\displaystyle=$ $\displaystyle
2B-4\operatorname{Re}\sqrt{1-(\Lambda-\mathrm{i}u)^{2}}+4u$
$\displaystyle+\delta(2a_{1}(\Lambda)\kappa^{0}(\pi))$
$\displaystyle+\delta^{2}(2a_{1}(\Lambda)\kappa^{1}(\pi)-a_{1}(\Lambda)\kappa^{0^{\prime}}(\pi))$
$\displaystyle+\delta^{3}(2a_{1}(\Lambda)\kappa^{2}(\pi)-a_{1}(\Lambda)\kappa^{1^{\prime}}(\pi)+\frac{1}{3}(a_{1}^{{}^{\prime\prime}}(\Lambda)-a_{1}(\Lambda))\kappa^{0}(\pi)+\frac{1}{3}a_{1}(\Lambda)\kappa^{0^{\prime\prime}}(\pi)).$

Taking an expansion on equation (69) and comparing order by order with respect
to the small $A$ on both sides, we obtain the coefficients of the first three
orders in $\varepsilon_{1}(\Lambda)$

$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$
$\displaystyle\varepsilon^{0}(\Lambda)+A(-2a_{2}(\Lambda)\varepsilon^{0}(0))$
$\displaystyle+A^{2}(-2a_{2}(\Lambda)\varepsilon^{1}(0)-a_{2}(\Lambda)\varepsilon^{0^{\prime}}(0))$
$\displaystyle+A^{3}(-a_{2}(\Lambda)\varepsilon^{1^{\prime}}(0)-2a_{2}(\Lambda)\varepsilon^{2}(0)-\frac{1}{3}(a_{2}^{{}^{\prime\prime}}(\Lambda)\varepsilon^{0}(0)+a_{2}(\Lambda)\varepsilon^{0^{\prime\prime}}(0))),$

with

$\displaystyle\varepsilon^{0}(0)$ $\displaystyle=$ $\displaystyle
2B+4(u-\sqrt{1+u^{2}})+\frac{2}{\pi u}\kappa^{0}(\pi)\delta+\frac{1}{\pi
u}(2\kappa^{1}(\pi)-\kappa^{0^{\prime}}(\pi))\delta^{2},$ (72)
$\displaystyle\varepsilon^{0^{\prime}}(0)$ $\displaystyle=$ $\displaystyle
0,\quad\varepsilon^{0^{\prime\prime}}(0)=\frac{4}{(1+u^{2})^{\frac{3}{2}}}.$
(73)

The quantities associated with the charge degrees of freedom can be calculated
using the same manner. Consequently, expressions of charge and spin equations
can be written in terms of the orders of $\delta$ and $A$

$\displaystyle\kappa(k)$ $\displaystyle=$
$\displaystyle\bar{\delta}\bm{M_{c}}\bar{A}^{T},$ (74)
$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$
$\displaystyle\bar{\delta}\bm{M_{s}}\bar{A}^{T},$ (75)

where
$\bar{\delta}=\left(\delta^{0},\delta^{1},\delta^{2},\delta^{3}\right),\bar{A}=\left(A^{0},A^{1},A^{2},A^{3}\right)$,
$T$ denotes transpose and $\bm{M_{c}},\bm{M_{s}}$ are two matrices which take
the following forms

$\displaystyle\bm{M_{c}}=$ $\displaystyle\left(\begin{matrix}-2\cos
k-\mu-2u-B&0&0&0\\\ 2a_{1}(\sin{k})\beta_{2}&\frac{4}{\pi
u}a_{1}(\sin{k})\beta_{1}&0&0\\\ -\frac{2}{\pi
u}a_{1}(\sin{k})\beta_{2}&-\frac{4}{\pi^{2}u^{2}}a_{1}(\sin{k})\beta_{1}+\frac{8}{\pi^{2}u^{2}}a_{1}(\sin{k})\beta_{2}&0&0\\\
\left(\begin{matrix}\frac{2}{\pi^{2}u^{2}}a_{1}(\sin{k})\beta_{2}\\\
+\frac{1}{3}a_{1}^{{}^{\prime\prime}}(\sin{k})\beta_{2}\\\
+\frac{4}{3(1+u^{2})^{\frac{3}{2}}}a_{1}(\sin{k})\end{matrix}\right)&0&0&0\end{matrix}\right),$
(76) $\displaystyle\bm{M_{s}}=$
$\displaystyle\begin{pmatrix}\left(\begin{matrix}2B+4u\\\
-4\operatorname{Re}\sqrt{1-(\Lambda-\mathrm{i}u)^{2}}\end{matrix}\right)&2a_{1}(\Lambda)\beta_{1}&0&\left(\begin{matrix}\frac{1}{3}(a_{1}^{{}^{\prime\prime}}(\Lambda)\\\
-a_{1}(\Lambda))\beta_{1}\\\ -\frac{2}{3}a_{1}(\Lambda)\end{matrix}\right)\\\
-2a_{2}(\Lambda)\beta_{2}&\left(\begin{matrix}\frac{4}{\pi
u}(a_{1}(\Lambda)\beta_{2}\\\
-a_{2}(\Lambda)\beta_{1})\end{matrix}\right)&\frac{8}{\pi^{2}u^{2}}a_{1}(\Lambda)\beta_{1}&0\\\
\frac{2}{\pi
u}a_{2}(\Lambda)\beta_{2}&\left(\begin{matrix}-\frac{4}{\pi^{2}u^{2}}a_{1}(\Lambda)\beta_{2}\\\
-\frac{8}{\pi^{2}u^{2}}a_{2}(\Lambda)\beta_{2}\\\
+\frac{4}{\pi^{2}u^{2}}a_{2}(\Lambda)\beta_{1}\end{matrix}\right)&0&0\\\
\left(\begin{matrix}-\frac{2}{\pi^{2}u^{2}}a_{2}(\Lambda)\beta_{2}\\\
-\frac{1}{3}a_{2}^{{}^{\prime\prime}}(\Lambda)\beta_{2}\\\
-\frac{4}{3(1+u^{2})^{\frac{3}{2}}}a_{2}(\Lambda)\end{matrix}\right)&0&0&0\end{pmatrix}.$
(77)

To simplify our notations in the following analysis, we define

$\beta_{1}=2-\mu-2u-B,\qquad\beta_{2}=2B+4(u-\sqrt{1+u^{2}}),$ (78)

which characterize the leading contributions in the charge and spin dressed
energies, respectively. By analysing the construction of the matrices
$\bm{M_{c}},\bm{M_{s}}$, it can be deduced that $\beta_{1}$ and $\beta_{2}$
are at least the second orders of $A$ or $\delta$, namely we have

$\displaystyle\kappa(k)$ $\displaystyle=$ $\displaystyle-2\cos
k-\mu-2u-B+2a_{1}(\sin{k})\beta_{2}A+\frac{4}{3(1+u^{2})^{\frac{3}{2}}}a_{1}(\sin{k})A^{3},$
(79) $\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$ $\displaystyle
2B+4(u-\operatorname{Re}\sqrt{1-(\Lambda-\mathrm{i}u)^{2}})$ (80)
$\displaystyle+2a_{1}(\Lambda)\beta_{1}\delta-\frac{2}{3}a_{1}(\Lambda)\delta^{3}-2a_{2}(\Lambda)\beta_{2}A-\frac{4}{3(1+u^{2})^{\frac{3}{2}}}a_{2}(\Lambda)A^{3}.$

In the above results, we only retain terms up to the third order. By applying
the Fermi point conditions $\kappa(Q)=0$ and $\varepsilon_{1}(A)=0$, we obtain
$\beta_{1}$ and $\beta_{2}$ as

$\displaystyle\beta_{1}$ $\displaystyle=$
$\displaystyle\delta^{2}+\frac{4\eta_{1}}{3\pi u}A^{3},$ (81)
$\displaystyle\beta_{2}$ $\displaystyle=$ $\displaystyle-\frac{4}{3\pi
u}\delta^{3}-\eta_{1}\left(A^{2}+\frac{2}{3\pi u}A^{3}\right)$ (82)

with $\eta_{1}=\frac{2}{(1+u^{2})^{3/2}}$. Utilizing iteration repeatedly, the
Fermi point conditions give rise to (up to the order
$O(\beta^{2}_{1},\beta^{2}_{2})$)

$\displaystyle\delta^{2}$ $\displaystyle=$
$\displaystyle\beta_{1}-\frac{4}{3\pi
u\eta^{\frac{1}{2}}_{1}}(-\beta_{2})^{\frac{3}{2}}\left\\{1-\frac{1}{\pi
u}\left[(\frac{-\beta_{2}}{\eta_{1}})^{\frac{1}{2}}+\frac{2\beta^{\frac{3}{2}}_{1}}{(-\beta_{2})}\right]\right\\},$
(83) $\displaystyle\eta_{1}A^{2}$ $\displaystyle=$
$\displaystyle-\beta_{2}-\frac{2}{3\pi
u}\left[\frac{(-\beta_{2})^{\frac{3}{2}}}{\eta^{\frac{1}{2}}_{1}}+2\beta_{1}^{\frac{3}{2}}\right]+\frac{2}{3\pi^{2}u^{2}}\frac{(-\beta_{2})^{2}}{\eta_{1}}$
(84)
$\displaystyle+\frac{8}{3\pi^{2}u^{2}\eta^{\frac{1}{2}}_{1}}\left[\frac{1}{2}(-\beta_{2})^{\frac{1}{2}}\beta_{1}^{\frac{3}{2}}+(-\beta_{2})^{\frac{3}{2}}\beta_{1}^{\frac{1}{2}}\right].$

Similarly, we apply the above technique to deal with the density Bethe ansatz
equations (20) and (21) and find the matrix form:

$\displaystyle\rho(k)$ $\displaystyle=$
$\displaystyle\bar{\delta}\bm{M_{dc}}\bar{A}^{T},$ (85)
$\displaystyle\sigma_{1}(\Lambda)$ $\displaystyle=$
$\displaystyle\bar{\delta}\bm{M_{ds}}\bar{A}^{T},$ (86)

where the matrices $\bm{M_{dc}},\bm{M_{ds}}$ are given by

$\displaystyle\bm{M_{dc}}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\frac{1}{2\pi}&0&0&0\\\ 2\cos
ka_{1}(\sin{k})\lambda_{1}&-\frac{2}{\pi^{2}u}\cos ka_{1}(\sin{k})&0&0\\\
-\frac{2}{\pi u}\cos ka_{1}(\sin{k})\lambda_{1}&\frac{4}{\pi u}\cos
ka_{1}(\sin{k})\left(\frac{2}{\pi
u}\lambda_{1}+\frac{1}{2\pi^{2}u}\right)&0&0\\\
\left(\begin{matrix}\frac{2}{\pi^{2}u^{2}}\cos ka_{1}(\sin{k})\lambda_{1}\\\
+\frac{1}{3}a_{1}^{{}^{\prime\prime}}(\sin{k})\cos k\lambda_{1}\\\
+\frac{1}{3}\cos
ka_{1}(\sin{k})\lambda_{2}\end{matrix}\right)&0&0&0\end{pmatrix},$ (87)
$\displaystyle\bm{M_{ds}}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\frac{1}{\pi}\operatorname{Re}\frac{1}{\sqrt{1-(\Lambda+\mathrm{i}u)^{2}}}&-\frac{1}{\pi}a_{1}(\Lambda)&0&-\frac{1}{6\pi}a_{1}^{{}^{\prime\prime}}(\Lambda)\\\
-2a_{2}(\Lambda)\lambda_{1}&\frac{4}{\pi
u}(a_{1}(\Lambda)\lambda_{1}+\frac{1}{2\pi}a_{2}(\Lambda))&-\frac{4}{\pi^{3}u^{2}}a_{1}(\Lambda)&0\\\
\frac{2}{\pi
u}a_{2}(\Lambda)\lambda_{1}&\left(\begin{matrix}-\frac{4}{\pi^{2}u^{2}}a_{1}(\Lambda)\lambda_{1}\\\
-\frac{8}{\pi^{2}u^{2}}a_{2}(\Lambda)\lambda_{1}\\\
-\frac{2}{\pi^{3}u^{2}}a_{2}(\Lambda)\end{matrix}\right)&0&0\\\
\left(\begin{matrix}-\frac{2}{\pi^{2}u^{2}}a_{2}(\Lambda)\lambda_{1}\\\
-\frac{1}{3}a_{2}^{{}^{\prime\prime}}(\Lambda)\lambda_{1}\\\
-\frac{1}{3}a_{2}(\Lambda)\lambda_{2}\end{matrix}\right)&0&0&0\end{pmatrix},$
(88)

where the constant values $\lambda_{1}=\frac{1}{\pi\sqrt{1+u^{2}}}$,
$\lambda_{2}=\frac{1-2u^{2}}{\pi(1+u^{2})^{5/2}}$ related to the interaction
scale. By integrating the densities within the Fermi points (22), (23), the
particle number $n_{c}$ and down-spin number $n_{\downarrow}$ are given by

$\displaystyle n_{c}$ $\displaystyle=$ $\displaystyle
1-2\int_{0}^{\pi-Q}\mathrm{d}k\rho(\pi-k)=1-2\delta\left(\frac{1}{2\pi}-\frac{2\lambda_{1}}{\pi
u}A+\frac{2\lambda_{1}}{\pi^{2}u^{2}}A^{2}\right)-\frac{4}{\pi^{3}u^{2}}\delta^{2}A,$
(89) $\displaystyle n_{\downarrow}$ $\displaystyle=$ $\displaystyle
2\left[A\left(\lambda_{1}-\frac{2}{\pi
u}\delta\left(\frac{1}{2\pi}-\frac{2\lambda_{1}}{\pi
u}A\right)\right)-A^{2}\frac{1}{\pi
u}\left(\lambda_{1}-\frac{1}{\pi^{2}u}\delta\right)+A^{3}\left(\frac{\lambda_{1}}{\pi^{2}u^{2}}+\frac{\lambda_{2}}{6}\right)\right].$
(90)

It is more intuitive to express the Fermi points in terms of particle numbers.
Based on the above relations, we have the small parameters near the Fermi
points $\delta=\pi-Q$ and $A$ in terms of doping parameter
$\hat{n}_{c}=1-n_{c}$ and $n_{\downarrow}$,

$\displaystyle\delta$ $\displaystyle=$
$\displaystyle\pi\hat{n}_{c}\left(1+\frac{2}{u}n_{\downarrow}+\frac{4}{u^{2}}n_{\downarrow}^{2}\right),$
(91) $\displaystyle A$ $\displaystyle=$ $\displaystyle\frac{\pi
u}{2}\zeta_{1}n_{\downarrow}\left[1+\zeta_{1}\left(\hat{n}_{c}+\frac{n_{\downarrow}}{2}\right)+\zeta^{2}_{1}\left(\hat{n}_{c}+\frac{n_{\downarrow}}{2}\right)^{2}-\zeta_{2}n_{\downarrow}^{2}\right]$
(92)

with $\zeta_{1}=\frac{1}{\pi u\lambda_{1}}$, and
$\zeta_{2}=\frac{\lambda_{2}}{24\lambda^{3}_{1}}$. Using the free energy (17)
with $T=0$ and the above expression (79) for $\kappa(k)$, we may obtain

$f=-\mu-
u-B-\frac{8\lambda_{1}}{3(1+u^{2})^{\frac{3}{2}}}A^{3}-\frac{2}{3\pi}\delta^{3}.$
(93)

On the other hand, in terms of the discrete symmetries of the repulsive
Hubbard model essler2005one , the free energy of the repulsive system can be
transformed to that of the attractive case via the relation
$f_{r}(\mu,B,T,u)=f_{a}(-B,-\mu,T,-u)-\mu-B$, where subscript $a(r)$ means
attractive (repulsive) interaction. It turns out that the results obtained not
only cover that of the attractive model Cheng:2018A  in the low density and
strong coupling limits, but also are valid for the arbitrary interaction
strength away from the strong coupling limit. In general, it is
straightforward to show that there is a one-to-one correspondence from the
high density area for repulsive case to low density regime for attractive case
Cheng:2018A . Take the derivative of free energy, one can obtain the
relationship between their thermodynamic quantities:

$\displaystyle\frac{\partial f_{r}}{\partial\mu}$ $\displaystyle=$
$\displaystyle-\frac{\partial f_{a}}{\partial B}-1\longrightarrow
n_{r,\mu}=1-2m_{a,B},$ (94) $\displaystyle\frac{\partial f_{r}}{\partial B}$
$\displaystyle=$ $\displaystyle-\frac{\partial
f_{a}}{\partial\mu}-1\longrightarrow 2m_{r,B}=1-n_{a,\mu}.$ (95)

Using these results, and by the definitions of spin and charge characteristic
velocities $v_{c}=\frac{\kappa^{{}^{\prime}}(k)}{2\pi\rho(k)}|_{k=Q}$,
$v_{s}=\frac{\varepsilon_{1}^{{}^{\prime}}(\Lambda)}{2\pi\sigma(\Lambda)}|_{\Lambda=A}$,
we obtain the following quantities

$\displaystyle\kappa^{{}^{\prime}}(Q)$ $\displaystyle=$ $\displaystyle
2\delta-\frac{1}{3}\delta^{3},$ (96) $\displaystyle\rho(Q)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}-\frac{2\lambda_{1}}{\pi
u}A+\frac{2\lambda_{1}}{\pi^{2}u^{2}}A^{2}-\frac{2\lambda_{1}}{\pi^{3}u^{3}}A^{3}+\frac{2\lambda_{1}}{\pi
u^{3}}\delta^{2}A+\frac{2\lambda_{1}}{3\pi u^{3}}A^{3}$ (97)
$\displaystyle-\frac{\lambda_{2}}{3\pi u}A^{3}+\frac{\lambda_{1}}{\pi
u}\delta^{2}A+\frac{2}{\pi^{3}u^{2}}\delta
A-\frac{2}{\pi^{4}u^{3}}\left(1+4\pi\lambda_{1}\right)\delta A^{2},$
$\displaystyle\varepsilon_{1}^{{}^{\prime}}(A)$ $\displaystyle=$
$\displaystyle 2\eta_{1}A+\frac{\eta_{2}}{6}A^{3},$ (98)
$\displaystyle\sigma_{1}(A)$ $\displaystyle=$
$\displaystyle\lambda_{1}+\frac{\lambda_{2}}{2}A^{2}-\frac{1}{\pi^{2}u}\delta+\frac{4\lambda_{1}}{\pi^{2}u^{2}}\delta
A-\frac{8\lambda_{1}}{\pi^{3}u^{3}}\delta A^{2}$ (99)
$\displaystyle+\frac{1}{\pi^{2}u^{3}}\delta
A^{2}-\frac{4}{\pi^{4}u^{3}}\delta^{2}A+\frac{1}{3\pi^{2}u^{3}}\delta^{3}-\frac{\lambda_{1}}{\pi
u}A+\frac{1}{\pi^{3}u^{2}}\delta A$ $\displaystyle+\frac{\lambda_{1}}{3\pi
u^{3}}A^{3}+\frac{\lambda_{1}}{\pi^{2}u^{2}}A^{2}-\frac{1}{\pi^{4}u^{3}}\delta
A^{2}-\frac{\lambda_{1}}{\pi^{3}u^{3}}A^{3}-\frac{\lambda_{2}}{6\pi u}A^{3}$

with $\eta_{2}=12\frac{1-4u^{2}}{\pi(1+u^{2})^{7/2}}$. The analytic results
and numerical simulation of sound velocities and thermodynamic properties by
the solution of TBA equations are showed in figure 10. From the figure 10 (a),
we observe that $v_{c}$ approaches zero at the critical point for the phase
transition from magnetized phase IV to the Mott phase V. On the other hand,
$v_{s}$ goes to zero that signifies an approach to the spin fully polarized
phase II. Thus the velocities can serve as a signature of quantum phase
trasnition. In figure 10 (b), we show the specific heat in terms of velocities
at various temperatures: demonstrating the analytical results in comparing
with numerical calculations. The analytical results of the specific heat are
given by

$\displaystyle C_{v}$ $\displaystyle=$ $\displaystyle\frac{\pi
T}{3}\left(\frac{1}{v_{c}}+\frac{1}{v_{s}}\right)\quad\text{phase IV},$ (100)
$\displaystyle C_{v}$ $\displaystyle=$ $\displaystyle\frac{\pi
T}{3}\frac{1}{v_{c}}\quad\text{phase II},$ (101) $\displaystyle C_{v}$
$\displaystyle=$ $\displaystyle\frac{\pi T}{3}\frac{1}{v_{s}}\quad\text{phase
V}.$ (102)

The quantity $C_{v}/T$ becomes temperature independent for the TLL states.

Figure 10: (a) Spin and charge velocities as function of chemical potential in
phase II, IV and V, in which analytic results are given by the equations
(96)-(99) for IV, equation (122) for II and equation (127) for V,
respectively. (b) shows the analytical and numerical results of the specific
heat in terms of chemical potential. The curves of $C_{v}/T$ at different
temperatures collapse to a single curve of the analytical result in each
phase. (c) shows the magnetic susceptibility: agreement between the numerical
and analytical results which are given by equations (104), (122), (126). In
the TLL phases, the susceptibility does not depend on temperature. The whole
numerical setting reads $B=0.825,\,u=1$.

###  III.4 III.4 Additivity rules of charge and spin susceptibilities

Now we calculate the charge and spin susceptibilities $\chi_{c}=\frac{\partial
n_{c}}{\partial\mu}|_{B(m)}$ and $\chi_{s}=\frac{\partial m}{\partial
B}|_{\mu(n_{c})}$ with $m=\frac{n_{c}}{2}-n_{\downarrow}$. In the grand
canonical ensemble driven by external potentials $B$ and $\mu$, we have the
following expressions of spin and charge susceptibilities

$\displaystyle\chi_{c}$ $\displaystyle=$ $\displaystyle\frac{\partial
n}{\partial\mu}|_{B}=\left(-\frac{1}{\pi}+\frac{4\lambda_{1}}{\pi
u}A-\frac{4\lambda_{1}}{\pi^{2}u^{2}}A^{2}-\frac{8}{\pi^{3}u^{2}}\delta
A\right)\frac{\partial\delta}{\partial\mu}$ (103)
$\displaystyle+\left(\frac{4\lambda_{1}}{\pi
u}\delta-\frac{8\lambda_{1}}{\pi^{2}u^{2}}\delta
A-\frac{4}{\pi^{3}u^{2}}\delta^{2}\right)\frac{\partial
A}{\partial\mu}=\chi_{c}^{(1)}+\chi_{c}^{(2)},$ $\displaystyle\chi_{s}$
$\displaystyle=$ $\displaystyle\frac{\partial m}{\partial
B}|_{\mu}=\left[-\frac{1}{2\pi}-\frac{4}{\pi^{3}u^{2}}\delta
A+\frac{2}{\pi^{2}u}\left(1+\pi\lambda_{1}\right)A-\frac{2}{\pi^{3}u^{2}}\left(1+5\pi\lambda_{1}\right)A^{2}\right]\frac{\partial\delta}{\partial
B}$ (104) $\displaystyle+\left[-2\lambda_{1}\left(1-\frac{2}{\pi
u}A+\frac{3}{\pi^{2}u^{2}}A^{2}\right)-\frac{2}{\pi^{3}u^{2}}\delta^{2}\right.$
$\displaystyle+\left.\frac{2}{\pi^{2}u}\left(1+\pi\lambda_{1}\right)\delta-\frac{4}{\pi^{3}u^{2}}\left(1+5\pi\lambda_{1}\right)\delta
A-\lambda_{2}A^{2}\right]\frac{\partial A}{\partial B}$ $\displaystyle=$
$\displaystyle\chi_{s}^{(1)}+\chi_{s}^{(2)}.$

We observe that the charge susceptibilities $\chi_{c}^{(1)},\,\chi_{c}^{(2)}$
(the spin susceptibility $\chi_{s}^{(1)},\,\chi_{s}^{(2)}$) denote the
contributions from different resources $\delta$ and $A$, respectively, see
(103) and (104). The susceptibilities can be split into two parts which are
reminiscent of the additivity rules found in the FFLO state of attractive
situation Cheng:2018A . It should be noted that the decomposition terms stem
from two sources, i.e., the changes of charge and spin with respect to the
chemical potential and magnetic field. In the (103) and (104),
$\frac{\partial\delta}{\partial\mu}$, $\frac{\partial A}{\partial\mu}$,
$\frac{\partial\delta}{\partial B}$ and $\frac{\partial A}{\partial B}$ can be
determined by taking the derivative of both sides of (81) and (82) (note:
$\beta_{1},\,\beta_{2}$ depend on $\mu$ and $B$):

$\displaystyle\frac{\partial A}{\partial\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{\eta_{1}\left(A+\frac{A^{2}}{\pi u}\right)\frac{\pi
u}{\delta}-\frac{4\eta_{1}}{\pi u}A^{2}},$ (105)
$\displaystyle\frac{\partial\delta}{\partial\mu}$ $\displaystyle=$
$\displaystyle-\frac{\left(A+\frac{A^{2}}{\pi
u}\right)}{2\delta\left(A+\frac{A^{2}}{\pi
u}\right)-\frac{8}{\pi^{2}u^{2}}\delta^{2}A^{2}},$ (106)
$\displaystyle\frac{\partial A}{\partial B}$ $\displaystyle=$
$\displaystyle-\frac{1-\frac{1}{\pi u}\delta}{\eta_{1}\left(A+\frac{A^{2}}{\pi
u}\right)-\frac{4\eta_{1}}{\pi^{2}u^{2}}\delta A^{2}},$ (107)
$\displaystyle\frac{\partial\delta}{\partial B}$ $\displaystyle=$
$\displaystyle-\frac{1-\frac{3A}{\pi u}}{2\delta\left(1+\frac{A}{\pi
u}\right)-\frac{8}{\pi^{2}u^{2}}\delta^{2}A}.$ (108)

We plot the magnetic susceptibility in figure (10) (c) and observe that the
susceptibility does not depend on temperature in the TLL regime, but becomes
divergent in the vicinities of the critical points.

In the canonical ensemble with fixed density and magnetization, we take the
total derivative of (91) and (92) with respect to $n_{c}$ and $m$ under the
condition $\mathrm{d}n_{\downarrow}=\mathrm{d}n_{c}/2-\mathrm{d}m$

$\displaystyle\mathrm{d}\delta$ $\displaystyle=$
$\displaystyle\left[\pi\left(1+\frac{2}{u}n_{\downarrow}+\frac{4}{u^{2}}n_{\downarrow}^{2}\right)-\pi\hat{n}_{c}\left(\frac{1}{u}+\frac{4}{u^{2}}n_{\downarrow}\right)\right]\mathrm{d}\hat{n}_{c}-\pi\hat{n}_{c}\left(\frac{2}{u}+\frac{8}{u^{2}}n_{\downarrow}\right)\mathrm{d}m,$
(109) $\displaystyle\mathrm{d}A$ $\displaystyle=$
$\displaystyle\mathrm{d}\hat{n}_{c}\left\\{-\frac{\pi
u\zeta_{1}}{4}\left[1+\zeta_{1}\left(\hat{n}_{c}+\frac{n_{\downarrow}}{2}\right)+\zeta^{2}_{1}\left(\hat{n}_{c}+\frac{n_{\downarrow}}{2}\right)^{2}-\zeta_{2}n_{\downarrow}^{2}\right]\right.$
(110) $\displaystyle\left.+\frac{\pi
u\zeta_{1}}{2}n_{\downarrow}\left[\frac{3\zeta_{1}}{4}+\frac{3\zeta^{2}_{1}}{2}\left(\hat{n}_{c}+\frac{n_{\downarrow}}{2}\right)+\zeta_{2}n_{\downarrow}\right]\right\\}$
$\displaystyle+\mathrm{d}m\left\\{-\frac{\pi
u\zeta_{1}}{2}\left[1+\zeta_{1}\left(\hat{n}_{c}+\frac{n_{\downarrow}}{2}\right)+\zeta^{2}_{1}\left(\hat{n}_{c}+\frac{n_{\downarrow}}{2}\right)^{2}-\zeta_{2}n_{\downarrow}^{2}\right]\right.$
$\displaystyle\left.+\frac{\pi
u\zeta_{1}}{2}n_{\downarrow}\left[-\frac{\zeta_{1}}{2}-\zeta^{2}_{1}\left(\hat{n}_{c}+\frac{n_{\downarrow}}{2}\right)+2\zeta_{2}n_{\downarrow}\right]\right\\}.$

Here $\frac{\partial\mu}{\partial n_{c}}|_{m},\frac{\partial B}{\partial
m}|_{n_{c}}$ can be similarly extracted from (81) (82). Similar to the
(105)-(108) with the difference in the quantity being differentiated, we have

$\displaystyle-\frac{\partial\mu}{\partial\hat{n}_{c}(m)}|_{m(\hat{n}_{c})}-\frac{\partial
B}{\partial\hat{n}_{c}(m)}|_{m(\hat{n}_{c})}=2\delta\frac{\partial\delta}{\partial\hat{n}_{c}(m)}|_{m(\hat{n}_{c})}+\frac{4\eta_{1}}{\pi
u}A^{2}\frac{\partial A}{\partial\hat{n}_{c}(m)}|_{m(\hat{n}_{c})},$ (111)
$\displaystyle\frac{\partial
B}{\partial\hat{n}_{c}(m)}|_{m(\hat{n}_{c})}=-\frac{2}{\pi
u}\delta^{2}\frac{\partial\delta}{\partial\hat{n}_{c}(m)}|_{m(\hat{n}_{c})}-\eta_{1}\left(A+\frac{A^{2}}{\pi
u}\right)\frac{\partial A}{\partial\hat{n}_{c}(m)}|_{m(\hat{n}_{c})}.$ (112)

Therefore, by solving this set of equations together with total differential
expressions (109) and (110), we obtain the charge and spin susceptibilities in
terms of $n_{c},\,m$ in the canonical ensemble

$\displaystyle\frac{1}{\chi_{c}}$ $\displaystyle=$
$\displaystyle\frac{\partial\mu}{\partial
n_{c}}|_{m}=2\delta\left(1-\frac{\delta}{\pi
u}\right)\frac{\partial\delta}{\partial\hat{n}_{c}}|_{m}-\eta_{1}A\left(1-\frac{3A}{\pi
u}\right)\frac{\partial A}{\partial\hat{n}_{c}}|_{m}$ (113) $\displaystyle=$
$\displaystyle\frac{1}{\chi_{c}^{(1)}}+\frac{1}{\chi_{c}^{(2)}},$
$\displaystyle\frac{1}{\chi_{s}}$ $\displaystyle=$
$\displaystyle\frac{\partial B}{\partial m}|_{\hat{n}_{c}}=-\frac{2}{\pi
u}\delta^{2}\frac{\partial\delta}{\partial
m}|_{\hat{n}_{c}}-\eta_{1}A\left(1+\frac{A}{\pi u}\right)\frac{\partial
A}{\partial m}|_{\hat{n}_{c}}$ (114) $\displaystyle=$
$\displaystyle\frac{1}{\chi_{s}^{(1)}}+\frac{1}{\chi_{s}^{(2)}}.$

In contrast to the the case for fixed external potentials, these formulas
expose reciprocal additivity relation at a fixed magnetization and density,
see (113) and (114). The two ensembles are related to each other by the
Jacobian determinant evaluated by total differential of $\mu,\,B$ with respect
to $n,\,m$. Apart from the numerical arithmetic and approximation method as we
have conducted above, we would like to discuss dressed charge matrix
$\mathrm{Z}$frahm1990critical ; frahm1991correlation ; essler2005one

$\mathrm{Z}=\left(\begin{matrix}\xi_{cc}(Q)&\xi_{cs}(A)\\\
\xi_{sc}(Q)&\xi_{ss}(A)\end{matrix}\right),$ (115)

whose elements are determined by

$\xi_{ab}\left(x_{b}\right)=\delta_{ab}+\sum_{d}\int_{-X_{d}}^{X_{d}}\mathrm{d}x_{d}\xi_{ad}\left(x_{d}\right)K_{db}\left(x_{d},x_{b}\right).$
(116)

Here the kernels are given by

$\displaystyle K_{cc}\left(x_{c},y_{c}\right)$ $\displaystyle=$ $\displaystyle
0,$ $\displaystyle K_{sc}\left(x_{c},x_{s}\right)$ $\displaystyle=$
$\displaystyle a_{1}\left(\sin\left(x_{c}\right)-x_{s}\right),$ $\displaystyle
K_{cs}\left(x_{c},x_{s}\right)$ $\displaystyle=$
$\displaystyle\cos\left(x_{c}\right)a_{1}\left(\sin\left(x_{c}\right)-x_{s}\right),$
$\displaystyle K_{ss}\left(x_{s},y_{s}\right)$ $\displaystyle=$
$\displaystyle-a_{2}\left(x_{s}-y_{s}\right).$ (117)

These dressed charges can be used to calculate rigorous solutions for
susceptibilities and conformal dimensions in asymptotics of correlation
functions frahm1990critical ; frahm1991correlation . We will discuss about
this study later.

One advantage of this method is that the iteration process is greatly
simplified both analytically and numerically, and the final result is clear
and unambiguous. In the grand canonical ensemble, some exact relations are
relevant to the dressed charge matrix (115), for example essler2005one ,

$\displaystyle\chi_{c}|_{B}$ $\displaystyle=$ $\displaystyle\frac{\partial
n_{c}}{\partial\mu}|_{B}=\frac{\mathrm{Z}^{2}_{cc}}{\pi
v_{c}}+\frac{\mathrm{Z}^{2}_{cs}}{\pi v_{s}}=\chi_{c}^{(1)}+\chi_{c}^{(2)},$
(118) $\displaystyle\chi_{s}|_{\mu}$ $\displaystyle=$
$\displaystyle\frac{\partial m}{\partial
B}|_{\mu}=\frac{(\mathrm{Z}_{cc}-2\mathrm{Z}_{sc})^{2}}{2\pi
v_{c}}+\frac{(\mathrm{Z}_{cs}-2\mathrm{Z}_{ss})^{2}}{2\pi
v_{s}}=\chi_{s}^{(1)}+\chi_{s}^{(2)}.$ (119)

Regarding the circumstance of variable changes, we can work out the Jacobian
determinant from the equations (118) and (119), namely,
$J=2(\mathrm{det}\mathrm{Z})^{2}/(\pi^{2}v_{c}v_{s})$. As a result, we obtain
111The second equality in the equation (6.79) of the book essler2005one
misses a factor of 2. for canonical ensemble

$\displaystyle\frac{1}{\chi_{c}}|_{m}$ $\displaystyle=$
$\displaystyle\frac{\partial\mu}{\partial
n_{c}}|_{m}=\frac{\pi}{4}\frac{v_{c}(\mathrm{Z}_{cs}-2\mathrm{Z}_{ss})^{2}+v_{s}(\mathrm{Z}_{cc}-2\mathrm{Z}_{sc})^{2}}{(\mathrm{det}\mathrm{Z})^{2}}=\frac{1}{\chi_{c}^{(1)}}+\frac{1}{\chi_{c}^{(2)}},$
(120) $\displaystyle\frac{1}{\chi_{s}}|_{n_{c}}$ $\displaystyle=$
$\displaystyle\frac{\partial B}{\partial
m}|_{n_{c}}=\frac{\pi}{2}\frac{v_{c}\mathrm{Z}_{cs}^{2}+v_{s}\mathrm{Z}_{cc}^{2}}{(\mathrm{det}\mathrm{Z})^{2}}=\frac{1}{\chi_{s}^{(1)}}+\frac{1}{\chi_{s}^{(2)}}.$
(121)

Comparing our analytic results (103) and (104) (or (113), (114)) with the
dressed charge formula (118) and (119) (or (120), (121)), we observe that the
derivative term with respect to $\delta$ in the (103) is equivalent to the
first term associated with the charge velocity $v_{c}$ in (118), whereas the
derivative term with respect to $A$ in (103) is equivalent to the second term
associated with the spin velocity $v_{s}$ in (118), respectively. Similar
correspondences can be found from (104) and (119) in the grand canonical
ensemble, as well as in these equations (113), (114) and (120), (121),
respectively.

We remark that the results we obtained can also be applied to the phase II and
the Mott phase V, i.e., omitting the terms related to $A$ in II or $\delta$ in
phase V, respectively. Explicitly, the corresponding thermodynamic quantities
and velocities are given by
phase II:

$n_{c}=1-\frac{\delta}{\pi},n_{\downarrow}=0,\chi_{c}=\frac{1}{2\pi\delta},\chi_{s}=\frac{1}{4\pi\delta},v_{c}=2\delta-\frac{1}{3}\delta^{3},$
(122)

phase V:

$\displaystyle n_{c}$ $\displaystyle=$ $\displaystyle 1,$ (123) $\displaystyle
n_{\downarrow}$ $\displaystyle=$
$\displaystyle\int_{-A}^{A}\mathrm{d}\Lambda\sigma_{1}(\Lambda)=2\left[\lambda_{1}A-\frac{\lambda_{1}}{\pi
u}A^{2}+A^{3}\left(\frac{\lambda_{1}}{\pi^{2}u^{2}}+\frac{\lambda_{2}}{6}\right)\right],$
(124) $\displaystyle\chi_{c}$ $\displaystyle=$ $\displaystyle 0,$ (125)
$\displaystyle\chi_{s}$ $\displaystyle=$
$\displaystyle\left[2\lambda_{1}\left(1-\frac{2}{\pi
u}A+\frac{3}{\pi^{2}u^{2}}A^{2}\right)+\lambda_{2}A^{2}\right]\frac{1}{\eta_{1}A\left(1+\frac{A}{\pi
u}\right)},$ (126) $\displaystyle v_{s}$ $\displaystyle=$
$\displaystyle\frac{\eta_{1}A+\frac{\eta_{2}}{12}A^{3}}{\pi\left(\lambda_{1}+\frac{\lambda_{2}}{2}A^{2}-\frac{\lambda_{1}}{\pi
u}A+\frac{\lambda_{1}}{3\pi
u^{3}}A^{3}+\frac{\lambda_{1}}{\pi^{2}u^{2}}A^{2}-\frac{\lambda_{1}}{\pi^{3}u^{3}}A^{3}-\frac{\lambda_{2}}{6\pi
u}A^{3}\right)}.$ (127)

Finally, we make direct comparison between our analytical results and
numerical results from the TBA equations. In figure 11, the first row shows
density (a), compressibility (b) in the grand canonical ensemble and (c) in
the canonical ensemble. Whereas the second row shows the magnetization (d),
susceptibility (e) in the grand canonical ensemble and (f) in the canonical
ensemble. The density and magnetization in figure 11 (a)(d) cross quantum
phases from II, to IV and V as chemical potential varies. In figure 11 (b),
(c), (e), (f), we observe additivity rules in both the grand canonical and the
canonical ensemble. All analytic results agree well with the corresponding
numerical results in this figure.

Figure 11:  Charge density (a) and (d) magnetization vs chemical potential in
phases II, IV and V: agreement between analytic results given by (122) for
phase II, (89) (90) for phase IV, (123)(124) for phase V and the numerical
calculations from the TBA equations. The compressibility (b) and
susceptibility (e) in grand canonical ensemble show agreement between the
analytic results given by (103) (104) and numerical simulations based on (118)
(119). While compressibility (c) and susceptibility (f) in canonical ensemble
show agreement between the analytic results given by 113 (114) and numerical
simulations based on (120) (121). All figures are set to $B=0.825,\,u=1$.

###  III.5 III.5 Quantum criticality near quadruple critical point

The results obtained in the previous subsections present an elementary
understanding of the ground state properties and the behaviour of the TLLs.
Although the criticality induced by the variation of some external potentials
such as magnetic field and chemical potential has been well studied in
literature, this is not the case for the interaction-driven quantum critical
behaviour, even though interaction plays a central role in many-body systems.
In this subsection, we study quantum phase transitions and universal scaling
functions of properties of the 1D Hubbard model in terms of external fields.
And in the last section we focus on the interaction driven quantum
transitions. From the phase diagram figure (8), it can be observed that a
phase transition occurs in the repulsive Hubbard model when a degree of
freedom appears, disappears or become gapped. Although phase transition occurs
at zero temperature, thermal and quantum fluctuations can reach the same level
of the energy scale in the V-shaped critical regions at finite temperatures. A
natural question is whether quantum critical region has its own set of
universal laws similar to the TLL. The obvious answer, of course, is that the
quantum criticality for the same universality class of models is insensitive
to the microscopic details of the systems and shares general and universal
critical phenomena characterised by the critical exponents sachdev2011quantum
; hazzard2011techniques . For example, the universal scaling laws of
magnetization and susceptibility satisfy sachdev2011quantum ;
hazzard2011techniques

$\displaystyle m$ $\displaystyle=$ $\displaystyle
m_{0}+T^{d/z+1-1/(vz)}O_{1}\left[\frac{\mu-\mu_{c}}{T^{1/vz}},\frac{B-B_{c}}{T^{1/vz}},\frac{u-u_{c}}{T^{1/vz}}\right],$
(128) $\displaystyle\chi_{s}$ $\displaystyle=$
$\displaystyle\chi_{s0}+T^{d/z+1-2/(vz)}O_{2}\left[\frac{\mu-\mu_{c}}{T^{1/vz}},\frac{B-B_{c}}{T^{1/vz}},\frac{u-u_{c}}{T^{1/vz}}\right],$
(129)

respectively. Where $m,\,\chi_{s}$ represent the first and second orders of
thermodynamic quantities respectively, $\mu_{c},\,B_{c},\,u_{c}$ denotes the
critical fields at the critical point and $z,\,v$ stand for the dynamical and
the correlation critical exponent, respectively. The scaling functions (128)
and (129) consist of two parts: the first term denotes the background
contributions stemming from the unchanged degrees of freedom, and the second
term accounts for the singular part stemming from the sudden change of the
density of state of one degree of freedom. Another salient feature of critical
region is that the characteristic length scale diverges, providing us with a
feasible opportunity to capture system information by fractional exclusive
statistics (FES) sachdev2011quantum ; fukui1995haldane ;
vitoriano2009fractional ; zhang2022interaction .

In what follows, we embark on the analytical derivation of the singular
behaviours of thermodynamical properties involved solely charge quasimomenta
$k$ and length-$1$ $\Lambda$ strings. Without losing generality, we first
consider the dressed energy equations for the magnetized phase IV

$\displaystyle\kappa(k)$ $\displaystyle=$
$\displaystyle-2\cos{k}-\mu-2u-B-T\int_{-\infty}^{\infty}\mathrm{d}\Lambda
a_{1}(\sin{k}-\Lambda)\ln(1+e^{-\frac{\varepsilon_{1}(\Lambda)}{T}}),$ (130)
$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$ $\displaystyle
2B-T\int_{-\pi}^{\pi}\mathrm{d}k\cos
ka_{1}(\sin{k}-\Lambda)\ln(1+e^{-\frac{\kappa(k)}{T}})$ (131)
$\displaystyle+T\int_{-\infty}^{\infty}\mathrm{d}\Lambda^{\prime}a_{2}\left(\Lambda-\Lambda^{\prime}\right)\ln(1+e^{-\frac{\varepsilon_{1}(\Lambda^{\prime})}{T}}),$

where $a_{n}(x)=\frac{1}{2\pi}\frac{2nu}{\left(nu\right)^{2}+x^{2}}$. There
are five phase transitions in total, I-II, II-III, III-V, II-IV and IV-V,
among which the four cases can be treated uniformly around the quadruple
point.

I-II phase transition: To this end, we first study the phase transition from
an empty lattice phase I to partially filled II with the absence of down-spin.
In this transition, the phase transition boundary is simply expressed as
$2+2u+B_{c}+\mu_{c}=0$ with subscript $B_{c}$ and $\mu_{c}$ stand for the
critical fields. Thus the charge dressed energy is given by $\kappa(k)=-2\cos
k-\mu-2u-B$, showing free fermions on a 1D lattice. This gives the free energy
$f=u+\frac{T^{3/2}}{2\pi^{1/2}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{\Delta
B+\Delta\mu}{T}}\right)$. Using fundamental thermodynamic relations, we
directly obtain the first and second order thermodynamic quantities

$\displaystyle m$ $\displaystyle=$
$\displaystyle\frac{n}{2}=-\frac{T^{\frac{1}{2}}}{4\pi^{\frac{1}{2}}}f_{\frac{1}{2}},$
(132) $\displaystyle\chi_{s}$ $\displaystyle=$
$\displaystyle\frac{\chi_{c}}{2}=-\frac{T^{-\frac{1}{2}}}{4\pi^{\frac{1}{2}}}f_{-\frac{1}{2}},$
(133) $\displaystyle\frac{c_{v}}{T}$ $\displaystyle=$
$\displaystyle-\frac{3T^{-\frac{1}{2}}}{8\pi^{\frac{1}{2}}}f_{\frac{3}{2}}+\frac{T^{-\frac{3}{2}}\left(\Delta
B+\Delta\mu\right)}{2\pi^{\frac{1}{2}}}f_{\frac{1}{2}}-\frac{T^{-\frac{5}{2}}\left(\Delta
B+\Delta\mu\right)^{2}}{2\pi^{\frac{1}{2}}}f_{-\frac{1}{2}}$ (134)

with $f_{n}=\mathrm{Li}_{n}\left(-\mathrm{e}^{\frac{\Delta
B+\Delta\mu}{T}}\right)$.

For the remaining four phase transitions, what we need to prepare before
obtaining the scaling forms is to express the dressed energy equations and the
free energy in terms of polylog functions. Under the assumption of
$n\rightarrow 1,\,n_{\downarrow}\rightarrow 0$, we expand the kernal function
$a_{n}$ around $k=\pi$ and $\Lambda=0$, the coupled equations become as

$\displaystyle\kappa(k)$ $\displaystyle=$
$\displaystyle-2\cos{k}+2\left(\frac{1}{\pi u^{3}}I_{1}-\frac{6}{\pi
u^{5}}I_{2}\right)\sin^{2}k+\left(-\mu-2u-B-\frac{2}{\pi u}I_{1}+\frac{2}{\pi
u^{3}}I_{2}\right)$ (135) $\displaystyle=$
$\displaystyle-2\cos{k}+2C_{1}\sin^{2}k+C_{2},$
$\displaystyle\varepsilon_{1}(\Lambda)$ $\displaystyle=$
$\displaystyle\Lambda^{2}\left(\eta_{1}+\frac{2J_{1}}{\pi
u^{3}}-\frac{12J_{2}}{\pi u^{5}}-\frac{I_{1}}{4\pi u^{3}}+\frac{3I_{2}}{8\pi
u^{5}}\right)$ (136)
$\displaystyle+\left(2B+4(u-\sqrt{1+u^{2}})-\frac{2J_{1}}{\pi
u}+\frac{2J_{2}}{\pi u^{3}}+\frac{I_{1}}{\pi u}-\frac{I_{2}}{4\pi
u^{3}}\right)$ $\displaystyle=$ $\displaystyle D_{1}\Lambda^{2}+D_{2},$

where the functions $C_{1},\,C_{2},\,D_{1},\,D_{2}$ denote the corresponding
coefficients

$\displaystyle C_{1}$ $\displaystyle=$ $\displaystyle\frac{1}{\pi
u^{3}}I_{1}-\frac{6}{\pi u^{5}}I_{2},$ (137) $\displaystyle C_{2}$
$\displaystyle=$ $\displaystyle-\mu-2u-B-\frac{2}{\pi u}I_{1}+\frac{2}{\pi
u^{3}}I_{2},$ (138) $\displaystyle D_{1}$ $\displaystyle=$
$\displaystyle\eta_{1}+\frac{2J_{1}}{\pi u^{3}}-\frac{12J_{2}}{\pi
u^{5}}-\frac{I_{1}}{4\pi u^{3}}+\frac{3I_{2}}{8\pi u^{5}},$ (139)
$\displaystyle D_{2}$ $\displaystyle=$ $\displaystyle
2B+4(u-\sqrt{1+u^{2}})-\frac{2J_{1}}{\pi u}+\frac{2J_{2}}{\pi
u^{3}}+\frac{I_{1}}{\pi u}-\frac{I_{2}}{4\pi u^{3}}.$ (140)

while the integrals $(I_{1},I_{2})$ and $(J_{1},J_{2})$ are related to spin
and charge degrees of freedom and make contributions only near zero point.
Integrating by parts, the formal solutions in terms of polylog functions read

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{d}\Lambda
T\ln(1+\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}})=-\frac{T^{\frac{3}{2}}\pi^{\frac{1}{2}}}{2D^{\frac{1}{2}}_{1}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{-D_{2}}{T}}\right),$
(141) $\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{d}\Lambda\Lambda^{2}T\ln(1+\mathrm{e}^{-\frac{\varepsilon_{1}(\Lambda)}{T}})=-\frac{T^{\frac{5}{2}}\pi^{\frac{1}{2}}}{4D^{\frac{3}{2}}_{1}}\mathrm{Li}_{\frac{5}{2}}\left(-\mathrm{e}^{\frac{-D_{2}}{T}}\right),$
(142) $\displaystyle J_{1}$ $\displaystyle=$ $\displaystyle
T\int_{0}^{\pi}\mathrm{d}k\cos k\ln(1+\mathrm{e}^{\frac{\kappa(k)}{T}})$ (143)
$\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),$
$\displaystyle J_{2}$ $\displaystyle=$ $\displaystyle
T\int_{0}^{\pi}\mathrm{d}k\cos
k\sin^{2}k\ln(1+\mathrm{e}^{\frac{\kappa(k)}{T}})=\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).$
(144)

With the help of the presentations of these four integrals, the Gibbs free
energy is given by

$f=-\mu-
u-B-2\lambda_{1}I_{1}-\lambda_{2}I_{2}+\frac{1}{\pi}J_{1}+\frac{1}{2\pi}J_{2}.$
(145)

Although we have greatly simplified the dressed equations, the four integrals
$(I_{1},I_{2})$ and $(J_{1},J_{2})$ are still coupled to each other. These
four quantities are intertwined and are needed to carefully distinguish the
primary and secondary contributions according to the critical field
conditions. Based on the above equations (135)-(144), we proceed to evaluate
critical behaviours from three aspects: critical fields, polylog functions and
free energy. The polylog functions contain the sources of criticality, from
which we can derive the scaling function of the free energy for quantum
critical region. For the part that causes criticality, we take expansions
under the limit of $T\gg\Delta B$ or $T\gg\Delta\mu$.

II-III phase transition: For the phase transition II-III, the critical field
is determined by $2-2u-B_{c}-\mu_{c}=0$ and $\varepsilon_{1}(\Lambda)$ is
gapped. The free energy is given by

$\displaystyle f=-\mu-
u-B+\frac{T^{\frac{3}{2}}}{2\pi^{\frac{1}{2}}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{-\Delta
B-\Delta\mu}{T}}\right),$ (146)

leading to the scattering forms of thermodynamics

$\displaystyle m$ $\displaystyle=$
$\displaystyle\frac{n}{2}=\frac{1}{2}+\frac{T^{\frac{1}{2}}}{4\pi^{\frac{1}{2}}}f_{\frac{1}{2}},$
(147) $\displaystyle\chi_{s}$ $\displaystyle=$
$\displaystyle\frac{\chi_{c}}{2}=-\frac{T^{-\frac{1}{2}}}{4\pi^{\frac{1}{2}}}f_{-\frac{1}{2}},$
(148) $\displaystyle\frac{C_{v}}{T}$ $\displaystyle=$
$\displaystyle-\frac{3T^{-\frac{1}{2}}}{8\pi^{\frac{1}{2}}}f_{\frac{3}{2}}-\frac{T^{-\frac{3}{2}}\left(\Delta
B+\Delta\mu\right)}{2\pi^{\frac{1}{2}}}f_{\frac{1}{2}}-\frac{T^{-\frac{5}{2}}\left(\Delta
B+\Delta\mu\right)^{2}}{2\pi^{\frac{1}{2}}}f_{-\frac{1}{2}}$ (149)

with $f_{n}=\mathrm{Li}_{n}\left(-\mathrm{e}^{\frac{-\Delta
B-\Delta\mu}{T}}\right)$. These present universal scaling behaviour for the
charge gapped phase transition in a lattice. The situation is very similar to
the case of I-II transition where only charge sector exists.

III-V phase transition: For the phase transition III-V, the emergence of spin
degree of freedom on a half-filled lattice generates the criticality with
respect to the critical field $B_{c}=2\sqrt{1+u^{2}}-2u$, independent of
chemical potential. This gives the integrable
$I_{1}\approx-\frac{T^{3/2}\pi^{1/2}}{2\eta^{1/2}_{1}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{-\frac{2\Delta
B}{T}}\right)$. In this case $\kappa(k)$ is always less than zero in charge
momentum space. For the charge degree, $|2+C_{2}|\gg T$, thus the $J_{1}$ with
respect to temperature magnitude can be evaluated by

$J_{1}=\frac{1}{2}\pi^{\frac{1}{2}}T^{\frac{3}{2}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{2+C_{2}}{T}}\right)\approx\frac{1}{2}\pi^{\frac{1}{2}}T^{\frac{3}{2}}\left[-\mathrm{e}^{\frac{2+C_{2}}{T}}+\frac{1}{2^{\frac{3}{2}}}\mathrm{e}^{\frac{2(2+C_{2})}{T}}\cdots\right]\ll\frac{1}{2}\pi^{\frac{1}{2}}T^{\frac{3}{2}}\approx
I_{1}.$ (150)

For this reason, $J_{1},J_{2}$ can be safely neglected, equivalent to the case
of XXX spin chain. Whereas for spin degrees of freedom, the term $I_{1}$ is
relevant to the low temperature criticality, in contrast the integral $I_{2}$
has higher order of power of the temperature than that of $I_{1}$. Under this
circumstance, we determine the free energy near the phase transition III-V as

$\displaystyle f=-\mu-
u-B+\frac{T^{\frac{3}{2}}\pi^{\frac{1}{2}}\lambda_{1}}{\eta^{\frac{1}{2}}_{1}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{-\frac{2\Delta
B}{T}}\right).$ (151)

Using the standard thermodynamic relation, we give the following forms of
scaling functions for density, magnetization, compressibility, susceptibility
and specific heat

$\displaystyle n$ $\displaystyle=$ $\displaystyle 1,\quad
m=\frac{1}{2}+\frac{T^{\frac{1}{2}}\pi^{\frac{1}{2}}\lambda_{1}}{\eta^{\frac{1}{2}}_{1}}f_{\frac{1}{2}},$
(152) $\displaystyle\chi_{c}$ $\displaystyle=$ $\displaystyle
0,\quad\chi_{s}=-\frac{2T^{-\frac{1}{2}}\pi^{\frac{1}{2}}\lambda_{1}}{\eta^{\frac{1}{2}}_{1}}f_{-\frac{1}{2}},$
(153) $\displaystyle\frac{C_{v}}{T}$ $\displaystyle=$
$\displaystyle-\frac{3T^{-\frac{1}{2}}\pi^{\frac{1}{2}}\lambda_{1}}{4\eta^{\frac{1}{2}}_{1}}f_{\frac{3}{2}}-\frac{2T^{-\frac{3}{2}}\pi^{\frac{1}{2}}\lambda_{1}\Delta
B}{\eta^{\frac{1}{2}}_{1}}f_{\frac{1}{2}}-\frac{4T^{-\frac{5}{2}}\pi^{\frac{1}{2}}\lambda_{1}\Delta
B^{2}}{\eta^{\frac{1}{2}}_{1}}f_{-\frac{1}{2}}$ (154)

with $f_{n}=\mathrm{Li}_{n}\left(-\mathrm{e}^{-\frac{2\Delta B}{T}}\right)$.

II-IV phase transition: For quantum phase transition II-IV, the phase
transition occurs at $\varepsilon_{1}(0)=0$ with charge dispersion
$\kappa(k)=-2\cos k-\mu-2u-B$. Let $Q$ denote Fermi point of $\kappa(k)$,
which satisfies $2\cos Q=-\mu-2u-B$. Via (81) corresponding to $\kappa(Q)=0$,
the critical field is obtained

$\displaystyle B_{c}$ $\displaystyle=$ $\displaystyle
2\sqrt{1+u^{2}}-2u-\frac{2}{3\pi u}(\pi-Q)^{3},$ $\displaystyle\mu_{c}$
$\displaystyle=$ $\displaystyle-2\cos
Q-2u-B_{c}=2-2\sqrt{1+u^{2}}-(\pi-Q)^{2}+\frac{2}{3\pi u}(\pi-Q)^{3}.$ (155)

From (135) and (136) at criticality, we observe that the term $C_{1}$ is a
small quantity, while $\frac{2+C_{2}}{T}=\frac{\kappa(\pi)}{T}$ in polylog
function in $J_{1},\,J_{2}$ is a large quantity. Using the property of the
polylog function

$\displaystyle\mathrm{Li}_{p}(-\mathrm{e}^{w})\stackrel{{\scriptstyle
w>>1}}{{\longrightarrow}}-2\sum_{k=0}^{\infty}\eta(2k)\frac{w^{p-2k}}{\Gamma(p-2k+1)}+O(-\mathrm{e}^{w}),$
(156)

we have

$\displaystyle J_{1}$ $\displaystyle\approx$
$\displaystyle\frac{1}{2}(1+C_{1})\left[-\frac{4}{3}\left(\kappa(\pi)\right)^{\frac{3}{2}}-\frac{\pi^{2}T^{2}}{6}\left(\frac{1}{\kappa(\pi)}\right)^{\frac{1}{2}}\right]$
(157)
$\displaystyle+\frac{3}{32}(1+5C_{1})\left[\frac{8}{15}\left(\kappa(\pi)\right)^{\frac{5}{2}}+\frac{\pi^{2}T^{2}}{3}\left(\kappa(\pi)\right)^{\frac{1}{2}}\right],$
$\displaystyle J_{2}$ $\displaystyle\approx$
$\displaystyle-\frac{1}{4}(1+3C_{1}+\frac{15}{2}C^{2}_{1})\left[\frac{8}{15}\left(\kappa(\pi)\right)^{\frac{5}{2}}+\frac{\pi^{2}T^{2}}{3}\left(\kappa(\pi)\right)^{\frac{1}{2}}\right].$
(158)

While $I_{1},\,I_{2}$ indicate that they are the orders of
$O(T^{\frac{3}{2}})$ and $O(T^{\frac{5}{2}})$, respectively. Therefore we
reasonably ignore the effect of $I_{2}$. The leading term $I_{1}$ is given by

$\displaystyle
I_{1}\approx-\frac{T^{\frac{3}{2}}\pi^{\frac{1}{2}}}{2D^{\frac{1}{2}}_{1}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{-\frac{D_{2}}{T}}\right).$
(159)

In the above equations (157) and (158), we have an approximation
$C_{1}\approx\frac{1}{\pi u^{3}}I_{1},C_{2}\approx-\mu-2u-B-\frac{2}{\pi
u}I_{1}$. It follows that

$\displaystyle\kappa(\pi)$ $\displaystyle=$ $\displaystyle
2+C_{2}=2-\mu-2u-B-\frac{2}{\pi u}I_{1}=2+2\cos Q-\frac{2}{\pi
u}I_{1}=\beta_{1}-\frac{2}{\pi u}I_{1},$ (160)

where $\beta_{1}=2-\mu-2u-B$ which is defined previously. Follow these
approximations, we rewrite the functions $J_{1}$ and $J_{2}$

$\displaystyle J_{1}$ $\displaystyle\approx$
$\displaystyle-\frac{2}{3}\beta_{1}^{\frac{3}{2}}+\frac{2\beta_{1}^{\frac{1}{2}}}{\pi
u}I_{1}-\frac{\pi^{2}T^{2}}{12\beta_{1}^{\frac{1}{2}}}-\frac{2\beta_{1}^{\frac{3}{2}}}{3\pi
u^{3}}I_{1}+\frac{1}{20}\beta_{1}^{\frac{5}{2}}-\frac{\beta_{1}^{\frac{3}{2}}}{4\pi
u}I_{1}+\frac{\pi^{2}T^{2}\beta_{1}^{\frac{1}{2}}}{32},$ (161) $\displaystyle
J_{2}$ $\displaystyle\approx$
$\displaystyle-\frac{2}{15}\beta_{1}^{\frac{5}{2}}+\frac{2\beta_{1}^{\frac{3}{2}}}{3\pi
u}I_{1}-\frac{\pi^{2}T^{2}\beta_{1}^{\frac{1}{2}}}{12}.$ (162)

Substituting (161) to the definition of $D_{2}$ in (140) which reflects the
distance away from QCP, $D_{2}$ can then be represented with $\Delta
B,\,\Delta\mu$ as

$\displaystyle\Delta t\equiv:D_{2}$ $\displaystyle\approx$ $\displaystyle
2B+4(u-\sqrt{1+u^{2}})-\frac{2J_{1}}{\pi u}\approx
2B+4(u-\sqrt{1+u^{2}})+\frac{4}{3\pi u}\beta_{1}^{\frac{3}{2}}$ (163)
$\displaystyle=$ $\displaystyle 2B+4(u-\sqrt{1+u^{2}})+\frac{4}{3\pi
u}\left(\beta_{1c}^{\frac{3}{2}}-\frac{3\beta_{1c}^{\frac{3}{2}}}{2}(\Delta
B+\Delta\mu)\right)-2B_{c}+2B_{c}$ $\displaystyle=$ $\displaystyle
2\left(1-\frac{\beta_{1c}^{\frac{1}{2}}}{\pi u}\right)\Delta
B-\frac{2\beta_{1c}^{\frac{1}{2}}}{\pi u}\Delta\mu,$

where $\beta_{1c}=2-\mu_{c}-2u-B_{c}$. In the above equation, we only inserted
the leading term $-\frac{2}{3}\beta_{1}^{\frac{3}{2}}$ of $J_{1}$ that is
irrelevant of $T$. While $J_{2}$ has been neglected and the identity
$2-\beta_{1}-2u-B-\mu=0$ has been used. Substituting the results of
$J_{1},\,J_{2},\,I_{1}$ into the free energy (145), we have

$f\approx-\mu-
u-B-\frac{2}{3\pi}\beta_{1}^{\frac{3}{2}}-\frac{1}{60\pi}\beta_{1}^{\frac{5}{2}}-\frac{\pi
T^{2}}{12\beta_{1}^{\frac{1}{2}}}-\frac{\pi
T^{2}\beta_{1}^{\frac{1}{2}}}{96}+\frac{T^{\frac{3}{2}}\pi^{\frac{1}{2}}a_{0}}{\eta^{\frac{1}{2}}_{1}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{-\frac{\Delta
t}{T}}\right),$ (164)

here
$a_{0}=\lambda_{1}-\frac{\beta_{1}^{\frac{1}{2}}}{\pi^{2}u}+\beta_{1}^{\frac{3}{2}}\left(\frac{2\lambda_{1}}{3\pi
u^{3}\eta_{1}}-\frac{1}{24\pi^{2}u}+\frac{1}{3\pi^{2}u^{3}}\right)$ and
$\Delta t=2\left(1-\frac{\beta_{1c}^{\frac{1}{2}}}{\pi u}\right)\Delta
B-\frac{2\beta_{1c}^{\frac{1}{2}}}{\pi u}\Delta\mu$. Taking derivatives of
free energy with respect to $B,\,\mu$ and $T$, we obtain the following scaling
forms of the thermal and magnetic properties

$\displaystyle n$ $\displaystyle=$ $\displaystyle
1-\frac{\beta_{1}^{\frac{1}{2}}}{\pi}-\frac{\beta_{1}^{\frac{3}{2}}}{24\pi}-\frac{2T^{\frac{1}{2}}\beta_{1c}^{\frac{1}{2}}a_{0}}{\pi^{\frac{1}{2}}u\eta^{\frac{1}{2}}_{1}}f_{\frac{1}{2}},$
(165) $\displaystyle m$ $\displaystyle=$
$\displaystyle\frac{1}{2}-\frac{\beta_{1}^{\frac{1}{2}}}{2\pi}-\frac{\beta_{1}^{\frac{3}{2}}}{48\pi}+\frac{T^{\frac{1}{2}}\left(1-\frac{\beta_{1c}^{\frac{1}{2}}}{\pi
u}\right)\pi^{\frac{1}{2}}a_{0}}{\eta^{\frac{1}{2}}_{1}}f_{\frac{1}{2}},$
(166) $\displaystyle\chi_{c}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\beta_{1}^{\frac{1}{2}}}+\frac{\beta_{1}^{\frac{1}{2}}}{16\pi}-\frac{4T^{-\frac{1}{2}}\beta_{1c}a_{0}}{\pi^{\frac{3}{2}}u^{2}\eta^{\frac{1}{2}}_{1}}f_{-\frac{1}{2}},$
(167) $\displaystyle\chi_{s}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi\beta_{1}^{\frac{1}{2}}}+\frac{\beta_{1}^{\frac{1}{2}}}{32\pi}-\frac{2T^{-\frac{1}{2}}\left(1-\frac{\beta_{1c}^{\frac{1}{2}}}{\pi
u}\right)^{2}\pi^{\frac{1}{2}}a_{0}}{\eta^{\frac{1}{2}}_{1}}f_{-\frac{1}{2}},$
(168) $\displaystyle\frac{C_{v}}{T}$ $\displaystyle=$
$\displaystyle\frac{\pi}{6\beta_{1}^{\frac{1}{2}}}+\frac{\pi\beta_{1}^{\frac{1}{2}}}{48}-\frac{3T^{-\frac{1}{2}}\pi^{\frac{1}{2}}a_{0}}{4\eta^{\frac{1}{2}}_{1}}f_{\frac{3}{2}}-\frac{2T^{-\frac{3}{2}}\pi^{\frac{1}{2}}a_{0}}{\eta^{\frac{1}{2}}_{1}}\left[\left(1-\frac{\beta_{1c}^{\frac{1}{2}}}{\pi
u}\right)\Delta B-\frac{\beta_{1c}^{\frac{1}{2}}}{\pi
u}\Delta\mu\right]f_{\frac{1}{2}}$ (169)
$\displaystyle-\frac{4T^{-\frac{5}{2}}\pi^{\frac{1}{2}}a_{0}}{\eta^{\frac{1}{2}}_{1}}\left[\left(1-\frac{\beta_{1c}^{\frac{1}{2}}}{\pi
u}\right)\Delta B-\frac{\beta_{1c}^{\frac{1}{2}}}{\pi
u}\Delta\mu\right]^{2}f_{-\frac{1}{2}}$

with $f_{n}=\mathrm{Li}_{n}\left(-\mathrm{e}^{-\frac{\Delta t}{T}}\right)$.

It is further noticed from (164) that the background of free energy involves
two types of terms, the temperature-dependent term $-\frac{\pi
T^{2}}{12\beta_{1}^{\frac{1}{2}}}-\frac{\pi T^{2}\beta_{1}^{\frac{1}{2}}}{96}$
comes from the contribution of TLL, whereas the remining terms $-\mu-
u-B-\frac{2}{3\pi}\beta_{1}^{\frac{3}{2}}-\frac{1}{60\pi}\beta_{1}^{\frac{5}{2}}$
is equivalent to the energy of non-interacting electrons on a lattice. For
simplicity and without losing generality, we only keep the order of
$\beta_{1}^{\frac{3}{2}}$ in the free energy. For the 1D repulsive Hubbard
model, zero-temperature background is given by

$\displaystyle f_{\mathrm{H}}\approx-\mu-
u-B-\frac{2}{3\pi}\beta_{1}^{\frac{3}{2}}\approx\beta_{1}-2-\frac{2}{3\pi}\beta_{1}^{\frac{3}{2}}.$
(170)

While the ground state free energy of a non-interacting lattice system is
given by

$\displaystyle
f_{\mathrm{NI}}=-\frac{1}{\pi}\left[\sum_{\uparrow,\downarrow}2\sin(k_{F,\sigma})+\mu_{\sigma}\sin(k_{F,\sigma})\right],$
(171)

in which $k_{F,\sigma}=\arccos(-\mu_{\sigma}/2)=\pi
n_{\sigma},\mu_{\uparrow}=\mu+B,\mu_{\downarrow}=\mu-B$. In the limit of high
density at spin-polarized phase, these two free energies (170) and (171) are
equivalent

$f_{\mathrm{NI}}\approx-\frac{2}{\pi}\left(\beta_{1}^{\frac{1}{2}}-\frac{1}{6}\beta_{1}^{\frac{3}{2}}\right)-(\mu+B)n\approx-\frac{2}{\pi}\left(\beta_{1}^{\frac{1}{2}}-\frac{1}{6}\beta_{1}^{\frac{3}{2}}\right)-(2-\beta_{1})\left(1-\frac{\beta_{1}^{\frac{1}{2}}}{\pi}\right)=f_{\mathrm{H}}.$
(172)

In light of quantum criticality, we refer to the phase II as the background.
Thus the terms that do not contribute to the singular part in thermodynamic
quantities can be readily recognised, namely,

$\displaystyle\chi^{\mathrm{II}}_{c}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\beta_{1}^{\frac{1}{2}}}+\frac{\beta_{1}^{\frac{1}{2}}}{16\pi},$
(173) $\displaystyle\chi^{\mathrm{II}}_{s}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi\beta_{1}^{\frac{1}{2}}}+\frac{\beta_{1}^{\frac{1}{2}}}{32\pi},$
(174) $\displaystyle\frac{C^{\mathrm{II}}_{v}}{T}$ $\displaystyle=$
$\displaystyle\frac{\pi}{6\beta_{1}^{\frac{1}{2}}}+\frac{\pi\beta_{1}^{\frac{1}{2}}}{48}.$
(175)

These background results are in agreement with (122) with
$\delta=\beta_{1}^{\frac{1}{2}}\left(1+\frac{1}{24}\beta_{1}\right)$, here
$\beta_{1}=2-\mu-2u-B\approx(\pi-Q)^{2}-\frac{1}{12}(\pi-Q)^{4}=\delta^{2}-\frac{1}{12}\delta^{4}$.
Furthermore the Wilson ratios for the phase II are given by

$R^{\chi_{s}}_{w}=\frac{4\pi^{2}}{3}\frac{\chi_{s}}{C_{v}/T}\approx 2,\qquad
R^{\chi_{c}}_{w}=\frac{\pi^{2}}{3}\frac{\chi_{c}}{C_{v}/T}\approx 1.$ (176)

IV-V phase transition: The phase transition from IV to V displays a novel
subtlety of quantum criticality with charge. At this phase transition, the
charge degree of freedom becomes gapped and the dressed energy $\kappa(\pi)$
approaches zero. We can directly use the dressed energy equations with the
help of (82) and $\kappa(\pi)=0$ under the condition $\delta=0$, giving the
critical fields

$\displaystyle\mu_{c}$ $\displaystyle=$ $\displaystyle
2-2\sqrt{1+u^{2}}+\frac{1}{(1+u^{2})^{\frac{3}{2}}}A^{2}_{c}\left(1-\frac{2}{\pi
u}A_{c}\right),$ (177) $\displaystyle B_{c}$ $\displaystyle=$
$\displaystyle-2u+2\sqrt{1+u^{2}}-\frac{1}{(1+u^{2})^{\frac{3}{2}}}A^{2}_{c}\left(1+\frac{2}{3\pi
u}A_{c}\right).$ (178)

From (78) and (84), it is convenient to represent $A_{c}$ in term of $B_{c}$

$A_{c}^{2}=-\frac{\beta_{2}}{\eta_{1}}-\frac{2}{3\pi
u}\frac{(-\beta_{2})^{\frac{3}{2}}}{\eta^{\frac{3}{2}}_{1}}+\frac{2}{3\pi^{2}u^{2}}\frac{(-\beta_{2})^{2}}{\eta^{2}_{1}},$
(179)

with $\eta_{1}=\frac{2}{(1+u^{2})^{3/2}}$. Due to the criticality induced from
charge degree, precisely opposite to the transition from phase II to IV, we
rewrite the (178) as

$\displaystyle 2B+4(u-\sqrt{1+u^{2}})+\frac{2\eta_{1}}{3\pi
u}A^{3}+\eta_{1}A^{2}=0.$ (180)

This equation represents the boundary condition. Using the relation $-D_{2}$
at zero temperature is approximate by $\eta_{1}A^{2}$, we get the main
contribution $I_{1}=2\eta_{1}A^{3}/3$ of the first iteration. Comparing this
equation with the boundary condition and (139)(140), up to the leading
contribution, we have $-D_{2}\approx\eta_{1}A^{2}+2J_{1}/(\pi u)$, which is
the leading order in the arguments of the polylog functions in $I_{1}$ and
$I_{2}$. The integral $I_{1}$ and $I_{2}$ serve as the background
contributions resulting from of the spin degrees of freedom of the TLL. Up to
the leading order of the temperature contributions and using the approximation
(156), we obtain the function for the second iteration

$\displaystyle I_{1}$ $\displaystyle\approx$
$\displaystyle\frac{2}{3}\eta_{1}A^{3}+\frac{\pi^{2}T^{2}}{12\eta_{1}A}+\frac{2A}{\pi
u}J_{1}-\frac{2A^{3}}{3\pi u^{3}}J_{1}+\frac{\pi
T^{2}}{144u^{3}\eta_{1}}A^{2},$ (181) $\displaystyle I_{2}$
$\displaystyle\approx$
$\displaystyle\frac{2}{15}\eta_{1}A^{5}+\frac{\pi^{2}T^{2}}{12\eta_{1}}A+\frac{2A^{3}}{3\pi
u}J_{1}.$ (182)

The term $2+C_{2}$ in the charge dressed energy determines the criticality in
the charge degrees of freedom. Substituting the above expressions of $I_{1}$
and $I_{2}$ into the $C_{2}$ in (138) and using the critical fields (177) and
(178), we give

$\Delta t\equiv:2+C_{2}=-\Delta\mu-\left(1-\frac{4A_{c}}{\pi
u+A_{c}}\right)\Delta B.$ (183)

Subsequently, the scaling function $J_{1}$ is given by

$J_{1}=\frac{1}{2}T^{\frac{3}{2}}\pi^{\frac{1}{2}}\left(1+\frac{2\eta_{1}}{3\pi
u^{3}}A^{3}\right)\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{\Delta
t}{T}}\right).$ (184)

Using the leading order contribution in $J_{1}$ and $I_{1},\,I_{2}$, we
finally obtain the following scaling form of the Gibbs free energy

$f\approx-\mu-
u-B-\frac{4}{3}\lambda_{1}\eta_{1}A^{3}-\frac{\pi^{2}T^{2}\lambda_{1}}{6\eta_{1}A}-\frac{2}{15}\lambda_{2}\eta_{1}A^{5}-\frac{\pi^{2}T^{2}\lambda_{2}}{12\eta_{1}}A+\frac{T^{\frac{3}{2}}b_{0}}{2\pi^{\frac{1}{2}}}\mathrm{Li}_{\frac{3}{2}}\left(-\mathrm{e}^{\frac{\Delta
t}{T}}\right)$ (185)

with
$b_{0}=1-\frac{4\lambda_{1}}{u}A+\left(\frac{\eta_{1}}{\pi}+2\lambda_{1}-u^{2}\lambda_{2}\right)\frac{2A^{3}}{3u^{3}}$.
Using the thermodynamic relations, the scaling forms of thermal and magnetic
properties are given explicitly by

$\displaystyle n$ $\displaystyle=$ $\displaystyle
1+\frac{T^{\frac{1}{2}}b_{0}}{2\pi^{\frac{1}{2}}}f_{\frac{1}{2}},$ (186)
$\displaystyle m$ $\displaystyle=$
$\displaystyle\frac{1}{2}-\frac{2\lambda_{1}A}{1+\frac{A}{\pi
u}}-\frac{1}{3}\lambda_{2}A^{3}+\frac{T^{\frac{1}{2}}\left(1-\frac{4A_{c}}{\pi
u+A_{c}}\right)b_{0}}{4\pi^{\frac{1}{2}}}f_{\frac{1}{2}},$ (187)
$\displaystyle\chi_{c}$ $\displaystyle=$
$\displaystyle-\frac{T^{-\frac{1}{2}}b_{0}}{2\pi^{\frac{1}{2}}}f_{-\frac{1}{2}},$
(188) $\displaystyle\chi_{s}$ $\displaystyle=$
$\displaystyle\frac{2\lambda_{1}}{\eta_{1}A\left(1+\frac{A}{\pi
u}\right)^{3}}+\frac{\lambda_{2}}{\eta_{1}}\frac{A}{1+\frac{A}{\pi
u}}-\frac{T^{-\frac{1}{2}}\left(1-\frac{4A_{c}}{\pi
u+A_{c}}\right)^{2}b_{0}}{4\pi^{\frac{1}{2}}}f_{-\frac{1}{2}},$ (189)
$\displaystyle\frac{C_{v}}{T}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}\lambda_{1}}{3\eta_{1}A}+\frac{\pi^{2}\lambda_{2}A}{6\eta_{1}}-\frac{3T^{-\frac{1}{2}}b_{0}}{8\pi^{\frac{1}{2}}}f_{\frac{3}{2}}-\frac{T^{-\frac{3}{2}}b_{0}}{2\pi^{\frac{1}{2}}}\left[\Delta\mu+\left(1-\frac{4A_{c}}{\pi
u+A_{c}}\right)\Delta B\right]f_{\frac{1}{2}}$ (190)
$\displaystyle-\frac{T^{-\frac{5}{2}}b_{0}}{2\pi^{\frac{1}{2}}}\left[\Delta\mu+\left(1-\frac{4A_{c}}{\pi
u+A_{c}}\right)\Delta B\right]^{2}f_{-\frac{1}{2}},$

where we denoted $f_{n}=\mathrm{Li}_{n}\left(-\mathrm{e}^{\frac{\Delta
t}{T}}\right)$.

Moreover, we observe that the phase V represents the source of the background
for the spin degrees of freedom. Thus we may easily identify the contributions
from Mott phase V in the above scaling functions, namely

$\displaystyle\chi^{\mathrm{V}}_{c}$ $\displaystyle=$ $\displaystyle 0,$ (191)
$\displaystyle\quad\chi^{\mathrm{V}}_{s}$ $\displaystyle=$
$\displaystyle\frac{2\lambda_{1}}{\eta_{1}A\left(1+\frac{A}{\pi
u}\right)^{3}}+\frac{\lambda_{2}}{\eta_{1}}\frac{A}{1+\frac{A}{\pi u}},$ (192)
$\displaystyle\frac{C^{\mathrm{V}}_{v}}{T}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}\lambda_{1}}{3\eta_{1}A}+\frac{\pi^{2}\lambda_{2}A}{6\eta_{1}}.$
(193)

Comparing above results with equations (123) to (127), we find these results
consistent with each other. Consequently, the Wilson ratios are given by

$R^{\chi_{s}}_{w}=\frac{4\pi^{2}}{3}\frac{\chi_{s}}{C_{v}/T}\approx\frac{8}{\left(1+\frac{A}{\pi
u}\right)^{3}}+\frac{8\lambda_{2}A^{3}}{\pi u\lambda_{1}},\quad
R^{\chi_{c}}_{w}=\frac{\pi^{2}}{3}\frac{\chi_{c}}{C_{v}/T}=0,$ (194)

showing the nature of the TLL at low temperature.

Figure 12: Row (a) ((b)) shows universal scaling behaviour of electron
density, spin susceptibility as well as specific heat v.s. magnetic field $B$
(the argument $(B-B_{c})/T$) for phase transitions II-IV. The analytical
results (solid lines) (165)-(169) agree well with the numerical solutions of
the TBA equations. The parameter setting is $u=1,\,\mu=-0.8382$. Row (c) ((d))
shows universal scaling behaviour of magnetization, charge susceptibility as
well as specific heat v.s. magnetic field $B$ (the argument $(B-B_{c})/T$) for
phase transitions IV-V. The analytical results (solid lines) (186)-(190) agree
well with the numerical solution of the TBA equations. The parameter setting
is $u=1,\,\mu=-0.8272$.

In figure 12 (a) (b) and figure 12 (c) (d), we plot universal scaling
behaviour of electron density and magnetization, spin and charge
susceptibilities as well as specific heat for phase transitions II-IV and
IV-V, respectively. They show that analytic expressions of their scaling
functions (165)-(169) and (186)-(190) are in good agreement with numerical
results obtained from the TBA equations. We note that in figure 12 (a), (c)
all lines at different temperatures intersect at QCP, while figure 12 (b), (d)
show the scaling function invariant in terms of $\frac{\Delta B}{T}$. It is
essential to note that all these scaling functions read off the dynamical
exponent $z=2$ and correlation critical exponent $\nu=1/2$, for example, the
susceptibility (129), also see sachdev2011quantum . In the critical region,
the polylog function represents the free fermion type of generating function
associated with dynamic critical exponent $z=2$ and correlation critical
exponent $v=\frac{1}{2}$ sachdev2011quantum ; continentino2017quantum . In the
zero temperature limit, the singular part of the susceptibility can be