# Experimentally tunable QED in dipolar-octupolar quantum spin ice

Han Yan (闫寒)  Institute for Solid State Physics, University of Tokyo, Kashiwa,
277-8581 Chiba, Japan  Department of Physics and Astronomy, Rice University,
Houston, TX 77005, USA  Smalley-Curl Institute, Rice University, Houston, TX
77005, USA  Alaric Sanders  TCM Group, Cavendish Laboratory, University of
Cambridge, Cambridge CB3 0HE, United Kingdom  Claudio Castelnovo  TCM Group,
Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United
Kingdom  Andriy H. Nevidomskyy  Department of Physics and Astronomy, Rice
University, Houston, TX 77005, USA  Rice Center for Quantum Materials, Rice
University, Houston, TX 77005, USA

###### Abstract

We propose a readily achievable experimental setting where an external
magnetic field is used to tune the emergent quantum electrodynamics (eQED) of
dipolar-octupolar quantum spin ice (DO-QSI). When the octupolar interactions
dominate, we find that a weak field tunes the DO-QSI’s quantum dynamics
without lifting the degeneracy of the low-energy ice-state manifold. This is a
unique advantage compared to conventional quantum spin ice. In $U(1)_{\pi}$
DO-QSI – the proposed ground state of QSI candidates Ce2Zr2O7, Ce2Sn2O7 and
Ce2Hf2O7 – the field tunes the emergent speed of light (and, consequently, the
emergent fine structure constant). Depending on the field’s alignment with the
crystal, one may induce different speeds of light for the two polarizations of
the emergent photons, in a fascinating analogue of the electro-optic Kerr
effect. In $U(1)_{0}$ DO-QSI – yet to be uncovered experimentally – we find a
number of unusual field-induced transitions, including a transition between
$0$\- and $\pi$-flux QSI phases, as well as phases with frustrated flux
configurations. Our proposal opens the gate to a plethora of experimentally
accessible, engineerable eQED phenomena in the emergent universes of quantum
spin ice.

Introduction — Quantum spin liquids (QSL) – exotic phases of matter where
strong quantum fluctuations prevent long-range order down to absolute zero
temperature – are some of the most exquisitely complex and rich behaviours
observed in quantum many-body physics [1]. Recent years have witnessed an
impressive experimental effort that brought to the fore several candidate
materials. While definitive experimental confirmation of a QSL state has
remained elusive thus far, concerted theoretical and experimental progress has
brought us ever so close to their detection.

Gauge field theory has emerged as a fundamental unifying principle to describe
spin liquid phases. A case in point is quantum spin ice (QSI) [2, 3], whose
behaviour is most naturally understood in terms of a $3+1D$ compact $U(1)$
emergent quantum electrodynamics (eQED) picture [4], featuring gapless
photons, gapped electrically charged matter (spinons) and Dirac monopoles
(visons) as quasiparticles. Various theoretical and numerical studies [5, 6,
7, 8] have inspected different aspects of QSI, including photons and their
speed of light [9, 10], spinon dynamics [11, 12, 13, 14] and Cerenkov
radiation [15], and vison properties [16, 17, 18, 19, 20], many of which go
far beyond the QED accessible in fundamental particle physics.

The recent theoretical prediction [21] of dipolar-octupolar quantum spin ice
(DO-QSI) and experimental discovery of candidate materials Ce2Zr2O7 [22, 23],
Ce2Hf2O7 [24, 25], and Ce2Sn2O7 [26] has further expanded the playground for
QSI physics. In these systems, the ground state doublets of the crystalline
electric fields of the magnetic ions are dominated by the
$\ket{J=5/2,m_{J}=\pm 3/2}$ wave-functions. As a consequence, the effective
pseudospin-$1/2$ operators $S_{i}^{\alpha}$ ($\alpha\in\\{x,y,z\\}$) do not
transform as conventional dipole moments under lattice and time-reversal
symmetries [21], leading to important consequences for their spectroscopic
signatures [27, 28, 11]. Most notably, the octupolar nature of the $S^{y}_{i}$
operators means that they do not couple, to linear order, to an applied
magnetic field.

In this Letter, we show how the properties of the DO-QSI offer a unique
opportunity to tune the parameters of the emergent QED via straightforward
application of external magnetic fields. More specifically, the applied fields
enable us to change the refractive index of photons, leading to the magnetic
analogue of the electro-optic Kerr effect. This process also tunes the speed
of emergent photons and therefore, the fine-structure constant. The effect is
also generically anisotropic, resulting in birefringence, and is directly
applicable to Ce2Zr2O7, Ce2Sn2O7 and Ce2Hf2O7, which are believed to reside in
the $\pi$-flux QSI phase. For systems in the $0$-flux ($U(1)_{0}$) QSI phase,
we further find that the application of external magnetic fields results in
novel phase transitions into the $\pi$-flux ($U(1)_{\pi}$) and various mixed-
flux phases.

Crucially, the dipolar-octupolar nature of the DO-QSI prevents the magnetic
field from directly polarizing the pseudo-spins, making the control of its
dynamics significantly easier than in conventional QSIs. The effect we predict
is achievable using readily available table-top (sub-Tesla) magnetic fields,
in contrast to the previously proposed electric field tuning that requires
very large ($\sim 10^{4}$ V/mm) field strength for similar effects to be
observed [29, 30].

The ability to precisely tune the dynamics of eQED using modest ($\lesssim 1$
T) magnetic fields is an important experimental handle that gives us access to
a vast range of eQED dynamics and phases without the need for challenging
(oftentimes impossible) fine tuning of material parameters. This will
drastically improve our efficiency in creating such new topological states of
matter and studying their properties.

Figure 1:  Schematic illustration of the pyrochlore lattice, and example of
the perturbative generation of the hexagon-flip term, via both the spin-spin
interaction and spin-field coupling. Notice that in our work, pertinent to DO-
QSI candidate materials, the arrows represent the local $S^{y}$ components of
the spins (plotted as pointing ‘into’ or ‘out of’ each tetrahedron for ease of
visualisation).

The model — We consider DO-QSI on the pyrochlore lattice (Fig. 1), described
by a Hamiltonian with nearest neighbor interactions and coupling to magnetic
field [31]

