{'fraction_non_alphanumeric': 0.034786976861009074, 'fraction_numerical': 0.01343756867916235, 'mean_word_length': 4.664591854881736, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unconventional spin glass state in elemental neodymium in the absence of extrinsic disorder

Spin glasses are a fascinating but highly complex magnetic state of matter, traditionally investigated in materials with a random distribution of exchange interactions. Intricately linked to spin frustration and disorder, spin glasses exhibit no long-range order and exude aging phenomena, distinguishing them from quantum spin liquids. Here we report a new type of non-collinear spin glass state, a so-called spin-Q glass, surprisingly observable in crystalline metallic neodymium. Using spin-polarized scanning tunnelling microscopy combined with ab initio calculations and atomistic spin-dynamics simulations, we visualize the variations in atomic-scale non-collinear order and its response to magnetic field and temperature. We quantify the aging phenomena relating the glassy behaviour to the crystalline symmetry and the energy landscape. This not only resolves the long-standing debate of the magnetism of neodymium, but also suggests that glassiness may arise in other magnetic solids lacking extrinsic disorder, providing new platforms for brain-inspired hardware and topological quantum excitations.

Spin glasses are one of the more intriguing, but least understood magnetic states of matter [1][2][3][4][5] . In stark contrast to ferromagnets or antiferromagnets, which form a long-range ordered state when cooled, spin glasses form a state characterized by seemingly random and uncorrelated magnetic patterns. In this way, the amorphous-like magnetization pattern in spin glasses may be compared to the amorphous structure of glasses, like silicon dioxide, which exhibit local structural correlations but lack a long-range ordered state. Spin glass behaviour has recently spurred keen interest in many fields, ranging from iron-based superconductors 6 to theoretical machine learning 4,5,7 , and it has been suggested to be relevant in quantum topological excitations [8][9][10] . Moreover, the concept of glassiness has impact far beyond condensed matter, involving e.g. the theory of biological evolution 11,12 , proteins and polymers 13 , and cosmology 14 . Spin glasses are characterized by a glass transition temperature and are furthermore distinguished by the phenomenon known as aging. Namely, the magnetic state depends on its history, driven by a distribution of distinctive spin-relaxation processes with time scales spanning many orders of magnitude 4,5 . Aging also distinguishes spin glasses from so-called quantum spin liquids 15 or spin ices 16 , which remain disordered down to zero temperature due to quantum fluctuations and lack of memory. The paradigm of such complex magnets involves magnetic frustration, derived from geometry or competing interactions, but unlike for spin liquids, disorder is traditionally also considered necessary to drive nonergodic behaviour in spin glasses.
There is still no clear understanding when spin glass behaviour can arise in magnetic materials. The most commonly debated models 1,3-5 describe randomly distributed spins with long-range magnetic interactions of alternating ferromagnetic and antiferromagnetic coupling, like that seen in prototypical materials such as dilute magnetic alloys 1,17 . In the thermodynamic limit, this is characterized by a hierarchical energy landscape with infinitely many local energy minima separated by energy barriers of multiple heights so that there is a broad distribution of transition times between different minima [3][4][5]18 . This results in an absence of local and long-range order. It is commonly believed that disorder is a key ingredient, in addition to competing magnetic interactions, to realize spin glass behaviour. There is, however, a concept of self-induced glassiness introduced initially for the stripe glass behaviour in high-temperature superconductors 19,20 and later developed for magnets 21,22 . Within this picture, it has been recently theorized that competing interactions alone can lead to the glassy state, even in the absence of external disorder.
Here we show that single-crystalline elemental neodymium exhibits a new type of spin glass behaviour in the absence of extrinsic disorder. Utilizing spin-polarized scanning tunnelling microscopy of the Nd (0001) surface, the magnetic state is found to exhibit local non-collinear magnetic order while lacking a longrange ordered state. This local order is defined by a spectral distribution of degenerate magnetic wave vectors, or Q states, which varies spatially. We probe the response of this so-called spin-Q glass to applied magnetic fields and variable temperature, to quantify its aging behaviour and energy landscape.
Harnessing ab initio methods and simulations, we quantify the competing long-range magnetic interactions and the favourable Q states, illustrating that this unconventional glassy behaviour results from valley-like pockets of degenerate Q states as proposed for self-induced spin glasses 21,22 . Our findings not only suggest that glassy behaviour and aging can be found in systems with crystalline order, but also unravels an unresolved debate about the magnetic ground state of elemental neodymium that has challenged scientists for several decades.

Atomic-scale visualization of the spin-Q glass state
To characterize the atomic-scale magnetization, we utilize spin-polarized scanning tunnelling microscopy and spectroscopy (SP-STM/SP-STS) 23 and probe the surface of thick films of Nd(0001) grown on W(110) (Supplementary Fig. 1). Here, we image the surface, at low temperature, in constant-current mode at two characteristic voltages that are related to the energies of the exchange-split surface state 24,25 , and we consider the subtracted image ( Supplementary Fig. 2). We utilize an out-of-plane sensitive antiferromagnetic Cr probe, which relates the contrast variations of the image directly to the z-projection (i.e., the c-axis in the dhcp structure) of the magnetization, which we refer to as the magnetization image ( Fig. 1a,d, Supplementary Fig. 2d). The magnetization images of the surface reveal strong and clear short-range magnetic order with periodicities varying from = 0.9 -4.5 nm, oriented along or near the high-symmetry axes. These atomic-scale variations are directly related to local non-collinear magnetic order, with varying periodicities depending on spatial location, defined by a superposition of local magnetic wave vectors 2 / (Fig. 1d,e; Supplementary Fig. 3-4). Remarkably, while clear shortrange order can be seen, defined by a local multi-Q state, there is no observable long-range ordered state found for any of the probed experimental conditions.
In order to better visualize the variations of local multi-Q order, we consider reciprocal-space images via Fast Fourier Transform (FFT) of the real-space magnetization images, which we refer to as Q-space images. A signature of the lack of long-range order and competing short-range order can be directly visualized by the smeared and broadly distributed spectral weight in various regions in Q space (Fig. 1b).
Similarly, Q-space images can be produced from smaller spatial regions defined by a sharpened and characteristic spectral weight of Q vectors (Fig. 1d, Supplementary Fig. 4). Various line-cuts along the Γ -M direction of multiple Q-space images taken from different spatial regions of Fig. 1a are illustrated in Fig.   1c. The resultant plots reveal a spatial variation in spectral weight in at least three distinct regions, or Qpockets, with significant spectral weight along the high-symmetry axes. In the ensuing discussion, we focus particularly on the three Q-pockets with wave vectors at around QA = 1.1-2.0 nm -1 , QB = 2.7-3.5 nm -1 and QC = 4.6-5.3 nm -1 . Utilizing this information, we plot a schematic of out-of-plane projected magnetization with respect to the atomic lattice ( Fig. 1e) for the two regions shown in Fig. 1d with the defined Q vectors from these Q-pockets. The corresponding Q-space images are a direct visualization of the multi-Q nature acquired over the spatial area of the given images, and they provide a quantitative comparison to previous neutron diffraction studies [26][27][28][29][30][31][32] . The spectral weight around these pockets, as well as blurring of the intensity of the Q states (Fig. 1b) is a first indication of the glassy nature of the magnetic state, and very reminiscent of spin-based analogues of stripe or checkerboard order in strongly correlated compounds 33 .
In order to illuminate the concept of a spin-Q glass, we first qualitatively illustrate the energy landscape in Q-space images for a spin-Q glass in comparison to a ferromagnet (Fig. 2). A long-range ordered state can be related to a global minimum in Q space 23,34,35 , where a single-domain ferromagnetic state is equivalent to a Q = 0 global minimum (Fig. 2a). In contrast, a spin-Q glass is distinguished by the existence of flat valleys defined by a distribution of many local minima, i.e. Q-pockets, at finite Q values ( Fig. 2b). Broad Q-pockets lead to a lack of a preferentially long-range ordered state. Instead, there are only local regions defined by a spectral weight of mixed Q vectors within the given pockets, and different regions exhibit random distributions of this spectral weight (see colour image in Fig. 2b). This is distinct from the behaviour observed in multi-Q states for thin 3d transition metal films, which show well-ordered domains 36 . Within the concept of self-induced glassiness for spins, such pockets may result from a strong competition of magnetic interactions, leading to strongly degenerate states 21,22 .

Theoretical analysis of the magnetic landscape of Nd(0001)
In order to analyse the origin of the spin-Q glass state in Nd(0001) and its unique and unexpected magnetic patterns not found in any other magnetic metal 37 , we utilize ab initio calculations to quantify and understand the exchange interactions and the energy landscape in Q space of bulk Nd. This element resides in a dhcp structure, whose importance on the magnetic exchange interactions is visualized in Fig.   3a. The calculated exchange interactions for bulk Nd are shown for both the dhcp and a hypothetical hcp structure. We note the minute energy difference between these crystal structures, resulting from the large similarities in atomic arrangement. The distance dependence in the hcp structure illustrates a prototypical behaviour with ferromagnetic nearest neighbour interactions and an oscillating RKKY-like interaction at large distances, as seen for e.g. for Gd 38,39 . In contrast, the dhcp structure has a significantly different behaviour: at shorter range the interactions are significantly weaker and primarily antiferromagnetic, while at larger distances RKKY interaction sets in. This leads to strong competing magnetic interactions in Nd, which creates conditions for frustrated magnetism and/or spin-glass behaviour.
In order to clarify the impact of the calculated exchange interactions in the dhcp structure of Nd, we evaluate the magnetic energy landscape by means of single-Q spin spirals, i.e. magnetic structures that can be parametrized by a single wave vector. For this purpose, we calculate the energy of helical spin spirals. Using the calculated magnetic exchange interactions from Fig. 3a, we parametrize an effective Heisenberg spin Hamiltonian from which the energy E(Q) of the single-Q spirals is then calculated (Fig.   3b). The energy-landscape exploration was performed by fixing the z-component Qz and then sweeping over Qx and Qy in the first Brillouin zone (BZ) (Supplementary Fig. 5-6). In Fig. 3b, we present the single-Q energy landscape for all possible Qx and Qy combinations in the cell spanned by the dhcp reciprocal lattice vectors for the case when Qz = 2/c, which corresponds to the configuration with the lowest single-Q energy. This choice of Qz corresponds to a 90 rotation of the moments in adjacent atomic layers within the dhcp unit cell. The Q-dependent energy is colour-coded such that red (blue) regions correspond to spin spirals with low (high) energy. This visualization illustrates the complex energy landscape of Nd, as exemplified by the broad and flat dark-red ring-like structure. Instead of a distinct, six-fold degenerate set of strong energy minima that could be expected for a spin-spiral magnet on a hexagonal lattice, the red ring structure shows that the energy barriers between global and local minima in this region of the BZ are very small. In addition, high-energy local minima, or pockets (in red) are distributed with hexagonal symmetry, just inside the BZ, while another low-energy valley-like structure goes along the BZ boundary. We note that magnetic anisotropy is not considered here, and we expect that it will bias the ring structure toward the high-symmetry directions, creating pockets akin to those seen in the experiments. Thus, Fig. 3b illustrates that the energy landscape of Nd has several broad Qpockets, which supports the formation of the experimentally observed spin-Q glass structure. The presence of strongly competing interactions leading to a glass-like energy landscape in Q space is a key manifestation of the concept of self-induced glassines [19][20][21] .
In addition to the single-Q energy-landscape explorations, we have also employed Monte Carlo and Atomistic Spin Dynamics (ASD) simulations in order to find the ground-state magnetic structure. From these simulations, we obtain real-space spin structures (Supplementary Figs. 7-8) as well as the static correlation function S(Q), which can be compared directly with the experimental Q-space images. In Fig.   3c, we present the simulated peaks of S(Q) in comparison with (i) the observed range of spectral weight from the Q-pockets QA, QB and QC in the SP-STM experiments (Bz = 0), for both pristine zero-field cooled samples and samples after being exposed to magnetic field), and (ii) reported neutron diffraction data 29,31 .
In comparison, the SP-STM and neutron diffraction data agree well with the simulated results in regions of lower Q values. This agreement indicates that the surface-derived measurements are significantly coupled to the bulk magnetic properties. The simulations also provide distinct correlation functions at the Brillouin zone boundary, which so far have not been detected experimentally.

