Supplemental Material for “How frustrated magnets hide a phase transition and challenge quantum spin liquids”

The data presented in Figs. 1 and 2 were obtained by digitising data from published work, referenced in the text. Any error introduced in digitising is less than the size of the symbols in our plots. Converting susceptibility units from the conventional emu/mole to dimensionless units requires knowledge of the mass density. For the end-members of the series, we used either published density values or obtained values from the Materials Project (https://materialsproject.org/). For compounds that are solid-solutions of two end-members, we assumed a Vegard’s law (linear) interpolation of the lattice constant and thus density. Values of Tg were taken from the position of the peak in zero-field-cooled susceptibility.

pseudobrookite structure. Randomness arises from mixing of the F e 3+ and T i 4+ ions among their respective crystallographic sites. Thus, due to the different oxidation states of F e and T i, increasing F e 3+ content beyond 2/3 is not possible within the same structure without destroying the insulating phase.

II. GLASS TRANSITION IN A WEAKLY DISORDERED SYSTEM
In this section, we estimate the glass transition temperature in a system of spins with exchange interactions with randomly located generic defects and demonstrate its vanishing in the limit of the vanishing density n imp of the defects. Such defects may be represented, for example, by vacancies in the lattice of the system, spinless foreign atoms or perturbations in the exchange interaction between the spins. In order to analyse the glass transition, we employ the replica trick. The replicated Hamiltonian of the system is given by where we consider, for simplicity, classical Ising spins s r = ±1 labelled by their locations r and r in the lattice and the replica index α; the Hamiltonian H 0 describes the pure systems; H dis accounts for the effect of the defects; J rr is the coupling between spins at locations r and r and is a function of only the relative position r − r . Quenched disorder is represented by the randomness of the locations ρ i of the defects in Eq. (1c); the matrix K r 1 ,r 2 is non-random and describes the modifications of the couplings J r,r caused by one defect located at ρ = 0.
The defects thus create identical perturbations near different locations ρ i . For example, vacancies in a lattice with nearest-neighbour exchange coupling J r,r = J are described by the elements K r−ρ,0 = K 0,r−ρ = −J at all sites r adjacent to the location ρ of the vacancy, effectively cancelling the interactions between the spin at the location of the vacancy and the neighbouring spins, and to vanishing elements of K r 1 ,r 2 for all other positions. If the defect is represented by a foreign atom with the same spin as the atoms in the compound, the elements of the matrix K r 1 ,r 2 may represent small perturbations of the exchange couplings around a defect. If quenched disorder is represented by random-bond disorder on all sites coming from, e.g., random strain fields, it can be described by the same model with the defect density set to the density of the sites, n imp ∼ n.
We perform averaging of the replicated partition function of the system over the locations ρ i of all where N is the total number of sites in the system and N imp is the number of defects. Keeping the first two cumulants of the impurity-induced perturbation in the action of the disorder-averaged partition function Z gives the action where we have introduced the ratio n imp = N imp /N of the number of defects to the total number of sites in the lattice and the quantities Taking into account higher-order cumulants of the disorder-induced perturbation may change coefficients of order unity in the mean-field phase boundary of the glass phase but will not lead to qualitatively new effects.
Decoupling the variables s α r s β r by the Hubbard-Stratonovich field Q αβ rr gives whereJ r,r = J r,r +n imp k rr andF r 1 r 2 r 1 r 2 = F r 1 r 2 r 1 r 2 −n imp k r 1 r 1 k r 2 r 2 and summation over repeated indices is implied. We note that in the limit of low concentrations n imp of the defects, their effect on the quantitiesJ r,r ≈ J r,r andF r 1 r 2 r 1 r 2 ≈ F r 1 r 2 r 1 r 2 in Eq. (5) may be neglected.
Summing over the spin degrees of freedom leads to the free energy as a function of the fields Q given by where G rr = s α r s α r is the correlator of spins in the same replica in the pure system. To obtain the phase boundary of the glass phase, we look for the mean-field solution for the fields Q αβ This solution has the same symmetry Q αβ r 1 r 2 = Q αβ r 2 r 1 = Q βα r 1 r 2 as the averaged correlators s α r 1 s β r 2 .
Furthermore, although the glass phase breaks replica symmetry [10][11][12] , it is sufficient to look for the replica-symmetric solution to obtain the phase boundary.
Taking into account the symmetry Q αβ r 1 r 2 = Q αβ r 2 r 1 = Q βα r 1 r 2 , the free energy (6) can be rewritten as whereQ is our convention for the matrix Q αβ r 1 r 2 andF andĜ ⊗Ĝ are superoperators acting on matrices of the respective dimension;QFQ ≡ Q αβ r 1 r 2Fr1r2r 1 ,r 2 Q αβ The mean-field transition temperature is given by where λ is the minimum eigenvalue of the (positive) operatorF (Ĝ ⊗Ĝ). Because this operator is independent of the the defect concentration in the limit of low defect density n imp , the glass transition temperature vanishes in the limit of vanishing defect concentration.
We emphasise that the result (8)  N iGa 2 S 4 ) are layered with weak interlayer coupling and, therefore, may be considered as effectively two-dimensional (2D) systems at short distances. In this section, we argue that the 3D nature of those systems, which accounts for the interlayer coupling, is essential to observe the glass transition in them.
There is abundant numerical evidence (see, for example, Refs. 13-18) of the absence of a glass state at a finite temperature T > 0 in 2D spin systems with bond disorder. This changes in the presence of field disorder, for a system described by the Hamiltonian where b i is a random field. In this case, the system has a finite glass order parameter due to the identity in the limit of weak disorder strength κ = b 2 j dis . The presence of a finite glass order parameter in a system with random-field disorder is similar, in some sense, to the existence of a finite magnetisation in a magnet in a magnetic field. Because the order parameter q αβ is nonzero at all temperatures, there is no phase transition between the glass (q αβ = 0) and non-glass (q αβ = 0) phases.
Both in the case of bond-and random-field disorder, the glass transition is absent; the system either has the glass order parameter at all temperatures or never has it. Because the glass transition is observed in all geometrically frustrated compounds discussed in the main text of this paper, the three-dimensional nature of the discussed layered compounds ( and N iGa 2 S 4 ) is, therefore, important for the existence of the transition.

