††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Edge zeros and boundary spinons in topological Mott insulators

Niklas Wagner  Institut für Theoretische Physik und Astrophysik and Würzburg-
Dresden Cluster of Excellence ct.qmat, Universität Würzburg, 97074 Würzburg,
Germany  Daniele Guerci  Center for Computational Quantum Physics, Flatiron
Institute, New York, NY 10010, USA  Andrew J. Millis  Department of Physics,
Columbia University, New York, NY 10027, USA  Center for Computational Quantum
Physics, Flatiron Institute, New York, NY 10010, USA  Giorgio Sangiovanni
<EMAIL_ADDRESS>Institut für Theoretische Physik und
Astrophysik and Würzburg-Dresden Cluster of Excellence ct.qmat, Universität
Würzburg, 97074 Würzburg, Germany

(August 30, 2024)

###### Abstract

We introduce a real-space slave rotor theory of the physics of topological
Mott insulators, using the Kane-Mele-Hubbard model as an example, and use it
to show that a topological gap in the Green function zeros corresponds to a
gap in the bulk spinon spectrum and that a zero edge mode corresponds to a
spinon edge mode. We then consider an interface between a topological Mott
insulator and a conventional topological insulator showing how the spinon edge
mode of the topological Mott insulator combines with the spin part of the
conventional electron topological edge state leaving a non-Fermi liquid edge
mode described by a gapless propagating holon and gapped spinon state. Our
work demonstrates the physical meaning of Green function zeros and shows that
interfaces between conventional and Mott topological insulators are a rich
source of new physics.

Introduction – A conventional metal is characterized by gapless low energy
electronic quasiparticles with charge $e$ and spin $1/2$ and an electronic
Green function with poles that define the Fermi surface. A conventional
topological insulator is characterized by a gap to fermionic (or other)
excitations and a gapless fermionic edge state. A strongly correlated Mott
insulator may exhibit fractionalized quasiparticles (spinons) with charge zero
but spin $1/2$ and also an electronic Green function characterized with a
surface of zeros that defines a “Luttinger surface” that generalizes the
concept of the Fermi surface to the Mott state [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12].

Recently it has been shown that Green function zeros play a crucial role also
in strongly interacting topological Mott insulators [13, 14, 15, 16, 17, 18,
19]. In particular, it is possible to define an interacting topological
invariant [20, 21, 22, 23, 24, 25], in terms of the fully interacting single
particle Green function in a manner that treats poles and zeros of the Green
function on equal footing and, moreover, extends the bulk-boundary
correspondence to the strongly interacting regime simply replacing edge poles
with zeros [23, 26, 15]. Recently, some of us and others conjectured [15],
based on analysis of one dimensional models, that at the interface between a
Mott and conventional topological insulator the gapless edge zero band could
annihilate the edge state of the conventional topological insulator. However,
the physical meaning of the edge zeros has remained unclear because the zeros
do not directly correspond to any physical excitation and the topological
invariant is not directly connected to the Hall conductance of the strongly
interacting system [17, 27, 16, 28].

Figure 1: Illustration of the spin and charge separation at an interface
between a topological Mott insulator (TMI) and weakly-interacting topological
insulator (TI). The edge state of the TI consists of a spin and charge
contribution whereas the TMI exhibits only a spinon edge state. Along the
interface the spinon channels compensate leaving only a charge channel. Note
that only one spin channel is shown.

In this work we use a real-space auxiliary particle method to connect the
properties of Green function zeros of a Mott insulator to the properties of
spinon excitations, in particular demonstrating a one-to-one correspondence
between bulk gaps in the Green function zeros and in the spinon spectrum and
between gapless boundary zeros and gapless boundary spinons. Moreover we give
a new perspective on the physics emerging at the interface [29] between a
topological Mott insulator [30, 31] and a conventional topological insulator,
showing (see Fig. 1) how the gapless spinon [32] (Green function zero) of the
topological Mott insulator cancels the spinon portion of the conventional
topological edge state, leaving a spinless holon excitation on the
conventional topological insulator side.

Model and Methods – We consider the interacting Kane-Mele Hamiltonian
$H=H_{t}+H_{\rm SOC}+H_{int}$ [33], which consists of electrons moving on
sites of a lattice with two sublattices per unit cell and subject to a local
on-site interaction. This model exhibits both metallic, conventional insulator
and Mott insulating phases, which are topological for $\lambda\neq 0$. This
model is of current interest as a model of transition metal dichalcogenides
moiré semiconductors [34, 35]. The Hamiltonian terms are:

$\begin{split}H_{t}&=-t\sum_{\langle\mathbf{r}\in A,\mathbf{r}^{\prime}\in
B\rangle}f^{A\dagger}_{\mathbf{r}}f^{B}_{\mathbf{r}^{\prime}}+h.c.,\\\ H_{\rm
SOC}&=\lambda\sum_{s=A,B}\sum_{\langle\mathbf{r},\mathbf{r}^{\prime}\rangle\in
s}f^{s\dagger}_{\mathbf{r}}e^{i\frac{\pi}{2}\nu_{\mathbf{r},\mathbf{r}^{\prime}}\sigma^{z}}f^{s}_{\mathbf{r}^{\prime}},\end{split}$
(1)

where
$f^{s}_{\mathbf{r}}=[f^{s}_{\mathbf{r}\uparrow},f^{s}_{\mathbf{r}\downarrow}]^{T}$
is a two-dimensional spinor, ${s=A,B}$ refers to the two inequivalent
sublattices, $\sigma^{z}$ is the Pauli matrix in the spin space and
$\nu_{\mathbf{r},\mathbf{r}^{\prime}}=+1(-1)$ if the path connecting
$\mathbf{r}$ to $\mathbf{r}^{\prime}$ winds counterclockwise (clockwise). In
this study we do not consider the effect of a staggered potential
distinguishing the two sublattices. The Coulomb interaction is given by the
on-site Hubbard term:

$H_{\rm int}=\frac{U}{2}\sum_{s=A,B}\sum_{\mathbf{r}\in
s}\left(\sum_{\sigma}n^{s}_{\mathbf{r}\sigma}-1\right)^{2},$ (2)

where $n_{\mathbf{r}\sigma}$ is the electron number at site $\mathbf{r}$ with
spin $\sigma$. At $\lambda=0$ the model has two important symmetries: the
discrete rotation $C_{3z}$ and the product of spinless time reversal and
inversion $C_{2z}\mathcal{T}$ which imply a gapless (Dirac) spectrum for the
weakly correlated half filled metallic phase and a gapless (Dirac) spinon
spectrum in the Mott insulating phase. $\lambda\neq 0$ breaks
$C_{2z}\mathcal{T}$, gapping the bulk electron (weakly correlated) or spinon
(strongly correlated) spectra in the half filled case [36, 33].

We are interested in studying the model for large interaction $U\gg t$ in the
paramagnetic Mott insulating regime considered also in Ref. [37, 38, 39]. We
employ the slave-rotor technique which consists of introducing the O(2) rotor
[40, 41, 42]:

$f^{s}_{\mathbf{r}\sigma}=e^{i\theta^{s}_{\mathbf{r}}}\psi^{s}_{\mathbf{r}\sigma},$
(3)

where the phase degree of freedom $\theta^{s}_{\mathbf{r}}$ is conjugate to
the local charge and $\psi^{s}_{\mathbf{r}\sigma}$ is an auxiliary fermion
(spinon). The enlarged Hilbert space of rotors and spinons is constrained to
$L^{s}_{\mathbf{r}}-\sum_{\sigma}\psi^{s\dagger}_{\mathbf{r}\sigma}\psi^{s}_{\mathbf{r}\sigma}+1=0$.
We treat the theory in a saddle point approximation justified in a large spin
degeneracy limit. Technical details of the calculations are given in the
Supplementary Material [43].

Although more accurate methods exist, this approach captures non perturbative
phenomena such as the Mott transition [42, 44] and offers a transparent
physical interpretation of the correlated state. Moreover, the method is
simple enough to allow simulations with different boundary conditions, that is
crucial for establishing the connection between bulk topological invariants
and edge modes of topological Mott insulators. Within the self-consistent
gaussian approximation, formally derived in the large-$M$ limit [42], the
slave rotor method gives a non-trivial description of the Mott regime,
including dispersive Hubbard bands, spinon excitations with hopping scaling as
$\sim t^{2}/U$ and preserving many qualitative properties of the electronic
spectral function [40, 41, 42, 44]. This both allows the study of systems with
an unprecedented level of reliability and enables studies of the finite-size
geometries needed for a detailed analysis of the bulk-edge correspondence in
the strongly correlated Mott regime.

Results –  We focus on the electron Green function
$G_{\sigma}(\mathbf{k},\omega+i0^{+})$ which in the slave rotor representation
is a composite correlation function involving the spinon ($\psi^{s}$) and
rotor ($\theta$) fields as
$G^{ss^{\prime}}_{\sigma}(\mathbf{r}-\mathbf{r}^{\prime},\tau)=-\langle
T_{\tau}\left(e^{i(\theta^{s}_{\mathbf{r}}(\tau)-\theta^{s^{\prime}}_{\mathbf{r}^{\prime}})}\psi^{s}_{\mathbf{r}\sigma}(\tau)\psi^{s^{\prime}\dagger}_{\mathbf{r}^{\prime}\sigma}\right)\rangle$.
In the saddle-point approximation the Green function takes a simple form with
an analytical expression detailed in the Supplementary Material [43] (see also
Ref. [45, 46] for similar derivations):

$\begin{split}G^{ss^{\prime}}_{\sigma}(\mathbf{k},i\omega)&=Z^{ss^{\prime}}G^{ss^{\prime}}_{\sigma,\psi}(\mathbf{k},i\omega)\\\
+\frac{T}{N}&\sum_{\mathbf{q},i\nu}G^{ss^{\prime}}_{\sigma,\psi}(\mathbf{k}-\mathbf{q},i\omega-i\nu)G^{ss^{\prime}}_{X}(\mathbf{q},i\nu),\end{split}$
(4)

where $\omega=i\pi(2n+1)T$, $\nu=i\pi 2nT$ with $T$ temperature and $G_{\psi}$
and $G_{X}$ are the Green function of the spinon and the rotor, respectively.
Finally, $Z^{ss^{\prime}}$ is the quasiparticle weight which vanishes in the
Mott insulator.

The left panel of Fig. 2 a) shows $|\det
G_{\sigma}(\mathbf{k},\omega+i0^{+})|$ for a large interaction and a non-zero
spin orbit coupling parameter $\lambda$, indicating the poles (blue) and zeros
(red) of the interacting single-particle Green function, along with the
spinons (black). We notice that the zeros reside in the middle of the Hubbard
gap, the two different branches are gapped and their dispersion resembles that
of the spinons.

