December 2020

Troubles with Global Monopoles in Quantum Gravity

Nemanja<EMAIL_ADDRESS>
QMAP, Department of Physics, University of California, Davis, CA 95616, USA

ABSTRACT

A theory with a global $O(3)$ symmetry broken at a scale $\eta$ admits
topological configurations: global monopoles realized by Goldstone fields
winding around the core of false vacuum. One may expect them to behave as
heavy, big composite objects, difficult to make and therefore mostly harmless.
However, after gravity is turned on, as long as Equivalence Principle holds,
one finds that global monopoles have a negative mass,
$M\sim-\eta/\sqrt{\lambda}$, where $\lambda$ is the field theory coupling.
This could catalyze an instability in the space of a global monopole, by the
production of additional global monopole-antimonopole pairs along with normal
particles, leading to energy production ab nihilo, if the pair energy is
dominated by their negative rest masses. In a theory with unbroken local
Poincaré symmetry, this could lead to a divergent ‘decay rate’ of the global
monopole configurations.

##  1 Introduction

It has been argued for a long time that quantum gravity does not permit exact
global symmetries (see, e.g. [1, 2, 3, 4, 5, 6, 7, 8]). An argument often
heard is that since black holes do not conserve global charges, and they count
among the Hilbert space states in quantum gravity, in processes where they
appear as virtual states the global charges cannot be conserved. Instead these
processes generate global symmetry violating operators in effective theory.
This argument is appealing, but it seems very difficult to implement it
directly, for example, as a diagrammatic calculation of a global symmetry
breaking operator generated by a loop including a virtual black hole.

More recently, different arguments emerged, invoking weak gravity conjecture
(WGC) [9] and AdS/CFT [10, 11]. These arguments rest on ideas on how the
theory of gravity may behave in the UV. In particular WGC is interesting since
it seems to block trying to realize global symmetries as a zero charge limit
of gauge symmetries. Yet while there is evidence supporting these ideas, we
still do not know if such specific frameworks for UV completions are unique or
generic. Hence it is interesting to pursue possible IR avenues of conflict
between global symmetries and gravity, to see if anything can be shaken out
from the mix of the two.

An example is provided by an interesting exploration of (in)compatibility of
non-compact, continuous shift symmetry with gravity in [12]. There, the
positivity of energy (combined with the implicit assumption that the geometry
is static) yields a bound on the total field variation exterior to a source.
This is despite the fact that in flat space the field variations can be
arbitrarily large. Further examples were provided in, for example, [13, 14,
15], where it was found that gravity appears to obstruct large field changes,
which indirectly implies that the naïve shift symmetry of the flat space limit
cannot coexist with gravity.

Here we will pursue a different approach, which will use only very basic IR
properties of gravity, employing very few ideas about the UV. We will consider
theories which have a spontaneously broken global $O(3)$ symmetry. Theories
with such structures seem legitimate in flat space without gravity. They have
a nontrivial vacuum manifold, with low energy Goldstone modes described by a
nonlinear $\sigma$-model. In addition to local excitations they also admit
topological global monopoles, which correspond to nontrivial Goldstone fields
winding around the region of false vacuum in the core. They exist when the
vacuum manifold is isomorphic to the sphere at infinity. In flat space, global
monopoles seem to be nonperturbative structures which are heavy, big composite
objects, that tend to be mostly harmless [16]. However even this is not
without controversy, since the issue of their stability has already been
debated vigorously [17, 18, 19, 20, 21, 22, 23, 24, 25]. We will briefly
review the status of this debate and outline the features relevant for our
argument. We will work with the symmetry breaking scale $\eta\ll M_{Pl}$, to
avoid any additional issues concerning large field variations111It is in fact
known that super-Planckian field values induce topological inflation [28, 29],
which in light of the results of [30] might even be viewed as a form of
instability [31].. Further we will work in weak coupling.

Turning gravity on lets global monopoles spring several surprises. First, in
the case of isolated global monopoles, their Goldstone fields support a solid
deficit angle [26]. This appears to regulate the flat space divergent energy
density distribution outside of a global monopole. Second, thanks to this,
global monopoles have a negative gravitational mass,
$M\sim-\eta/\sqrt{\lambda}$ [27], where $\lambda$ is the field theory
coupling. Finally, we will invoke a simple argument using Equivalence
Principle to infer that as a result the inertial mass of global monopoles is
also negative. Combining this with unbroken local Poincaré symmetry, and using
continuity of dispersion relations, implies that global monopoles should be
described in a local theory as ghosts, since they reside on mass shells which
are inverted relative to normal matter. To a good approximation they seem to
be ‘hard’ at low energies.

This immediately suggests possible catastrophic consequences, just like in the
case of theories with hard fundamental ghosts [32, 33, 34, 35]: if they are at
least metastable in local theory, and their rest energy is dominated by mass
terms, global monopole-antimonopole pairs may provide means to catalyze an
instability of the background spacetime of a global monopole, due to their
pair production along with normal positive energy particles. In this case, at
least, the positive energy theorem of [36, 37] is not directly applicable
since the background has a deficit angle. As long as the process amplitude is
not exactly zero, and the net energy of the pair is negative, the ‘decay rate’
of a global monopole would diverge when the theory has unbroken Poincaré
symmetry in the UV. This would be reminiscent to how the flat space stable
vacuum appears to reject remnants [38].

Hence we may infer that a consistent theory of quantum gravity with unbroken
local Poincaré symmetry in the UV may not support stable global monopoles
which arise in field theories with spontaneously broken global symmetries. If
this is true, our example, albeit somewhat special, suggests a glimpse into
what may go wrong with quantum gravity if we try to push global symmetries
onto it. An alternative, which would prevent this conclusion, is that the pair
energy is never negative despite the fact that their local rest masses are.
For this to be the case, a generalization of the theorems [36, 37] to some
spaces which are not asymptotically flat could be true, although there is no
proof at this time.

##  2 Global Monopole Solutions and Their Properties

Global monopoles are topological solitons of field theories with spontaneous
symmetry breaking, which contain an $O(3)$ symmetry, with a potential

$V(\vec{\Phi})=\frac{\lambda}{4}\,\Bigl{(}\vec{\Phi}^{2}-\eta^{2}\Bigr{)}^{2}\,.$
(1)

Because any vacuum at infinity spontaneously breaks $O(3)$ to $O(2)$, the
theory has massless Goldstone bosons. They are crucial for the existence of a
global monopole because they give long range fields which wind around the
global monopole core, providing the topological charge distinguishing a global
monopole from the vacuum. In such theories there are configurations where the
field in the target space develops a hedgehog configuration around the core of
the false vacuum,

$\vec{\Phi}=f(r)\,\frac{\vec{r}}{r}\,,$ (2)

wrapping around it (in the example above, exactly once) as it interpolates to
the degenerate manifold of true vacua at infinity (for a review see [16]). As
the Goldstones range out to infinity in a global monopole configuration,
interpolating between the different true vacua far away, they yield stress
energy tensor which vanishes only as [16]

$T^{0}{}_{0}=-T^{r}{}_{r}\sim\frac{\eta^{2}}{r^{2}}\,,$ (3)

while the angular components vanish much faster. As a result, individual
global monopoles formally seem to have a divergent field energy,

$E=4\pi\int^{R}dr\,r^{2}\,T^{0}{}_{0}\sim 4\pi\eta^{2}R\,.$ (4)

Hence individual global monopoles cannot be produced in flat space unless
infinite energy initial conditions are used. Thus one may imagine that due to
such huge energy in the global monopole field, its mass, and therefore its
inertia, should be very large. However this single source solution can be
regulated by adding a global antimonopole, which is a parity image of a global
monopole, at a distance $\sim R$ from the global monopole [16]. In that case
the energy (4) will be cut off at infinity, since the antimonopole serves as a
sink of field lines, at a finite distance from the source. This yields a
faster dilution of stress energy far away.

Early on it has been asserted that a global monopole may be unstable [17] (see
also [18, 19]) since it can unwind by itself, thanks to an angular
perturbation of the spherical isopotential shells around it. A small
perturbation of the global monopole on, say, the North Pole has a tendency to
grow in the axial direction, and pull in the Goldstone field lines towards the
vertical axis, where they can unwind.

A more detailed analysis has shown that this instability really induces a
force acting on the global monopole, accelerating it in the direction of the
perturbation [18, 20, 21, 22]. We could imagine this configuration of a
“leashed” global monopole as a global monopole-antimonopole pair separated by
a large distance. The array breaks spherical symmetry, since the presence of a
global antimonopole introduces a direction, and so the field lines tend to
align with it and form a global ‘flux tube’ [20, 21]. The initial realignment
starts slowly due to a large separation, but then the constant force pulls the
two together [22].

Thus we assert that the instability observed in [17] can be viewed as a signal
of the onset of the process of global monopole-antimonopole annihilation in
the large initial distance limit [18, 20, 22]. Moreover, although the
topological charge is accounted for by the boundary conditions at infinity,
the energy barrier between configurations with a different charge seems to be
finite [23]. This suggests that the number of global monopoles could be
‘rearranged’ with a finite energy investment [23] despite the apparent
divergence of field energy (4).

More careful examinations with gravity turned on reveal interesting surprises.
To see how the energy (4) behaves, we may use General Relativity as a probe.
It turns out that the leading order sources in $T^{\mu}{}_{\nu}$ do not yield
any local gravitational forces. Instead, the leading order $T^{\mu}{}_{\nu}$
which had yielded the divergence (4) in flat space now sources a global solid
deficit angle of the space exterior to a global monopole [26]. The metric of a
global monopole outside of the core is

$ds^{2}=-\Bigl{(}1-8\pi
G_{N}\eta^{2}-\frac{2G_{N}M}{r}\Bigr{)}\,dt^{2}+\Bigl{(}1-8\pi
G_{N}\eta^{2}-\frac{2G_{N}M}{r}\Bigr{)}^{-1}dr^{2}+r^{2}d\Omega_{2}\,.$ (5)

The solid angle surrounding a global monopole at infinity is
$\Delta\phi=4\pi(1-8\pi G_{N}\eta^{2})$, which follows from the leading order
behavior of the metric as $r\rightarrow\infty$. Thus the leading order energy
is effectively degravitated. Locally, all it does is enhance the value of
Newton’s constant in the space outside of a global monopole. This is obvious
from Gauss’ law at infinity, where the angular deficit implies that there is
stronger fields for a fixed source mass since there is less spatial volume for
the field lines to dilute [39]. The long range effects of $T^{\mu}{}_{\nu}$
persist only as the deficit angle [26, 27]. They are important for the
understanding of a global monopole-antimonopole pair creation as we will argue
later; they are crucial for understanding their local gravitational and
inertial properties.

At the subleading level, gravity reveals an even greater surprise [27]: a
global monopole has a negative gravitational mass! Technically, this rather
odd feature can be understood as follows. Inside a global monopole core the
field configuration is near the false vacuum $\Phi=0$, where the energy
density is $\rho\sim V\sim\lambda\eta^{4}/4$, which behaves like a positive
cosmological constant, exhibiting repulsion of normal matter due to the
negative pressure $p\sim-\rho$, which induce forces $\propto(\rho+3p)<0$. The
interior geometry is approximately given by de Sitter static patch metric,

$ds^{2}=-\Bigl{(}1-H_{core}^{2}r^{2}\Bigr{)}\,dt^{2}+\Bigl{(}1-H_{core}^{2}r^{2}\Bigr{)}^{-1}dr^{2}+r^{2}d\Omega_{2}\,,$
(6)

where

$H_{core}^{2}\sim\frac{8\pi G_{N}V}{3}\sim\frac{2\pi
G_{N}\lambda\eta^{4}}{3}\,.$ (7)

