# Magnetic excitations, phase diagram and order-by-disorder in the extended
triangular-lattice Hubbard model

Josef Willsher  0000-0002-6895-6039<EMAIL_ADDRESS>Department of Physics
TQM, Technische Universität München, James-Franck-Straße 1, D-85748 Garching,
Germany  Hui-Ke Jin  0000-0002-0880-6557 Department of Physics TQM, Technische
Universität München, James-Franck-Straße 1, D-85748 Garching, Germany
Johannes Knolle  0000-0002-0956-2419 Department of Physics TQM, Technische
Universität München, James-Franck-Straße 1, D-85748 Garching, Germany  Munich
Center for Quantum Science and Technology (MCQST), 80799 Munich, Germany
Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom

###### Abstract

The dynamical structure factor is an important observable of quantum magnets
but due to numerical and theoretical limitations, it remains a challenge to
make predictions for Hubbard-like models beyond one dimension. In this work,
we study the magnetic excitations of the triangular lattice Hubbard model
including next-nearest neighbor hopping. Starting from the 120∘ and stripe
magnetic orders we compute the magnon spectra within a self-consistent random
phase approximation. In the stripe phase, we generically find accidental zero
modes related to a classical degeneracy known from the corresponding
$J_{1}$–$J_{2}$ Heisenberg model. We extend the order-by-disorder mechanism to
Hubbard systems and show how quantum fluctuations stabilize the stripe order.
In addition, the frustration-induced condensation of magnon modes allows us to
map out the entire phase diagram which is in remarkable agreement with recent
numerical works. We discuss connections to experiments on triangular lattice
compounds and the relation of our results to the proposed chiral spin liquid
phase.

##  I Introduction

The Hubbard model has played a central role in condensed matter research as
one of the simplest effective models describing interacting electrons [1, 2].
It displays a whole variety of fascinating many-body phenomena such as high-
temperature superconductivity, Mott metal-insulator transitions, and different
forms of quantum magnetism [3]. In one dimension, the ground state and
excitations of the Hubbard model have been understood in great detail [4, 5,
6]. However, in the absence of exact methods, in two dimension and beyond the
ground state phase diagram remains a challenge even for state-of-the-art
numerical methods [7, 8, 9], with access to the excitaion spectrum proving
even more difficult. In addition, on geometrically frustrated lattices the
complexity is even higher [10, 11, 12, 13] since competing magnetic
interactions provide a potential for realizing quantum spin liquid (QSL)
phases [14, 15, 16, 17, 18, 19, 20, 21]. In particular, the successful
synthesis of possible QSL candidates [22] on the triangular lattice, including
organic compounds $\kappa$(ET)2Cu2(CN)3 and Me3EtSb[Pd(dmit)2]2 [23, 24], or
rare-earth compounds YbMgGaO4 [25, 26] and NaYbO2 [27] has triggered a lot of
activity. More recently, the rapid progress in the field of twisted moiré
materials [28, 29, 30, 31, 32] provides a versatile platform for the study of
Hubbard physics on various lattices and with different band structures [33,
34]. In particular, the triangular-lattice Hubbard model can be realized in
twisted bilayer boron nitride [35, 36], twisted WSe2/WS2 moiré superlattices
[37, 38], and twisted double-bilayer WSe2 [39].

Much progress has been made in the study of the Hubbard model on the
triangular lattice by using combinations of theory and numerics. In the
large-$U$ limit and with only nearest neighbor (NN) hopping $t$, the Hubbard
model at half filling is reduced to the antiferromagnetic Heisenberg model for
which the celebrated resonating valence bond (RVB) state was proposed as a
trial wave function by Anderson in 1973 [14] — nowadays considered to be the
first QSL state. However, the following numerical studies provided clear
evidence for a long-range 120∘ antiferromagnetic (AFM) order in the strong
coupling limit [40, 41, 42, 43, 44]. Away from strong coupling, a spin
stiffness analysis based on a mean-field approximation indicates a loss of
120∘ magnetic order at $U\sim{}6t$ [45]. An intermediate QSL phase (the so-
called weak-Mott insulator) might additionally be stabilized for moderate
Hubbard interaction $U\sim{}8t$, as indicated by the variational cluster
approximation (VCA) [46], path-integral renormalization group analysis [47],
and series expansions [48]. Moreover, recent density matrix renormalization
group (DMRG) studies give strong evidence for a gapped chiral QSL in this
regime [49, 50, 51, 52, 53, 54, 55].

The Hubbard model with additional next-nearest-neighbor (NNN) hoppings
$t^{\prime}$ is further frustrated and studies have mostly been restricted to
the large-$U$ model, i.e., the $J_{1}$–$J_{2}$ AFM Heisenberg model on the
triangular lattice. In the classical limit where the spins are treated as
$O(3)$ vectors, it is predicted that the coplanar 120∘ order is destabilized
and a highly degenerate four-sublattice state appears between
$1/8\leq{}J_{2}/J_{1}\leq{}1$ [56]; for larger values of $J_{2}/J_{1}>1$
spiral states dominate. The intermediate regime has been one of the
paradigmatic examples of the order-by-disorder mechanism [57, 58] where
quantum fluctuations treated within spin-wave theory lift the classical
degeneracy such that a collinear stripe ordered phase is stabilized [56, 59].
Near the classical critical point at $J_{2}/J_{1}=1/8$ an exciting possibility
is the emergence of a magnetically disordered phase, such as the $U(1)$ Dirac
QSL state which has been proposed within the variational Monte Carlo (VMC)
approach [60, 61]. DMRG works have suggested gapped [62, 63] and more recently
gapless [64, 65, 66] QSL ground states.

Only few works have studied the full Hubbard model with both the NN and NNN
hopping at finite $U$ including charge and spin fluctuations. A ground-state
phase diagram has been obtained in the $t^{\prime}$-$U$ plane by using a VCA
method with a 12-site cluster [67] and VMC [68], which together point to a
rich phase diagram including insulating spin disordered phases, similar to
other extensions of the Hubbard model [69, 70, 71, 72, 73]. However,
theoretical studies of the dynamical spin excitations as well as an
understanding of the metal-insulator and magnetic phase transitions are
missing, despite being of particular importance in light of the recent above-
mentioned candidate materials.

Despite recent breakthroughs in numerical techniques which allow ever more
detailed studies of Hubbard-like models, strong limitations on system size and
frequency resolution remain for the calculation of dynamical properties in two
dimensions. Specifically, DMRG studies of two-dimensional magnets are limited
in their prediction of dynamical correlation functions due to entanglement
growth in time. Recent effort has begun to shine light on the dynamical
properties of the strong coupling Heisenberg limits [66, 74], but is currently
unable to probe Hubbard physics in the full $t/U$ regime. In this work we
provide calculations of dynamical properties of the triangular lattice
$t$–$t^{\prime}$–$U$ Hubbard model. Working within a self-consistent mean
field theory we are able to calculate the dynamical spin-spin correlation
function within the random phase approximation (RPA). Our method, which
includes harmonic spin and charge fluctuations around the magnetically ordered
states, has been pioneered for the square lattice Hubbard model [75, 76, 77]
and successfully applied to describe inelastic neutron scattering experiments
for cuprates [78] and iron-based superconductor [79, 80, 81, 82] parent
compounds.