Figure 2: a) Left panel: color plot of the determinant of the slave-rotor
Green function for $\lambda/t=0.5$ and $U/t=10$ (blue and red indicate poles
and zeros of the Green function respectively) along with spinon dispersion
multiplied by 10 (black). Right panel: momentum-summed spectral weight
$A_{\uparrow}(\omega)=\sum_{\mathbf{k}}\Tr
A_{\uparrow}(\mathbf{k},\omega)/N_{s}$ with $N_{s}$ number of unit cells. b)
Zoom on the $\mathbf{K}$ point for different $\lambda$. Black lines show the
spinon dispersion. The Chern number $N_{\uparrow}$ for spin up associated to
the zeros is given in each plot title.

The Kane-Mele Hubbard model exhibits a topological transition as $\lambda$ is
tuned through $\lambda=0$. The two phases may be distinguished via the spin-
dependent topological invariant [20, 21, 47, 22, 23, 24, 17]:

$N_{\sigma}=\frac{\varepsilon_{\mu\nu\rho}}{24\pi^{2}}\int_{\omega\mathbf{k}}\Tr\left[G_{\sigma}\partial_{\mu}G^{-1}_{\sigma}G_{\sigma}\partial_{\nu}G^{-1}_{\sigma}G_{\sigma}\partial_{\rho}G^{-1}_{\sigma}\right],$
(5)

where $\int_{\omega\mathbf{k}}\equiv\int d\omega\int_{\rm BZ}d^{2}\mathbf{k}$
and $G$ is the interacting Green function. In the non-interacting limit
$N_{\sigma}$ corresponds to $C_{\sigma}=\sum_{n}\int_{\rm
BZ}d^{2}\mathbf{k}\theta(-\epsilon_{\mathbf{k}n\sigma})\Omega_{n\sigma}(\mathbf{k})/(2\pi)$
with
$\Omega(\mathbf{k})=i\innerproduct{\partial_{k_{x}}u_{\mathbf{k}n}}{\partial_{k_{y}}u_{\mathbf{k}n}}+h.c.$
the Berry curvature in the non-interacting limit and $u_{\mathbf{k}n}$ being
the eigenstates. In the interacting case $N_{\sigma}$ can be different from
the spin Hall conductance [25], as recently outlined in Refs. [17, 27, 16,
28]. In the Mott insulating phase $N_{\sigma}$ is determined by the Green
function zeros. Interestingly, applying the non-interacting electron formulae
to the mean field spinon Hamiltonian gives the same result as Eq. 5,
indicating a connection between the topology of zeros and spinons that will be
explored more fully below.

Fig. 2b) presents the evolution of the spectrum of zeros and spinons across
the topological transition at $\lambda=0$. At $\lambda/t=-0.5$ (left panel) we
observe that both the locus of Green function zeros and the spinon dispersion
are gapped. At the transition point $\lambda=0$ we observe gapless (Dirac)
behavior near $\mathbf{K}$ and $\mathbf{K}^{\prime}$ in the spinon spectrum
and in the locus of Green function zeros. The gapless behavior is protected by
the combination of $C_{3z}$ and $C_{2z}\mathcal{T}$ symmetries. For
$\lambda>0$ the gaps in the spinon dispersion and zero locus re-open and the
sign of the Chern number changes. More explicitly, we found that for each of
the two spin blocks both spinons and zeros are described by an effective
Hamiltonian which around $\mathbf{K},\mathbf{K}^{\prime}$ ($\tau=\pm)$ reads:

$\mathcal{H}_{\tau}(\mathbf{k})=v(\tau
k_{x}\sigma^{x}+k_{y}\sigma^{y})+\tau\Delta\sigma^{z}+O(k^{2}),$ (6)

where both $v$ and $\Delta$ depend on the interaction strength $U$ and on
whether the spinon or Green function zeros are considered. At the topological
transition $\Delta\rightarrow 0$. This behavior is further evidence of a deep
connection between Green function zeros and the spinon spectrum in topological
Mott insulators.

Edge zeros and spinons – We now consider edge states by generalizing the
slave-rotor method to allow for site/bond dependent mean field parameters as
detailed in the Supplementary Material [43]. In Fig. 3 we compare results
obtained for cylindrical geometries periodic in the direction parallel to one
of the lattice vectors, $\mathbf{a}_{1}$, and either open or periodic in the
other direction (see insets to Fig 3). The periodic case (along
${\mathbf{a}_{2}}$) shows gapped spinons as well as gapped zeros. The results
for open boundary condition reveal instead gapless modes in both the spinon
dispersion and the zeros. As shown in the Supplementary Material (Fig. S3) the
gapless modes are localized on the boundaries, mirroring the behavior of the
non-interacting case but with one important difference: the direction of the
edge zeros is reversed compared to the edge spinons and to the non-interacting
case.

Figure 3: Slave-rotor results for the Kane-Mele-Hubbard model with
$\lambda/t=0.2$ and $U/t=5.5$ on a cylindrical geometry periodic along the
crystal $\mathbf{a}_{1}$ direction and 12 sites along the orthogonal direction
($\mathbf{a}_{2}$) with either open (OBC, left panel) or periodic (PBC, right
panel) boundary conditions. $\overline{\mathbf{K}}$,
$\overline{\mathbf{K}}^{\prime}$, $\overline{\mathbf{M}}^{\prime}$ are the
projections of the high symmetry Brillouin zone points
$\mathbf{K},\mathbf{K}^{\prime},\mathbf{M}^{\prime}$ on the $\mathbf{a}_{1}$
direction. The color corresponds to the determinant of the Green function and
solid black lines indicate the spinon dispersion. In the case with OBC, we
find gapless modes with linear dispersion. For each spin block there is one
spinon/zero eigenvalue localized at either edge of the cylinder (see Fig. S3).

Annihilation of edge zeros and poles –  Next, we study an interface of a
strongly interacting slab which consequently hosts edge zeros and a weakly
correlated slab with edge poles. We consider a cylindrical geometry with
periodic boundary conditions along ${\mathbf{a}_{1}}$ and open boundary
conditions along ${\mathbf{a}_{2}}$ such that we have two (identical)
interfaces (see Fig. 4 a,b). Because of the simultaneous presence of both
poles and zeros, the analysis is best made by directly inspecting the
eigenvalues of the Green function individually, rather than looking at its
determinant. In Fig. 4e-h we show the minimal and maximal eigenvalues of the
Green’s function, thus resolving zeros and poles respectively. When the two
halves of the torus are disconnected (left column), gapless edge poles and
edge zeros are evident. Upon putting them in contact with each other (right
column), the gapless modes disappear, consistent with a mutual hybridization
of the two. This behavior is mirrored for the spinons: in the uncoupled case
both sides of the slab host gapless edge spinons (see Fig. 4c), although with
different renormalization. When coupled, the gapless edge spinons disappear
(Fig. 4d).

It is possible to understand the annihilation using two different pictures –
either looking at the spinons or the zeros. Considering the spinons, the TMI
has a spinon edge mode which at the interface has the opposite direction
compared to the edge state of the TI. Consequently, the edge modes cancel each
other (as sketched in Fig. 1). From the perspective of the zeros the situation
is slightly different: the zeros have the same direction at the interface as
the TI edge states (see Fig. S3 and also Ref.[15]) and instead of compensating
the edge state, they annihilate it.

Figure 4: a) and b): Sketch of the (un-)coupled Kane-Mele-Hubbard slab
geometry. One side, shown in orange, is in the Mott phase ($U/t=5.5$), the
other (green) is weakly correlated ($U/t=2.5$). The spin-orbit coupling is set
to $\lambda=0.2t$. Along ${\mathbf{a}_{1}}$ we consider periodic boundary
conditions, i.e. $k_{1}$ is a good quantum number. Instead, along
${\mathbf{a}_{2}}$ we consider a finite number of sites, i.e. 24 sites in
total. c) and d): spinon eigenvalues for the uncoupled and coupled case. The
color indicates on which side of the interface the eigenvalues are localized.
e) and f): maximal absolute value eigenvalue of the slave-rotor Green function
for the two sketched cases. g) and h): minimal absolute value eigenvalue. In
the disconnected case, c), e) and g), gapless edge spinons, edge poles and
edge zeros appear. Upon coupling the two slabs, i.e. going from c), e) and g)
to d), f) and h), the edge modes annihilate.

Spin and charge edge modes separation –  To provide further understanding of
the transport properties of the interface between a 2D TI and a 2D Mott TI we
consider the low energy theory of the edge modes of the two systems implied by
the slave rotor results. The edge modes of the conventional TI are chiral
fermions with chirality defined by the spin and described by the fermionic
Hamiltonian $H_{\rm TI}=\sum_{\sigma}\int
dxf^{\dagger}_{L\sigma}\left(-iv_{F}\sigma\partial_{x}\right)f_{L\sigma}$
where $\sigma=\pm$ for spin $\uparrow/\downarrow$. However, the edge modes of
the topological Mott Insulator are chiral spinons described by the Hamiltonian
$H_{\rm TMI}=\sum_{\sigma}\int
dx\psi^{\dagger}_{R\sigma}\left(+iv^{*}_{F}\sigma\partial_{x}\right)\psi_{R\sigma}$
with Fermi velocity $v^{*}_{F}$ determined by the exchange constants in the
Mott TI.

To understand the coupling at the interface we fractionalize the fermionic
edge modes of the conventional TI, writing separately the charge and spinon
parts and we then consider that at the interface the spinon modes are coupled.
The time-reversal symmetry and the spin $S^{z}$ conservation constrain the
coupling between the two different edge spinons to $\Delta H=m\int
dx\,\sum_{\sigma}\psi^{\dagger}_{L\sigma}\psi_{R\sigma}+h.c.$ so that the
resulting 1D spinon Hamiltonian,
$H_{\sigma}=\sigma(v_{F}-v^{*}_{F})k/2+\sigma(v_{F}+v^{*}_{F})k\mu^{z}/2+m\mu^{x}$
with $\mu$ Pauli matrices in the $L,R$ degrees of freedom, predicts a gapped
edge spectrum in agreement with the results obtained in our explicit slave
rotor calculation. The charge modes of the topological insulator, which remain
ungapped, provide charge transport confined at the interface between the TI
and the TMI, leading to the edge state pattern shown in Fig. 1, where a
conventional electronic edge state flows along the free sides of the
conventional TI and a conventional spinon edge state flows along the free
edges of the Mott TI, but at the interface the spinon modes are gapped,
leaving a gapless charge mode to carry the current through the interface. Note
that the gap in the spin spectrum means that a physical electron cannot
propagate along the Mott TI/TI interface, explaining why the wave packet
evolution of a physical electron presented in Ref. [15] decays as it
propagates into the interface.

The picture just presented is based on a chiral edge mode analysis. In the
physical case of unbroken time reversal symmetry with two counterpropagating
modes in each system, additional possibilities may occur: in particular
residual interactions may lead to gapping of the interface charge modes and
charge backscattering at the corner where the free edge of the conventional TI
meets the TI/Mott TI interface. Study of this situation is left to a future
paper.