$\begin{split}\mathcal{H}=&\sum_{\langle
ij\rangle}\left[\sum_{\alpha=x,y,z}J_{\alpha}S_{i}^{\alpha}S_{j}^{\alpha}+J_{xz}(S_{i}^{x}S_{j}^{z}+S_{i}^{z}S_{j}^{x})\right]\\\
&-\sum_{i}(S_{i}^{z}\tilde{g}_{zz}+S_{i}^{x}\tilde{g}_{xz}){\bf\it\hat{z}}_{i}\cdot{\bf\it
B}\,.\end{split}$ (1)

Here, the spin operators are in the local basis, and ${\bf\it\hat{z}}_{i}$ is
the unit vector of the local $z$-axis pointing towards/away from the
tetrahedron center. An important feature of this Hamiltonian, as mentioned
earlier, is that the $S^{y}_{i}$ operators do not couple to the external
magnetic field to linear order, due to their octupolar symmetry [21].

We study the effect of tuning the magnetic field ${\bf\it B}$ in DO-QSI with
dominant octupolar interaction $J_{y}$. For simplicity, we take $J_{\pm}\equiv
J_{z}=J_{x}\ll J_{y}$, $J_{xz}=0=\tilde{g}_{xz}$, and $\tilde{g}_{zz}=1$.
Varying these parameters may change the details of the phase diagram but not
the qualitative physics. The Hamiltonian is then more conveniently written as
a dominant part plus a perturbation,
$\mathcal{H}=\mathcal{H}_{I}+\mathcal{H}^{\prime}$, where

$\displaystyle\mathcal{H}_{I}$ $\displaystyle=\sum_{\langle
ij\rangle}J_{y}\sigma_{i}^{z}\sigma_{j}^{z},$ (2)
$\displaystyle\mathcal{H}^{\prime}$ $\displaystyle=\sum_{\langle
ij\rangle}\frac{J_{\pm}}{2}(\sigma_{i}^{+}\sigma_{j}^{-}+\sigma_{i}^{-}\sigma_{j}^{+})-\sum_{i}\frac{{\bf\it\hat{z}}_{i}\cdot{\bf\it
B}}{2}(\sigma_{i}^{+}+\sigma_{i}^{-})\,.$

Here, we relabeled $S^{y}=\sigma_{i}^{z}$, $S^{z}\pm iS^{x}=\sigma^{\pm}_{i}$,
to bring the Hamiltonian to a more customary notation.

The low energy Hilbert space is approximately spanned by the ice states, i.e.,
product states that satisfy the 2-in-2-out rules in the $\sigma_{i}^{z}$
basis, which are ground states of $\mathcal{H}_{I}$. In the presence of the
non-commuting $\mathcal{H}^{\prime}$, the ground states and low-energy
excitations (i.e., photons and visons of the eQED) are coherent quantum
superpositions of the ice states.

The aforementioned absence of coupling to the $\sigma^{z}_{i}=S^{y}_{i}$
pseudo-spin components means that applied magnetic field action is purely
transverse, without biasing a preferred polarisation of the ice states. One
therefore expects the field to tune the DO-QSI phase but not altogether
disrupt it up to $|{\bf\it B}|\sim J^{y}$ (i.e., well in excess of the one
Tesla scale) [27].

The low-energy dynamics is determined by matrix elements between ice states
generated perturbatively. Note that individual terms in $\mathcal{H}^{\prime}$
take the system out of the ice state manifold by creating pairs of spinons. It
is only at higher orders in perturbation theory in $\mathcal{H}^{\prime}$ –
creating, moving and annihilating spinons around closed, non-self-retracing
paths – that one finds quantum processes connecting different ice states. The
simplest such process takes place around an elementary hexagonal loop on the
lattice, illustrated in Fig. 1 and embodied by a perturbative Hamiltonian term
operating within the ice state sector,

$\mathcal{H}_{\text{eff}}=\sum_{p\in\\{\text{all
}\hexagon\\}}g_{p}\left[\mathcal{O}_{p}+\mathcal{O}_{p}^{\dagger}\right]+\ldots\,,$
(3)

where
$\mathcal{O}_{p}=\sigma_{1}^{+}\sigma_{2}^{-}\sigma_{3}^{+}\sigma_{4}^{-}\sigma_{5}^{+}\sigma_{6}^{-}$
is the hexagon-flip operator, and $g_{p}$ is a function of $J_{\pm}$ and
${\bf\it B}$. In the following, we shall consider only these simplest (single
hexagon) terms and ignore higher order processes.

There are four types of hexagons facing different directions. They generally
have different $g_{p}$ values due to the non-uniform magnetic field coupling
proportional to ${\bf\it\hat{z}}_{i}\cdot{\bf\it B}$ on different sublattice
sites. To leading order, we find [32]

$g_{p}=\frac{3J_{\pm}^{3}}{2J_{y}^{2}}+\frac{5J_{\pm}^{2}}{4J_{y}^{3}}({\bf\it\hat{n}}_{p}\cdot{\bf\it
B})^{2}\,,$ (4)

where ${\bf\it\hat{n}}_{p}$ is the unit vector orthogonal to the hexagon $p$.

We will often discuss the four hexagons surrounding the 4-tetrahedron object
shown in Fig. 1. Notice that each hexagon contains $3$ of the $4$ sublattice
sites; therefore we can conveniently label the hexagon types by the sublattice
index $\mu\in\\{0,1,2,3\\}$ that they do not contain. For example, Fig. 1
highlights a hexagon containing sublattice sites $0,1,2$, and we label it
being of type $\mu=3$.

In order to understand the role of this coupling constant in the Hamiltonian
$\mathcal{H}_{\text{eff}}$ in Eq. (3), it is convenient to use the Villain
expansion [33],

$\sigma^{+}_{j}=e^{ia_{j}/2}\sqrt{(S+1/2)^{2}-(\sigma_{j}^{z})^{2}}\;e^{ia_{j}/2}\,,$
(5)

