Herbertsmithite samples
ZnCu3(OH)6Cl2 single crystals were synthesized as described in ref. 54, using a recrystallization method. Powders of ZnCl2, CuO and H2O were mixed in a quartz tube with a ratio of 2.015 g:0.235 g:4.5 ml. The tube was sealed under vacuum and laid horizontally in a three-zone gradient furnace, with the temperature of hot and cold ends set at 180 °C and 160 °C, respectively. Millimetre-scale single crystals were obtained after 3 months. Extended Data Fig. 1a shows the photographs of the three ZnCu3(OH)6Cl2 single crystals studied. The lattice structure was confirmed by X-ray Laue diffraction, as exemplified by the clear Bragg peaks of Sample 1′ (Extended Data Fig. 1b). Measurements were performed along the c axis for Samples 1 and 2 and the a axis for Sample 3.
The stoichiometry of Zn:Cu ratio is found to be 0.97:3.03 for the samples reported here, by using inductively coupled plasma mass spectrometry. The refinement of a single-crystal X-ray diffraction measurement indicates that 32.5% of the Zn2+ sites and 10.8% of the Cu2+ sites are inter-substituted24. Extended Data Fig. 1c shows the d.c. susceptibility (Sample 1′) in SI units measured by a SQUID magnetic property measurement system (Quantum Design) and its Curie–Weiss fitting by \(\chi =\,{\chi }_{0}+\frac{{C}_{{\rm{Curie}}}}{T-{\theta }_{{\rm{CW}}}}\). Fitting in a temperature range of 150 K \(\le T\le\) 320 K yields \({\chi }_{0}=-5\times {10}^{-6}\), \({\theta }_{{\rm{CW}}}=-280\) K and \({C}_{{\rm{Curie}}}=0.165\) K. Fitting in a temperature range of 2 K \(\le T\le\) 6 K yields \({\chi }_{0}=(4.3\pm 0.2)\times {10}^{-4}\), \({\theta }_{{\rm{CW}}}=-1.07\pm 0.03\;{\rm{K}}\) and \({C}_{{\rm{Curie}}}=0.0134\pm 0.0002\;{\rm{K}}\) corresponding independently to S = 1/2 at 32.5% ± 0.5% of the Zn sites. The Curie–Weiss fitting of the low-temperature d.c. susceptibility is stable as long as the fitting range is within 2 K \(\le T\le 10\) K, although it becomes sharply fitting range dependent below 1 K at which the susceptibility starts diverting from the Curie–Weiss behaviour (Fig. 3c)14. These sample characterization results are comparable with past stoichiometry studies11,19,44, neutron diffraction studies11,19,55 and d.c. susceptibility studies14,34,56.
As the Cu occupation probability of witness spin sites in our single crystals, we take 33%. This value is based on the coincidence of estimates from the single-crystal X-ray diffraction measurement24 and d.c. susceptibility Curie–Weiss fitting (Extended Data Fig. 1c), both of which are performed on our single crystal (Sample 1′). We note that the precise nature of site disorder in herbertsmithite has not yet been fully determined11,22,23,57, with reported Cu substitution percentage at Zn sites ranging from 12% (ref. 17) to 36% (ref. 55), and of the Zn substitution percentage at Cu ranging from 0% (refs. 11,57) to 10% (refs. 22,55).
Measurements
Magnetic flux \(\varPhi (t)\) noise measurements were performed with a 19-turn single superconducting pickup coil connected to a d.c. SQUID SQ1200 (Star Cryoelectronics), which is designed to maximize the noise measurement sensitivity. Susceptibility measurements were performed with a 10-turn-each in-series counter-wound superconducting pickup coil connected to a SQUID SP550 (Quantum Design) (Fig. 1c), using a solenoid whose magnetic field was calibrated by using an indium cylinder in a superconducting state (\(\chi =\,-1\)). In both setups, three 0.2-mm-diameter silver wires were directly attached to the sample for thermalization, and measurements were taken at least 20 min after the target temperature was reached. Both setups were shielded by multiple nested niobium and mu-metal cylinders.
Flux noise data at each temperature was recorded for 1,000 s at 20 kSa s−1. For susceptibility measurements at each temperature, the sample magnetic response was recorded as the magnetic field \({\mu }_{0}H\) was swept over \(0\,{\rm{\mu }}{\rm{T}}\to \,-4\,{\rm{\mu }}{\rm{T}}\to 4\,{\rm{\mu }}{\rm{T}}\to -4\,{\rm{\mu }}{\rm{T}}\to 0\,{\rm{\mu }}{\rm{T}}\) in steps of \(\sim 0.05\,{\rm{\mu }}{\rm{T}}\) (zero-field cooling). The slope of these \(({\mu }_{0}{{M}})/({\mu }_{0}H)\) data, representing the magnetic susceptibility, is extracted by a linear fitting with standard error bars from the linear fit. The micro-Tesla d.c. susceptibility of Sample 1 (Fig. 3c) is obtained by subtracting an offset constant. This offset value is determined so that the measured micro-Tesla susceptibility smoothly connects to the measurement result in the magnetic property measurement system in the overlapping temperature range of 2 K ≤ T ≤ 3 K (Extended Data Fig. 1d). For the long-term spin evolution under a 2-μT field, a sample was thermalized at \({T}_{1}=400\,{\rm{mK}}\) for 1 h and the temperature was rapidly dropped to a lower temperature \({T}_{2}\) in less than 5 min. After 20 min of thermalization at \({T}_{2}\) (that is, starting from \(t=\) 1,200 s), the spin evolution was recorded for 80,000 s (~1 day) at 1 kSa s−1. In Fig. 3d (inset), the signal is averaged for every 100 s.
Flux noise data are processed in a similar method as that in ref. 32. The distribution of \(\varPhi (t)\) is Gaussian with the expected statistical fluctuations (Extended Data Fig. 2a). The PSD with frequency resolution \(\Delta \omega =\left(2{\rm{\pi }}\,{\rm{rad}}\right)\times \left(0.1\,{\rm{Hz}}\right)=0.6\) rad Hz is first calculated from 100 split segments32, with its error bars determined by the standard error of segment averaging. The empty-coil measurement result is subtracted as a background contribution. This PSD is plotted in Fig. 2c,d and Extended Data Fig. 2b. To increase the signal-to-noise ratio, the PSD is averaged over a 10\(\Delta \omega\) or \(100\Delta \omega\) window at high frequencies. The power-law index shown in Fig. 3a is obtained by fitting the \(\Delta \omega =0.6\) rad Hz PSD with \({S}_{\varPhi }(\omega ,T)\propto {\omega }^{-\alpha (T)}\) in the range of 0.6 rad Hz \(\le \omega \le\) 600 rad Hz (Extended Data Fig. 2b). The error bars are the standard error from fitting. The variance \({\sigma }_{\varPhi }^{2}\) shown in Fig. 3b is calculated by integrating the \(\Delta \omega =0.6\) rad Hz PSD in the range of 0.6 rad Hz \(\le \omega \le\) 600 rad Hz (error bars are propagated from the PSD). This is demonstrably equivalent to the variance directly calculated from the time sequence of \(\varPhi (t)\) (appropriately filtered to the corresponding frequency range), and the variance peak temperature remains close to \({T}^{* }\) for all different frequency-integration ranges (Extended Data Fig. 2c).
The EA spin glass order parameter \({q}_{{\rm{EA}}}\) shown in Fig. 3d is extracted from the micro-Tesla d.c. susceptibility by solving the formula58,59
$${\chi} (T)=\frac{C\left(1-{q}_{{\rm{EA}}}(T)\right)}{T-\theta \left(1-{q}_{{\rm{EA}}}(T)\right)}.$$
(8)
The constants \(C\) and \(\theta\) are obtained by fitting in the temperature range above \({T}^{* }\) for which \({q}_{{\rm{EA}}}(T)\) vanishes: 270 mK \(\le T\le\) 500 mK. The error bars of \({q}_{{\rm{EA}}}(T)\) propagate from \(\chi (T)\) and the standard error of the fitted parameters \(C\) and \(\theta .\)
The spin noise data shown in Figs. 2 and 3 were measured in Sample 1. The equivalent measurements were performed for Sample 2 and Sample 3. As shown in Extended Data Fig. 3, the transition at \(T^{*} \approx 260\,\mathrm{mK}\) in the witness spin noise power index \(\alpha\) from the PSD \({S}_{\varPhi }\left(\omega ,T\right)\propto {\omega }^{-\alpha (T)}\), the witness spin noise variance \({\sigma }_{\varPhi }^{2}\) and the micro-Tesla susceptibility \(\chi\); and the \(-\mathrm{ln}(t)\) relaxation of the sample flux \(\varPhi \left(t\right)\) below \({T}^{* }\) are reproduced in multiple samples. The Sample 3 response was smaller compared with the other two samples because the small crystal size made it difficult to fill the full length of the pickup coil. Accordingly, the PSD fitting range is limited to 0.6 rad Hz \(\le \omega \le\) 60 rad Hz and the susceptibility measurement result in Extended Data Fig. 3c is scaled for comparison with Sample 1.
The witness spin dynamics and the associated transition that we observe in all herbertsmithite samples had not been previously measured in either a.c. or d.c. susceptibility studies14,44. One possible reason could be a difference in measurement conditions. Pioneering d.c. susceptibility measurements14 were performed at magnetic fields near \(B=\) 0.05 T. Although that field is small compared with the energy scale of the observed transition temperature \({T}^{* }\approx 260\,\mathrm{mK}\), it is empirically known that fairly a small field can considerably suppress a sharp peak signature of spin glass transition, making it difficult to detect. For example, a d.c. field of \(B=\) 0.04 T is capable of suppressing the sharp peak signature of a 21.5-K spin glass transition in Fe0.5Mn0.5TiO3 (ref. 60), 0.05 T for a 17-K transition in CdCr1.7In0.3S4 (ref. 61), 0.06 T for a 15.5-K transition in Gd0.37Al0.63 (ref. 62) and 0.04 T for a 0.2-K transition in Gd3Ga5O12 (ref. 63). Another possible reason is the difference between powder samples and single crystals. The herbertsmithite susceptibility studies14,44 were performed on powder samples before the establishment of single-crystal growth54. Differences between single crystal and such powder samples might possibly be caused by different chemical compositions, enhanced surface effects and randomized direction of applied field, which may have prevented the observation of a spin glass transition therein. All of the d.c. susceptibility measurements reported in the present work were carried out at \(B\le 5\,{\rm{\mu }}{\rm{T}}\) on single crystals, and all of them yield a sharp transition at a virtually identical \({T}^{* }\approx 260\,\mathrm{mK}\), supporting the plausible conclusion that this phenomenon is intrinsic to herbertsmithite single crystals in ambient magnetic fields of \({|B|}\le 5\,{\rm{\mu }}{\rm{T}}.\)
Spinon-mediated interactions via Z
2 QSL
The Hamiltonian for mutual witness spin interactions in equations (1) and (2)40,64 is derived from the coupling between a witness spin \({{\boldsymbol{s}}}_{i}\) and a kagome spin \({{\boldsymbol{s}}}_{l}^{{\rm{Kagome}}}\) as
$${H}_{{\rm{coupling}}}=\gamma {{\boldsymbol{s}}}_{i}\cdot {{\boldsymbol{s}}}_{l}^{{\rm{Kagome}}}.$$
(9)
The intra-kagome spin susceptibility \({\zeta }_{{lm}}\) is calculated using linear response equations (3) and (4) from the spinon band structure of a Z2\(\left[0,\,{\rm{\pi }}\right]\beta\) QSL. A Z2\(\left[0,\,{\rm{\pi }}\right]\beta\) QSL is the only gapped QSL that is compatible with lattice symmetries at the mean field level, and is in the neighbourhood of the U(1)\(\left[0,\,{\rm{\pi }}\right]\) state whose energy is the lowest among different U(1) QSLs28.
In units of the nearest-neighbour spinon hopping energy \({t}_{1}=0.4{J}_{{\rm{K}}}\approx 76\,{\rm{K}}\) (ref. 28), the Z2\(\left[0,\,{\rm{\pi }}\right]\beta\) QSL Hamiltonian contains real parameters for the second-neighbour hopping \({t}_{2}\), gap \({\Delta }_{2}\), and two Lagrange multipliers \({\lambda }_{1}\) and \({\lambda }_{3}\) which enforce the physical Hilbert space constraint of half-filling. We solve for all four using the standard self-consistent mean field approach. Reference 36 gives the following mean field spinon Hamiltonian for the kagome planes:
$${{\boldsymbol{s}}}_{l}^{\mathrm{Kagome}}=\frac{1}{2}\mathop{\sum }\limits_{\alpha ,\beta =\left\{\uparrow ,\downarrow \right\}}{f}_{i\alpha }^{\,\dagger }{\sigma }_{\alpha \beta }\,{f}_{i\beta },$$
(10)
$$\begin{array}{l}{\hat{H}}_{\mathrm{QSL}}=\left\{{\mathop{\sum}\limits_{i}^{N}}{\lambda }_{3}\left({f}_{i\uparrow }^{\,\dagger }{f}_{i\uparrow }+{f}_{i\downarrow }^{\,\dagger }{f}_{i\downarrow }\right)+{\lambda }_{1}\left({f}_{i\uparrow }^{\,\dagger }{f}_{i\downarrow }^{\,\dagger }+{f}_{i\downarrow }{f}_{i\uparrow }\right)\right\}\\\qquad\quad+\left\{\mathop{\sum }\limits_{ij}\left({t}_{1}{\nu }_{ij}^{\left(1\right)}+{t}_{2}{\nu }_{ij}^{\left(2\right)}\right)\left({f}_{i\uparrow }^{\,\dagger }{f}_{j\uparrow }+{f}_{i\downarrow }^{\,\dagger }{f}_{j\downarrow }-{f}_{i\uparrow }{f}_{j\uparrow }^{\,\dagger }-{f}_{i\downarrow }{f}_{j\downarrow }^{\,\dagger }\right)\right.\\\left.\qquad\quad\qquad\quad+{\Delta }_{2}{\nu }_{ij}^{\left(2\right)}\left({f}_{i\uparrow }^{\,\dagger }{f}_{j\downarrow }^{\,\dagger }-{f}_{i\downarrow }^{\,\dagger }{f}_{j\uparrow }^{\,\dagger }-{f}_{i\uparrow }{f}_{j\downarrow }+{f}_{i\downarrow }\,{f}_{j\uparrow }\right)\vphantom{\mathop{\sum }\limits_{ij}}\right\},\end{array}$$
(11)
where \({\nu }_{{ij}}^{\left(1\right)}\) is non-zero only for first-nearest neighbours (and is \(1\) or \(-1\) as defined in ref. 36), and \({\nu }_{{ij}}^{\left(2\right)}\) is non-zero only for second-nearest neighbours. There are \(N\) sites in the system. It is convenient to rewrite the diagonal terms using
$$\left\{{f}_{i},{f}_{j}^{\,\dagger }\right\}={\delta }_{{ij}},$$
(12)
$$\left\{{f}_{i},{f}_{j}\right\}=0,$$
(13)
which gives
$$\begin{array}{l}{\hat{H}}_{\mathrm{QSL}}=N{\lambda }_{3}+\left\{\mathop{\sum }\limits_{i}\frac{{\lambda }_{3}}{2}\left({f}_{i\uparrow }^{\,\dagger }{f}_{i\uparrow }-{f}_{i\uparrow }{f}_{i\uparrow }^{\,\dagger }+{f}_{i\downarrow }^{\,\dagger }{f}_{i\downarrow }-{f}_{i\downarrow }{f}_{i\downarrow }^{\,\dagger }\right)\right.\\\left.\qquad\qquad\qquad\qquad\quad\;+\frac{{\lambda }_{1}}{2}\left({f}_{i\uparrow }^{\,\dagger }{f}_{i\downarrow }^{\,\dagger }+{f}_{i\downarrow }{f}_{i\uparrow }-{f}_{i\downarrow }^{\,\dagger }{f}_{i\uparrow }^{\,\dagger }-{f}_{i\uparrow }{f}_{i\downarrow }\right)\vphantom{\mathop{\sum }\limits_{i}}\right\}\\\qquad\qquad\quad\quad\;+\left\{\mathop{\sum }\limits_{ij}\left({t}_{1}{\nu }_{ij}^{\left(1\right)}+{t}_{2}{\nu }_{ij}^{\left(2\right)}\right)\left({f}_{i\uparrow }^{\,\dagger }{f}_{j\uparrow }+{f}_{i\downarrow }^{\,\dagger }{f}_{j\downarrow }-{f}_{i\uparrow }{f}_{j\uparrow }^{\,\dagger }-{f}_{i\downarrow }{f}_{j\downarrow }^{\,\dagger }\right)\right.\\\left.\qquad\qquad\qquad\qquad\quad\;+{\Delta }_{2}{\nu }_{ij}^{\left(2\right)}\left({f}_{i\uparrow }^{\,\dagger }{f}_{j\downarrow }^{\,\dagger }-{f}_{i\downarrow }^{\,\dagger }{f}_{j\uparrow }^{\,\dagger }-{f}_{i\uparrow }{f}_{j\downarrow }+{f}_{i\downarrow }{f}_{j\uparrow }\right)\vphantom{\mathop{\sum }\limits_{ij}}\right\}.\end{array}$$
(14)
The initial \(N{\lambda }_{3}\) acts as an overall chemical potential and can be dropped. The following basis is then block diagonal:
$${\hat{H}}_{\mathrm{QSL}}=\mathop{\sum }\limits_{ij}\left(\begin{array}{cc}\left(\begin{array}{cc}{f}_{i\uparrow }^{\,\dagger } & {f}_{i\downarrow }\end{array}\right) & \left(\begin{array}{cc}{f}_{i\uparrow } & {f}_{i\downarrow }^{\,\dagger }\end{array}\right)\end{array}\right)\left(\begin{array}{cc}({h}_{{ij}}) & (0)\\ (0) & (-{h}_{{ij}})\end{array}\right)\left(\begin{array}{l}\left(\begin{array}{l}{f}_{j\uparrow }\\ {f}_{j\downarrow }^{\,\dagger }\end{array}\right)\\ \left(\begin{array}{l}{f}_{j\uparrow }^{\,\dagger }\\ {f}_{j\downarrow }\end{array}\right)\end{array}\right),$$
(15)
where
$${h}_{{ij}}=\left(\begin{array}{cc}\frac{{\lambda }_{3}}{2}{\delta }_{{ij}}+{t}_{a}{\nu }_{{ij}}^{a} & \frac{{\lambda }_{1}}{2}{\delta }_{{ij}}+{\Delta }_{2}{\nu }_{{ij}}^{\left(2\right)}\\ \frac{{\lambda }_{1}}{2}{\delta }_{{ij}}+{\Delta }_{2}{\nu }_{{ij}}^{\left(2\right)} & -\frac{{\lambda }_{3}}{2}{\delta }_{{ij}}-{t}_{a}{\nu }_{{ij}}^{a}\end{array}\right)$$
(16)
(a sum over \(a=1,\,2\) is implicit). Hence, all the information is contained in the upper matrix:
$${\hat{H}}_{\mathrm{QSL}}^{U}=\mathop{\sum }\limits_{ij}{\left(\begin{array}{cc}{f}_{i\uparrow }^{\,\dagger } & {f}_{i\downarrow }\end{array}\right)}_{\alpha }{h}_{{ij}}^{\alpha \beta }{\left(\begin{array}{l}{f}_{j\uparrow }\\ {f}_{j\downarrow }^{\,\dagger }\end{array}\right)}_{\beta }.$$
(17)
The self-consistency conditions are
$${\Delta }_{{ij}}=-2\left\langle {f}_{i\uparrow }{f}_{j\downarrow }\right\rangle =2\left\langle {f}_{i\downarrow }\,{f}_{j\uparrow }\right\rangle ,$$
(18)
$${t}_{{ij}}=2\left\langle {f}_{i\uparrow }^{\,\dagger }{f}_{j\uparrow }\right\rangle =2\left\langle {f}_{i\downarrow }^{\,\dagger }{f}_{j\downarrow }\right\rangle ,$$
(19)
$$0=\left\langle {f}_{i\uparrow }{f}_{j\uparrow }\right\rangle =\left\langle {f}_{i\downarrow }{f}_{j\downarrow }\right\rangle =\left\langle {f}_{i\uparrow }^{\,\dagger }{f}_{j\downarrow }\right\rangle =\left\langle {f}_{i\downarrow }^{\,\dagger }{f}_{j\uparrow }\right\rangle .$$
(20)
The global half-filling constraint on the physical Hilbert space is enforced by the Lagrange multipliers \({\lambda }_{1}\) and \({\lambda }_{3}\):
$${\lambda }_{1}:0=\mathop{\sum }\limits_{i}\left\langle {f}_{i\uparrow }{f}_{i\downarrow }\right\rangle -\left\langle {f}_{i\downarrow }{f}_{i\uparrow }\right\rangle ,$$
(21)
$${\lambda }_{3}:1=\mathop{\sum }\limits_{i}\left\langle {f}_{i\uparrow }^{\,\dagger }{f}_{i\uparrow }\right\rangle +\left\langle {f}_{i\downarrow }^{\,\dagger }{f}_{i\downarrow }\right\rangle .$$
(22)
This must now be diagonalized:
$${\hat{H}}_{\mathrm{QSL}}^{U}=\mathop{\sum }\limits_{ij}\left(\begin{array}{cc}{\gamma }_{i1}^{\dagger } & {\gamma }_{i2}^{\dagger }\end{array}\right){D}_{{ij}}\left(\begin{array}{l}{\gamma }_{j1}\\ {\gamma }_{j2}\end{array}\right)$$
(23)
with diagonal \(D\), and
$$h={UD}{U}^{\dagger },$$
(24)
$$\left(\begin{array}{c}{f}_{i\uparrow }\\ {f}_{i\downarrow }^{\dagger }\end{array}\right)={U}_{{ij}}\left(\begin{array}{c}{\gamma }_{j1}\\ {\gamma }_{j2}\end{array}\right)=\left(\begin{array}{c}{U}_{{ij}}^{11}{\gamma }_{j1}+{U}_{{ij}}^{12}{\gamma }_{j2}\\ {U}_{{ij}}^{21}{\gamma }_{j1}+{U}_{{ij}}^{22}{\gamma }_{j2}\end{array}\right),$$
(25)
and the Hermitian conjugate gives the other required terms:
$$\left(\begin{array}{cc}{f}_{i\uparrow }^{\,\dagger } & {f}_{i\downarrow }\end{array}\right)=\left(\begin{array}{cc}{U}_{{ij}}^{11* }{\gamma }_{j1}^{\dagger }+{U}_{{ij}}^{12* }{\gamma }_{j2}^{\dagger } & {U}_{{ij}}^{21* }{\gamma }_{j1}^{\dagger }+{U}_{{ij}}^{22* }{\gamma }_{j2}^{\dagger }\end{array}\right).$$
(26)
In this basis,
$$\left\{{\gamma }_{i\alpha },{\gamma }_{j\beta }\right\}=\left\{{\gamma }_{i\alpha }^{\dagger },{\gamma }_{j\beta }^{\dagger }\right\}=0,$$
(27)
$$\left\{{\gamma }_{i\alpha },{\gamma }_{j\beta }^{\dagger }\right\}={\delta }_{{ij}}{\delta }_{\alpha \beta },$$
(28)
$$\left\langle {\gamma }_{i\alpha }^{\dagger }{\gamma }_{j\beta }\right\rangle ={\delta }_{{ij}}{\delta }_{\alpha \beta }{n}_{{\rm{D}}}\left({D}_{{ii}}^{\alpha \alpha }\right),$$
(29)
where \({n}_{{\rm{D}}}\) is the Fermi–Dirac distribution. However, note that the eigenvalues \({D}_{{ii}}\) are ordered low to high, and the spectrum is symmetric about zero. Hence, at \(T\approx 0\) (since \({T}^{* }/{J}_{{\rm{K}}}=0.26\,{\rm{K}}/190\,{\rm{K}}\ll 1\)), \({n}_{{\rm{D}}}\left({D}_{{mm}}^{22}\right)=0\) and \({n}_{{\rm{D}}}\left({D}_{{mm}}^{11}\right)=1\). Feeding these expressions into the self-consistency conditions gives
$${\Delta }_{2}{\nu }_{{ij}}^{\left(2\right)}=2\mathop{\sum }\limits_{{{m}}}{U}_{{im}}^{11}{U}_{{jm}}^{21* },$$
(30)
$${t}_{a}{\nu }_{{ij}}^{\left(a\right)}=2\mathop{\sum }\limits_{{{m}}}{U}_{{im}}^{11* }{U}_{{jm}}^{11},$$
(31)
$${\lambda }_{1}:0=\mathop{\sum }\limits_{im}{U}_{{im}}^{12}{U}_{{im}}^{22* }-{U}_{{im}}^{21* }{U}_{{im}}^{11},$$
(32)
$${\lambda }_{3}:1=\mathop{\sum }\limits_{i}{\left|{U}_{{ii}}^{11}\right|}^{2}+{\left|{U}_{{ii}}^{22}\right|}^{2}.$$
(33)
We set \({t}_{1}=1\), defining the energy scale. Working in \(q\) space at \({q}=\,0\) (since the gap should be constant), we found a self-consistent solution with
$${\Delta }_{2}=0.4583,$$
(34)
$${t}_{2}=-0.2849,$$
(35)
$${\lambda }_{1}=0.4327,$$
(36)
$${\lambda }_{3}=1.500.$$
(37)
We used a tolerance of \({10}^{-3}\) in finding the constraints with the Lagrange multipliers:
$${\lambda }_{1}\Rightarrow \mathop{\sum }\limits_{im}{U}_{{im}}^{12}{U}_{{im}}^{22* }-{U}_{{im}}^{21* }{U}_{{im}}^{11}=-8.4\times {10}^{-4}\,\left( \sim 0\right),$$
(38)
$${\lambda }_{3}\Rightarrow \mathop{\sum }\limits_{i}{|{U}_{{ii}}^{11}|}^{2}+{|{U}_{{ii}}^{22}|}^{2}=1.001\left( \sim 1\right).$$
(39)
The overall energy gap (identified from the density of states),
$$2\Delta =0.44{t}_{1}=33\,{\rm{K}},$$
(40)
is essentially equal to the gap (\(0.43{t}_{1}\)) identified previously using exact diagonalization39.
In the witness–witness spin interactions \({J}_{{ij}}\) in equations (1)–(4), the only free parameter is \(\gamma\). We constrain \(|\gamma |=60\,{\rm{K}}\approx {J}_{{\rm{K}}}/3\) by requiring a match to the widely reported experimental value of the Curie–Weiss temperature \({\theta }_{{\rm{CW}}}\left(1\,{\rm{K}} < T\right)=-1.1\,{\rm{K}}\).
$${r}_{0}=\frac{\hslash {v}_{{\rm{F}}}}{2\Delta },$$
(41)
where the spinon band structure enters via the gap \(2\Delta\), and the spinon Fermi velocity \({v}_{{\rm{F}}}\) of the parent U(1)\(\left[0,\,{\rm{\pi }}\right]\) gapless QSL from which the Z2\(\left[0,\,{\rm{\pi }}\right]\beta\) forms28 with
$${v}_{{\rm{F}}}=\frac{{\sqrt{2}t}_{1}d}{\hslash },$$
(42)
where \(d\) is the nearest-neighbour kagome Cu spacing. Equations (41) and (42) lead to equation (6).
Witness spin Monte Carlo simulations
In herbertsmithite, the witness spin sites (that is, Zn2+ sites) form a triangular lattice on the \(ab\) plane, staggered along the \(c\) axis with a period of 3. This witness spin lattice effectively connects as a simple cubic lattice17. We simulate a witness spin lattice with the size of \(X\times X\times Z=45\times 45\times 4.\) The smallest cell containing one witness spin, which is \(1\times 1\times 1\) under this notation, is a rhombic prism with the side length \(a/\sqrt{3}=3.95\,\mathring{\rm A}\) and height \(c=14.09\,\mathring{\rm A}\). The direction of its rhombic base is rotated by 90° around the c axis, compared with the rhombic base of the conventional unit cell of herbertsmithite. To satisfy periodic boundary conditions, \(X\) has to be a multiple of three and \(Z\) has to be an even number.
We created a Monte Carlo simulation of the witness spins using the Metropolis–Hasting algorithm. We modelled the witness spins as classical Ising spins \({s}_{i}=\pm 1/2\). Although witness spins in herbertsmithite are not Ising like, they are not Heisenberg like, either. Electron spin resonance demonstrates a strong Dzyaloshinskii–Moriya interaction (\(D/{J}=\,0.08\)), leading to spin anisotropy45. Moreover, there is evidence for additional easy-axis anisotropy beyond the Dzyaloshinskii–Moriya interaction65. A key consequence of this magnetic anisotropy is that using a pure Heisenberg representation to model witness spin dynamics would be incorrect and that using an Ising representation is a reasonable approximation choice. The use of Ising spins has a further pragmatic justification in the context of spin glass. Adding a tiny amount of anisotropy to a Heisenberg spin glass can lead to a spin glass in the Ising universality class, making the Ising-spin model effective in reproducing experimental observations66. We initialized the system with the size 45 × 45 × 4 and periodic boundary conditions with 33% of potential witness spin sites occupied (\(N=2,673\) spins). The configuration of occupied sites is randomly assigned using a seeded random number generator. We average all of our results over 128 different configurations, keeping the same sets of seed across all runs. For each seeded configuration, we also average our results over three simulation runs.
In all cases, the initial state of the system has each spin in random uncorrelated states, corresponding to infinite temperature. We then run our simulation starting at \({T}=\,\) 400 mK—above the freezing transition temperature so as to avoid quenching the glass—and ending at 50 mK at intervals of 10 mK. In addition, we also conduct a separate run going from 10 K to 1.6 K at intervals of 0.2 K to ensure that the Curie–Weiss temperature is \({\theta }_{{\rm{CW}}}=-1.1\,\) K, and from 1.5 K to 0.5 K at intervals of 0.1 K for completeness.
At each temperature point, we first equilibrate the system by updating the system over 1,000 sweeps, with each sweep consisting of \(N=2,673\) update steps. We then sampled the spin-per-site noise \(s\left(t\right)=\frac{1}{N}{\sum }_{i}{s}_{i}(t)\), EA spin glass order parameter \({q}_{{\rm{EA}}}=\frac{1}{N}{\sum }_{i}{\overline{2{s}_{i}\left(t\right)}}^{\,2}\) and AF order parameter \({\phi }_{\mathrm{AF}}=\overline{\frac{1}{N}{\left({\sum }_{i}{\left(-1\right)}^{k}(2{s}_{i}(t))\right)}^{2}}\) (\(k=\mathrm{0,1}\) for each sublattice of the bipartite witness spin sites) over 100,000 sweeps. The bar represents an average over the Monte Carlo sweep time. From \(s\left(t\right)\), the magnetization noise \({{M}}(t)\) and d.c. magnetic susceptibility \(\chi\) are estimated using
$${{M}}\left(t\right)={\rho }_{V}{\mu }_{0}g{\mu }_{{\rm{B}}}s\left(t\right)\sqrt{\frac{N}{{N}_{{\rm{EXP}}}}},$$
(43)
$$\chi ={\rho }_{V}{\mu }_{0}{\left(g{\mu }_{{\rm{B}}}\right)}^{2}N\frac{\overline{{s}^{2}\left(t\right)}-{\left(\overline{s\left(t\right)}\right)}^{2}}{{k}_{{\rm{B}}}T},$$
(44)
where \({\rho }_{V},{\mu }_{0},g\,=\,2,\,{\mu }_{{\rm{B}}}\,\mathrm{and}\,{k}_{{\rm{B}}}\) are the number density per volume of witness spins in herbertsmithite (33% per Zn sites), vacuum permeability, electron \(g\)-factor, Bohr magneton and Boltzmann constant, respectively. The factor \(\sqrt{N/{N}_{{\rm{EXP}}}}\), where \({N}_{{\rm{EXP}}}\) is the number of witness spins in the volume of herbertsmithite Sample 1 (~3 mm3), is required to approximately estimate the order of the magnetization noise magnitude that generally scales as \({{M}}(t)\propto 1/\sqrt{{N}_{{\rm{EXP}}}}\) (ref. 32). The error bars of \(\chi\), \({\phi }_{{\rm{AF}}}\) and \({q}_{{\rm{EA}}}\) are the standard error of averaging.
The predicted witness spin magnetization noise \({{M}}(t)\) is then processed in the same method as the experimental spin noise in ZnCu3(OH)6Cl2 (see the ‘Measurements’ section). The distribution is Gaussian with small statistical fluctuations (Extended Data Fig. 4a). The PSD in Fig. 4c,d and Extended Data Fig. 4b has a frequency resolution of \(\Delta \omega =\left(2{\rm{\pi }}\,{\rm{rad}}\right)\times \left({10}^{-5}/{\rm{MCS}}\right)=6\times {10}^{-5}\) rad/MCS and is averaged over a 10\(\Delta \omega\) or \(100\Delta \omega\) window at high frequencies; the error bars are the standard error of averaging. The power index in Fig. 5a is obtained by fitting in \(6\times {10}^{-5}\) rad/MCS \(\le \omega \le\) \(1\times {10}^{-3}\) rad/MCS (Extended Data Fig. 4b); the error bars are the standard error from fitting. In Fig. 5a, open symbol points are used at temperatures above \({T}^{* }\) at which the power-law fitting is challenging (R < 0.98). The variance in Fig. 5b is calculated by integrating the PSD from \(6\times {10}^{-5}\) rad/MCS \(\le \omega \le\) \(6\times {10}^{-2}\) rad/MCS. The variance peak temperature remains at \({T}^{* }\) for different integration ranges (Extended Data Fig. 4c). When 1 MCS = 100 μs, the simulated PSD and the measured experimental PSD roughly correspond in the same frequency window (Extended Data Fig. 5). They remain consistent with each other for any values in the range of 1 MCS < 100 μs, as long as the PSD continues to be scale invariant down to lower frequency in both simulation and experiment. The challenges in fitting PSD at temperatures above \({T}^{* }\), which were not seen in the experiment, may be resolved by simulating the PSD in a lower Monte Carlo frequency range or by simulating the fully quantum mechanical theory.
When we perform the equivalent simulations for different witness spin concentrations from 15% to 60%, the transition temperature \({T}^{* }\) changes from 200 mK to 100 mK, and the nearest-neighbour witness spin interaction energy scale (equation (7)) changes from 2.5 K to 0.5 K. Even if these different witness spin concentrations are used in the model, these quantitative changes do not alter the conclusion of this work.
We finally note that whether the predicted transition is of true spin glass type has not been examined in detail. Direct theoretical investigations into this point, such as finite-size scaling of the spin glass susceptibility, are left for future work.
Neutron scattering structure factor
Extended Data Fig. 6 shows the witness spin structure factor \({\Sigma} ({\boldsymbol{q}})\) that is calculated from a spin configuration snapshot in our Monte Carlo simulation at 2 K over 1,000 sweeps and averaging over 128 configurations:
$${\Sigma} \left({\boldsymbol{q}}\right)=\left|F\left({\boldsymbol{q}}\right)\right|^{2}\,\left\langle \left|\sum \limits_{i}{s}_{i}{{\rm{e}}}^{{{i}}{\boldsymbol{q}}\cdot {{\boldsymbol{r}}}_{i}}\right|^{2}\right\rangle ,$$
(45)
where \(F({\boldsymbol{q}})\) is the magnetic form factor of Cu2+. It reasonably matches the low-energy neutron scattering structure factor in ZnCu3(OH)6Cl2, which shows diffuse scattering without a sharp peak17,26, although with a fairly low signal-to-noise ratio. Comparable features are the shape of the mid-intensity contribution (green) that extends throughout the in-plane and out-of-plane directions, the high-intensity contribution (red) at the out-of-plane peak at (00\(\frac{3}{2}\)), and the high-intensity in-plane circular shape contribution with correct \(|{\boldsymbol{q}}|\) and approximately equally distributed intensity. Further comparison of the precise shape of the experimental neutron scattering intensity pattern in the (HK0) plane requires improvement both in the precision of the neutron scattering experiment and in the model.
Considering alternatives to spinon-mediated witness spin interactions
Although we considered a large number of alternative hypotheses, the only cases capable of explaining the full range of experimental data involved spinon-mediated couplings via spin liquids. In reviewing these hypotheses, there are two strong constraints.
First, any theoretical model with sufficiently rapid variation with distance of the witness spin to witness spin interaction decay (local couplings) will result in a sizable population of isolated witness spins. With Zn site occupation probability \(p\), the percentage of isolated witness spins having no nearest neighbour is \({\left(1-p\right)}^{6}\); with \(p=\,0.33\), this gives 9% (3% of Zn sites). These isolated witness spins must contribute a d.c. magnetic susceptibility, which diverges as \(1/T\) as \(T\to 0\). As shown in Extended Data Fig. 7a, if only 0.7% of Zn sites are occupied by isolated witness spins, this would be enough to show a d.c. susceptibility evolving as \(1/T\) as \(T\to 0\), which is not observed in any of our experiments.
Second, witness spin to witness spin interactions evolving too slowly with distance (\(1/{r}^{2}\) or slower) will lead to an unphysical divergence in the sum over spins forming the structure factor. As well as being unphysical, this situation is incompatible with the experimental inelastic neutron scattering structure factor that shows broad features in momentum space at 2 K, consistent with dominant AF nearest-neighbour correlations17. An example of the structure factor for a slowly decaying spin-wave mediated interaction is shown in Extended Data Fig. 7b,c.
In the context of these constraints, we have considered and ruled out a variety of alternatives to spinon-mediated witness spin interactions including nearest-neighbour local exchange, next-nearest-neighbour local exchange, direct dipolar witness spin interactions, dimer correlations mediating the witness spin interaction, random singlets and random spin clusters, and spin-wave mediating the witness spin interaction (Supplementary Discussion).
The model in the main text discusses the Z2\(\left[0,\,{\rm{\pi }}\right]\beta\) QSL scenario, whereas another candidate kagome QSL in herbertsmithite is the U(1)[0, \({\rm{\pi }}\)] state with a Dirac nodal spinon Fermi surface. We calculated its spinon band structure using the mean field decoupling of ref. 28 and then the witness spin interaction using the same methods we used for Z2: in this case, the calculation becomes that of an Ruderman-Kittel-Kasuya-Yosida interaction between witness spins mediated by the spinon Dirac nodal Fermi surface. We find that all witness couplings are AF, decaying approximately as \(1/{r}^{3}\). The predictions of the witness spin noise, susceptibility and order parameters via these U(1) QSL spinon-mediated witness spin interactions are shown in Extended Data Fig. 8. Each panel can be compared with Figs. 4 and 5. The U(1) model predicts a transition at 110 mK, which is smaller than the value of 150 mK of the Z2 QSL prediction. Thus, our Z2 QSL model is more consistent with the experiment. However, the U(1) QSL model also reproduces the qualitative features of the experiment in Figs. 2 and 3, and cannot be fully excluded using the existing data.