Because the exterior metric has a deficit angle, the interior and exterior can
be smoothly matched at the “surface” of the global monopole, where the metrics
and their derivatives continuously connect. This yields [27], up to
corrections $\propto 8\pi G_{N}\eta^{2}\ll 1$,

$r_{0}\simeq\frac{2}{\sqrt{\lambda}\eta}\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}M\simeq-\frac{16\pi\eta}{3\sqrt{\lambda}}\,.$
(8)

The mass in the exterior metric must be negative by Birkhoff’s theorem, since
the interior mass [27, 40] is $M_{core}\sim(\rho_{core}+3p_{core})\,{\rm
Vol}_{core}\sim-2\rho_{core}\,{\rm Vol}_{core}$. Importantly, also note that
as long as

$\frac{2G_{N}|M|}{r_{0}}\sim\frac{2}{3}8\pi
G_{N}\eta^{2}~{}\mbox{\raisebox{-2.58334pt}{$\stackrel{{\scriptstyle<}}{{\sim}}$}}~{}1\,,$
(9)

or so as long as $8\pi
G_{N}\eta^{2}~{}\mbox{\raisebox{-2.58334pt}{$\stackrel{{\scriptstyle<}}{{\sim}}$}}~{}1$,
the description of a monopole using flat space physics should be reliable,
since gravity is too weak in its core. Clearly, we are in this limit when
$8\pi G_{N}\eta^{2}\ll 1$. Further note that by combining Eqs. (8),

$r_{0}|M|=\frac{r_{0}}{\lambda_{compton}}\simeq\frac{32\pi}{3\lambda}\gg 1\,,$
(10)

which shows that in perturbation theory, $8\pi
G_{N}\eta^{2}~{}\mbox{\raisebox{-2.58334pt}{$\stackrel{{\scriptstyle<}}{{\sim}}$}}~{}1$,
$\lambda\ll 1$, the global monopole core size is much larger than its Compton
wavelength. So, global monopoles are classical, extended objects222As a
technical aside, to fully describe them we need to keep the radial mode in the
theory. Since the mass of the radial mode of $\vec{\Phi}$ is
$m_{\Phi}\sim\sqrt{\lambda}\eta$, we find $r_{0}m_{\Phi}={\cal O}(1)$,
implying that global monopoles are about the same size as the Compton
wavelength of the radial scalar modes, integrated out in the true vacuum.
Hence global monopoles explore the physics above the cutoff of the
perturbative $\sigma$-model describing only the Goldstones, since the cutoff
of the low energy effective theory is set by $m_{\Phi}$., which as long as
$8\pi
G_{N}\eta^{2}~{}\mbox{\raisebox{-2.58334pt}{$\stackrel{{\scriptstyle<}}{{\sim}}$}}~{}1$,
$\lambda\ll 1$ have gravitational properties that can be consistently
described using linearized General Relativity.

But this is not all: as long as gravity is governed by General Relativity,
which obeys Equivalence Principle, the global monopole inertial mass
$M_{inertial}$ must also be negative since it must be equal to its
gravitational mass $M$. To see this, imagine a probe particle near a global
monopole, in free fall in the geometry with the metric (5), following the
geodesic equation

$\ddot{x}^{\mu}+\Gamma^{\mu}_{\nu\lambda}\dot{x}^{\nu}\dot{x}^{\lambda}=0\,.$
(11)

Let the particle stand to the right of the monopole, moving slowly and being
separated a distance $r$. Because the global monopole gravitational mass $M$
is negative, the particle experiences a repulsive acceleration controlled by
the global monopole’s negative gravitational mass, as the particle mass
cancels out of the geodesic equation. So the force pushing the particle away
from the global monopole is $\vec{F}\sim-
mM\vec{r}_{0}/r^{2}=m|M|\vec{r}_{0}/r^{2}$.

Now: change coordinates from the rest frame of the global monopole to the rest
frame of the particle, and look at the motion of the global monopole governed
by the gravitational field of the particle, with a Schwarzschild metric of the
mass $m$. To the leading order, the global monopole follows geodesic equation
in the new coordinate system, experiencing an attractive acceleration induced
by the positive particle mass – since this time the global monopole’s inertial
mass cancels from the geodesic equation, the gravitational field and the
acceleration imparted by it is the same for all probes, controlled only by the
particle’s positive mass $m$, as follows from Equivalence principle. The force
pulling the global monopole towards the particle (on the right of the
monopole) is $\hat{F}\sim mM_{inertial}\,\vec{r}_{0}/r^{2}$. Again by
Equivalence Principle – since the motion of the particle in the global
monopole field must be the same as the motion of the global monopole in the
particle’s field – the forces $F$ and $\hat{F}$ must be equal and opposite
[41, 42], and hence $M_{inertial}=-|M|=M\equiv M_{gravitational}$, as claimed.
Therefore

$M_{inertial}\equiv M\simeq-\frac{16\pi\eta}{3\sqrt{\lambda}}\,.$ (12)

The calculation of the mass $M$ can be made more precise by generalizing the
notion of the ADM mass in the presence of the deficit angle, which changes the
boundary conditions at infinity. The generalization of the ADM mass for
spacetimes with solid deficit angle [43, 44] confirms the intuitive results
after the renormalization of $G_{N}$ by the deficit angle. This leaves us
precisely with the result of (8) to the leading order in the $8\pi
G_{N}\eta^{2}$ expansion [27, 43, 44], allowing the mass of a global monopole
to be negative.

The price, of course, is the deficit angle, amounting to the alteration of the
boundary conditions at infinity. In flat space, the net global monopole number
– ie the total winding number – enumerates superselection sectors. To create a
single global monopole would presumably be obstructed by the need to ‘bring’
the charge from infinity, and so spend an infinite amount of energy, since
Goldstone fields dilute too slowly far away. However, these superselection
sectors seem to only be separated by finite energy barriers [23]. This, at
least, suggests that a global monopole-antimonopole pair can be formed
readily. In flat space without gravity, if the initial state is the Minkowski
vacuum, or even a configuration with a global monopole on shell, a global
monopole-antimonopole pair could be created quantum-mechanically if the
conservation laws permit it. In this limit all the energies are dominated by
positive rest masses, by the theorems [36, 37] in asymptotically flat spaces.
So the pair will quickly fall in and annihilate – with a zero net energy
release – existing a mere moment only because of the uncertainty principle. If
so, the flat space vacuum would be stable, since these events would merely be
one of the many quantum fluctuations that come and go. The diagrams
representing their contributions to any observables would be disconnected
bubbles, which factor out. In this argument the gravitational effects are
completely ignored.

Turning gravity on brings out the negative global monopole masses and deficit
angles. Having objects with negative masses, with all local energy densities
being positive, is a very unsettling proposition. Already [27] note that if
these configurations really existed, they could lead to the formation of Bondi
dipoles [45]. After noting this option, [27] step back and dismiss it on the
grounds that global monopoles are topological objects which are not ‘isolated
point particles, but extended sources.’ A single global monopole is indeed an
extended object in the strict geometric sense: it affects the boundary
conditions at infinity. Yet at high energies a global monopole can unwind as
the symmetry is restored. The charge and the field energy would then flow away
to infinity.

On the other hand, global monopole-antimonopole pairs extend far less, with
their long distance effects diluting away more quickly. This naturally
regulates the energy divergences far away. At distances much larger than the
global monopole-antimonopole separation, the net charge of the pair will look
to be zero, which means that the field strengths of the Goldstones must
diminish faster than in the case of a single global monopole. Therefore in a
quantum mechanical process global monopoles could be created as pairs and so a
net global monopole number should not even change333Such events could be
extremely rare, but if the rates are not zero the problems will arise..
However the numbers of global monopoles and global antimonopoles may
fluctuate.

This should survive the onset of gravity; in pair creation, the curvature
deformation mimicking deficit angle may be induced locally near individual
sources after the pair emerges, but will not extend to infinity as long as the
field energy between the pair is finite. This local distortion of the geometry
may subsequently propagate farther – or not – depending on the fate of the
pair. However in practice, in the linearized gravity regime applicable to our
analysis due to choosing $8\pi G_{N}\eta^{2}\ll 1$, $\lambda\ll 1$, inside the
region of space where the pair forms, and close to the global monopoles, they
will behave as negative mass objects when probed by normal matter for all
practical intents and purposes.

We contend that this may allow the global monopole-antimonopole pair creation
to drive an instability of a configuration initially containing a global
monopole on shell through particle production. As we stressed above, the
formation of a global monopole-antimonopole pair out of the flat space vacuum
is not prohibited a priori if all conservation laws are satisfied. Since
Goldstone fields are cut off, this is a finite energy configuration. As we
stated above, this configuration should be interpreted as a pair of global
monopoles connected with a ‘flux tube’. This is motivated by the studies of
Goldhaber’s instability [17, 18], whereby a global monopole unwinds
directionally before annihilating with a global antimonopole [18, 20, 21, 22].
Since the global antimonopole is needed to regulate the Goldstone field energy
(4) anyway, it is natural to interpret Goldhaber’s instability as the tendency
of Goldstone fields to collapse to a flux tube [20, 21]. We maximize its
energy by assuming no energy loss in the field rearrangement, so that the
interaction energy mediated by the flux tube is given by the regulated
expression (4). This way we reproduce the formula for the constant force
acting on a global monopole (as already noted in [21]).

We can then view the global monopoles at the end points of the tube as locally
negative mass objects, when we get near them (see also [46] for a similar
point of view). The question then is, what is the energy of the pair, and
could it be negative. This question is subtle444We thank A. Vilenkin for
useful communications about this issue.. Although it may appear that we could
make the energy negative for a pair even in vacuum, when they are initially
sufficiently close, so the flux tube between them is short, and they are
moving, so their kinetic energy may add to their negative rest masses, the
overall global monopole-antimonopole resonance energy in asymptotically flat
space is prohibited from being negative by the theorems of [36, 37]. This
would block the instability in the vacuum.

However if the background already contains global monopoles, and therefore has
a deficit angle, the case is not so clear cut. The theorems [36, 37] do not
automatically guarantee that the pair energy must be positive. The total
energy of the system might not be positive definite due to different boundary
conditions at infinity. In this case, if a global monopole-antimonopole pair
pops out in a background already containing a global monopole, with an overall
negative energy, it must be accompanied by at least one positive mass
particle-antiparticle pair, such that the total system satisfies energy-
momentum conservation enforced by local Poincaré invariance – that the net
energy change is zero. Hence the production rate may be kinematically allowed.
Once the positive mass particles get on shell, they will fly away, since they
can be light – or even massless, being photons or even gravitons themselves.

In many cases the global monopole-antimonopole pairs which are created will
fall in and unwind away, absorbing some of the positive mass particles
floating aroud; in some, they may even end up wiggling around for a while and
increasing their (negative) kinetic energy by emitting positive mass
particles. But their ultimate fate past the fact that they could be created
seems largely irrelevant for the instability. If they are on shell, the
positive mass particles will escape regardless of the fate of the global
monopole-antimonopole pair, and stress-energy may flow out of the space which
is asymptotically locally the same as the vacuum, except for the deficit
angle. When local Poincaré symmetry is unbroken in the UV, the phase space for
the outgoing decay products has divergent volume just like similar processes
for the production of hard elementary ghosts, where the energy-momentum
conservation holds [32, 33, 34]. As a result, as long as the decay rate is not
exactly zero, it should diverge. The suppression due to the compositeness of
the global monopoles may make their pair production rate per unit phase space
volume exponentially small, but as long as it is not zero the disaster is
lurking. All it takes is a pair with a net negative energy and local Poincaré
invariance in the UV.

An alternative outcome is that even though the space has a deficit angle, and
so is not asymptotically flat, the results of [36, 37] still hold and exclude
negative energy pairs. If this is true then our discussion below may be viewed
as an explanation why [36, 37] need to extend to this case, despite the global
redefinition of the ADM mass [43, 44]. In this case reconciling the negative
mass of an isolated global monopole with the positive energy of a pair would
be quite interesting to see.