We find broadly overlapping self-consistent solutions for both the 120∘\- and
stripe-ordered phases. Instead of comparing the respective ground state
energies to produce a pure mean-field phase diagram, we use the calculation of
spin excitations to examine the stability of the mean-field magnetization in
the presence of fluctuations. Most of the phase boundaries of the 120∘ and
stripe phases can be mapped out by observing the condensation of _accidental
soft modes_ when the magnon spectrum comes down to zero-energy at a wavevector
incommensurate with the existing magnetic order. We additionally find a novel
phase boundary at small positive $t^{\prime}$ and $U$ where magnetic order is
instead destabilized directly by a vanishing of the spin wave velocity. From
the components of the structure factor we observe that the corresponding
fluctuations are out-of-plane indicating an instability to a phase with non-
collinear spin correlations. Of course, our method is only able to describe
excitations within a long-range magnetically ordered phase but we argue that
the unusual out-of-plane instability is an indication of the nearby chiral
spin liquid phase.

In addition to the phase diagram, we provide a detailed understanding of the
magnetic excitations in the stripe phase. Similar to the corresponding
$J_{1}$–$J_{2}$ model, for a stripe state with wave vector $M$ we find within
our Hubbard model approach accidental soft modes at $M^{\prime}$ as a direct
signature of the underlying magnetic degeneracy for sizable $t^{\prime}$. We
then show that the inclusion of quantum fluctuations via an approximate
treatment of self-energy effects gaps out the accidental mode at $M^{\prime}$
stabilizing the collinear stripe phase. Thus, we extend the order-by-disorder
mechanism to the triangular lattice Hubbard model.

Our manuscript is structured as follows: In Section II we present the model
and its mean-field solution. In Section III we explain the calculation of the
magnetic phases’ spin fluctuations in the RPA scheme. In Section IV we discuss
in detail the stability of the stripe ordered phase beyond harmonic
fluctuations; and in Section V discuss how phase transitions of the model are
diagnosed and used to map out the phase diagram. We conclude with a discussion
of our results and open questions.

##  II Mean field theory

The Hubbard model Hamiltonian is given by

$H=-\sum_{i,\delta}t_{\delta}\,a^{\dagger}_{i\sigma}a_{i+\delta\sigma}+U\sum_{i}a^{\dagger}_{i\uparrow}a_{i\uparrow}\,a^{\dagger}_{i\downarrow}a_{i\downarrow},$
(1)

where $\delta$ is a nearest or next-nearest neighbor displacement vector on
the triangular lattice, as shown in Fig. 1(a), and $a^{\dagger}_{i\sigma}$ is
a creation operator of a fermion with spin $\sigma$ at site $i$. The (next)
nearest neighbor hopping is henceforth referred to as $t$ ($t^{\prime}$). The
on-site repulsion is $U$ and the chemical potential is fixed to half filling.

Figure 1: Model and Phases: (a) Triangular lattice with nearest neighbor $t$
(solid) and next-nearest neighbor $t^{\prime}$ hoppings (dashed line)
highlighted. (b) Brillouin zone showing the highlighted magnetic Brillouin
zones for each ordered phase. The path drawn is used to plot spectral
functions. (c) 120∘ (left) and stripe (right) ordered states. Unit cells
highlighted in grey boxes and magnetic moments $\Delta_{\alpha}$ highlighted
for each order. (d) Phase diagram of the triangular lattice Hubbard model,
calculated using mean field theory plus fluctuations with on-site repulsion
$U$. The phases highlighted are: central 120∘ order (blue); peripheral stripe
order (green); at low $U$ a metallic phase (gray); intermediate insulating
spin-disordered phases (white); and a potential noncoplanar insulating state
(pink). Red phase boundary signifies the vanishing spin stiffness at the
transition due to out-of-plane fluctuations.

First, we review the standard self-consistent mean-field theory for the to two
different magnetic orderings on the triangular lattice: 120∘ order at the $K$
point in the Brillouin zone (BZ); and stripe order at the $M$ point. Fig. 1(b)
shows the first BZ of a triangular lattice with the ordering vectors $K$ and
$M$ highlighted, alongside other high-symmetry points and a path between them
which is used for our analysis of dynamical spectral functions across the BZ.
In the presence of symmetry breaking and magnetic order formation, the
magnetic unit cell is enlarged; it contains $N_{\mathrm{UC}}=3$ sites for the
120∘ and $N_{\mathrm{UC}}=2$ sites for stripe order with the magnetic ordering
pattern highlighted in Fig. 1(c). The respective magnetic BZ (MBZ) is shrunk
covering 1/3 and 1/2 of the BZ area [Fig. 1(b)]. One may consider equivalent
stripe orderings at the other midpoints of the BZ edge (labeled $M^{\prime}$),
which are related by a sublattice-dependent spin rotation and are part of the
underlying classical manifold of magnetic states [83].

With these magnetic states we perform a Hartree–Fock (HF) decoupling of the
Hamiltonian Eq. (1) which is implemented by introducing _sublattice fermions_
$a_{\alpha}$, where the sublattice index $\alpha=A,B,C$ for the 120∘ Ansatz
and $\alpha=A,B$ for the stripe Ansatz. The mean-field Hamiltonian elements
are hence labeled by sublattice index $\alpha$ in addition to spin indices
$\sigma=\uparrow,\downarrow$.

The ordering is encoded in this formalism with the sublattice vectors
$\vec{\Delta}_{\alpha}$ (components $(\Delta_{\alpha})_{\mu}$, for
$\mu=x,y,z$). In our calculations, we choose the following basis for mean-
field vectors as highlighted in real space in Fig. 1(c): the 120∘ state has
$\vec{\Delta}_{A}=(0,\Delta,0)$ with $\vec{\Delta}_{B,C}$ rotated clockwise by
$120$ and $240$ degrees; the stripe state has
$\vec{\Delta}_{A,B}=\pm(0,0,\Delta)$, where we have used the fact that the
mean-field collinear state is independent of the choice of ordering direction
to fix it out of the plane. Written in this compact form, the decoupled
Hamiltonian is

$H_{\mathrm{MF}}=\sum_{i,\alpha}a_{i,\alpha\sigma}^{\dagger}[-(\Delta_{\alpha})_{\mu}(\sigma^{\mu})_{\sigma\sigma^{\prime}}]a_{i,\alpha\sigma^{\prime}}+\sum_{i,\delta}t_{\delta}\,a^{\dagger}_{i,\sigma}a_{i+\delta,\sigma},$
(2)

where the repeated spin indices are implicitly summed over. In addition to
this, the sublattice magnetization vector is defined self-consistently

$\Delta_{\alpha,\mu}=U\left<(S_{\alpha})_{\mu}\right>_{\mathrm{GS}}=\frac{U}{2}\left<a_{\alpha\sigma}^{\dagger}(\sigma_{\mu})_{\sigma\sigma^{\prime}}a_{\alpha\sigma}\right>_{\mathrm{GS}}$
(3)

where $\vec{S}_{\alpha}$ is the sublattice spin vector and expectation value
$\left<~{}\cdot~{}\right>_{\mathrm{GS}}$ is with respect to the ground state.
The sublattice magnetization $m=2\Delta/U$ is calculated self consistently for
each state.

We next block-diagonalize the system by Fourier-transformation; introducing
the sublattice spinor
$\Psi_{\mathbf{k}}=([a_{\mathbf{k}}]_{\uparrow,A},[a_{\mathbf{k}}]_{\downarrow,A},\dots)^{T}$
with 6 (or 4) compound indices $m=(\alpha,\sigma)$ for 120∘ (stripe)
sublattice fermions, we can rewrite the Hamiltonian as
$H_{\mathrm{MF}}=\sum_{\mathbf{k}}\Psi^{\dagger}_{\mathbf{k}}H(\mathbf{k})\Psi_{\mathbf{k}}$.
Explicitly the mean-field Hamiltonian in the stripe phase is therefore

