# Confinement-Deconfinement transition and $Z_{2}$ symmetry in $Z_{2}+$Higgs
theory

Minati Biswal<EMAIL_ADDRESS>Indian Institute of Science Education and
Research, Mohali 140306, India  Sanatan Digal<EMAIL_ADDRESS>The Institute
of Mathematical Sciences, Chennai 600113, India  Homi Bhabha National
Institute, Training School Complex, Anushakti Nagar, Mumbai 400085, India
Vinod Mamale<EMAIL_ADDRESS>The Institute of Mathematical Sciences,
Chennai 600113, India  Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar, Mumbai 400085, India  Sabiar Shaikh
<EMAIL_ADDRESS>The Institute of Mathematical Sciences, Chennai
600113, India  Homi Bhabha National Institute, Training School Complex,
Anushakti Nagar, Mumbai 400085, India

###### Abstract

We study the Polyakov loop and the $Z_{2}$ symmetry in the lattice
$Z_{2}+$Higgs theory in 4D Euclidean space using Monte Carlo simulations. The
results show that this symmetry is realised in the Higgs symmetric phase for
large number of temporal lattice sites. To understand the dependence on the
number of temporal sites, we consider a one dimensional model by keeping terms
of the original action corresponding to a single spatial site. In this
approximation the partition function can be calculated exactly as a function
of the Polyakov loop. The resulting free energy is found to have the $Z_{2}$
symmetry in the limit of large temporal sites. We argue that this is due to
$Z_{2}$ invariance as well as dominance of the distribution or density of
states corresponding to the action.

††preprint: IMSc/2021/02/01

##  I Introduction

$Z_{N}$ symmetry plays an important role in the confinement-deconfinement (CD)
transition in pure $SU(N)$ gauge theories tHooft:1977nqb ; McLerran:1981pb ;
Belyaev:1991gh . In these theories, at finite temperature, the allowed gauge
transformations are classified by the centre of the gauge group, i.e $Z_{N}$.
Under these $Z_{N}$ gauge transformations, i.e $Z_{N}$ symmetry, the Polyakov
loop ($L$) transforms like magnetisation, in spin models Kogut:1979wt . In the
confinement and deconfinement phases the Polyakov loop acquires vanishing and
non-zero thermal average values respectively, hence plays the role of an order
parameter for the confinement-deconfinement (CD) transition, Kuti:1980gh ;
McLerran:1980pk ; Weiss:1980rj ; Creutz:1980zw ; Svetitsky:1982gs . In the
deconfined phase, the $Z_{N}$ symmetry is spontaneously broken which leads to
$N$-degenerate states Yaffe:1982qf ; Svetitsky:1985ye ; Celik:1983wz .

The $Z_{N}$ symmetry of pure $SU(N)$ gauge theory is spoiled when matter
fields are included. Gauge transformations which are not periodic in temporal
directions can not act on the matter fields. These may act only on the gauge
fields but in the process the action does not remain invariant. There are many
studies on the effect of matter fields on this symmetry. Perturbative loop
calculations of the Polyakov loop effective potential show that this symmetry
is explicitly broken by matter fields in the fundamental representation
Weiss:1981ev ; Belyaev:1991np ; Guo:2018scp . The mean-field approximations of
lattice partition functions in the strong coupling limit also show the
explicit breaking of the $Z_{N}$ symmetry Green:1983sd ; Karsch:2000zv . Non-
perturbative studies of $CD$ transition in $2-$colour QCD show a sharp
transition suggesting small explicit breaking of $Z_{2}$ symmetry Satz:1985js
.

Recent non-perturbative Monte Carlo simulations of $SU(2)+$Higgs theory show
that the strength of $Z_{2}$ explicit breaking depends on the Higgs condensate
Biswal:2015rul . These studies find that the $CD$ transition exhibits critical
behaviour in the Higgs symmetric phase for large number of temporal sites
($N_{\tau}$) Biswal:2015rul . The distributions of the Polyakov loop are found
to be $Z_{2}$ symmetric, albeit within statistical errors, suggesting the
realisation of $Z_{2}$ symmetry in the Higgs symmetric phase. In reference
Biswal:2016xyq , it was argued that the emergence of $Z_{2}$ symmetry is due
to enhancement of the configuration/ensemble space with $N_{\tau}$. This
enhancement makes it possible that the change in the Euclidean action due to
$Z_{2}$ “rotation” of gauge links can be compensated by appropriately changing
the Higgs field. This was numerically tested by updating the Higgs field using
Monte Carlo steps after $Z_{2}$ rotating the gauge fields.

The non-invariance of the action under $Z_{2}$ gauge transformation which are
not periodic in temporal directions does not necessarily imply the explicit
breaking of $Z_{2}$ symmetry. The presence of $Z_{2}$ symmetry or it’s
explicit breaking can only be inferred from the free energy of the Polyakov
loop. In the free energy or the partition function calculations, two factors
play an important role. The distribution of the action also known as the
density of states ($DoS$) and the Boltzmann factor. The latter clearly does
not respect the $Z_{2}$ symmetry. So the realisation of the $Z_{2}$ symmetry
must come from the $DoS$ and it’s dominance over the Boltzmann factor.
Computing the $DoS$ in $SU(N)+$Higgs theory is a difficult task as the
configuration space is infinite. In this situation, the $Z_{2}+$Higgs theory
in four dimensions provides a suitable alternative. Since the field variables
take values $\pm 1$, it is possible to calculate the $DoS$ with some
simplifications.

The $Z_{2}+$Higgs theory has been extensively studied in literature
Balian:1974ts ; Balian:1974xw ; Creutz:1979he ; Fradkin:1978dv ;
Callaway:1981rt ; Bhanot:1981ug ; Bhanot:1980ga . The phase diagram of this
theory is found to be similar to that of $SU(N)+$Higgs theories in $3$ and
$4-$dimensions Jongeward:1980wx ; Creutz:1983ev . This theory has no continuum
limit and the temporal direction of $SU(N)$ gauge theories has no analog in
this case. For the purpose of understanding the role of $DoS$ on the $Z_{2}$
symmetry it is sufficient to consider the Higgs field to be periodic in any
one of the four dimensions. This will affect the gauge transformations, which
are not periodic in this direction but preserve the invariance of the pure
gauge action, exactly the same way matter fields affect the corresponding
gauge transformations in $SU(N)$ gauge theories. The Polyakov loop along this
“temporal direction” behaves as an order parameter for the $CD$ transition.

In this paper, we consider one of the four dimensions to define the “temporal
direction” and the Polyakov loop. The $Z_{2}$ symmetry of this observable and
the nature of $CD$ transition is then studied by varying the number of lattice
points, $N_{\tau}$, along the temporal direction. The computations are mostly
done on the Higgs symmetric side of the Higgs transition line. Our results
show that the $Z_{2}$ symmetry is realised for large $N_{\tau}$. Also the
behaviour of the $CD$ transition is found to be similar to the pure gauge case
apart from the location of the critical point. To understand the role of
$N_{\tau}$ a $0+1$ dimensional model is considered by keeping temporal
component of the gauge Higgs interaction corresponding to a given spatial
coordinate. The reason for this choice is the fact that only the temporal
component of the gauge Higgs interaction is sensitive to the $Z_{2}$ gauge
transformations. For the one dimensional model the Polyakov loop can take
values $\pm 1$. For each of these cases the free energy can be calculated
exactly. The free energy calculations show the emergence of $Z_{2}$ symmetry
in the large $N_{\tau}$ limit for arbitrary interaction coupling. Further the
Monte Carlo results for the distribution of the interaction term is reproduced
well by $0+1$ dimensional $DoS$ with a simple Boltzmann factor, though with a
different value of the interaction strength coupling. The $DoS$ for both
values of the Polyakov loop is sharply peaked at zero. $Z_{2}$ symmetry is
clearly observed near the peak, the differences appear when the action takes
the limiting values. Since the peak height grows with $N_{\tau}$, the $DoS$
will dominate the thermodynamics in the $N_{\tau}\to\infty$, restoring the
$Z_{2}$ symmetry.

This paper is organised as follows. In section II, we discuss the $Z_{2}$
symmetry in $Z_{2}+$Higgs theory. This is followed by numerical simulations of
$CD$ transition and the $Z_{2}$ symmetry in pure gauge theory and in the
presence of Higgs in section III. In section IV, we derive the free energy of
the Polyakov loop in a $0+1$ dimensional model, and relate the results to
$4-$dimensional Monte Carlo simulations. In section V, discussions and
conclusions are presented.

##  II  $Z_{2}$ symmetry in $Z_{2}+$Higgs gauge theory.

The action for the $Z_{2}+$Higgs theory in four dimensional Euclidean lattice
$(N_{s}^{3}\times N_{\tau})$ is given by,

$S=-\beta_{g}\sum_{P}U_{P}-\kappa\sum_{{n},\hat{\mu}}\Phi_{{n}+\hat{\mu}}U_{n,\hat{\mu}}\Phi_{n}.$
(1)

Here $n=(n_{1},n_{2},n_{3},n_{4})$ represents a point on the lattice with
$1\leq n_{1},n_{2},n_{3}\leq N_{s}$ and $1\leq n_{4}\leq N_{\tau}$. As
mentioned above we assume that the fourth direction is the temporal direction.
$U_{n,\hat{\mu}}$ represents the gauge links in $\hat{\mu}$ direction between
the lattice point $n$ and $n+\hat{\mu}$. The Higgs field $\Phi_{n}$ lives at
the site $n$. Both $U_{n,\hat{\mu}}$ and $\Phi_{n}$ take values $\pm 1$.
$\beta_{g}$ is the gauge coupling and $\kappa$ is the gauge Higgs interaction
strength. Figure. 1 shows a schematic layout of the gauge links and Higgs
variable on the lattice.

Figure 1: Position of gauge links U and Higgs fields $\Phi$ on lattice

The plaquette $U_{P}$ which is path ordered product of the links along an
elementary square on the $\mu-\nu$ plane, i.e

$U_{P}=U_{n,\hat{\mu}}U_{n+\hat{\mu},\hat{\nu}}U_{n+\hat{\nu},\hat{\mu}}U_{n,\hat{\nu}}.$
(2)

Figure. 2 shows the sketch of an elementary plaquette.

nn+$\hat{\mu}$n+$\hat{\mu}$+$\hat{\nu}$n+$\hat{\nu}$$U_{n,\hat{\mu}}$$U_{n+\hat{\mu},\hat{\nu}}$$U_{n+\hat{\nu},\hat{\mu}}$$U_{n,\hat{\nu}}$
Figure 2: Sketch of an elementary plaquette $U_{P}$

The pure gauge part of the action, first term in Eq.1, is invariant under the
$Z_{2}$ gauge transformations,

$U_{n,\hat{\mu}}\rightarrow V_{n}U_{n,\hat{\mu}}V^{-1}_{n+\hat{\mu}}$ (3)

where $V_{n}=\pm 1\in Z_{2}$. The $V_{n}$’s satisfy the following boundary
condition,

$V(\vec{n},n_{4}=1)=zV(\vec{n},n_{4}=N_{\tau}).$ (4)

$z=\pm 1\in Z_{2}$. So the gauge transformations can be classified by the
group $Z_{2}$. For $z=-1$ the gauge transformations are anti-periodic in the
temporal direction.

The Polyakov loop, which is defined as the product of links along the temporal
direction, i.e,

$L(\vec{n})=\displaystyle\prod_{n_{4}=1}^{N_{\tau}}U_{{(\vec{n},n_{4})},\hat{4}}$
(5)

transforms non-trivially under $Z_{2}$ gauge transformations Svetitsky:1982gs
. It is easy to see that the Polyakov loop transforms as,

$L(\vec{n})\rightarrow zL(\vec{n}).$ (6)

This transformation property of the Polyakov loop under $Z_{2}$ (or $Z_{N}$ in
general) gauge transformation is similar to that of magnetisation in the Ising
model. The partition function in the pure gauge case ($\kappa=0$) is given by,

${\cal Z}=\int DUe^{-S}.$ (7)

Since the action for $\kappa=0$ is invariant under $Z_{2}$ gauge
transformations, any configuration and it’s gauge rotated counterpart will
contribute equally to the partition function. Therefore the distribution of
the Polyakov loop exhibits $Z_{2}$ symmetry in this case. Equivalently the
free energy of the Polyakov loop will have $Z_{2}$ symmetry.

The presence of the Higgs field changes the space of allowed gauge
transformations. The reason being that the Higgs field is required to be
periodic in the temporal direction. Under a gauge transformation, $\Phi_{n}$
transforms as,

$\Phi_{n}\rightarrow V_{n}\Phi_{n}.$ (8)

Now the periodic boundary condition of $\Phi$ would be spoiled if non-periodic
gauge transformations, characterised by $z=-1$ are allowed. In this case given
a configuration, one can define a $Z_{2}$ counterpart in which only the gauge
links are $Z_{2}$ rotated. Obviously these pair of configurations will not
contribute equally to the partition function for $\kappa\neq 0$. So according
to the Boltzmann factor, $\sum_{\vec{n}}L(\vec{n})$ and
$-\sum_{\vec{n}}L(\vec{n})$ are non degenerate. This situation is similar to
the presence of an external field in the Ising model. However, the status of
$Z_{2}$ symmetry in the free energy can be answered only after integrating out
the Higgs field for a given $L(\vec{n})$ and it’s $Z_{2}$ rotated
configurations.