##  3 The Global Monopoles’ Fizz

Now we turn to the details of how global monopoles could trigger the
instability in the space surrounding them. As we noted in the beginning, we
will be working with safely subcritical global monopoles, where $8\pi
G_{N}\eta^{2}\ll 1$, $\lambda\ll 1$ so that we can for the most part ignore
the deficit angle. Locally, it only renormalizes $G_{N}$ in this limit;
however it is crucial for the existence of negative masses.

Consider now a global monopole-antimonopole pair, separated by a distance $R$,
in a locally asymptotically flat space of a distant global monopole. The
global monopole far away sets up the deficit angle at infinity. In the least
favorable case to our argument, when the force is attractive, since it is
constant [16, 22], the time it takes for the pair to annihilate is controlled
by the initial separation,

$\Delta T\sim\sqrt{\frac{R|M|}{4\pi\eta^{2}}}\,.$ (13)

At distances $r>r_{0}$ we will treat individual global monopoles comprising
the pair as hard particles, with negative masses $M$ of Eq. (12). To the
leading order, when they are close enough their total energy may be dominated
by the kinetic terms due to negative masses. The negativity of the mass
combined with unbroken local Poincaré symmetry implies that global monopoles
would behave as ghosts, with negative energy

$E=-\sqrt{\vec{p}^{2}+M^{2}}\,,$ (14)

for a given momentum $\vec{p}$. To see this, imagine an observer boosted
relative to a global monopole at rest. In the rest frame of the moving
observer, it appears that the global monopole is moving in the opposite
direction, with the relativistic energy $E\neq M$. By continuity of boosts, we
can dial down the relative velocity $\vec{v}$ to zero smoothly, which means
that in the limit $\vec{v}=0$ we must reproduce $E=-M$. Therefore
$E=-\sqrt{\vec{p}^{2}+M^{2}}$. This means that the mass shells of the global
monopoles are inverted relative to normal particles.

Next, we estimate how far the monopoles can be so that the net energy is
negative, using the energy balance formula based on (4). Note that if the
initial state was empty Minkowski space the net energy would have to be non-
negative by [36, 37]. However if there is a deficit angle due to a global
monopole far away, this need not be true. Note also that the formula (4) is at
best only an approximation in a conical space, but we will nevertheless use it
to illustrate the issues. So for a global monopole-antimonopole at rest, with
negative masses (8), to the leading order we take the excess energy over the
energy of the background monopole far away to be roughy given by the flat
space formula for the flux tube energy, but with negative monopole masses,

$E_{total}\simeq
4\pi\eta^{2}R+2M=4\pi\eta^{2}R-\frac{32\pi\eta}{3\sqrt{\lambda}}\,,$ (15)

and allowing for it to be negative if the monopole-antimonopole separation $R$
is small enough,

$R_{cr}~{}\mbox{\raisebox{-2.58334pt}{$\stackrel{{\scriptstyle<}}{{\sim}}$}}~{}\frac{8}{3\sqrt{\lambda}\eta}\simeq\frac{4}{3}r_{0}\,.$
(16)

The distance $R_{cr}$ is just-so larger than the global monopole core size,
meaning the global monopole-antimonopole pair are in contact. They will not
unwind right away, but may do so, at a rate depending on $r_{0}$. Still, even
this may suffice to trigger a cataclysm. Note that as we take
$\lambda\rightarrow 0$ the danger from long lived negative energy pairs seems
to grow.

Moreover, the monopoles are initially moving, having additional kinetic energy
so that $-M\rightarrow-M/\sqrt{1-v^{2}}$, so the critical separation might be
larger, and the net energy of the pair might still be negative. This opens a
window for a global monopole-antimonople separation to be at least a few times
larger than $r_{0}$, improving the approximation where the monopoles are
treated as hard ghosts. Thus when the pair are moving initially, by causality
they will survive for at least a time $\Delta
T\sim({R_{cr}}/{\sqrt{\lambda}\eta})^{1/2}>{\rm few}\times r_{0}$. Such global
monopole-antimonopole pairs would be approximately on shell: quantum
production of photons will be faster than global monopole-antimonopole
annihilation. Hence these global monopole-antimonopole pairs ought to be
included among those states which the space that contains a global monopole
can decay into, at least as metastable resonances.

Since in locally Poincaré invariant theories we must enforce energy
conservation, additional positive energy particles would be produced, such
that the total energy change is zero [32, 33]. A simple cartoon for a process
describing this is a spontaneous creation of a bound system of scalar quanta
forming the global monopole-antimonople pair along with a swarm of free
streaming ‘photons’ mediated by gravitational interactions (or any other
particle-antiparticle pair, due to the universality of gravity, including
gravitons or even the $O(3)$ Goldstones). The global monopole setting the
background can be taken to be far enough that its only effect on local physics
is via the deficit angle, which invalidates positive energy condition.
Otherwise the dynamics is as in any locally flat space. Treating the
background as a ‘vacuum’, the total process would involve many different
internal scalar excitations, whose number would be at least as large as
$M/m_{\Phi}\sim Mr_{0}\gg 1$, with an arbitrary number of external photons
popping out, as long as the net energy is zero. A channel with a given number
of scalars and two photons is sketched in Fig. 1.

Figure 1: One of many channels for producing a global monopole-antimonopole
pair + two photons, mediated by gravity, far from an initial global monopole
on shell. Since global monopoles are composite, at least $\sim 2n\sim 2Mr_{0}$
scalar ‘quanta’ are required to make up the pair.

Clearly, there is an exponentially large number of such channels. The total
amplitude of the process is obtained by summing amplitudes for each channel. A
large number of channels means that they would sum up mostly out of phase,
yielding destructive interference and suppressing the amplitude. The rate
would be effectively suppressed by compositeness of the global monopole-
antimonopole pair. However even an exponential suppression $\sim\exp(-S)$ is
not enough to stop the disaster, as long as the pair energy could be negative.

Indeed, this could trigger a catastrophe because the Lorentz-invariant phase
space for such processes is infinite. This is regardless of whether global
monopole-antimonopole pair end up going to infinity, or not. The important
point is, there is no stopping the photons, which would carry energy away with
them. If realized, this leads to asymptotic space instability, resulting in
energy production ab nihilo.

Although the initial states which contain a global monopole and deficit angle
break the local Poincaré spontaneously, this breaking is in the IR and does
not help with UV divergences. The reason that the phase space diverges for a
process like

$|0_{asymptotic}\rangle\rightarrow|M\bar{M}\,\gamma\gamma\ldots\,\rangle$

is that the global monopole masses are negative, due to deficit angle. The
divergence can be readily inferred from the transition amplitude following the
analyses in [32, 33, 34, 35]. Let’s consider the simplest divergent
contribution coming from the process $0_{asymptotic}\rightarrow
M\bar{M}\,\gamma\gamma$ whose one channel is depicted in Fig. 1. The
contribution to the decay rate of the asymptotically flat space far from the
global monopole source via the production of a global monopole-antimonopole
pair plus two photons is

$\displaystyle\Gamma$ $\displaystyle\ni$ $\displaystyle\int
d^{4}p_{1}d^{4}p_{2}d^{4}k_{1}d^{4}k_{2}\,|A|^{2}\,\delta(p_{1}^{2}-M^{2})\delta(p_{2}^{2}-M^{2})\delta(k_{1}^{2})\delta(k_{2}^{2})\times~{}~{}~{}~{}$
(17)
$\displaystyle~{}~{}\times\,\theta(-p_{10})\theta(-p_{20})\theta(k_{10})\theta(k_{20})\delta^{4}(p_{1}+p_{2}+k_{1}+k_{2})\,,$

Here $A$ is the transition amplitude for the process $0\rightarrow
M\bar{M}\,\gamma\gamma$, defined in the standard way by $\langle
p_{1},p_{2},k_{1},k_{2}|0_{asymptotic}\rangle={A}(p_{1},p_{2},k_{1},k_{2})\
\delta^{(4)}(p_{1}+p_{2}+k_{1}+k_{2})$ (thanks to the local Poincaré
symmetry), $p_{j}$ are the global monopole $4$-momenta and $k_{j}$ are the
photon $4$-momenta, and the signs in the step functions
$\theta(-p_{j0}),\theta(k_{j0})$ take care of the mutual reversal of the mass
shells, accounting for the negative inertial mass (12).

To proceed, it is convenient [34] to rewrite the integral in (17) by
introducing an auxiliary $4$-momentum $P$, and integrate over it to enforce
$P=p_{1}+p_{2}=-k_{1}-k_{2}$, which executes $4$-momentum conservation in two
steps. This yields [34]

$\Gamma\ni\int d^{4}P\,\gamma(P)\,,$ (18)

where

$\displaystyle\gamma(P)$ $\displaystyle=$ $\displaystyle\int
d^{4}p_{1}d^{4}p_{2}d^{4}k_{1}d^{4}k_{2}\,|A|^{2}\,\delta(p_{1}^{2}-M^{2})\delta(p_{2}^{2}-M^{2})\delta(k_{1}^{2})\delta(k_{2}^{2})\times~{}~{}~{}~{}$
(19)
$\displaystyle~{}~{}\times\,\theta(-p_{10})\theta(-p_{20})\theta(k_{10})\theta(k_{20})\delta^{(4)}(P+k_{1}+k_{2})\
\delta^{(4)}(P-p_{1}-p_{2})\,.$

In this integral, the auxiliary variable $P$ appears as an external momentum.
The local Poincaré symmetry then asserts that – since there are no preferred
directions in the UV – $\gamma$ can depend only on the magnitude
$P^{\mu}P_{\mu}=-s$. Note that this is a consequence of the Lagrangian of the
theory. Then changing the integration measure in (18) to $s$ and
$\vec{v}=\vec{P}/\sqrt{s}$ yields

$\Gamma\ni\int\frac{1}{2}\frac{d^{3}\vec{v}}{\sqrt{1+\vec{v}^{2}}}\int ds\ s\
\gamma(s)\,.$ (20)

This change of variables is analogous to spherical polar coordinates in $4D$.
However the difference is that $P$ is a vector in a $4D$ Minkowski space, and
thus in the decomposition $P^{\mu}P_{\mu}=P^{2}(1-\vec{v}^{2})=-s$ the
manifold of $\vec{v}$ is not compact. Indeed, using
$d^{4}P=dPP^{3}dH_{3}=dss\,dH_{3}/2$, and realizing that $dH_{3}$ is the
measure on the mass shells defined by $v_{0}^{2}=1+\vec{v}^{2}$, we see that
$\int dH_{3}/2=\frac{1}{2}\int
dv_{0}\,d^{3}\vec{v}\,\delta(v_{0}^{2}-1-\vec{v}^{2})=\int\frac{1}{2}\frac{d^{3}\vec{v}}{\sqrt{1+\vec{v}^{2}}}\sim
2\pi\int dv\,v$. If this were an integration over euclidean vectors, it would
correspond to integrating over a hemisphere from the pole to the equator,
which would give a finite answer since the sphere is compact555In that case,
the sign under the square root in the denominator of (20) would be negative..
However, the Lorentzian signature implies that the integration is over an
infinite range of the lengths of $\vec{v}$. Thus as long as $|A|^{2}\neq 0$
the integral diverges.

We can understand this by considering $4$-momentum conservation
$p_{1}+p_{2}+k_{1}+k_{2}=0$. Since the photons are null,
$k_{1}^{2}=k_{2}^{2}=0$. For the global monopoles, we have
$p_{1}^{2}=p_{2}^{2}=M^{2}$. Rewriting the conservation equation as the
expression for $p_{1}$ in terms of the rest then yields, after a simple
manipulation, $p_{2}\cdot k_{1}+p_{2}\cdot k_{2}+k_{1}\cdot k_{2}=0$. Next, we
can perform a boost to the rest frame of $p_{2}$, where $p_{2}=(M,0)$. In this
frame $k_{j}\cdot p_{2}=ME_{j}$. Using the projection constraint and defining
the angle $\theta$ between the two photons via
$\vec{k}_{1}\cdot\vec{k}_{2}=E_{1}E_{2}\cos(\theta)$, we finally obtain the
equation