Magnetic field evolution of the spin-Q glass state
Next, we experimentally characterize the response of the magnetization to external magnetic fields. Outof-plane field dependence is illustrated in Fig. 4 for a few chosen fields up to Bz = 7 T at T = 1.3 K (Supplementary Figs. 9-10). The application of magnetic field should favour states with smaller Q, eventually leading to a preferential Q state near the zone center 34,40 , when the Zeeman energy exceeds the local exchange energy. However, magnetic imaging of a given area at variable magnetic field is contrary to this intuition. At increasing magnetic fields, no distinct Q state becomes favourable, with the spectral weight strongly broadening along high-symmetry directions toward higher Q (Fig. 4g). The magnetic order is sensitive to fields on the order of Bz = 0.5 T, illustrating the degeneracy driven by the Qpockets, while at the highest applied fields (Bz = 7 T;) there is no distinct Q state, demonstrating the strong local exchange energy. Moreover, the application of in-plane fields reveals similar behaviour: favourable Q states collapse onto an axis related to the direction of the applied field, but with a significant smearing of the spectral weight along one particular axis ( Supplementary Fig. 11). We note that the fielddependent behaviour seen here cannot be attributed to the picture of local domains, as there are neither clear domain boundaries that can be traced in magnetic field, nor a repeating or favourable Q-structure.
We note that magnetostriction effects should saturate beyond 1 T, at the measured temperature, as previously reported 41 .