The Polyakov loop and Ising spins are similar in how they transform under
respective transformations. However there is an important difference between
them. This becomes clear when one compares $L(\vec{n})$ and an Ising spin at a
spatial point $\vec{n}=\\{n_{1},n_{2},n_{3}\\}$. A given value of $L(\vec{n})$
is associated with an entropy factor. This is because there are many different
combinations of $U_{(\vec{n},n_{4}),\hat{4}}$ and $\Phi_{\vec{n},n_{4}}$ are
possible for a given value of $L(\vec{n})$. Larger the $N_{\tau}$, larger is
the corresponding entropy. This aspect of the Polyakov loop needs to be taken
into account to understand the explicit breaking or realisation of $Z_{2}$
($Z_{N}$) symmetry, which is done in section IV. In the following section, we
discuss the algorithm of the Monte Carlo simulations Creutz:1979kf , present
simulation results for the phase diagram in the $\beta_{g}-\kappa$ plane,
distribution of the Polyakov loop and $CD$ transition in the Higgs symmetric
phase etc.

##  III Numerical technique and Monte Carlo simulation results.

In the Monte Carlo simulations, the Metropolis algorithm is used for sampling
the statistically significant configurations Hastings:1970aa . To update a
particular gauge link $U_{n,\mu}$, we consider the change in the action by
flipping it. If the action decreases then the flipped gauge link is accepted
for the new configuration. If the action increases by $\Delta S$ then the new
link is accepted with probability $Exp(-\Delta S)$. The same procedure is
adopted for $\Phi_{n}$. The process of updating is carried out over all $n$
and $\mu$ in multiple sweeps. Configurations separated by $10$ sweeps are used
in our analysis, which brings down the autocorrelation between successive
configurations to an acceptable level. For this simulations, $N_{\tau}=4-24$
and $N_{s}=16-84$ with $N_{s}/N_{\tau}=4$ lattices have been considered
Engels:1981ab .

The pure gauge simulations are initially performed to understand the nature of
$CD$ transition and $Z_{2}$ symmetry of the Polyakov loop. The simulations
were repeated in the presence of $\Phi$ to study its effects. The pure gauge
transition has been studied previously in the mean-field approximations
Balian:1974ts , which finds the transition is first order in four dimensions.
Also using duality transformations it can be shown that the critical
$\beta_{g}\sim 0.4407$ for $\kappa=0$ Balian:1974xw . These results are
supported by Monte Carlo simulations of smaller lattices Creutz:1979he . The
simulations carried out in this work are also consistent with these results.
In figure. 3 the average of the Polyakov loop is plotted $vs$ $\beta_{g}$ for
$N_{\tau}=4,8$. There is a range in $\beta_{g}$ for which clearly separated
peaks in the distribution of the Polyakov loop has been observed. We take
average of the Polyakov loop values corresponding to each peak separately.
Therefore we have two points in the figure for a given $\beta_{g}$. The two
peaks also suggest that the transition is first order. For larger lattice
sizes the range of $\beta_{g}$ over which two states are observed increases
Damgaard:1986jg . This is expected as strength of fluctuations relatively
decrease with volume (when correlation length is smaller than the spatial size
of the system), making it difficult for the field to climb over the barrier
and cross to the other side.

Figure 3: The average of the Polyakov loop vs $\beta_{g}$ for $N_{\tau}=$ 4
and 8.

The effect of the $\Phi$ field on the $CD$ transition and $Z_{2}$ symmetry is
expected to depend on $\kappa$. To relate these two aspects of pure gauge
theory to the phases of the Higgs field, simulations were performed to obtain
the Higgs transition line. For a given $\beta_{g}$, $\kappa>\kappa_{c}$
corresponds to the Higgs broken phase. In this phase the action term
dominates. For $\kappa<\kappa_{c}$ the fluctuations of the Higgs rather than
the action dominate the thermodynamic properties. This situation is similar to
the Ising model at high temperatures. In Fig.4 the Higgs transition line is
plotted in the $\beta_{g}-\kappa$ plane. The location of the phase boundary is
obtained by studying the $\kappa$ dependence of the interaction term and it’s
fluctuations for different values of $\beta_{g}$. In our simulations the Higgs
transition is found to be first order for intermediate range of $\beta$ and
crossover for both small and large $\beta$, as observed in previous studies
Jongeward:1980wx ; Creutz:1983ev . For large $\beta_{g}$ critical $\kappa_{c}$
remains flat and increases with $\beta_{g}$ in the small $\beta_{g}$ range. As
expected the critical values ($\beta_{c},\kappa_{c}$) depend on $N_{\tau}$,
however for $N_{\tau}$ values used in our simulations the dependence is very
mild.

Figure 4: Phase diagram

In the Higgs broken phase, i.e large $\kappa$, the interaction term dominates
over the entropy. The action takes the largest value when all the temporal
links are $+1$. So it is expected that in the Higgs phase $Z_{2}$ symmetry is
badly broken, also observed in our simulations. In the Higgs symmetric phase,
it is the fluctuations of Higgs in other words the distribution of the
interaction term dominate. In this phase there is a possibility for
realisation of $Z_{2}$ symmetry. In figure. 5 we show $CD$ transition in the
Higgs symmetric phase $(\kappa=.13)$. For comparison, $\kappa=0$ results also
have been included. The $CD$ transition is first order even in the presence of
$\Phi$, though the transition point shifts to lower values of $\beta_{g}$.

Figure 5: The average of the Polyakov loop vs $\beta_{g}$ for $N_{\tau}=$ 4
and 8.

To check the $N_{\tau}$ dependence of the $Z_{2}$ symmetry at $\kappa=.13$,
the distribution of Polyakov loop is computed both in the confined and the
deconfined phases for $N_{\tau}=2,3,$ and 8. In the deconfined phase, $L<0$
data is $Z_{2}$ rotated and then compared with $L>0$ data. The
distributions/histograms are shown in figures. 7-11. For $N_{\tau}=2$ the
histograms clearly show there is no $Z_{2}$ symmetry. In the deconfinement
side there is no $Z_{2}$ symmetry as the two Polyakov loop sectors do not
overlap. For $N_{\tau}=3$ the two peaks corresponding to the two sectors are
approaching towards each other. For $N_{\tau}=8$, the histogram of Polyakov
loop for two $Z_{2}$ sectors agree well with each other.

Figure 6: Histogram of L in the confined phase.

Figure 7: Histogram of L in the deconfined phase.

Figure 8: Histogram of L in the confined phase.

Figure 9: Histogram of L in the deconfined phase.

Figure 10: Histogram of L in the confined phase.

Figure 11: Histogram of L in the deconfined phase.

The $\kappa$ dependence of the $Z_{2}$ symmetry is studied by computing the
thermal average of the temporal part of the interaction, i.e
$sk_{4}=\sum_{n}\Phi_{n}U_{n,\hat{4}}\Phi_{n+\hat{4}}^{\dagger}$ and the
corresponding susceptibility $\chi_{sk_{4}}$. These simulations are carried
out in the deconfined phase, as there are two $Z_{2}$ states corresponding to
each sector of the Polyakov loop. The results for
($\left<sk_{4}\right>,\chi_{sk_{4}}$) are shown in figures. 15-15. For all
$N_{\tau}$ values the difference in ($\left<sk_{4}\right>,\chi_{sk_{4}}$) for
these two sectors is vanishingly small for small enough $\kappa$. For larger
$N_{\tau}$, the kappa value at which the two polyakov loop sectors differ
significantly in $sk_{4}$ and $\chi_{sk_{4}}$ is higher. For the largest
considered, $N_{\tau}=24$ , the two sectors agree in
($\left<sk_{4}\right>,\chi_{sk_{4}}$) up to the Higgs crossover point
$\kappa<\kappa_{c}$. When Higgs transition is first order the $Z_{2}$ symmetry
is observed in the Higgs symmetric phase even for $\kappa>\kappa_{c}$. Note
that for $\kappa>\kappa_{c}$ the Higgs symmetric phase is meta-stable.

Figure 12: $sk_{4}$ average vs $\kappa$ for $\beta_{g}=0.435$ on $16^{3}\times
4$ lattice

Figure 13: $sk_{4}$ fluctuation vs $\kappa$ for $\beta_{g}=0.435$ on
$16^{3}\times 4$ lattice

Figure 14: $sk_{4}$ average vs $\kappa$ for $\beta_{g}=0.435$ on $64^{3}\times
16$ lattice

Figure 15: $sk_{4}$ fluctuation vs $\kappa$ for $\beta_{g}=0.435$ on
$64^{3}\times 16$ lattice

It is clear from our $3+1$ dimensional simulations that the $Z_{2}$ symmetry
is realised in the Higgs symmetric phase for large $N_{\tau}$. To understand
these results, a $0+1$ dimensional model with temporal component of the
interaction term corresponding to a single spatial coordinates has been
considered..

##  IV The partition function and density of states in $0+1$ dimensions

The temporal component of the gauge Higgs interaction corresponding to a
particular spatial site can be written as,

$S_{1D}=-\kappa
sk_{4},~{}~{}~{}sk_{4}=\sum_{n=1}^{N_{\tau}}\Phi_{n}U_{n}\Phi_{n+1}.$ (9)

$n$ denotes the temporal lattice site, i.e $1\leq n\leq N_{\tau}$. $\Phi_{n}$
satisfies the periodic boundary condition $\Phi_{N_{\tau}+1}=\Phi_{1}$. Since
the action will not be invariant if a $z=-1$ gauge transformation is made on
$U_{i}$’s, the action breaks the $Z_{2}$ symmetry explicitly. For this model
the Polyakov loop can take values $\pm 1$. To see the $N_{\tau}$ dependence of
the $Z_{2}$ symmetry we calculate the free energy $V(L,N_{\tau})$. To simplify
the calculations we set $U_{i}=1$, for $i=1,2,...N_{\tau}-1$ and
$U_{N_{\tau}}=L$. All other configurations of $U_{i}$ corresponding to a given
value of $L$ are gauge equivalent. Now the partition function for $L=1$ is
nothing but that of the one dimensional Ising chain. For $L=-1$ the only
difference is that the coupling between $\Phi_{N_{\tau}}$ and $\Phi_{1}$ is
anti-ferromagnetic. For each choice of $L$ the partition function can be
calculated exactly, i.e,

${\cal Z}(L=1)=\lambda_{1}^{N_{\tau}}+\lambda_{2}^{N_{\tau}},~{}~{}{\cal
Z}(L=-1)=\lambda_{1}^{N_{\tau}}-\lambda_{2}^{N_{\tau}},$ (10)

where $\lambda_{1}=e^{\kappa}+e^{-\kappa}$ and
$\lambda_{2}=e^{\kappa}-e^{-\kappa}$. The corresponding free energies in large
$N_{\tau}$ are given by,

$V(L=1)=V(L=-1)=N_{\tau}{\rm log}(\lambda_{1}).$ (11)

This results show that there is $Z_{2}$ symmetry in $0+1$ dimensions in the
limit of $N_{\tau}\to\infty$.

Figure 16: $\rho(sk_{4})$ for $\kappa=0$ in 0+1

Figure 17: $\rho(sk_{4})$ for $\kappa=0$ in 0+1

Figure 18: $\rho(sk_{4})$ for $\kappa=0$ in 0+1

Figure 19: $\rho(sk_{4})$ for $\kappa=0$ in 0+1

As noted previously the restoration of the $Z_{2}$ symmetry must come from the
$Z_{2}$ symmetry of the entropy or the $DoS$. For $L=1$ the sequence of
allowed value of $sk_{4}$ is $\\{N_{\tau},N_{\tau}-4,......\geq-N_{\tau}\\}$.
On the other hand for $L=-1$ the corresponding sequence is
$\\{N_{\tau}-2,N_{\tau}-6,........\geq 2-N_{\tau}\\}$. The $DoS$ or
$\rho(sk_{4})$ for $N_{\tau}=4,8,12$ and $16$ are shown in figures.19-19. For
small $N_{\tau}$ there are clear difference for $L=\pm 1$. The difference
persists for the largest as well as smallest values of $sk_{4}$. For large
$N_{\tau}$, $\rho(sk_{4})$’s for both $L=\pm 1$ are well described by a
gaussian centred at $sk_{4}=0$, with $\sqrt{N_{\tau}}$ as standard deviation.
The logarithm of the peak hight is given by $~{}\simeq{\rm log}N_{\tau}!-2{\rm
log}(N_{\tau}/2)!+{\rm log}2$ for $N_{\tau}$ even. For $N_{\tau}=2n+1$ the
same can be approximated by ${\rm log}N_{\tau}!-{\rm log}(n^{2}+n)+{\rm
log2}$. The thermodynamics in the $N_{\tau}\to\infty$ limit will be dominated
by peak height and distribution of $\rho(sk_{4})$ around the peak, which is
$Z_{2}$ symmetric, for all finite $\kappa$. Interestingly this situation is
similar to one dimensional Ising chain where entropy dominates for any non-
zero finite temperature.

In order to take into account the effect of nearest neighbour coupling along
the spatial direction we consider $1+1$ dimensional model with $N_{s}=2$ and
vary $N_{\tau}$. In this case the Polyakov loop can take value $L=0,\pm 2$.
The exact calculation of $\rho(sk)$ get increasingly difficult with
$N_{\tau}$. One can however consider generating configurations randomly by
giving equal probability for each allowed value of a given variable. The
results for the distribution of the total action for $N_{\tau}=4$ and
$N_{\tau}=16$ are shown in Figs.21-21. As one can see that for higher
$N_{\tau}$, $\rho(sk)$ around the peak $sk=0$ do not depend on $L$.

Figure 20: $\rho(sk)$ for $\kappa=0$ in 0+1

Figure 21: $\rho(sk)$ for $\kappa=0$ in 0+1