which recreates the spin algebra using the $U(1)$ operators $a_{i}$
canonically conjugate to $\sigma^{z}_{i}$,
$[a_{i},\sigma_{j}^{z}]=i\delta_{ij}$. If we take the limit $S\to\infty$, we
recognize that $\mathcal{O}_{p}\simeq(S+1/2)^{6}e^{i\Phi_{p}}$, to leading
order, where $\Phi_{p}=\nabla_{p}\times a$ are the gauge fluxes of the
emergent QED magnetic field through the hexagon [4, 34]. Eq. (3) thus becomes
$\mathcal{H}_{\text{eff}}\propto\sum_{\mu}g_{\mu}\cos(\Phi_{\mu})$. This
highlights the importance of the sign of $J_{\pm}$ in Eq. (4). If $J_{\pm}>0$,
we have $g>0$ regardless of $|{\bf\it B}|$, and the $U(1)_{\pi}$ state with
$\Phi_{p}\equiv\pi$ remains always lowest in energy. If $J_{\pm}<0$, however,
a competition arises between $J_{\pm}$ and $|{\bf\it B}|$ capable of changing
the sign of $g$ independently on the four sublattices, depending on the
direction of the applied field. In this case, the ground state flux
configuration must be determined by minimizing the semiclassical energy
$2\sum_{\mu=0}^{3}g_{\mu}\cos(\Phi_{\mu})$ over $\Phi_{\mu}\in[-\pi,\pi]^{4}$
subjected to the quantization condition $\sum_{\mu=0}^{3}\Phi_{\mu}=0$
($\mathrm{mod}\;2\pi$) around the four-tetrahedron structure shown in Fig. 1.
We plot the resulting phase diagram in Fig. 2.

Figure 2:  Panels (a-c): $B-J_{\pm}$ phase diagrams for [111], [110] and [100]
applied fields. Zero-flux $U(1)_{0}$ and $\pi$-flux $U(1)_{\pi}$ states
compete with the mixed-flux $\pi\pi 00$ and a continuously-varying frustrated
flux phase, see Fig. 4 for further details. Phase transitions are shown as
solid black lines. Panels (d-h): raw effective hexagon-flip operator
coefficient strength from the leading-order expansion, Eq. (4). The color
scale is in units of $J_{y}$. Note that $J_{\pm}=0$ marks the line in which
all standard QSI hexagon-flips vanish, i.e., the classical limit.

This procedure yields four possible ground states: the standard $U(1)_{\pi}$
and $U(1)_{0}$ states, in addition to a mixed-flux $\pi\pi 00$ state (in which
precisely two hexagon sublattices carry flux $\pi$, and two have vanishing
flux); and a continuously variable, non-quantized flux arrangement which we
dub ‘frustrated flux’ (FF) state. Photon dynamics corresponds to quadratic
fluctuations about these ground states.

Figure 3:  Numerical (color map) and analytical (red traces) dynamical
$S^{y}S^{y}$ correlations in the $\pi\pi 00$, $U(1)_{0}$, $U(1)_{\pi}$,
$g_{\mu}=(\pm 1/2,\pm 1/2,[\pm]1,[\pm]1)$ phases (left column); and in the
frustrated-flux, $g_{\mu}=(-2,1,1,1)$ phase (right column). As discussed in
the main text, the photon dispersion is not affected by changing pairs of
signs of $g$, hence the identical correlations in all three phases in the left
panels. White arrows in panel d) mark shoulders in the spectral weight due to
three-photon up-conversion processes, which do not occur in the non-frustrated
phases on the left. Both $T$ and $g_{\mu}$ are in units of
$J_{\pm}^{3}/J_{y}^{2}$.

Photon dispersion and emergent Kerr effect — We calculated the spectrum of the
emergent photon by two methods. First, we performed a semiclassical path
integral simulation using the method of Ref. 16; second, we calculated the
photon dispersion using the large-$S$ formalism of Refs. 9, 34, 17. The
details of these calculations may be found in the Supplemental Material [32].
We see from Fig. 3 that the strongest features in the simulations are
precisely the modes of the quadratic theory.

First, we remark on the rich physics achievable in the $J_{\pm}>0$ case
relevant to the known cerium pyrochlore materials. Generically, magnetic field
couples differently to the four types of hexagons, thus breaking the
crystalline symmetry. As a result, the degeneracy of the photon branches is
lifted or, in the language of optics, the medium becomes birefringent, as seen
in Fig. 3. Near the $\Gamma$ point, the splitting – and therefore the change
in speed of light $c$ of eQED for different polarizations – scales as $\Delta
c\propto B^{2}/J_{y}$, in a fascinating emergent analogue of the electro-
optical Kerr effect (except with the magnetic, rather than electric field
being applied here). For instance, using the parameters of Ce2Hf2O7 given in
Ref. 25, $\Delta c/c\sim 1$ can be achieved at the easily accessible field of
$\sim 0.3$ T.

In the special case of ${\bf\it B}\parallel$ [100], all hexagon types have the
same value of $g_{\mu}$: $g_{\mu}\propto(1,1,1,1)$. In this case, no
birefringence ensues but the speed of light gets altered. This in turn affects
the fine structure constant $\alpha=q^{2}/(\hbar c)$ ($q$ being the spinon
charge), illustrating the fascinating level of control achievable with even
modest magnetic fields.

Tuning eQED phases with applied fields — We next consider the case
$J_{\pm}<0$, in which the QSI is in the $U(1)_{0}$ phase in the absence of
magnetic field. In addition to the birefringent physics discussed above, the
applied magnetic field now also tunes the sign of $g_{p}$, allowing one to
change the ground state flux configuration from $U(1)_{0}$ to $\pi-$flux,
$\pi\pi 00$ flux or FF, as shown in Fig. 4. The photon dynamics for the
$\pi\pi 00$ and $\pi$-flux phases, depicted in Fig. 3, are equivalent to the
standard 0-flux $U(1)_{0}$ spin liquid (see SM [32] for details). However,
this is not the case for the FF phase, which does not exhibit quantized flux,
as can be seen in Fig. 4 – this is a crucial distinction from the other three
phases, as it reflects frustration of the emergent magnetic fluxes that gives
rise to new behaviours.

Figure 4:  Field direction dependence of the ground state flux phase, in the
interesting case $\|{\bf\it B}\|=\sqrt{4.2J_{y}J_{\pm}}$, where the effect of
the magnetic field is comparable to the intrinsic hexagon-flip terms. Top
panel: Phase diagram at constant magnetic field, labeled by high-symmetry
directions on the unit sphere. Bottom panel: variation of the four fluxes as
the magnetic field is rotated (namely, it is moved along the green path in the
top panel).