The charge-spin separation taking place at the boundary can be also inferred
employing the Ioffe-Larkin composition rules. To this aim we consider the
interface of two semi-infinite planes, where the TI edge electrons interact
with the edge spinon of the TMI. The longitudinal conductivity along the
interface can be readily obtained as
$\sigma=\sigma_{f}+1/(\sigma^{-1}_{\psi}+\sigma^{-1}_{X})$ [48, 49] where
$\sigma_{f}$ is the contribution of the TI edge state, and the contribution
from the TMI is comprised by the spinon part $\sigma_{\psi}$ and the rotor
contribution $\sigma_{X}$. In the topological Mott insulator the rotor degree
of freedom $X$ is gapped and, as a result, $\sigma_{X}=0$. Thus, the
conductivity reduces to $\sigma=\sigma_{f}$ where each spin species
contributes with the quantized value $\sigma=e^{2}/h$. In other words, at the
interface only the edge modes of the TI couples to the electromagnetic field
and contribute to the charge transport. Finally, assuming to reduce the value
of the interaction below the Mott transition, when the slave rotor is
condensed $\sigma_{X}=\infty$, the conductivity becomes
$\sigma=\sigma_{f}+\sigma_{\psi}$. Both edge modes have finite charge, spin
and charge are bound in the electronic quasiparticle and the spin-resolved
conductivity vanishes.

Conclusions – The results of this paper indicate a deep connection between the
dispersion of Green function zeros (a mathematical concept introduced in the
context of generalizations of the Luttinger theorem) and the spinon spectrum
of a topological Mott insulator. Using slave rotor calculations to study the
Kane-Mele Hubbard model we found that in the bulk both the Green function
zeros and the spinon spectrum exhibit gaps, that close at the topological
transition where the spin-dependent topological index $N_{\sigma}$ changes
sign. In the topological phases the Green function zeros and the spinon
spectrum exhibit gapless boundary modes. At an interface between a topological
Mott insulator and conventional topological insulator the coupling between the
Green function zeros/gapless spinons and the fermionic edge state of the
conventional topological insulator leads to a novel spin-charge separated
state in which the spin degrees of freedom are gapped and charge transport
along the boundary is provided by charge $e$ spin zero holon states.

Our results raise many questions and suggest many directions for future
research. We investigated the large-$U$ regime of a Kane-Mele-Hubbard model
focusing on the paramagnetic Mott insulator without discussing the competition
with magnetic phases [50, 51, 37] or charge order/sublattice potential;
extending our results to these cases is of great interest. Further, we
restricted ourselves to the particle-hole symmetric case. Breaking particle-
hole symmetry by varying an external potential leads to interesting interface
structures now under investigation. The theory of backscattering at the
junction where the non-topological and topological edge modes meet is of both
theoretical and experimental interest, as is the extension of our
considerations to chiral Mott states and fractional Chern states. Going beyond
the approximate slave-rotor method used here, and studying a wider range of
models is also of great interest; and many-body physics of the coupled edge
states may yield further novel interface phases.

Note: During the completion of this manuscript we became aware of a related
work [52]; where there is overlap our results agree.

Acknowledgments – We acknowledge useful discussions with R. Raimondi, J. H.
Pixley, J. Cano, L. Crippa, E. König, S. Bollmann, M. Fabrizio, S. Sachdev and
B. Trauzettel. N.W. and G.S. are supported by the SFB 1170 Tocotronics, funded
by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
Project-ID 258499086. G.S. acknowledges financial support from the DFG through
the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in
Quantum Matter ct.qmat (EXC 2147, project-id 390858490). AJM acknowledges
support from the Columbia Materials Science and Engineering Research Center on
Precision Assembled Quantum Materials through NSF grant DMR-2011738. The
Flatiron Institute is a division of the Simons Foundation. The numerical
calculations have been carried out using the Julia language [53]. N.W and G.S
gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-
centre.eu) for funding this project by providing computing time on the GCS
Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre (www.lrz.de). This
research was supported in part by grants NSF PHY-1748958 and PHY-2309135 to
the Kavli Institute for Theoretical Physics (KITP).

## References

  * Oshikawa [2000] M. Oshikawa, Topological Approach to Luttinger’s Theorem and the Fermi Surface of a Kondo Lattice, Phys. Rev. Lett. 84, 3370 (2000), publisher: American Physical Society.
  * Senthil  _et al._ [2003] T. Senthil, S. Sachdev, and M. Vojta, Fractionalized Fermi Liquids, Physical Review Letters 90, 216403 (2003).
  * Paramekanti and Vishwanath [2004] A. Paramekanti and A. Vishwanath, Extending Luttinger’s theorem to Z_2 fractionalized phases of matter, Phys. Rev. B 70, 245118 (2004), publisher: American Physical Society.
  * Rosch [2007] A. Rosch, Breakdown of Luttinger’s theorem in two-orbital Mott insulators, The European Physical Journal B 59, 495 (2007).
  * Sakai  _et al._ [2009] S. Sakai, Y. Motome, and M. Imada, Evolution of Electronic Structure of Doped Mott Insulators: Reconstruction of Poles and Zeros of Green’s Function, Physical Review Letters 102, 056404 (2009).
  * Sakai  _et al._ [2010] S. Sakai, Y. Motome, and M. Imada, Doped high- T c cuprate superconductors elucidated in the light of zeros and poles of the electronic Green’s function, Physical Review B 82, 134505 (2010).
  * Volovik [2012] G. E. Volovik, Topology of quantum vacuum (2012), arXiv:1111.4627 [cond-mat, physics:gr-qc, physics:hep-ph].
  * Seki and Yunoki [2017] K. Seki and S. Yunoki, Topological interpretation of the Luttinger theorem, Phys. Rev. B 96, 085124 (2017), publisher: American Physical Society.
  * Heath and Bedell [2020] J. T. Heath and K. S. Bedell, Necessary and sufficient conditions for the validity of Luttinger’s theorem, New Journal of Physics 22, 063011 (2020), publisher: IOP Publishing.
  * Wen [2021] X.-G. Wen, Low-energy effective field theories of fermion liquids and the mixed $u(1)\times{}\mathbb{R}^{d}$ anomaly, Phys. Rev. B 103, 165126 (2021), publisher: American Physical Society.
  * Else  _et al._ [2021] D. V. Else, R. Thorngren, and T. Senthil, Non-Fermi Liquids as Ersatz Fermi Liquids: General Constraints on Compressible Metals, Phys. Rev. X 11, 021005 (2021), publisher: American Physical Society.
  * Skolimowski and Fabrizio [2022] J. Skolimowski and M. Fabrizio, Luttinger’s theorem in the presence of Luttinger surfaces, Physical Review B 106, 10.1103/physrevb.106.045109 (2022), publisher: American Physical Society (APS).
  * Fabrizio [2020] M. Fabrizio, Landau-Fermi liquids without quasiparticles, Phys. Rev. B 102, 155122 (2020), publisher: American Physical Society.
  * Fabrizio [2022] M. Fabrizio, Emergent quasiparticles at Luttinger surfaces, Nature communications 13, 1561 (2022), publisher: Nature Publishing Group UK London.
  * Wagner  _et al._ [2023] N. Wagner, L. Crippa, A. Amaricci, P. Hansmann, M. Klett, E. J. König, T. Schäfer, D. Di Sante, J. Cano, A. J. Millis, A. Georges, and G. Sangiovanni, Mott insulators with boundary zeros, Nature Communications 14, 7531 (2023).
  * Setty  _et al._ [2023a] C. Setty, F. Xie, S. Sur, L. Chen, M. G. Vergniory, and Q. Si, Electronic properties, correlated topology and Green’s function zeros (2023a), arXiv:2309.14340 [cond-mat].
  * Blason and Fabrizio [2023] A. Blason and M. Fabrizio, Unified role of Green’s function poles and zeros in topological insulators (2023), _eprint: 2304.08180.
  * Setty  _et al._ [2023b] C. Setty, S. Sur, L. Chen, F. Xie, H. Hu, S. Paschen, J. Cano, and Q. Si, Symmetry constraints and spectral crossing in a Mott insulator with Green’s function zeros (2023b), arXiv:2301.13870 [cond-mat].
  * Setty  _et al._ [2023c] C. Setty, F. Xie, S. Sur, L. Chen, S. Paschen, M. G. Vergniory, J. Cano, and Q. Si, Topological Diagnosis of Strongly Correlated Electron Systems (2023c), arXiv:2311.12031 [cond-mat].
  * Ishikawa and Matsuyama [1986] K. Ishikawa and T. Matsuyama, Magnetic field induced multi-component QED 3 and quantum Hall effect, Zeitschrift für Physik C Particles and Fields 33, 41 (1986), publisher: Springer.
  * Ishikawa and Matsuyama [1987] K. Ishikawa and T. Matsuyama, A microscopic theory of the quantum Hall effect, Nuclear Physics B 280, 523 (1987).
  * Gurarie [2011] V. Gurarie, Single-particle Green’s functions and interacting topological insulators, Physical Review B 83, 085426 (2011).
  * Essin and Gurarie [2011] A. M. Essin and V. Gurarie, Bulk-boundary correspondence of topological insulators from their respective Green’s functions, Physical Review B 84, 125132 (2011).
  * Wang and Zhang [2012] Z. Wang and S.-C. Zhang, Simplified Topological Invariants for Interacting Insulators, Physical Review X 2, 031008 (2012).
  * Slagle  _et al._ [2015] K. Slagle, Y.-Z. You, and C. Xu, Exotic quantum phase transitions of strongly interacting topological insulators, Phys. Rev. B 91, 115121 (2015).
  * Manmana  _et al._ [2012] S. R. Manmana, A. M. Essin, R. M. Noack, and V. Gurarie, Topological invariants and interacting one-dimensional fermionic systems, Physical Review B 86, 205119 (2012).
  * Gavensky  _et al._ [2023] L. P. Gavensky, S. Sachdev, and N. Goldman, Connecting the many-body Chern number to Luttinger’s theorem through St\v{r}eda’s formula (2023), arXiv:2309.02483 [cond-mat].
  * Zhao  _et al._ [2023] J. Zhao, P. Mai, B. Bradlyn, and P. Phillips, Failure of Topological Invariants in Strongly Correlated Matter, Physical Review Letters 131, 106601 (2023).
  * Aasen  _et al._ [2020] D. Aasen, R. S. Mong, B. M. Hunt, D. Mandrus, and J. Alicea, Electrical Probes of the Non-Abelian Spin Liquid in Kitaev Materials, Physical Review X 10, 031014 (2020).
  * Raghu  _et al._ [2008] S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang, Topological Mott Insulators, Phys. Rev. Lett. 100, 156401 (2008), publisher: American Physical Society.
  * Maciejko and Fiete [2015] J. Maciejko and G. A. Fiete, Fractionalized topological insulators, Nature Physics 11, 385 (2015), publisher: Nature Publishing Group UK London.
  * 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 (2010), publisher: Springer Science and Business Media LLC.
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Z 2 Topological Order and the Quantum Spin Hall Effect, Physical Review Letters 95, 146802 (2005).
  * Zhao  _et al._ [2022] W. Zhao, K. Kang, L. Li, C. Tschirhart, E. Redekop, K. Watanabe, T. Taniguchi, A. Young, J. Shan, and K. F. Mak, Realization of the Haldane Chern insulator in a moiré lattice (2022).
  * Mak and Shan [2022] K. F. Mak and J. Shan, Semiconductor moiré materials, Nature Nanotechnology 17, 686 (2022), publisher: Nature Publishing Group UK London.
  * Haldane [1988] F. D. M. Haldane, Model for a Quantum Hall Effect without Landau Levels: Condensed-Matter Realization of the “Parity Anomaly”, Phys. Rev. Lett. 61, 2015 (1988), publisher: American Physical Society.
  * Rachel and Le Hur [2010] S. Rachel and K. Le Hur, Topological insulators and Mott physics from the Hubbard interaction, Phys. Rev. B 82, 075106 (2010), publisher: American Physical Society.
  * Rachel [2018] S. Rachel, Interacting topological insulators: a review, Reports on Progress in Physics 81, 116501 (2018), publisher: IOP Publishing.
  * Fernández López and Merino [2022] M. Fernández López and J. Merino, Bad topological semimetals in layered honeycomb compounds, Physical Review B 105, 115138 (2022).
  * Florens and Georges [2002] S. Florens and A. Georges, Quantum impurity solvers using a slave rotor representation, Phys. Rev. B 66, 165111 (2002), publisher: American Physical Society.
  * Florens  _et al._ [2003] S. Florens, P. San José, F. Guinea, and A. Georges, Coherence and Coulomb blockade in single-electron devices: A unified treatment of interaction effects, Phys. Rev. B 68, 245311 (2003), publisher: American Physical Society.
  * Florens and Georges [2004] S. Florens and A. Georges, Slave-rotor mean-field theories of strongly correlated systems and the Mott transition in finite dimensions, Phys. Rev. B 70, 035114 (2004), publisher: American Physical Society.
  * [43] See supplementary materials at URL.
  * Lee and Lee [2005] S.-S. Lee and P. A. Lee, U(1) gauge theory of the hubbard model: Spin liquid states and possible application to $\kappa\mathrm{\text{$-$}}(\mathrm{BEDT}\mathrm{\text{$-$}}\mathrm{TTF}{)}_{2}{\mathrm{cu}}_{2}(\mathrm{CN}{)}_{3}$, Phys. Rev. Lett. 95, 036403 (2005).
  * He and Lee [2022] W.-Y. He and P. A. Lee, Magnetic impurity as a local probe of the U(1) quantum spin liquid with spinon Fermi surface, Phys. Rev. B 105, 195156 (2022), publisher: American Physical Society.
  * He and Lee [2023] W.-Y. He and P. A. Lee, Electronic density of states of a U(1) quantum spin liquid with spinon Fermi surface. II. Zeeman magnetic field effects, Phys. Rev. B 107, 195156 (2023), publisher: American Physical Society.
  * Haldane [2004] F. D. M. Haldane, Berry Curvature on the Fermi Surface: Anomalous Hall Effect as a Topological Fermi-Liquid Property, Phys. Rev. Lett. 93, 206602 (2004), publisher: American Physical Society.
  * Ioffe and Larkin [1989] L. B. Ioffe and A. I. Larkin, Gapless fermions and gauge fields in dielectrics, Physical Review B 39, 8988 (1989).
  * Aldape  _et al._ [2022] E. E. Aldape, T. Cookmeyer, A. A. Patel, and E. Altman, Solvable theory of a strange metal at the breakdown of a heavy Fermi liquid, Physical Review B 105, 235111 (2022).
  * Ran  _et al._ [2008] Y. Ran, A. Vishwanath, and D.-H. Lee, A direct transition between a Neel ordered Mott insulator and a $d_{x^{2}-y^{2}}$ superconductor on the square lattice (2008), arXiv:0806.2321 [cond-mat].
  * Hermele  _et al._ [2004] M. Hermele, T. Senthil, M. P. A. Fisher, P. A. Lee, N. Nagaosa, and X.-G. Wen, Stability of U ( 1 ) spin liquids in two dimensions, Physical Review B 70, 214437 (2004).
  * [52] S. Bollmann, C. Setty, U. Seifert, and E. König, Topological green’s function zeros in an exactly solved model and beyond, in preparation.
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Review 59, 65 (2017).
  * Moshe and Zinn-Justin [2003] M. Moshe and J. Zinn-Justin, Quantum field theory in the large N limit: a review, Physics Reports 385, 69 (2003), publisher: Elsevier BV.
  * Coleman [1984] P. Coleman, New approach to the mixed-valence problem, Phys. Rev. B 29, 3035 (1984), publisher: American Physical Society.
  * Kotliar and Ruckenstein [1986] G. Kotliar and A. E. Ruckenstein, New Functional Integral Approach to Strongly Correlated Fermi Systems: The Gutzwiller Approximation as a Saddle Point, Physical Review Letters 57, 1362 (1986).
  * dem [2005] Orbital-selective Mott transition in multiband systems: Slave-spin representation and dynamical mean-field theory, Phys. Rev. B 72, 205124 (2005), publisher: American Physical Society.
  * Lechermann  _et al._ [2007] F. Lechermann, A. Georges, G. Kotliar, and O. Parcollet, Rotationally invariant slave-boson formalism and momentum dependence of the quasiparticle weight, Physical Review B 76, 10.1103/physrevb.76.155102 (2007), publisher: American Physical Society (APS).
  * Guerci and Fabrizio [2017] D. Guerci and M. Fabrizio, Unbinding slave spins in the Anderson impurity model, Phys. Rev. B 96, 201106 (2017), publisher: American Physical Society.
  * Raimondi and Castellani [1993] R. Raimondi and C. Castellani, Lower and upper Hubbard bands: A slave-boson treatment, Phys. Rev. B 48, 11453 (1993), publisher: American Physical Society.
  * Rüegg  _et al._ [2010] A. Rüegg, S. D. Huber, and M. Sigrist, \mathsfZ_2-slave-spin theory for strongly correlated fermions, Phys. Rev. B 81, 155118 (2010), publisher: American Physical Society.
  * Riegler  _et al._ [2020] D. Riegler, M. Klett, T. Neupert, R. Thomale, and P. Wölfle, Slave-boson analysis of the two-dimensional Hubbard model, Phys. Rev. B 101, 235137 (2020), publisher: American Physical Society.
  * Seufert  _et al._ [2021] J. Seufert, D. Riegler, M. Klett, R. Thomale, and P. Wölfle, Breakdown of charge homogeneity in the two-dimensional Hubbard model: Slave-boson study of magnetic order, Phys. Rev. B 103, 165117 (2021), publisher: American Physical Society.