To find out how well the $\rho(sk_{4})$ describe the Monte Carlo simulations
of the $4D$ partition function, the thermal average of the distribution
function $H(sk_{4})$ of $sk_{4}$ has been computed. For each configuration
$H(sk_{4})$ is given by the number of spatial sites with a given value of
$sk_{4}$. Note that the distribution of $sk_{4}$ takes into account the
Boltzmann factor which shifts the peak of $\rho(sk_{4})$ to the right.

Figure 22: $H(sk_{4})$ for $\kappa=0.1,~{}\beta_{g}=0.435$ for $3+1$ dimension

Figure 23: $H(sk_{4})$ fitted with $0+1$ density of states with a Boltzmann
factor

The figure. 23 shows the distribution $H(sk_{4})$ for $N_{\tau}=16$ at
$\kappa=0.1$ and $\beta_{g}=0.435$. For these values of $\kappa$ and
$\beta_{g}$, the system is found to be in the deconfined and Higgs symmetric
phase. The thermal average of the Polaykov loop for the two sectors are found
to be $\left<L\right>=0.5896\pm 0.002$ and $-0.5897\pm 0.00199$. Since the
$\left<L\right>\neq 1$ there is a smaller but finite fraction of spatial site
where the Polyakov loop takes opposite value. This results in the lower
envelope in $H(sk_{4})$. The results clearly show that $H(sk_{4})$ for both
the Polyakov loop sectors can be approximately described by single function in
other words the presence of $Z_{2}$ symmetry.

In figure. 23, we try to fit the $3+1$ dimensional simulation result with
$0+1$ dimensional $DoS$ by including an extra Boltzmann factor, i.e
exp$(\kappa^{\prime}sk_{4})$. The resulting fit agree very well with
$H(sk_{4})$. We expect that the $0+1$ results can describe the $3+1$ Monte
Carlo simulations in most of the phase diagram except for critical points.
Note here, $H(sk_{4})$ values correspond to $\kappa=0.1$, however to fit $DoS$
one needs a $\kappa$ value which is higher. This is due to the fact that in
$3+1$ dimensions $sk_{4}$ at a given spatial point interacts with $sk_{4}$ at
the nearest neighbour sites. Considering a mean-field approximation one can
compute the free energy difference between $L=1$ and $L=-1$ at
$\kappa=\kappa^{\prime}$ for the $3+1$ dimensional system at $\kappa=0.1$,
which turns out to be $10^{-10}$.

##  V Conclusions

In this paper the CD transition and $Z_{2}$ symmetry are studied in
$Z_{2}$+Higgs theory for $3+1$ dimensional Eulidean space. The results show
that for large $N_{\tau}$ the $Z_{2}$ symmetry is realized in the Higgs
symmetric phase within statistical errors. To understand the mechanism of
emergence of the $Z_{2}$ symmetry a simplified one dimension model of
$Z_{2}+$Higgs is considered by keeping only the temporal interaction terms at
a given spatial site. The partition function and the corresponding free energy
for each of the two Polyakov loop sectors is exactly calculated. It is shown
that the free energy difference between the two Polyakov loop sectors vanishes
in the large $N_{\tau}$ limit, which leads to $Z_{2}$ symmetry purely due to
dominance of entropy. The $DoS$ for finite $N_{\tau}$ are calculated exactly
where the asymmetry between the different Polyakov loop sectors rapidly
decreases with $N_{\tau}$. The effect of nearest neighbour interaction along
the spatial directions in a simple model shows the persistence of $Z_{2}$
symmetry in the $DoS$. Further it is shown that the $3+1$ Monte Carlo
simulations can be reproduced using the $DoS$ of the one dimensional model.

REFERENCES

## References

  * (1) G. ’t Hooft, Nucl. Phys. B 138, 1-25 (1978) doi:10.1016/0550-3213(78)90153-0
  * (2) L. D. McLerran and B. Svetitsky, Phys. Rev. D 24, 450 (1981) doi:10.1103/PhysRevD.24.450
  * (3) V. M. Belyaev, Phys. Lett. B 254, 153-157 (1991) doi:10.1016/0370-2693(91)90412-J
  * (4) J. B. Kogut, Rev. Mod. Phys. 51, 659 (1979) doi:10.1103/RevModPhys.51.659
  * (5) J. Kuti, J. Polonyi and K. Szlachanyi, doi:10.1016/0370-2693(81)90987-4
  * (6) L. D. McLerran and B. Svetitsky, doi:10.1016/0370-2693(81)90986-2
  * (7) N. Weiss, Phys. Rev. D 24, 475 (1981) doi:10.1103/PhysRevD.24.475
  * (8) M. Creutz, Phys. Rev. D 21, 2308-2315 (1980) doi:10.1103/PhysRevD.21.2308
  * (9) B. Svetitsky and L. G. Yaffe, Nucl. Phys. B 210, 423-447 (1982) doi:10.1016/0550-3213(82)90172-9
  * (10) L. G. Yaffe and B. Svetitsky, Phys. Rev. D 26, 963 (1982) doi:10.1103/PhysRevD.26.963
  * (11) B. Svetitsky, Phys. Rept. 132, 1-53 (1986) doi:10.1016/0370-1573(86)90014-1
  * (12) T. Celik, J. Engels and H. Satz, Phys. Lett. B 125, 411-414 (1983) doi:10.1016/0370-2693(83)91314-X
  * (13) N. Weiss, Phys. Rev. D 25, 2667 (1982) doi:10.1103/PhysRevD.25.2667
  * (14) V. M. Belyaev, I. I. Kogan, G. W. Semenoff and N. Weiss, Phys. Lett. B 277, 331-336 (1992) doi:10.1016/0370-2693(92)90754-R
  * (15) Y. Guo and Q. Du, JHEP 05, 042 (2019) doi:10.1007/JHEP05(2019)042 [arXiv:1810.13090 [hep-ph]].
  * (16) F. Green and F. Karsch, Nucl. Phys. B 238, 297-306 (1984) doi:10.1016/0550-3213(84)90452-8
  * (17) F. Karsch, E. Laermann, A. Peikert, C. Schmidt and S. Stickan, Nucl. Phys. B Proc. Suppl. 94, 411-414 (2001) doi:10.1016/S0920-5632(01)00988-4 [arXiv:hep-lat/0010040 [hep-lat]].
  * (18) H. Satz, Phys. Lett. B 157, 65-69 (1985) doi:10.1016/0370-2693(85)91213-4
  * (19) M. Biswal, S. Digal and P. S. Saumia, Nucl. Phys. B 910, 30-39 (2016) doi:10.1016/j.nuclphysb.2016.06.025 [arXiv:1511.08295 [hep-lat]].
  * (20) M. Biswal, M. Deka, S. Digal and P. S. Saumia, Phys. Rev. D 96, no.1, 014503 (2017) doi:10.1103/PhysRevD.96.014503 [arXiv:1610.08265 [hep-lat]].
  * (21) R. Balian, J. M. Drouffe and C. Itzykson, Phys. Rev. D 10, 3376 (1974) doi:10.1103/PhysRevD.10.3376
  * (22) R. Balian, J. M. Drouffe and C. Itzykson, Phys. Rev. D 11, 2104 (1975) [erratum: Phys. Rev. D 19, 2514 (1979)] doi:10.1103/PhysRevD.11.2104
  * (23) M. Creutz, doi:10.1103/PhysRevD.21.1006
  * (24) E. H. Fradkin and S. H. Shenker, Phys. Rev. D 19, 3682-3697 (1979) doi:10.1103/PhysRevD.19.3682
  * (25) D. J. E. Callaway and L. J. Carson, Phys. Rev. D 25, 531-537 (1982) doi:10.1103/PhysRevD.25.531
  * (26) G. Bhanot and B. A. Freedman, Nucl. Phys. B 190, 357-364 (1981) doi:10.1016/0550-3213(81)90566-6
  * (27) G. Bhanot and M. Creutz, BNL-27833.
  * (28) G. A. Jongeward and J. D. Stack, Phys. Rev. D 21, 3360 (1980) doi:10.1103/PhysRevD.21.3360
  * (29) M. Creutz, L. Jacobs and C. Rebbi, Phys. Rept. 95, 201-282 (1983) doi:10.1016/0370-1573(83)90016-9
  * (30) M. Creutz, L. Jacobs and C. Rebbi, Phys. Rev. Lett. 42, 1390 (1979) doi:10.1103/PhysRevLett.42.1390
  * (31) W. K. Hastings, Biometrika 57, 97-109 (1970) doi:10.1093/biomet/57.1.97
  * (32) J. Engels, F. Karsch and H. Satz, Nucl. Phys. B 205, 239-252 (1982) doi:10.1016/0550-3213(82)90387-X
  * (33) P. H. Damgaard and U. M. Heller, Phys. Lett. B 171, 442-448 (1986) doi:10.1016/0370-2693(86)91436-X


11institutetext: IISlab, Technical University of Kosice, Kosice, Slovakia
22institutetext: 2nd Department of Internal Medicine, Pavol Jozef Safarik
University and Louis Pasteur University Hospital, Kosice, Slovakia

# Ensembled Autoencoder Regularization for Multi-Structure Segmentation for
Kidney Cancer Treatment

David Jozef Hresko 11 Marek Kurej Jakub Gazda 1122 Peter Drotar 11

###### Abstract

The kidney cancer is one of the most common cancer types. The treatment
frequently include surgical intervention. However, surgery is in this case
particularly challenging due to regional anatomical relations. Organ
delineation can significantly improve surgical planning and execution. In this
contribution, we propose ensemble of two fully convolutional networks for
segmentation of kidney, tumor, veins and arteries. While SegResNet
architecture achieved better performance on tumor, the nnU-Net provided more
precise segmentation for kidneys, arteries and veins. So in our proposed
approach we combine these two networks, and further boost the performance by
mixup augmentation.

###### Keywords:

segmentation kidney SegResNet nnU-Net ensemble mixup

##  1 Introduction

Kidney cancer (KC) has become one of the ten most common cancer types in the
general population, while its incidence steadily increased from the 1970s
until the mid-1990s. This trend was attributed to the improved diagnosis
secured by the broad introduction of advanced radiological imaging to clinical
practice [1]. However, KC is strongly associated with risk factors such as
cigarette smoking, obesity, and arterial hypertension. Therefore, even though
the incidence has recently levelled off, it is not expected to drop anytime
soon. On the contrary, it will pose a significant threat in industrialized
countries. The treatment of KC is based on ablation, chemotherapy, and surgery
[2], [11]. Luckily, most KCs are detected early and only incidentally by
radiologic imaging for other diseases when the curative treatment - surgery -
is still possible. However, surgery is specifically challenging due to complex
regional anatomical relations. Kidney parsing can help tackle this issue and
improve preoperative planning and perioperative decisions, leading to a higher
chance of successful tumour resection.

Abdominal organ segmentation from medical imaging is a topic that has been
researched for several years. Particularly, the computed tomography can
provide detailed information about organs and structures in abdominal area.
This information can be used for surgical planning, diagnosis and various
clinical decisions. In recent years the contour delineation of anatomical and
pathological structures became important for applications such as visual
support during surgery. Since the manual delineation is tedious and time
consuming the automated highly accurate methods are strongly desired.

The advent of convolutional neural networks brought many efficient solution
also for medical images segmentation. The fully convolutional neural network,
particularly U-Net architecture [9] became de-facto golden standard for
medical image segmentation. Moreover, the pipeline optimisation process
implemented in nnU-Net [5] boosted segmentation performance and made the U-net
easily applicable for any segmentation task. Considering all these recent
advances, the abdominal organ segmentation appears to be almost solved problem
[6]. However, even though this may be partially true for a big organs such as
liver, kidneys and spleen, the smaller structures still represent challenge
for automated segmentation approaches.

Challenges, like LiTS, KiTS, FLARE did not consider smaller structure namely
arteries and veins. However segmentation of arteries and veins is crucial for
urological surgeries namely laparoscopic nephrectomy [10],[8]. In this case
the precise information about target segmental arteries is required in order
to avoid insufficient clamping.

In this paper we propose ensemble of two encoder-decoder based convolutional
neural networks: nnU-net and SegResNet [7]. To further boost the performance
of nnU-net, we additionally applied manifold mixup augumentation [12]. Similar
improvement was already proposed by authors of [3] on KiTS21 challenge, which
improved overall performance of their solution for kidney and kidney tumor
segmentation task.

The rest of the paper is organized as follows. In the next section we provide
detailed data description. In the methodology section, the proposed solution
is explained. Finally we present the results and discuss different aspects of
our submission.

##  2 Data

The data were acquired by Siemens dual-source 64-slice CT scanner and the
contrast media was injected during CT image acquisition. The further technical
settings of CT are: X-ray tube current is 480 mA, B25F convolution kernel,
exposure time equal to 500 ms and voltage 120KV. The slice thickness is
0.5mm/pixel and spacing of images is from 0.47 mm/pixel to 0.74 mm/pixel and
0.75 mm/pixel for z-direction, respectively.

Altogether there are 130 3D abdominal CT images. The ground-truth
corresponding to four classes, kidney, tumor, vein and artery, was determined
by three medically trained experts and validated by experienced radiologist.
All images were cropped to the same size of $150\times 150\times 200$ to focus
on the four structures of interests

From 130 cases, 70 are used for training and validation of the model, 30 for
open test phase evaluation and 30 for closed test phase evaluation.

##  3 Methodology

Besides the state-of-art architecture nnU-Net, which is commonly used to solve
medical segmentation tasks, we decided to also adapt and verify not so well-
known architectures. To be more concrete, we chose autoencoder based
architecture with additional regularization named SegResNet. Based on the
obtained results from training phase we decided to ensemble the output from
nnU-Net and SegResNet, which were later additionally fine-tuned for closed
test phase of KiPA22 challenge.

###  3.1 Preprocessing

