# Finite time effects in single and double Compton scattering

V. Dubrovich1, T. Zalialiutdinov2,3 1 Special Astrophysical Observatory, St.
Petersburg Branch, Russian Academy of Sciences, 196140, St. Petersburg, Russia
2 Department of Physics, St. Petersburg State University, Petrodvorets,
Oulianovskaya 1, 198504, St. Petersburg, Russia
3 Petersburg Nuclear Physics Institute named by B.P. Konstantinov of National
Research Centre ”Kurchatov Institut”, St. Petersburg, Gatchina 188300, Russia
<EMAIL_ADDRESS>
###### Abstract

The process of Compton scattering by a free electron with subsequent
reemission of one or two photons is considered in the assumption of finite
interaction time. The corresponding cross sections are obtained in the
framework of relativistic quantum electrodynamics using a modified form of
fermion propagator with complex transmitted momentum. It is shown that finite
time effects can be observable at sufficiently low energies of scattered
photons. The proposed method also regularizes arising infrared divergence in
the cross section of the double Compton effect. Possible experimental
verification of considered theoretical approach is discussed.

##  I Introduction

In this paper, we consider the derivation of single and double Compton
scattering cross section, assuming that the incoming photon interacts with the
electron for only a finite time $\tau$. The explicit introduction of time-
dependent adiabatic factors at the interaction vertex immediately makes the
theory no longer Lorentz-invariant [1]. It complicates the construction of
quantities such as the transition probabilities and the scattering cross
sections. However, the approach developed for interactions with unstable
particles proves helpful here. Recently, it was shown that for unstable
fermion particles, the following modified Feynman propagator might be
introduced [2, 3, 4, 5]

$\displaystyle
S_{F}(p)=\frac{\hat{p}+m_{F}-\mathrm{i}\Gamma/2}{p^{2}-(m_{F}-\mathrm{i}\Gamma/2)^{2}},$
(1)

where $p$ is the four-vector of fermion momentum, $m_{F}$ is the mass of
fermion, $\hat{p}=\gamma^{\mu}p_{\mu}$, $\Gamma$ corresponds to the decay
width of the particle and, in general, could depend on transmitted momentum
$q$. In contrast to the naive introduction of factor $\Gamma$ in the numerator
of ordinary Feynman propagator [6]

$\displaystyle
S_{F}(p)=\frac{\hat{p}+m_{F}}{p^{2}-m_{F}^{2}-\mathrm{i}m_{F}\Gamma},$ (2)

the representation Eq. (1) satisfy the electromagnetic Ward identity, which
provides the gauge-invariant description of processes [7, 5].

Parameter $\Gamma$ introduced in Eq. (1) and associated with the particle
decay implies the finiteness of time interaction in the scattering processes.
The corresponding interaction interval is defined by the fermion decay time
$\tau=1/\Gamma$. In principle, the localaized in time effects for scattering
of stable particles can be introduced in a similar manner. Then the
corresponding $\Gamma$ in Eq. (1) can be associated with the duration time of
laser pulse $\tau_{L}$. This picture can be true for a single quasi-
monochromatic photon with incident energies much higher than the width of
photon non-monochromaticity. Without loss of generality, the imaginary
addition to the momentum transferred from a photon to a fermion can be
attributed either to the lifetime of an unstable particle, or to the
interaction time (inverse width) of a quasi-monochromatic photon with a stable
fermion.

It should be noted that the exactly similar situation arises in atomic physics
for the scattering of photons by an atomic electron. In this case, parameter
$\Gamma$ is associated with the natural width of atomic level, which can be
rigorously introduced within bound state QED theory of line profile or within
the quantum mechanical description of quasi-stationary states [8, 9]. However,
in contrast to QED in external fields the small imaginary addition to the
energy of free stable fermion can be introduced only phenomenologically. With
such restrictions and with the use of propagator in the form of Eq. (1), we
try to construct the cross section for the scattering of quasi-monochromatic
photon by a free fermion.

The paper is organised as follows. In section II fully relativistic QED
derivation of single Compton scattering of quasi-monochromatic photon is
presented. The section III is devoted to the similar consideration of double
Compton effect. The results and discussions are presented in section IV.

Throughout this paper relativistic units are employed, where the velocity of
light $c=1$, Planck’s constant $\hbar=1$, respectively. The electron mass
$m_{e}$ is written explicitly. The electron charge $e=-\left|e\right|$ is
related to the fine structure constant $\alpha$ via $\alpha=e^{2}/4\pi$. For
$4-$vectors and tensor a standard notations are used for covariant (lower
index) and contravariant (upper index) components, which are related to each
other by metric tensor with Minkowski metrics
$g_{\mu\nu}=g^{\mu\nu}=(1,-1,-1,-1)$, so that, for example, for vector
components $a_{\mu}$: $a_{\mu}=g_{\mu\nu}a^{\nu}$ (implying here and further
Einstein’s sum convention). Greek indices run over set $(0,\,1,\,2,\,3)$ and
Latin indices run over set $(1,\,2,\,3)$. The contraction of Dirac
$\gamma^{\mu}$-matrices ($\mu=0,\,1,\,2,\,3$) with four-vector
$a=(a_{0},\textbf{a})$ is denoted as $\hat{a}=\gamma^{\mu}a_{\mu}$.

##  II Single compton scattering

We start from the standard QED derivation of Klein-Nishina formula [10] for
Compton process but with the modified propagator in the form given by Eq. (1).
The corresponding S-matrix element for the single photon scattering
$\gamma(k_{1})+e^{-}\rightarrow\gamma(k_{2})+e^{-}$ is [11, 12]

$\displaystyle S^{(2)}_{fi}=(2\pi)^{4}\delta^{(4)}(p_{f}+k_{2}-p_{i}-k_{1})$
(3)
$\displaystyle\times\frac{m_{e}e^{2}}{\sqrt{8V^{4}E_{i}E_{f}\omega_{1}\omega_{2}}}{\cal
M}(k_{1},k_{2}),$

where $p_{i}=(E_{i},\textbf{p}_{i})$, $p_{f}=(E_{f},\textbf{p}_{f})$ are the
four-vectors of initial and final electron momentum, respectively,
$k_{1}=(\omega_{1},\textbf{k}_{1})$ is the four vector of incident photon with
frequency $\omega_{1}$ and wave-vector $\textbf{k}_{1}$,
$k_{2}=(\omega_{2},\textbf{k}_{2})$ is the four-vectors of outgoing photon,
and $\cal{M}$ is the Feynman amplitude of the process

$\displaystyle{\cal
M}(k_{1},k_{2})=\overline{u}(p_{f})\left\\{\hat{\varepsilon}(k_{2})S_{F}(p_{i}+k_{1})\hat{\varepsilon}(k_{1})\right.$
(4)
$\displaystyle+\left.\hat{\varepsilon}(k_{1})S_{F}(p_{i}-k_{2})\hat{\varepsilon}(k_{2})\right\\}u(p_{i}).$

In Eq. (4) $u(p)$ is the Dirac spinor for free electron with its Dirac adjoint
defined as $\overline{u}(p)=u^{\dagger}(p)\gamma_{0}$, $\gamma^{\mu}$ are the
Dirac matrices, $\varepsilon(k)$ is the polarization four-vector of a photon
with momentum $k$ and $S_{F}$ is the Feynman propagator given by the
representation Eq. (1).

The transition rate per unit time to one define state can be found according
to the definition [11]

$\displaystyle
w=\frac{|S^{(2)}_{fi}|^{2}}{T}=V(2\pi)^{4}\delta^{(4)}(p_{f}+k_{2}-p_{i}-k_{1})$
(5)
$\displaystyle\times\frac{m_{e}^{2}e^{4}}{4V^{4}E_{i}E_{f}\omega_{1}\omega_{2}}|{\cal
M}(k_{1},k_{2})|^{2},$

where $T\rightarrow\infty$ is the observation time and $V$ is the phase
volume. Since we are interested to a transition rate $dw$ to a group of final
states with momenta in the intervals
$(\textbf{p}_{f},\textbf{p}_{f}+d\textbf{p}_{f})$ and
$(\textbf{k}_{2},\textbf{k}_{2}+d\textbf{k}_{2})$ we must multiply Eq. (5) by
the number of these states which is

$\displaystyle\frac{V^{2}d^{3}\textbf{p}_{f}d^{3}\textbf{k}_{2}}{(2\pi)^{6}}.$
(6)

With our choice of normalization for the states, the volume $V$ contains one
scattering center and the incident photon flux is $F=c/V$ ($c$ is the speed of
light) [13]. Then the corresponding differential cross section can be found as
follows

$\displaystyle
d\sigma_{\mathrm{sc}}=\frac{dw}{F}=\delta^{(4)}(p_{f}+k_{2}-p_{i}-k_{1})$ (7)
$\displaystyle\times\frac{m_{e}^{2}e^{4}}{4E_{i}E_{f}\omega_{1}\omega_{2}}|{\cal
M}(k_{1},k_{2})|^{2}\frac{d^{3}\textbf{p}_{f}d^{3}\textbf{k}_{2}}{(2\pi)^{2}}.$

Integration over $d\textbf{p}_{f}$ in Eq. (7) can be performed with the use of
equality [11, 14]

$\displaystyle\frac{d^{3}\textbf{p}_{f}}{2E_{f}}=\int
d^{4}p_{f}\delta(p_{f}^{2}-m_{e}^{2})\theta(p_{f}^{0}),$ (8)

where $\theta$ is the Heaviside step-function. After that Eq. (7) takes the
form

$\displaystyle
d\sigma_{\mathrm{sc}}=\delta((p_{i}+k_{1}-k_{2})^{2}-m_{e}^{2})\theta(m_{e}+\omega_{1}-\omega_{2})$
(9) $\displaystyle\times\frac{m_{e}e^{4}}{2\omega_{1}\omega_{2}}|{\cal
M}(k_{1},k_{2})|^{2}\frac{d^{3}\textbf{k}_{2}}{(2\pi)^{2}}.$

We now specialize to a reference frame where the electron initially at rest,
i.e. $\textbf{p}_{i}=0$ and $E_{i}=m_{e}$. Then the energy conservation law
given by the argument of delta-function in Eq. (9) is

$\displaystyle(p_{i}+k_{1}-k_{2})^{2}-m_{e}^{2}$ (10)
$\displaystyle=2m_{e}(\omega_{1}-\omega_{2})-2\omega_{1}\omega_{2}(1-\cos\theta_{12})=0,$

and the frequency of outgoing photon $\omega_{2}$ can be expressed via
$\omega_{1}$ and the angle $\theta_{12}$ between corresponding photon wave-
vectors $\textbf{k}_{1}$ and $\textbf{k}_{2}$ as follows

$\displaystyle\omega_{2}=\frac{\omega_{1}}{1+\frac{\omega_{1}}{m_{e}}(1-\cos\theta_{12})}.$
(11)

Finally, performing summation over the spin and polarizations of final states
and averaging over the spins and polarizations of initial states the
differential cross section takes the form

$\displaystyle\frac{d\sigma_{\mathrm{sc}}}{d\Omega_{2}}=\frac{e^{4}}{16\pi^{2}}\left(\frac{\omega_{2}}{\omega_{1}}\right)^{2}X(\omega_{1},\omega_{2}),$
(12)

together with notations

$\displaystyle
X(\omega_{1},\omega_{2})=\frac{1}{2m_{e}^{2}}F(\omega_{1},\omega_{2})$ (13)
$\displaystyle\times\left(\Gamma^{4}+16m_{e}^{2}\left(\Gamma^{2}+4\omega_{1}^{2}\right)+16\Gamma^{2}m_{e}\omega_{1}\right)^{-1}$
$\displaystyle\times\left(\Gamma^{4}+16m_{e}^{2}\left(\Gamma^{2}+4\omega_{2}^{2}\right)-16\Gamma^{2}m_{e}\omega_{2}\right)^{-1},$

and

$\displaystyle
F(\omega_{1},\omega_{2})=64m_{e}\left(\Gamma^{6}(\omega_{1}-\omega_{2})\right.$
(14)
$\displaystyle+4\Gamma^{2}m_{e}^{2}(\omega_{1}-\omega_{2})\left(9\Gamma^{2}+4\omega_{1}^{2}+\omega_{2}^{2}\right)$
$\displaystyle+8m_{e}^{3}\left(4\Gamma^{2}\left(\omega_{1}^{2}-4\omega_{1}\omega_{2}+\omega_{2}^{2}\right)+5\Gamma^{4}\right.$
$\displaystyle\left.+16\omega_{1}\omega_{2}\left(\omega_{1}^{2}+\omega_{2}^{2}\right)\right)+128m_{e}^{5}\left(\Gamma^{2}+(\omega_{1}-\omega_{2})^{2}\right)$
$\displaystyle+128m_{e}^{4}(\omega_{1}-\omega_{2})\left(\Gamma^{2}-2\omega_{1}\omega_{2}\right)$
$\displaystyle\left.+2\Gamma^{4}m_{e}\left(\Gamma^{2}+4\omega_{1}^{2}-6\omega_{1}\omega_{2}+4\omega_{2}^{2}\right)\right).$