Aging of the spin-Q glass state
While field-dependent imaging illustrates the degeneracy within the Q-pockets, the distinguishing property of a spin glass is the observation of aging. Aging can be described by the existence of multiple relaxation time scales, leading to a magnetization state that never fully relaxes 17 . To monitor the evolution of the ground state, we repeatedly apply the following procedure: (i) sweep up the magnetic field (max |Bz| = 7 T) for a finite time, (ii) reduce the magnetic field to zero, (iii) image the zero-field state of the same area.
After field-dependent cycling, the zero-field magnetic state at T = 1.3 K illustrates sequential redistributions in the spectral weight within and between the various Q-pockets ( The overall tendency is that the angular distribution of the spectral weight sharpens along the high-symmetry axes, while it broadens along the Q axis within all three Q-pockets, and there is a redistribution between QA, QB and QC (Fig. 5g). Moreover, the system never reverts to the initial zero-field cooled state or a singularly defined Q state, nor does it favour a periodicity defined by a region with a particular Q-state distribution. These points rule out that fielddependent cycling can be understood as the evolution of a favourable domain. To rule out hysteretic effects, we also performed field sweeps at positive and negative fields, but there is no clear correlation between such subsequent conditions, nor any correlation with local defects on the surface. We note that similar aging effects can be seen at different temperatures ( Fig. 5h), as we detail below, as well as with the application of an in-plane field ( Supplementary Fig. 13).
The evolution of intermittent spectral weight distributions frozen after subsequent field-sweep cycles is a signature of slow aging dynamics as well as a hierarchy of states separated by small energies. In contrast, the application of temperature illustrates the presence of fast dynamics concomitant with slow relaxation dynamics. At T = 4.2 K the spectral weight within the QA pockets "melts" away ( Fig. 5h), and the associated long-wavelength pattern is no longer visible in comparison to images of the same area at Supplementary Fig. 14). The loss of spectral intensity in the QA pockets is most likely due to the fact that our imaging method is too slow to resolve the increased fluctuations of the magnetization states within the QA pockets, which are thermally activated at this temperature. Nevertheless, we observe similar slow aging behaviour after application of magnetic field at T = 4.2 K for the other Q-pockets as seen at lower temperatures for exactly the same sample area ( Fig. 5h and Supplementary Fig. 14). While it is unclear how to describe the dynamical behaviour of the magnetization, the presence of multiple relaxation times clearly exemplifies the glassy behaviour in neodymium, analogous to aging effects seen in metallic alloys 3,17 . Counterintuitively, the application of higher temperature leads to more well-defined Q states (e.g. for T = 7 K); the same area imaged at T = 40 mK in comparison shows similar glassy distributions to the previously shown data taken at T = 1.3 K ( Supplementary Fig. 15). This rules out transitions into various long-range ordered multi-Q states 29,31 with decreasing temperature, but rather exemplifies the emergence of a spin-Q glass state. Finally, we observe that cleaner samples with reduced defect density