$H_{\mathrm{stripe}}(\mathbf{k})=\begin{pmatrix}\zeta_{\mathbf{k}}\mathbb{I}-\sigma\cdot\Delta_{A}&\delta_{\mathbf{k}}\mathbb{I}\\\
\delta_{\mathbf{k}}^{*}\mathbb{I}&\zeta_{\mathbf{k}}\mathbb{I}-\sigma\cdot\Delta_{B}\end{pmatrix},$
(4)

with $\mathbb{I}$ the identity matrix and $AA,BB$ and $AB,BA$ hoppings written
as $\zeta_{\mathbf{k}}$ and $\delta_{\mathbf{k}}$, respectively. These are
calculated by summing over all $AA$ or $AB$ hopping vectors $\delta$ for both
NN and NNN bonds, defined in the following way:

$\delta_{\mathbf{k}}=\sum_{\delta=A\to
B}t_{\delta}e^{i\mathbf{k}\cdot\delta},\quad\zeta_{\mathbf{k}}=\sum_{\delta=A\to
A}t_{\delta}e^{i\mathbf{k}\cdot\delta}.$ (5)

Explicit expressions for these matrix elements are given in the Appendix
alongside the equivalent formulation for the 120∘ mean-field Hamiltonian.

Diagonalizing the HF Hamiltonian produces eigenstates $\ket{\mathbf{k},m}$ of
the Hamiltonian with spectrum $\omega_{m}(\mathbf{k})$, where $m=1\dots
2N_{\mathrm{UC}}$. Mean field expectation values can be computed in this
basis, including the ordering vector via Eq. (3); in a system of linear
dimension $L$, its expectation value is given

$\displaystyle\Delta_{\alpha,\mu}$
$\displaystyle=\frac{U}{2}\left<\Psi^{\dagger}\sigma^{\mu}_{\alpha}\Psi\right>_{\mathrm{GS}}$
$\displaystyle=\frac{U}{2N}\sum_{\mathbf{k},m}n_{f}[\omega_{m}(\mathbf{k})-\mu]\bra{\mathbf{k},m}\sigma_{\alpha}^{\mu}\ket{\mathbf{k},m}$
(6)

where $N=L^{2}N_{\mathrm{UC}}$ is the number of sites and
$n_{f}[\omega_{m}(\mathbf{k})-\mu]$ is the Fermi–Dirac distribution at zero
temperature. The sublattice-$\alpha$ Pauli matrix $\sigma_{\alpha}^{\mu}$ is
defined as acting on a $\Psi^{\sigma^{\prime}}_{\beta}$ object through the
Pauli matrix $(\sigma^{\mu})_{\sigma\sigma^{\prime}}$ on its spin index and on
its sublattice index diagonally via $\delta^{\alpha\beta}$. At half-filling,
the chemical potential $\mu(U,\Delta)$ is always fixed such that the fermion
occupation

$n=\left<\Psi^{\dagger}\Psi\right>_{\mathrm{GS}}=\frac{1}{N}\sum_{\mathbf{k},m}n_{f}[\omega_{m}(\mathbf{k})-\mu]$
(7)

is equal to $1$ per site.

Numerically, the magnetic order $\Delta$ is computed self consistently as
follows: for a given choice $\Delta$ the mean-field Hamiltonian is found, then
the chemical potential $\mu$ is fixed by minimizing $n-1$ numerically, and
finally the expectation of magnetic ordering [found via Eq. (II)] is computed.
The difference
$\left<\Psi^{\dagger}\sigma^{\mu}_{\alpha}\Psi\right>_{\mathrm{GS}}-2\Delta/U$
is then minimized numerically to self-consistently compute $\Delta$.

Figure 2: Dynamical structure factor $S(\mathbf{q},\omega)$ of (a) the 120∘
and (b) stripe magnetic ordered phases, taken along the path in the Brillouin
zone highlighted in the text. Dotted line represents the zero of
$[1-U\chi^{(0)}(\mathbf{q},\omega)]$. Systems $U=30$, linear size $L=720$
($N=N_{\mathrm{UC}}L^{2}$ sites), with (a) $t^{\prime}=0$ and (b)
$t^{\prime}=0.6$.

##  III Spin fluctuations

The spin susceptibility is one of the prime observables for understanding the
physics of Hubbard systems: This quantity gives access to the _dynamical
structure factor_ , an experimentally observable quantity which is routinely
probed in inelastic neutron scattering experiments [84]. As Goldstone modes in
an AFM system, the spin waves extend down to low energy with linear
dispersion, and understanding the velocity of these modes at low energy allows
the prediction of quantum corrections to the sublattice magnetization. Access
to the dynamical spin susceptibility will hereby provide information about
phase boundaries beyond the basic mean-field approximation. The high-energy
spectrum of spin excitations may in particular show strong deviations from
Heisenberg-like spin-only physics, due to the effect of charge fluctuations
especially in the intermediate-$U$ regime.

The dynamical spin susceptibility $\chi(\mathbf{q},\omega)$ can be expressed
in terms of the momentum-space spin-spin correlation function [85],

$[\chi(\mathbf{q},\omega)]^{\mu\nu}_{\alpha\beta}=\frac{i}{2N}\int\differential{t}\,e^{i\omega
t}\left<TS_{\mathbf{q},\alpha}^{\mu}(t)S_{-\mathbf{q},\beta}^{\nu}(0)\right>_{\mathrm{GS}},$
(8)

defined with sublattice indices $\alpha,\beta$ and spin indices $\mu,\nu$.
Explicitly, the spin operators can be written as
$S_{\mathbf{q},\alpha}^{\mu}=\sum_{\mathbf{k}}\Psi_{\mathbf{k}+\mathbf{q}}^{\dagger}[\sigma^{\mu}_{\alpha}]\Psi_{\mathbf{k}}$.

In order to calculate this susceptibility, we focus on the magnetic
interaction channel and apply the RPA formalism [75, 76, 82]. The
susceptibility tensor is defined through the Dyson equation

$[\chi]^{\mu\nu}_{\alpha\beta}=[\chi]^{\mu\rho}_{\alpha\gamma}U^{\rho\sigma}_{\gamma\lambda}[\chi^{(0)}]^{\sigma\nu}_{\lambda\beta}+[\chi^{(0)}]^{\mu\nu}_{\alpha\beta},$
(9)

where the interactions are contained in the tensor $U$. This is
diagrammatically expressed as

(10)

in terms of the non-interacting correlation function
$\chi^{(0)}(\mathbf{q},\omega)$. The interaction vertex is taken to be
diagonal in both indices

$U^{\rho\sigma}_{\gamma\lambda}=U\mathbb{I}=\begin{gathered}\includegraphics[scale={0.66}]{vertex.pdf}\,.\end{gathered}$
(11)

The solution to the equation can be conveniently written by summing up all
bubble diagrams [77, 86]; in matrix form this is

$\chi(\mathbf{q},\omega)=\left[\mathbb{I}-U\chi^{(0)}(\mathbf{q},\omega)\right]^{-1}\chi^{(0)}(\mathbf{q},\omega).$
(12)

Our task therefore becomes to evaluate this bare correlation function using
the mean-field state; it can be expressed as