In these equations the frequency $\omega_{2}$ of outgoing photon is given by
Eq. (11).

It is easy to check that in the limit $\Gamma\rightarrow 0$ Eq. (12) turns to
the ordinary Klein-Nishina formula for monocromatic Compton scattering cross
section

$\displaystyle\frac{d\sigma_{\mathrm{sc}}}{d\Omega_{2}}\equiv\lim\limits_{\Gamma\rightarrow
0}\frac{d\sigma_{\mathrm{sc}}^{\Gamma}}{d\Omega_{2}}=\frac{r_{e}^{2}}{2}\left(\frac{\omega_{2}}{\omega_{1}}\right)^{2}$
(15)
$\displaystyle\times\left(\frac{\omega_{2}}{\omega_{1}}+\frac{\omega_{1}}{\omega_{2}}-\sin^{2}\theta_{12}\right),$

where $r_{e}=e^{2}/(4\pi m_{e})$ is the classical electron radius.

The angular dependence of ordinary Compton cross section Eq. (16) for the
different energies of incident photon is presented in Fig. 1. These results
can be compared with the equation (12) for finite parameter $\Gamma$. Within
the considered approach which is based on the modified form of fermion
propagator given by Eq. (1), the scattering process of a quasi-monochromatic
photon with a width $\Gamma$ on stable particle is equivalent to the
scattering of a monocromatic photon by an unstable fermion with lifetime
$\tau=1/\Gamma$. As an example we consider the single Compton process for the
muon with mass $m_{\mu}\approx 207m_{e}$ and decay time
$\tau_{\mu}=2.197\times 10^{-6}$ s [15]. In Fig. 2 the relative difference
between ordinary and modified Compton cross sections is presented.

Figure 1: The dependence of differential Compton scattering cross section on
the angle $\theta_{12}$ between incident and outgoing photon (see Eq. (15)).
Here $r_{e}=e^{2}/(4\pi m_{e})$ and the photon energy is measured in units of
electron rest energy $\varepsilon_{e}\equiv\omega_{1}/m_{e}$. Figure 2:
Relative difference $\Delta\mathrm{x}_{\mathrm{rel}}$ in percent between
modified and ordinary Compton cross sections for photon scattering on free
muon (see Eqs. (12) and (15), respectively). Here
$\varepsilon_{\mu}\equiv\omega_{1}/m_{\mu}$ and parameter $\Gamma$ in Eqs.
(13), (14) is $\Gamma_{\mu}=2.996\times 10^{-10}$ eV.

It is important to consider the nonrelativistic limit of the above equations,
which corresponds to the Thomson scattering. In this case $\omega_{1}/m_{e}\ll
1$ and the energy conservation law Eq. (11) is $\omega_{1}\approx\omega_{2}$.
Then setting $\omega_{1}=\omega_{2}$ in Eq. (15) and performing the angular
integration we arrive at [6, 16]

$\displaystyle\sigma_{\mathrm{Th}}=\frac{8\pi}{3}r_{e}^{2}.$ (16)

We will be interested in a similar expression but for a finite parameter
$\Gamma$. Keeping the terms of the oder $\Gamma^{2}$ in the Eq. (15),
integrating over the angles and taking the nonreltivistic limit we have

$\displaystyle\sigma_{\mathrm{Th}}^{\mathrm{\Gamma}}=\frac{\sigma_{\mathrm{Th}}}{1+\Sigma^{\Gamma}(\omega_{1})},$
(17)

where

$\displaystyle\Sigma^{\Gamma}(\omega_{1})=2\pi\Gamma^{2}\left(70m_{e}^{2}\omega_{1}^{2}+120m_{e}^{3}\omega_{1}+60m_{e}^{4}\right.$
(18)
$\displaystyle\left.+10m_{e}\omega_{1}^{3}-13\omega_{1}^{4}\right)(15m_{e}^{2}\omega_{1}^{4})^{-1}.$

This result Eq. (17) can be compared with one obtained from the classical
electrodynamics with a radiation damping effects taking into account which is
[17, 18]

$\displaystyle\sigma_{\mathrm{Th}}=\frac{\sigma_{\mathrm{Th}}}{1+(\gamma/\omega_{1})^{2}}.$
(19)

From comparison of Eqs. (17) and (19) it follows that QED derivation leads to
a different power-law dependence of Eq. (17) on the incident radiation
frequency $\omega_{1}$. It is important to note that the finite time effects
in Thomson scattering have not been previously considered within the framework
of the QED formalism. Existing approaches to solving this problem refer only
to classical electrodynamics.

##  III Double compton scattering

In this section, we will carry out similar calculations for the double Compton
effect. To begin with, we recall a brief derivation of the scattering cross
section of the process
$\gamma(k)+e^{-}\rightarrow\gamma(k_{1})+\gamma(k_{2})+e^{-}$ [19, 12].
Following [12] the double Compton (DC) effect is described by the third order
S-matrix element

$\displaystyle S^{(3)}_{fi}=(2\pi)^{4}\delta^{(4)}(p_{f}+k_{1}+k_{2}-p_{i}-k)$
(20)
$\displaystyle\times\frac{m_{e}e^{3}}{\sqrt{8V^{5}E_{i}E_{f}\omega\omega_{1}\omega_{2}}}{\cal
M}(k,k_{1},k_{2}),$

where $\cal{M}$ is the Feynman amplitude of the process

$\displaystyle{\cal
M}(k,k_{1},k_{2})=\overline{u}(p_{f})\left\\{\hat{\varepsilon}(k)S_{F}(k-p_{f})\hat{\varepsilon}^{*}(k_{1})S_{F}(k-k_{1}-p_{f})\hat{\varepsilon}^{*}(k_{2})+(1\leftrightarrow
2)\right.$ (21)
$\displaystyle\left.+\hat{\varepsilon}^{*}(k_{1})S_{F}(-k_{1}-p_{f})\hat{\varepsilon}(k)S_{F}(k-k_{1}-p_{f})\hat{\varepsilon}^{*}(k_{2})+(1\leftrightarrow
2)\right.$
$\displaystyle\left.+\hat{\varepsilon}^{*}(k_{1})S_{F}(-k_{1}-p_{f})\hat{\varepsilon}^{*}(k_{2})S_{F}(-k_{1}-k_{2}-p_{f})\hat{\varepsilon}(k)+(1\leftrightarrow
2)\right\\}u(p_{i}).$

In the reference frame where the electron initially at rest, i.e.
$\textbf{p}_{i}=0$ and $E_{i}=m_{e}$, the energy conservation law takes the
form

$\displaystyle\omega_{2}=\frac{m_{e}(\omega-\omega_{1})-\omega\omega_{1}(1-\cos\theta_{1})}{m_{e}+\omega(1-\cos\theta_{2})-\omega_{1}(1-\cos\theta_{12})}.$
(22)

After the summation over the spin and polarizations of final states, averaging
over the spins and polarization of initial states and integrating over
$\omega_{2}$ in Eq. (20) the differential cross section is

$\displaystyle\frac{d\sigma_{\mathrm{dc}}}{d\omega_{1}d\Omega_{1}d\Omega_{2}}=\frac{m_{e}e^{6}}{2^{8}\pi^{5}}\left(\frac{\omega_{1}}{\omega}\right)\frac{m_{e}(\omega-\omega_{1})-\omega\omega_{1}(1-\cos\theta_{1})}{(m_{e}+\omega(1-\cos\theta_{2})-\omega_{1}(1-\cos\theta_{12}))^{2}}X(\omega,\omega_{1}),$
(23)

where

$\displaystyle
X(\omega,\omega_{1})=\frac{1}{4}\mathrm{Tr}\sum\limits_{\mathrm{polarizations}}|{\cal
M}(\omega,\omega_{1})|^{2}.$ (24)

Summation in Eq. (23) and algebra with Dirac matrices can be performed in a
fully analytical way with the use of FeynCalc software [20, 21]. To
demonstrate the dependence of Eq. (23) on parameter $\Gamma$ we follow the
examples considered in [19, 12]. The results of numerical evaluation of double
Compton scattering differential cross section for the energy
$\varepsilon_{e}=\omega/m_{e}=2$ (1 MeV) of incident photon and $\Gamma=0$ are
presented in Fig. 3. In this case the following infrared singularities occurs:
two for the angles $\theta_{1}=\theta_{2}=\frac{\pi}{2}$,
$\theta_{12}=\frac{\pi}{2}$ at $\varepsilon_{e}=0$ and $\varepsilon_{e}=2/3$,
and for the angles $\theta_{1}=\frac{\pi}{2}$, $\theta_{2}=0$,
$\theta_{12}=\frac{\pi}{2}$ at $\varepsilon_{e}=0$.

Figure 3: Energy spectra of double Compton scattering cross section for given
directions of both photons. The incident energy is
$\varepsilon_{e}=\omega/m_{e}=2$ (1 MeV). The solid line corresponds to the
angles $\theta_{1}=\theta_{2}=\frac{\pi}{2}$, $\theta_{12}=\frac{\pi}{2}$,
while the dashed line to $\theta_{1}=\theta_{2}=\frac{\pi}{2}$,
$\theta_{12}=0$.

To consistently cure the issue, the DC process has to be treated together with
the next to leading order (NLO) radiative corrections to the Compton process
[22, 23]. In this context, one is usually interested in the total Compton
scattering cross-section at order $\alpha^{3}$. NLO corrections display a
logarithmic divergence which can be shown to cancel with the corresponding DC
one. Employing a regularization - introducing a photon mass in the standard
approach [22] \- in the calculation of both cross-sections and summing the
two, the dependence on the regularization parameters drops out, leaving a
finite radiative correction at order $\alpha^{3}$. The drawback of this
approach is that the radiative correction cross-section now depends on the
energy resolution of the experiment, $\omega_{\mathrm{res}}$. The argument to
justify summing the two processes is that below some $\omega_{\mathrm{res}}\ll
1$ the experiment is unable to distinguish contributions from virtual photon
emission, relevant to the computation of the NLO correction, and from real
photon emission by DC. Adding all contributions to the total CS scattering
cross-section it was found that the correction can exceed the naive
$\alpha/\pi$ level at sufficiently high energies [24].

However, all of the regularization procedures are related to the consideration
of the final state of the system [22]. Alternative physical approach for
solving this problem based on the analysis of the input photon - taking into
account its finiteness in time. Indeed, the standard consideration assumes
that the incident photon is taken in the form of a monochromatic plane wave
infinite in time and space [14, 11, 25, 16]. Only in this case there is
possibility to form an arbitrarily low-frequency photon at the end of the
process. However, in reality, there are no such waves in nature. There are
only real photons, which are formed by some radiation mechanism during a
finite time interval [26]. As a result, the scattering of such a photon by an
electron lead to a finite spectrum that is completely finite in the entire
frequency range. Recently in [27] it was shown that the scattering cross
section for inverse double Compton (IDC) process have no ordinary infrared
divergence related to the resolution of photon detector in final states.
Instead of this the singularity is moved to the incident photon energies
$\omega_{1}$ and $\omega_{2}$ which are the input parameters for differential
IDC cross section. From consideration of the inverted process, it becomes
obvious that restrictions must be imposed not on the resolution of the
detector $\omega_{\mathrm{res}}$ (which hypothetically can be made arbitrarily
sensitive to the low-energy part of the DC spectrum) but on the photon source.
The latter means that the incident photon should actually be considered as
non-monochromatic one.

It is useful to compare the behaviour of DC differential cross section near
the singularity for finite parameter $\Gamma\neq 0$ and limiting case
$\Gamma=0$. The corresponding results are presented in Fig. 4 for
$\varepsilon_{e}=2$ (1 MeV) and angles $\theta_{1}=\theta_{2}=\frac{\pi}{2}$,
$\theta_{12}=0$. From Fig. 4 it is seen that in the case $\Gamma\neq 0$ the
infinite peak arising as a result of infrared divergence is greatly smoothed
out.

Figure 4: Energy spectra of double Compton scattering cross section for given
directions of both photons and different values of width parameter in units of
electron mass $y=\Gamma/m_{e}$. The incident energy is
$\varepsilon_{e}=\omega/m_{e}=2$ (1 MeV) and the angles are
$\theta_{1}=\theta_{2}=\frac{\pi}{2}$, $\theta_{12}=0$. The solid blue line
goes to infinity at $\varepsilon_{1e}=2/3$ (or $\varepsilon_{2e}=0$ according
to energy conservation law Eq. (11)).