At higher temperatures, the nonlinearities of the system are exposed by the
stronger fluctuations. Working in terms of the emergent large-$S$ magnetic
flux $\Phi$, the hexagon-flip operator $\mathcal{O}\sim e^{i\Phi}$, and we see
two contributions from the order of $\Phi^{4}$ terms: (i) a Hartree-Fock
correction to the dispersion, and (ii) a lifetime broadening of the photon
modes. More subtly, we see a high temperature effect in the FF phase (Fig. 3d)
that is absent in the non-frustrated phases (panel c): the appearance of a
structured continuum of excitations above the photon bandwidth. In the
language of nonlinear optics, this behaviour is referred to as the three-
photon up-conversion (see SM [32] for details).

Magnetic field tuning in DO-QSI vs QSI — We emphasize the differences between
DO-QSI and QSI in regards to the rich physics discussed above. In conventional
QSI, applied magnetic fields unavoidably couple also to the component of the
spins subjected to the ice rules 111Note that, due to the 3D non-collinear
nature of the local easy axes, it is not possible to apply a purely transverse
external field. and favour a polarised, classical ice-state as a result. This
selection is only countered by the kinetic terms in the Hamiltonian, of the
order of $g|_{{\bf\it B}=0}\sim J_{\pm}^{3}/J_{z}^{2}$ and much smaller than
the dominant exchange coupling. One would therefore be limited in the use of
applied fields to tune the eQED behaviour to very small field strengths. While
some precursors of the behaviours discussed in our work may be achievable in
conventional QSI, it is likely that they would remain very subtle; most likely
one would not be able to observe the rich phase diagram and, for instance, the
exotic physics of the FF phase at large fields ($|{\bf\it
B}|\gtrsim\sqrt{|J_{\pm}|J_{z}}\gg g|_{{\bf\it B}=0}$).

In contrast, DO-QSI is free from such a problem, since the field does not
couple to the dominant octupolar component $\sigma^{z}$ subjected to the ice-
rule, and acts as a purely transverse field. Hence, a phase transition out of
the eQED phase is only expected for $|{\bf\it B}|\gtrsim J_{y}$ [27].

This is crucial for its ability to access the predicted tunable eQED regime in
experiments. Among the DO-QSI candidates uncovered to date (summarised for
convenience in Table 1 in the SM [32]), it has been suggested that the
compounds $\mathrm{Ce}_{2}\mathrm{Hf}_{2}\mathrm{O}_{7}$ and
$\mathrm{Ce}_{2}\mathrm{Sn}_{2}\mathrm{O}_{7}$ may have dominant $J_{y}$
interactions, while $\mathrm{Ce}_{2}\mathrm{Zr}_{2}\mathrm{O}_{7}$ has
comparable $J_{x}$ and $J_{y}$. Our proposal is thus directly applicable to
the former two materials, for which single crystals exist.

Discussion — In this work, we demonstrated that various interesting eQED
physics in DO-QSI systems can be accurately tuned using experimentally
accessible applied magnetic fields. Modest field strengths
($\lesssim\sqrt{|J_{\pm}|J_{y}}$) generically cause sizable, anisotropic
changes to the photon dispersion due to an emergent analog of the electro-
optic Kerr effect. Strong fields in the $U(1)_{0}$ phase of DO-QSI have the
even more dramatic effect of driving the system into eQED phases with
different background flux configurations. The frustrated-flux phase is
particularly unusual, as it allows for the presence of nonlinear three-photon
processes that qualitatively alter the photon dispersion. Experimentally, our
proposal can be straightforwardly applied to DO-QSI candidate materials with
dominant $J_{y}$ interactions. The by-design emergent-photon dynamics also
opens the gate to other possible interesting photonic phenomena including
lensing effects and wave guides based on total internal reflection.

While this work focuses on the photon sector of the physics, many other
important aspects of the tunable eQED remain to be studied. One of them is the
spinon spectrum at energy scale $J_{y}$ [36, 11, 28, 19, 27]. It is already
known that the mean-field spinon bands can show signatures of the underlying
flux structure. In the $U(1)_{\pi}$ and $\pi\pi 00$ phases, the spinon unit
cell is effectively doubled [28, 36, 27]. Spinon bands in the FF phase are
believed to have finite Chern number, manifested in the thermal Hall effect
[37]. Furthermore, the spinon spectrum is easier to measure in neutron
scattering at small wave-vectors, so characterizing the spectrum in different
phases is crucial to identifying eQED phase transitions.

In addition, the Zeeman coupling ${\bf\it\hat{z}}_{i}\cdot{\bf\it
B}\sigma^{x}_{i}$ contributes to spinon propogation with the subtle difference
that it hops spinons between nearest-neighbour tetrahedra, instead of the next
nearest neighbour ones (as due to the $J_{\pm}\sigma^{+}\sigma^{-}$ terms).
The latter preserve the spinon number on the ‘up’ and ‘down’ tetrahedral
sublattice separately, while the former couples the two sublattices and lifts
the spinon band degeneracy [27]. This suggests that by applying a magnetic
field, one can introduce net spinon charge imbalances between the two
tetrahedral sublattices. Later quenching the field to zero, one can abruptly
decouple the populations of ‘up’ and ‘down’ spinons, locking in the charge
imbalance. Such a procedure may be used to generate QSI states hosting spinons
and antispinons (on different tetrahedral sublattices) that are unable to
annihilate and therefore survive down to vanishing temperatures, given proper
boundary conditions.

Another interesting perspective comes from the eQED field theory discussed in
Refs. 10, 38, 15. With the photon speed controlled by the applied magnetic
field, we expect the fine-structure constant to be tunable too. This provides
an exciting avenue to probe both weakly and strongly coupled eQED,
ultrarelativistic effects, and transitions between different regimes. Emergent
light by-design may soon become a reality in DO-QSI, presenting us with a new
landscape of QED beyond elementary particle physics to be explored.

###### Acknowledgements.

The authors would like to thank Gang Chen and Yamashita Minoru for helpful
discussions. This work was funded in part by the Engineering and Physical
Sciences Research Council (EPSRC) grants No. EP/T028580/1 and No. EP/V062654/1
(C.C.). A.H.N. and H.Y. are grateful for the hospitality of the Aspen Center
for Physics, supported by the U.S. National Science Foundation grant
PHY-2210452.