$\displaystyle[\chi(\mathbf{q},\omega)^{(0)}]^{\mu\nu}_{\alpha\beta}=\frac{1}{N}\sum_{\mathbf{k}}$
$\displaystyle\sum_{m,l}^{2N_{\mathrm{UC}}}[n_{f}(\omega_{\mathbf{k},m})-n_{f}(\omega_{\mathbf{k}+\mathbf{q},l})]$
$\displaystyle\times$
$\displaystyle\frac{[U_{\mathbf{k},m}]^{\nu}_{\alpha}[U_{\mathbf{k}+\mathbf{q},l}^{*}]^{\mu}_{\beta}}{\omega+\omega_{\mathbf{k},m}-\omega_{\mathbf{k}+\mathbf{q},l}-i\delta},$
(13)

where $[U_{\mathbf{k},m}]^{\mu}_{\alpha}$ is the $m$ eigenvector of the
Hamiltonian $H^{\mu\alpha}_{\mathbf{k}}$. The divergence of the susceptibility
is regulated by a width $i\delta$, which we taken proportional to $1/(NU)$ so
as to scale with the level spacing in a finite system and to smoothly take the
thermodynamic and strong coupling limits.

The bare correlation function has poles corresponding to particle-hole
excitations, which are gapped for half-filling in all magnetic phases. The
collective modes are contained in the zeros of
$[1-U\chi^{(0)}(\mathbf{q},\omega)]$ corresponding to the spin wave
excitations in a fully self-consistent treatment. In order to reproduce the
dynamical spectral function, we calculate the full matrix
$\chi(\mathbf{q},\omega)$ and then evaluate the following sum over spins and
sublattices

$S(\mathbf{q},\omega)=\operatorname{Im}\sum_{\alpha,\beta}^{N_{\mathrm{UC}}}\,e^{-\mathbf{q}\cdot(r_{\alpha}-r_{\beta})}\sum_{\mu,\nu=x,y,z}[\chi(\mathbf{q},\omega)]^{\mu\nu}_{\alpha\beta}.$
(14)

Moreover, the eigenvalues of $1-U\chi^{(0)}(\mathbf{q},\omega)$ are evaluated
for each wavevector $\mathbf{q}$, and the zeros are plotted to show the spin
wave dispersion.

Fig. 2(a) shows a typical structure factor $S(\mathbf{q},\omega)$ evaluated in
the 120∘ ordered phase away from the phase boundaries. The ordering wavevector
is at the $K$ points in the BZ, and spectral weight vanishes at the $\Gamma$
point in agreement with previous results [45], in which only the
$t^{\prime}=0$ case was studied. This coplanar but non-collinear state has
three spin excitation modes; at low energies the susceptibility matrix is
approximately diagonal in global coordinates $\mu=x,y,z$, allowing us to
characterize each of these modes as having either purely in-plane
($S^{x},S^{y}$) or out-of-plane ($S^{z}$) character. Magnon excitations, which
are transverse to the local ordering vector in each sublattice, hence contain
components in all three spin components. The modes are distinguished by
calculating $\operatorname{Im}\chi^{zz}$ for out-of-plane modes, and
$\operatorname{Im}\chi^{+-}$ for in-plane modes (interpreted as being
$\left<S^{+}S^{-}\right>$ and depending on the matrix elements
$\chi^{xx}+\chi^{yy}+i\chi^{xy}$). These Goldstone modes have a linear
dispersion at low $\omega$, where the single in-plane mode has a higher
velocity than the two degenerate transverse modes $c_{\parallel}>c_{\perp}$.

Looking to the high-energy behavior, the in- and out-of-plane modes generally
mix character but remain degenerate across the $M\to Y\to\Gamma$ line and at
the other high symmetry points $Y_{1},M^{\prime}$ where the highest-energy
bands meet and where the fluctuations are predominantly in-plane.

Fig. 2(b) shows a typical structure factor $S(\mathbf{q},\omega)$ evaluated in
the stripe ordered phase away from the phase boundaries. The ordering
wavevector here is at the $M$ point in the BZ with linear dispersion; a
similar linear mode also appears at the $\Gamma$ point but with vanishing
spectral weight. There is only one excitation mode in this collinear state
which is twofold degenerate everywhere in the BZ.

A further advantage of the self-consistent RPA approach is that it allows the
calculation of dynamical response functions away from half filling. The method
naturally includes charge fluctuations from particle-hole excitations which
cannot be captured within a spin-only treatment. We have checked that both the
120∘ and stripe phases are stable against hole doping; the two-particle
continuum is visible as weak spectral weight, and there is an increase in the
spin-wave velocity for low doping. We also confirmed a strong particle-hole
asymmetry, as was seen for the $t^{\prime}=0$ model in Ref. [45], but since a
full study of the doping dependent phase diagram is beyond the scope of this
work we concentrate on the case of half-filling for the remainder.

##  IV Order-by-Disorder in the Stripe Phase

Turning back to the stripe phase spectrum in Fig. 2(b), one can observe a
peculiar mode at the point $M^{\prime}$ which goes quadratically to zero
despite the magnetic ordering wavevector being at the inequivalent $M$ point.
In the large-$U$ limit, the $t,t^{\prime}$ Hubbard model becomes the
corresponding frustrated Heisenberg $J_{1}$–$J_{2}$ model, whose stripe order
has an exact _accidental zero mode_ according to linear spin wave theory
(LSWT) [59, 56]. Taken at face value, the zero-energy spin excitations at this
point with zero stiffness would produce a divergent quantum correction to the
sublattice magnetization, destroying the stripe state, thus, at the harmonic
level no long-range order is selected. However, this is an artifact of the
quadratic-fluctuation approximation; in the Heisenberg model it is known that
quantum fluctuations generically gap the accidental zero mode at $M^{\prime}$
to higher energies — the gap becomes of the order $J$ — leading to a stable
phase while the Goldstone modes at the $M$ point remain gapless [59]. The
general scenario that quantum (or thermal or quenched disorder) fluctuations
stabilize a long-range ordered state is known as _order-by-disorder_.

Figure 3: Stripe phase instability in the RPA scheme of the accidental zero
mode around the $M^{\prime}$ point. Subfigures (a)–(c) plot the transverse
susceptibility $\operatorname{Im}\chi^{+-}(\mathbf{q},\omega)$ for $U=30$,
$20$, and $10$. The white dashed line is the mode $\omega_{\mathbf{q}}$
calculated from the zeros of $1-U\chi^{(0)}$. Calculations taken
$t^{\prime}=0.8$ and linear size $L=720$ ($N=2L^{2}$ sites). Subfigure (d)
shows the energy of the mode minimum $\omega_{M^{\prime}}$ for varying values
of $t^{\prime}$ across the phase diagram has a cubic dependence on $1/U$.

Turning back to the Hubbard model’s stripe phase in the RPA approximation, we
find that reducing $U$ actually tends to lower the minimum of the mode
$\omega_{M^{\prime}}$ below zero energy. Fig. 3(a)–(c) plots the transverse
susceptibility $\operatorname{Im}\chi^{+-}$ for varying $U$, showing that the
minimum is significantly below zero energy at $U\approx 15$ with positive
$t^{\prime}$ (and happens sooner for negative $t^{\prime}$). In Fig. 3(d) this
is put quantitatively, showing the $U^{-3}$ dependence of this minimum energy
in the RPA scheme.

The presence of the accidental zero mode at $M^{\prime}$ is again an artifact
of the harmonic fluctuation approximation scheme, and will not hold up in the
presence of quantum fluctuations. Indeed, the square lattice antiferromagnetic
and ferromagnetic phases were studied in Refs. [86, 87], where quantum
fluctuations were calculated systematically for the Hubbard model by extending
it to contain $\mathcal{N}$ orbitals. By evaluating the fermion self-energy to
order $1/\mathcal{N}$, the effect on the susceptibility was shown to be a
shifted pole with renormalized weight. In the following, we provide an
approximate treatment of quantum fluctuations around the $M^{\prime}$ point on
the triangular lattice stripe phase and show that the accidental mode remains
gapped down to low $U$. We use the fact that the effect of quantum
fluctuations on $\chi$ can be calculated by replacing $\chi^{(0)}$ in Eq. (12)
with $\phi=\chi^{(0)}+\delta\phi$, where the additional term originates from
diagrams which appear alongside the free loop $\chi^{(0)}$ in the Dyson
equation [86].