Discussion
We have demonstrated that the magnetic ground state of elemental neodymium is a spin-Q glass. We have observed aging, that is, glassiness, in an elemental magnetic solid that lacks chemical or extrinsic disorder. This represents the first direct experimental confirmation of the concept of self-induced glassiness 21         The UppASD code, http://www.physics.uu.se/uppasd.

Growth of Nd(0001) islands and spin contrast
Measurements were performed on epitaxial islands of Nd(0001) grown on W(110). The W(110) substrate was cleaned by repeated cycles of annealing at T = 1250 °C in a gradually reduced oxygen atmosphere (from p = 1 x 10 -7 mbar in the first cycle down to p = 2 x 10 -8 ) and flashing at T = 2400 °C, as described in Ref. 1   surface (Fig. S1c) reveals the expected surface state of Nd(0001), which is characterized by an exchange splitting into a majority peak visible below EF and a minority peak above EF, with an additional narrow peak at EF 2 . Previous studies have measured the temperature dependence of the exchange splitting ∆ , illustrating that the splitting persists far above TN = 19.9 K 3 . The slowly diminishing value of ∆ above TN is unusual compared to other lanthanide metals 6 , which is another signature of strong local correlations in Nd. 7 We note that the narrow peak at EF has been seen on other lanthanides as well 2,3,5 , i.e., there is no connections to the spin-Q glass behaviour only found in Nd so far. Comparing spectra at T = 1.3 K and T = 40 mK, this resonance narrows at lower temperature. Also, we do not see this feature strongly change in magnetic field. Therefore, we believe this feature is most likely related to a van Hove singularity.
In order to acquire spin contrast, we utilize a Cr bulk tip 8 , which we prepare by tip pulsing until achieving an optimal out-of-plane contrast (ref. 9 ). We note that we can sometimes observe probes with a canted magnetization, but we have never observed probes with only a pure in-plane and stable magnetization.
Constant-current SP-STM images near either peak of the exchange-split surface state using out-of-plane tips show the expected magnetic patterns in real space (Fig S2b, c) as well as in the FFT (Fig. S2f, g), as described in the main paper. Topographic effects due to subsurface defects are nearly identical at both imaging biases, while the magnetic contrast is nearly fully inverted. This is due to the fact that in the tunneling energy window defined by the applied bias, there is a bulk symmetry gap in the centre of the surface-projected Brillouin zone. Thus, tunneling almost exclusively happens between the tip and the highly spin-polarized surface state, i.e., out of the majority state at negative sample bias and into the minority state at positive sample bias, respectively. 10,11 For the magnetization images in the main manuscript, we acquired SP-STM images at VS = 200 mV and VS = -150 mV, respectively (Fig S2b, c), and subsequently subtracted them from each other (Fig S2d). Below, we illustrate that the overall changes of the FFT patterns in magnetic field for the individual SP-STM images taken at each surface state energy (Fig. S2f-h), are nearly identical to the subtracted magnetization images (cf. Fig. S10). We