## References

  * Moessner and Moore [2021] R. Moessner and J. E. Moore,  _Topological Phases of Matter_ (Cambridge University Press, 2021).
  * Gingras and McClarty [2014] M. J. P. Gingras and P. A. McClarty, Reports on Progress in Physics 77, 056501 (2014).
  * Udagawa and Jaubert [2021] M. Udagawa and L. Jaubert, eds.,  _Spin Ice_, Springer Series in Solid-State Sciences, Vol. 197 (Springer International Publishing, 2021).
  * Hermele  _et al._ [2004] M. Hermele, M. P. A. Fisher, and L. Balents, Phys. Rev. B 69, 064404 (2004).
  * Huang  _et al._ [2018] C.-J. Huang, Y. Deng, Y. Wan, and Z. Y. Meng, Physical Review Letters 120, 167202 (2018).
  * Shannon  _et al._ [2012] N. Shannon, O. Sikora, F. Pollmann, K. Penc, and P. Fulde, Physical Review Letters 108, 067204 (2012).
  * Udagawa and Moessner [2019] M. Udagawa and R. Moessner, Physical Review Letters 122, 117201 (2019).
  * Hosoi  _et al._ [2022] M. Hosoi, E. Z. Zhang, A. S. Patri, and Y. B. Kim, Physical Review Letters 129, 097202 (2022).
  * Benton  _et al._ [2012] O. Benton, O. Sikora, and N. Shannon, Phys. Rev. B 86, 075154 (2012).
  * Pace  _et al._ [2021] S. D. Pace, S. C. Morampudi, R. Moessner, and C. R. Laumann, Physical Review Letters 127, 117205 (2021).
  * Chen [2023] G. Chen, Physical Review Research 5, 033169 (2023).
  * Desrochers and Kim [2023] F. Desrochers and Y. B. Kim, Spectroscopic signatures of fractionalization in octupolar quantum spin ice (2023).
  * Lee  _et al._ [2012] S. Lee, S. Onoda, and L. Balents, Physical Review B 86, 104412 (2012).
  * Savary and Balents [2012] L. Savary and L. Balents, Physical Review Letters 108, 037202 (2012).
  * Morampudi  _et al._ [2020] S. C. Morampudi, F. Wilczek, and C. R. Laumann, Physical Review Letters 124, 097204 (2020).
  * Szabó and Castelnovo [2019] A. Szabó and C. Castelnovo, Physical Review B 100, 014417 (2019).
  * Kwasigroch [2020] M. P. Kwasigroch, Physical Review B 102, 125113 (2020).
  * [18] A. Sanders and C. Castelnovo,  arXiv:2310.17076 .
  * Chen [2017a] G. Chen, Physical Review B 96, 195127 (2017a).
  * Bergman  _et al._ [2006] D. L. Bergman, G. A. Fiete, and L. Balents, Physical Review B 73, 134402 (2006).
  * Huang  _et al._ [2014] Y.-P. Huang, G. Chen, and M. Hermele, Phys. Rev. Lett. 112, 167203 (2014).
  * Gao  _et al._ [2019] B. Gao, T. Chen, D. W. Tam, C.-L. Huang, K. Sasmal, D. T. Adroja, F. Ye, H. Cao, G. Sala, M. B. Stone, C. Baines, J. A. T. Verezhak, H. Hu, J.-H. Chung, X. Xu, S.-W. Cheong, M. Nallaiyan, S. Spagna, M. B. Maple, A. H. Nevidomskyy, E. Morosan, G. Chen, and P. Dai, Nature Physics 15, 1052 (2019).
  * Gaudet  _et al._ [2019] J. Gaudet, E. M. Smith, J. Dudemaine, J. Beare, C. R. C. Buhariwalla, N. P. Butch, M. B. Stone, A. I. Kolesnikov, G. Xu, D. R. Yahne, K. A. Ross, C. A. Marjerrison, J. D. Garrett, G. M. Luke, A. D. Bianchi, and B. D. Gaulin, Phys. Rev. Lett. 122, 187201 (2019).
  * Porée  _et al._ [2022] V. Porée, E. Lhotel, S. Petit, A. Krajewska, P. Puphal, A. H. Clark, V. Pomjakushin, H. C. Walker, N. Gauthier, D. J. Gawryluk, and R. Sibille, Physical Review Materials 6, 044406 (2022), 2203.16295 .
  * Porée  _et al._ [2023a] V. Porée, A. Bhardwaj, E. Lhotel, S. Petit, N. Gauthier, H. Yan, V. Pomjakushin, J. Ollivier, J. A. Quilliam, A. H. Nevidomskyy, H. J. Changlani, and R. Sibille, (2023a), arXiv:2305.08261 .
  * Sibille  _et al._ [2020] R. Sibille, N. Gauthier, E. Lhotel, V. Porée, V. Pomjakushin, R. A. Ewings, T. G. Perring, J. Ollivier, A. Wildes, C. Ritter, T. C. Hansen, D. A. Keen, G. J. Nilsen, L. Keller, S. Petit, and T. Fennell, Nature Physics 16, 546 (2020).
  * Li and Chen [2017] Y.-D. Li and G. Chen, Physical Review B 95, 041106 (2017).
  * Yao  _et al._ [2020] X.-P. Yao, Y.-D. Li, and G. Chen, Physical Review Research 2, 013334 (2020).
  * Lantagne-Hurtubise  _et al._ [2017] E. Lantagne-Hurtubise, S. Bhattacharjee, and R. Moessner, Physical Review B 96, 125145 (2017).
  * Mandal [2019] I. Mandal, The European Physical Journal B 92, 187 (2019).
  * Patri  _et al._ [2020] A. S. Patri, M. Hosoi, and Y. B. Kim, Phys. Rev. Research 2, 023253 (2020).
  * [32] See Supplemental Material at PlaceholderURL for details of perturbation theory, numerics and large-$S$ formalism.
  * Villain [1979] J. Villain, Zeitschrift für Physik B Condensed Matter and Quanta 33, 31 (1979).
  * Kwasigroch  _et al._ [2017] M. P. Kwasigroch, B. Douçot, and C. Castelnovo, Physical Review B 95, 134439 (2017).
  * Note [1] Note that, due to the 3D non-collinear nature of the local easy axes, it is not possible to apply a purely transverse external field.
  * Chen [2017b] G. Chen, Physical Review B 96, 085136 (2017b).
  * Yang  _et al._ [2020] H.-J. Yang, H. S. Kim, and S. Lee, Physical Review B 102, 060405 (2020).
  * Laumann and Moessner [2023] C. R. Laumann and R. Moessner,  (2023), arXiv:2302.06635 .
  * Yahne  _et al._ [2022] D. R. Yahne, B. Placke, R. Schäfer, O. Benton, R. Moessner, M. Powell, J. W. Kolis, C. M. Pasco, A. F. May, M. D. Frontzek, E. M. Smith, B. D. Gaulin, S. Calder, and K. A. Ross, (2022), arXiv:arXiv:2211.15140 .
  * Porée  _et al._ [2023b] V. Porée, H. Yan, F. Desrochers, S. Petit, E. Lhotel, M. Appel, J. Ollivier, Y. B. Kim, A. H. Nevidomskyy, and R. Sibille, (2023b), arXiv:2304.05452 .
  * Smith  _et al._ [2023] E. M. Smith, J. Dudemaine, B. Placke, R. Schäfer, D. R. Yahne, T. DeLazzer, A. Fitterman, J. Beare, J. Gaudet, C. R. C. Buhariwalla, A. Podlesnyak, G. Xu, J. P. Clancy, R. Movshovich, G. M. Luke, K. A. Ross, R. Moessner, O. Benton, A. D. Bianchi, and B. D. Gaulin, Phys. Rev. B 108, 054438 (2023).
  * Bhardwaj  _et al._ [2022] A. Bhardwaj, S. Zhang, H. Yan, R. Moessner, A. H. Nevidomskyy, and H. J. Changlani, npj Quantum Materials 7, 51 (2022).