##  IV Results and discussion

In this work, we analyzed finite-time effects in scattering one or two photons
on a free fermion particle. Localized in time interaction was introduced by a
small complex addition to the photon momentum (which means going beyond the
monochromatic plane-wave approximation [28]). To do this, we used the theory
of propagation functions of unstable fermions, which also contain a complex
addition to the momentum related to the particle’s lifetime [6, 5]. At
sufficiently low incident photon energies but still much larger than the
corresponding pulse width, the cross sections of the considered processes can
differ significantly compared to the standard approach. The expressions
obtained in this work can also be applied to the scattering of photons by
unstable particles, where the particle’s lifetime plays the role of inverse
width. The considered approach also leads to the regularization of the
infrared divergence in the cross section of the double Compton effect in
complete analogy with bound-state QED, where the natural width of atomic level
plays the role of regularising parameter.

It is worth noting that other approaches for introducing a finite interaction
time have also been discussed in a recent work [29]. In particular, in [29] it
was proposed to multiply the amplitude of the process, Eq. (21), by a
regularizing factor containing the parameter of a pulse duration from the
source.

Considering the above, it becomes evident that with an increase in the
detectors sensitivity to the low energy part of the spectrum, the cross
section of double Compton effect in the soft photon region can be studied
experimentally. In particular, direct observation of peak flattening from
infrared divergence, see Fig. 4, could verify the proposed theoretical
approach. An experimental study of the considered effects can be carried out
with optical photons. For the double Compton effect, this will result in the
reemitted soft photon being in the radio-frequency range, where sufficiently
sensitive detectors are available.

##  V Acknowledgements

The authors are grateful to Dr. V. Zalipaev (Laboratory of Quantum Processes
and Measurements, ITMO University, St. Petersburg, Russia) and Dr. D. Solovyev
(Department of Physics, St. Petersburg State University, Russia) for valuable
discussion.

## References

  * Sucher [1957] J. Sucher, S-matrix formalism for level-shift calculations, Phys. Rev. 107, 1448 (1957).
  * LÓPEZ CASTRO  _et al._ [1991] G. LÓPEZ CASTRO, J. L. M. LUCIO, and J. PESTIEAU, W± propagator on the resonance region, Modern Physics Letters A 06, 3679 (1991), https://doi.org/10.1142/S0217732391004255 .
  * CASTRO  _et al._ [1996] G. L. CASTRO, J. LUCIO M., and J. PESTIEAU, Remarks on the w propagator at the resonance, International Journal of Modern Physics A 11, 563 (1996), https://doi.org/10.1142/S0217751X96000262 .
  * Nowakowski and Pilaftsis [1993] M. Nowakowski and A. Pilaftsis, On gauge invariance of breit-wigner propagators, Zeitschrift für Physik C Particles and Fields 60, 121 (1993).
  * Kuksa [2015] V. Kuksa, Complex-mass definition and the structure of unstable particle’s propagator, Advances in High Energy Physics 2015, 490238 (2015).
  * Weinberg  _et al._ [1995] S. Weinberg, W. S, and T. de campos,  _The Quantum Theory of Fields_, Quantum Theory of Fields, Vol. 2: Modern Applications No. v. 1 (Cambridge University Press, 1995).
  * Kuksa [2014] V. I. Kuksa, Unstable states in quantum theory, Physics of Particles and Nuclei 45, 568 (2014).
  * Andreev  _et al._ [2008] O. Y. Andreev, L. N. Labzowsky, G. Plunien, and D. A. Solovyev, Qed theory of the spectral line profile and its applications to atoms and ions, Phys. Rep. 455, 135 (2008).
  * Zalialiutdinov  _et al._ [2018] T. A. Zalialiutdinov, D. A. Solovyev, L. N. Labzowsky, and G. Plunien, Qed theory of multiphoton transitions in atoms and ions, Phys. Rep. 737, 1 (2018).
  * Klein and Nishina [1929] O. Klein and Y. Nishina, Über die streuung von strahlung durch freie elektronen nach der neuen relativistischen quantendynamik von dirac, Zeitschrift für Physik 52, 853 (1929).
  * Bjorken and Drell [1964] J. D. Bjorken and S. D. Drell,  _Relativistic quantum mechanics_ (Mcgraw-Hill College, 1964).
  * Jauch and Rohrlich [1976] J. M. Jauch and F. Rohrlich,  _The theory of photons and electrons. The relativistic quantum field theory of charged particles with spin one-half_, 2nd ed., Texts and Monographs in Physics (Springer, Berlin, 1976).
  * Rapoport  _et al._ [1978] L. Rapoport, B. Zon, and N. Manakov, Theory of multiphoton processes in atoms (1978).
  * Akhiezer and Berestetskii [1965] A. I. Akhiezer and V. B. Berestetskii,  _Quantum Electrodynamics_ (Wiley-Interscience, New York, 1965).
  * Milotti [1998] E. Milotti, Energy levels and radiative transitions in muonic hydrogen, Atomic Data and Nuclear Data Tables 70, 137 (1998).
  * Berestetskii  _et al._ [1982] V. Berestetskii, E. Lifshits, and L. Pitaevskii,  _Quantum Electrodynamics_ (Oxford Butterworth-Heinemann, 1982).
  * Heinzl and Ilderton [2013] T. Heinzl and A. Ilderton, Corrections to laser electron thomson scattering, (2013).
  * Schwinger  _et al._ [1998] J. Schwinger, L. Deraad, K. Milton, W. Tsai, and J. Norton,  _Classical Electrodynamics_, Advanced book program (Avalon Publishing, 1998).
  * Mandl and Skyrme [1952] F. Mandl and T. Skyrme, The theory of the double compton effect, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 215, 497 (1952).
  * Mertig  _et al._ [1991] R. Mertig, M. Böhm, and A. Denner, Feyn calc - computer-algebraic calculation of feynman amplitudes, Computer Physics Communications 64, 345 (1991).
  * Shtabovenko  _et al._ [2016] V. Shtabovenko, R. Mertig, and F. Orellana, New developments in feyncalc 9.0, Computer Physics Communications 207, 432 (2016).
  * Brown and Feynman [1952] L. M. Brown and R. P. Feynman, Radiative corrections to compton scattering, Phys. Rev. 85, 231 (1952).
  * Ravenni and Chluba [2020] A. Ravenni and J. Chluba, The double compton process in astrophysical plasmas, Journal of Cosmology and Astroparticle Physics 2020 (10), 025.
  * Mork [1971] K. J. Mork, Radiative corrections. ii. compton effect, Phys. Rev. A 4, 917 (1971).
  * Landau and Lifshitz [1965] L. D. Landau and E. M. Lifshitz,  _Quantum Mechanics: Non-Relativistic Theory_ (Pergamon Press, 1965).
  * Dawson and Fried [1970] J. F. Dawson and Z. Fried, Plane-wave packets and their limitation in nonlinear compton scattering, Phys. Rev. D 1, 3363 (1970).
  * Dubrovich and Zalialiutdinov [2021] V. Dubrovich and T. Zalialiutdinov, Some features of the direct and inverse double-compton effect as applied to astrophysics, Physics 3, 1167 (2021).
  * Zalialiutdinov  _et al._ [2014] T. Zalialiutdinov, Y. Baukina, D. Solovyev, and L. Labzowsky, Theory of the multiphoton cascade transitions with two photon links: comparison of quantum electrodynamical and quantum mechanical approaches, Journal of Physics B: Atomic, Molecular and Optical Physics 47, 115007 (2014).
  * Lötstedt and Jentschura [2009] E. Lötstedt and U. D. Jentschura, Correlated two-photon emission by transitions of dirac-volkov states in intense laser fields: Qed predictions, Phys. Rev. A 80, 053419 (2009).


# Unusual low-temperature behavior in the half-filled band of the one-
dimensional extended Hubbard model in atomic limit

Onofre Rojas  Departamento de Física, Universidade Federal de Lavras, CP 3037,
37200-900, Lavras-MG, Brazil  S. M. de Souza  Departamento de Física,
Universidade Federal de Lavras, CP 3037, 37200-900, Lavras-MG, Brazil  J.
Torrico  Departamento de Física, Instituto de Cíencias Exatas, Universidade
Federal de Alfenas, 37133-840 Alfenas, Minas Gerais, Brazil  L. M. Verissimo
Donostia International Physics Center, Paseo Manuel de Lardizabal 4, E-20018
San Sebastián, Spain  Instituto de Física, Universidade Federal de Alagoas,
57072-970 Maceió, Alagoas, Brazil  M. S. S. Pereira  Instituto de Física,
Universidade Federal de Alagoas, 57072-970 Maceió, Alagoas, Brazil  M. L. Lyra
Instituto de Física, Universidade Federal de Alagoas, 57072-970 Maceió,
Alagoas, Brazil  Oleg Derzhko  Institute for Condensed Matter Physics,
National Academy of Sciences of Ukraine, Svientsitskii Street 1, 79011 L’viv,
Ukraine

(August 28, 2024)

###### Abstract

Recently, a kind of finite-temperature pseudo-transition was observed in
several quasi-one-dimensional models. In this work, we consider a genuine one-
dimensional extended Hubbard model in the atomic limit, influenced by an
external magnetic field and with the arbitrary number of particles controlled
by the chemical potential. The one-dimensional extended Hubbard model in the
atomic limit was initially studied in the seventies and has been investigated
over the past decades, but it still surprises us today with its fascinating
properties. We rigorously analyze its low-temperature behavior using the
transfer matrix technique and provide accurate numerical results. Our analysis
confirms that there is an anomalous behavior in the half-filled band,
specifically occurring between the alternating pair ($AP$) and paramagnetic
($PM$) phases at zero temperature. Previous investigations did not deeply
identify this anomalous behavior, maybe due to the numerical simplicity of the
model, but from analytical point of view this is not so easy to manipulate
algebraically because one needs to solve an algebraic cubic equation. In this
study, we explore this behavior and clearly distinguish the pseudo-transition,
which could easily be mistaken with a real phase transition. This anomalous
behavior mimics features of both first- and second-order phase transitions.
However, due to its nature, we cannot expect a finite-temperature phase
transition in this model.

##  I Introduction

Recent studies on various effective-one-dimensional models with short-range
interactions have revealed intriguing thermal behaviors, resembling first and
second-order phase transitions [1]. This peculiar behavior, especially around
pseudo-transitions, has been further analyzed in reference [2], focusing on
spin correlation functions. The simplest models where this unusual behavior
arises are in the Ising diamond chain [3, 4, 5]. Another similar decorated
Ising chain discussed is the Ising sawtooth-like chain model [5], along with
two and three-leg Ising ladder models [6, 7, 8]. Pseudo-transitions also
appear in other models, such as the Ising-Heisenberg diamond chain [9, 10],
and the one-dimensional double-tetrahedral model with alternating Ising spins
and delocalized electrons [11]. Similar phenomena are observed in ladder
models with Ising-Heisenberg coupling [12] and triangular tube models [13],
highlighting a pattern of pseudo-transition. Further investigations [6, 14,
15] have shown some kind of universality of pseudo-critical exponents, while
still satisfying the Rushbrooke inequality. All the above models are related
to classical and classical-quantum spin models; here by classical we mean that
the Hamiltonian does not contain non-commuting terms. However, models of other
natures, like the extended Hubbard model in a diamond chain structure [16],
and the Potts model on a diamond chain structure, have also exhibited this
unusual behavior. On the other hand, exploring pseudo-transitions in genuine
one-dimensional systems without decoration couplings is rather interesting. In
this sense, the Zimm-Bragg-Potts model was recently explored and found to
exhibit this anomalous behavior [17].

On the other hand, the Hubbard model [18] stands as a foundational model in
the modern theory of strongly correlated electrons, playing a paradigmatic
role in the study of electronic correlations in quantum materials. This model
is particularly significant in contexts where interactions are crucial. The
Bethe ansatz method helps to understand this model better, especially in
figuring out how electrons behave in different conditions and how they
interact with each other [19]. The one-dimensional extended Hubbard model is a
simplified theoretical model that describes the behavior of electrons in a
one-dimensional chain, which additionally includes nearest-neighbor
interaction energy terms. The nearest-neighbor interaction energy term
describes the Coulomb repulsion between electrons occupying neighboring sites.
The one-dimensional extended Hubbard model has been widely studied in the
literature due to its relevance in understanding the electronic properties of
one-dimensional materials such as carbon nanotubes and organic conductors
[20]. It also serves as a simple prototype model for studying strongly
correlated electron systems [21]. On the other hand, several theoretical
investigations have focused on the extended Hubbard model in the one-
dimensional case. Numerous investigations regarding the half-filled ground
states of the extended Hubbard model have been conducted. Tsuchiizu and
Furusaki [22] revisited the ground-state phase diagram of the one-dimensional
half-filled extended Hubbard model, revealing new phases and transitions using
a renormalization-group approach. Glocke et al. [23] utilized density-matrix
renormalization group methods on transfer matrices to study the thermodynamics
of the one-dimensional extended Hubbard model at half-filling. This highlights
the detection of phase transitions through standard thermodynamic measures
like isothermal compressibility and magnetic susceptibility. The study
identifies a unique phase with long-range dimer order and delineates the phase
diagram, comparing it with quantum Monte Carlo studies [24, 25]. The phase
diagram and critical exponents of the one-dimensional $U-V$ model at quarter-
filling were identified using exact diagonalization and various limit results,
identifying a transition from Luttinger liquid to charge density wave
insulator, noting dominant superconducting or spin density wave fluctuations
depending on the value of $V$ [26]. Further investigation was also performed
in reference [27].