In case of the nnU-Net architecture, preprocessing methods including
transformation, re-sampling, normalization and scaling were automatically
handled by the nnU-Net pipeline. This includes, the resampling strategy for
anisotropic data was third order spline interpolation for in-plane and nearest
neighbours for out-of-plane. Global dataset percentile clipping along with
z-score with global foreground mean and standard deviation was chosen as a
normalization strategy. The clipping of HU values were automatically preset
with default setting (0.5 and 99.5 percentile).

In case of the SegResNet architecture, as the first step we applied regular
channel-wise normalization and scaling to clip CT scan values. After that
several augmentations were applied, each with the probability of 0.5.
Concretely, we randomly noised original scan with gaussian noise, followed by
random rotation, zooming, axis flipping and elastic deformations using
bilinear interpolation to calculate output values.

###  3.2 Post-processing

Based on our observations during training phase, we noticed that some faulty
segmentations of kidney and tumor structures were mainly caused by small
redundant segments, which were located far from the expected structure
location. To prevent this phenomena, we decided to apply keep largest
connected component method as the post-processing step during separate
inference of both models.

###  3.3 Proposed approach

Our proposed method relies on ensemble technique and consists of two different
encoder-decoder based variations of U-Net architecture. The first one is the
nnU-Net architecture, which already dominated several segmentation challenges.
Following the approach of [3] we also took advantege of mixup augmentation.
Instead of regular version of mixup [13], we utilized extended version named
manifold mixup [12]. The key difference here is that, the mixing process of
batch samples can be also applied on inner layers of the network, not only on
the input layer. The mixed minibatch is then passed forward regularly from
$k$-th inner layer to the output layer, which is used to compute the loss
value and gradients.

During extensive experiments performed on KiPA22 public dataset we discovered
that, there is an architecture, which is capable of better results on some of
the segmented classes. Concretely, performance of the nnU-Net was outperformed
by SegResNet on tumor class. To benefit from both architectures, we decided to
perform ensemble of these two architectures to better localize targeted
structures. Based on the obtained results from individual architectures, we
keep predictions of kidney, renal vein and renal artery class from nnU-Net,
while tumor were produced by SegResNet. Overall concept of proposed method is
depicted in Figure 1.

Figure 1: Overview of proposed ensemble method

##  4 Results

The evaluation of the model performance is based on two different approaches.
The first one relies on the area based metric, concretely Dice Similarity
Coefficient (DSC) is used to evaluate the area-based overlap index. The second
one is based on distance between predicted and ground truth segmentation mask.
Here, to evaluate the coincidence of the surface for less sensitive to
outliers, Average Hausdorff Distance (AHD) is used. Additionally, outliers
sensitive Hausdorff Distance (HD) is also used to further evaluate the
segmentation quality.

To train the final version of nnU-Net, we used stochastic gradient descent
(SGD) optimizer with the initial learning rate of 0.01. The length of training
was set to 1000 epochs. For manifold mixup augmentation we set the
hyperparameter value of mixing coefficient to $\alpha$ = 0.1. In case of the
final version of SegResNet, we used AdamW optimizer with the initial learning
rate of 0.0001. The length of training was set to 4000 epochs. For both models
we used combined Dice and Cross Entropy loss, batch size equal to two.

The results from the training phase, along with comparisons to other
experiments with different models and configurations are presented in Table 1.
As can be seen, proposed manifold mixup improved performance of the nnU-Net
for each class. SegResNet outperformed nnU-Net for tumor class segmentation
task with notable differences. Based on the [4] we were also curious if
transformer based architecture is truly capable of better results than
standard convolutional network, so we also trained UNETR model on this
dataset. Our results proved that this vision transformer was not able to reach
the performance of any baseline version of tested models.

Table 1: Overall performance of examined models on KiPA22 public train data | Kidney | Tumor | Vein | Artery |
---|---|---|---|---|---
Network | DSC | HD | AHD | DSC | HD | AHD | DSC | HD | AHD | DSC | HD | AHD |
| nnU-Net
---
(baseline)
0.960 | 17.108 | 0.484 | 0.856 | 12.432 | 1.412 | 0.805 | 12.002 | 1.180 | 0.839 | 15.995 | 0.512 |
| nnU-Net
---
+
mixup
0.963 | 17.022 | 0.424 | 0.893 | 10.114 | 1.225 | 0.823 | 12.112 | 0.830 | 0.849 | 16.395 | 0.449 |
UNETR | 0.951 | 19.434 | 0.874 | 0.838 | 14.768 | 1.992 | 0.773 | 18.109 | 3.556 | 0.822 | 20.542 | 1.987 |
| SegResNet
---
0.961 | 17.045 | 0.512 | 0.901 | 9.100 | 0.970 | 0.818 | 11.224 | 1.952 | 0.873 | 19.156 | 2.443 |
Figure 2: Example segmentation for the 61 case. Individual contours denote
ground truth label for specific class. Predicted segmentation classes are:
yellow=kidney, blue=tumor, green=renal vein, red=renal artery.

The Fig.2 shows the example prediction of our proposed method. We randomly
selected case 61 from public dataset, which contains all four labels. As can
be seen, the segmentation of all structures is almost identical with the
ground truth regions. The yellow colour represents healthy kidney tissue
located laterally to renal vessels. The red and green colours medially
represent the renal artery and renal vein, respectively. Finally, the blue
colour dorsally represents kidney cancer. This picture informs the surgeon
that the kidney cancer has not invaded renal vessels yet and has a high chance
of successful resection. However, a surgeon must evaluate all sections, which
may present different spatial relation between kidney cancer and renal
vessels. In more complicated cases, human visual evaluation might not be
enough to evaluate the ground truth correctly.

Overall performance of our proposed method was also independently measured on
KiPA22 challenge test data during open test phase and closed test phase of the
challenge. The obtained results from open test phase can be seen in the table
2. The results from the closed test phase were not available during the time
of writing this paper.

Table 2: Overall performance of proposed method on KiPA22 open test data | Kidney | Tumor | Vein | Artery |
---|---|---|---|---|---
Network | DSC | HD | AHD | DSC | HD | AHD | DSC | HD | AHD | DSC | HD | AHD |
| Ensembled
---
model
0.957 | 17.008 | 0.464 | 0.880 | 9.104 | 1.616 | 0.835 | 12.733 | 0.322 | 0.878 | 16.195 | 0.322 |

##  5 Discussion

Here we investigate in detail some of the segmentation cases that reached
lower scores. Fig. 3 represents ”a double incorrect case” because neither the
ground truth nor the prediction is correct. First, the ground truth (yellow
contour, representing healthy renal parenchyma) contains a pathological
(probably cystic) lesion (it is impossible to tell without considering other
phases of the contrast enhancement). Second, the prediction ignores this
lesion and goes around it without rendering and classifying the pathological
lesion. Furthermore, Fig. 4 represents ”an incorrect prediction case”. The
ground truth (blue contour) depicts a tumour located in the dorsal part of the
left kidney. However, the segmentation fails to render the very dorsal part of
the tumour, which again appears as its cystic component. There are no density
differences between the cystic component of the tumour and the healthy
perirenal fat here, which could be why the prediction failed. Altogether,
while it seems that the renal parenchyma, renal vessels and solid tumours have
been segmented correctly, the mistakes are associated mostly with cystic
lesions or cystic components of kidney cancers.

Figure 3: Example fault segmentation for the 6 case. Individual contours
denote ground truth label for specific class. Predicted segmentation classes
are: yellow=kidney, green=renal vein, red=renal artery. Figure 4: Example
fault segmentation for the 32 case. Individual contours denote ground truth
label for specific class. Predicted segmentation classes are: yellow=kidney,
blue=tumor, green=renal vein, red=renal artery.

##  6 Conclusions

In this paper we proposed ensemble method for segmentation of structures in
abdominal area. The ensemble is build on combination of nnU-Net and SegResNet
architecture. The method achieves highly competitive score on all four
structures.

## References

  * [1] Chow, W.H., Dong, L.M., Devesa, S.S.: Epidemiology and risk factors for kidney cancer. Nature Reviews Urology 7(5), 245–257 (2010)
  * [2] Dahle, D.O., Skauby, M., Langberg, C.W., Brabrand, K., Wessel, N., Midtvedt, K.: Renal cell carcinoma and kidney transplantation: a narrative review. Transplantation 106(1), e52–e63 (2022)
  * [3] Gazda, M., Bugata, P., Gazda, J., Hubacek, D., Hresko, D.J., Drotar, P.: Mixup augmentation for kidney and kidney tumor segmentation. In: Heller, N., Isensee, F., Trofimova, D., Tejpaul, R., Papanikolopoulos, N., Weight, C. (eds.) Kidney and Kidney Tumor Segmentation. pp. 90–97. Springer International Publishing, Cham (2022)
  * [4] Hatamizadeh, A., Tang, Y., Nath, V., Yang, D., Myronenko, A., Landman, B., Roth, H., Xu, D.: Unetr: Transformers for 3d medical image segmentation (2021). https://doi.org/10.48550/ARXIV.2103.10504, https://arxiv.org/abs/2103.10504
  * [5] Isensee, F., Jaeger, P.F., Kohl, S.A.A., Petersen, J., Maier-Hein, K.H.: nnu-net: a self-configuring method for deep learning-based biomedical image segmentation. Nature Methods 18(2), 203–211 (Feb 2021). https://doi.org/10.1038/s41592-020-01008-z, https://doi.org/10.1038/s41592-020-01008-z
  * [6] Ma, J., Zhang, Y., Gu, S., Zhu, C., Ge, C., Zhang, Y., An, X., Wang, C., Wang, Q., Liu, X., Cao, S., Zhang, Q., Liu, S., Wang, Y., Li, Y., He, J., Yang, X.: Abdomenct-1k: Is abdominal organ segmentation a solved problem? (2020). https://doi.org/10.48550/ARXIV.2010.14808, https://arxiv.org/abs/2010.14808
  * [7] Myronenko, A.: 3d mri brain tumor segmentation using autoencoder regularization (2018). https://doi.org/10.48550/ARXIV.1810.11654, https://arxiv.org/abs/1810.11654
  * [8] Porpiglia, F., Fiori, C., Checcucci, E., Amparore, D., Bertolo, R.: Hyperaccuracy three-dimensional reconstruction is able to maximize the efficacy of selective clamping during robot-assisted partial nephrectomy for complex renal masses. European Urology 74(5), 651–660 (2018). https://doi.org/https://doi.org/10.1016/j.eururo.2017.12.027, https://www.sciencedirect.com/science/article/pii/S0302283817310953
  * [9] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Navab, N., Hornegger, J., Wells, W.M., Frangi, A.F. (eds.) Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015. pp. 234–241. Springer International Publishing, Cham (2015)
  * [10] Shao, P., Qin, C., Yin, C., Meng, X., Ju, X., Li, J., Lv, Q., Zhang, W., Xu, Z.: Laparoscopic partial nephrectomy with segmental renal artery clamping: Technique and clinical outcomes. European Urology 59(5), 849–855 (2011). https://doi.org/https://doi.org/10.1016/j.eururo.2010.11.037, https://www.sciencedirect.com/science/article/pii/S0302283810011395
  * [11] Shao, P., Tang, L., Li, P., Xu, Y., Qin, C., Cao, Q., Ju, X., Meng, X., Lv, Q., Li, J., Zhang, W., Yin, C.: Precise segmental renal artery clamping under the guidance of dual-source computed tomography angiography during laparoscopic partial nephrectomy. European Urology 62(6), 1001–1008 (2012). https://doi.org/https://doi.org/10.1016/j.eururo.2012.05.056, https://www.sciencedirect.com/science/article/pii/S0302283812006409
  * [12] Verma, V., Lamb, A., Beckham, C., Najafi, A., Mitliagkas, I., Courville, A., Lopez-Paz, D., Bengio, Y.: Manifold mixup: Better representations by interpolating hidden states (2018). https://doi.org/10.48550/ARXIV.1806.05236, https://arxiv.org/abs/1806.05236
  * [13] Zhang, H., Cisse, M., Dauphin, Y.N., Lopez-Paz, D.: mixup: Beyond empirical risk minimization (2017). https://doi.org/10.48550/ARXIV.1710.09412, https://arxiv.org/abs/1710.09412


# The impact of Rashba spin-orbit coupling in charge-ordered systems

Rodrigo A. Fontenele  Instituto de Física, Universidade Federal do Rio de
Janeiro Cx.P. 68.528, 21941-972 Rio de Janeiro RJ, Brazil  Sebastião dos Anjos
Sousa Júnior  Instituto de Física, Universidade Federal do Rio de Janeiro
Cx.P. 68.528, 21941-972 Rio de Janeiro RJ, Brazil  Tarik P. Cysne  Instituto
de Física, Universidade Federal Fluminense, 24210-346 Niterói RJ, Brazil
Natanael C. Costa  Instituto de Física, Universidade Federal do Rio de Janeiro
Cx.P. 68.528, 21941-972 Rio de Janeiro RJ, Brazil

###### Abstract

We study the impact of the Rashba spin-orbit coupling (RSOC) on the stability
of charge-density wave (CDW) in systems with large electron-phonon coupling
(EPC). Here, the EPC is considered in the framework of the Holstein model at
the half-filled square lattice. We start obtaining the phase diagram of the
Rashba-Holstein model using the Hartree-Fock mean-field theory, and
identifying the boundaries of the CDW and Rashba metal phases. As our main
result, we notice that the RSOC disfavors the CDW phase, driving the system to
a correlated Rashba metal. Proceeding, we employ a cluster perturbation theory
(CPT) approach to investigate the phase diagram beyond the Hartree-Fock
approximation. The quantum correlations captured by CPT indicate that the RSOC
is even more detrimental to CDW than previously anticipated. That is, the
Rashba metal region is observed to be expanded in comparison to the mean-field
case. Additionally, we investigate pairing correlations, and the results
further strengthen the identification of critical points.