Variation in local order and Q-state distribution
Magnetization images show a spatially dependent spectral weight of Q states, as illustrated in the FFT in Fig. 1b of the main manuscript. In Fig. S3, we utilize inverse FFT filtering at each identified Q-pocket and plot the resultant magnetization images, where the overall intensity is related to the contribution of the selected Q states. The images illustrate that the large variation in spectral weight of each of the Qpockets and that the underlying local order has various contributions from the named Q-pockets as schematically illustrated in Fig. 2b.
In Fig. 1c  where  is the Gilbert damping parameter that determines the rate of dissipation of energy from the spin system. Thermal effects are included in the ASD method by means of Langevin dynamics where a stochastic noise term is added to the effective magnetic field Bi. All simulations were done with the UppASD software 13 and the algorithms used are given in full detail in Ref. 14 .
For the correlation functions displayed in Fig. 3c in the main manuscript we followed a simulated cooling protocol where the temperature was decreased from T = 20 K (which is above the ordering temperature) with steps of ΔT = 5 K while performing 10 5 Monte Carlo sweeps at each temperature. The correlation functions were then sampled using 10 5 ASD steps with a Gilbert damping parameter of 0.01 and a time step of 10 -16 s. The Q-space correlation function S(Q) is obtained as the Fourier transform of the realspace trajectories of the system, which where sampled every 100 th time step.
The energetics of the single-Q spirals, displayed in Fig. 3b in the main text were evaluated by constructing helical spin spirals for given wave vectors Q, where the plane of rotation of the spins was defined to be perpendicular to the wave vector, i.e. the magnetization of the spin-spiral has both in-plane and out-of-plane components. The spatial magnetization profile for a single-Q spin spiral is constructed as:

̅ , sin ̅ • cos ̅ • ̂
where ̂ so that the in-plane component of the magnetization is perpendicular to the Q-vector.
The energy of each single-Q spiral was then obtained by evaluating the spin Hamiltonian defined above.
We note that since only isotropic Heisenberg exchange interactions are included in the Hamiltonian, the energy of a single-Q spin spiral does not depend on the orientation of the plane where the spins rotate, i.e., cycloidal and helical spin-states are degenerate for the same Q vector.
The schematic magnetization images shown in Fig. 1e in the main manuscript was obtained as linear combinations of spin-spirals with distinct Q vectors having different wave vectors (QA=1.5 nm -1 , QB=3.0 nm -1 , and QC=5.0 nm -1 ). This is analogous to an inverse Fourier transform of the selected Q vectors, with the difference that for the construction of the spirals used for the images in Fig. 1e, the magnetization of the resulting multi-Q spin-spirals were rescaled so that the norm of the magnetization was kept constant.
A proper linear combination of single-Q spirals would otherwise result in longitudinal fluctuations, i.e.
changes of the local magnetic moment magnitudes, but in order to stay consistent with our unit-vector based Hamiltonian we here enforced constant moment magnitudes.

Periodicity of the spin-spiral energy landscape
The single-Q spin spiral energy landscape E(Q), which is presented in Fig. 3b

Decomposition of the simulated magnetic structure for various Q vectors
The schematic magnetization images shown in Fig. 1e  stay consistent with our unit-vector based Hamiltonian, we here enforced constant moment magnitudes.
In Fig. S7, we illustrate the resultant site-dependent magnetization for various Q vectors taken from the calculation illustrated in Fig. 3 in the main manuscript.
As a supplement to the schematic picture in Fig. 2b in the main manuscript, we have simulated a quenching, i.e. a rapid cooling of the system going from T = 7 K to T = 0 K, using Metropolis MC simulations (Fig. S8)  magnetic fields applied close to one of the basal high-symmetry directions. [16][17][18] At intermediate fields, a variation in spectral weight is observed, similar to the out-of-plane experiments. This spectral weight broadens in higher applied fields, as illustrated by the line cuts shown in Fig. S11. We note here that due to the high defect density, the FFTs of this region are noisier in comparison to the out-of-plane data.