Recent experiments have successfully demonstrated the realization of an
extended fermionic Hubbard model using a two-dimensional lattice composed of
dopant-based quantum dots. Quantum dots, often likened to artificial atoms,
can be accurately arranged into structures resembling artificial molecules and
lattices. These arrangements offer adjustable hopping amplitudes and
interaction strengths, as well as the ability to design specific point
symmetries. This advancement marks a significant step in the exploration of
complex quantum systems and could pave the way for new insights into the
behaviors of electronic systems [28]. Recent angle-resolved photoemission
spectroscopy (ARPES) studies on the one-dimensional extended Hubbard model,
employing bosonization and time-dependent calculations, reveal insights into
electron-phonon coupling and interactions in one-dimensional systems [29].
Earlier, Epstein et al. [30] discussed the metal-insulator transition of
N-methyl phenazinium (NMP) tetracyanoquinodimethane (TCNQ) based on the
strongly correlated Hubbard model ($t\ll U$). A considerable number of
theoretical studies of this model have been undertaken, and many of its
properties are now well known. However, this simple model still surprises us
with unexpected features, which are the focus in the following sections.

Although the simplest version of the extended Hubbard model has been
considerably studied and applied to various physical systems, here we consider
a typical one-dimensional extended Hubbard model in the atomic limit
(neglecting the hopping term),

$\displaystyle\boldsymbol{H}=$
$\displaystyle\sum_{i=1}^{N}\left[U\boldsymbol{n}_{i,\uparrow}\boldsymbol{n}_{i,\downarrow}+V\boldsymbol{n}_{i}\boldsymbol{n}_{i+1}\right.$
$\displaystyle\left.-\mu(\boldsymbol{n}_{i,\uparrow}+\boldsymbol{n}_{i,\downarrow})-h(\boldsymbol{n}_{i,\uparrow}-\boldsymbol{n}_{i,\downarrow})\right],$
(1)

where $U$ is the on-site Coulomb interaction, $V$ is the Coulomb interaction
between electrons on the neighboring sites, $\mu$ is the chemical potential,
and $\boldsymbol{n}_{i,\sigma}$ is the corresponding number operator at site
$i$, with spin $\sigma=\\{\uparrow,\downarrow\\}$ and
$\boldsymbol{n}_{i}=\boldsymbol{n}_{i,\uparrow}+\boldsymbol{n}_{i,\downarrow}$.
And the last term reports the contribution of external magnetic field $h$.
Despite its simplicity, the present model still astonishes us by providing
further interesting anomalous properties not previously elucidated elsewhere,
through careful analysis.

The present work is organized as follows: In Sec. 2, we give the
thermodynamics of the one-dimensional extended Hubbard model in the atomic
limit, identifying each eigenvalue of the cubic root and determining which one
is the largest, a topic not previously elucidated. In Sec. 3, we analyze a
peculiar property in the low-temperature region and explore the so-called
pseudo-critical temperature transition, along with the corresponding region
where the main properties undergo significant changes at this anomalous
temperature. This phenomenon is what we refer to as the quasi-phase diagram
[1]. In Sec. 4, we explore additional physical quantities, reporting the
influence of pseudo-transition. Finally, in Sec. 5, we present our
conclusions.

##  II Thermodynamics of the Model

In the 1970s, the thermodynamics of the one-dimensional extended Hubbard model
in the atomic limit garnered significant interest, with pioneering analyses by
Bari [31], Beni and Pincus [32], and Gallinar [33] employing the transfer
matrix approach. These studies examined specific heat, static magnetic
susceptibility, and density-density correlation functions at various
temperatures, particularly focusing on half-filled band and briefly on the
quarter-filled band case with infinite intra-atomic Coulomb repulsion. Later,
Mancini and Mancini [34, 35] advanced this research using Green’s function and
equations of motion formalism, finding four distinct phases and diverse charge
orderings at zero temperature. Their work also considered the influence of
external magnetic fields [36] on thermodynamic properties like magnetization
and specific heat, identifying critical fields associated with polarization
levels.

###  II.1 Transfer matrix

In order to express the transfer matrix of the model, we use the following
natural basis
$\\{|0\rangle,|\uparrow\rangle,|\downarrow\rangle,|\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$}\rangle\\}$.
The first state corresponds to the vacuum state, the second state denotes the
spin-up state, the third state is the spin-down state, and the fourth state
corresponds to two spins with opposite spins on same site.

In principle, this model can be solved using the transfer matrix technique
[32, 37], and the transfer matrix is given by

$\mathbf{W}=\left[\begin{array}[]{cccc}1&yw_{0,1}&y^{-1}w_{0,1}&w_{0,2}\\\
yw_{0,1}&y^{2}w_{1,1}&w_{1,1}&yw_{1,2}\\\
y^{-1}w_{0,1}&w_{1,1}&y^{-2}w_{1,1}&y^{-1}w_{1,2}\\\
w_{0,2}&yw_{1,2}&y^{-1}w_{1,2}&w_{2,2}\end{array}\right],$ (2)

where $w_{0,1}={\rm e}^{\beta\mu/2}$, $w_{0,2}={\rm e}^{\beta(\mu-U/2)}$,
$w_{1,1}={\rm e}^{\beta(\mu-V)}$, $w_{1,2}={\rm e}^{\beta(3\mu/2-2V-U/2)}$,
$w_{2,2}={\rm e}^{\beta(2\mu-4V-U)}$, $y={\rm e}^{\beta h/2}$, and $\beta$ is
the inverse temperature, $\beta=1/(k_{B}T)$.

Certainly, we can proceed by calculating
$\text{det}(\mathbf{W}-\lambda\mathbf{1})=0$ to obtain the eigenvalues of the
transfer matrix, as was previously done by Beni and Pincus [32].
Alternatively, we can manipulate the transfer matrix for convenience,
considering the spin inversion symmetry, in order to simplify our perturbative
approach calculations later on. Thus, we can employ a new set of basis vectors
$\\{|0\rangle,|\updownarrow\rangle,|\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$}\rangle,|\leftrightarrow\rangle\\}$,
where
$|\updownarrow\rangle=\left(y|\uparrow\rangle+y^{-1}|\downarrow\rangle\right)/\sqrt{z}$
corresponds to the symmetric state and
$|\leftrightarrow\rangle=\left(y|\downarrow\rangle-y^{-1}|\uparrow\rangle\right)/\sqrt{z}$
denotes the antisymmetric state, with $z=y^{2}+y^{-2}=2\cosh(\beta h)$. Note
that these states are invariant under simultaneous spin inversion and magnetic
field inversion. In this new basis, the transfer matrix (2) simply becomes:

$\mathbf{W}=\left[\begin{array}[]{cccc}1&w_{0,1}\sqrt{z}&w_{0,2}&0\\\
w_{0,1}\sqrt{z}&w_{1,1}z&w_{1,2}\sqrt{z}&0\\\
w_{0,2}&w_{1,2}\sqrt{z}&w_{2,2}&0\\\
0&0&0&0\end{array}\right]=\left[\begin{array}[]{cc}\mathbf{V}&\mathbf{0}\\\
\mathbf{0}&0\end{array}\right],$ (3)

where the matrix $\mathbf{V}$, expressed in the basis
$\\{|0\rangle,|\updownarrow\rangle,|\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$}\rangle\\}$,
results in

$\mathbf{V}=\left[\begin{array}[]{ccc}1&w_{0,1}\sqrt{z}&w_{0,2}\\\
w_{0,1}\sqrt{z}&w_{1,1}z&w_{1,2}\sqrt{z}\\\
w_{0,2}&w_{1,2}\sqrt{z}&w_{2,2}\end{array}\right].$ (4)

Obviously, the eigenvalue corresponding to the state $|\leftrightarrow\rangle$
is null.

It is straightforward to diagonalize the transfer matrix by solving the
determinant equation ${\rm det}(\mathbf{V}-\lambda\mathbf{1})=0$, which leads
to the following secular equation

$\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}=0,$ (5)

where the coefficients result in

$\displaystyle a_{0}=$ $\displaystyle
z\left(w_{1,2}^{2}+w_{0,2}^{2}w_{1,1}+w_{0,1}^{2}w_{2,2}\right)$
$\displaystyle-z\left(w_{1,1}w_{2,2}-2w_{0,2}w_{0,1}w_{1,2}\right),$
$\displaystyle a_{1}=$ $\displaystyle
z\left(w_{1,1}+w_{1,1}w_{2,2}-w_{0,1}^{2}-w_{1,2}^{2}\right)+w_{2,2}-w_{0,2}^{2},$
$\displaystyle a_{2}=$ $\displaystyle-1-zw_{1,1}-w_{2,2}.$ (6)

Therefore, the roots of the algebraic cubic equation (5) can be expressed
conveniently using trigonometric function, i.e.,

$\lambda_{j}=2\sqrt{Q}\cos\left(\tfrac{\phi-2\pi
j}{3}\right)-\frac{1}{3}a_{2},\quad j=0,1,2,$ (7)

with

$\displaystyle\phi$
$\displaystyle=\arccos\left(\tfrac{R}{\sqrt{Q^{3}}}\right),$ (8)
$\displaystyle Q$ $\displaystyle={\displaystyle\frac{a_{2}^{2}-3a_{1}}{9},}$
(9) $\displaystyle R$
$\displaystyle={\displaystyle\frac{9a_{1}a_{2}-27a_{0}-2a_{2}^{3}}{54}.}$ (10)

Note that, it is sufficient to restrict the cubic equation solution without
loss of generality to the interval of $0<\phi<\pi$; thus the eigenvalues are
ordered as follows: $\lambda_{0}>\lambda_{1}>\lambda_{2}$. This criterion is
discussed in more detail in reference [16]. To analyze the characteristics of
the eigenvalues in the interval $0<\phi<\pi$, after some trigonometric
manipulation, we have

$\displaystyle\sqrt{Q}-\tfrac{a_{2}}{3}<\lambda_{0}<2\sqrt{Q}-\tfrac{a_{2}}{3},$
(11)

obviously $\lambda_{0}$ is definitely positive, because $a_{2}<0$. Similarly,
for the second eigenvalue, we can express

$\displaystyle-\sqrt{Q}-\tfrac{a_{2}}{3}<\lambda_{1}<\sqrt{Q}-\tfrac{a_{2}}{3},$
(12)

in this case, $\lambda_{1}$ can be positive or negative depending on the
Hamiltonian parameters. For the last eigenvalue, the corresponding interval
can be expressed as

$\displaystyle-2\sqrt{Q}-\tfrac{a_{2}}{3}<\lambda_{2}<-\sqrt{Q}-\tfrac{a_{2}}{3},$
(13)

similarly, $\lambda_{2}$ can be positive or negative. Particularly to satisfy
$\lambda_{2}<0$, the following condition must be met:
$\sqrt{Q}>-\tfrac{a_{2}}{3}$, which implies that $a_{1}<0$. It is noteworthy
that for other intervals, the arrangement of the cubic root solutions merely
exchanges; for details, see the Table 2 of reference [16].

As said above, the fourth eigenvalue $\lambda_{3}$ of the transfer matrix (2)
becomes null.

###  II.2 Thermodynamic quantities and correlators

To analyze thermodynamic quantities, we use the grand partition function for a
chain consisting of $N$ site:

$\displaystyle\Xi(T,\mu,h,N)=\lambda_{0}^{N}+\lambda_{1}^{N}+\lambda_{2}^{N}.$
(14)

Given the hierarchy of the eigenvalues
($\lambda_{0}>\lambda_{1}>\lambda_{2}$), one can determine the grand potential
per site in the thermodynamic limit $N\to\infty$, which is solely dictated by
the largest eigenvalue of the transfer matrix:

$\displaystyle\Omega(T,\mu,h)=$ $\displaystyle-
k_{B}T\lim_{N\rightarrow\infty}\frac{\ln\left[\Xi(T,\mu,h,N)\right]}{N}$
$\displaystyle=$ $\displaystyle-k_{B}T\ln\lambda_{0}.$ (15)