##  I Introduction

A growing class of two-dimensional (2D) materials offers opportunities to
explore condensed-matter phenomena in a new and flexible platform Mas-Ballesté
_et al._ (2011); Gupta  _et al._ (2015); Novoselov  _et al._ (2016); Ahn
(2020). In addition to their technological applications, the prospect of
solving some old puzzles in solid-state physics motivates much of the
fundamental studies in this field. Graphene was the first genuinely 2D
material synthesized, but apart from specific situations Kotov  _et al._
(2012), electronic correlations do not play a significant role on the observed
properties, due to their highly dispersive band structure Zhang  _et al._
(2019); Paiva  _et al._ (2005). Therefore, part of the research on $2$D
materials was developed within the single-particle perspective and neglecting
many-body effects Castro Neto  _et al._ (2009). With the progress in the
quality of samples and synthesis of novel families of materials, a renewed
interest in the phases of matter driven by electronic correlation reappeared
in the $2$D materials community. Among those compounds of interest, we mention
the superconducting phase in twisted bilayer graphene Cao  _et al._ (2018),
the magnetism in Cr2Ge2Te6 Gong  _et al._ (2017) and CrI3 Huang  _et al._
(2017), and the charge-density wave (CDW) phase in some transition metal
dichalcogenide (TMDs) family Manzeli  _et al._ (2017).

For a CDW phase, the key feature is the occurrence of an electron-phonon
coupling (EPC), which creates instabilities in the metallic system, leading to
a ground state with a spatially-modulated charge order Grüner (1988); Gorḱov
and Grüner (1989). Such a state is characterized by a gapped single-particle
charge spectrum and a gapless collective mode formed by phonon excitations.
Interestingly, CDW and superconducting phases have been observed in TMDs
compounds which stabilize in the centrosymmetric (1T) Sugawara  _et al._
(2015) and the non-centrosymmetric (2H) Lin  _et al._ (2020) structures. In
addition, due to the high atomic number of transition metal atoms, the
compounds in the TMD family typically have a strong _intrinsic_ spin-orbit
coupling (SOC). However, the interplay between the leading effects of a SOC
and an EPC, regarding in particular how the SOC affects CDW formation, is
still an open issue. Indeed, the study of the interplay between EPC and SOC in
TMDs is challenging due to its complex orbital nature and tri-layer structure,
with some works suggesting a SOC-induced renormalization of the electron-
phonon coupling Ge and Liu (2012); Lian  _et al._ (2022).

Different types of spin-orbit interaction may also appear in materials
depending on the point-group symmetries of their lattice structure. In
addition to the intrinsic SOC characteristic of materials with high atomic
numbers, other types of SOCs can be induced in materials by _extrinsic_
mechanisms. For instance, the Rashba SOC (RSOC) occurs in systems with broken
$(z\rightarrow-z)$-symmetry, frequently appearing in surface-confined
electrons of semiconductor heterostructures Nitta  _et al._ (1997), alloys Ast
_et al._ (2007); Rotenberg  _et al._ (1999), and surface-states of metallic
systems LaShell  _et al._ (1996). It can also be induced in 2D materials by
proximity to suitable substrates Min  _et al._ (2006); Cysne  _et al._ (2018a,
b). In this context, calculations using density functional theory indicate
that pristine lead chalcogenides monolayers [Pb$X$, with $X$= S, Se, Te] are
2D materials with strong RSOC Hanakata  _et al._ (2017). Its monolayers
possess a buckled square lattice that breaks $(z\rightarrow-z)$-symmetry,
allowing the occurrence of the Rashba effect, whose intensity is adjusted by
changing the buckling angle, via application of strain Hanakata  _et al._
(2018). The RSOC may also occur in other 2D materials such as 2D perovskites,
and Janus monolayer of transition metals dichalcogenides Chen  _et al._
(2021).

One significant impact of the RSOC is the emergence of a spin-dependent
momentum shift, resulting in spin-splitting effects in the band structure.
This spin splitting can have profound consequences on the electronic
properties of the system, particularly on electron correlations, directly
affecting the emergence of new phases. For instance, in certain families of
iridates and pyrochlores – materials that demonstrate strong electron-electron
and spin-orbit interactions – the interplay between RSOC and electron-electron
interactions may give rise to topologically nontrivial states, Weyl
semimetals, or even Kitaev spin liquids Pesin and Balents (2010); Witczak-
Krempa  _et al._ (2014); Rau  _et al._ (2016). Despite recent advancements in
this field, the effects of RSOC on electron-phonon interacting systems, such
as those compounds mentioned earlier, remain less clear.

In view of these stimulating results, in this paper we investigate the impact
of RSOC in a charge-ordered system. Here, we consider dispersionless phonon
degrees of freedom coupled to the electrons, within the framework of the
Holstein model Holstein (1959): a simplified model whose ground state at the
half-filled square lattice exhibits a staggered CDW phase for any electron-
phonon coupling Costa  _et al._ (2020). This charge-ordered phase may undergo
significant changes when external parameters are introduced – such as
electronic doping Dee  _et al._ (2019); Bradley  _et al._ (2021), pressure or
strain Cohen-Stead  _et al._ (2019); Araújo  _et al._ (2022), nonlinear
electron-phonon couplings Li  _et al._ (2015); Dee  _et al._ (2020); Paleari
_et al._ (2021), disorder Xiao  _et al._ (2021), doping with magnetic
impurities Sousa-Júnior  _et al._ (2023), or phonon dispersion Costa  _et al._
(2018) –, with the enhancement of pairing correlations for most of cases.
Despite this scrutinization, the impact of a SOC on the Holstein CDW phase is
barely unknown. We bridge this gap by examining the Rashba-Holstein
Hamiltonian through two methodologies: mean-field (MF) and cluster
perturbation theory (CPT) approaches.

The paper is organized as follows: we present the Rashba-Holstein model in
Sec. II, with the results for the MF and CPT approaches being discussed in
Secs. III and IV, respectively. In Section V we present our conclusions.

##  II The Rashba-Holstein Hamiltonian

We consider electrons hopping in a single orbital square lattice. Here, we
assume the system has a broken $(z\rightarrow-z)-$symmetry, as occurs in lead
chalcogenides monolayers Hanakata  _et al._ (2017, 2018) and in confined
electron gases Nitta  _et al._ (1997); Ast  _et al._ (2007); Rotenberg  _et
al._ (1999); LaShell  _et al._ (1996), which allows the occurrence of RSOC. In
addition, we use the Holstein model to describe the electron-phonon coupling,
which captures the local interaction of electrons with single-Einstein phonon
degrees of freedom Holstein (1959).

For the Rashba-Holstein model, the Hamiltonian reads

$\displaystyle\mathcal{H}=$
$\displaystyle-t\sum_{\langle\textbf{i},\textbf{j}\rangle,\sigma}(c_{\textbf{i},\sigma}^{\dagger}c_{\textbf{j},\sigma}+h.c.)-\mu\sum_{\textbf{i},\sigma}n_{\textbf{i},\sigma}$
(1)
$\displaystyle+\hbar\omega_{0}\sum_{\textbf{i}}a_{\textbf{i}}^{\dagger}a_{\textbf{i}}-g\sum_{\textbf{i},\sigma}n_{\textbf{i},\sigma}(a_{\textbf{i}}+a_{\textbf{i}}^{\dagger})$
$\displaystyle-
t_{R}\sum_{\langle\textbf{i},\textbf{j}\rangle}c_{\textbf{i}}^{\dagger}\left({\overleftrightarrow{\alpha}}_{\textbf{i},\textbf{j}}\times{\overleftrightarrow{\sigma}}\right)_{z}c_{\textbf{j}}~{},$

where the sums run over a square lattice, with
$\langle\textbf{i},\textbf{j}\rangle$ denoting nearest neighbors. Here we use
the usual second quantization formalism, with $c_{\textbf{i},\sigma}$
($c_{\textbf{i},\sigma}^{\dagger}$) being a fermionic operator that
annihilates (creates) an electron with spin $\sigma=\uparrow,\downarrow$ on a
given site $\mathbf{i}$;
$n_{\textbf{i},\sigma}=c_{\textbf{i},\sigma}^{\dagger}c_{\textbf{i},\sigma}$
is the fermion number operator. Similarly, $a_{\textbf{i}}$
($a_{\textbf{i}}^{\dagger}$) is a bosonic operator that annihilates (creates)
a phonon with energy $\hbar\omega_{0}$ on a given site $\mathbf{i}$. In the
last term of Eq. (1),
${\overleftrightarrow{\alpha}}_{\textbf{i},\textbf{j}}\equiv\left(\alpha_{\textbf{i},\textbf{j}}^{x},\alpha_{\textbf{i},\textbf{j}}^{y},0\right)$
is the displacement tensor $\alpha_{\textbf{i},\textbf{j}}^{\nu}\equiv
i(\delta_{\textbf{i},\textbf{j}+\alpha_{\nu}}-\delta_{\textbf{i},\textbf{j}-\alpha_{\nu}})$,
with $\alpha_{\nu}$ being unitary translation in the $\nu$ direction, and
${\overleftrightarrow{\sigma}}\equiv\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)$
are Pauli matrices related to electronic spin. The first two terms on the
right-hand side of Eq. (1) correspond to the hopping of the electrons and
their chemical potential, respectively. The third term describes the local
harmonic oscillators, while the fourth represents the EPC, with strength $g$.
Finally, the last term is the RSOC, with strength $t_{R}$. We set the energy
scale in units of the hopping integral $t$.

At this point, it is also worth mentioning that the ground state of the half-
filled pure Holstein model ($t_{R}=0$) is CDW for any EPC $g>0$, due to
instabilities from the van Hove singularity (vHs) and the Fermi surface
nesting (FSN) Hohenadler and Batrouni (2019); Costa  _et al._ (2020). On the
other hand, for a pure Rashba metal, i.e. $t_{R}\neq 0$ and $g=0$, both FSN
and vHs are absent at the half-filling, due to the splitting of the spin
textures – similar to those presented in Fig. 1 (a) and (c). Such feature
leads to no divergence in the density of states at the Fermi level
$N(\omega=0)$ – consequently, no divergence in the noninteracting electronic
susceptibility. Therefore, within a weak coupling analysis of these two cases,
one should expect that the RSOC must suppress CDW formation. Although this may
be true for the weak coupling case, the scenario is less clear in the
intermediate and strong coupling regimes, where the Fermi surface does not
play a relevant role. In what follows, we investigate this problem at the
ground state, by performing analyses within both MF and CPT approaches, fixing
$\omega_{0}/t=1$ while varying $t_{R}$ and $g$.

##  III Mean-field theory

We start by discussing the mean-field approach to the Rashba-Holstein
Hamiltonian. To this end, we change the notation of Eq. (1), rewriting the
phonon creation and annihilation operators in terms of position and momentum
ones. For instance, by considering the pure Holstein model ($t_{R}=0$), its
Hamiltonian is

$\displaystyle\mathcal{H}_{\text{H}}=$
$\displaystyle-t\sum_{\langle\textbf{i,j}\rangle,\sigma}(c_{\textbf{i},\sigma}^{\dagger}c_{\textbf{j},\sigma}+h.c.)-\mu\sum_{\textbf{i},\sigma}n_{\textbf{i},\sigma}$
$\displaystyle-g\sqrt{\frac{2M\omega_{0}}{\hbar}}\sum_{\textbf{i},\sigma}n_{\textbf{i}}^{\sigma}\hat{X}_{\textbf{i}}+\sum_{\textbf{i}}\left(\dfrac{\hat{P}_{\textbf{i}}^{2}}{2M}+\dfrac{M\omega^{2}_{0}}{2}\hat{X}_{\textbf{i}}^{2}\right),$

with
$\hat{X}_{\bf{i}}=\sqrt{\dfrac{\hbar}{2M\omega_{0}}}(a_{\textbf{i}}+a_{\textbf{i}}^{\dagger})$
and $\hat{P}_{\bf{i}}=\mathit{i}\sqrt{\dfrac{\hbar
M\omega_{0}}{2}}(a_{\textbf{i}}^{\dagger}-a_{\textbf{i}})$ being the position
and momentum operators, respectively, of a harmonic oscillator in a given site
$\mathbf{i}$. Here, we define the mass of the harmonic oscillators ($M$) and
the lattice constant as unities.

Within a static mean-field approach, we assume that the kinetic energy of
lattice vibrations is negligible if compared with potential energy
($\hat{P}_{\textbf{i}}^{2}/2M\ll M\omega_{0}^{2}\hat{X}_{\textbf{i}}^{2}/2$),
i.e., we deal with a permanent distortion that breaks the translation symmetry
of the lattice. This assumption corresponds to the adiabatic regime Li  _et
al._ (2010); Alexandrov  _et al._ (1994), i.e. $\omega_{0}\to 0$ while keeping
$M\omega^{2}_{0}\equiv K$ finite. In this regime, the Hamiltonian is a
quadratic form for electron creation/annihilation operators, therefore its
ground state is exact by a mean-field approach.

Proceeding, we approximate the position operator by its expected value,
assuming the ansatz $\hat{X}_{\textbf{i}}\rightarrow\langle
X_{\bf{i}}\rangle=X_{0}+\cos({\bf{Q}}\cdot{\bf{i}})X_{1}$, with
$\textbf{Q}=(\pi,\pi)$. Notice that other ansatze involving arbitrary
wavevectors $\textbf{Q}=(q_{x},q_{y})$ can be explored, even in incommensurate
cases Lopes  _et al._ (2023). However, owing to the four-fold rotational
symmetry inherent to the square lattice, along with the nesting properties and
the commensurate filling, it is expected that the $(\pi,\pi)$ configuration is
the one that most significantly reduces the internal energy. Therefore, the MF
Holstein Hamiltonian becomes

