We assume throughout this work that the operators $S_i^{\alpha }$ transform as dipoles. As we discuss in Sec. 9, in the rare earth pyrochlores, it is frequently useful to use pseudospins to address the low-energy physics. These have mixed dipolar-quadrupolar character in non-Kramers ions [4].

Because crystalline symmetries are discrete, the largest irreducible representations of point groups are three-dimensional, so all multipoles higher than dipoles must become reducible.

We use lowercase $t_{2g}$ and $e_g$ to denote the decomposition of a single quadrupole into irreps, whereas we use the capital letter names to distinguish the collective decomposition of all four quadrupoles in a tetrahedron.

One can check that $|{\epsilon }_{\alpha \beta \gamma }|$ is an invariant symbol of $T_d$, meaning that $$|{\epsilon }_{\alpha \beta \gamma }|=|{\epsilon }_{\rho \delta \sigma }|R(g{)}_{\rho \alpha }R(g{)}_{\delta \beta }R(g{)}_{\sigma \gamma }$$for all $g\in T_d$, where $R(g)$ is the corresponding O(3) rotation.

The tensor $\mathcal{R}$ is one of three rank-2 “residue tensors” defined by pairwise contraction of the antisymmetric Levi-Civita symbol with a general rank-3 tensor [93]. In our case, two are identical and the third is zero, because ${\mathcal{M}}^{\alpha \beta \gamma }={\mathcal{M}}^{\alpha \gamma \beta }$ following from the symmetry of the quadrupole tensors $Q$.

In fact, its antisymmetric components are proportional to the trace vector ${\chi }_t^{\alpha }$ [Eq. (27)].

This description is redundant and covers the four-sphere multiple times, for example, $q(0,\pi /2,0,\mathcal{S},-\mathcal{T})=q(0,0,0,\mathcal{S},\mathcal{T})$. In practice, one can simply work with a unit five-component vector and reconstruct the symmetric trace-free matrix without any redundancy.

Note that, for uniaxial tensors with $\mathcal{T}=0$ in Eq. (54), $\mathrm{Tr}[q^3]={\mathcal{S}}^3/\sqrt{6}$, so the sign of $\mathrm{Tr}[q^3]$ is equivalent to the sign of $S$.

For $S>1$ it is sufficient to evaluate $⟨S|Q^{\alpha \beta }|S⟩$, but for $S=1$ there is also a nontrivial cross term coming from $⟨-1|Q^{\alpha \beta }|+1⟩$.

We thank Karlo Penc for pointing out that such a constraint exists for $S=1$.

The only difference between the Monte Carlo and projector calculations is the slight broadening of the intensity of the peak at the (331) point. This peak has an intensity of approximately 2.93 in the Luttinger-Tisza calculation and approximately 3.45 in the Monte Carlo calculation—it is not a Bragg peak; rather, it is a consequence of the long-wavelength fluctuations associated to the $\mathit{q}=\mathbf{0}$ band touching between the two flat bands and five dispersing bands, where the number of low-energy modes is seven rather than two. The broadening of the peak at finite temperature can be reproduced in a finite-temperature self-consistent Gaussian approximation calculation.