Here, we introduce several useful thermodynamic quantities. The entropy is
calculated as ${\cal S}=-\frac{\partial\Omega}{\partial T}$; the internal
energy is given by
$\mathcal{E}=-\frac{\partial\ln\left(\lambda_{0}\right)}{\partial\beta}=k_{B}T^{2}\frac{\partial\ln\left(\lambda_{0}\right)}{\partial
T}$; analogously we can obtain specific heat at constant chemical potential
$C=T\frac{\partial{\cal S}}{\partial T}=\frac{\partial{\cal E}}{\partial T}$
[38]; the magnetization can be expressed as $M=-\frac{\partial\Omega}{\partial
h}$; magnetic susceptibility $\chi=\frac{\partial M}{\partial h}$; the
electron density $\rho$ is determined by
$\rho=-\frac{\partial\Omega}{\partial\mu}$; and the isothermal compressibility
is derived from $\kappa=\tfrac{1}{\rho^{2}}\frac{\partial\rho}{\partial\mu}$.

It is also feasible to determine other quantities using $\rho$ and $M$, which
can be expressed as
$\rho=\langle\boldsymbol{n}_{\uparrow}\rangle+\langle\boldsymbol{n}_{\downarrow}\rangle$
and
$M=\langle\boldsymbol{n}_{\uparrow}\rangle-\langle\boldsymbol{n}_{\downarrow}\rangle$,
leading to the expressions

$\displaystyle\langle\boldsymbol{n}_{\uparrow}\rangle=$
$\displaystyle\frac{\rho+M}{2}\quad\text{and}\quad$
$\displaystyle\langle\boldsymbol{n}_{\downarrow}\rangle=$
$\displaystyle\frac{\rho-M}{2}.$ (16)

Similarly, we can also obtain the following quantities

$\displaystyle\langle\boldsymbol{n}_{\uparrow}\boldsymbol{n}_{\downarrow}\rangle=$
$\displaystyle\frac{\partial\Omega}{\partial U},\quad$
$\displaystyle\langle\boldsymbol{n}_{i}\boldsymbol{n}_{i+1}\rangle=$
$\displaystyle\frac{\partial\Omega}{\partial V},$ (17)

where we define

$\langle\boldsymbol{n}_{i}\boldsymbol{n}_{i+1}\rangle\equiv\langle\left(\boldsymbol{n}_{i,\uparrow}+\boldsymbol{n}_{i,\downarrow}\right)\left(\boldsymbol{n}_{i+1,\uparrow}+\boldsymbol{n}_{i+1,\downarrow}\right)\rangle.$
(18)

Furthermore, for a transfer matrix with a non-degenerate and positively
defined spectrum, one can easily define the correlation length, $\xi$, as the
inverse logarithm of the ratio of the largest and the second-largest
eigenvalue. Although our transfer matrix eigenvalues spectrum is non-
degenerate, some eigenvalues could become negative. Therefore, we need to
compare the magnitude of each eigenvalue. For instance, focusing near the
pseudo-transition, the eigenvalues become $\lambda_{0}>0$ and, similarly,
$\lambda_{1}>0$ is also positive. Nevertheless, $\lambda_{2}<0$. However, in
terms of magnitude, we cannot determine whether $|\lambda_{1}|$ or
$|\lambda_{2}|$ is the second largest. Depending on the Hamiltonian
parameters, we may have $|\lambda_{1}|>|\lambda_{2}|$ or
$|\lambda_{1}|<|\lambda_{2}|$. Therefore, for our case, we define the
correlation length [39, 40] as follows:

$\xi=\left[\ln\left(\frac{\lambda_{0}}{\max(|\lambda_{1}|,|\lambda_{2}|)}\right)\right]^{-1}.$
(19)

Note that $\max(|\lambda_{1}|,|\lambda_{2}|)$ has nothing to do with
exchanging cubic root solutions but depends solely on the Hamiltonian
parameters.

Evidently the correlation length, $\xi$, becomes more intricate due to the
competition between the magnitudes of the second-largest eigenvalues. This
complexity in $\xi$ is more pronounced near the point where the magnitudes of
these eigenvalues compete. Far from this point of competition, the behavior of
$\xi$ could be more straightforwardly described.

##  III Phase diagrams

In what follows, we discuss the ground-state energy, focusing on peculiar
regions where anomalous behavior appears. Then we consider the low-temperature
case to discuss quasi-phase diagrams and define a pseudo-critical temperature.

###  III.1 Zero-temperature phase diagram

In the absence of magnetic field $h$, the first ground state to consider is
the frustrated phase $FR_{1}$

$\displaystyle|FR_{1}\rangle=$
$\displaystyle\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N/2}|0,\sigma_{2j}\rangle\quad\text{or}\quad\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N/2}|\sigma_{2j},0\rangle$
(20)

with the respective ground-state energy

$\displaystyle E_{FR_{1}}=$ $\displaystyle-\frac{h}{2}-\frac{\mu}{2}.$ (21)

Note that the frustrated ground-state energy becomes just at $h=0$. For a non-
null magnetic field, the state aligns with the magnetic field, thereby
becoming non-frustrated. The electron density for this phase is given by
$\rho=1/2$, commonly known as quarter-filled band. The residual entropy (per
site) of this phase is expressed as $\mathcal{S}=k_{B}\ln(2)/2$.

Another phase is the frustrated phase $FR_{2}$ given by

$\displaystyle|FR_{2}\rangle=$
$\displaystyle\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N}|\updownarrow\rangle,$
$\displaystyle=$
$\displaystyle\tfrac{1}{\sqrt{z}^{N}}\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N}\left(y|\uparrow\rangle+y^{-1}|\downarrow\rangle\right).$
(22)

The ground-state energy for this phase is

$\displaystyle E_{FR_{2}}=$ $\displaystyle V-h-\mu.$ (23)

Here again, the frustrated ground state energy becomes only when $h=0$.
However, when $h\neq 0$, the state gradually aligns with the magnetic field,
losing its frustration. The residual entropy of this phase is
$\mathcal{S}=k_{B}\ln(2)$, and the particle density is $\rho=1$ (half
occupancy).

At a half-filling band or electron density $\rho=1$, we observe the
alternation pair phase $AP$

$\displaystyle|AP\rangle=$
$\displaystyle\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N/2}|0,\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$}\rangle\quad\text{or}\quad\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N/2}|\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$},0\rangle.$
(24)

The corresponding ground-state energy is

$\displaystyle E_{AP}=$ $\displaystyle\frac{U}{2}-\mu.$ (25)

Notably, the $AP$ phase has no residual entropy and can be identified as
charge density wave ($CDW$) [23].

The third frustrated phase $FR_{3}$, with an electron density of $\rho=3/2$,
in the absence of a magnetic field, is described as

$\displaystyle|FR_{3}\rangle=$
$\displaystyle\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N/2}|\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$},\sigma_{2j}\rangle\quad\text{or}\quad\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N/2}|\sigma_{2j},\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$}\rangle.$
(26)

The ground-state energy for this phase is

$\displaystyle E_{FR_{3}}=$
$\displaystyle\frac{U}{2}+2V-\frac{h}{2}-\frac{3\mu}{2}.$ (27)

Again the ground state is frustrated at $h=0$, whereas for $h\neq 0$ the
system loses its frustration. For $h=0$, the corresponding residual entropy
leads to $\mathcal{S}=k_{B}\ln(2)/2$.

Lastly, in the fully filled phase $FF$, the ground state results in

$\displaystyle|FF\rangle=$
$\displaystyle\operatorname*{\text{\raisebox{1.07639pt}{\scalebox{0.6}{$\bigotimes$}}}}_{j=1}^{N}|\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$}\rangle.$
(28)

The corresponding ground-state energy is

$\displaystyle E_{FF}=$ $\displaystyle U-4V-2\mu.$ (29)

The electron density in this phase is $\rho=2$, and there is no residual
entropy.

It is important to note that in all the aforementioned states, the
magnetization obviously becomes $M=0$ at $h=0$. Furthermore, including a
magnetic field $h$ into the phase diagrams generally results in properties
similar to those observed in the absence of a magnetic field (not
illustrated). The main distinction is that the regions previously identified
as frustrated ($FR$) become unfrustrated in the presence of a magnetic field.
An equivalent analysis was earlier performed by Mancini and Mancini [35, 36],
offering a different perspective on the phase diagram. For example, in the
$FR_{1}$ (quarter filled) and $FR_{3}$ (three quarter filled) regions, under a
sufficiently strong magnetic field, the spin arrangements align parallel to
the external magnetic field. In contrast, the $FR_{2}$ phase, characterized by
one particle per site with randomly oriented spins, begins to show alignment
under the influence of a magnetic field. As the magnetic field strength
increases, the spins start aligning with the external field, transitioning
into a fully polarized or paramagnetic phase $PM$.

Figure 1:  Phase diagram in the $\mu/U-h/U$ plane, with solid lines, describes
the zero-temperature phase diagram for the fixed $V=0.55$. Meanwhile, the
background color density plot corresponds to electron density (left) and
magnetization (right), assuming the low temperature $T=0.03$.

For a clearer illustration, we present in Fig. 1 the phase diagram in the
$h/U-\mu/U$ plane for fixed $V=0.55$. Black solid lines describe the zero-
temperature phase diagram, while the red solid line delineates an unusual
phase boundary where anomalous behavior emerges. In the same plot, we also
incorporate a color density plot corresponding to electron density $\rho$
(left) and magnetization $M$ (right) in the low temperature region $T=0.03$.
In the left panel, we observe a phase transition between alternating pair
($AP$) and paramagnetic ($PM$) phases in the half-filled region, where a
constant density is clearly visible. In contrast, we note a distinct change in
magnetization from $M=0$ in the $AP$ phase to $M=1$ in the $PM$ phase. The
half-filling band region be our focus from now on.

Table 1: The second column provides zero-temperature phase boundary conditions, the third column gives the chemical potential at the phase boundary, while the fourth column reports the associated residual boundary entropy. Boundary | $h$ | $\mu$ | $\mathcal{S}/k_{B}$
---|---|---|---
$FF-FR_{3}$ | $\mu-U-2V$ | $h+U+2V$ | $\ln(2)$
$FR_{3}-AP$ | $4V-\mu$ | $4V-h$ | $\ln(3)/2$
$FR_{3}-FR_{2}$ | $\mu-U-4V$ | $h+U+4V$ | $\ln(1+\sqrt{3})$
$FR_{2}-AP$ | $V-\tfrac{U}{2}$ |  | $\ln(2)$
$FR_{1}-AP$ | $\mu-U$ | $\mu+U$ | $\ln(3)/2$
$FR_{1}-FR_{2}$ | $2V-\mu$ | $2V-h$ | $\ln(1+\sqrt{3})$

Table 1 presents the behavior of the phase boundaries: The second column lists
the magnetic field $h$ at each phase boundary, and the fourth column details
the residual entropy at the different boundaries. Typically, at the interface
between two phases, the residual entropy is higher than in the adjacent
phases. However, in Table 1, the phase boundary between $FR_{2}$ and $AP$
exhibits an unusual behavior, the residual boundary entropy does not exceed
the residual entropy of $FR_{2}$ [41]. This anomaly leads to unexpected
behavior at finite temperatures, as discussed in Sec. III.2.

Figure 2:  Phase diagram in the $U-V$ plane under the assumption of zero
magnetic field and $\mu=1.58$. The solid line represents the phase diagram at
zero temperature, while the background density plot illustrates the entropy at
$T=0.01$.

The phase diagram in the $U-V$ plane, assuming a zero magnetic field and
considering a suitable chemical potential ($\mu=1.58$) to enhance the
anomalous behavior at a half-filled band, is illustrated in Fig. 2. This phase
diagram features a fully filled band phase $FF$, an alternating pair phase
$AP$, and three frustrated phases. The first frustrated phase $FR_{1}$
corresponds to a quarter-filled band of electrons, with spins that can point
either up or down, an electron density of $\rho=1/2$, and a residual entropy
of $\mathcal{S}=k_{B}\ln(2)/2$. In the half-filled region, there is another
frustrated phase $FR_{2}$ with one electron per site and spins that can also
randomly point up or down, having a residual entropy of
$\mathcal{S}=k_{B}\ln(2)$. The third frustrated phase $FR_{3}$ corresponds to
a $3/4$-filled band of electrons or a quarter-filled band of holes, with an
electron density of $\rho=3/2$ and a residual entropy of
$\mathcal{S}=k_{B}\ln(2)/2$. Additionally, the diagram shows a fully-filled
band of electrons. All these curves can be obtained from Table 1, second
column, assuming $h=0$. Solid lines delineate standard phase boundaries, while
the dashed line indicates an anomalous boundary between two half-filled
regions.