We can calculate an approximate form of $\phi$, generally given by

$\delta\phi(\mathbf{k},\omega)=\frac{1}{\mathcal{N}}\frac{t^{2}}{U^{3}}\begin{pmatrix}K_{\mathbf{k}}+J_{\mathbf{k}}\,U\omega/2t^{2}&K_{\mathbf{k}}\gamma_{\mathbf{k}}\\\
K_{\mathbf{k}}\gamma_{\mathbf{k}}&K_{\mathbf{k}}-J_{\mathbf{k}}\,U\omega/2t^{2}\end{pmatrix}$
(15)

where $K_{\mathbf{k}}$, $J_{\mathbf{k}}$ are factors coming from loop diagrams
[86]. These are momentum dependent and of the order $1/\mathcal{N}$ in the
controlled expansion parameter. The term $\gamma_{\mathbf{k}}$ is the off-
diagonal element of $\chi^{(0)}$, which is zero when evaluated at the
$M^{\prime}$ point in the RPA. We are therefore able to simply approximate the
effect of the additional diagrams for small momentum around $M^{\prime}$ by
taking $K_{\mathbf{k}}=J_{\mathbf{k}}$ as an order $1$ constant, and
$\gamma_{\mathbf{k}}=0$. In Fig. 4(a) we show that indeed the spectrum is
corrected to be gapped at the $M^{\prime}$ point (for concreteness we used
$U=20$, $t=0.8$). Thus, we find that similar to the Heisenberg limit,
fluctuations also stabilize the triangular lattice Hubbard model’s stripe
phase.

Figure 4: Stripe phase stability including quantum fluctuations of the
accidental zero mode around the $M^{\prime}$ point. Subfigure (a) plots the
transverse susceptibility $\operatorname{Im}\chi^{+-}(\mathbf{q},\omega)$ for
$U=20$. The white dashed line is the mode $\omega_{\mathbf{k}}$ calculated
from the zeros of $1-U\phi$. Calculations taken $t^{\prime}=0.8$ and linear
size $L=720$ ($N=2L^{2}$ sites). Subfigure (b) shows the energy of the mode
minimum $\omega_{M^{\prime}}$ for varying values of $t^{\prime}$ including the
phenomenological self-energy correction. Scaling at high $U$ is linear, with
the dashed line $\propto t^{2}/U$.

In this corrected picture and for large on-site repulsion $U$, we find that
the gap magnitude scales as $\omega_{M^{\prime}}\propto t^{2}/U$ in agreement
with the Heisenberg limit [59] (the exact constant relating the gap to $J$ in
this limit depends on the values $J_{M^{\prime}}$, $K_{M^{\prime}}$). For
positive $t^{\prime}$ this stability extends down to below the phase
transition to a spin-disordered phase, suggesting the order-by-disorder effect
stabilizes the whole stripe phase here.

For negative $t^{\prime}$, the $\omega_{M^{\prime}}$ gap appears to vanish
again at around $U\approx 20$, where the correction $\propto t^{2}/U$ becomes
comparable to the $\omega_{M^{\prime}}$ of the quadratic RPA theory, see Fig.
4 [and the dashed line in Fig. 3(d)]. We suspect that in this regime either
the phase becomes disordered, or that diagrams involving charge fluctuations
become significant and may further stabilize the phase. One could further
improve our approximate treatment by including higher-order self-energy
corrections and by calculating these self-consistently [88] which comes at a
considerable numerical cost. This is beyond the scope of our work but we argue
that our results provide good evidence that for much of the coupling regime
$U\gtrsim 10$ supporting stripe ordering, the magnon spectrum at the
$M^{\prime}$ point is gapped. In the following, we will therefore ignore the
artificial zero mode at the $M^{\prime}$ point within the harmonic
approximation and instead focus on the softening of collective modes at other
points in the BZ for mapping out the phase diagram.

##  V Phase Transitions & Discussion

We can now map out the phase diagram of the next-nearest neighbor triangular
lattice Hubbard model, going beyond mean-field theory by requiring stability
against fluctuations. Our main result is shown in Fig. 1(d) which is in
remarkable in agreement with numerical works based on the variational cluster
approximation [67] and more recently on the variational Monte Carlo method
[68]. We find the presence of 120∘ order at low $|t^{\prime}|$, stripe order
for increased $|t^{\prime}|$, and instabilities towards spin disordered phases
in the intermediate-$U$ regime. At sufficiently low $U$ there is a metallic
phase characterized by no self-consistent magnetic order and a gapless
spectrum.

Despite only considering quadratic fluctuations in the RPA scheme, we observe
two types of phase transition out of the ordered phases: Firstly, the magnon
modes come down to zero energy with a linear dispersion around some
incommensurate critical wavevector; secondly, the spin-wave velocity around
the ordering wavevector vanishes. We will discuss the physical implications of
these two transitions in the following. These two types of transitions can be
easily characterized, in contrast to the descending accidental zero mode at
$M^{\prime}$ discussed in the previous section, which does not show any
critical behaviour (remaining quadratic as it descends down in energy at the
RPA level).

The first type of transition occurs when the ordered phase is destabilized
once a magnon branch reaches zero energy at a wavevector incommensurate with
the ordering vector. In this regime, despite a self-consistent mean-field
solution, the Ansatz is not stable against fluctuations. We find that at the
vertical boundaries of the 120∘ phase, as drawn in the phase diagram Fig.
1(d), the transition is of this type whereby a transverse spin fluctuation
mode closes at the $M$ point (and equivalently the $Y_{1}$ point in the MBZ)
see Fig. 5(a,b). This was previously observed along the $t^{\prime}=0$ line
for varying $U$ [45], but we find the behaviour persists along the whole
extent of the phase boundary.

Numerically the boundary is ascertained by examining two features, which
happen simultaneously: the diverging matrix elements
$\chi(\mathbf{q}_{\mathrm{closing}},0)$ at zero energy which become larger
than the same elements at the ordering wavevector, as well as the lowest zero
eigenvalue of $[1-U\chi(\mathbf{q}_{\mathrm{closing}},\omega)^{(0)}]$ (the
mode energy) reaching zero $\omega_{\mathrm{closing}}=0$. Passing through the
phase boundary, these modes become linear and then square-root like, with no
pole in the imaginary part of the susceptibility matrix at the critical
closing point.

The RPA scheme can easily be extended to the $U\to\infty$ limit where it
predicts the direct transition between 120∘ and stripe ordered phases at
$t^{\prime}_{c}\approx 0.35$. At strong coupling the system also appears
symmetric in $t^{\prime}$, and from the 120∘ phase the closing wavevector is
at the $M$ point. In connection to the classical results for the
$J_{1}$–$J_{2}$ Heisenberg model, the exchanges to leading order in
perturbation theory are $J_{1}=4t^{2}/U$ and $J_{2}=4t^{\prime 2}/U$. Our
results predict in this limit $J_{2,c}\approx 0.123$, in good agreement with
the classical $J_{2,c}=1/8$ [56, 59].