Supplementary material for
“”

Niklas Wagner1, Daniele Guerci2, Andrew J. Millis3,2 and Giorgio Sangiovanni1

1Institut für Theoretische Physik und Astrophysik and Würzburg-Dresden Cluster
of Excellence ct.qmat, Universität Würzburg, 97074 Würzburg, Germany, Germany

2 Center for Computational Quantum Physics, Flatiron Institute, 162 5th
Avenue, NY 10010, USA
3 Department of Physics, Columbia University, New York, NY 10027, USA

These supplementary materials contain the details of analytic calculations as
well as additional numerical details supporting the results presented in the
main text.

##  Appendix A Conventions

We choose a system of coordinates in which the scalar product of the
$\pm{\bm{K}}$ vectors with the nearest neighbor vectors ${\bm{u}}_{j=1,2,3}$
and ${\bm{\gamma}}_{j=1,2,3}$ are given by

$\begin{array}[]{c||c|c|c|c|c|c}{\bm{r}}&{\bm{\gamma}}_{1}&{\bm{\gamma}}_{2}&{\bm{\gamma}}_{3}&{\bm{u}}_{1}&{\bm{u}}_{2}&{\bm{u}}_{3}\\\
\hline\cr\hline\cr
e^{i\mathbf{K}\cdot{\bm{r}}}&\omega^{*}&\omega^{*}&\omega^{*}&1&\omega&\omega^{*}\\\
\hline\cr
e^{i\mathbf{K}^{\prime}\cdot{\bm{r}}}&\omega&\omega&\omega&1&\omega^{*}&\omega\end{array}\,\,,\quad\omega=e^{2i\pi/3}.$
(1)

This can be realized, for instance, with
${\bm{u}}_{j}=C^{j-1}_{3z}{\bm{u}}_{1}$ and ${\bm{u}}_{1}=a(1,0)/\sqrt{3}$,
${\bm{a}}_{1}=a(\sqrt{3},1)/2$, ${\bm{a}}_{2}=a(-\sqrt{3},1)/2$,
$\bm{K}=-\bm{K}^{\prime}=(4\pi/3a)(0,1)$ and
${\bm{\gamma}}_{1}={\bm{a}}_{1}+{\bm{a}}_{2}$,
${\bm{\gamma}}_{j}=C^{j-1}_{3z}{\bm{\gamma}}_{1}$. We also introduce the
Fourier transform of the hopping process:

$\begin{split}&\epsilon_{\mathbf{k}}=2\lambda\sum^{3}_{j=1}\cos{\bm{\gamma}}_{j}\cdot\mathbf{k},\\\
&\epsilon_{\mathbf{k}A/B\sigma}=2\lambda\sum^{3}_{j=1}\cos\left({\bm{\gamma}}_{j}\cdot\mathbf{k}\pm\pi\sigma/2\right),\\\
&V_{\mathbf{k}}=\sum^{3}_{j=1}e^{i\mathbf{k}\cdot{\bm{u}_{j}}}.\end{split}$
(2)

##  Appendix B Details on the slave rotor mapping and on the large-$M$ theory

In this section we detail the slave rotor mapping, providing the Hamiltonian
expressed in terms of the auxiliary representation, the saddle-point equations
and the self-consistent expression of the interacting Green function.
Employing the slave rotor mapping introduced in the main text Eq. (3) the
interacting part of the Hamiltonian takes the form:

$H_{\rm int}=\frac{U}{2}\sum_{s,\mathbf{r}\in s}L^{s\,2}_{\mathbf{r}}.$ (3)

Hoppings are dressed by a dynamical Peierls phase:

$H_{t}=-t\sum_{\langle\mathbf{r}\in A,\mathbf{r}^{\prime}\in
B\rangle}e^{-i\left(\theta^{A}_{\mathbf{r}}-\theta^{B}_{\mathbf{r}^{\prime}}\right)}\psi^{A\dagger}_{\mathbf{r}}\psi^{B}_{\mathbf{r}^{\prime}}+h.c.,$
(4)

and

$H_{\rm
SOC}=\lambda\sum_{s=A,B}\sum_{\langle\mathbf{r},\mathbf{r}^{\prime}\rangle\in
s}e^{-i(\theta^{s}_{\mathbf{r}}-\theta^{s}_{\mathbf{r}^{\prime}})}\psi^{s\dagger}_{\mathbf{r}}e^{i\frac{\pi}{2}\nu_{\mathbf{r},\mathbf{r}^{\prime}}\sigma^{z}}\psi^{s}_{\mathbf{r}^{\prime}}.$
(5)

The mapping is supplemented by the constraint:

$L^{s}_{\mathbf{r}}-\sum_{\sigma}\psi^{s\dagger}_{\mathbf{r}\sigma}\psi^{s}_{\mathbf{r}\sigma}+1=0.$
(6)