IV. CROSSOVER IN THE MAGNETIC PERMEABILITY IN THE COULOMB PHASE
In this section, we describe the mechanism of a crossover in the behaviour of a GF medium that drives the glass transition in a spin ice or a system of classical vector spins. Such a crossover is a property of the clean system and, albeit insufficient to cause a transition on its own, determines the conditions favourable for the onset of the glass phase.
We assume that the medium is described by the Coulomb phase [19][20][21][22][23] . The role of the variables in this description is played by a coarse-grained classical field B associated with the polarisations of the spins. In the absence of impurities and interactions other than the exchange interactions between the spins, the classical ground states of the system (states that minimise the total exchange energy) are hugely degenerate, corresponding to all configurations of the field B that satisfy the condition divB = 0.
Flipping a spin creates a pair of monopoles with opposite charges that can be further separated by local rotations of the spins. Creating an isolated monopole requires a large energy of the order of the AF coupling and will not be central to our arguments here. Defects with spins that mismatch the spins of the parent compound also correspond to monopole charges Q i at the locations of the defects.
Interactions other than the antiferromagnetic exchange interactions between the spins may lift the degeneracy of the abovementioned classical ground states. Here, we assume that such interactions are present and, in the limit of small B, lead to the temperature-independent contribution 1 8π B 2 dr (measured in appropriate units) to the free energy of the system. It may come, for example, from the magnetic dipole-dipole interactions between the spins or from interactions beyond nearest neighbours 24 . Another possible mechanism for generating this energy is quantum tunnelling [25][26][27] between classical ground states of the regions with divB = 0. While it does not require the existence of an energy scale other than the exchange energy, the amplitude of the tunnelling is in general exponentially suppressed by the number of spins that need to be rotated in order to arrive from one classical ground state to another [17][18][19]. The exact nature of the energy of the field B is not important for our consideration and cannot be uniquely inferred from the available data, which is why we leave the investigation of the origin of this energy for future studies.
The free energy of the systems also includes the "entropic" contribution 22  Combining all the discussed contributions gives the free energy of low-energy excitations in a GF material in the form where Q i and r i are the charges and the locations of the monopoles introduced by the defects; E Q is the energy of the defect, which may depend on the state of the defect spin; h(r) and h 0 (r) are weak quenched random fields that account for disorder other than monopole charges created by magnetic vacancies, e.g. random strain. In Eq. (S1), we have introduced the field φ, minimisation with respect to which enforces the Gauss theorem in the form divB = 4π i Q i δ(r − r i ). For states that minimise the free energy (12), the scalar field φ plays the role of the potential of the Coulomb field, ∇φ = − T +T T B + 4πh + h 0 B, where the random fields 4πh and h 0 B come from the fluctuations of the effective magnetisation and magnetic permeability. The free energy (S1) describes a system of randomly located magnetic charges (defects) in a disordered medium with the average magnetic permeability determined by the properties of the disorder-free system and given by Eq. (1) in the main text.
We emphasise that the permeability µ(T ) given by Eq. (1) characterises the interaction between excitations in the system and is in general unrelated to the magnetic susceptibility χ(T ), which characterises response to the external magnetic field. The energy E ≈ E 0 + 1 of a system of monopoles in the material consists of a temperature-independent contribution E 0 determined by the exchange coupling between spins and the interaction contribution ∝ [µ(T )] −1 .
The temperature dependence of the permeability µ(T ) can thus be independently verified via spin resonance spectroscopy of multi-charge, e.g, dipole, excitations in the system.
The glass state in the system described by the free energy (S1) requires sufficiently strong disorder represented by randomly located impurities, which create monopole charges Q i , or other sources of randomness represented by the quenched parameters h and h 0 in Eq. (1). We emphasise that the hidden energy scale T * , the glass transition temperature at zero impurity concentration, cannot be explained solely by the residual disorder distinct from impurities because the glass transition temperature grows with decreasing disorder. This trend, the exact microscopic mechanism of which we leave for future studies, suggests also a significant role of the disorder-free medium for the "hidden" energy scale T * discussed in this paper.