From the stripe ordered phase however, the closing wavevector in the RPA
scheme is dependent on the position along the phase boundary; at high-$U$ the
closing wavevector is at the corner of the three-site 120∘ unit cell MBZ, see
Fig. 5(c,d). At reduced $U$ towards the bottom of the lobe of stripe order in
the phase diagram Fig. 1(d), the closing vector moves towards $Y$ and then to
$K$. This is in the region where the RPA scheme results in a negative-energy
accidental mode, and we expect that the position of this closing vector could
be moved in a self-consistent calculation which considers higher-order
fluctuations.

Figure 5: Critical modes closing for $U=30$ with (a,b) $t^{\prime}=0.358$ in
the 120∘ phase and (c,d) $t^{\prime}=0.365$ in the stripe ordered phase. The
figures (a,c) show the dynamical structure factor $S(\mathbf{q},\omega)$ along
the paths from $\Gamma$ to the respective ordering wavevector. The dashed
lines represent the zero of $[1-U\chi^{(0)}(\mathbf{q},\omega)]$. The static
structure factor $S(\mathbf{q})$ is plotted in (b,d). Linear size $L=720$
($N=N_{\mathrm{UC}}L^{2}$ sites) for spectral functions and $L=72$ for static
structure factor.

The second type of transition is triggered by long-wavelength spin
excitations. When the spin-wave velocity of excitations about the ordering
wavevector goes to zero, the quantum correction to the sublattice
magnetization are expected to diverge [59]. We find that this occurs at finite
$t^{\prime}\approx 0.2$–$0.4$ in the 120∘ ordered phase where it preempts the
closing of the mode at $M$. Therefore, the order is directly destroyed by
long-wavelength transverse fluctuations in this regime, marked as a red line
in Fig. 1(d). Fig. 6 shows the transverse susceptibility components, and
highlights the massless out-of-plane Goldstone excitations at this transition.

It is tempting to relate the second unusual transition to the recent numerical
works observing a chiral QSL for $t^{\prime}=0$ and small $U$ [49, 50, 51].
Indeed, the strong out-of plane fluctuations would be naturally related to an
instability towards a non-coplanar magnetically ordered state, which has been
proposed as the parent state of the chiral QSL [89, 90, 91] and the soft
structure factor at this point is also observed numerically in the CSL phase
[55]. Of course, we cannot make definite claims about the nature of the
disordered phases because our RPA scheme assumes the presence of a magnetic
state. Nevertheless, our results certainly highlight the different nature of
the order-disorder transition from the 120∘ phase along the bottom and sides
of the phase boundary. Along the sides there is a clear tendency to melting
because of short-wavelength fluctuations, and along the bottom there is a
regime where out-of-plane fluctuations destabilize the order. We also note
that the actual transition is expected to occur at a higher $U_{c}$ (and
possibly different $t^{\prime}$) because the harmonic fluctuations are further
renormalized by higher-order corrections.

Figure 6: Vanishing of the spin wave velocity along the red transition line of
the phase diagram. (a) The transverse susceptibility
$\operatorname{Im}\chi^{+-}(\mathbf{q},\omega)$. The dashed lines represent
the zero of $[1-U\chi^{(0)}(\mathbf{q},\omega)]$. (b) The static structure
factor $S(\mathbf{q})$. Calculations for $U=4.6$, $t^{\prime}=0.289$; linear
size $L=720$ ($N=3L^{2}$ sites) for dynamical susceptibility and $L=72$ for
static structure factor.

##  VI Conclusion

We have studied the spin fluctuations of the $t$–$t^{\prime}$–$U$ Hubbard
model on the triangular lattice and mapped out the phase diagram at half
filling. Employing a self-consistent mean-field plus RPA approximation allowed
us to calculate the spin structure factor in the thermodynamic limit.
Therefore, the method is ideally suited to provide predictions for inelastic
neutron scattering experiments complementary to the usual spin wave
calculations applicable only in the insulating large interaction limit.

Despite its limitations, we have shown that the method provides a
comprehensive picture of the spin fluctuations in the ordered stripe and 120∘
phases. Thereby, it allows the delineation of phase boundaries according to
the stability of magnon modes. We find that the phase diagram shown in Fig.
1(d) is in remarkable agreement with recent numerical works. In general, the
self consistent RPA theory is well suited to the study of itinerant magnetic
systems, including away from half-filling. In the future it would be
worthwhile to explore the entire doping dependence of the triangular lattice
Hubbard model.

We have provided a first study of spin excitations within the stripe
magnetically ordered phase. We showed how the inclusion of quantum
fluctuations removes an accidental zero mode from the underlying magnetic
degeneracy, which provides an extension of the order-by-disorder mechanism to
the finite-$U$ regime [56, 59]. This was done by modifying the bare
susceptibility to account for a fluctuation-renormalized self-energy [86, 87].
In the future it would be interesting to include higher-order corrections and
perform a fully self-consistent calculation. While technically and numerically
challenging, this would provide a more complete picture of the effect of
charge-fluctuations and the stability of magnetic phases for hole and electron
doping.

The phase boundaries in the large $U$ limit are in agreement with a classical
transition between ordered phases at $J_{2}/J_{1}=1/8$ [56, 59]. In the finite
$U$ regime we also find spin disordered regimes which have been at the centre
of recent numerical works both for the corresponding Heisenberg limit [43, 63,
61] and for the Hubbard model [49, 50, 51]. Given the difficulty of accessing
with purely numerical methods the dynamical behavior of two-dimensional
frustrated models with QSL ground states, a complementary approach is provided
by the self-consistent RPA method which has been extended to parton-mean-field
descriptions of Heisenberg models with fractionalized excitations [92, 93]. An
exciting direction for future research is to explore similar schemes for QSL
phases of Hubbard models in order to understand the nature of spin excitations
in the presence of spin fractionalization and charge fluctuations.

## Acknowledgements

We would like to thank Philipp Gegenwart, Markus Drescher and Frank Pollmann
for useful discussions. J.K. acknowledges support from the Imperial-TUM
flagship partnership. The research is part of the Munich Quantum Valley, which
is supported by the Bavarian state government with funds from the Hightech
Agenda Bayern Plus. H.-K.J. is funded by the European Research Council (ERC)
under the European Unions Horizon 2020 research and innovation program (Grant
agreement No. 771537).

Data and materials availability – Code and data are available on Zenodo [94].