Aging in Nd(0001)
The distinguishing feature of spin glasses is the presence of aging dynamics. The most traditional definition of aging, as inspired by experimental evidence, is characterized by magnetization states which never fully relax resulting from the presence of multiple relaxation time scales, regardless of the waiting time. This relaxation behaviour can span many orders of magnitude in time, as exemplified by the magnetic alloys such as Cu-Mn 19,20 , which exhibits glassy dynamics below a particular freezing temperature above which the material is paramagnetic. illustrates a seemingly random change in the spectral weight of Q states within a given spatial region.
We also see aging behaviour for applied in-plane magnetic fields. Fig. S13 illustrates a region measured at T = 40 mK. We note that in the pristine state (Fig. S13a), there are predominately Q-states along one orientation. As this area was imaged near an edge of the island, we suggest that strain might have locked the Q states into this direction during the initial zero-field cool-down of the sample. However, after applying an in-plane field sweep, the Q states completely randomize their orientation (Fig. S13b), uninfluenced by the step edge. This illustrates the negligible impact of the island edges to the energy landscape. After subsequent in-plane field sweeps, we see similar behaviour as in out-of-plane fields, in which the spectral weight randomizes and shows no preferential distribution. Moreover, to compare the behaviour at T = 40 mK with the above data taken at T = 1.3 K, we also performed a field sweep in an (d-f, j-k) Corresponding Q-space images of the images above (scale bar = 3 nm -1 ). Warming up the surface from 1.3 K to 4.2 K results in depopulation of the QA pocket. Similar out-of-plane magnetic field dependence and aging behaviour is observed at T = 4.2 K for the other Q-pockets. The magnetization image in (i), measured after magnetic field-dependent measurements, is the aged state of the pristine state in (b).
out-of-plane field, on the same area (Fig. S13e), and observe the same aging behaviour that is illustrated in the main manuscript. Therefore, we conclude that the aging behaviour is independent of the applied field direction.
Temperature dependence of spin-Q glass state Fig. S14 illustrates some of the temperature-dependent data, whose line-cuts are shown in Fig. 5h of the main manuscript, including the magnetic field dependence and resultant aging at T = 4.2 K. Comparison of magnetization images taken at T = 1.3 K and T = 4.2 K, in the pristine state, shows that the spectral weight in the QA region is strongly reduced when going to higher temperature. The reduced intensity at low Q values is most likely due to the bandwidth of our measurement (~7 kHz), suggesting that the order around the QA pockets fluctuates faster than our time resolution in STM imaging, leading to a reduced intensity. Imaging the same region at T = 4.2 K in out-of-plane magnetic fields illustrates the same trends in the spectral weight in other regions of Q-space similar to the observations at T = 1.3 K. Namely, the spectral weight broadens along the other pockets without the appearance of any favourable Q state.
Likewise, aging behaviour is seen for these Q-pockets, similar to the behaviour seen at lower temperature, illustrating that there are multiple relaxation time` scales present ( Fig. S14b and Fig. S14i).
Previous neutron diffraction data illustrates the emergence of new Q states with decreasing temperature, starting at T < 19.9 K, down to T = 1.8 K. Similarly, we performed temperature-dependent measurements.
At T = 7 K (Fig. S15), we observe primarily Q states at high Q values with an absence of longer-range ordered states as observed at T = 1.3 K. Images of the same region at T = 40 mK illustrate a behaviour similar to that observed at T = 1.3 K, with the emergence of spectral weight in the QA pockets which are absent in images at T = 7 K. Moreover, we note that we do not observe any new Q states at T = 40 mK, in comparison to images taken at T = 1.3 K, suggesting that either the system is completely frozen into a finite number of Q-states, or that fluctuations persist down to lower temperature, inhibiting certain Qstates to freeze.

Effects of bulk defects
In Fig. 5i of the main manuscript, we illustrated the overall effect of defects on the spectral weight of Q states, where the magnetization images utilized are displayed in Fig. S13. The comparison between the dirty (Fig. S16a) and the clean (Fig. S16b)