To make progress we turn to the action formalism:

$\begin{split}\mathcal{S}=&\int^{\beta}_{0}d\tau\sum_{s=a,b}\sum_{\mathbf{r}\in
s}\left[-iL^{s}_{\mathbf{r}}\partial_{\tau}\theta^{s}_{\mathbf{r}}+\psi^{s\dagger}_{\mathbf{r}}\left(\partial_{\tau}-h^{s}_{\mathbf{r}}\right)\psi^{s}_{\mathbf{r}}\right.\\\
&\left.+h^{s}_{\mathbf{r}}\left(L^{s}_{\mathbf{r}}+1\right)\right]+\frac{U}{2}\sum_{s,r\in
s}L^{s\,2}_{r}+H_{0}[\theta,\psi],\end{split}$ (7)

where $\tau\in[0,\beta]$ with $\beta=1/T$, $h^{s}_{\mathbf{r}}(\tau)$ is the
field enforcing the constraint in Eq. (6),
$H_{0}[\theta,\psi]=H_{\Delta}+H_{t}+H_{\rm SOC}$ is expressed in terms of the
spinon $\psi$ and the phase variable $e^{i\theta}$ see Eqs. (3),(4) and (5),
respectively. Integrating out $L^{s}_{\mathbf{r}}$ we find:

$\begin{split}\mathcal{L}=&\sum_{s,\mathbf{r}\in
s}\psi^{s\dagger}_{\mathbf{r}}\left(\partial_{\tau}-h^{s}_{\mathbf{r}}\right)\psi^{s}_{\mathbf{r}}+\frac{(\partial_{\tau}\theta^{s}_{\mathbf{r}}+ih^{s}_{\mathbf{r}})^{2}}{2U}+h^{s}_{\mathbf{r}}\\\
&+H_{0}[\theta,\psi].\end{split}$ (8)

The simplest yet non-trivial treatment of the fluctuations of the theory (8)
was developed in Ref. [40, 42] and consists of introducing the complex boson
$X^{s}_{\mathbf{r}}$:

$X^{s}_{\mathbf{r}}=e^{i\theta^{s}_{\mathbf{r}}},$ (9)

with the local constraint $|X^{s}_{\mathbf{r}}(\tau)|^{2}=1$ which is enforced
by the additional Lagrangian multiplier $\lambda^{s}_{\mathbf{r}}$.
Introducing $X$ in the Lagrangian (8) we find:

$\begin{split}\mathcal{L}=&\sum_{s,\mathbf{r}\in
s}\psi^{s\dagger}_{\mathbf{r}}\left(\partial_{\tau}-h^{s}_{\mathbf{r}}\right)\psi^{s}_{\mathbf{r}}+\frac{|\partial_{\tau}X^{s}_{\mathbf{r}}|^{2}-h^{s2}_{\mathbf{r}}}{2U}+h^{s}_{\mathbf{r}}\\\
&+\sum_{s,\mathbf{r}\in
s}\rho^{s}_{\mathbf{r}}\left(|X^{s}_{\mathbf{r}}|^{2}-1\right)+\frac{h^{s}_{\mathbf{r}}}{2U}\left(X^{s*}_{\mathbf{r}}\partial_{\tau}X^{s}_{\mathbf{r}}-h.c.\right)\\\
&+H_{0}[X,\psi].\end{split}$ (10)

Then, the mean field theory can be rigorously derived as the saddle-point of
the non-linear sigma model in Eq. (10) in the large-$M$ limit [54]. The limit
consists of replacing O(2) with O(2$M$) or, in other words, extending
$X^{s}_{\mathbf{r}}$ to an $M$ components field $X^{s,\alpha}_{\mathbf{r}}$
with $\alpha=1,\cdots,M$ and
$\sum^{M}_{\alpha=1}|X^{s,\alpha}_{\mathbf{r}}(\tau)|^{2}=1$. We also
introduce $\alpha=1,\cdots,M$ Fermi fields $\psi^{s,\alpha}_{\mathbf{r}}$ for
each spin and sublattice degree of freedom. Performing the Hubbard-
Stratonovich transformation which introduces the bond variables $Q^{X}$ and
$Q^{\psi}$ we have:

$\begin{split}&\mathcal{L}_{M}=\sum_{s,\mathbf{r}\in
s}\sum^{M}_{\alpha=1}\psi^{s,\alpha\dagger}_{\mathbf{r}}\left(\partial_{\tau}-h^{s}_{\mathbf{r}}\right)\psi^{s,\alpha}_{\mathbf{r}}-\frac{h^{s2}_{\mathbf{r}}}{2U}+h^{s}_{\mathbf{r}}+\sum_{s,\mathbf{r}\in
s}\rho^{s}_{\mathbf{r}}\left(\sum^{M}_{\alpha=1}|X^{s,\alpha}_{\mathbf{r}}|^{2}-M\right)\\\
&+\sum_{s,\mathbf{r}\in
s}\sum^{M}_{\alpha=1}\frac{|\partial_{\tau}X^{s,\alpha}_{\mathbf{r}}|^{2}}{2U}+\frac{h^{s}_{\mathbf{r}}}{2U}\left(X^{s,\alpha*}_{\mathbf{r}}\partial_{\tau}X^{s,\alpha}_{\mathbf{r}}-h.c.\right)-\lambda
M\sum_{s,\langle\mathbf{r},\mathbf{r}^{\prime}\rangle\in
s}Q^{X}_{ss}(\mathbf{r},\mathbf{r}^{\prime})Q^{\psi}_{ss}(\mathbf{r},\mathbf{r}^{\prime})\\\
&+\sum_{s,\langle\mathbf{r},\mathbf{r}^{\prime}\rangle\in
s}\sum^{M}_{\alpha=1}\left[Q^{\psi}_{ss}(\mathbf{r},\mathbf{r}^{\prime})\psi^{s,\alpha\dagger}_{\mathbf{r}}e^{i\frac{\pi}{2}\nu_{\mathbf{r},\mathbf{r}^{\prime}}\sigma^{z}}\psi^{s,\alpha}_{\mathbf{r}^{\prime}}+Q^{X}_{ss}(\mathbf{r},\mathbf{r}^{\prime})X^{s,\alpha*}_{\mathbf{r}}X^{s,\alpha}_{\mathbf{r}^{\prime}}\right]+tM\sum_{\langle\mathbf{r}\in
a,\mathbf{r}^{\prime}\in
b\rangle}Q^{X}_{ab}(\mathbf{r},\mathbf{r}^{\prime})Q^{\psi}_{ab}(\mathbf{r},\mathbf{r}^{\prime})+h.c.\\\
&-t\sum_{\langle\mathbf{r}\in a,\mathbf{r}^{\prime}\in
b\rangle}\sum^{M}_{\alpha=1}\left[Q^{\psi}_{ab}(\mathbf{r},\mathbf{r}^{\prime})\psi^{a,\alpha\dagger}_{\mathbf{r}}\psi^{b,\alpha}_{\mathbf{r}^{\prime}}+Q^{X}_{ab}(\mathbf{r},\mathbf{r}^{\prime})X^{a,\alpha*}_{\mathbf{r}}X^{b,\alpha}_{\mathbf{r}^{\prime}}+h.c.\right].\end{split}$
(11)

In the limit $M\to\infty$ the solution is obtained by taking the saddle point
of the previous Lagrangian. Assuming that translational invariance and the
point group symmetries of the lattice model are preserved, the self-
consistency equations associated with $h^{s}_{\mathbf{r}}$ and
$\rho^{s}_{\mathbf{r}}$ read:

$\begin{split}&\frac{T}{N_{s}}\sum_{\mathbf{k}\neq
0,i\nu}G^{ss}_{X}(\mathbf{k},i\nu)+Z^{ss}=1,\\\
&\frac{T}{N_{s}}\sum_{\sigma,\mathbf{k}}\sum_{i\omega}G^{ss}_{\sigma,\psi}(\mathbf{k},i\omega)=1-\frac{2h^{s}}{U}-\frac{T}{UN_{s}}\sum_{\mathbf{k}\neq
0,i\nu}i\nu G^{ss}_{X}(\mathbf{k},i\nu)\left[e^{i\nu 0^{+}}+e^{-i\nu
0^{+}}\right],\end{split}$ (12)

with $\omega=(2n+1)\pi T$, $\nu=2n\pi T$, $N_{s}$ number of unit cells. For
all calculations we used a temperature $T/t=0.01$. $Z^{ss}\equiv Z^{s}$ is the
quasiparticle weight which is determined self-consistently. Note that we do
not consider a staggered potential or chemical potential shift and thus the
second constraint is automatically fulfilled for $h^{s}=0$. In addition, we
also have the saddle-point equations for the bond variables $Q^{X}$:

$\begin{split}Q^{X}_{ss}&=\frac{T}{6\lambda
N}\sum_{\sigma,\mathbf{k}}\sum_{i\omega}\epsilon_{\mathbf{k}s\sigma}G^{ss}_{\sigma,\psi}(\mathbf{k},i\omega)e^{i\omega
0^{+}},\\\
Q^{X}_{ab}&=\frac{T}{3N}\sum_{\sigma,\mathbf{k}}\sum_{i\omega}V_{\mathbf{k}}G^{ba}_{\sigma,\psi}(\mathbf{k},i\omega)e^{i\omega
0^{+}},\end{split}$ (13)

and $Q^{\psi}$:

$\begin{split}Q^{\psi}_{ss}&=Z^{ss}+\frac{T}{6\lambda N}\sum_{\mathbf{k}\neq
0,i\nu}\epsilon_{\mathbf{k}}G^{ss}_{X}(\mathbf{k},i\nu),\\\
Q^{\psi}_{ab}&=Z^{ab}+\frac{T}{3N}\sum_{\mathbf{k}\neq
0,i\nu}V_{\mathbf{k}}G^{ba}_{X}(\mathbf{k},i\nu),\end{split}$ (14)

where the $\mathbf{k}=0$ condensate term is isolated from the finite
$\mathbf{k}$ contribution of the fluctuations. The quantities
$\epsilon_{\mathbf{k}}$, $\epsilon_{\mathbf{k}s\sigma}$ and $V_{\mathbf{k}}$
are the Fourier transform of the hopping process and are given in Eq.(2). In
Eq. (12) and (14) the condensate part [44] is given by
$Z^{ss^{\prime}}=\langle
X^{s\dagger}(\mathbf{k}=0,\tau)X^{s^{\prime}}(\mathbf{k}=0,\tau)\rangle/N=\sqrt{Z^{s}Z^{s^{\prime}}}$
with $s,s^{\prime}=A,B$. It is finite in the weakly correlated regime and
vanishes in the Mott insulator. The Green functions of the spinon and of the
rotor are given by:

$G^{-1}_{\sigma,\psi}(\mathbf{k},i\omega)=i\omega+h-\begin{bmatrix}Q^{\psi}_{aa}\epsilon_{\mathbf{k}a\sigma}&-tQ^{\psi}_{ab}V_{\mathbf{k}}\\\
-tQ^{\psi}_{ba}V^{*}_{\mathbf{k}}&Q^{\psi}_{bb}\epsilon_{\mathbf{k}b\sigma}\end{bmatrix},$
(15)