$\cos(\theta)=1+M(\frac{1}{E_{1}}+\frac{1}{E_{2}})\,.$ (21)

Since the photons are on the positive mass shell, the term in the parenthesis
is positive definite. Thus if $M>0$, there would be no solutions for $\theta$,
since the RHS would exceed unity. However since the global monopole mass is
negative, $M<0$, the RHS will be smaller than unity at energies
$~{}\mbox{\raisebox{-2.58334pt}{$\stackrel{{\scriptstyle>}}{{\sim}}$}}~{}{\rm
few}\times|M|$, and the kinematical constraint (21) will be satisfied for some
scattering angle $\theta$. This will happen not only in the rest frame of
$p_{2}$, but in any frame boosted relative to it, with the same value of
$s=-(p_{1}+p_{2})^{2}$, which must also be included among the final decay
states. Since there is infinitely many such frames the rate will diverge. This
last step is precisely the integration over $\vec{v}$ above.

As long as we maintain that the local Poincaré symmetry is unbroken, the only
way to reverse this conclusion seems to be to have $|A|^{2}\equiv 0$, or to
block the global monopole-antimonopole pair from negative values. But this is
precisely the point we want to make: as long as a global monopole-antimonopole
pair with negative energy can be formed in the background with an initial
global monopole, however unlikely this outcome might be, the decay rate will
blow up due to unbroken local Poincaré invariance. It appears that one way to
accomplish this switch-off of $A$ is to decouple gravity completely, sending
$M_{Pl}\rightarrow\infty$. The world without gravity is not what we are after,
however.

We would be remiss to ignore the possibility that the positive energy theorems
from asymptotically flat space [36, 37] somehow extend to the conical space of
a global monopole. This could preclude the situation that the global monopole
rest masses ever dominate the flux tube energy. We cannot exclude this,
although it is difficult to see how this could occur666Note that if we could
take the limit $\lambda\rightarrow 0$, we might think the negative mass
monopoles would decouple (they would in flat space without gravity). However
this would enhance the global symmetry from $O(3)$ to the full symmetry group
of the $3D$ Euclidean space by adding three commuting continuous translation
shift symmetries. From the point of view of quantum gravity this would seem to
be an even worse outcome. for all values of $\lambda$. Further if this is the
case, then reconciling such a statement with the negative mass of an isolated
global monopole which follows from the redefinition of the ADM mass to conical
spaces [43, 44] would be quite interesting to see.

Alternatively, if the theory has a UV completion which is not locally Poincaré
invariant, the instability rate can be cut off and global monopoles may
coexist with gravity. An example might be global monopole-like configurations
in ‘emergent gravity’ models involving superfluid Helium [47]. Since the UV
completion is just the standard nonrelativistic quantum mechanics, the absence
of local Poincaré invariance in the UV introduces a cutoff to the integrals
contributing to the decay rate.

Another alternative is to declare that stable global monopole states cannot
exist in a world with gravity. E.g. those states must include an arbitrary
number of excitations and hence may not be normalizable in the standard way.
This would be in line with the lore that global symmetries cannot coexist with
quantum gravity. Thus our example, while special to an extent since it
requires a global spontaneously broken symmetry, may be a piece of evidence
supporting this lore.

Simply put, in the presence of a spontaneously broken global symmetry in a
world with General Relativity the solitonic structures supported by
spontaneous breaking of global symmetry provide long range Goldstone fields
sourced by the topological charge. They backreact on the geometry and induce a
local distortion around the global monopole cores, which behaves like a
‘cavity’ that traps a section of de Sitter space in the core of the monopole.
This region is nonsingular, smooth and completely ‘naked’ to the exterior: it
is not hiding behind a horizon, and is not covered by a shell of material.
Thus the relativistic pressures of the vacuum energy, which also gravitate,
allow the mass of the configuration to be negative. A deficit angle at
infinity may permit configurations of multiple such objects where the negative
masses dominate. If so, without a strict cutoff, the disaster looms. A space
of a single global monopole becomes crowded quickly.

##  4 Summary

In summary, we have argued that global monopole configurations in General
Relativity with unbroken local Poincaré invariance might be unstable to
quantum effects. They have a deficit angle and as a consequence, a negative
mass, both gravitational and inertial [27, 41]. Although their production rate
may be very small, if they can be on shell when paired up with global
antimonopoles, and have the pair energy which is negative, which might happen
when an initial global monopole is already present, obstructing the theorems
of [36, 37], normal particles which conserve total $4$-momentum could be
emitted. The rate for these processes is divergent because the unbroken local
Poincaré symmetry in the UV implies that the phase space accessible to the out
states is infinite.

One possible ‘obstruction’ to our argument is that the result of [36, 37]
might still be valid despite the fact that a global monopole sets up a conical
space. If so, our discussion could be taken as an indirect argument to this
extent, since otherwise a rapid instability will occur. A more general
positive energy theorem would exclude the pair production instability. If on
the other hand even a single pair might pop up with a negative energy, this
would raise an obstacle to having global symmetries in a QFT coupled to
universal gravity, with the only invocation of UV being to require that the
single-excitation states in asymptotically locally flat space exist
forever777A pragmatic demand might only be that the locally flat space, and
isolated individual monopoles, exist for
$~{}\mbox{\raisebox{-2.58334pt}{$\stackrel{{\scriptstyle>}}{{\sim}}$}}~{}10^{10}$
years, to fit the observations. If so, the finiteness of the observable
universe might help regulate divergence of the decay rate. and that the theory
is locally Poincaré invariant. Everything else rests on the IR properties of
gravity.

If the latter option is the outcome, then how does gravity exorcise global
monopoles, and a global symmetry needed to produce them, out of the spectrum
of the theory? At this point, this is a wide open question. One option is that
those states simply do not belong to the theory since they would not be
normalizable in the usual sense, as the number of particles diverges. Yet
another possibility is that gravity further destabilizes global monopole-
antimonopole pairs (unlike individual global monopoles which seem to be more
stable in classical gravity than without it [24]), and hence they never get on
shell with negative energy. Either way, the global monopoles merely flag
disconnected superselection sectors determined by boundary conditions at
infinity. This is unlike what happens with topological couplings in gauge
theories, where gauge symmetry allows sectors with fixed winding numbers to
mix, since after all there is no gauge symmetry here. Restricting the theory
to a single sector would effectively break global symmetry.

A question then is how this could manifest locally. Perhaps quantum gravity,
via, e.g. wormhole effects, such as those proposed long ago in [1, 3, 4, 5],
induces explicit global symmetry breaking terms such as
$\sum_{j}(\vec{A}^{j}\cdot\vec{\Phi})^{2}$ that would lift the Goldstones, and
preclude global monopoles altogether even in the classical limit. It seems
interesting to explore this further.

Acknowledgments: We thank Guido D’Amico, Roberto Emparan, David E. Kaplan,
Surjeet Rajendran, Alexander Vilenkin, Grigory Volovik and Alexander Westphal
for interesting conversations and criticism. We are grateful to CERN
Department of Theoretical Physics and to DESY Theory Group for kind
hospitality during the course of the work reported here. This work is
supported in part by the DOE Grant DE-SC0009999.

## References

  * [1] L. F. Abbott and M. B. Wise, Nucl. Phys. B 325, 687 (1989).
  * [2] T. Banks, M. Dine and N. Seiberg, Phys. Lett. B 273, 105 (1991).
  * [3] M. Kamionkowski and J. March-Russell, Phys. Lett. B 282, 137 (1992).
  * [4] R. Holman, S. D. H. Hsu, T. W. Kephart, E. W. Kolb, R. Watkins and L. M. Widrow, Phys. Lett. B 282, 132 (1992).
  * [5] R. Kallosh, A. D. Linde, D. A. Linde and L. Susskind, Phys. Rev. D 52, 912 (1995).
  * [6] J. Polchinski, Int. J. Mod. Phys. A 19S1, 145-156 (2004).
  * [7] T. Banks and N. Seiberg, Phys. Rev. D 83, 084019 (2011).
  * [8] T. Daus, A. Hebecker, S. Leonhardt and J. March-Russell, Nucl. Phys. B 960, 115167 (2020).
  * [9] N. Arkani-Hamed, L. Motl, A. Nicolis and C. Vafa, JHEP 0706, 060 (2007).
  * [10] D. Harlow and H. Ooguri, Phys. Rev. Lett. 122, no. 19, 191601 (2019).
  * [11] D. Harlow and H. Ooguri, arXiv:1810.05338 [hep-th].
  * [12] A. Nicolis, JHEP 0807, 023 (2008).
  * [13] A. Hebecker, P. Henkenjohann and L. T. Witkowski, Fortsch. Phys. 65, no.3-4, 1700011 (2017).
  * [14] M. J. Dolan, P. Draper, J. Kozaczuk and H. Patel, JHEP 1704, 133 (2017).
  * [15] P. Draper and S. Farkas, JHEP 1905, 158 (2019).
  * [16] A. Vilenkin and E. P. S. Shellard, “Cosmic Strings and Other Topological Defects,” Cambridge University Press, 2000.
  * [17] A. S. Goldhaber, Phys. Rev. Lett. 63, 2158 (1989).
  * [18] S. H. Rhie and D. P. Bennett, Phys. Rev. Lett. 67, 1173 (1991).
  * [19] D. P. Bennett and S. H. Rhie, Phys. Rev. Lett. 65, 1709-1712 (1990).
  * [20] I. Chuang, B. Yurke, R. Durrer and N. Turok, Science 251, 1336-1342 (1991).
  * [21] G. W. Gibbons, Lect. Notes Phys. 383, 110-138 (1991).
  * [22] L. Perivolaropoulos, Nucl. Phys. B 375, 665 (1992).
  * [23] A. Achucarro and J. Urrestilla, Phys. Rev. Lett. 85, 3091-3094 (2000).
  * [24] H. Watabe and T. Torii, Phys. Rev. D 66, 085019 (2002).
  * [25] A. Achucarro, B. Hartmann and J. Urrestilla, JHEP 07, 006 (2005).
  * [26] M. Barriola and A. Vilenkin, Phys. Rev. Lett. 63, 341 (1989).
  * [27] D. Harari and C. Lousto, Phys. Rev. D 42, 2626 (1990).
  * [28] A. D. Linde, Phys. Lett. B 327, 208 (1994).
  * [29] A. Vilenkin, Phys. Rev. Lett. 72, 3137 (1994).
  * [30] I. Cho and A. Vilenkin, Phys. Rev. D 56, 7621 (1997).
  * [31] I. Cho and J. Guven, Phys. Rev. D 58, 063502 (1998).
  * [32] J. M. Cline, S. Jeon and G. D. Moore, Phys. Rev. D 70, 043543 (2004).
  * [33] D. E. Kaplan and R. Sundrum, JHEP 0607, 042 (2006).
  * [34] R. Emparan and J. Garriga, JHEP 03, 028 (2006).
  * [35] G. Dvali, [arXiv:1107.0956 [hep-th]].
  * [36] R. Schon and S. T. Yau, Commun. Math. Phys. 79, 231-260 (1981).
  * [37] E. Witten, Commun. Math. Phys. 80, 381 (1981).
  * [38] L. Susskind, [arXiv:hep-th/9501106 [hep-th]].
  * [39] N. Kaloper and D. Kiley, JHEP 0603, 077 (2006).
  * [40] J. Ehlers, I. Ozsvath, E. L. Schucking and Y. Shang, Phys. Rev. D 72, 124003 (2005).
  * [41] A. Einstein, L. Infeld and B. Hoffmann, Annals Math. 39, 65 (1938).
  * [42] G. W. Gibbons and P. J. Ruback, Phys. Rev. Lett. 57, 1492 (1986).
  * [43] U. Nucamendi and D. Sudarsky, Class. Quant. Grav. 14, 1309-1327 (1997).
  * [44] U. Nucamendi and D. Sudarsky, Class. Quant. Grav. 17, 4051-4058 (2000).
  * [45] H. Bondi, Rev. Mod. Phys. 29, 423 (1957).
  * [46] G. R. Dvali, [arXiv:hep-th/0004057 [hep-th]].
  * [47] G. Volovik, JETP Letters 112, 505 (2020).