###  III.2 Anomalous behavior in low-temperature region and quasi-phase
diagram

In our exhaustive analysis, we explore a unique property emerging in the low-
temperature region, a phenomenon not extensively observed or detailed in
previous studies [31, 32, 33, 35, 36]. Remarkably, our findings reveal
intriguing phenomena within such a simplistic model.

The background of Fig. 2 presents a density plot of entropy in the $V-U$ plane
at a low temperature ($T=0.01$). Our focus is on the zero-temperature boundary
between the $FR_{2}$ phase and the $AP$ phase. Given the absence of true phase
transitions, with only crossover lines occurring at finite temperatures, we
refer to the zero-temperature phases at finite temperatures as quasi-phases.
In this low-temperature region, these phases are termed the quasi-frustrated
$qFR_{2}$ region and the quasi-alternating pair $qAP$ region or quasi-$CDW$
region, where the arrangement of most spins remains similar to their
configurations at zero temperature.

Figure 3: Entropy as a function of $V$ for several values of low temperatures,
assuming fixed $U=1$ and $\mu=1.58$.

Another perspective of the background density plot illustrating the entropy,
as depicted in Fig. 2, can be seen in Fig. 3 for fixed values of $U=1$ and
$\mu=1.58$, plotted against $V$. Here, it is evident how residual entropy
influences thermal entropy. The peaks at $V\approx 0.145$ and $V\approx 0.29$
correspond to the standard interphases between $FF-FR_{3}$ (with
$\mathcal{S}=k_{B}\ln(2)$) and $FR_{3}-FR_{2}$ (with
$\mathcal{S}=k_{B}\ln(1+\sqrt{3})$), respectively. However, in the interphase
between $FR_{2}-AP$ occurring at $V=0.5$, no peaks are observed, indicating
anomalous behavior. Here, residual entropy leads to $\mathcal{S}=k_{B}\ln(2)$.

Figure 4:  (a) Entropy against $T$ for several values of magnetic fields,
assuming $U=1$ and $V=0.56$. (b) Magnetization as a function of $T$. (c)
Specific heat as a function of $T$. (d) Magnetic susceptibility as a function
of $T$. All curves correspond to the legend given in (a).

Further evidence of this anomaly is explored in Fig. 4. Panel (a) presents the
entropy (${\cal S}$) as a function of temperature ($T$) for several values of
the external magnetic field, as indicated within the panel. Here, one can
clearly observe a swift change in entropy at a specific low temperature
region. Similarly, in panel (b), we illustrate the magnetization ($M$) as a
function of temperature ($T$); the colored curves refer to the same set of
parameters as in panel (a). Once again, we observe that the magnetization for
all sets of parameters is essentially null up to a certain temperature where
an unusual feature arises, followed by a rapid change in magnetization. Almost
full polarization is achieved above the anomalous temperature with increasing
magnetic fields. For moderately high temperatures, standard thermal
magnetization governs. Another quantity we examine is the specific heat as a
function of temperature under the same conditions as in panel (a). Here, one
can clearly observe very sharp peaks, with the sharpness increasing at lower
temperatures, akin to a second-order phase transition, although it is merely a
sharp peak with no divergence. The magnetic susceptibility as a function of
temperature is depicted in panel (d), where we also observe a very sharp peak
at a certain temperature. Although it mimics a second-order phase transition,
there is no divergence, that would indicate a genuine phase transition.
Therefore, we will refer to these anomalous peaks as a “false” phase
transition or pseudo-transition at a finite temperature.

Figure 5:  Quasi-phase diagram depicting the half-filled band in the $h/U-T/U$
plane, with fixed parameters $V/U=0.54$ and $\mu/U=1.58$. The dashed curve
delineates the pseudo-critical-temperature $T_{p}$. Panels (a) and (b)
showcase density plots of magnetization and entropy, respectively.

In exploring another intriguing aspect of the low-temperature region, Fig. 5a
presents a density plot of magnetization $M$ in the $h/U-T/U$ plane, assuming
$V/U=0.54$ and $\mu/U=1.58$. Here, the dashed line represents the pseudo-
critical temperature $T_{p}$. Below a magnetic field of $h/U=0.04$ at
sufficiently low temperatures, the system resides in a quasi alternating-pair
$qAP$ region, characterized by particles predominantly arranged in alternating
pairs. As the magnetic field increases within this low temperature regime, the
system goes swiftly into a quasi-paramagnetic $qPM$ region, where most spins,
one per site, align parallel to the magnetic field. Conversely, at very weak
magnetic fields and increasing temperatures, the $FR_{2}$ frustrated system
leads into a quasi-frustrated $qFR_{2}$ region. The boundary between $qFR_{2}$
and $qPM$ is marked by a standard crossover, or as visualized in the
magnetization plot, a diffuse boundary, appearing as a fuzzy region. In panel
(a), one can observe the magnetization shift between $qAP$ and $qPM$ regions,
although the boundary between $qAP$ and $qFR_{2}$ remains indistinct due to
null magnetization in both regions. Similarly, Fig. 5b displays the entropy
$\mathcal{S}$, using the same parameters as panel (a). Here, the crossover
boundary, denoted by a dashed line, outlines the pseudo-critical temperature
as a function of magnetic field. In contrast, panel (b) shows that the
boundary between $qAP$ and $qPM$ becomes indistinguishable with almost null
entropy, while the boundary between $qAP$ and $qFR_{2}$ is sharply defined.
This distinction is due to the residual entropy of $FR_{2}$ being
$\mathcal{S}=k_{B}\ln(2)$. However, at temperatures above $T>0.06$, the
entropy significantly increases, becoming clearly larger than
$\mathcal{S}\gtrsim 0.7k_{B}$.

###  III.3 Pseudo-critical temperature condition

This detailed analysis further explores the interesting behavior in the low-
temperature region of the model. Previously, the pseudo-critical temperature
was defined for models where the free energy involved a simple square root
expression [1, 2, 14]. Defining the pseudo-critical temperature generally
poses a challenge, as the only available indicator is an anomalous behavior in
the low-temperature region, marked by significant changes in entropy and
magnetization, along with sharp peaks in correlation length, specific heat,
and magnetic susceptibility. At very low temperatures, these peaks occur at
roughly the same temperature. However, as the anomalous behavior shifts to
higher temperatures $T_{p}$, these peaks appear at slightly different
temperatures for each quantity, complicating the definition of the pseudo-
critical temperature.

An interesting observation is that this anomalous behavior is evident in
fundamental quantities such as the ratio between the second largest and the
largest eigenvalues. To analyze this more deeply, consider the temperature and
magnetic field dependence of the eigenvalues. At a fixed magnetic field
$h_{1}$, there exists a special temperature, the pseudo-critical temperature
$T_{p}$, which satisfies the relation

$\left.\frac{\partial}{\partial
T}\left(\tfrac{\lambda_{1}\left(T,h_{1}\right)}{\lambda_{0}\left(T,h_{1}\right)}\right)\right|_{T_{p}}=0.$
(30)

A similar analysis can be conducted by taking derivatives with respect to the
magnetic field of the ratio of the second largest to the largest eigenvalues,
assuming a fixed temperature $T_{p^{*}}$:

$\left.\frac{\partial}{\partial
h}\left(\tfrac{\lambda_{1}\left(T_{p^{*}},h\right)}{\lambda_{0}\left(T_{p^{*}},h\right)}\right)\right|_{h_{1}}=0.$
(31)

Comparing $T_{p^{*}}$ and $T_{p}$, we observe that this quantity becomes
slightly different as the temperature increases.

In order to explore further properties of the eigenvalues given in (7), let us
define $\hat{\lambda}_{j}(T,h)=\lambda_{j}{\rm
e}^{\varepsilon_{0}(h)/k_{B}T}$, where $\varepsilon_{0}(h)$ represents the
ground-state energy. This formulation, as discussed in reference [42],
primarily aims to handle the ground-state energy explicitly, thereby
circumventing the issue of dealing with extremely large numbers. Thus, Fig. 6a
depicts the three eigenvalues $\hat{\lambda}_{j}$ as functions of temperature,
under the conditions of a null magnetic field, $\mu=1.58$, $U=1$, and
$V=0.55$. Below the temperature $T_{1}$, $\hat{\lambda}_{2}$ approaches
$-\hat{\lambda}_{0}$, whereas above $T_{1}$, $\hat{\lambda}_{2}$ approaches
$-\hat{\lambda}_{1}$ but quickly diverges due to thermal fluctuations.

Following this, we introduce the function

$g(T,h)=\hat{\lambda}_{0}(T,h)+\hat{\lambda}_{1}(T,h)+2\hat{\lambda}_{2}(T,h).$
(32)

This function, $g(T,h)$, can yield either positive or negative values.
Therefore, the condition

$g(T_{1},h_{1})=0,$ (33)

for a fixed $h_{1}$ leads to the identification of a specific temperature,
$T_{1}$. This approach provides a practical method for examining the
temperature-dependent behavior of the system under study.

Panel (b) of the Fig. 6 illustrates the function $g(T,h)$ as it varies with
temperature, under the conditions of a null magnetic field, $\mu=1.58$, $U=1$,
and $V=0.55$. In this representation, one can observe that the function
$g(T,h)$ passes through zero at a specific temperature, denoted as $T_{1}$.
This provides a clear visual indication of the pseudo-critical point in the
behavior of the system.

Alternatively, the condition (33) can be simplified through some algebraic
manipulation, even in more general cases and with generic coefficients of a
cubic polynomial as outlined in (5). This leads to a simple yet interesting
result:

$2\tilde{a}_{2}^{3}+\tilde{a}_{1}\tilde{a}_{2}+\tilde{a}_{0}=0,$ (34)

where $\tilde{a}_{j}$ is defined as $a_{j}(T_{1},h_{1})$.

Furthermore, we can use the expressions defined in (6) within the context
provided by (34) to obtain after some algebraic manipulation the following
condition

$z^{2}w_{1,1}^{2}+3zw_{1,1}-w_{0,2}^{2}+2=0.$ (35)

Given that $zw_{1,1}\sim w_{0,2}\gg 1$ the dominant terms are
$z^{2}w_{1,1}^{2}$ and $w_{0,2}^{2}$, while linear term $3zw_{1,1}$ and $2$
become negligible and are thus omitted for the simplification. Consequently,
the above expression reduces simply to:

$w_{1,1}z-w_{0,2}=0.$ (36)

Figure 6:  (a) Cubic root solutions $\hat{\lambda}_{j}$ as a function of
temperature, under the condition of a null magnetic field ($h=0$), with fixed
values of chemical potential ($\mu=1.58$), on-site interaction strength
($U=1$), and nearest-neighbor interaction strength ($V=0.55$). (b) Plot of
$g(T,0)$ against temperature $T$, using the same parameter set as in panel
(a).

Furthermore, by employing a perturbative approach, as elaborated in Appendix A
and illustrated through Eqs. (41) and (42) for the unperturbed result, we
reaffirm the critical equation derived earlier, simplifying our findings into
this coherent expression.

In regions of low temperature, $T_{p}\approx T_{p^{*}}\approx T_{1}$, which
are very close to each other. When this occurs, we can observe anomalous
behavior, and we simply define it as the pseudo-critical temperature, $T_{p}$.

For the present model the condition (36) provides a sufficiently accurate
condition for determining the pseudo-critical temperature, where $T_{p}\approx
T_{1}$ holds. Specifically, the pseudo-critical temperature for a fixed
magnetic field $h_{1}$ can be represented using the results $w_{0,2}={\rm
e}^{\beta(\mu-U/2)}$, $w_{1,1}={\rm e}^{\beta(\mu-V)}$, and $z=2\cosh(\beta
h)$, leading to the condition

$2\cosh(\beta_{1}h_{1})={\rm e}^{\beta_{1}(V-U/2)}.$ (37)

Assuming $h_{1}>0$ and at very low-temperatures, we deduce that
$h_{1}=V-\frac{U}{2}$. Under these conditions, a pseudo-transition is observed
at low temperatures.

Eventually, to refine our previous result, we can apply a first-order
perturbation approximation with respect to the non-half-filled contribution,
as detailed in Appendix A. By using the results given in Eqs. (47) and (48) we
arrive at the modified condition:

$w_{1,1}z-w_{0,2}+\frac{3}{2}(1+w_{2,2})=0.$ (38)

Moreover, further corrections can be applied up to a second-order
approximation, as detailed in the Appendix A.

Figure 7:  Plot of $T_{1}$ (thin black line), $T_{p^{*}}$ (thick green line),
and $T_{p}$ (thin red line) as functions of $h$. The curves are obtained under
the assumption of fixed parameters $\mu=1.58$, $U=1$, and varying $V$ values:
$V=0.54$ (inner curve), $V=0.56$ (middle curve), and $V=0.58$ (outer curve).