$\displaystyle\mathcal{H}_{\text{H}}=$
$\displaystyle-t\sum_{\langle\textbf{i,j}\rangle,\sigma}(c_{\textbf{i},\sigma}^{\dagger}c_{\textbf{j},\sigma}+\text{H.c.})-\mu\sum_{\textbf{i},\sigma}n_{\textbf{i},\sigma}$
(3)
$\displaystyle-g\sqrt{\frac{2M\omega_{0}}{\hbar}}\sum_{\textbf{i},\sigma}n_{\textbf{i},\sigma}\big{[}X_{0}+X_{1}\cos({\bf{Q}}\cdot{\bf{i}})\big{]}$
$\displaystyle+\dfrac{NM\omega^{2}_{0}}{2}\left(X_{0}^{2}+X_{1}^{2}\right),$

where $N$ is the number of lattice sites.

Figure 1: The mean-field result for the density of states of the Rashba-
Holstein model for EPC (a) $g/t=0.2$, and (b) $g/t=0.8$, and its electronic
dispersion in (c) and (d), respectively. The results are for fixed RSOC
$t_{\rm R}/t=0.5$, with the color codes indicating the chiral spin-texture
polarization of bands.

The inclusion of the Rashba term changes the electron dispersion, and the
eventual mean-field Hamiltonian in reciprocal space reads

$\displaystyle\mathcal{H}_{\text{MF}}$
$\displaystyle=\sum_{\mathbf{k},\sigma}(\varepsilon_{\mathbf{k}}-\mu^{\prime})c_{\mathbf{k}\sigma}^{\dagger}c_{\mathbf{k}\sigma}^{\phantom{\dagger}}+\sum_{\mathbf{k},\sigma,\sigma^{\prime}}(\widehat{V}_{\mathbf{k}})_{\sigma\sigma^{\prime}}^{\phantom{\dagger}}c_{\mathbf{k}\sigma}^{\dagger}c_{\mathbf{k}\sigma^{\prime}}^{\phantom{\dagger}},$
(4)
$\displaystyle-g\sqrt{\frac{M\omega_{0}}{2\hbar}}X_{1}\sum_{\mathbf{k},\sigma}(c_{\mathbf{k}+\mathbf{Q}\sigma}^{\dagger}c_{\mathbf{k}\sigma}^{\phantom{\dagger}}+\text{H.c.})$
$\displaystyle+\dfrac{NM\omega^{2}_{0}}{2}\left(X_{0}^{2}+X_{1}^{2}\right),$

with $\varepsilon_{\mathbf{k}}=-2t(\cos{k_{x}}+\cos{k_{y}})$,
$\widehat{V}_{\mathbf{k}}\equiv
2t_{\text{R}}(\sigma_{x}\sin{k_{y}}-\sigma_{y}\sin{k_{x}})$, and
$\mu^{\prime}=\mu+g\sqrt{2M\omega_{0}/\hbar}X_{0}$. In a basis of Nambu
spinors
$\Phi^{\dagger}_{\mathbf{k}}=~{}(c_{\mathbf{k}\uparrow}^{\dagger},c_{\mathbf{k}\downarrow}^{\dagger},c_{\mathbf{k}+\mathbf{Q}\uparrow}^{\dagger},c_{\mathbf{k}+\mathbf{Q}\downarrow}^{\dagger})$,
we obtain a quadratic form representation of $\mathcal{H}_{MF}$,

$\displaystyle\mathcal{H}_{\text{MF}}=\frac{1}{2}\sum_{\mathbf{k}}\Phi^{\dagger}_{\mathbf{k}}\widehat{H}_{\mathbf{k}}\Phi_{\mathbf{k}}+\dfrac{NM\omega^{2}_{0}}{2}\left(X_{0}^{2}+X_{1}^{2}\right)$
(5)

where

$\displaystyle\widehat{H}_{\mathbf{k}}=\left(\begin{array}[]{cccc}\varepsilon_{\mathbf{k}}-\mu^{\prime}&(\widehat{V}_{\mathbf{k}})_{\uparrow\downarrow}&-\Delta
X&0\\\
(\widehat{V}_{\mathbf{k}})_{\downarrow\uparrow}&\varepsilon_{\mathbf{k}}-\mu^{\prime}&0&-\Delta
X\\\ -\Delta
X&0&\varepsilon_{\mathbf{k}+\mathbf{Q}}-\mu^{\prime}&(\widehat{V}_{\mathbf{k}+\mathbf{Q}})_{\uparrow\downarrow}\\\
0&-\Delta
X&(\widehat{V}_{\mathbf{k}+\mathbf{Q}})_{\uparrow\downarrow}&\varepsilon_{\mathbf{k}+\mathbf{Q}}-\mu^{\prime}\\\
\end{array}\right)$ (10)

with $\Delta X=g\sqrt{2M\omega_{0}/\hbar}X_{1}$. Notice that we added $X_{0}$
to the chemical potential $\mu^{\prime}$, since this field shifts the energy
level, while the $X_{1}$ field modulates the electronic distribution,
therefore being proportional to the amplitude of the CDW order parameter.
111Indeed, if we define $\langle
n_{i}\rangle=n+\cos({\bf{Q}}\cdot{\bf{i}})\delta n$, from Eq. (3) one obtains
$\delta n=\omega_{0}\sqrt{\frac{M\hbar\omega_{0}}{2}}\frac{X_{1}}{g}$.

Figure 2: The behavior of the ground state mean-field parameter $X_{1}$
(proportional to the CDW order parameter) as a function of the electron-phonon
coupling $g$, for different values of RSOC $t_{R}$’s.

The parameters $X_{0(1)}$ are determined self-consistently by minimizing the
Helmholtz free energy, through the diagonalization of Eq. (5). We use the
minimization solutions to describe the ground state properties of the system,
such as the density of state (DOS) and the electron dispersion, e.g., as
displayed in Fig. 1 for small and large values of EPC, while keeping
$t_{R}/t=0.5$ fixed. As presented in Figs. 1 (a) and (c), when $g/t=0.2$, the
RSOC splits the electronic spectra into two branches with well-defined
chiralities, $\langle\overleftrightarrow{\sigma}\rangle.({\bf
k}\times\hat{z})/|{\bf k}|=\pm 1$ for right-handed and left-handed. Such a
split shifts the energy levels and, consequently, the van Hove singularities,
leading to a reduction in the DOS at the half-filling, as shown in Fig. 1 (a).
On the other hand, for large EPC, a CDW phase emerges and creates a charge gap
at the half-filling, as displayed in Figs. 1 (b) and (d), for $g/t=0.8$. This
result provides a few hints that the RSOC is harmful to the CDW phase.

To further emphasize the effects of the RSOC on the CDW phase, Fig. 2 displays
the behavior of $X_{1}$ as a function of $g/t$, for different spin-orbit
couplings, in the ground state. For fixed values of $g/t$, where the
minimization of the energy leads to CDW, it is clear that the main effect of
the RSOC is to reduce $X_{1}$, which is deleterious to charge modulation. At
this point, a remark should be made: the asymptotic behavior to $X_{1}$ as a
function of $g/t$ may be an artifact of the mean-field solution. Then, we
define $X_{1}=10^{-3}$ as the lower bound limit for the CDW phase, which is
enough to reproduce $g_{c}\to 0$ when $t_{R}\to 0$, as expected for the pure
Holstein model at the half-filled square lattice. This lower bound limit
defines the critical region where the CDW phase is destroyed (or emerges).

Figure 3: The ground-state phase diagram of Holstein-Rashba model in the half-
filled square lattice. The black square symbols correspond to the mean-field
solution, while the red star ones denote the CPT results. The blue heatmap
indicates the amplitude of charge modulation in the CDW phase for the MFT
solution. The red curve is just a guide to the eye.

Repeating similar analyses as those of Fig. 2, we are able to present the
mean-field phase diagram of the Rashba-Holstein model at half-filling, as
shown in Fig. 3. Here, the CDW and Rashba metal phases are separated by the
curve going through the black squares, while the color-map scheme describes
the amplitude of the CDW charge modulation, $\delta n$. As complementary
information, Fig. 4 displays the critical CDW temperature $T_{\rm CDW}$, below
which the system opens a gap 222Again, notice that tail-like behavior of
$T_{\rm CDW}$ for large SOC is a mean-field artifact.. In line with the
previous discussion, the RSOC tends to suppress the CDW phase, decreasing both
the critical temperatures and the order parameter, eventually leading to a
Rashba metal.

Despite being obtained by a static mean-field approach, the phase diagram of
Fig. 3 is expected to provide the correct transition lines when $\omega_{0}\to
0$. However, $\omega_{0}/t=1$ is clearly away from this limit, which may
require less biased methodologies to investigate the Hamiltonian properties,
as discussed in the next Section.

##  IV Cluster Perturbation Theory

Now, we use the cluster perturbation theory, which is a convenient technique
to extract spectral weights and the electronic DOS of a correlated system.
Briefly speaking, the CPT consists of partitioning an infinite lattice into
clusters of finite size where the Hamiltonian is exactly solved, while
coupling the clusters by perturbation theory approaches. A pedagogical
introduction to this method can be found in Ref. Hohenadler  _et al._ , 2003;
Sénéchal  _et al._ , 2002, 2000; Sénéchal, 2010.

Figure 4: CDW critical temperature as a function of $t_{R}/t$ for different
values of $g/t$.

Within this approach, we decompose the Hamiltonian as

$\displaystyle H=H_{\text{cluster}}+T~{},$ (11)

where $H_{\text{cluster}}=K+H_{\text{int}}$ describes the operations inside a
given finite-sized cluster, with
$K=\sum_{\alpha,\beta}t_{\alpha,\beta}c_{\alpha}^{\dagger}c_{\beta}$ and
$H_{\text{int}}=-g\sum_{\textbf{i},\sigma}n_{\textbf{i},\sigma}(a_{\textbf{i}}+a_{\textbf{i}}^{\dagger})$
being the hopping and the EPC terms, respectively. The term
$T=\sum_{\alpha,\beta}V_{\alpha,\beta}c_{\alpha}^{\dagger}c_{\beta}$ is the
one that connects different clusters of the superlattice (i.e., intercluster
sites). Given this, one can write the Green’s functions as

$\displaystyle\textbf{G}(\omega)^{-1}=\mathbf{G}^{\prime}(\omega)^{-1}-\mathbf{V},$
(12)

where $\textbf{G}(\omega)$ and $\mathbf{G}^{\prime}(\omega)$ are the infinite
lattice and the (exact) finite-sized cluster Green’s functions, respectively.
As the cluster superlattice has a reduced Brillouin zone, one can connect any
wavevector ${\bf k}$ of the original Brillouin zone to a wavevector
${\bf\tilde{k}}$ in the reduced one by ${\bf k}={\bf\tilde{k}}+{\bf K}$, where
${\bf K}$ is an element of the reciprocal superlattice. Therefore, in the
reduced Brillouin zone, Eq. (12) is given by

$\displaystyle\textbf{G}(\textbf{\~{k}},\omega)^{-1}=\mathbf{G}^{\prime}(\omega)^{-1}-\mathbf{V}(\textbf{\~{k}})~{},$
(13)

which eventually leads to the CPT Green’s functions

$\displaystyle
G_{\text{CPT}}(\mathbf{k},\omega)=\frac{1}{L}\sum_{a,b}e^{-i\mathbf{k}(\mathbf{r}_{a}-\mathbf{r}_{b})}G_{a,b}(\mathbf{\tilde{k}},\omega)~{},$
(14)

with ${\bf r}_{a(b)}$ being the position of the intra-cluster sites.

The CPT Green’s functions take into account the leading order in the
interaction term, and allow one to obtain the spectral weight

$\displaystyle
A_{0}(\mathbf{k},\omega)=-\frac{1}{2\pi}\text{Im}[G_{\text{CPT}}(\mathbf{k},\omega)],$
(15)

and the DOS

$\displaystyle
N(\omega)=\frac{1}{N}\sum_{\mathbf{k}}A_{0}(\mathbf{k},\omega).$ (16)

In what follows, we use the DOS to identify the transition from a gapped CDW
phase to a Rashba metal.

Figure 5: The ground state energy of the Holstein model for $g/t=1.5$, in a
$2\times 2$ plaquette, as a function of the cutoff number of phonons per site
$N_{p}$. The red squares (blue circles) symbols denote the ground state
energies for the case with (without) a shift to symmetrize the electron-phonon
term.

At this point, it is worth discussing a few technical issues about our CPT
implementation. Due to the phonon degrees of freedom, the Hilbert space is
infinite, which demands a cutoff in the total number of phonons to perform the
exact diagonalization procedure. Among the different ways of establishing this
cutoff, we have defined a maximum limit, $N_{p}$, for the number of phonons
per site. Its value is determined through a systematic analysis of the ground
state energy $E_{g}$ as a function of $N_{p}$, as displayed in Fig. 5, for
fixed $g/t=1.5$ and $\omega_{0}/t=1$, in a $2\times 2$ plaquette. Notice that
$E_{g}(N_{p})$ goes asymptotically to the limit $N_{p}\to\infty$, so the error
from the cut-off may be disregarded for a given value of $N_{p}$. However, the
electron-phonon term described in Eq. (1) requires a substantial number of
phonons to attain the correct ground state, as evidenced by the behavior of
the blue circle symbols in Fig. 5. To circumvent this difficulty and keep
$N_{p}$ as small as possible, we perform a shift in the phonon operators
$a_{\mathbf{i}}\to\big{(}a_{\mathbf{i}}+g\langle n\rangle\big{)}$, with
$\langle n\rangle$ being the average electron density. We emphasize that such
a shift is equivalent to a change in the origin of the phonon position
operators $\hat{X}_{\mathbf{i}}$ to deal just with excitations around its mean
value. This drastically reduces the number of required phonons, as shown by
the behavior of the red square symbols in Fig. 5. Despite this improvement,
the greater the EPC, the larger the values of $N_{p}$, therefore we limit our
analyses to intermediate coupling strengths for clusters of $2\times 2$ sites.