## References

  * Hubbard [1963] J. Hubbard, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 276, 238 (1963).
  * Arovas  _et al._ [2022] D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, Annual review of condensed matter physics 13, 239 (2022).
  * Auerbach [2012] A. Auerbach,  _Interacting electrons and quantum magnetism_ (Springer Science & Business Media, 2012).
  * Lieb and Wu [1968] E. H. Lieb and F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968).
  * Giamarchi [2003] T. Giamarchi,  _Quantum physics in one dimension_ , Vol. 121 (Clarendon press, 2003).
  * Essler  _et al._ [2005] F. H. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin,  _The one-dimensional Hubbard model_ (Cambridge University Press, 2005).
  * Hirsch [1985] J. E. Hirsch, Physical Review B 31, 4403 (1985).
  * White  _et al._ [1989] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar, Phys. Rev. B 40, 506 (1989).
  * LeBlanc  _et al._ [2015] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 5, 041041 (2015).
  * Vannimenus and Toulouse [1977] J. Vannimenus and G. Toulouse, J. Phys. C 10, L537 (1977).
  * Toulouse  _et al._ [1987] G. Toulouse _et al._ , Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications 9, 99 (1987).
  * Diep [2004] H. T. Diep,  _Frustrated spin systems_ (World Scientific, 2004).
  * Lacroix  _et al._ [2011] C. Lacroix, P. Mendels, and F. Mila,  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_ , Vol. 164 (Springer, 2011).
  * Anderson [1973] P. W. Anderson, Mater. Res. Bull. 8, 153 (1973).
  * Anderson [1987] P. W. Anderson, Science 235, 1196 (1987).
  * Lee [2008] P. A. Lee, Science 321, 1306 (2008).
  * Balents [2010] L. Balents, Nature 464, 199 (2010).
  * Savary and Balents [2016] L. Savary and L. Balents, Rep. Prog. Phys. 80, 016502 (2016).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
  * Knolle and Moessner [2019] J. Knolle and R. Moessner, Annu. Rev. Condens. Matter Phys. 10, 451 (2019).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Science 367 (2020).
  * Li  _et al._ [2020] Y. Li, P. Gegenwart, and A. A. Tsirlin, Journal of Physics: Condensed Matter 32, 224004 (2020).
  * Kanoda and Kato [2011] K. Kanoda and R. Kato, Annual Review of Condensed Matter Physics 2, 167 (2011).
  * Powell and McKenzie [2011] B. J. Powell and R. H. McKenzie, Reports on Progress in Physics 74, 056501 (2011).
  * Li  _et al._ [2015] Y. Li, H. Liao, Z. Zhang, S. Li, F. Jin, L. Ling, L. Zhang, Y. Zou, L. Pi, Z. Yang, _et al._ , Scientific reports 5, 1 (2015).
  * Shen  _et al._ [2016] Y. Shen, Y.-D. Li, H. Wo, Y. Li, S. Shen, B. Pan, Q. Wang, H. Walker, P. Steffens, M. Boehm, _et al._ , Nature 540, 559 (2016).
  * Ding  _et al._ [2019] L. Ding, P. Manuel, S. Bachus, F. Grußler, P. Gegenwart, J. Singleton, R. D. Johnson, H. C. Walker, D. T. Adroja, A. D. Hillier, _et al._ , Physical Review B 100, 144432 (2019).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Wu  _et al._ [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. Lett. 121, 026402 (2018).
  * Wu  _et al._ [2019] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Phys. Rev. Lett. 122, 086402 (2019).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, _et al._ , Nature materials 19, 861 (2020).
  * Kennes  _et al._ [2021] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. Basov, A. N. Pasupathy, and A. Rubio, Nature Physics 17, 155 (2021).
  * Ni  _et al._ [2019] G. Ni, H. Wang, B.-Y. Jiang, L. Chen, Y. Du, Z. Sun, M. Goldflam, A. Frenzel, X. Xie, M. Fogler, _et al._ , Nature communications 10, 1 (2019).
  * Xian  _et al._ [2019] L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, Nano letters 19, 4934 (2019).
  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, _et al._ , Nature 579, 353 (2020).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, _et al._ , Nature 579, 359 (2020).
  * An  _et al._ [2020] L. An, X. Cai, D. Pei, M. Huang, Z. Wu, Z. Zhou, J. Lin, Z. Ying, Z. Ye, X. Feng, _et al._ , Nanoscale horizons 5, 1309 (2020).
  * Huse and Elser [1988] D. A. Huse and V. Elser, Phys. Rev. Lett. 60, 2531 (1988).
  * Bernu  _et al._ [1994] B. Bernu, P. Lecheminant, C. Lhuillier, and L. Pierre, Phys. Rev. B 50, 10048 (1994).
  * Capriotti  _et al._ [1999] L. Capriotti, A. E. Trumper, and S. Sorella, Phys. Rev. Lett. 82, 3899 (1999).
  * White and Chernyshev [2007] S. R. White and A. L. Chernyshev, Phys. Rev. Lett. 99, 127004 (2007).
  * Li  _et al._ [2022] Q. Li, H. Li, J. Zhao, H.-G. Luo, and Z. Y. Xie, Phys. Rev. B 105, 184418 (2022).
  * Singh [2005] A. Singh, Phys. Rev. B 71, 214406 (2005).
  * Sahebsara and Sénéchal [2008] P. Sahebsara and D. Sénéchal, Phys. Rev. Lett. 100, 136402 (2008).
  * Yoshioka  _et al._ [2009] T. Yoshioka, A. Koga, and N. Kawakami, Phys. Rev. Lett. 103, 036401 (2009).
  * Yang  _et al._ [2010] H.-Y. Yang, A. M. Läuchli, F. Mila, and K. P. Schmidt, Phys. Rev. Lett. 105, 267204 (2010).
  * Shirakawa  _et al._ [2017] T. Shirakawa, T. Tohyama, J. Kokalj, S. Sota, and S. Yunoki, Phys. Rev. B 96, 205130 (2017).
  * Szasz  _et al._ [2020] A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore, Phys. Rev. X 10, 021042 (2020).
  * Wietek  _et al._ [2021] A. Wietek, R. Rossi, F. Šimkovic, M. Klett, P. Hansmann, M. Ferrero, E. M. Stoudenmire, T. Schäfer, and A. Georges, Phys. Rev. X 11, 041013 (2021).
  * Szasz and Motruk [2021a] A. Szasz and J. Motruk, Phys. Rev. B 103, 235132 (2021a).
  * Zhu  _et al._ [2022] Z. Zhu, D. N. Sheng, and A. Vishwanath, Phys. Rev. B 105, 205110 (2022).
  * Chen  _et al._ [2022] B.-B. Chen, Z. Chen, S.-S. Gong, D. N. Sheng, W. Li, and A. Weichselbaum, Phys. Rev. B 106, 094420 (2022).
  * Cookmeyer  _et al._ [2021] T. Cookmeyer, J. Motruk, and J. E. Moore, Phys. Rev. Lett. 127, 087201 (2021).
  * Jolicoeur  _et al._ [1990] T. Jolicoeur, E. Dagotto, E. Gagliano, and S. Bacci, Phys. Rev. B 42, 4800 (1990).
  * Villain  _et al._ [1980] J. Villain, R. Bidaux, J.-P. Carton, and R. Conte, Journal de Physique 41, 1263 (1980).
  * Henley [1989] C. L. Henley, Physical review letters 62, 2056 (1989).
  * Chubukov and Jolicoeur [1992] A. V. Chubukov and T. Jolicoeur, Phys. Rev. B 46, 11137 (1992).
  * Kaneko  _et al._ [2014] R. Kaneko, S. Morita, and M. Imada, Journal of the Physical Society of Japan 83, 093707 (2014).
  * Iqbal  _et al._ [2016] Y. Iqbal, W.-J. Hu, R. Thomale, D. Poilblanc, and F. Becca, Phys. Rev. B 93, 144411 (2016).
  * Zhu and White [2015] Z. Zhu and S. R. White, Phys. Rev. B 92, 041105 (2015).
  * Hu  _et al._ [2015] W.-J. Hu, S.-S. Gong, W. Zhu, and D. N. Sheng, Phys. Rev. B 92, 140403 (2015).
  * Gong  _et al._ [2017a] S.-S. Gong, W. Zhu, J.-X. Zhu, D. N. Sheng, and K. Yang, Phys. Rev. B 96, 075116 (2017a).
  * Hu  _et al._ [2019] S. Hu, W. Zhu, S. Eggert, and Y.-C. He, Phys. Rev. Lett. 123, 207203 (2019).
  * Drescher  _et al._ [2022] M. Drescher, L. Vanderstraeten, R. Moessner, and F. Pollmann, Dynamical signatures of symmetry broken and liquid phases in an $s=1/2$ heisenberg antiferromagnet on the triangular lattice (2022), arXiv:2209.03344 [cond-mat.str-el] .
  * Misumi  _et al._ [2017] K. Misumi, T. Kaneko, and Y. Ohta, Phys. Rev. B 95, 075124 (2017).
  * Tocchio  _et al._ [2020] L. F. Tocchio, A. Montorsi, and F. Becca, Phys. Rev. B 102, 115150 (2020).
  * Morita  _et al._ [2002] H. Morita, S. Watanabe, and M. Imada, Journal of the Physical Society of Japan 71, 2109 (2002).
  * Watanabe  _et al._ [2008] T. Watanabe, H. Yokoyama, Y. Tanaka, and J. Inoue, Phys. Rev. B 77, 214505 (2008).
  * Tocchio  _et al._ [2009] L. F. Tocchio, A. Parola, C. Gros, and F. Becca, Phys. Rev. B 80, 064419 (2009).
  * Tocchio  _et al._ [2013] L. F. Tocchio, H. Feldner, F. Becca, R. Valentí, and C. Gros, Phys. Rev. B 87, 035143 (2013).
  * Szasz and Motruk [2021b] A. Szasz and J. Motruk, Phys. Rev. B 103, 235132 (2021b).
  * Sherman  _et al._ [2022] N. E. Sherman, M. Dupont, and J. E. Moore, Spectral function of the $j_{1}-j_{2}$ heisenberg model on the triangular lattice (2022), arXiv:2209.00739 [cond-mat.str-el] .
  * Schrieffer  _et al._ [1989] J. R. Schrieffer, X. G. Wen, and S. C. Zhang, Phys. Rev. B 39, 11663 (1989).
  * Chubukov and Frenkel [1992] A. V. Chubukov and D. M. Frenkel, Phys. Rev. B 46, 11884 (1992).
  * Singh and Tešanović [1990] A. Singh and Z. Tešanović, Phys. Rev. B 41, 614 (1990).
  * Peres and Araújo [2002] N. M. R. Peres and M. A. N. Araújo, Phys. Rev. B 65, 132404 (2002).
  * Knolle  _et al._ [2010] J. Knolle, I. Eremin, A. V. Chubukov, and R. Moessner, Phys. Rev. B 81, 140506 (2010).
  * Brydon and Timm [2009] P. M. R. Brydon and C. Timm, Phys. Rev. B 80, 174401 (2009).
  * Kaneshita and Tohyama [2010] E. Kaneshita and T. Tohyama, Phys. Rev. B 82, 094441 (2010).
  * Knolle  _et al._ [2011] J. Knolle, I. Eremin, and R. Moessner, Phys. Rev. B 83, 224503 (2011).
  * Jolicoeur and Le Guillou [1989] T. Jolicoeur and J. C. Le Guillou, Phys. Rev. B 40, 2727 (1989).
  * Lovesey [1984] S. Lovesey,  _Theory of Neutron Scattering from Condensed Matter_, International series of monographs on physics No. Bd. 2 (Clarendon Press, 1984).
  * Moriya [2012] T. Moriya,  _Spin fluctuations in itinerant electron magnetism_ , Vol. 56 (Springer Science & Business Media, 2012).
  * Singh [1991] A. Singh, Phys. Rev. B 43, 3617 (1991).
  * Edwards and Hertz [1973] D. Edwards and J. Hertz, Journal of Physics F: Metal Physics 3, 2191 (1973).
  * Ghosh and Singh [2008] S. Ghosh and A. Singh, Phys. Rev. B 77, 094430 (2008).
  * Hickey  _et al._ [2017] C. Hickey, L. Cincio, Z. Papić, and A. Paramekanti, Phys. Rev. B 96, 115115 (2017).
  * Saadatmand and McCulloch [2017] S. N. Saadatmand and I. P. McCulloch, Phys. Rev. B 96, 075117 (2017).
  * Gong  _et al._ [2017b] S.-S. Gong, W. Zhu, J.-X. Zhu, D. N. Sheng, and K. Yang, Phys. Rev. B 96, 075116 (2017b).
  * Ho  _et al._ [2001] C.-M. Ho, V. N. Muthukumar, M. Ogata, and P. W. Anderson, Phys. Rev. Lett. 86, 1626 (2001).
  * Zhang and Li [2020] C. Zhang and T. Li, Phys. Rev. B 102, 075108 (2020).
  * Willsher  _et al._ [2022] J. Willsher, H.-K. Jin, and J. Knolle, Magnetic excitations, phase diagram and order-by-disorder in the extended triangular-lattice Hubbard model (2022).

## Appendix

The stripe Hamiltonian is

$H_{\mathrm{stripe}}(\mathbf{k})=\begin{pmatrix}\zeta_{\mathbf{k}}\mathbb{I}-\sigma\cdot\Delta_{A}&\delta_{\mathbf{k}}\mathbb{I}\\\
\delta_{\mathbf{k}}^{*}\mathbb{I}&\zeta_{\mathbf{k}}\mathbb{I}-\sigma\cdot\Delta_{B}\end{pmatrix}.$
(16)

In order to define and work with this matrix as a function of momentum, we
take the inverse lattice as a basis of momentum $(k_{1},k_{2})$ such that if
$\mathbf{r}_{i}$ are the lattice unit vectors,
$\mathbf{k}\cdot\mathbf{r}_{i}=k_{i}$. Written as a function of $k_{1,2}$, the
$AB$ hopping element is explicitly given by

$\delta_{\mathbf{k}}=t\left[1+e^{-ik_{1}}+e^{-ik_{2}}+e^{-i(k_{1}+k_{2})}\right]+t^{\prime}\left[e^{ik_{1}}+e^{i(k_{1}-k_{2})}+e^{-i(2k_{1}+k_{2})}+e^{-2ik_{1}}\right]$
(17)

and the diagonal $AA/BB$ hopping is given by

$\zeta_{\mathbf{k}}=2t\cos(k_{1})+2t^{\prime}\cos(k_{2}).$ (18)

Now for the 120∘ phase, the $6\times 6$ Hamiltonian reads

$H_{120}(\mathbf{k})=\begin{pmatrix}\zeta_{\mathbf{k}}\mathbb{I}-\sigma\cdot\Delta_{A}&\delta^{AB}_{\mathbf{k}}\mathbb{I}&\delta^{AC}_{\mathbf{k}}\mathbb{I}\\\
\delta^{*\,AB}_{\mathbf{k}}\mathbb{I}&\zeta_{\mathbf{k}}\mathbb{I}-\sigma\cdot\Delta_{B}&\delta^{BC}_{\mathbf{k}}\mathbb{I}\\\
\delta^{*\,AC}_{\mathbf{k}}\mathbb{I}&\delta^{*\,BC}_{\mathbf{k}}\mathbb{I}&\zeta_{\mathbf{k}}\mathbb{I}-\sigma\cdot\Delta_{C}\end{pmatrix},$
(19)

with off-diagonal elements given by

$\delta^{AB}_{\mathbf{k}}=t\left[1+e^{-ik_{2}}+e^{i(k_{1}-k_{2})}\right],\quad\delta^{AC}_{\mathbf{k}}=t\left[1+e^{-ik_{2}}+e^{-ik_{1}}\right],\quad\delta^{BC}_{\mathbf{k}}=t\left[1+e^{-ik_{1}}+e^{i(k_{2}-k1)}\right].$
(20)

and diagonal elements given by

$\zeta_{\mathbf{k}}=2t^{\prime}\left[\cos(k_{1})+\cos(k_{2})+\cos(k_{1}-k_{2})\right].$
(21)