Supplementary Materials for “Experimentally tunable QED in dipolar-octupolar
quantum spin ice”

Han Yan, Alaric Sanders, Claudio Castelnovo, Andriy H. Nevidomskyy

##  I Perturbation Theory

The starting point is the Hamiltonian $\mathcal{H}^{\prime}$ in Eq. (2) of the
main text:

$\displaystyle\mathcal{H}^{\prime}$ $\displaystyle=\sum_{\langle
ij\rangle}\frac{J_{\pm}}{2}(\sigma_{i}^{+}\sigma_{j}^{-}+\sigma_{i}^{-}\sigma_{j}^{+})-\sum_{i}\frac{{\bf\it\hat{z}}_{i}\cdot{\bf\it
B}}{2}(\sigma_{i}^{+}+\sigma_{i}^{-})$ (S1)

By performing a perturbative expansion with respect to $\mathcal{H}^{\prime}$,
relative to the Hamiltonian $\mathcal{H}_{I}=\sum_{\langle
ij\rangle}J_{y}\sigma_{i}^{z}\sigma_{j}^{z}$ enforcing the spin-ice rules, we
obtain an effective Hamiltonian of the general form

$\mathcal{H}_{\text{eff}}=(1-\mathcal{P})\bigg{[}-\mathcal{H}^{\prime}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}^{\prime}+\mathcal{H}^{\prime}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}^{\prime}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}^{\prime}+\dots\bigg{]}(1-\mathcal{P})\,,$
(S2)

where $1-\mathcal{P}$ is the projection operator onto the ice states. The only
non-vanishing terms in the above expansion are therefore those that map spin
ice configurations to spin ice configurations, with the affected spins forming
closed loops (or trivial self-retracing paths).

It is useful to split $\mathcal{H}^{\prime}$ into two contributions, the XX
part $\mathcal{H}_{\rm XX}=\sum_{\langle
ij\rangle}\frac{J_{\pm}}{2}(\sigma^{+}_{i}\sigma^{-}_{j}+\sigma^{-}_{i}\sigma^{+}_{j})$
common to all QSI models and the Zeeman term $\mathcal{H}_{\rm
Z}=-\sum_{i}\frac{{\bf\it z}_{i}\cdot{\bf\it B}}{2}\sigma^{x}_{i}$ unique to
DO-QSI.

The ice rule constraint immediately excludes any contribution at first order.
At second order, we have two non-vanishing contributions. From
$\mathcal{H}_{\rm XX}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}_{\rm XX}$,
we obtain only a term proportional to
$J_{\pm}^{2}\\{\sigma_{i}^{+}\sigma_{j}^{-},\sigma_{i}^{-}\sigma_{j}^{+}\\}$
arising from two factors of $\sigma^{+}\sigma^{-}$ occurring on the same
pyrochlore bond. This is a term from standard spin ice [4], which counts the
number of nearest-neighbour, oppositely aligned spins. Within the ice
manifold, each tetrahedron has exactly six antialigned nearest-neighbour
spins, and this term is therefore a trivial constant. The second contribution,
arising from $\mathcal{H}_{\rm
Z}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}_{\rm Z}$, is proportional to
$B^{2}\\{\sigma_{i}^{+},\sigma_{i}^{-}\\}$ and so constant.

At third order, we find non-vanishing terms from (I) $\mathcal{H}_{\rm
XX}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}_{\rm
XX}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}_{\rm XX}$ and (II)
permutations of $\mathcal{H}_{\rm
Z}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}_{\rm
Z}\frac{\mathcal{P}}{\mathcal{H}_{I}}\mathcal{H}_{\rm XX}$. From (I), we
obtain further trivial constant terms and the well known hexagon-flip term
$\frac{3J_{\pm}^{3}}{2J_{y}^{2}}(\mathcal{O}+\mathcal{O}^{\dagger})$, while a
short calculation reveals (II) to be a constant.

The leading order non-trivial effect of the magnetic field comes at fourth
order, from processes similar to those depicted in Fig. 1 in the main text.
Two channels contribute to this process, corresponding to different
arrangements of the field-mediated single-spin flips, see Fig. S1a,b).

Figure S1:  Hexagon-flip processes occurring in the perturbative expansion.
Following the convention of Refs. 4, 18, we represent spins as hard-core
bosons, with down (up) spins written as empty (filled) circles, all in the
$S^{y}$ basis. Gray circles are spectator spins that are decoupled from the
process. The two leading order processes responsible for field-mediated
hexagon-flip terms are show in panels a) and b). Panel c) shows the fourth
order pure-XX correction to the hexagon-flip Hamiltonian, note that spin 7 is
flipped twice. See Ref. 18 for further remarks on the combinatorics associated
with these diagrams.

It may be show from these considerations that the effective Hamiltonian at
order $J_{\pm}^{3}/J_{y}^{2}$ and $B^{2}J_{\pm}^{2}/J_{y}^{3}$ is of the form