Figure 6: The DOS of the Rashba-Holstein model for EPC (a) $g/t=1.2$ and (b)
$g/t=1.5$, for different values of RSOC $t_{R}$.

Now, using the CPT method we turn to discuss the effects of the RSOC on the
spectral properties of our charge-ordered system, for fixed $\omega_{0}/t=1$,
and at half-filling. Figure 6 (a) displays the behavior of $N(\omega)$ for
$g/t=1.2$ and different values of $t_{R}$. When $t_{R}=0$, the system has a
large charge gap around $\omega=0$, in line with our expectation of a CDW
phase for the Holstein model in the half-filled square lattice. Interestingly,
the size of the gap remains nearly unchanged for $t_{R}/t\lesssim 0.5$, which
shows that the CDW phase is robust against spin-orbit perturbations. However,
for larger RSOC, the charge gap drastically decreases, until $N(\omega\approx
0)$ exhibits a peak instead of a gap for $t_{R}/t\approx 1.25$. This specific
change in the behavior of $N(\omega)$ around $\omega\approx 0$ (and
$t_{R}/t\approx 1.25$) is a clear indication of the occurrence of a phase
transition from the CDW phase to a Rashba metal.

A very similar behavior occurs for different values of $g$, as displayed in
Fig. 6 (b), for $g/t=1.5$. The main difference between this case and the
previous one is the need for a larger RSOC to destroy the charge gap. In view
of this, here we define as critical points the values of $t_{R}$ at which
$N(\omega)$ exhibits a peak instead of a gap at the half-filling. These CPT
critical points are presented as red star symbols in the phase diagram of Fig.
3, with the red (dashed and solid) lines being just guides to the eye.

In order to further emphasize the phase transitions, we extract the charge gap
directly from Fig. 6 by defining the CDW order parameter, $\Delta_{\rm CDW}$,
as the energy difference between the two peaks away from $\omega=0$. The
behavior of $\Delta_{\rm CDW}$ as a function of $t_{R}$ is displayed in Fig.
7, for $g/t=1.2$ (red diamond symbols) and 1.5 (blue triangle symbols), from
which one may notice that the gap closing is consistent with a second order
phase transition. In particular, the suppression of the CDW phase seems to
occur at $t_{R(c)}=1.1\pm 0.1$ for $g/t=1.2$, and at $t_{R(c)}=2.1\pm 0.1$ for
$g/t=1.5$.

Figure 7: The CDW charge gap, $\Delta_{\rm CDW}$ as a function of the RSOC,
$t_{R}/t$, for fixed $\omega_{0}/t=1$.

There are clear differences between the CPT and MFT results, which may be
understood through two key remarks. First, and foremost, the CPT method takes
into account the effects of short-range electronic correlations, while the MFT
deals with local ones. That is, the MFT exhibits the points where the Rashba
metal becomes on-site correlated, while the CPT provides insights when these
correlations become long-ranged. Therefore, the region enclosed between the
MFT and CPT lines in Fig. 3 defines a correlated Rashba metal, i.e. metallic
(from CPT) but with on-site correlations (from MFT approach). As a second
remark, we recall that the mean-field solution is exact at the ground state
when $\omega_{0}\to 0$, thus the MFT line in Fig. 3 can also be interpreted as
as the boundary of the adiabatic limit. That is, if one reduces $\omega_{0}$,
the CPT transition line must be pushed into the MFT ones. Then, it provides
the behavior of the critical points when $\omega_{0}$ changes.

At this point, it is worth discussing the emergence of pairing correlations in
the absence of the CDW ones. As the occurrence of superconductivity is a
challenge for both methodologies, we analyze the s-wave pair correlation
functions Paiva  _et al._ (2004)

$\displaystyle P_{s}({\bf i},{\bf j})=\frac{1}{2}\langle b^{\dagger}_{{\bf
i}}b_{{\bf j}}+\text{H.c.}\rangle$ (17)

with $b^{\dagger}_{{\bf i}}=c^{\dagger}_{{\bf i}\uparrow}c^{\dagger}_{{\bf
i}\downarrow}$ ($b_{{\bf j}}=c_{{\bf j}\downarrow}c_{{\bf j}\uparrow}$),
restricting our analysis to nearest-neighbor sites. This function is computed
in real space, after partitioning the system in the cluster scheme described
above. Figure 8 (a) presents $P_{s}(|{\bf i}-{\bf j}|=1)$ as a function of
$t_{\rm R}$ for $g/t=1.2$ and 1.5, from which we see an enhancement of the
pairing correlations. This result provides us with an interesting feature: the
derivative of the nearest-neighbor pair correlation function ($dP_{s}/dt_{R}$)
exhibits peaks for values very close to the CDW critical points, as shown in
Fig. 8 (b). Such a behavior is not a coincidence; as charge-charge
correlations compete with pairing ones, then it is expected that $P_{s}({\bf
i},{\bf j})$ should be enhanced around the points where the CDW gap closes,
independently whether the ground state is superconducting or not. Therefore,
although this result is not enough to support the existence of
superconductivity, it reinforces our previous findings about the nature of the
charge gap being a CDW phase, and about the location of the CPT critical
points presented in the phase diagram of Fig. 3.

Figure 8: The s-wave pair correlation function $P_{s}({\bf i},{\bf j})$ as a
function of RSOC and its numerical derivative $\text{d}P_{s}({\bf i},{\bf
j})/\text{d}t_{\rm R}$ for two values of EPC. The correlation function is
calculated between pairs at a site ${\bf i}$ and its nearest-neighbor ${\bf
j}$ symbolized here by $({\bf i},{\bf j})=(R_{1})$.

Last, but not least, we recall that the spin-orbit coupling has been
investigated in different contexts of strongly correlated systems, in
particular for those involving electron-electron interactions. Then, it is
imperative to compare these cases with ours, highlighting both their
similarities and differences. For instance, for the repulsive Hubbard model in
the half-filled square lattice, the inclusion of the RSOC leads to the
emergence of novel magnetic phases, which may include spiral and striped ones;
away from the half-filling, ferromagnetism may appear Zhang  _et al._ (2015);
Brosco and Capone (2020); Kennedy  _et al._ (2022); Wang and Zhu (2023);
Kawano and Hotta (2023); Hodt  _et al._ (2023). This occurs due to the spin-
flip term in the RSOC, which disfavors the staggered antiferromagnetic ground
state of the half-filling, suppressing the charge gap, while keeping gapless
spin excitations. Due to the spin dependence in the interaction term in the
Hubbard model, the CPT spectral properties for this case exhibit a clear
dependence of the chiral bands and DOS as a function of the Hubbard
interaction strength $U$ Brosco and Capone (2020). It is substantially
different from our present study, since the interaction term in the Rashba-
Holstein model is not spin-dependent, with the spin and charge gaps
emerging/disappearing together. Furthermore, noteworthy changes manifest in
the attractive Hubbard model when subjected to the RSOC. This model exhibits
s-wave superconductivity for any filling Paiva  _et al._ (2004); Moreo and
Scalapino (1991); Fontenele  _et al._ (2022), whose critical temperatures are
significantly increased with the inclusion of the RSOC Tang  _et al._ (2014);
Rosenberg  _et al._ (2017); Ptok  _et al._ (2018). This enhancement results
from the suppression of charge-charge correlations, in line with our findings
for the Rashba-Holstein model. Given that the antiadiabatic regime
($\omega_{0}\to\infty$) of the Holstein one can be mapped onto the attractive
Hubbard model, we anticipate the emergence of superconductivity in our
scenario. However, achieving this would require the application of unbiased
methodologies, such as quantum Monte Carlo methods, which is beyond the scope
of the present study.

##  V Conclusions

In this work, we have studied the impact of the Rashba spin-orbit coupling to
a charge-ordered system, as the Holstein model. To this end, we investigated
the Rashba-Holstein model using two methodologies: (i) a static mean-field
approach, and (ii) the cluster perturbation theory – the latter provides the
spectral functions of an infinite system with short-range correlations, i.e.
beyond a Hartree-Fock approximation. Given this, our main findings are
summarized in the phase diagram of Fig. 3, which presents the CPT and MFT
results, from which one notices that the RSOC is harmful to the CDW phase.
This deleterious feature of the RSOC could be understood from the diminishing
DOS at the half-filling when $t_{R}$ increases, from a weak coupling point-of-
view; or by a large energetic price to pay to avoid spin-flip processes
(within the CDW phase), from a strong coupling point-of-view. Furthermore, as
shown in Figs. 6 and 7, the transition from the CDW phase to a correlated
Rashba metal seems to be continuous, with the transition line being pushed to
larger values of $t_{R}$ when $\omega_{0}\to 0$. Even though the critical
phenomena properties (such as critical exponents) are beyond the scope of this
work, our phase diagram $g\times t_{\rm R}$ provides accurate critical points,
whose fundamental properties may shed light into the growing class of Rashba
materials Hanakata  _et al._ (2017, 2018); Chen  _et al._ (2021).

## ACKNOWLEDGMENTS

The authors are grateful to R.R. dos Santos for enlightening discussions and
suggestions. Financial support from Fundação Carlos Chagas Filho de Amparo à
Pesquisa do Estado do Rio de Janeiro, grant number E-26/200.258/2023 -
SEI-260003/000623/2023 (N.C.C.); and from CNPq grant number 313065/2021-7
(N.C.C.) is gratefully acknowledged. We also acknowledge support from INCT-IQ.
T.P.C. acknowledge support from CAPES (Brazil). S.A.S.Jr. acknowledges support
from CNPq (Brazil) and Fundação Carlos Chagas Filho de Amparo à Pesquisa do
Estado do Rio de Janeiro process SEI-260003/006232/2023. R.A.F. acknowledges
support from Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do
Rio de Janeiro and CNPq (Brazil).