Supplementary Information to accompany the article

“Complete field-induced spectral response of the spin-1/2 triangular-lattice
antiferromagnet CsYbSe2”

Tao Xie, A. A. Eberharter, Jie Xing, S. Nishimoto, M. Brando, P. Khaneko, J.
Sichelschmidt, A. A. Turrini, D. G. Mazzone,

P. G. Naumov, L. D. Sanjeewa, N. Harrison, A. S. Sefat, B. Normand, A. M.
Läuchli, A. Podlesnyak and S. E. Nikitin

##  S1 Characterization

###  S1.1 Single-Crystal X-Ray Diffraction

We have performed a careful investigation of the lattice structure of our
CsYbSe2 crystals by analysing several batches of single crystals using a
Bruker Quest D8 single-crystal X-ray diffractometer. The structure was refined
by the Rietveld method using the FullProf software package Rodríguez-Carvajal
(1993), which delivered excellent structural solutions with no evidence for
site mixing. In the refinements, we treated the anisotropic displacement
parameters and site occupancies as free variables, and indeed observed no site
disorder. The crystallographic data are presented in Table S1, and the result
of a Rietveld refinement is shown in Fig. S1; we have also reported this
crystal structure in the Cambridge Crystallographic Data Centre (CCDC)
(1952075) Xing  _et al._ (2020). We also found no significant residual
electron density, which would suggest the presence of some interstitial atoms.
All of these results confirmed the high quality of our single crystals.

Figure S1: Rietveld refinement of single-crystal X-ray diffraction data at 100 K. $F^{2}_{\mathrm{obs}}$ and $F^{2}_{\mathrm{cal}}$ denote respectively the observed and calculated structure factors. Table S1: Crystallographic data for CsYbSe2 determined by single-crystal X-ray diffraction. empirical formula | CsYbSe2
---|---
formula weight (g/mol) | 463.87
$T$, K | 100
Crystal habit | red plates
Crystal dimensions, mm | $0.14\times 0.10\times 0.03$
crystal system | hexagonal
space group | P63/mmc (No. 194)
_a_ , Å | 4.1466(2)
_c_ , Å | 16.5050(1)
volume, Å3 | 245.77(3)
Z | 2
Density (calc), g/cm3 | 6.268
$\mu(\mathrm{MoK_{\alpha}})$, mm-1 | 40.932
_F_(000) | 386
$T_{\rm{max}}$, $T_{\rm{min}}$ | 0.4415, 1.0000
$\theta$ range for data collection | 2.47-29.93
reflections collected | 3895
final $R~{}{[}I\textgreater 2\sigma(I){]}R_{1},R_{w2}$ | 0.0274/0.0833
final $R$ (all data) $R_{1},R_{w2}$ | 0.0284/0.0843
Goodness of fit, $F^{2}$ | 1.343

###  S1.2 Crystal Structure and Next-Neighbour Magnetic Interaction

One of the key questions in the analysis of every compound in the Yb-
delafossite family is whether the ground state of the $S=1/2$ system can be a
quantum spin liquid (QSL) at zero applied magnetic field. In our work we have
confirmed the presence of 120∘ order in CsYbSe2, at temperatures below 0.4 K
and at least over a spatial range far exceeding the lattice constant (Fig. 1c
of the main text and Fig. S6(e) below). The same order is found in KYbSe2
Scheie  _et al._ (2021), but notably not in NaYbSe2 Dai  _et al._ (2021).
Factors destabilizing this order include a possible frustrated coupling
between triangular lattice (TL) planes and a possible next-neighbour
interaction within the planes, where (as noted in the main text) recent
numerical studies have achieved partial agreement that the ground state is a
QSL in the regime $0.06\lesssim J_{2}/J_{1}\lesssim 0.15$ Kaneko  _et al._
(2014); Li  _et al._ (2015); Zhu and White (2015); Hu  _et al._ (2015); Iqbal
_et al._ (2016); Hu  _et al._ (2019).

Table S2: Crystallographic information for selected Yb-selenide delafossites Xing  _et al._ (2021).  Chemical formula | NaYbSe2 | KYbSe2 | RbYbSe2 | CsYbSe2
---|---|---|---|---
Crystal system | trigonal | trigonal | trigonal | hexagonal
Space group | R$\overline{3}m$ (No. 166) | R$\overline{3}m$ (No. 166) | R$\overline{3}m$ (No. 166) |  $P6_{3}/mmc$ (No. 194)
Stacking type | ABC | ABC | ABC | AA
Ionic radius, Å | 1.02 (Na+) | 1.38 (K+) | 1.52 (Rb+) | 1.67 (Cs+)
Interlayer distance, Å | 6.92 | 7.56 | 7.88 | 8.25

Addressing first the issue of layer stacking, the majority of Yb delafossites
have a structure described by the R$\overline{3}m$ space group, as listed in
Table. S2, which has an ABC stacking that does indeed suggest frustration of
antiferromagnetic (AF) interlayer couplings. However, some members of the
series with larger alkali-metal ions, including Cs (Table. S2), have a
structure with space group P$6_{3}/mmc$ that has AA stacking, and hence no
interlayer frustration. Thus the available materials examples tend to suggest
that the stacking of TL layers is not a relevant factor. Turning to $J_{2}$,
in our work we have deduced a weak next-neighbour interaction in CsYbSe2,
$J_{2}=0.03J_{1}$. The authors of Ref. Scheie  _et al._ (2021) deduced a value
$J_{2}=0.05J_{1}$ in KYbSe2, placing it closer to the QSL regime. While the
authors of Ref. Dai  _et al._ (2021) were not able to determine a $J_{2}$
value, one may certainly suggest that the larger lattice constant inherent to
the members with larger alkali-metal ions (Table. S2) causes a reduction of
$J_{2}$ and hence an increasing stabilization of the 120∘-ordered ground state
at base temperature in the Yb-delafossite materials.

Figure S2: Isothermal magnetization measured in a pulsed magnetic field. The
black curve shows the raw magnetization data, the red dashed line the
estimated van Vleck contribution and the blue curve the intrinsic
magnetization obtained by subtraction of the van Vleck part. The orange
symbols are the low-field (up to 7 T) magnetization data measured in a MPMS-7.
Open circles show the bulk moment extracted from the field-dependence of the
$(0,0,4)$ Bragg peak measured by neutron diffraction at $T<0.05$ K. The solid
red curve in the lower panel shows the first derivative of the magnetization,
$dM/dB$.

###  S1.3 Magnetization in Pulsed Magnetic Field

We measured the isothermal magnetization at $T=0.4$ K in pulsed magnetic
fields up to 60 T, and for calibration in a MPMS-7 with a 3He insert at fields
up to 7 T. We note that the magnetic field in all our experiments is applied
in the $ab$ plane, unless otherwise stated, and that we do not distinguish
between in-plane directions. In the raw magnetization data shown by the black
solid line in Fig. S2, we observe a clear and continuous increase above a
saturation field of approximately 10 T. This contribution results from van
Vleck paramagnetism and can be approximated in lowest order by a linear
function (the dashed red line), from which the van Vleck susceptibility may be
estimated as $\chi^{VV}\approx 0.0152~{}\mu_{B}$/T. Subtracting this
contribution from the raw data Ranjith  _et al._ (2019a, b); Ding  _et al._
(2019) leaves an intrinsic magnetization (solid blue line in Fig. S2, shown
also in Fig. 1b of the main text) displaying near-perfect saturation above a
field $B_{\mathrm{Sat}}$ that we estimate most accurately from the excitation
spectrum of the fully polarized phase (Fig. 2b of the main text). In Fig. S2
we show also the bulk magnetization obtained from the field-induced
enhancement of the (0 0 4) Bragg peak measured on TASP at the Paul Scherrer
Institute (PSI), and these results agree quantitatively with the direct
measurements. We defer the accurate modelling of these data, which we
performed to obtain the orange line in Fig. 1b of the main text, to Sec. S4.

Figure S3: Temperature-dependence of the specific heat measured in a range of
applied in-plane magnetic fields ($\mathbf{B}\perp{\hat{c}}$).

###  S1.4 Specific Heat

A sharp peak in $C(T)$ can be used to establish the presence of long-range
order in the system. In Fig. S3(a) we group the applied fields, namely those
below 3 T and above 6 T, in which no clear peak appears in $C(T)$ at all, only
a broad hump at the low fields and a weak shoulder at the high fields. In Fig.
S3(b) we show that relatively sharp, $\lambda$-shaped peaks can be found at 3,
4 and 5 T, in partial correspondence with the long-range order of the UUD
phase found by neutron diffraction. The temperatures of these sharp peaks are
indicated as solid circles in the phase diagram shown in Fig. 1d of the main
text and the temperatures of the features found at $B\leq 2$ T and $B\geq 6$ T
as open circles.

In more detail, the specific heat at zero field in Fig. S3(a) shows a broad
hump with no trace of a phase transition. This should be contrasted with the
neutron diffraction measurements in Fig. 1c of the main text, which display a
series of weak magnetic Bragg peaks with a clear onset at $T\simeq 0.4$ K.
However, the analysis of the correlation lengths presented in Sec. S2.2 leads
to the result that these are finite at zero field, signalling a lack of true,
long-ranged AF order even at $T=0.02$ K, which is consistent with the absence
of a sharp peak in the corresponding $C(T)$ data. The specific heat does show
sharp, $\lambda$-shaped peaks in the field range 3-5 T [Fig. S3(b)] that
reflect phase transitions into a long-range-ordered phase on, and apparently
near, the 1/3-plateau state. These thermodynamic data are confirmed by the
field-dependence of the (1/3, 1/3, 1) Bragg peak presented in Fig. 1b of the
main text and analysed in detail in Fig. S7 below.

Figure S4: (a) Photographs of the CsYbSe2 crystals coaligned and assembled for
neutron scattering experiments. (b) Rocking curve of the coaligned sample
measured by neutron diffraction, with a Gaussian fit shown by the solid line.

##  S2 Neutron Scattering Experiments

###  S2.1 Sample preparation

We coaligned around 200 single crystallites on copper plates, as shown in Fig.
S4(a), to obtain a mosaic sample with a total mass of approximately 0.5 g. The
rocking curve obtained by a neutron diffraction measurement of one selected
Bragg peak [Fig. S4(b)] shows that the mosaicity of this sample is
approximately $2\degree$ FWHM.

Figure S5: Constant-energy slices in the $(H~{}H~{}L)$ plane taken at zero
energy transfer under different magnetic fields, with the data symmetrized
about the $[H~{}H~{}0]$ axis. The strong spots at $Q=(0,~{}0,~{}L)$ for even-
integer $L$ are nuclear Bragg peaks.

###  S2.2 Elastic Neutron Scattering at CNCS