and

$G^{-1}_{X}(\mathbf{k},i\nu)=\frac{\nu^{2}}{U}-\frac{2i\nu
h}{U}+\rho+\begin{bmatrix}Q^{X}_{aa}\epsilon_{\mathbf{k}}&-tQ^{X}_{ab}V_{\mathbf{k}}\\\
-tQ^{X}_{ba}V^{*}_{\mathbf{k}}&Q^{X}_{bb}\epsilon_{\mathbf{k}}\end{bmatrix}.$
(16)

For the sake of completeness, we recall that in imaginary time and real space
the spinon and rotor Green functions are defined as
$G^{ss^{\prime}}_{\sigma,\psi}(\mathbf{r}-\mathbf{r}^{\prime},\tau)=-\langle
T_{\tau}\psi^{s}_{\mathbf{r}\sigma}(\tau)\psi^{s^{\prime}\dagger}_{\mathbf{r}^{\prime}\sigma}\rangle$
and $G^{ss^{\prime}}_{X}(\mathbf{r}-\mathbf{r}^{\prime},\tau)=\langle
T_{\tau}X^{s}_{\mathbf{r}}(\tau)X^{s^{\prime}\dagger}_{\mathbf{r}^{\prime}}\rangle$,
respectively. In the rotor Green’s function $G_{X}$ we only keep the
fluctuation part, the contribution of the condensed fraction has been singled
out in Eq. 12 and 14. Finally, within the large-$M$ limit it is
straightforward to find that the slave rotor expression of the interacting
single-particle Green function reads:

$\begin{split}G^{ss^{\prime}}_{\sigma}(\mathbf{k},i\omega)=&Z^{ss^{\prime}}G^{ss^{\prime}}_{\sigma,\psi}(\mathbf{k},i\omega)\\\
+\frac{T}{N}&\sum_{\mathbf{q},i\nu}G^{ss^{\prime}}_{\sigma,\psi}(\mathbf{k}-\mathbf{q},i\omega-i\nu)G^{ss^{\prime}}_{X}(\mathbf{q},i\nu),\end{split}$
(17)

with $i\nu$ bosonic Matsubara frequency $\nu=2nT$ and $\mathbf{q}$ momentum
transfer between the rotor and the spinon. We emphasize that the first term in
Eq. (17) is absent in the Mott insulating regime $Z^{ss^{\prime}}=0$ while the
second gives rise to the incoherent Hubbard bands and zeros discussed in the
main text.

###  B.1 Numerical solution of the self-consistency equations

In this section we provide additional results of the self-consistency
equations for periodic and cylindrical boundary conditions. The results for
the former are in agreement with earlier studies in Ref. [37, 39]. Fig. S1
shows the evolution of $Q^{X}$, $Q^{\psi}$, quasiparticle weight $Z$ and the
rotor gap as a function of $U$. The latter is related to the Lagrange
multiplier $\rho$ and given by $\Delta
X=\rho+\operatorname{min}(\epsilon_{X})$ where $\epsilon_{X}$ is the
dispersion of the rotors. We emphasize that the rotor gap grow linearly with
$(U-U_{c})$ after the Mott transition. Interestingly, the dispersion of the
spinon is finite and goes like $t/U$ even for $Z=0$ differently from the
results obtained at mean-field level with other slave-particle approaches [55,
56, 57, 58, 59]. This property results from the inclusion of Gaussian
fluctuations [60, 61, 62, 63] in the self-consistency equation which are taken
into account by the saddle-point equations obtained in the large-$M$ limit.

Turning to cylindrical boundary conditions, Figure S2 shows the site-dependent
values of the self-consistent parameters for a slab geometry with interface
between a TMI and a TI. As shown there the renormalization parameters vary
slightly within the bulk but are mostly unaffected by the interface – apart
from those directly at the interface. Note that the results shown in Fig.4 are
for a smaller slab, hence the spinon dispersion is more strongly influenced by
coupling the two sides. Finally, Fig. S3 shows the localization of edge
spinons and egde zeros, revealing the opposite direction of the two edge
modes.

Figure S1: Self-consistent slave-rotor parameters as a function of interaction
strength $U$. $\Delta X$ is the gap of the rotors, $Z$ is the quasiparticle
weight and $Q_{\psi}^{t}$ and $Q_{X}^{t}$ are the renormalization of the
$t$-hopping term for spinon and rotor respectively. The calculations are done
with $\lambda/t=0.5$.  Figure S2: Self-consistent slave-rotor parameters as a
function of site for a cylindrical geometry (see Fig. 4) with interface
between a TMI ($U/t=12$, orange shade) and a TI ( $U/t=4$, green shade).
$\Delta X$ is the gap of the rotors, $Z$ is the quasiparticle weight and
$Q_{\psi}^{t}$ and $Q_{X}^{t}$ are the renormalization of the $t$-hopping term
for spinon and rotor respectively. The calculations are done with
$\lambda/t=0.5$.

(a) Spinons

(b) Zeros

Figure S3:  Localization of edge spinons (a) and edge zeros (b) for the
results shown in Fig. 3 . The red-blue color scheme in (b) shows the
determinant of the slave-rotor Green function, zoomed-in on the edge state.
Black and green dots show the corresponding edge character.


value $a_{0}$. To simplify the dispersion laws, as before, one introduces the
plasma frequency of the moving particles, $\omega_{p}$, with
$\omega_{p}^{2}=q_{+}^{2}n_{0}/m\varepsilon_{0}$, and the normalized
wavevector $\vec{q}_{0}=\vec{k}a_{0}$ [24, 25, 27].

The computation of the term in the longitudinal dispersion law involving the
temperature derivatives at constant entropy is the same as in both previous
approaches. One uses results on the temperature derivative, as shown in Eq.
(32), and on the derivative at constant entropy of a function of $\Gamma$, as
demonstrated in Eq. (33).

To analyze terms with integrals of $g-1$ in both of the dispersion laws, in
addition to the previous identities, use the following strategy, which is
slightly different from the strategy used in the two previous approaches.
First, take number density derivatives outside of the integrals. Then, change
variables inside the integrals to $x=r/a_{0}$ and introduce the normalized
wavevector $\vec{q}_{0}=\vec{k}a_{0}$. In this case, integrals are rewritten
in terms of functions of $\Gamma(n_{0},T(n,s))$ and $|\vec{q}_{0}|$. Next, one
can use the following identity for an arbitrary function
$f(\Gamma(n_{0},T(n,s)),|\vec{q}_{0}|)$ on the derivative of number density at
constant entropy, where we use the chain rule, Eq. (32), and the definition of
$\Gamma$.

$\begin{gathered}n\frac{\partial f}{\partial
n}\Big{|}_{s}\big{(}\Gamma(n_{0},T(n,s)),|\vec{q}_{0}|\big{)}=-\frac{2}{3}\Gamma(n_{0},T(n,s))\frac{\partial
f}{\partial\Gamma}\Big{|}_{|\vec{q}_{0}|}f_{1}\big{(}\Gamma(n,T)\big{)}\end{gathered}$
(81)

The transverse dispersion relation in the normalized variables is the same as
in the Lagrangian approach, as given by Eq. (74), as it does not depend on
number density derivatives at constant entropy.

However, the longitudinal dispersion relation in the normalized variables is
now different. For each equilibrium value $\Gamma_{0}$, it is provided in Eq.
(82). The terms appearing in the dispersion relation are given by Eqs. (32),
(36), (37), (76), (77), (83), (84), (85).

$\begin{gathered}\left(\frac{\omega_{L}}{\omega_{p}}\right)^{2}(\Gamma_{0},|\vec{q}_{0}|)=1+\frac{|\vec{q}_{0}|^{2}}{\Gamma_{0}}f_{2}(\Gamma_{0})+j_{1}(\Gamma_{0},|\vec{q}_{0}|)+j_{2}(\Gamma_{0},|\vec{q}_{0}|)\\\
+b(\Gamma_{0},|\vec{q}_{0}|)+\ell_{1}(\Gamma_{0},|\vec{q}_{0}|),\end{gathered}$
(82) $\begin{gathered}j_{1}(\Gamma,|\vec{q}|)=-\frac{2\Gamma}{3}\frac{\partial
j}{\partial\Gamma}\Big{|}_{|\vec{q}|}(\Gamma,|\vec{q}|)\left(f_{1}(\Gamma)-\frac{2}{3}f_{1}^{2}(\Gamma)+\frac{\Gamma}{3}\frac{\mathrm{d}f_{1}}{\mathrm{d}\Gamma}(\Gamma)\Big{(}1-2f_{1}(\Gamma)\Big{)}\right),\end{gathered}$
(83)
$\begin{gathered}j_{2}(\Gamma,|\vec{q}|)=\frac{4\Gamma^{2}}{9}f_{1}^{2}(\Gamma)\frac{\partial^{2}j}{\partial\Gamma^{2}}\Big{|}_{|\vec{q}|}(\Gamma,|\vec{q}|),\end{gathered}$
(84)
$\begin{gathered}\ell_{1}(\Gamma,|\vec{q}|)=-\frac{2\Gamma}{3}f_{1}(\Gamma)\frac{\partial\ell}{\partial\Gamma}\Big{|}_{|\vec{q}|}(\Gamma,|\vec{q}|).\end{gathered}$
(85)

As discussed in both previous variational approaches in detail, due to
insufficient numerical data regarding how the pair distribution function
$g(\Gamma,|\vec{r}|/a)$ depends on $\Gamma$, we use a step function
approximation to simplify integrals and make calculations more practical. The
step function approximation for the pair distribution function is such that it
has a value of $0$ when $|\vec{r}|/a<(R/a)(\Gamma)$, and $1$ when
$|\vec{r}|/a>(R/a)(\Gamma)$.

The dependence of $R/a$ on $\Gamma$ is determined by computing the equilibrium
energy, as given in Eq. (20), for the given approximation. We then use Eq.
(30) to relate it to $u(\Gamma)$. The result of this computation is provided
in Eq. (40).

In this case, the transverse dispersion relation is provided in Eq. (78),
where $R/a$ is the correct value for a given $\Gamma_{0}$. This result agrees
with the QLCA result for the same pair distribution function approximation
[28]. The integrals found in the longitudinal dispersion relation, as given by
Eqs. (37), (76), (77), can be simplified as before, as described by Eqs. (41),
(79), where $R/a$ is the correct value for a given $\Gamma$.

Now, let us consider different fits for the nonlocal contribution to the
thermodynamic energy, $u(\Gamma)$, as defined in Eq. (30). First, let us
examine a simple fit $u(\Gamma)=-0.9\Gamma$, which is accurate for very strong
coupling where $\Gamma\gg 1$ [15]. As discussed previously, for such a fit,
$R/a=\sqrt{6/5}$, independent of the value of $\Gamma$. The results for the
transverse dispersion relation for different values of $\Gamma_{0}$ are
presented as before in Fig. 4(a).