## References

  * Mas-Ballesté  _et al._ (2011) Rubén Mas-Ballesté, Cristina Gómez-Navarro, Julio Gómez-Herrero, and Félix Zamora, “2D materials: to graphene and beyond,” Nanoscale 3, 20–30 (2011).
  * Gupta  _et al._ (2015) Ankur Gupta, Tamilselvan Sakthivel, and Sudipta Seal, “Recent development in 2D materials beyond graphene,” Progress in Materials Science 73, 44–126 (2015).
  * Novoselov  _et al._ (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. Castro Neto, “2D materials and van der Waals heterostructures,” Science 353, aac9439 (2016), https://www.science.org/doi/pdf/10.1126/science.aac9439 .
  * Ahn (2020) Ethan C Ahn, “2D materials for spintronic devices,” npj 2D Materials and Applications 4, 17 (2020).
  * Kotov  _et al._ (2012) Valeri N. Kotov, Bruno Uchoa, Vitor M. Pereira, F. Guinea, and A. H. Castro Neto, “Electron-Electron Interactions in Graphene: Current Status and Perspectives,” Rev. Mod. Phys. 84, 1067–1125 (2012).
  * Zhang  _et al._ (2019) Y.-X. Zhang, W.-T. Chiu, N. C. Costa, G. G. Batrouni, and R. T. Scalettar, “Charge Order in the Holstein Model on a Honeycomb Lattice,” Phys. Rev. Lett. 122, 077602 (2019).
  * Paiva  _et al._ (2005) Thereza Paiva, R. T. Scalettar, W. Zheng, R. R. P. Singh, and J. Oitmaa, “Ground-state and finite-temperature signatures of quantum phase transitions in the half-filled Hubbard model on a honeycomb lattice,” Phys. Rev. B 72, 085123 (2005).
  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. 81, 109–162 (2009).
  * Cao  _et al._ (2018) Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43–50 (2018).
  * Gong  _et al._ (2017) Cheng Gong, Lin Li, Zhenglu Li, Huiwen Ji, Alex Stern, Yang Xia, Ting Cao, Wei Bao, Chenzhe Wang, Yuan Wang, Z. Q. Qiu, R. J. Cava, Steven G. Louie, Jing Xia, and Xiang Zhang, “Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals,” Nature 546, 265–269 (2017).
  * Huang  _et al._ (2017) Bevin Huang, Genevieve Clark, Efrén Navarro-Moratalla, Dahlia R. Klein, Ran Cheng, Kyle L. Seyler, Ding Zhong, Emma Schmidgall, Michael A. McGuire, David H. Cobden, Wang Yao, Di Xiao, Pablo Jarillo-Herrero, and Xiaodong Xu, “Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit,” Nature 546, 270–273 (2017).
  * Manzeli  _et al._ (2017) S. Manzeli, D. Ovchinnikov, D. Pasquier, O.V. Yazyev, and A. Kis, “2D transition metal dichalcogenides,” Nature Reviews Materials 2, 17033 (2017).
  * Grüner (1988) G. Grüner, “The dynamics of charge-density waves,” Rev. Mod. Phys. 60, 1129–1181 (1988).
  * Gorḱov and Grüner (1989) L. P. Gorḱov and G. Grüner,  _Charge density waves in solids_ (Elsevier, Amsterdam, 1989).
  * Sugawara  _et al._ (2015) Katsuaki Sugawara, Yuki Nakata, Ryota Shimizu, Patrick Han, Taro Hitosugi, Takafumi Sato, and Takashi Takahashi, “Unconventional Charge-Density-Wave Transition in Monolayer 1T-TiSe2,” ACS Nano 10, 1341–1345 (2015).
  * Lin  _et al._ (2020) Dongjing Lin, Shichao Li, Jinsheng Wen, Helmuth Berger, László Forró, Huibin Zhou, Shuang Jia, Takashi Taniguchi, Kenji Watanabe, Xiaoxiang Xi, and Mohammad Saeed Bahramy, “Patterns and driving forces of dimensionality-dependent charge density waves in 2H-type transition metal dichalcogenides,” Nature Communications 11, 2406 (2020).
  * Ge and Liu (2012) Yizhi Ge and Amy Y. Liu, “Effect of dimensionality and spin-orbit coupling on charge-density-wave transition in 2H-TaSe2,” Phys. Rev. B 86, 104101 (2012).
  * Lian  _et al._ (2022) Chao-Sheng Lian, Christoph Heil, Xiaoyu Liu, Chen Si, Feliciano Giustino, and Wenhui Duan, “Intrinsic and doping-enhanced superconductivity in monolayer $1H\text{$-$}{\mathrm{TaS}}_{2}$: Critical role of charge ordering and spin-orbit coupling,” Phys. Rev. B 105, L180505 (2022).
  * Nitta  _et al._ (1997) Junsaku Nitta, Tatsushi Akazaki, Hideaki Takayanagi, and Takatomo Enoki, “Gate Control of Spin-Orbit Interaction in an Inverted I${\mathrm{n}}_{0.53}$G${\mathrm{a}}_{0.47}$As/I${\mathrm{n}}_{0.52}$A${\mathrm{l}}_{0.48}$As Heterostructure,” Phys. Rev. Lett. 78, 1335–1338 (1997).
  * Ast  _et al._ (2007) Christian R. Ast, Jürgen Henk, Arthur Ernst, Luca Moreschini, Mihaela C. Falub, Daniela Pacilé, Patrick Bruno, Klaus Kern, and Marco Grioni, “Giant Spin Splitting through Surface Alloying,” Phys. Rev. Lett. 98, 186807 (2007).
  * Rotenberg  _et al._ (1999) Eli Rotenberg, J. W. Chung, and S. D. Kevan, “Spin-Orbit Coupling Induced Surface Band Splitting in Li/W(110) and Li/Mo(110),” Phys. Rev. Lett. 82, 4066–4069 (1999).
  * LaShell  _et al._ (1996) S. LaShell, B. A. McDougall, and E. Jensen, “Spin Splitting of an Au(111) Surface State Band Observed with Angle Resolved Photoelectron Spectroscopy,” Phys. Rev. Lett. 77, 3419–3422 (1996).
  * Min  _et al._ (2006) Hongki Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, Leonard Kleinman, and A. H. MacDonald, “Intrinsic and Rashba spin-orbit interactions in graphene sheets,” Phys. Rev. B 74, 165310 (2006).
  * Cysne  _et al._ (2018a) Tarik P. Cysne, Aires Ferreira, and Tatiana G. Rappoport, “Crystal-field effects in graphene with interface-induced spin-orbit coupling,” Phys. Rev. B 98, 045407 (2018a).
  * Cysne  _et al._ (2018b) Tarik P. Cysne, Jose H. Garcia, Alexandre R. Rocha, and Tatiana G. Rappoport, “Quantum Hall effect in graphene with interface-induced spin-orbit coupling,” Phys. Rev. B 97, 085413 (2018b).
  * Hanakata  _et al._ (2017) Paul Z. Hanakata, A. S. Rodin, Alexandra Carvalho, Harold S. Park, David K. Campbell, and A. H. Castro Neto, “Two-dimensional square buckled Rashba lead chalcogenides,” Phys. Rev. B 96, 161401 (2017).
  * Hanakata  _et al._ (2018) Paul Z. Hanakata, A. S. Rodin, Harold S. Park, David K. Campbell, and A. H. Castro Neto, “Strain-induced gauge and Rashba fields in ferroelectric Rashba lead chalcogenide $\mathrm{Pb}X$ monolayers ($X=\mathrm{S}$, Se, Te),” Phys. Rev. B 97, 235312 (2018).
  * Chen  _et al._ (2021) Jiajia Chen, Kai Wu, Wei Hu, and Jinlong Yang, “Spin–Orbit Coupling in 2D Semiconductors: A Theoretical Perspective,” The Journal of Physical Chemistry Letters 12, 12256–12268 (2021).
  * Pesin and Balents (2010) D. Pesin and L. Balents, “Mott physics and band topology in materials with strong spin-orbit interaction,” Nature Physics 6, 376–381 (2010).
  * Witczak-Krempa  _et al._ (2014) W. Witczak-Krempa, G. Chen, Y.B. Kim, and L. Balents, “Correlated quantum phenomena in the strong spin-orbit regime,” Annual Review of Condensed Matter Physics 5, 57–82 (2014).
  * Rau  _et al._ (2016) J.G. Rau, E.K.-H. Lee, and H.-Y. Kee, “Spin-orbit physics giving rise to novel phases in correlated systems: Iridates and related materials,” Annual Review of Condensed Matter Physics 7, 195–221 (2016).
  * Holstein (1959) T Holstein, “Studies of polaron motion,” Annals of Physics 8, 325–342 (1959).
  * Costa  _et al._ (2020) Natanael C. Costa, Kazuhiro Seki, Seiji Yunoki, and Sandro Sorella, “Phase diagram of the two-dimensional Hubbard-Holstein model,” Communications Physics 3, 80 (2020).
  * Dee  _et al._ (2019) P. M. Dee, K. Nakatsukasa, Y. Wang, and S. Johnston, “Temperature-filling phase diagram of the two-dimensional Holstein model in the thermodynamic limit by self-consistent Migdal approximation,” Phys. Rev. B 99, 024514 (2019).
  * Bradley  _et al._ (2021) Owen Bradley, George G. Batrouni, and Richard T. Scalettar, “Superconductivity and charge density wave order in the two-dimensional Holstein model,” Phys. Rev. B 103, 235104 (2021).
  * Cohen-Stead  _et al._ (2019) B. Cohen-Stead, N. C. Costa, E. Khatami, and R. T. Scalettar, “Effect of strain on charge density wave order in the Holstein model,” Phys. Rev. B 100, 045125 (2019).
  * Araújo  _et al._ (2022) Maykon V. Araújo, José P. de Lima, Sandro Sorella, and Natanael C. Costa, “Two-dimensional $t-{t}^{{}^{\prime}}$ Holstein model,” Phys. Rev. B 105, 165103 (2022).
  * Li  _et al._ (2015) Shaozhi Li, E. A. Nowadnick, and S. Johnston, “Quasiparticle properties of the nonlinear Holstein model at finite doping and temperature,” Phys. Rev. B 92, 064301 (2015).
  * Dee  _et al._ (2020) Philip M. Dee, Jennifer Coulter, Kevin G. Kleiner, and Steven Johnston, “Relative importance of nonlinear electron-phonon coupling and vertex corrections in the Holstein model,” Communications Physics 3, 145 (2020).
  * Paleari  _et al._ (2021) G. Paleari, F. Hébert, B. Cohen-Stead, K. Barros, RT. Scalettar, and G. G. Batrouni, “Quantum Monte Carlo study of an anharmonic Holstein model,” Phys. Rev. B 103, 195117 (2021).
  * Xiao  _et al._ (2021) B. Xiao, N. C. Costa, E. Khatami, G. G. Batrouni, and R. T. Scalettar, “Charge density wave and superconductivity in the disordered Holstein model,” Phys. Rev. B 103, L060501 (2021).
  * Sousa-Júnior  _et al._ (2023) Sebastião dos Anjos Sousa-Júnior, Raimundo R. dos Santos, and Natanael C. Costa, “Magnetic impurities in a charge-ordered background,” Phys. Rev. B 107, 075140 (2023).
  * Costa  _et al._ (2018) N. C. Costa, T. Blommel, W.-T. Chiu, G. Batrouni, and R. T. Scalettar, “Phonon Dispersion and the Competition between Pairing and Charge Order,” Phys. Rev. Lett. 120, 187003 (2018).
  * Hohenadler and Batrouni (2019) M. Hohenadler and G. G. Batrouni, “Dominant charge density wave correlations in the Holstein model on the half-filled square lattice,” Phys. Rev. B 100, 165114 (2019).
  * Li  _et al._ (2010) Zhou Li, D. Baillie, C. Blois, and F. Marsiglio, “Ground-state properties of the Holstein model near the adiabatic limit,” Phys. Rev. B 81, 115114 (2010).
  * Alexandrov  _et al._ (1994) A. S. Alexandrov, V. V. Kabanov, and D. K. Ray, “From electron to small polaron: An exact cluster solution,” Phys. Rev. B 49, 9915–9923 (1994).
  * Lopes  _et al._ (2023) Nei Lopes, Daniel Reyes, Natanael C. Costa, Mucio A. Continentino, and Christopher Thomas, “Incommensurate charge density wave in multiband intermetallic systems exhibiting competing orders,” Phys. Rev. B 107, 205141 (2023).
  * Note (1) Indeed, if we define $\langle n_{i}\rangle=n+\mathop{cos}\nolimits({\bf{Q}}\cdot{\bf{i}})\delta n$, from Eq. (3) one obtains $\delta n=\omega_{0}\sqrt{\frac{M\hbar\omega_{0}}{2}}\frac{X_{1}}{g}$.
  * Note (2) Again, notice that tail-like behavior of $T_{\rm CDW}$ for large SOC is a mean-field artifact.
  * Hohenadler  _et al._ (2003) Martin Hohenadler, Markus Aichhorn, and Wolfgang von der Linden, “Spectral function of electron-phonon models by cluster perturbation theory,” Phys. Rev. B 68, 184304 (2003).
  * Sénéchal  _et al._ (2002) David Sénéchal, Danny Perez, and Dany Plouffe, “Cluster perturbation theory for Hubbard models,” Phys. Rev. B 66, 075129 (2002).
  * Sénéchal  _et al._ (2000) D. Sénéchal, D. Perez, and M. Pioro-Ladrière, “Spectral Weight of the Hubbard Model through Cluster Perturbation Theory,” Phys. Rev. Lett. 84, 522–525 (2000).
  * Sénéchal (2010) David Sénéchal, “An introduction to quantum cluster methods,” (2010), arXiv:0806.2690 [cond-mat.str-el] .
  * Paiva  _et al._ (2004) Thereza Paiva, Raimundo R. dos Santos, R. T. Scalettar, and P. J. H. Denteneer, “Critical temperature for the two-dimensional attractive Hubbard model,” Phys. Rev. B 69, 184501 (2004).
  * Zhang  _et al._ (2015) Xin Zhang, Wei Wu, Gang Li, Lin Wen, Qing Sun, and An-Chun Ji, “Phase diagram of interacting Fermi gas in spin–orbit coupled square lattices,” New Journal of Physics 17, 073036 (2015).
  * Brosco and Capone (2020) Valentina Brosco and Massimo Capone, “Rashba-metal to Mott-insulator transition,” Phys. Rev. B 101, 235149 (2020).
  * Kennedy  _et al._ (2022) Welberth Kennedy, Sebastião dos Anjos Sousa-Júnior, Natanael C. Costa, and Raimundo R. dos Santos, “Magnetism and metal-insulator transitions in the Rashba-Hubbard model,” Phys. Rev. B 106, 165121 (2022).
  * Wang and Zhu (2023) Xuanlan Wang and Wei Zhu, “Intertwined charge and spin orders through the exchange-interaction on a square lattice,” Physics Letters A 461, 128652 (2023).
  * Kawano and Hotta (2023) Masataka Kawano and Chisa Hotta, “Phase diagram of the square-lattice Hubbard model with Rashba-type antisymmetric spin-orbit coupling,” Phys. Rev. B 107, 045123 (2023).
  * Hodt  _et al._ (2023) Erik Wegner Hodt, Jabir Ali Ouassou, and Jacob Linder, “Transient dynamics and quantum phase diagram for the square lattice Rashba-Hubbard model at arbitrary hole doping,” Phys. Rev. B 107, 224427 (2023).
  * Moreo and Scalapino (1991) A. Moreo and D. J. Scalapino, “Two-dimensional negative-U Hubbard model,” Phys. Rev. Lett. 66, 946–948 (1991).
  * Fontenele  _et al._ (2022) Rodrigo A. Fontenele, Natanael C. Costa, Raimundo R. dos Santos, and Thereza Paiva, “Two-dimensional attractive Hubbard model and the BCS-BEC crossover,” Phys. Rev. B 105, 184502 (2022).
  * Tang  _et al._ (2014) Ho-Kin Tang, Xiaosen Yang, Jinhua Sun, and Hai-Qing Lin, “Berezinskii-Kosterlitz-Thoules phase transition of spin-orbit coupled Fermi gas in optical lattice,” Europhysics Letters 107, 40003 (2014).
  * Rosenberg  _et al._ (2017) Peter Rosenberg, Hao Shi, and Shiwei Zhang, “Ultracold Atoms in a Square Lattice with Spin-Orbit Coupling: Charge Order, Superfluidity, and Topological Signatures,” Phys. Rev. Lett. 119, 265301 (2017).
  * Ptok  _et al._ (2018) Andrzej Ptok, Karen Rodríguez, and Konrad Jerzy Kapcia, “Superconducting monolayer deposited on substrate: Effects of the spin-orbit coupling induced by proximity effects,” Phys. Rev. Mater. 2, 024801 (2018).