Figure S5 shows two-dimensional (2D) slices through the measured intensity
dataset at zero energy in the ($H~{}H~{}L$) plane at fields from 0 to 5 T. At
0 T we observe weak but nonetheless clear intensity peaks at
$Q=(1/3,~{}1/3,~{}L)$ for odd-integer $L$, which become weaker and almost
invisible at 1 T. For $B\geq 2$ T we find rods of magnetic intensity elongated
along the $[1/3~{}1/3~{}L]$ direction, with peaks at the same $Q$ values.
These results confirm the 2D nature of the magnetic system and the persistence
of threefold periodicity in the plane of the triangular lattice at all applied
fields.

For a specific magnetic Bragg peak, the spin correlation length can be
estimated using the formula $\xi=2\pi/\sqrt{w^{2}-R^{2}}$ Young  _et al._
(2013), where $w$ is the full width at half maximum (FWHM) height of the
magnetic peak and $R$ is the instrumental momentum resolution at this peak. To
estimate this resolution, we first prepared 1D intensity cuts along the [$H$
$H$ 0] and [0 0 $L$] directions at the nuclear peaks (002) [Figs. S6(a) and
S6(b)], (004) and (006), then fitted a Gaussian function to them. The FWHM
values of these three nuclear peaks reflect the momentum resolution of the
instrument perpendicular (for [$H$ $H$ 0]) and parallel (for [0 0 $L$]) to
${\hat{Q}}$ at each reciprocal lattice point. Polynomial fits to the FWHM as a
function of $L$ [Figs. S6(c)] yields an approximate estimate of the resolution
as a function of $\left|\bf{Q}\right|$ for an arbitrarily chosen momentum
transfer, $\mathbf{Q}$, in the perpendicular and parallel directions [Fig.
S6(d)].

To obtain the FWHM of the magnetic Bragg peaks for a quantitative
determination of the correlation length, we fitted 1D cuts through the
magnetic peaks to the Voigt function, which is defined as

$\displaystyle y(x)$ $\displaystyle=y_{0}+(f_{1}*f_{2})(x),$ (S1)

where the $*$ denotes a convolution. In Eq. (S1),
$f_{1}(x)=\frac{2A}{\pi}\frac{w_{L}}{{4(x-x_{c}})^{2}+w_{L}^{2}}$ is the
Lorentz function, with $x_{c}$ as the peak centre and $w_{L}$ as the FWHM, and
$f_{2}(x)=\sqrt{\frac{4\ln 2}{\pi}}\frac{1}{w_{G}}\exp(-\frac{4\ln
2}{w_{G}^{2}}x^{2})$ is a Gaussian with peak centre $x=0$, unit peak area and
FWHM $w_{G}$, which should be fixed to the instrumental resolution,
$R(\left|\bf{Q}\right|)$, for a specific $\bf{Q}$. The FWHM of the Voigt
function is then $w_{V}$ = $0.5346w_{L}$ \+ $\sqrt{0.2166\cdot
w_{L}^{2}+w_{G}^{2}}$. By inserting $w_{V}$ for $w$ and
$R(\left|\bf{Q}\right|)$ for $R$ into the above expression for $\xi$, we
estimate the correlation lengths at $\mathbf{Q}$ = (1/3, 1/3, 1) for each
different field. Figure S6(e) shows the correlation lengths we deduce from all
of our CNCS data, which from our choice of the (0 0 $L$) series of nuclear
peaks correspond to the in- and out-of-plane correlations. At $B=0$ we obtain
the values $\xi_{ab}=60(7)$ Å and $\xi_{c}=23(5)$ Å quoted in the main text.
We observe that the correlation lengths dip at 1 T before $\xi_{ab}$ rises
strongly, to values in excess of 200 Å, on and just above the 1/3 plateau. To
gauge the meaning of this number, in Fig. S6(f) we show a 1D cut along
$[H~{}H~{}1]$ at $B=4$ T: the horizontal blue bar is the $|Q|$-resolution at
the magnetic Bragg peak (1/3, 1/3, 1), which is clearly resolution-limited,
and thus the intrinsic $\xi$ values diverge as expected for the long-range-
ordered UUD state.

We also state for completeness that the experimental data in the upper panel
of Fig. 1c of the main text were integrated over the interval
$K=[-0.05,0.05]$, symmetrized according to the crystal symmetry and unfolded
for visual clarity, meaning that the intensities at $\pm L$ are equivalent.

Figure S6: (a,b) 1D intensity cuts at zero energy along the $[0~{}0~{}L]$ and
$[H~{}H~{}2]$ directions at nuclear (0, 0, 2) Bragg peak at 0 T. (c) FWHM of
nuclear peaks obtained from cuts along the [$H~{}H$ 0] direction (black) and
$[0~{}0~{}L]$ direction (red) at $Q$ = (0, 0, 2), (0, 0, 4) and (0, 0, 6).
Solid lines in both panels show polynomial fits. (d) Perpendicular and
parallel $|Q|$-resolution estimated from the FWHMs shown in panel (c). (e)
Perpendicular ($\xi_{\perp}$) and parallel ($\xi_{\parallel}$) correlation
lengths calculated from the magnetic Bragg peaks shown in Fig. S5. For this
geometry $\xi_{\perp}\equiv\xi_{ab}$, the in-plane correlation length, and
$\xi_{\parallel}\equiv\xi_{c}$, the out-of-plane correlation length. (f) An
example 1D cut along the [$H~{}H$ 1] direction through the magnetic Bragg peak
($1/3$, $1/3$, 1) at 4 T. The horizontal blue bar represents the instrumental
$|Q|$-resolution, which indicates a resolution-limited peak at 4 T.

###  S2.3 Elastic Neutron Scattering at TASP

The dependence of the integrated intensity of the magnetic Bragg peak (1/3,
1/3, 1) on magnetic field is summarized in Fig. 1b of the main text and the
dependence on temperature in Fig. 1c. Figure S7(a) shows the raw elastic
neutron scattering data measured on TASP in scans along the $[H~{}H~{}1]$
direction at base temperature for a series of magnetic fields. The peak
intensity is clearly suppressed by a weak field, with a minimum around 1 T,
before becoming sharper and stronger as the field is increases, reaching a
maximum at approximately 4 T. This intensity then decreases beyond 5 T to very
small values higher in the V regime. Figure S7(b) shows the data obtained for
the same scan at zero field for multiple temperatures below 1 K, where it is
clear that a robust peak is present at $T=0.35$ K and becomes sharper and
stronger towards base temperature.

Figure S7: (a) Intensity of the ($1/3$, $1/3$, 1) magnetic Bragg peak measured
at $T<0.05$ K for a series of magnetic fields, which are offset by 35 units
for clarity. Solid lines and shading show Gaussian fits. (b) Intensity of this
peak measured at zero field for a series of temperatures below 1 K, which are
offset by 10 units for clarity.

###  S2.4 Two-Dimensional Excitation Character

The spin excitations we observe in our measurements retain their highly 2D
nature over the full energy range and under all applied magnetic fields. In
the representative 2D constant-energy slices presented in Fig. S8, the
intensity distribution in all cases takes the form of multiple rods extending
largely unchanged along the $L$ direction. This definitive proof of 2D spin
excitations is fully consistent with previous reports on CsYbSe2 Xing  _et
al._ (2019) and related delafossites Dai  _et al._ (2021); Scheie  _et al._
(2021). As a consequence we can integrate our intensity data over a wide $L$
range when analysing the in-plane spin dynamics, and for the data shown both
below and in Fig. 3 of the main text, this integration range, $1.2\leq L\leq
3.8$, is denoted by the $L$ label 2.5.

Figure S8: Constant-energy 2D intensity slices at different energies and
magnetic fields. The data have been symmetrized about the [0 0 $L$] axis.
Figure S9: Definition of the subtracted background in the INS spectra and
illustration of selected INS spectra before and after background subtraction.
(a) Left: raw data for the spin excitation spectrum along $[H~{}H~{}2.5]$ at 5
T and 70 mK; right: corresponding spectrum obtained after background
subtraction. This background is extracted from the shaded area $0\leq H\leq
0.05$). (b) Open black squares show the 1D energy cut obtained by integrating
over the shaded region in panel (a). The solid black line shows a Lorentz fit
to this cut with a linear background, the blue line shows the Lorentz peak fit
alone and the open red circles show the difference, which is used as the real
background for subtraction purposes. (c) INS spectra before and after
background subtraction for $B=3$ T. (d) Spectra before and after subtraction
for $B=4$ T.

###  S2.5 Background Definition and Subtraction for INS Spectra

The INS spectrum is always contaminated by a background contribution that
arises from incoherent neutron scattering and from scattering due to the
sample environment. Raw spectral data obtained on CNCS for the $[H~{}H~{}2.5]$
direction (meaning with the broad $L$ integration described above) are shown
in the left panels of Figs. S9(a,c,d) for three different magnetic fields at
our base temperature of 70 mK. For an accurate characterization of the
background, we make use of the fact that the 5 T spectrum at the $\Gamma$
point, $I(\mathbf{Q}=\mathbf{0},E)$, consists of a single and well-defined
inelastic peak near 1 meV. We assume that an appropriate 1D intensity cut,
prepared from the shaded region in Fig. S9(a) and shown by the black points in
Fig. S9(b), can be described within the energy window 0.5 meV $\leq E\leq$ 1.5
meV by a Lorentzian peak and a linear background. Thus we fitted the measured
signal to the form

$\displaystyle
I(E)=a_{0}+a_{1}E+\frac{I_{0}W_{0}^{2}}{(E-E_{0})^{2}+W_{0}^{2}},$ (S2)

where $I_{0}$, $W_{0}$ and $E_{0}$ characterize respectively the intensity,
width and centre of the inelastic peak. By subtracting the fitted Lorentzian,
shown by the blue line in Fig. S9(b), from the raw spectrum, we obtained a
residual intensity that we consider as being fully representative of the
background [shown by the red points in Fig. S9(b)]. The results of subtracting
this background are shown in the right panel of Fig. S9(a), which shows clear
limits to the extent of each continuum. Because our collected spectra provide
no evidence that the background varies with $\mathbf{Q}$, we subtracted the
same form from our INS data cuts over the full reciprocal space. Similarly,
the background has no dependence on the magnetic field and examples of full
background-subtracted spectra at $B=3$ and 4 T are shown in Figs. S9(c,d). We
stress that this method provides a minimal parameter-free model for describing
the background, and its simplicity far outweighs the disadvantage of a minor
(statistically insignificant) oversubtraction appearing at some high $H$
values.

##  S3 Linear Spin-Wave Theory

In our analysis we use linear spin-wave theory (SWT) to achieve two separate
goals. Working above the saturation field, $B_{\mathrm{sat}}$, we make use of
the fact that the excitations are guaranteed to be well defined spin waves to
obtain the most accurate available fit of the magnetic interaction parameters
of CsYbSe2 (Sec. S3.1). Working below $B_{\mathrm{sat}}$, we use linear SWT
for a preliminary indication of the locations and energy scales of putative
$\Delta S=1$ excitations in the TLHAF, and hence of departures from
semiclassical magnetism arising due to quantum corrections. In Sec. S3.2 we
outline how the SWT results shown in Fig. 4 of the main text were obtained.

###  S3.1 Fitting of Magnetic Interactions

Figure S10: (a) Fit quality shown as a function of assumed XXZ-model
parameters $J_{1}$ and $\Delta*J_{1}$ for the dominant nearest-neighbour
interaction. Here $J_{2}$ has been fixed to its optimal value of 0.011 meV.
(b) Fit quality shown as a function of $J_{1}$ and $J_{2}$, assuming $\Delta\
=1$. Red crosses in both panels indicate the position of the global minimum.

As described in the main text, we determine the magnetic interactions by using
the 11 T dataset from CAMEA, which shows a single, sharp magnon mode with its
maximum at the $\Gamma$ point (Fig. 2b of the main text). We quantified the
location of this mode at 17 points in reciprocal space and fitted these to the
Hamiltonian S3 using SpinW Toth and Lake (2015).