As before, from Eqs. (41), (79), it follows that the functions $j,b,\ell$
appearing in the longitudinal dispersion relation, as given in Eq. (82), are
also independent of $\Gamma$. This significantly simplifies the dispresion
relation, and from Eqs. (83), (84), (85), we deduce that
$j_{1},j_{2},\ell_{1}$ are all zero. Additionally, as previously discussed,
for such a fit, we have $f_{1}(\Gamma)=1$ according to Eq. (32). With these
results, the longitudinal dispersion law, as given in Eq. (82) and
supplemented with Eq. (36) and the simplified expression for $b$ in Eq. (79),
becomes as follows. It is worth noting that this result is exactly the ideal
gas speed of sound term combined with the QLCA expression for the same pair
distribution function approximation [28].

$\begin{gathered}\left(\frac{\omega_{L}}{\omega_{p}}\right)^{2}(\Gamma_{0},|\vec{q}_{0}|)=\frac{1}{3}+\frac{5}{9}\frac{|\vec{q}_{0}|^{2}}{\Gamma_{0}}-2\frac{\cos\left(\sqrt{\frac{6}{5}}|\vec{q}_{0}|\right)}{\left(\sqrt{\frac{6}{5}}|\vec{q}_{0}|\right)^{2}}+2\frac{\sin\left(\sqrt{\frac{6}{5}}|\vec{q}_{0}|\right)}{\left(\sqrt{\frac{6}{5}}|\vec{q}_{0}|\right)^{3}}\end{gathered}$
(86)

The results for different values of $\Gamma_{0}$ are presented in Fig. 6(a),
and the following properties can be observed from the figure. As in both
previous approaches, for weak coupling where $\Gamma_{0}=1$, the dispersion
relation resembles the behavior of an ideal gas. However, as $\Gamma_{0}$
increases, the behavior undergoes a dramatic change. Similar to the Eulerian
approach, but unlike the original Lagrangian approach, there is a critical
value of $\Gamma_{0}=7.0$ at which we observe the onset of negative
dispersion. In other words, for small values of $|\vec{q}_{0}|$,
$\omega_{L}\leq\omega_{p}$ instead of $\omega_{L}\geq\omega_{p}$. However,
unlike the Eulerian approach but similar to the original Lagrangian approach,
there are no unstable regions. This means that for all $\Gamma_{0}$, we find
that $\omega_{L}^{2}>0$ for all $|\vec{q}_{0}|$.

Figure 6: The longitudinal dispersion law in the modified Lagrangian approach,
expressed in normalized variables, is shown for different values of
$\Gamma_{0}$ as given by Eq. (82). We assume the step function approximation
for the pair distribution function, and therefore, we use Eqs. (32), (36),
(40), (41), (79) for a given $u(\Gamma)$ defined in Eq. (30). In (a),
$u(\Gamma)=-0.9\Gamma$, and the result is provided by Eq. (86). In (b),
$u(\Gamma)=-0.9\Gamma+0.5944\Gamma^{1/3}-0.2786$.

The simpler fit for $u(\Gamma)$ does not account for thermal effects in the
nonlocal interaction. Therefore, we also consider a more precise fit,
$u(\Gamma)=-0.9\Gamma+0.5944\Gamma^{1/3}-0.2786$, which is accurate for
$\Gamma\geq 1$ and includes thermal effects [15].

The results for different values of $\Gamma_{0}$ for the transverse dispersion
law are presented in Fig. 4(b), and for the longitudinal dispersion law, they
are shown in Fig. 6(b). As in both previous approaches, one can observe that
the more complicated fit for $u(\Gamma)$ yields results that are qualitatively
the same as in the case of the simpler fit, with small quantitative
differences. For instance, the onset of negative dispersion is now observed at
$\Gamma_{0}=9.5$.

The results for the more precise fit are also compared to the results obtained
by QLCA [28] for the same approximation of the pair distribution function, as
well as to current fluctuation spectra obtained from molecular dynamics
simulations provided by the authors of Ref. 25, where dispersion laws are
identified by the peaks in the spectra.

The comparison for the longitudinal dispersion law is presented in Fig. 1. The
predicted dispersion law agrees with the molecular dynamics simulations. For
$\Gamma_{0}=1,5,10$, the agreement extends through the range of
$0<|\vec{q}_{0}|<5$. For the large value $\Gamma_{0}=80$, the agreement is
observed for the range of $0<|\vec{q}_{0}|<2$. Additionally, the predicted
onset of negative dispersion at $\Gamma_{0}=9.5$ aligns excellently with the
range of values from $\Gamma_{0}=9.5$ to $\Gamma_{0}=10.0$, estimated by the
molecular dynamics simulations [25, 27]. Moreover, the predicted longitudinal
dispersion law in the limit as $\Gamma_{0}\to\infty$ tends to the QLCA result,
which is considered accurate in such a limit [28]. This can be observed for a
simpler fit for $u(\Gamma)$ by taking the limit $\Gamma_{0}\to\infty$ directly
in Eq. (86), while for the more complicated fit for $u(\Gamma)$, it was
confirmed numerically. Consequently, in such a strong coupling limit, the
predicted longitudinal dispersion law, following the QLCA result, tends
towards the value $\omega_{p}/\sqrt{3}$ in the limit as
$|\vec{q}_{0}|=|\vec{k}|a_{0}\to\infty$.

The comparison for the transverse dispersion law is presented in Fig. 2.
Similar to the previous case of the Lagrangian approach, the predicted
dispersion law at all $\Gamma_{0}$ exactly agrees with the QLCA result and
aligns with the results of molecular dynamics simulations at large values of
wavevector. It correctly predicts that the transverse dispersion law tends to
the value $\omega_{p}/\sqrt{3}$ in the limit as
$|\vec{q}_{0}|=|\vec{k}|a_{0}\to\infty$. However, it fails to predict the
numerically observed disappearance of transverse waves at small values of the
wavevector [26]. Nonetheless, similar to the case of QLCA, this can be
remedied by considering relaxation in the context of generalized hydrodynamics
[26].

##  IV Discussion and Conclusion

We have explored variational principles for the hydrodynamics of the classical
OCP in both the Eulerian approach to hydrodynamics, where the hydrodynamic
functions depend on the position in the laboratory frame, and the Lagrangian
approach, where the hydrodynamic functions depend on some reference position,
such as the initial or equilibrium position of a fluid particle. We motivated
Lagrangian densities that are used in the variational principles and showed
how to obtain equations of motion and conservation laws for momentum and
energy when the variational principle is nonlocal. In the Eulerian approach,
one has to introduce constraint fields in the variational principle to obtain
continuity equations for number density and specific entropy. Because of that,
one has to be careful when eliminating the introduced constraint fields from
the resulting equations, as we have shown. In the Lagrangian approach, such
issues with constraint fields do not occur.

Consistency with equilibrium results from thermodynamics was ensured by
employing the energy conservation law and calculating the energy expression
from a given variational principle. This expression was evaluated in
equilibrium and matched with the known thermodynamic expression. Consequently,
all the considered variational principles yield the same result for
equilibrium energy, even though their out-of-equilibrium behavior differs due
to various assumptions regarding the pair distribution function. To fully
specify each of the considered variational approaches, we also assumed that
the out-of-equilibrium pair distribution function depends solely on number
densities and specific entropies in terms of $n,s$, rather than using more
complex models such as $(n+n^{\prime})/2,(s+s^{\prime})/2$ or
$\sqrt{nn^{\prime}},\sqrt{ss^{\prime}}$.

To analyze different proposed variational principles, we obtained linearized
equations of motion and calculated longitudinal and transverse dispersion
laws. In all of these approaches, the dispersion laws depend not only on the
pair distribution function but also on its adiabatic derivatives. This is
similar to the case of Euler equations, where dispersion laws depend on the
adiabatic derivatives of energy. However, in our case, we have included a
nonlocal contribution to the energy of the OCP that depends on the pair
distribution function. Nevertheless, we are not aware of precise numerical
data regarding the adiabatic derivatives of the pair distribution function.

Due to this lack of data, we used a simple step function approximation to the
pair distribution function. We chose the radius for the transition of the step
function to be consistent with the equilibrium energy expression. This
approach allows us to significantly simplify integrals that appear in the
dispersion laws. We hope that our accurate results will serve as motivation to
perform a more careful numerical analysis of the pair distribution function in
equilibrium in the future. With precise knowledge of its adiabatic
derivatives, we can avoid the simple step function approximation we used and
obtain even more accurate results.

In each of the considered variational approaches, dispersion laws were
computed using two fits for the nonlocal contribution to the energy of the
OCP. First, a simpler fit was considered, which is valid for large values of
the plasma parameter, $\Gamma\gg 1$, and for which simple analytical
expressions for the dispersion laws can be obtained, see Eqs. (42), (78),
(80), (86). Then, a more complicated fit that is valid for $\Gamma>1$ was
considered, for which dispersion laws were computed numerically. We found that
compared to the simpler fit with its simple analytical expressions, the
results are qualitatively the same, with some small quantitative differences.
This can be seen in Figs. 3, 4, 5, 6.

Finally, we compared the obtained dispersion laws for different variational
approaches across a wide range of equilibrium values of the plasma parameter,
$\Gamma_{0}$, to the results of theoretical calculations using QLCA and the
numerical results of molecular dynamics simulations. The results can be
observed in Figs. 1, 2. Indeed, we can see that different variational
approaches yield significantly different results for the dispersion laws,
emphasizing the importance of understanding how the pair distribution function
behaves out of equilibrium within the variational principles.

From the comparison of the transverse dispersion laws in Fig. 2, we can
observe that for small values of the wavevector, $|\vec{q}_{0}|$, the Eulerian
variational approach correctly predicts the absence of transverse waves in the
system. Meanwhile, for large values of the wavevector, the modified Lagrangian
approach agrees with the QLCA result and correctly predicts the presence of
transverse waves. As with QLCA, one can improve our models by introducing
relaxation, as done in generalized hydrodynamics, which allows for the
accurate prediction of the disappearance of transverse waves for small
wavevector values. The transverse dispersion law does not differentiate all of
the variational principles from the QLCA result. To understand the differences
between theoretical models, it is, therefore, better to examine the
longitudinal dispersion law, as shown in Fig. 1.

From the comparison of the longitudinal dispersion laws, we conclude that the
modified Lagrangian variational principle provides very accurate results for
all values of the wavevector, $|\vec{q}_{0}|$, and for all considered values
of $\Gamma_{0}$ when compared to the results of molecular dynamics
calculations. It is also more accurate than the other variational approaches
considered and QLCA. In the regime of large wavelengths and small wavevectors,
the modified Lagrangian variational approach correctly predicts the onset of
negative dispersion to occur at $\Gamma_{0}=9.5$, in agreement with the range
of values $9.5\leq\Gamma_{0}\leq 10.0$ predicted by the molecular dynamics
simulations. In the limit of strong coupling, where $\Gamma_{0}\to\infty$, the
modified Lagrangian variational approach tends towards the QLCA result. Thus,
in the regime of small wavelengths and large wavevectors, it predicts that
both longitudinal and transverse dispersion laws tend towards the finite value
$\omega_{p}/\sqrt{3}$ as $|\vec{q}_{0}|=|\vec{k}|a_{0}\to\infty$.