$\mathcal{H}_{\text{eff}}=\sum_{p\in\\{\text{all
}\hexagon\\}}g_{p}\left[\mathcal{O}_{p}+\mathcal{O}_{p}^{\dagger}\right]\ .$
(S3)

In practice, we calculate the values of $g_{p}$ on each type of hexagon by six
site exact diagonalisation. Starting from the fourth order series of Eq. (S2),
we project $\mathcal{H}_{I},\ \mathcal{H}^{\prime},$ and $\mathcal{P}$ on to
the basis of all possible hexagon states and so recast them as $2^{6}$ by
$2^{6}$ matrices. The coefficient $g_{p}$ given in Eq. (4) of the main text
then corresponds to the matrix element of $\mathcal{H}_{\rm eff}$ connecting
the state of clockwise pointing spins to the state of counter-clockwise
pointing spins. This calculation includes all the relevant spinon creation,
propagation and annihilation processes implicitly, automatically including all
relevant combinatorial factors.

Other higher order contributions, such as the pure XX contribution at fourth
order shown in Fig. S1c) can be generated when one includes higher order
perturbations and more neighbouring spins, but are neglected in our work since
they do not alter the qualitative features of the phase diagram in Fig. 2.

##  II Semiclassical simulations

We studied the effective Hamiltonian in Eq. (S3) using the semiclassical
method detailed in Refs. 16, 18, which effectively produces a zeroth-order
approximation to the imaginary-time path integral of the system.

In this approach, we simulate finite temperature dynamics by the following
procedure:

  1. 1.

Using a classical Metropolis algorithm, we generate an ensemble of equilibrium
configurations at the desired temperature.

  2. 2.

For each of these configurations, we perform classical Landau–Lifshitz time
evolution.

  3. 3.

Finally, we average any observables over the ensemble of thermal
configurations.

We simulate a periodic cubic supercell of $12^{3}\times 16=27\,648$ spins.
Starting from a temperature of $T=10$ (in units of $J_{\pm}^{3}/J_{y}^{2}$),
we anneal to the final temperature via $256$ logarithmically-spaced steps,
using Monte Carlo acceptance rates to monitor convergence. We average the time
evolution dynamics over an ensemble of $32$ thermal configurations. In all
other respects, the simulation parameters are the same as those used in Ref.
18.

##  III Relevance to DO-QSI candidate materials

In Table 1 we summarise for convenience some relevant details of the cerium-
based pyrochlore materials studied in the literature to date.

Material | $J_{y}$ | $J_{x}$ | $J_{z}$ | $J_{xz}$ | $\begin{matrix}\text{single}\\\ \text{crystal}\end{matrix}$
---|---|---|---|---|---
$\mathrm{Ce}_{2}\mathrm{Hf}_{2}\mathrm{O}_{7}$ [25] | $0.044$ | $0.011$ | $0.016$ | $-0.002$ | Yes 222private communication
$\mathrm{Ce}_{2}\mathrm{Hf}_{2}\mathrm{O}_{7}$ [25] | $0.022$ | $0.046$ | $0.011$ | $-0.001$ | Yes
$\mathrm{Ce}_{2}\mathrm{Sn}_{2}\mathrm{O}_{7}$ [39] | $-0.001$ | $0.045$ | $-0.012$ | $0$ | No
$\mathrm{Ce}_{2}\mathrm{Sn}_{2}\mathrm{O}_{7}$ [26, 40] | $0.069$ | $0.034$ | $0.034$ | $0$ | No
$\mathrm{Ce}_{2}\mathrm{Zr}_{2}\mathrm{O}_{7}$ [41] | $0.063$ | $0.062$ | $0.011$ | $0$ | Yes
$\mathrm{Ce}_{2}\mathrm{Zr}_{2}\mathrm{O}_{7}$ [42] | $0.069$ | $0.068$ | $0.013$ | $0$ | Yes
Table 1:  Estimated coupling parameters for DO-QSI candidates, quoted in meV.

##  IV Large-$S$ formalism

The spin Hamiltonian of Eq. (S3) admits an analytical solution via the Villain
expansion given in Eq. (5) in the main text. To quadratic order in spin
fluctuations, we obtain [34, 9]

$\mathcal{H}_{\mathrm{large-}S}\simeq\sum_{j\in\\{\rm{Spins}\\}}U_{\alpha(j)}(\sigma^{z}_{j}(t))^{2}-\sum_{p\in\\{\hexagon\\}}g_{\mu(p)}\cos(\Phi_{p}(t))\,,$
(S4)

where $\mu(p)\in\\{0,1,2,3\\}$ denotes the type of hexagon,
$\alpha(j)\in\\{0,1,2,3\\}$ denotes the sublattice of spin $j$, and
$U_{\alpha(j)}=2\sum_{\mu\neq\alpha}g_{\mu}$ represents a sum over the six
hexagon-flip coefficients bordering spin $j$.

Noting that in the ice-state manifold, $\langle\sigma^{z}\rangle=0$, we first
statically extremise this Hamiltonian and obtain the zero-temperature flux
configuration $\Phi^{\rm GS}_{p}$. We find that this always results in one of
the four possible ground states: 0-flux, $\pi$-flux, $\pi\pi 00$ flux or
frustrated flux (FF) state. We then expand in fluctuations, setting
$\Phi_{p}(t)=\Phi^{\rm GS}_{p}+\delta\Phi_{p}(t)$ and truncating the series at
quadratic order in $\delta\Phi$. Working in the radiation gauge, $\sigma^{z}$
is identified with $\partial_{t}a_{j}$ and the normal modes found in the usual
way.

In order to provide a concrete expression for the quadratic modes, we are now
forced to change index convention. The spins $j$ of the pyrochlore lattice may
be uniquely indexed in terms of an FCC lattice site ${\bf\it R}_{t}$,
corresponding to the centre of an ‘up’ tetrahedron, and a sublattice index
$\alpha\in\\{0,1,2,3\\}$. We find that the Fourier modes of $\sigma^{z}$ and
$a$,