$\displaystyle\mathcal{H}$ $\displaystyle=$ $\displaystyle J_{1}\sum_{\langle
i,j\rangle}({S}^{x}_{i}{S}^{x}_{j}+{S}^{z}_{i}{S}^{z}_{j}+\Delta{S}^{y}_{i}{S}^{y}_{j})$
(S3) $\displaystyle+J_{2}\sum_{\langle\langle
i,j\rangle\rangle}\mathbf{S}_{i}\cdot\mathbf{S}_{j}-\mu_{B}g_{ab}B\sum_{i}{S}^{z}_{i}.$

The optimal fit in the 3D parameter space of $J_{1}$, $\Delta$ and $J_{2}$ is
$J_{1}=0.395(7)$ meV, $\Delta\times J_{1}=0.39(2)$ meV and $J_{2}$ = 0.011(4)
meV. In Fig. S10 we illustrate the quality of this fit by showing two 2D
cross-sections: fixing $J_{2}$ to its optimal value [Fig. S10(a)] shows that
$J_{1}=\Delta\times J_{1}$, i.e. the nearest-neighbour interaction is
isotropic within the precision of our measurements; setting $\Delta\ =1$ [Fig.
S10(b)] allows us to optimize $J_{2}$, and hence to conclude that magnetic
Hamiltonian of CsYbSe2 is described by the $J_{1}$-$J_{2}$ Heisenberg model
with $J_{1}=0.395(8)$ meV and $J_{2}=0.011(4)$ meV ($J_{2}/J_{1}$ = 0.03).

###  S3.2 SWT Spectra

Linear SWT ceases to provide a complete description of the spectrum of the
TLHAF as soon as the field is lowered below $B_{\mathrm{sat}}$, when quantum
corrections become finite. The TLHAF in an in-plane magnetic field has long
been known Chubokov and Golosov (1991) to exhibit a deformed 120∘ phase (Y),
the 1/3 plateau phase (UUD), and a V (or 2:1) phase, with two of the spins
parallel, below saturation, as illustrated schematically in Fig. 1d of the
main text. At the mean-field level, the magnetization is perfectly linear in
field from 0 to $B_{\mathrm{sat}}=9.6$ T [Fig.S11], which is obtained exactly
due to the classical nature of the problem at $B\geq B_{\mathrm{sat}}$.
However, $M(B)$ has no plateau around $B_{\mathrm{sat}}/3$, as this phase is
stabilized only at higher order in $1/S$. In linear SWT, there is a large
family of classically degenerate states at any field $B<B_{\mathrm{sat}}$, and
the precise orientation of the Y, UUD or V state at any given field is set
manually in our analysis. The spin excitations around this fixed ground state
were then calculated using SpinW and the resulting spectra at all selected
fields in our experimental range ($B=0$-11 T) are shown as the orange lines in
Fig. 4 of the main text.

##  S4 Magnetization Response of the TLHAF

To interpret the magnetic response despite the effects of the finite-
temperature rounding, we estimate $B_{\mathrm{Sat}}$ from the TLHAF model
parameters obtained by high-field INS (Fig. S11 below) and perform a grand
canonical DMRG calculation of the full magnetization curves that allows us to
deduce the boundaries of the 1/3 magnetization plateau.

Figure S11: Magnetization, $M(B)$, of the $J_{1}$-$J_{2}$ TLHAF with
$J_{2}$/$J_{1}$ = 0.03 calculated by grand canonical DMRG and mean-field
method. The inset shows the 75-site open cluster used for the DMRG
calculation.

For a state-of-the-art numerical determination of the magnetization response
of the $S=1/2$ TLHAF, we apply the grand canonical DMRG method Hotta  _et al._
(2013). This technique computes the infinitesimally small magnetization
response to a change in the applied field, and thus the physical quantities it
provides mimic the thermodynamic limit to approximately 1 part in $10^{3}$ for
a 2D system. The method is based on a graded division of a finite-sized
cluster into a centre region and edge regions, such that the centre reproduces
the continuous bulk response by using the near-zero-energy states of the edge
as a type of buffer. In a system of fixed size and shape, the grading is
introduced by modulating the energy scale with the externally imposed function

$\displaystyle f(r)=\frac{1}{2}\left[1+\cos\Big{(}\frac{\pi
r}{R}\Big{)}\right],$ (S4)

which deforms the Hamiltonian smoothly from its standard energy at the centre
of the system ($r=0$) to zero at the open cluster edges ($r=R$). After
obtaining the lowest eigen-wavefunction of the deformed Hamiltonian, the
magnetization can be read as $M=\langle S^{z}(r=0)\rangle$, because this
wavefunction optimizes the expectation value to its thermodynamic limit at any
given magnetic field. To mimic the results of a bulk measurement,
appropriately weighted for the different possible directions of the magnetic
field with respect to the triangular lattice, we performed our calculations
using the hexagonal 75-site cluster shown in the inset of Fig. S11. The
calculations were performed for the $J_{1}$-$J_{2}$ TLHAF using the
interaction parameters deduced in Subsec. S3.1 and an in-plane $g$-factor of
$g_{ab}=3.2$, as deduced from ESR and INS in the main text. The result shown
in Fig. S11 is that reproduced in Fig. 1d of the main text, and was used for
our extraction of the lower and upper boundaries ($B_{\rm l}$ and $B_{\rm u}$)
of the 1/3 plateau.

Figure S12: (a) Example of the cylinder geometry used in our MPS calculations.
The three sublattices are indicated by different colours. Numbers denote the
index in the linear arrangement of the MPS and the three centre sites (used as
initial excitation sites for the time evolution) are marked by red indices.
Grey cells illustrate the wrapping specified by the periodic (XC) boundary
conditions. (b) Local magnetization $\braket{S_{i}^{z}}$ measured in the the
ground state in the $1/3$ magnetization plateau showing up-up-down order. The
three center sites are marked in green.

##  S5 MPS Calculations

The wave function of the spin system is represented by a matrix-product state
(MPS), which assigns a tensor to each lattice site. The accuracy of the MPS
approximation is controlled by the maximum bond dimension, $\chi$, of these
tensors. To study the TLHAF we tested values of this parameter up to
$\chi=1024$, and will illustrate its role below. We used the Python package
TenPy Hauschild and Pollmann (2018), specifically its algorithms for the
ground state and for time-evolution, as well as its support for preserving the
$U(1)$ symmetry of a model.

We start by describing the cylinder geometry used for MPS simulations. We
implement a cylinder with circumference $C$, length $L$ and the “XC” boundary
conditions shown in Fig. S12(a) Szasz  _et al._ (2020); Scheie  _et al._
(2021). This geometry provides a high momentum resolution along the Brillouin-
zone path $\Gamma$-K-M that is of primary interest for the comparison with
experiment. In our MPS calculations we determine the time-dependent spin-spin
correlation function

$\displaystyle\begin{split}C&{}^{\alpha\beta}_{\mathbf{r}}(\mathbf{x},t)\coloneqq\braket{\hat{\mathcal{S}}^{\alpha}_{\mathbf{r}+\mathbf{x}}(t)\hat{\mathcal{S}}^{\beta}_{\mathbf{r}}(0)}\\\
&=\braket{[\hat{S}^{\alpha}_{\mathbf{r}+\mathbf{x}}(t)-\braket{\hat{S}^{\alpha}_{\mathbf{r}+\mathbf{x}}(t)}][\hat{S}^{\beta}_{\mathbf{r}}(0)-\braket{\hat{S}^{\beta}_{\mathbf{r}}(0)}]}\\\
&=\braket{\hat{S}^{\alpha}_{\mathbf{r}+\mathbf{x}}(t)\hat{S}^{\beta}_{\mathbf{r}}(0)}-\braket{\hat{S}^{\alpha}_{\mathbf{r}+\mathbf{x}}(t)}\braket{\hat{S}^{\beta}_{\mathbf{r}}(0)},\\\
\end{split}$ (S5)

where the disconnected part is subtracted in order to remove the Bragg peaks
in $C^{zz}$. Using $\ket{0}$ to denote the ground state and $E_{0}$ for the
ground-state energy of $H$, we rewrite Eq. (S5) in the form

$\displaystyle C^{\alpha\beta}_{\mathbf{r}}(\mathbf{x},t)$
$\displaystyle=\braket{0}{\hat{S}^{\alpha}_{\mathbf{r}+\mathbf{x}}\
e^{-i(H-E_{0})t}\ \hat{S}^{\beta}_{\mathbf{r}}}{0}$
$\displaystyle-\braket{0}{\hat{S}^{\alpha}_{\mathbf{r}+\mathbf{x}}}{0}\braket{0}{\hat{S}^{\beta}_{\mathbf{r}}}{0},$
(S6)

where $\mathbf{r}$ is the site at which the initial spin operator is applied
and $\mathbf{x}$ is the vector separation in the two-point correlator. Because
of the residual $U(1)$ symmetry of the Heisenberg Hamiltonian in a field, the
total magnetization in the $z$ direction is conserved under time evolution. We
take advantage of this symmetry to increase the performance of the calculation
by restricting the MPS to a single magnetization sector. It follows that only
the correlators which satisfy $\alpha\beta\in\\{zz,+-,-+\\}$ are non-zero. The
other spin-spin correlator functions can be obtained based on the identity
$C^{xx}_{\mathbf{r}}(\mathbf{x},t)\equiv
C^{yy}_{\mathbf{r}}(\mathbf{x},t)\equiv\frac{1}{2}(C^{+-}_{\mathbf{r}}(\mathbf{x},t)+C^{-+}_{\mathbf{r}}(\mathbf{x},t))$.

In a typical time-evolution process one selects a single site, $\mathbf{r}$,
at the centre of the cylinder and uses a single time-evolved state to obtain
the spectral function. However, the physics of the TLHAF in a field is such
that there is a UUD-type symmetry breaking of the $z$-axis magnetization at
all finite fields. Our MPS set-up is therefore chosen such that the ground
state in a field is a representative symmetry-broken state with a pronounced
UUD-type local $z$-axis magnetization pattern, as shown in Fig. S12(b). In
fact it is more efficient to repeat the dynamical evolution using three
distinct starting sites than it is to use a single MPS with all three
symmetry-broken states superposed in a spatially uniform way. The reason for
this is that the approximate tripling of the bond dimension required to
represent this symmetric MPS leads to a much more expensive time evolution,
because of the expected $\chi^{3}$ scaling, than using three different
starting states at a lower bond dimension would do.

To obtain the spectral functions that correspond to the experiment, it is
necessary to restore the lost translation symmetry. To this end, the spectral
function $S_{\mathbf{r},\alpha\beta}(\mathbf{Q},\omega)$ is calculated for the
three sites $\mathbf{r}$ corresponding to the three sublattices of the central
unit cell of a single symmetry-broken groundstate. We then take the equally
weighted average to obtain the symmetric spectral function

$S_{\alpha\beta}(\mathbf{Q},\omega)=\frac{1}{3}\sum_{\mathbf{r}}S_{\mathbf{r},\alpha\beta}(\mathbf{Q},\omega).$
(S7)

This procedure is repeated for all $\alpha\beta\in\\{zz,+-,-+\\}$, meaning
that in total nine time-evolution runs are required at every value of the
magnetic field, $B$.

In addition, symmetries of the Heisenberg Hamiltonian allow us to obtain the
negative-time correlators at no additional computational cost as the complex
conjugate of the positive-time correlators, i.e.
$C^{\alpha\beta}_{\mathbf{r}}(\mathbf{x},-t)=\overline{C^{\alpha\beta}_{\mathbf{r}}(\mathbf{x},t)}$.
Putting these together, Eq. (5) of the main text simplifies to the form