The excellent agreement of the dispersion laws predicted by the modified
Lagrangian variational principle, as compared to the molecular dynamics
simulations, shows that our variational hydrodynamic approach can be used to
describe motion on small length scales where the wavelength is comparable to
the average distance between the particles. This suggests that one can apply
this general variational approach with reasonable confidence in the future to
nonlinear problems of the OCP, such as the motion of vortices and adiabatic
expansion. It also suggests that one can apply it to other long-range systems,
for example, Yukawa fluids, two-component plasma, two-dimensional OCP, and
fluids consisting of electric or magnetic dipoles.

Determining the correct action indeed requires experimental input, and we have
shown how the procedure of developing theory using the variational principle
works. But once the correct input has been achieved, and the complete set of
variables has been assumed, one gets the payoff that a self-consistent theory
that can be applied to nonlinear processes has been obtained.

In our analysis, we have ignored the effects of viscosity, relaxation, and
heat transfer. In the future, one can consider adding them by including
additional terms in the equations of motion produced by the variational
principle. Viscosity can be included by adding a term related to the velocity
gradient, as seen in the case of the Navier-Stokes equations. Relaxation can
be included by adding a linear in velocity friction force between moving
particles and the stationary neutralizing background, as in the two-fluid
plasma equations. Alternatively, relaxation can be introduced for the
hydrodynamic variables for their approach to local thermodynamic equilibrium
values, as in generalized hydrodynamics. Heat transfer can be accounted for by
adding a term related to the temperature gradient to the energy conservation
law obtained from the variational principle, similar to the case of the
Navier-Stokes equations.

###### Acknowledgements.

This research has been funded by the AFOSR (FA9550-21-1-0295). The views,
opinions, and/or findings expressed are those of the authors and should not be
interpreted as representing the official views or policies of the Department
of Defense or the U.S. Government. We thank T. Chuna, Z. Johnson, M. S.
Murillo, L. G. Silvestri for many valuable discussions and Z. Donkó, I.
Korolov for the provided numerical data from the molecular dynamics
simulations.

## References

  * Brush, Sahlin, and Teller [1966] S. G. Brush, H. L. Sahlin, and E. Teller, “Monte Carlo Study of a One-Component plasma. I,” J. Chem. Phys. 45, 2102–2118 (1966).
  * Hansen [1973] J. P. Hansen, “Statistical Mechanics of Dense Ionized Matter. I. Equilibrium Properties of the Classical One-Component Plasma,” Phys. Rev. A 8, 3096–3109 (1973).
  * van Horn [1968] H. M. van Horn, “Crystallization of White Dwarfs,” Astrophys. J. 151, 227 (1968).
  * van Horn [2019] H. M. van Horn, “The crystallization of white dwarf stars,” Nat. Astron. 3, 129–130 (2019).
  * Demmel  _et al._ [2005] F. Demmel, S. Hosokawa, W.-C. Pilgrim, and S. Tsutsui, “Evidences for optic modes in molten NaI,” Nucl. Instrum. Methods Phys. Res. B 238, 98–101 (2005).
  * Demmel, Hosokawa, and Pilgrim [2021] F. Demmel, S. Hosokawa, and W.-C. Pilgrim, “Collective particle dynamics of molten NaCl by inelastic x-ray scattering,” J. Phys. Condens. Matter 33, 375103 (2021).
  * Hansen and McDonald [1975] J. P. Hansen and I. R. McDonald, “Statistical mechanics of dense ionized matter. IV. Density and charge fluctuations in a simple molten salt,” Phys. Rev. A 11, 2111–2123 (1975).
  * Bataller  _et al._ [2019] A. Bataller, A. Latshaw, J. P. Koulakis, and S. Putterman, “Dynamics of strongly coupled two-component plasma via ultrafast spectroscopy,” Opt. Lett. 44, 5832–5835 (2019).
  * Killian  _et al._ [2007] T. Killian, T. Pattard, T. Pohl, and J. Rost, “Ultracold neutral plasmas,” Phys. Rep. 449, 77–130 (2007).
  * Maxson  _et al._ [2013] J. M. Maxson, I. V. Bazarov, W. Wan, H. A. Padmore, and C. E. Coleman-Smith, “Fundamental photoemission brightness limit from disorder induced heating,” New J. Phys. 15, 103024 (2013).
  * Morfill and Ivlev [2009] G. E. Morfill and A. V. Ivlev, “Complex plasmas: An interdisciplinary research field,” Rev. Mod. Phys. 81, 1353–1404 (2009).
  * Landau and Lifshitz [1987] L. D. Landau and E. M. Lifshitz,  _Fluid Mechanics: Volume 6 (Course of Theoretical Physics)_ , 2nd ed. (Butterworth-Heinemann, 1987).
  * Desbiens, Arnault, and Clérouin [2016] N. Desbiens, P. Arnault, and J. Clérouin, “Parametrization of pair correlation function and static structure factor of the one component plasma across coupling regimes,” Phys. Plasmas 23, 092120 (2016).
  * DeWitt [1978] H. E. DeWitt, “Statistical mechanics of dense plasmas: Numerical simulation and theory,” J. Phys. Colloq. 39, C1–173–C1–180 (1978).
  * Khrapak and Khrapak [2014] S. A. Khrapak and A. G. Khrapak, “Simple thermodynamics of strongly coupled one-component-plasma in two and three dimensions,” Phys. Plasmas 21, 104505 (2014).
  * Landau and Lifshitz [1976] L. D. Landau and E. M. Lifshitz,  _Mechanics: Volume 1 (Course of Theoretical Physics)_ , 3rd ed. (Butterworth-Heinemann, 1976).
  * Landau and Lifshitz [1980a] L. D. Landau and E. M. Lifshitz,  _The Classical Theory of Fields: Volume 2 (Course of Theoretical Physics)_ , 4th ed. (Butterworth-Heinemann, 1980).
  * Herivel [1955] J. W. Herivel, “The derivation of the equations of motion of an ideal fluid by Hamilton’s principle,” Math. Proc. Camb. Philos. Soc. 51, 344–349 (1955).
  * Seliger and Whitham [1968] R. L. Seliger and G. B. Whitham, “Variational principles in continuum mechanics,” Proc. R. Soc. A 305, 1–25 (1968).
  * Lin [1963] C. C. Lin, “Hydrodynamics of helium II,” in  _Liquid helium: Proceedings of the international school of physics “Enrico Fermi”_ , Course XXI (Academic Press, 1963) pp. 93–146.
  * Putterman [1982] S. Putterman, “Comments on the variational principle and superfluid mechanics,” Phys. Lett. A 89, 146–148 (1982).
  * Murillo [2006] M. S. Murillo, “Ultrafast Dynamics of Strongly Coupled Plasmas,” Phys. Rev. Lett. 96, 165001 (2006).
  * Acciarri, Moore, and Baalrud [2022] M. D. Acciarri, C. Moore, and S. D. Baalrud, “Strong Coulomb coupling influences ion and neutral temperatures in atmospheric pressure plasmas,” Plasma Sources Sci. Technol. 31, 125005 (2022).
  * Hansen, McDonald, and Pollock [1975] J. P. Hansen, I. R. McDonald, and E. L. Pollock, “Statistical mechanics of dense ionized matter. III. Dynamical properties of the classical one-component plasma,” Phys. Rev. A 11, 1025–1039 (1975).
  * Korolov  _et al._ [2015] I. Korolov, G. J. Kalman, L. Silvestri, and Z. Donkó, “The Dynamical Structure Function of the One-Component Plasma Revisited,” Contrib. Plasma Phys. 55, 421–427 (2015).
  * Khrapak  _et al._ [2019] S. A. Khrapak, A. G. Khrapak, N. P. Kryuchkov, and S. O. Yurchenko, “Onset of transverse (shear) waves in strongly-coupled Yukawa fluids,” J. Chem. Phys. 150, 104503 (2019).
  * Mithen, Daligault, and Gregori [2012] J. P. Mithen, J. Daligault, and G. Gregori, “Onset of negative dispersion in the one-component plasma,” AIP Conf. Proc. 1421, 68–72 (2012).
  * Khrapak  _et al._ [2016] S. A. Khrapak, B. Klumov, L. Couëdel, and H. M. Thomas, “On the long-waves dispersion in Yukawa systems,” Phys. Plasmas 23, 023702 (2016).
  * Hansen [1981] J. P. Hansen, “Plasmon dispersion of the strongly coupled one component plasma in two and three dimensions,” J. Phys. Lett. 42, 397–400 (1981).
  * Kalman and Golden [1990] G. Kalman and K. I. Golden, “Response function and plasmon dispersion for strongly coupled Coulomb liquids,” Phys. Rev. A 41, 5516–5527 (1990).
  * Golden and Kalman [2000] K. I. Golden and G. J. Kalman, “Quasilocalized charge approximation in strongly coupled plasma physics,” Phys. Plasmas 7, 14–32 (2000).
  * Murillo [1998] M. S. Murillo, “Static local field correction description of acoustic waves in strongly coupling dusty plasmas,” Phys. Plasmas 5, 3116–3121 (1998).
  * Kaw and Sen [1998] P. K. Kaw and A. Sen, “Low frequency modes in strongly coupled dusty plasmas,” Phys. Plasmas 5, 3552–3559 (1998).
  * Diaw and Murillo [2015] A. Diaw and M. S. Murillo, “Generalized hydrodynamics model for strongly coupled plasmas,” Phys. Rev. E 92, 013107 (2015).
  * Jou  _et al._ [1985] D. Jou, C. Pérez-García, L. S. García-Colín, M. López de Haro, and R. F. Rodríguez, “Generalized hydrodynamics and extended irreversible thermodynamics,” Phys. Rev. A 31, 2502–2508 (1985).
  * Chen [2019] F. F. Chen,  _Introduction to Plasma Physics and Controlled Fusion_ , 3rd ed. (Springer, 2019).
  * Goldstein, Poole, and Safko [2001] H. Goldstein, C. Poole, and J. Safko,  _Classical Mechanics_ , 3rd ed. (Pearson, 2001).
  * Berestetskii, Lifshitz, and Pitaevskii [1982] V. B. Berestetskii, E. M. Lifshitz, and L. P. Pitaevskii,  _Quantum Electrodynamics: Volume 4 (Course of Theoretical Physics)_ , 2nd ed. (Pergamon, 1982).
  * Landau and Lifshitz [1980b] L. D. Landau and E. M. Lifshitz,  _Statistical Physics: Volume 5 (Course of Theoretical Physics)_ , 3rd ed. (Butterworth-Heinemann, 1980).
  * Chorin and Marsden [1993] A. J. Chorin and J. E. Marsden,  _A Mathematical Introduction to Fluid Mechanics_ , 3rd ed. (Springer, 1993).
  * Frahm [1983] C. P. Frahm, “Some novel delta-function identities,” Am. J. Phys. 51, 826–829 (1983).