In Fig. 7, the temperatures $T_{1}$ (represented by a thin black line),
$T_{p^{*}}$ (depicted by a thick green line), and $T_{p}$ (illustrated with a
thin red line) are shown as functions of $h$, under fixed parameters
$\mu=1.58$, $U=1$, and varying $V$ values: $0.54$ (inner curve), $0.56$
(middle curve), and $0.58$ (outer curve). For the inner curve, $T_{1}$,
$T_{p^{*}}$, and $T_{p}$ are visually almost overlapping. However, as we
observe the middle curve, the discrepancy between $T_{1}\approx T_{p^{*}}$ and
$T_{p}$ becomes more pronounced at higher temperatures. This discrepancy is
even more significant for the outer curve. Based on this discrepancy, we can
conclude that the pseudo-critical transition is marked for temperatures below
$T_{p}\lesssim 0.06$, where all curves overlap (when $T_{1}\approx T_{p^{*}}$
and $T_{p}$ start to diverge). A practical approach to observe the pseudo-
critical temperature at higher temperatures is to consider when $\Delta
T_{p}=T_{p}-T_{p^{*}}\rightarrow 0$. This is because when $\Delta T_{p}$
becomes significant, it becomes less meaningful to refer to it as a pseudo-
critical temperature.

Table 2:  Pseudo-critical temperature for $\mu=1.58$, $U=1$, and $V=0.56$. The second column is computed using equation (30), the third column utilizes equation (31), and the fourth column is determined by employing equation (33). $h$ | $T_{p}$ | $T_{p^{*}}$ | $T_{1}$
---|---|---|---
$0.020$ | $0.0789151013282$ | $0.0831149115123$ | $0.0831316467916$
$0.040$ | $0.0688806941568$ | $0.0711199370000$ | $0.0711277471472$
$0.050$ | $0.0582435091646$ | $0.0591242044255$ | $0.0591267047363$
$0.055$ | $0.0489342571491$ | $0.0492080726555$ | $0.0492080726555$
$0.058$ | $0.0392459736812$ | $0.0392913669397$ | $0.03929143951989$
$0.059$ | $0.0336811084457$ | $0.0336914844331$ | $0.03369149663655$
$0.0599$ | $0.0221859758365$ | $0.0221860334134$ | $0.02218603343118$

##  IV Further Physical Quantities

Overall, physical quantities such as magnetic susceptibility or specific heat
play critical roles in thermodynamics, which are essential for understanding
the behavior of the Hubbard model. Therefore, we explore these physical
quantities, focusing particularly on the anomalous properties they exhibit.

Figure 8:  (a) Internal energy plotted against $T/T_{p}$ for various $T_{p}$
values and magnetic fields as specified in Table 2, assuming $U=1$ and
$V=0.56$. (b) Entropy shown as a function of $T/T_{p}$. (c) Magnetization
represented as a function of $T/T_{p}$. (d) Electron density depicted as a
function of $T/T_{p}$.

To provide a clearer depiction of this anomaly, we divide the temperature $T$
by the pseudo-critical temperature $T_{p}$ listed in Table 2, where the
anomaly manifests. Thus in Fig. 8a, the internal energy $\mathcal{E}$ is
presented as a function of temperature normalized by the pseudo-critical
temperature ($T/T_{p}$) for various fixed magnetic fields, as described within
the panel. For each fixed magnetic field, a corresponding pseudo-critical
temperature $T_{p}$ is given in Table 2. At a magnetic field of $h=0.0599$,
the internal energy $\mathcal{E}$ remains almost constant, but there is a
noticeable “jump” at $T/T_{p}=1$. This indicates a change from the quasi-
alternating pair $qAP$ region to a one-electron per site region (almost
polarized region) or simply a quasi-paramagnetic $qPM$ region. As the external
magnetic field decreases, this sharp boundary becomes more gradual, although
the “jump” in internal energy becomes larger. Similarly, in panel (b), the
entropy is observed under the same conditions as in (a), highlighting these
features. Notably, for a magnetic field $h\lesssim 0.06$, the boundary between
$qAP$ and $qPM$ regions becomes more distinct. As $h$ approaches $0.06$, this
characteristic resembles a first-order or discontinuous phase transition,
despite the absence of actual discontinuity. Panel (c) depicts the
magnetization $M$ as a function of $T/T_{p}$. This allows for the
corroboration of the spin orientation in the $qAP$ region, where the
magnetization is nearly null, while in the $qPM$ region, the magnetization is
almost saturated. This feature is more pronounced at weaker magnetic fields,
becoming smoother as the magnetic field decreases. Finally, panel (d)
illustrates the electron density $\rho$ as a function of $T/T_{p}$, using the
same set of magnetic fields considered in panel (a). Here, it is observed that
the electron density remains almost constant $\rho\approx 1$ (but always
smaller than $1$, $\rho<1$) for temperatures $T/T_{p}<1$, corresponding to a
half-filled band. At $T/T_{p}=1$, there is a small but distinct depression in
electron density, which is sharper and more pronounced at higher magnetic
fields. As expected, for lower magnetic fields, the electron density decreases
with increasing temperature.

Figure 9:  (a) Correlation length plotted against $T/T_{p}$, with $T_{p}$
values provided in Table 2, assuming $U=1$ and $V=0.56$. Dashed lines describe
the function $1/\ln(\lambda_{0}/\lambda_{1})$. (b) Specific heat variation
with $T/T_{p}$. (c) Magnetic susceptibility versus $T/T_{p}$. (d) Isothermal
electron compressibility as a function of $T/T_{p}$.

In Fig. 9a, the correlation length $\xi$ is illustrated as a function of
$T/T_{p}$ for a fixed external magnetic field, as specified inside the panel.
It is important to note that for each fixed magnetic field, there is a
corresponding pseudo-critical temperature $T_{p}$, as listed in Table 2. The
correlation length, $\xi$ (defined in Sec. IIB, Eq. (19), and represented by
solid lines), demonstrates a decreasing function with an inflection point
around $T/T_{p}=1$. As the magnetic field decreases, the inflection curvature
vanishes. Contrary to what has been previously reported in the literature [1,
14, 16, 42, 41], $\xi$ does not show the typical peak observed around the
pseudo-transition [1, 14, 16, 42, 41]. However, the function
$1/\ln(\lambda_{0}/\lambda_{1})$ (depicted by dashed lines) exhibits a sharp
peak as it approaches $h=0.06$ at $T/T_{p}=1$. The peak of this function
mimics a typical peak of $\xi$ around the pseudo-transition found in
references [1, 14, 16, 42, 41], indicative of a swift change between $qAP$ and
$qPM$ regions. However, as the magnetic field is decreased, this sharp peak
becomes more gradual. Similarly, panel (b) also illustrates the specific heat
$C$ as a function of $T/T_{p}$, using the same set of magnetic fields. Here,
we can observe a sharp peak, akin to that seen in a second-order phase
transition, although there is no actual divergence. Furthermore, panel (c)
displays the magnetic susceptibility plotted against $T/T_{p}$, clearly
showing the sharp peak at $T/T_{p}=1$, which mimics a second-order phase
transition. Finally, panel (d) demonstrates the behavior of isothermal
compressibility $\kappa$. As $h$ approaches $0.06$, $\kappa$ diminishes,
indicating that the system becomes less compressible under these conditions.
Conversely, for smaller magnetic fields, $\kappa$ increases, suggesting that
the system is more easily compressible. At relatively low magnetic fields,
$\kappa$ exhibits a minor peak with a maximum at $T/T_{p}=1$, though this peak
becomes less pronounced with lower magnetic fields.

It is worth mentioning that such a quantity as the entropy always tends to
increase as $T$ grows in accordance with the second law of thermodynamics,
which indicates the thermal stability of the system. Similarly, the positivity
of specific heat and isothermal compressibility around anomalous behavior also
serves as an indicator of stability.

Figure 10:  (a) Electron density for spin-up,
$\langle\boldsymbol{n}_{\uparrow}\rangle$, as a function of temperature
$T/T_{p}$ across various weak magnetic fields, assuming $U=1$ and $V=0.56$.
(b) Electron density for spin-down,
$\langle\boldsymbol{n}_{\downarrow}\rangle$, as a function of temperature
$T/T_{p}$, corresponding to the magnetic fields used in (a). (c) On-site
average, $\langle\boldsymbol{n}_{\uparrow}\boldsymbol{n}_{\downarrow}\rangle$,
as a function of temperature $T/T_{p}$, assuming the magnetic field values
from (a). (d) On-site correlation function,
$c(\boldsymbol{n}_{\uparrow},\boldsymbol{n}_{\downarrow})$, as a function of
temperature $T/T_{p}$, for the same set of parameters as in (a).

In Fig. 10a, the spin-up electron density
$\langle\boldsymbol{n}_{\uparrow}\rangle$ is illustrated as a function of
temperature $T/T_{p}$, for several weak magnetic fields. Notably, in the half-
filled band condition, it is clear how the spin-up averages are organized: for
$T/T_{p}<1$, the spin-up average alternates, given by
$\langle\boldsymbol{n}_{\uparrow}\rangle\approx 0.5$, with the remaining spins
arranged as spin-down. For $T/T_{p}>1$, most spins are upwardly arranged. The
higher the magnetic field (but $h<0.6$), the stronger the curvature change at
$T/T_{p}=1$. Similarly, in panel (b), the spin-down electron density
$\langle\boldsymbol{n}_{\downarrow}\rangle$ is reported as a function of
temperature $T/T_{p}$, for the same set of magnetic fields as in panel (a).
Here, for $T/T_{p}<1$, almost half of the spins are arranged with spin-down on
average, while the remaining are spin-up, as illustrated in panel (a). In
contrast, for $T/T_{p}>1$, there are nearly no spin-up electrons, as most
spins align with the external magnetic field. Another interesting quantity is
the on-site electron average
$\langle\boldsymbol{n}_{\uparrow}\boldsymbol{n}_{\downarrow}\rangle$ as a
function of temperature $T/T_{p}$, assuming the same magnetic field values as
in panel (a). This confirms the electron arrangement in the low-temperature
region: For $T/T_{p}<1$, the average leads to
$\langle\boldsymbol{n}_{\uparrow}\boldsymbol{n}_{\downarrow}\rangle\approx
0.5$, indicating that most spins are arranged in pairs alternating with empty
sites. For $T/T_{p}>1$,
$\langle\boldsymbol{n}_{\uparrow}\boldsymbol{n}_{\downarrow}\rangle\rightarrow
0$, as most spins are parallel and aligned with the external magnetic field.
For a very weak magnetic field, the curvature change is evident, and it
diminishes as the magnetic field increases. Similarly, in panel (d), the
correlation function
$c(\boldsymbol{n}_{\uparrow},\boldsymbol{n}_{\downarrow})=\langle\boldsymbol{n}_{\uparrow}\boldsymbol{n}_{\downarrow}\rangle-\langle\boldsymbol{n}_{\uparrow}\rangle\langle\boldsymbol{n}_{\downarrow}\rangle$
is depicted as a function of temperature $T/T_{p}$, again for the same set of
parameters as in panel (a). For $T/T_{p}<1$, the correlation function is
roughly $c(\boldsymbol{n}_{\uparrow},\boldsymbol{n}_{\downarrow})=0.25$, while
for $T/T_{p}>1$, it becomes negative, as
$\langle\boldsymbol{n}_{\uparrow}\boldsymbol{n}_{\downarrow}\rangle\approx 0$,
$\langle\boldsymbol{n}_{\uparrow}\rangle\approx 0.5$ and
$\langle\boldsymbol{n}_{\downarrow}\rangle\approx 0.5$, leading to a negative
correlation function of approximately $\approx-0.25$.

##  V Conclusion

Our study has provided a comprehensive re-examination of the one-dimensional
extended Hubbard model in the atomic limit, a subject initially explored in
the 1970s. This result not only confirms previous outcomes of the model but
also reveals novel insights, especially the occurrence of finite-temperature
pseudo-transitions in specific quasi-one-dimensional models.

Initially, our analysis focused on delineating the zero-temperature phase
diagram at zero magnetic field. Here, we identified three distinct types of
frustrated phases, enriching our understanding of the intricate phase
structure within the model. This phase characterization was facilitated by the
employment of the transfer matrix technique, followed by precise numerical
analyses, which unveiled anomalous behaviors within the half-filled band of
the model. These behaviors are particularly pronounced in the low-temperature
regime, especially during the gradual shift between alternating pair ($AP$)
and paramagnetic ($PM$) phases. The subtlety of these anomalous behaviors,
possibly overlooked in earlier studies due to the deceptively simple numerical
facade of the model, but from an algebraic perspective, requiring the solution
of an algebraic cubic equation introduces a level of slightly more elaborate
manipulation, revealing the intricate nature of the model.