$\displaystyle a(\boldsymbol{R}_{t}+\boldsymbol{e}_{\rho}/2)$
$\displaystyle=\sqrt{\frac{2}{N_{t}}}\sum_{\boldsymbol{k}}e^{i\boldsymbol{k}\cdot(\boldsymbol{R}_{t}+\boldsymbol{e}_{\rho}/2)}a_{\rho}(\boldsymbol{k})$
$\displaystyle\sigma^{z}(\boldsymbol{R}_{t}+\boldsymbol{e}_{\sigma}/2)$
$\displaystyle=\sqrt{\frac{2}{N_{t}}}\sum_{\boldsymbol{k}}e^{i\boldsymbol{k}\cdot(\boldsymbol{R}_{t}+\boldsymbol{e}_{\sigma}/2)}p_{\sigma}(\boldsymbol{k})\,,$

evolve as

$\displaystyle\partial_{t}a_{\alpha}({\bf\it k},t)$
$\displaystyle=\sum_{\mu\neq\alpha}2g_{\mu}\sigma^{z}_{\alpha}(-{\bf\it k},t)$
$\displaystyle\partial_{t}\sigma^{z}_{\alpha}({\bf\it k},t)$
$\displaystyle=-\sum_{\nu}g_{\nu}\sin({\bf\it h}_{\alpha\nu}\cdot{\bf\it
k})\sin({\bf\it h}_{\beta\nu}\cdot{\bf\it k})a_{\beta}(-{\bf\it k},t)$
$\displaystyle{\bf\it h}_{\mu\nu}$
$\displaystyle=\frac{a_{0}}{\sqrt{8}}\frac{\boldsymbol{e}_{\mu}\times\boldsymbol{e}_{\lambda}}{\|\boldsymbol{e}_{\mu}\times\boldsymbol{e}_{\lambda}\|}\,,$

where $\boldsymbol{e}_{\alpha}$ denotes the vector joining the centre of an
‘up’ tetrahedron to the spin on sublattice $\alpha$, and $a_{0}$ denotes the
cubic lattice parameter.

##  V Photon up-conversion in the FF case

As Fig. 3d in the main text shows, the FF phase is distinct from the non-
frustrated cases in that a structured continuum of excitations with several
step-like features appears. In optics, this phenomenon is known as three-
photon up-conversion. Fundamentally, it arises from the expansion of
$\cos\Phi_{p}$ about its zero-temperature expectation value: the series
expansions of $\cos(0+\delta\Phi)$ and $\cos(\pi+\delta\Phi)$ contain only
even powers of $\delta\Phi$, while no such symmetry applies in the frustrated
flux case. In the large-$S$ formalism, any beyond-quadratic terms in the
expansion are interaction vertices - it is only in the FF phase that the
$\delta\Phi^{3}$ vertex survives. It can readily be shown [18] that pairs of
contracted three-photon vertices have the same order in renormalization group
as the Hartree-Fock corrections in large-$S$ theory, suggesting that these
effects will be of similar importance for spin-$1/2$ moments (note that
Hartree-Fock vertices are believed to correct the speed of light by a factor
of $2.7$ [34]).

##  VI Dualities of the model

The semiclassical analysis we presented amounts to expanding the emergent
magnetic field flux $b(t)=\Phi_{0}+\delta b(t)$ to quadratic order in the
field fluctuations $\delta b(t)$, where $\Phi_{0}$ is the static flux
background as depicted in Figs. 2 and 4 in the main text. This allows us to
map the behaviour in the $U(1)_{\pi}$ and $\pi\pi 00$ phases onto the
$U(1)_{0}$ spin ice phase. Here we show that this mapping also works at the
level of the spin Hamiltonian.

We define the unitary operator

$\displaystyle U[a]=\prod_{j}\exp(ia_{j}S^{y}_{j})\,,$ (S5)

where the index $j$ runs over all spins, and $a_{j}\in[0,2\pi)$. It can be
shown that $U[a]\mathcal{O}_{p}U[a]^{\dagger}=e^{i\nabla_{p}\times
a}\mathcal{O}_{p}$. The lattice curl
$\nabla_{p}\times=a_{1}-a_{2}+a_{3}-a_{4}+a_{5}-a_{6}$ is the alternating sum
of emergent gauge fields $a$. In the large-$S$ limit, one may interpret this
as the statement that emergent magnetic fields may be added together, with the
caveat of $2\pi$ periodicity.

By ‘subtracting $\pi$’ from two of the four sublattices, one can map the
$\pi\pi 00$ phase to $U(1)_{0}$. Given a gauge field configuration $a_{\pi\pi
00}$ generating a $\pi\pi 00$ state, it follows that

$\begin{split}U[a_{\pi\pi
00}]\mathcal{H}_{\text{eff}}[g_{0},g_{1},g_{2},g_{3}]&U[a_{\pi\pi
00}]^{\dagger}\\\
&=\mathcal{H}_{\text{eff}}[-g_{0},-g_{1},g_{2},g_{3}]\,.\end{split}$ (S6)

The gauge behaviour of all anisotropic $g$ models is therefore symmetric with
respect to flipping the signs of any two $g$ values (though one may expect the
symmetries of the spinons to fractionalize due to the enlarged gauge unit
cell).

The frustrated-flux phase, by contrast, is not dual to $U(1)_{0}$, as the
frustration of the emergent fluxes cannot be subtracted away. For
concreteness, let $a_{\rm FF}$ generate the flux configuration
$[\pi,\pi/3,\pi/3,\pi/3]$, as an example. Then, letting $\sigma(p)$ denote the
sublattice of hexagon ${p}$, we find

$\begin{split}&U[a_{\rm FF}]H[g_{0},g_{1},g_{2},g_{3}]U[a_{\rm
FF}]^{\dagger}\\\
&=H\left[-g_{0},\frac{g_{1}}{2},\frac{g_{2}}{2},\frac{g_{3}}{2}\right]+\sum_{p,\sigma(p)\neq
0}ig_{p}\frac{\sqrt{3}}{2}(\mathcal{O}_{p}-\mathcal{O}^{\dagger}_{p})\,.\end{split}$
(S7)

This highly unusual complex hexagon-flip term has previously been identified
in spin-$1/2$ QSI on the breathing pyrochlore lattice [18]. It is associated
with the breaking of a) inversion symmetry and b) charge-conjugation symmetry,
both of which are facilitated here by a [111] magnetic field. The finding
represents the first time such a coupling has been suggested in the pyrochlore
QSI literature, to the best of our knowledge.

Following Ref. 18, we recognise that these terms are decoupled from the
fluctuating background to leading order, manifesting instead as the three-
photon terms.