$\displaystyle\begin{split}S_{\alpha\beta}(\mathbf{Q},\omega)=&\frac{2}{3}\sum_{\mathbf{r}}\int_{0}^{\infty}\mathrm{d}t\sum_{\mathbf{x}}e^{-i\,\mathbf{x}\cdot\mathbf{Q}}\\\
&\left[\cos(\omega t)\,{\rm
Re}\,C^{\alpha\beta}_{\mathbf{r}}(\mathbf{x},t)-\sin(\omega t)\,{\rm
Im}\,C^{\alpha\beta}_{\mathbf{r}}(\mathbf{x},t)\right],\end{split}$ (S8)

where the outer sum runs over the three centre sites and the correlators are
multiplied by a Gaussian envelope in space and time as described in the
Methods section. In Fig. S13 we show spectra obtained using this procedure in
the UUD (1/3-plateau) phase, at a magnetic field equivalent to $B=3$ T, to
illustrate the degree to which our calculations have converged to the spectral
function of the infinite system. We observe that increasing the cylinder
circumference from $C=6$ to 8, the length from $L=30$ to 60 or the bond
dimension from $\chi=256$ to 512 make only minor changes to the quality of the
spectra. As a result of this study we perform the majority of our calculations
with $C=6$ and $L=30$, using bond dimensions $\chi=512$ or 1024, as stated in
the Methods section.

Figure S13: Comparison of spectral functions conmputed using the MPS method
for different system sizes and bond dimensions. Here we show the quantity
$S_{xx}(\mathbf{Q},\omega)+S_{zz}(\mathbf{Q},\omega)$ at a magnetic field
corresponding to the experimental value $B=3$ T. The columns show increasing
values of the bond dimension, $\chi$, while the rows show cylinders of
different sizes, as labelled in the panels in the first column.

##  S6 Comparison of Measured and Calculated Spectral Functions

Figure S14: Constant-$\mathbf{Q}$ cuts through the experimental (symbols) and
calculated (solid lines) spectral functions for six different magnetic fields
and six $\mathbf{Q}$ values along the [$H$ $H$ $2.5$] line. For clarity the
cuts are shown with a vertical offset of 0.2.

For a quantitative comparison of the MPS calculations with the experimental
INS data, we present our results in two different ways. In Figs. 3b, 3d, 3f,
3h, 3j and 3l of the main text we show the calculated spectral functions for
fields of 0, 2, 3, 4, 5 and 8 T with the symbols from Figs. 3a, 3c, 3e, 3g, 3i
and 3k overlayed, which shows that the energies and wavevectors of all the
characteristic spectral features are reproduced with quantitative accuracy. In
Fig. S14 we show constant-$\mathbf{Q}$ cuts through the measured and
calculated spectral functions for a range of $\mathbf{Q}$ points and applied
magnetic fields. This extensive comparison reveals that the MPS calculations
lose their quantitative accuracy in locating the excitation features on
approaching the low-field limit ($B=0$ and 1 T), as noted in the main text,
but do an excellent job for all $\mathbf{Q}$ points at $B\geq 2$ T. Regarding
the calculated intensities, we find that the MPS calculations achieve
quantitative accuracy at low $|\mathbf{Q}|$ and between K and M, but clearly
experience some challenges with the low-energy features near the K and M
points.

##  S7 Bound States of Spin Waves in the UUD Phase

With the goal of understanding the excitation spectrum in the UUD phase of the
$S=1/2$ TLAF, we focus on excitations at the same total $S^{z}$ as the 1/3
magnetization plateau. In order to elucidate the spectral response at the
Heisenberg point, we consider an ad-hoc XXZ Hamiltonian with parallel
anisotropy and magnetic field axes:

$\displaystyle\mathcal{H}$ $\displaystyle=$ $\displaystyle\sum_{\langle
i,j\rangle}{\textstyle\frac{1}{2}}J_{xy}\left(S^{+}_{i}S_{j}^{-}+S^{-}_{i}S_{j}^{+}\right)+J_{z}\,S^{z}_{i}S^{z}_{j}$
(S9) $\displaystyle-h\sum_{i}S^{z}_{i},$

where $h=\mu_{B}g_{ab}B$. At $m=1/3$ the ground state in the Ising limit,
$J_{z}/J_{xy}\gg 1$ is the UUD product state illustrated in Fig. 4c of the
main text. An excitation at the same $S^{z}$ consists at least of one pair of
local spin-flips, $U\rightarrow D$ and $D\rightarrow U$. If the two flipped
spins are spatially well separated, they cost an energy $3\times J_{z}$ in the
Ising limit, but if they are nearest neighbors then they cost only $2\times
J_{z}$ in the same limit. We note that, because the total $S^{z}$ is unchanged
by two opposing spin flips, the $z$-axis magnetic field has no effect on the
energies of these two-spin excitations.

As the next step we consider the perturbative limit with a finite, but small,
$J_{xy}$. The nearest-neighbour bound states now show a spatially localized
yet partially mobile structure. A schematic representation is shown in Fig. 4c
of the main text, where the flipped red spin is confined, but can hop around
the flipped blue spin along a finite chain of six sites, i.e. on the hexagon
around the blue spin. The discrete energy spectrum of this six-site chain
shows a lifting of the bound-state degeneracy to 6 flat bands with energies
$\\{2J_{z}+J_{xy},2J_{z}+J_{xy}/2,2J_{z}-J_{xy}/2,2J_{z}-J_{xy}\\}$, with the
second and third energies each twofold degenerate. The manifold of excitations
arising from pairs of flipped spins far away from each other is also split at
first order in $J_{xy}$, and forms a two-particle continuum at an energy
around $3J$ whose support can be derived from the dispersion of the single
spin-flips.

Figure S15: Cylinder MPS calculations of the longitudinal spectral function
for the nearest-neighbour $S=1/2$ TLAF with XXZ spin interctions
$J_{z}/J_{xy}=100$ (a), 5 (b), 2(c) and 1 (d). Panels (e-h) show spectra on a
logarithmic intensity scale. One can see how the faint two magnon scattering
continuum starts to overlap with the bound states as $J_{z}$ is reduced.

The next step is to determine whether this structure of bound states is
visible in the longitudinal dynamical structure factor,
$S_{zz}(\mathbf{Q},\omega)$. For this we performed MPS calculations for the
nearest-neighbour model at different values of $J_{z}/J_{xy}$, and in Fig.
S15(a) we show the spectral function at $J_{z}/J_{xy}=100$ along the usual
$\Gamma$-K-M path in the Brillouin zone. Indeed we find that the spectral
weight becomes $\mathbf{Q}$-dependent, and that the excitation energies are
confined to the perturbative branches (horizontal dashed red lines) discussed
above. As $J_{z}/J_{xy}$ is lowered to $5$ [Fig. S15(b)], $2$ [Fig. S15(c)]
and the Heisenberg point [Fig. S15(d)], one observes that the bound-state
branches acquire an increasing dispersion and that they also broaden in
energy; in particular, while the lower bound state remains relatively sharp as
$J_{z}/J_{xy}\rightarrow 1$, the upper scattering states broaden into the bow-
tie structure of continuum IV.

The origin of this broadening can be seen in the lower part of Fig. S15, where
the same four spectral functions are shown on a logarithmic intensity scale.
In Fig. S15(e) we observe the appearance of the two-particle scattering
continuum centred at an energy of approximately $15J_{xy}$. As $J_{z}/J_{xy}$
is lowered, the two-particle continuum overlaps more and more with the energy
window of the bound states that were well separated in our perturbative
treatment near the Ising point. On proceeding to the Heisenberg point [Fig.
S15(h)], it appears that only a part of the lowest bound-state branch remains
as the well defined and isolated mode observed in experiment (Fig. 4a of the
main text).

## References

  * Rodríguez-Carvajal (1993) J. Rodríguez-Carvajal, Physica B: Condens. Matter 192, 55 (1993).
  * Xing  _et al._ (2020) J. Xing, L. D. Sanjeewa, J. Kim, G. R. Stewart, M.-H. Du, F. A. Reboredo, R. Custelcean, and A. S. Sefat,  (2020), 10.25505/fiz.icsd.cc23j954.
  * Scheie  _et al._ (2021) A. O. Scheie, E. A. Ghioldi, J. Xing, J. A. M. Paddison, N. E. Sherman, M. Dupont, D. Abernathy, D. M. Pajerowski, S.-S. Zhang, L. O. Manuel, A. E. Trumper, C. D. Pemmaraju, A. S. Sefat, D. S. Parker, T. P. Devereaux, J. E. Moore, C. D. Batista, and D. A. Tennant, “Witnessing quantum criticality and entanglement in the triangular antiferromagnet KYbSe2,” (2021), arXiv:2109.11527 .
  * Dai  _et al._ (2021) P.-L. Dai, G. Zhang, Y. Xie, C. Duan, Y. Gao, Z. Zhu, E. Feng, Z. Tao, C.-L. Huang, H. Cao, A. Podlesnyak, G. E. Granroth, M. S. Everett, J. C. Neuefeind, D. Voneshen, S. Wang, G. Tan, E. Morosan, X. Wang, H.-Q. Lin, L. Shu, G. Chen, Y. Guo, X. Lu, and P. Dai, Phys. Rev. X 11, 021044 (2021).
  * Kaneko  _et al._ (2014) R. Kaneko, S. Morita, and M. Imada, J. Phys. Soc. Jpn. 83, 093707 (2014).
  * Li  _et al._ (2015) P. H. Y. Li, R. F. Bishop, and C. E. Campbell, Phys. Rev. B 91, 014426 (2015).
  * 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).
  * Iqbal  _et al._ (2016) Y. Iqbal, W.-J. Hu, R. Thomale, D. Poilblanc, and F. Becca, Phys. Rev. B 93, 144411 (2016).
  * Hu  _et al._ (2019) S. Hu, S. Eggert, and Y.-C. He, Phys. Rev. Lett. 123, 207203 (2019).
  * Xing  _et al._ (2021) J. Xing, L. D. Sanjeewa, A. F. May, and A. S. Sefat, APL Materials 9, 111104 (2021).
  * Ranjith  _et al._ (2019a) K. M. Ranjith, D. Dmytriieva, S. Khim, J. Sichelschmidt, S. Luther, D. Ehlers, H. Yasuoka, J. Wosnitza, A. A. Tsirlin, H. Kühne, and M. Baenitz, Phys. Rev. B 99, 180401 (2019a).
  * Ranjith  _et al._ (2019b) K. M. Ranjith, S. Luther, T. Reimann, B. Schmidt, P. Schlender, J. Sichelschmidt, H. Yasuoka, A. M. Strydom, Y. Skourski, J. Wosnitza, H. Kühne, T. Doert, and M. Baenitz, Phys. Rev. B 100, 224417 (2019b).
  * 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, and A. A. Tsirlin, Phys. Rev. B 100, 144432 (2019).
  * Young  _et al._ (2013) O. Young, A. R. Wildes, P. Manuel, B. Ouladdiaf, D. D. Khalyavin, G. Balakrishnan, and O. A. Petrenko, Phys. Rev. B 88, 024411 (2013).
  * Xing  _et al._ (2019) J. Xing, L. D. Sanjeewa, J. Kim, G. R. Stewart, A. Podlesnyak, and A. S. Sefat, Phys. Rev. B 100, 220407 (2019).
  * Toth and Lake (2015) S. Toth and B. Lake, J. Phys. Condens. Matter 27, 166002 (2015).
  * Chubokov and Golosov (1991) A. V. Chubokov and D. I. Golosov, J. Phys. Condens. Matter 3, 69 (1991).
  * Hotta  _et al._ (2013) C. Hotta, S. Nishimoto, and N. Shibata, Phys. Rev. B 87, 115128 (2013).
  * Hauschild and Pollmann (2018) J. Hauschild and F. Pollmann, SciPost Phys. Lect. Notes , 5 (2018).
  * Szasz  _et al._ (2020) A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore, Phys. Rev. X 10, 021042 (2020).