Further investigations explored the low-temperature pseudo-phase transitions
at the half-filling band, where the aforementioned anomalous behavior is
predominantly observed. We observed that close to this anomalous region, a
pseudo-transition emerges that exhibits characteristics reminiscent of both
first- and second-order phase transitions. This pseudo-transition is
particularly fascinating, as it mimics the intricacies of real phase
transitions, adding a layer of sophistication to our understanding of the
model. Moreover, our exploration extended to diverse properties such as
electron density, on-site correlations, and nearest-site electron averages. In
examining these quantities, we discerned the underlying reasons for the
occurrence of pseudo-transitions, thereby deepening our grasp of the model
behavior. The intrinsic nature of the model, defined by its nearest-neighbor
interaction, precludes the occurrence of real phase transitions at finite
temperatures, thereby accentuating the distinctive properties and behaviors
inherent to the model within a low-dimensional system and underscoring its
peculiarity.

The property we have explored may not necessarily be exclusive to one-
dimensional models; however, uncovering this anomalous residual boundary
entropy could become a challenging task in higher dimensional systems.

###### Acknowledgements.

O. R. and S. M. S. thank CNPq and FAPEMIG for partial financial support. M. L.
L. and M. S. S. Pereira also acknowledge the financial support of FAPEAL. O.
D. acknowledges support through the EURIZON project (project #3025), which is
funded by the European Union under grant agreement No.871072.

##  Appendix A Transfer matrix around pseudo-transition

The phenomenon of pseudo-transition arises near the alternating pairs $AP$ and
paramagnetic $PM$ phases. Although the solution of (7) is exact and analytic,
the analytical expression for finding the pseudo-critical temperature could
becomes a cumbersome task due to involving cubic root expression. An
alternative approach to determine this condition is to split the transfer
matrix into two terms, namely,
$\mathbf{V}=\mathbf{V}_{0}+\varsigma\mathbf{V}_{1}$. Here, $\mathbf{V}_{0}$
represents the core structure matrix, whose matrix elements just includes the
half-filled case, while $\varsigma\mathbf{V}_{1}$ describes a small
perturbation applied to the matrix, corresponding to non-half-filled band. The
detailed formulations of these matrices are provided as follows:

$\mathbf{V}_{0}=\left[\begin{array}[]{ccc}0&0&w_{0,2}\\\ 0&w_{1,1}z&0\\\
w_{0,2}&0&0\end{array}\right]$ (39)

and

$\mathbf{V}_{1}=\left[\begin{array}[]{ccc}1&w_{0,1}\sqrt{z}&0\\\
w_{0,1}\sqrt{z}&0&w_{1,2}\sqrt{z}\\\
0&w_{1,2}\sqrt{z}&w_{2,2}\end{array}\right],$ (40)

respectively. In this context, $\varsigma$ acts as a small formal parameter
that, for our purposes, is set to 1 to measure how the system deviates from
the half-filled-band limit. This approach allows us to analyze the effects of
slight perturbations on the system by examining the changes introduced by
$\mathbf{V}_{1}$ in relation to the original matrix $\mathbf{V}_{0}$.

Therefore, diagonalizing $\mathbf{V}_{0}$, we have the following eigenvalues

$\displaystyle v_{0}^{(0)}=$ $\displaystyle w_{1,1}z,$ (41) $\displaystyle
v_{1}^{(0)}=$ $\displaystyle w_{0,2},$ (42) $\displaystyle v_{2}^{(0)}=$
$\displaystyle-w_{0,2}.$ (43)

From this solution one can clearly verify that $v_{1}^{(0)}>v_{2}^{(0)}$,
however, we cannot affirm any condition between $v_{0}^{(0)}$ and
$v_{1}^{(0)}$, since there is no restriction, so one can establish the
following condition $v_{1}^{(0)}=v_{0}^{(0)}$ this would be essential to find
the pseudo-critical temperature, as discussed in Sec. III.3. The corresponding
eigenvectors can be expressed by

$\displaystyle|u_{0}^{(0)}\rangle=$ $\displaystyle|\updownarrow\rangle,$ (44)
$\displaystyle|u_{1}^{(0)}\rangle=$
$\displaystyle\tfrac{1}{\sqrt{2}}\left(|0\rangle+|\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$}\rangle\right),$
(45) $\displaystyle|u_{2}^{(0)}\rangle=$
$\displaystyle\tfrac{1}{\sqrt{2}}\left(-|0\rangle+|\rotatebox[origin={c}]{90.0}{$\leftrightharpoons$}\rangle\right).$
(46)

Now let us improve our previous result, thus the first-order correction on the
transfer matrix eigenvalues can be obtained perturbatively, resulting in

$\displaystyle v_{0}^{(1)}=$ $\displaystyle 0,$ (47) $\displaystyle
v_{1}^{(1)}=$ $\displaystyle\tfrac{1}{2}\left(1+w_{2,2}\right),$ (48)
$\displaystyle v_{2}^{(1)}=$
$\displaystyle\tfrac{1}{2}\left(1+w_{2,2}\right).$ (49)

Here evidently we have the following condition
$v_{0}^{(1)}<v_{1}^{(1)}=v_{2}^{(1)}$, because $w_{2,2}>0$.

Similarly, the second-order corrections of eigenvalues becomes,

$\displaystyle v_{0}^{(2)}=$
$\displaystyle\tfrac{z\left(w_{0,1}+w_{1,2}\right)^{2}}{2\left[w_{1,1}z-w_{0,2}\right]}+\tfrac{z\left(w_{0,1}-w_{1,2}\right)^{2}}{2\left[w_{1,1}z+w_{0,2}\right]},$
(50) $\displaystyle v_{1}^{(2)}=$
$\displaystyle\tfrac{1}{8}\tfrac{\left(w_{2,2}-1\right)^{2}}{w_{0,2}}-\tfrac{z\left(w_{0,1}+w_{1,2}\right)^{2}}{2\left[w_{1,1}z-w_{0,2}\right]},$
(51) $\displaystyle v_{2}^{(2)}=$
$\displaystyle-\tfrac{1}{8}\tfrac{\left(w_{2,2}-1\right)^{2}}{w_{0,2}}-\tfrac{z\left(w_{0,1}-w_{1,2}\right)^{2}}{2\left[w_{1,1}z+w_{0,2}\right]}.$
(52)

It is evident that $u_{1}^{(2)}>u_{2}^{(2)}$ and $u_{2}^{(2)}<0$.

Therefore, the corresponding eigenvalues up to order ${\cal
O}(\varsigma^{3})$, are given approximately as follow

$\displaystyle\lambda_{j}=$ $\displaystyle v_{j}^{(0)}+\varsigma
v_{j}^{(1)}+\varsigma^{2}v_{j}^{(2)}+{\cal O}(\varsigma^{3}).$ (53)

It is evident that the eigenvalues up to second-order correction indicate that
$\lambda_{1}>\lambda_{2}$. However, we cannot affirm just by looking the
perturbative eigenvalues which eigenvalues is the largest one $\lambda_{0}$ or
$\lambda_{1}$. Although we have confirmed at the end of Sec. II.1 that
$\lambda_{0}$ must be the largest one. Nevertheless, our perturbative result
could be useful to find the condition of pseudo-critical temperature up to
second-order approximation. Therefore, by using the condition given in (33) we
can find the following relation

$w_{1,1}z-w_{0,2}+\tfrac{3\left(1+w_{2,2}\right)}{2}-\tfrac{1}{8}\tfrac{\left(w_{2,2}-1\right)^{2}}{w_{0,2}}-\tfrac{z\left(w_{0,1}-w_{1,2}\right)^{2}}{2\left[w_{1,1}z+w_{0,2}\right]}=0,$
(54)

this result would be relevant when we look for a more accurate value of the
pseudo-critical temperature.

## References

  * [1] S. M. de Souza and O. Rojas, Solid State Commun. 269, 131 (2018)
  * [2] I. M. Carvalho, J. Torrico, S. M. de Souza, O. Rojas, and O. Derzhko, Ann. Phys. 402, 45 (2019)
  * [3] J. Strečka, Acta Phys. Pol. A 137, 610 (2020)
  * [4] W. Yin, https://doi.org/10.48550/arXiv.2312.11722
  * [5] W. Yin, https://doi.org/10.48550/arXiv.2401.00948
  * [6] T. Hutak, T. Krokhmalskii, O. Rojas, S. M. de Souza, and O. Derzhko, Phys. Lett. A 387, 127020 (2021)
  * [7] W. Yin, Phys. Rev. Research 6, 013331 (2024)
  * [8] W. Yin, https://doi.org/10.48550/arXiv.2006.15087
  * [9] J. Torrico, M. Rojas, S. M. de Souza, O. Rojas, and N. S. Ananikian, EPL 108, 50007 (2014)
  * [10] J. Torrico, M. Rojas, S. M. de Souza, and O. Rojas, Phys. Lett. A 380, 3655 (2016)
  * [11] L. Galisova and J. Strečka, Phys. Rev. E 91, 022134 (2015)
  * [12] O. Rojas, J. Strečka, and S. M. de Souza, Solid State Commun. 246, 68 (2016)
  * [13] J. Strečka, R. C. Alecio, M. Lyra, and O. Rojas, J. Magn. Magn. Mater. 409, 124 (2016)
  * [14] O. Rojas, J. Strecka, M. L. Lyra, and S. M. de Souza, Phys. Rev. E 99, 042117 (2019)
  * [15] T. Krokhmalskii, T. Hutak, O. Rojas, S. M. de Souza, and O. Derzhko, Physica A 573, 125986 (2021)
  * [16] O. Rojas, S. M. de Souza, J. Torrico, L. M. Verissimo, M. S. S. Pereira, and M. L. Lyra, Phys. Rev. E 103, 042123 (2021)
  * [17] Y. Panov and O. Rojas, Phys. Rev. E 103, 062107 (2021)
  * [18] J. Hubbard, Proc. R. Soc. Lond. A 276, 238 (1963)
  * [19] F. H. L. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin, The One-Dimensional Hubbard Model (Cambridge University Press, 2009)
  * [20] Z. Wang, D. Psiachos, R. F. Badilla, and S. Mazumdar, J. Phys.: Condens. Matter 21, 095009 (2009)
  * [21] E. Dagotto and Y. Tokura, MRS Bulletin 33, 1037 (2008)
  * [22] M. Tsuchiizu and A. Furusaki, Phys. Rev. B 69, 035103 (2004)
  * [23] S. Glocke, A. Klümper, and J. Sirker, Phys. Rev. B 76, 155121 (2007)
  * [24] P. Sengupta, A. W. Sandvik, and D. K. Campbell, Phys. Rev. B 65, 155113 (2002)
  * [25] A. W. Sandvik, L. Balents, and D. K. Campbell, Phys. Rev. Lett. 92, 236401 (2004)
  * [26] F. Mila and X. Zotos, EPL 24, 133 (1993)
  * [27] H. Ding and J. Zhang, Chinese J. Phys. 54, 237 (2016)
  * [28] X. Wang, E. Khatami, F. Fei, J. Wyrick, P. Namboodiri, R. Kashid, A. F. Rigosi, G. Bryant, and R. Silver, Nat. Commun. 13, 6824 (2022)
  * [29] H.-X. Wang, Y.-M. Wu, Y.-F. Jiang, and H. Yao, Phys. Rev. B 109, 045102 (2024)
  * [30] A. J. Epstein, S. Etemad, A. F. Garito, and A. J. Heeger, Solid State Commun. 9, 1803 (1971)
  * [31] R. A. Bari, Phys. Rev. B 3, 2662 (1971)
  * [32] G. Beni and P. Pincus, Phys. Rev. B 9, 2963 (1974)
  * [33] J. P. Gallinar, Phys. Rev. B 11, 4421 (1975)
  * [34] F. Mancini, Eur. Phys. J. B 47, 527 (2005)
  * [35] F. Mancini and F. P. Mancini, Phys. Rev. E 77, 061120 (2008)
  * [36] F. Mancini and F. P. Mancini, Eur. Phys. J. B 68, 341 (2009)
  * [37] R. J. Baxter, Exactly Solved Models in Statistical Mechanics (Academic Press, New York, 1982)
  * [38] V. Derzhko and J. Jedrzejewski, https://doi.org/10.48550/arXiv.1004.2786
  * [39] N. Alves Jr. and C. S. O. Yokoi, Braz. Jour. Phys. 30, 667 (2000)
  * [40] Y. Wang, W. Guo, and H. W. J. Blöte, Phys. Rev. E 91, 032123 (2015)
  * [41] O. Rojas, Braz. J. Phys. 50, 675 (2020); Acta Phys. Pol. A 37, 933 (2019)
  * [42] R. A. Pimenta, O. Rojas, and S. M. de Souza, J. Magn. Magn. Mater. 550, 169070 (2022)