Emergence of the spin polarized domains in the kagome lattice Heisenberg antiferromagnet Zn-barlowite (Zn_0.95Cu_0.05)Cu₃(OD)₆FBr

Abstract

Kagome lattice Heisenberg antiferromagnets are known to be highly sensitive to perturbations caused by the structural disorder. NMR is a local probe ideally suited for investigating such disorder-induced effects, but in practice, large distributions in the conventional one-dimensional NMR data make it difficult to distinguish the intrinsic behavior expected for pristine kagome quantum spin liquids from disorder-induced effects. Here we report the development of a two-dimensional NMR data acquisition scheme applied to Zn-barlowite (Zn_0.95Cu_0.05)Cu₃(OD)₆FBr kagome lattice, and successfully correlate the distribution of the low energy spin excitations with that of the local spin susceptibility. We present evidence for the gradual growth of domains with a local spin polarization induced by 5% Cu²⁺ defect spins occupying the interlayer non-magnetic Zn²⁺ sites. These spin-polarized domains account for ~60% of the sample volume at 2 K, where gapless excitations induced by interlayer defects dominate the low-energy sector of spin excitations within the kagome planes.

Introduction

The quantum spin liquid (QSL) is a novel state of matter formed by entangled spin singlets, in which magnetic frustration effects prevent spins from undergoing magnetic long-range order^1,2,3,4,5. The last few decades have seen concerted efforts to identify model materials of the QSL. However, many of them undergo magnetic long-range order or spin freezing at low temperatures. Copper hydroxyhalide materials herbertsmithite ZnCu₃(OH)₆Cl₂^{6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23} and Zn-barlowite ZnCu₃(OH)₆FBr^{21,24,25,26,27,28,29,30} are among the few exceptions without magnetic symmetry breaking even near absolute zero. As shown in Fig. 1, Cu²⁺ ions with electron spin-1/2 in Zn-barlowite form a corner-sharing triangular lattice known as the kagome lattice. The nearest-neighbor Cu²⁺-Cu²⁺ antiferromagnetic super-exchange interaction is as large as J ~160 K^7,25 and geometrically frustrated. While earlier studies point toward the realization of a proximate kagome QSL state in these materials^{6,7,8,11,12,15,16,17,19,20,21,21,22,23,24,27}, the structural disorder in these materials^14,26 complicates the interpretation of the experimental results. Accordingly, the exact nature of the kagome planes has been highly controversial, necessitating a clear understanding of the influence of disorder on the kagome planes.

Fig. 1: Kagome lattice.

a The crystal structure of Zn-barlowite “Zn_0.95” (Zn_0.95Cu_0.05)Cu₃(OD)₆FBr. For clarity, D-sites attached to O-sites are not shown in any of the panels. b A c-axis view of the kagome plane of Zn_0.95, and the interlayer Zn (cyan), Br (blue), and F (green) sites above and below the kagome plane. c A wider field of view of the kagome plane in Zn_0.95, without D, O, and Br sites. The Zn²⁺ site surrounded by a black circle near the middle of the panel represents the Cu²⁺ interlayer defect occupying one of the twenty Zn sites within this field of view with a 5% probability. The circle with orange-to-yellow shading represents the spin-polarized domains, which encompasses 12 of the 20 ¹⁹F sites in this field of view, corresponding to 60% as observed at 2 K in panel d. All the kagome Cu sites in this domain have polarized spins or interact with them. d Temperature dependence of the fraction F_para of the ¹⁹F sites belonging to the spin-polarized domains estimated in this work. Blue shading depicts the freezing of the lattice distortion starting at ~160 K and completing at ~50 K²¹.

Recent work, including site-selective anomalous X-ray scattering experiments, established that the interlayer non-magnetic Zn²⁺ sites are occupied by extra Cu²⁺ defect spins with ~15% probability in herbertsmithite¹⁴ and ~5% probability in the most Zn-rich Zn-barlowite, synthesized as a powder²⁶. Accordingly, their actual chemical formula is (Zn_0.85Cu_0.15)Cu₃(OH)₆Cl₂ for herbertsmithite and (Zn_0.95Cu_0.05)Cu₃(OH)₆FBr for Zn-barlowite. We will use the acronym Zn_0.95 to represent Zn-barlowite hereafter. One important difference between Zn-barlowite and herbertsmithite is in the location of the interlayer Cu²⁺ defects: in herbertsmithite, they occupy the centered, octahedral Zn²⁺ site, while in Zn-barlowite they occupy an off-center, trigonal prismatic site that is distinct from the octahedral Zn²⁺ site^25,26.

The extra interlayer Cu²⁺ defect spins in both of these materials could directly contribute to the magnetic response of the sample, such as enhanced bulk spin susceptibility χ_bulk at low temperatures. In addition, each interlayer defect spin interacts with the three nearest-neighbor (nn) kagome Cu²⁺ sites located in each of the two adjacent kagome planes. Furthermore, the structural distortion in the vicinity of the interlayer defects^11,17,21,30 might enhance local spin susceptibility χ_local at low temperatures through Dzyaloshinskii-Moriya interaction^9,10,13,19 and/or randomness introduced to the exchange interaction^31,32,33. On the other hand, the upper bound of the concentration of the non-magnetic Zn²⁺ inter-site defects diluting the magnetic Cu²⁺ sites within the kagome planes is ~1% in both of these materials^14,26. This contradicts the persisting claim that as much as ~5% of kagome Cu sites may be diluted by non-magnetic inter-site Zn²⁺ defects, which cause local spin singlets accompanied by oscillatory spin density induced in their vicinity^12,20.

Our earlier NMR measurements in the isotope-enriched single crystal samples of herbertsmithite successfully separated the main NMR peaks from those associated with the nn and the next nearest-neighbor (nnn) ²D (nuclear spin I = 1)¹⁵ and ¹⁷O (nuclear spin I = 5/2)¹⁷ sites of the 15% interlayer Cu²⁺ defect spins. We demonstrated that the Knight shift of the main NMR peak arising from the sites far away from the interlayer defects decreases toward zero. The implication of the diminishing local spin susceptibility χ_local far from the defects has been the subject of an intense debate, i.e., whether the proximate kagome QSL realized in herbertsmithite and Zn_0.95 is gapped or gapless^{11,12,15,17,20,21,24}. On the other hand, the NMR Knight shift at the nn sites of the interlayer defects in herbertsmithite is negative and its magnitude grows with decreasing temperature roughly in proportion to the bulk averaged spin susceptibility χ_bulk^15,17. This means that the interlayer defect moments, which are polarized by the external magnetic field B_ext applied to conduct NMR measurements, induce a negative hyperfine magnetic field B_hyp pointing in the opposite direction from B_ext. In principle, one can probe the spatial extent of the influence of interlayer Cu²⁺ defect spins on the adjacent kagome planes by investigating the nnn and further ²D and ¹⁷O sites. But we were unable to resolve these peaks at low temperatures, because the splitting between the multiple quadrupole-split NMR peaks is obscured by extreme magnetic line broadening with both positive as well as negative B_hyp.

Very recently, we used inverse Laplace transform (ILT) T₁ analysis technique^{21,34,35,36,37,38,39,40,41} to successfully deduce the density distribution function P(1/T₁) of the nuclear spin-lattice relaxation rate 1/T₁ in herbertsmithite and Zn_0.95²¹. The integral of P(1/T₁) is normalized to 1, and P(1/T₁) represents the probability for the nuclear spin to relax with a given value of 1/T₁. In general, 1/T₁ probes low-energy spin excitations in the form of the wave vector q integral of the dynamical electron spin susceptibility at the NMR frequency ω( = 2πf),

$$\frac{1}{{T}_{1}}=\frac{{\gamma }_{n}^{2}{k}_{B}T}{{\mu }_{B}^{2}{\hslash }^{2}}\mathop{\sum}\limits_{{{{\bf{q}}}}}| {A}_{{{{\bf{q}}}}}{| }^{2}\frac{{\chi }^{{\prime\prime} }({{{\bf{q}}}},\omega )}{\omega },$$

(1)

where A_q is the q dependent hyperfine form factor, and χ^″(q, ω) is the imaginary part of the dynamical electron spin susceptibility (i.e., spin fluctuations at the resonance frequency ω).

It turned out that the ⁶³Cu nuclear spin-lattice relaxation rate 1/T₁ in Zn_0.95 develops a bimodal distribution below ~30 K due to the gradual emergence of spin singlets with inhomogeneous gap Δ, which has a large distribution ranging from a few K to as high as Δ ~30 K²¹. These Cu sites involved in spin-singlet formation exhibit diminishing 1/T₁ with decreasing temperature, and the absolute upper bound of their volume fraction is F_singlet ~50%. On the other hand, the remaining Cu sites with a fraction of 1 − F_singlet ~50% or greater do not participate in the spin-singlet formation and remain correlated but paramagnetic. These paramagnetic Cu spins exhibit a large and roughly constant 1/T₁ ~Tχ^″ ≃ 10³ s⁻¹, which is typical for gapless excitations of paramagnetic spins coupled by super-exchange interactions, but their origin was not clear.

The aim of this study is to probe the origin, nature, and fraction of these paramagnetic Cu sites based on the ¹⁹F NMR study of Zn_0.95. ¹⁹F has a nuclear spin I = 1/2, lacks nuclear quadrupole moment, and hence ¹⁹F NMR is immune from nuclear quadrupole interaction that splits and complicates ²D and ¹⁷O NMR lineshapes. ¹⁹F NMR, therefore, provides an avenue to probe the influence of interlayer defects on the kagome planes. Nonetheless, earlier ¹⁹F NMR studies based on the conventional one-dimensional NMR approach^21,24 did not provide a clear-cut picture of the nature of defects in Zn_0.95, because the line broadening obscured the distinction between the intrinsic and defect-induced phenomena.

In this paper, we report the development of a two-dimensional NMR data acquisition scheme and its application to ¹⁹F NMR measurements in Zn_0.95. Instead of separately measuring the position by position variations of χ_local and low-energy spin excitations from the distributions of the Knight shift and 1/T₁, respectively, we generate a two-dimensional correlation map C(f; 1/T₁) between them. In essence, this procedure allows us to count the number of ¹⁹F sites that have certain values of the Knight shift K and 1/T₁ (i.e., χ_local and local spin excitations). We demonstrate that interlayer Cu²⁺ defect spins induce distinct spin-polarized domains at low temperatures, in which ¹⁹F sites exhibit similar behavior as in antiferromagnetic barlowite Cu₄(OH)₆FBr with Néel temperature T_N = 15 K (abbreviated as Cu₄ hereafter)^{25,42,43,44,45,46,47,48}. The ¹⁹F sites in these domains of Zn_0.95 detect gapless spin excitations and enhanced χ_local induced by interlayer Cu²⁺ defects spins. The volume fraction of these domains reaches as much as F_para ~60% at 2 K despite the low 5% concentration of the interlayer defects in Zn_0.95, indicating the gradual spatial growth of these domains at low temperatures.

The rest of this article is organized as follows. We will begin the next RESULTS section with a presentation of the one-dimensional NMR data. We will explain how much one can learn from such conventional NMR data with the aid of ILT, and also its limitations. The following sub-section will develop the two-dimensional NMR data acquisition scheme, and discuss how one can de-convolute the ¹⁹F NMR peak into two components with different magnitudes of 1/T₁. In the DISCUSSION section, we will compare Zn_0.95 and Cu₄, and discuss the implications of our findings.

Results

Conventional one-dimensional NMR data and their limitations

We first summarize the conventional one-dimensional ¹⁹F NMR results, and their limitations in elucidating the influence of defects. In Fig. 2a, we present the representative ¹⁹F NMR lineshapes observed for a deuterated, polycrystalline Zn_0.95 sample of (Zn_0.95Cu_0.05)Cu₃(OD)₆FBr in an external magnetic field of B_ext = 2.4 T. The full width at half maximum (FWHM) of the lineshape is as narrow as 26 kHz at 250 K owing to the small nuclear dipole moment of ²D. Upon cooling, the peak frequency ¹⁹f_peak of the lineshape slightly increases down to ~20 K, then decreases toward the bare resonance frequency ¹⁹f₀ = ¹⁹γ_nB_ext ≃ 96.02 MHz marked by the vertical dashed line. ¹⁹f₀ is the resonance frequency expected in the absence of hyperfine interactions with electron spins in the lattice, and the nuclear gyromagnetic ratio of the ¹⁹F nucleus is ¹⁹γ_n/2π = 40.05 MHz/Tesla. We can summarize the change of ¹⁹f_peak using the NMR Knight shift ¹⁹K_peak defined as

$${}^{19}{K}_{{{{\rm{peak}}}}}=\frac{{}^{19}{f}_{{{{\rm{peak}}}}}-{}^{19}{f}_{0}}{{}^{19}{f}_{0}}=\mathop{\sum }\limits_{i=1}^{12}\frac{{A}_{{{{\rm{hf}}}}}}{{N}_{A}{\mu }_{B}}\,{\chi }_{{{{\rm{local}}}}}^{i}+{}^{19}{K}_{{{{\rm{chem}}}}},$$

(2)

where the summation goes over the twelve nn Cu sites of ¹⁹F nucleus with local spin susceptibility ${\chi }_{{{{\rm{local}}}}}^{i}$, N_A is Avogadro’s number, A_hf is the positive hyperfine coupling with these nn Cu sites²⁴, and ¹⁹K_chem(≃0.01%) is a very small and temperature-independent chemical shift²⁴.

Fig. 2: ¹⁹F NMR lineshapes.

Representative ¹⁹F NMR lineshapes of a"Zn_0.95” (Zn_0.95Cu_0.05)Cu₃(OD)₆FBr and b “Cu₄” Cu₄(OD)₆FBr, both measured in B_ext = 2.4 T and normalized for Boltzmann factor. For clarity, the vertical origin is shifted at different temperatures. The lineshapes of Cu₄ below T_N = 15 K is extremely broad due to the static hyperfine magnetic field from ordered Cu moments. The gray vertical dashed line represents the bare resonance frequency ¹⁹f₀ with no Knight shift.

In Fig. 3, we summarize ¹⁹K_peak (▴) in comparison to χ_bulk. Our ¹⁹K_peak results are similar to the earlier reports in ref. ²⁴ and ref. ²¹; ¹⁹K_peak begins to decrease below ~20 K toward zero. An earlier report fitted the results with an activation form ¹⁹K_peak ~exp(−Δ/k_BT) with a gap Δ²⁴, but we found that large experimental uncertainties of ¹⁹K_peak caused by extremely broad lineshapes below 20 K make it rather difficult to rule out a linear temperature dependence expected for the gapless spin liquid picture based on Dirac Fermions⁴⁹. We recall that ¹⁹K_peak does not necessarily decrease to zero, unless all the 12 nn Cu sites of a given ¹⁹F site form spin singlets; but less than half of Cu sites are actually involved in spin-singlet formation in Zn_0.95²¹. In fact, ¹⁹K_peak exhibits a temperature dependence analogous to that of χ_bulk only down to ~50 K, where the latter begins to grow. The broad ¹⁹F NMR lineshapes observed at low temperatures exhibit long tails stretching hundreds of kHz toward both higher and lower frequencies. This means that a majority of ¹⁹F sites sense large, distributed hyperfine magnetic fields B_hyp with both positive and negative signs, and only a small fraction of ¹⁹F sites exhibit diminishing ¹⁹K_peak below 20 K.

Fig. 3: ¹⁹F NMR Knight shifts ¹⁹K.

Summary of the Knight shift measured in B_ext = 2.4 T for Zn_0.95: ¹⁹K_peak as determined at the nominal peak frequency ¹⁹f_peak of the overall lineshape (▴), ¹⁹K_main () and ¹⁹K_para () for the main and paramagnetic side peaks as determined from the deconvolution of C(f; 1/T₁) based on double Gaussian fit. For comparison, we show ¹⁹K measured for the single peak observed for antiferromagnetic Cu₄ (). The vertical error bars specify the maximum possible range of uncertainties. Using the right axis, we also show the bulk magnetization data of χ_bulk measured in 2.4 T for Zn_0.95 (red) and Cu₄ (blue).

If there is no magnetic inhomogeneity, we expect ${\chi }_{{{{\rm{local}}}}}^{i}={\chi }_{{{{\rm{bulk}}}}}$ for all sites, and hence the NMR lineshape would remain narrow. In other words, the broadening of the lineshape indicates that ${\chi }_{{{{\rm{local}}}}}^{i}$ develops a large distribution below ~50 K. Interestingly, this is the same temperature range, where the structural disorder freezes³⁰ as evidenced by the broadening of the ⁷⁹Br NQR lineshapes; the oscillation on the transverse relaxation measurements also uncovered the signature of structural dimer formation below ~50 K³⁰.

We also measured 1/T₁ at the peak frequency ¹⁹f_peak. The magnitude of 1/T₁ develops a significant distribution below ~160 K, where the lattice freezing sets in²¹. Accordingly, we resorted to deduce the stretched-fit value of $1/{T}_{1,{{{\rm{str}}}}}$ based on the conventional stretched fit of the recovery curve M(t) as follows:

$$M(t)={M}_{0}-A{e}^{-{(t/{T}_{1,{{{\rm{str}}}}})}^{\beta }},$$

(3)

where M₀ is the saturated nuclear magnetization, A is the change in the nuclear magnetization immediately after the inversion π pulse or saturation comb pulses, and β is the phenomenological stretched exponent. See Fig. 1g in ref. ²¹ for the similar results of M(t) and the fit. If magnetic inhomogeneity is negligibly small in Zn_0.95, we expect β = 1. We emphasize that the stretched fit analysis of M(t) is convenient but largely empirical, because it implicitly assumes a specific functional form for the density distribution function P(1/T₁) of 1/T₁^{50,51,52,53,54}. Strictly speaking, it is justifiable only if P(1/T₁) may be represented with the modified Bessel function of the second kind⁵⁴. Empirically, $1/{T}_{1,{{{\rm{str}}}}}$ is often a reasonable approximation of the center of gravity 1/T_1,cg of the true distributed values of 1/T₁ as determined from ILT^{21,38,39,40,55}, including the present case.

We summarize the temperature dependences of $1/{T}_{1,{{{\rm{str}}}}}$ and β observed at ¹⁹f_peak in Fig. 4a, b. These results are very similar to our earlier work conducted in a lower magnetic field 0.72 T²¹. $1/{T}_{1,{{{\rm{str}}}}}$ smoothly decreases with temperature, and does not exhibit an activation behavior expected for fully and homogeneously gapped QSL. On the other hand, the power law behavior expected for the gapless QSL formed by Dirac Fermions⁴⁹ is not observed, either. A key consideration in interpreting the 1/T₁ results is that, again, ¹⁹F sites are subjected to the fluctuations of the transferred hyperfine fields B_hyp from 12 nn Cu²⁺ sites. Since the absolute upper bound of the fraction of spin singlets occupying the Cu sites is F_singlet ~50% in Zn_0.95²¹, many of the ¹⁹F sites are under the influence of paramagnetic Cu spins that are not involved in spin singlets. This explains why the results for $1/{T}_{1,{{{\rm{str}}}}}$ as well as χ_local do not meet the naive expectations for purely gapped or purely gapless QSL.

Fig. 4: Conventional stretched fit results of $1/{T}_{1,{{{\rm{str}}}}}$ in Zn_0.95.

a Temperature dependence of $1/{T}_{1,{{{\rm{str}}}}}$ observed at ¹⁹f_peak for Zn_0.95. For comparison, we also show the center of gravity 1/T_1,cg of the distributed 1/T₁ as determined from ILT. b Temperature dependence of the stretched exponent β observed at the peak. c Point-by-point ¹⁹F NMR lineshape at 30 and 4.2 K. d, e$1/{T}_{1,{{{\rm{str}}}}}$ and β measured at each data point in c. Notice that $1/{T}_{1,{{{\rm{str}}}}}$ is enhanced at both higher and lower frequencies. β deviates from one more significantly near the center of the lineshapes, because the fast and slow contributions are superimposed.

We also measured $1/{T}_{1,{{{\rm{str}}}}}$ at various frequencies within the extremely broad lineshape. We summarize the representative results of the frequency dependence of $1/{T}_{1,{{{\rm{str}}}}}$ and β in Fig. 4c–e. $1/{T}_{1,{{{\rm{str}}}}}$ is larger at both lower and higher frequency sides of the peak. This suggests that the 1/T₁ relaxation process at ¹⁹F sites in the tailed sections of the broad lineshape is enhanced by low-energy spin excitations associated with the source of the hyperfine magnetic fields that gives rise to the line broadening. But we are unable to pinpoint their origin based on these measurements.

Two-dimensional NMR data C(f; 1/T ₁)

Proximate quantum spin liquid materials that do not undergo magnetic long-range order are often structurally disordered and exhibit significantly distributed 1/T₁. In our recent work, we overcame the difficulties in conventional NMR data acquisition protocols outlined in the previous section by using the ILTT₁ analysis technique^38,39. In the ILT analysis of the distributed 1/T₁, one only assumes that 1/T₁ has a spatial distribution, and each nuclear spin relaxes exponentially with the jth value 1/T_1j with the probability density P(1/T_1j). Utilizing Tikhonov regularization, one can optimize the fit of the experimentally observed recovery curve M(t) with a summation of single exponentials with distributed 1/T₁,

$$M(t)=\mathop{\sum}\limits_{j}{M}_{0}\left[1-2{e}^{-t/{T}_{1j}}\right]P(1/{T}_{1j}),$$

(4)

where the summation of the probability P(1/T_1j) is normalized to 1. See refs. ^36,38,39 for the details of the numerical ILT procedures.

In Fig. 5a, we summarize representative P(1/T₁) results deduced for the distribution of 1/T₁ in Zn_0.95 measured at ¹⁹f_peak in B_ext = 2.4 T. These P(1/T₁) results are similar to our previous findings in B_ext = 0.72 T²¹, and indicate that 1/T₁ has a roughly symmetrical distribution above ~100 K. The center of gravity 1/T_1,cg of the distribution of P(1/T₁) summarized in Fig. 4a is similar to $1/{T}_{1,{{{\rm{str}}}}}$, as noted above. Below ~100 K, P(1/T₁) begins to develop a tail toward larger values of 1/T₁. This means that some of the ¹⁹F nuclear spins resonating around ¹⁹f_peak have enhanced 1/T₁, but the majority of ¹⁹F nuclear spins around ¹⁹f_peak still relax with the slower values comparable to $1/{T}_{1,{{{\rm{str}}}}}$. We emphasize that the conventional stretched fit analysis alone cannot clarify such qualitative changes in the distribution of 1/T₁.

Fig. 5: Density distribution function P(1/T₁).

a Representative results of P(1/T₁) deduced at the peak frequency ¹⁹f_peak of the ¹⁹F NMR linehapes observed for Zn_0.95. The integral of P(1/T₁) is normalized to 1 (one). b Corresponding results observed for Cu₄. Notice that the side peak of P(1/T₁) observed for Zn_0.95 with 1/T₁ ~300 s⁻¹ is absent for Cu₄. The large distribution observed at 4.2 K is caused by antiferromagnetic long-range order. We fixed the Tikhonov regularization parameter for ILT as α = 1 to achieve the same effective resolution. The distributions of 1/T₁ is negligibly small with β ≃ 1 above ~160 K for Zn_0.95 and above T_N = 15 K for Cu₄, hence the apparent width of P(1/T₁) in these regions is set primarily by the effective resolution of our ILT with fixed α = 1. (See Supplementary Materials of refs. ^21,38 for details about the effective resolution of the ILT analysis.).

Successful demonstration of the distribution P(1/T₁) of 1/T₁ in Zn_0.95²¹ was a major step toward clarifying the nature of the disorder in kagome materials. However, we were unable to identify the origin of the ¹⁹F sites with enhanced 1/T₁, nor can we estimate their population based on such a one-dimensional NMR data set alone, even with ILT. Notice that we were measuring two sets of one-dimensional distributions separately in these measurements: (a) the distribution of the resonant frequency ¹⁹f and hence ¹⁹K in the form of the NMR lineshape, which provides us with the histogram of the distribution of χ_local averaged over 12 Cu²⁺ sites nn to each ¹⁹F site; (b) the distribution of 1/T₁ in the form of P(1/T₁) at the peak frequency ¹⁹f_peak. What is still missing is the correlation between the distributions in static local spin susceptibility χ_local and local spin dynamics reflected on 1/T₁.

Here, we take one more step forward to clarify the nature of the disorder by advancing the ILTT₁ analysis technique. In addition to P(1/T₁) at the peak frequency ¹⁹f_peak of the lineshape, we also measure M(t) and use ILT to deduce P(1/T₁) at other frequencies within the broadened lineshape, and thereby correlate the histograms of ¹⁹K and 1/T₁ distributions. In Fig. 6, we illustrate how this two-dimensional NMR data acquisition scheme works. The top panel Fig. 6a shows the conventional overall NMR lineshape at 30 K. The broad lineshape peaked at ¹⁹f_peak ≃ 96.06 MHz represents the histogram of the total hyperfine magnetic field at ¹⁹F sites transferred from the 12 nn Cu sites. The result of P(1/T₁) measured at ¹⁹f_peak, which we now represent as C(¹⁹f_peak; 1/T₁) to show the specific frequency, is shown in the bottom panel Fig. 6c underneath the peak of the total lineshape in Fig. 6a. Note that this is the same P(1/T₁) curve at 30 K presented in Fig. 5a, rotated clockwise by 90^∘. The integrated area of the C(¹⁹f_peak; 1/T₁) as a function of 1/T₁ is normalized to 1, as before.

Fig. 6: Two-dimensional correlation map C(f; 1/T₁) of Zn_0.95 and 1/T₁-resolved NMR lineshapes.

a Conventional one-dimensional NMR lineshape of all ¹⁹F nuclear spins observed for Zn_0.95 at 30 K. b Example of the color contour map of the frequency f vs. 1/T₁ correlation function C(f; 1/T₁) at 30 K deduced from the ILT measurements of P(1/T₁) at various frequencies within the lineshape. A side peak marked as Para emerges at low temperatures centered around f = 96.02 MHz or below, with enhanced 1/T₁ ≃ 100 s⁻¹. c Vertical slice of the color contour map of C(f; 1/T₁) at a fixed frequency f is nothing but the one-dimensional P(1/T₁) result at that frequency. For example, C(f; 1/T₁) sliced vertically at ¹⁹f_peak = 96.06 MHz is P(1/T₁) at 30 K presented in Fig. 5a. d Horizontal slice of C(f; 1/T₁) at a fixed value of 1/T₁ yields the 1/T₁-resolved NMR lineshape with the contributions from only the ¹⁹F nuclear spins with that particular value of 1/T₁. Note that the ¹⁹F NMR 1/T₁-resolved lineshape for 1/T₁ = 300 s⁻¹ is peaked at 96.02 MHz, which is lower than the main peak frequency ¹⁹f_peak ≃ 96.06 MHz. The Knight shift ¹⁹K_para arises from the former.

Next, we repeat the measurement of the recovery curves M(t) at different frequencies f within the broad lineshape, and deduce the corresponding C(f; 1/T₁) between f = 95.9 and 96.15 MHz, as also shown in Fig. 6c. The integrated area of C(f; 1/T₁) at each frequency is normalized to the relative intensity of the one-dimensional NMR lineshape in Fig. 6a. For example, the overall ¹⁹F signal intensity at 96.02 MHz in Fig. 6a is only ~50% of the peak intensity at ¹⁹f_peak ≃ 96.06 MHz, hence the integrated area of C(96.02 MHz; 1/T₁) is normalized to ~0.5. Once we deduce multiple C(f; 1/T₁) curves across the broad NMR lineshape in panel Fig. 6c, we can generate the two-dimensional color contour plot of C(f; 1/T₁) as a function of both f and 1/T₁ and normalize its integral in the two-dimensional parameter space to 1. We present the normalized color contour plot of C(f; 1/T₁) in panel Fig. 6b. C(f; 1/T₁) represents the two-dimensional correlation map between the distributed resonance frequency f and the distributed low-energy spin excitations reflected on 1/T₁. Therefore, the magnitude of C(f; 1/T₁) at specific values of f and 1/T₁ represents the relative probability of a ¹⁹F nuclear spin to have particular values of f and 1/T₁.

These two-dimensional NMR data of C(f; 1/T₁) require an order of magnitude longer data acquisition time in comparison to the conventional one-dimensional NMR data, because we need to repeat 1/T₁ measurements at various frequencies with high precision required by ILT. But the results of the correlation map C(f; 1/T₁) provide much richer information. For example, we can slice the color contour map in Fig. 6b horizontally at a given value of $1/{T}_{1}^{{\prime} }$, and plot $C(f;1/{T}_{1}^{{\prime} })$ as a function of f, as summarized in the right panel Fig. 6d. Then the resulting curves represent the NMR lineshapes for nuclear spins with the specific values of $1/{T}_{1}^{{\prime} }$. Unlike the conventional NMR lineshape in panel Fig. 6a measured from the sum of all ¹⁹F nuclear spins with various 1/T₁ values, the NMR lineshapes in Fig. 6d are resolved for different 1/T₁ values. We opt to coin them as ILTT₁ resolved NMR lineshapes. Notice that the ILTT₁ resolved NMR lineshape with 1/T₁ ≃ 300 s⁻¹ is centered around 96.02 MHz (purple curve in Fig. 6d), whereas the majority of the ¹⁹F sites are concentrated around the main peak near 96.05 MHz with much slower values of 1/T₁ ≃ 30 s⁻¹ (green curve in Fig. 6d). To the best of our knowledge, this is the first example of successfully resolving the inhomogeneously broadened NMR lineshape for different values of 1/T₁, although narrow NMR lineshapes with homogeneous broadening were previously resolved with the aid of fast Fourier transform⁵⁶.

Temperature dependence of C(f; 1/T ₁)

In Fig. 7a–d, we summarize the representative results of C(f; 1/T₁) observed for Zn_0.95 at various temperatures. See Supplementary Fig. 1 for the complete set of C(f; 1/T₁) data and Supplementary Movie 1 for the movie summary of the evolution of C(f; 1/T₁) as a function of temperature. Notice that C(f; 1/T₁) has a nearly symmetrical oval shape around 250 K, because both f and 1/T₁ have only a minor distribution centered around their central value. This main peak remains fairly narrow along both the f and 1/T₁ axes, at least down to 10 K. But a distinct component emerges below ~120 K on the upper left-hand side of the main peak, and gradually splits off. In what follows, we call the emergent side peak with larger values of 1/T₁ as the paramagnetic peak for the reasons to become clear below. The two-dimensional distribution of the paramagnetic peak continuously grows along the f axis toward 2 K, but the growth of the distribution along the 1/T₁ axis appears less significant at low temperatures.

Fig. 7: Evolution of C(f; 1/T₁) for Zn_0.95.

a–d Representative C(f; 1/T₁) results at various temperatures observed for Zn_0.95. The vertical dashed line marks the zero Knight shift frequency ¹⁹f₀. See Supplementary Materials for the complete set of data as well as the animation (movie) of the temperature evolution of C(f; 1/T₁). e–h C(f; 1/T₁) at corresponding temperatures, deduced from the summation of the two Gaussian functions in the two-dimensional f and 1/T₁ space that best fit the experimental results in a–d. Note that above 100 K, the two-dimensional gaussian does not capture the behavior of the side peak as it is merged into the main peak.

In order to separate the contributions of the main and paramagnetic peaks, we deconvoluted C(f; 1/T₁) into two separate two-dimensional Gaussian functions of f and 1/T₁ by conducting the least χ² fit in the two-dimensional parameter space. See Supplementary Fig. 2 and Supplementary Fig. 3 for the details of the fit. We show the color contour map of the summation of these two Gaussian components in Fig. 7e–h. The two-component fit reproduces the experimental data in Fig. 7a–d fairly well. The paramagnetic side peak merges into the main peak at higher temperatures, and we were unable to clearly resolve them above 100 K. In Fig. 1d, we summarize the temperature dependence of the fraction F_para of the ¹⁹F sites involved in the paramagnetic side peak estimated from the integral of the two-dimensional Gaussian peak in the f-1/T₁ space. F_para gradually grows with decreasing temperature, and reaches as large as F_para ~0.6 at 2 K. We emphasize that the faint signature of the paramagnetic side peak seems to persist even at 150 K or above, as shown in Supplementary Fig. 1. It is consistent with our earlier finding that the ²D and ¹⁷O sites nn to the interlayer defects are distinct even at 295 K in herbertsmithite^15,17. On the other hand, the estimation of F_para was difficult above 100 K due to the unstable nature of the two-dimensional fit of the small side peak involving a large number of free parameters.

From the central values of the two-dimensional Gaussian function for the paramagnetic peak, we estimated the central values of the distributed ¹⁹K_para and $1/{T}_{1}^{{{{\rm{para}}}}}$, whereas ¹⁹K_main and $1/{T}_{1}^{{{{\rm{main}}}}}$ were estimated from the two-dimensional Gaussian function of the main peak; we summarize these results in Figs. 3, 8. The paramagnetic peak is somewhat asymmetrical along the frequency axis with slightly more weight along the lower frequency side, but this does not significantly affect our estimation of the central value of ¹⁹K_para. The central values of ¹⁹K_main and $1/{T}_{1}^{{{{\rm{main}}}}}$ for the main peak are close to the nominal peak values of the one-dimensional data discussed in the previous section.

Fig. 8: Temperature dependences of 1/T₁.

The central value $1/{T}_{1}^{{{{\rm{main}}}}}$ () of the main peak and $1/{T}_{1}^{{{{\rm{para}}}}}$ () of the paramagnetic peak observed for Zn_0.95, as estimated from the peak of the double Gaussian fit of C(f; 1/T₁). The vertical error bars specify the maximum possible range of uncertainties. For comparison, we also show $1/{T}_{1,{{{\rm{str}}}}}$ observed at the nominal peak frequency ¹⁹f_peak of Zn_0.95 (▴, same stretched fit results as in Fig. 4a) and Cu₄ ().

Discussions

Our two-dimensional correlation maps C(f; 1/T₁) and ILTT₁ resolved ¹⁹F NMR lineshapes revealed that two distinct types of ¹⁹F sites exist in Zn_0.95 at low temperatures: (i) ¹⁹F sites with characteristic signatures of paramagnetic spins, i.e. enhanced Knight shift ¹⁹K_para due to growing spin polarization, and large temperature-independent $1/{T}_{1}^{{{{\rm{para}}}}}$ induced by gapless low-energy spin excitations obeying χ^″ ~1/T. The volume fraction F_para of these ¹⁹F sites reaches as much as ~60% at 2 K, indicating that spin-polarized domains gradually extends in space, as schematically shown in Fig. 1c. (ii) The main intrinsic ¹⁹F sites that are mostly immune from these effects, suggesting that they are spatially more distanced from the source of the spin polarization and more reflective of the intrinsic behavior of the kagome planes. These two different types of ¹⁹F sites become indistinguishable above ~120 K within our experimental resolutions.

Interestingly, the negative Knight shift ¹⁹K_para of the paramagnetic peak shows qualitatively the same behavior as the negative Knight shift observed at the nn ²D and ¹⁷O sites of the 15% interlayer Cu²⁺ defect spins occupying the Zn²⁺ sites in herbertsmithite^15,17. This strongly suggests that the paramagnetic ¹⁹F sites are also located in the vicinity of the 5% interlayer Cu²⁺ defects of Zn_0.95. In order to test this scenario, we repeated two-dimensional ¹⁹F NMR measurements in a deuterated powder sample of barlowite Cu₄(OD)₆FBr. Cu₄ is the parent antiferromagnetic phase of Zn_0.95 with Cu²⁺ spins occupying all the interlayer Zn²⁺ sites in Zn_0.95, and undergoes a Néel transition into an antiferromagnetically ordered state at T_N = 15 K^25,47,48.

We summarize the ¹⁹F NMR results of Cu₄ in Figs. 2b, 3, 5b, 8, 9 in comparison to the results for Zn_0.95. See Supplementary Fig. 4 for the complete set of C(f; 1/T₁) data of Cu₄ and Supplementary Movie 2 for the movie summary of the evolution of C(f; 1/T₁) as a function of temperature. Our one-dimensional NMR data for Cu₄ are similar to earlier ¹⁹F NMR measurements conducted in a limited temperature range below 100 K for a protonated powder sample of Cu₄⁴⁷. Notice that two-dimensional correlation maps C(f; 1/T₁) of Cu₄ in Fig. 9 show no signs of a split-off peak. Instead, Cu₄ exhibits only one type of ¹⁹F NMR peak in the entire paramagnetic state above T_N; this is consistent with the crystal structure of Cu₄ with a crystallographically unique ¹⁹F site.

Fig. 9: Two-dimensional correlation map C(f; 1/T₁) of Cu₄.

a–d The two-dimensional correlation map C(f; 1/T₁) observed for Cu₄ at representative temperatures. The ¹⁹F NMR lineshape gradually broadens below 250 K, as shown in Fig. 2b, toward T_N = 15 K, but the distribution of 1/T₁ is always minimal and roughly symmetrical at all frequencies in the paramagnetic state. Unlike Zn_0.95, there is no hint of a distinct split-off peak above T_N. The drastic change of C(f; 1/T₁) at 4.2 K is due to antiferromagnetic long-range order.

A striking finding here is that ¹⁹K_para observed for Zn_0.95 shows nearly identical behavior as the unique ¹⁹F sites in antiferromagnetic Cu₄ above T_N, as shown in Fig. 3. This contrasts with the behavior of ¹⁹K_main and reinforces our conclusion that the paramagnetic ¹⁹F sites observed for Zn_0.95 are in the close proximity with interlayer Cu²⁺ defects occupying the Zn²⁺ sites, where the structural and magnetic environment is locally similar to Cu₄’s. Since χ_bulk of Cu₄ shows strong enhancement below ~100 K, naturally, one expects analogous enhancement of the local spin susceptibility χ_local in the vicinity of Cu²⁺ interlayer defects of Zn_0.95. The negative sign of the hyperfine magnetic field B_hyp leads to a strongly negative ¹⁹K_para in Zn_0.95. This further suggests that the conclusions from prior one-dimensional ¹⁹F NMR studies that neglect the effects of impurities should be reassessed²⁴.

While our ¹⁹F NMR measurements establish the gradual growth of spin-polarized domains induced in the vicinity of the Cu²⁺ interlayer defects, our measurements do not reveal the exact spin configuration in these domains. The growth of staggering spin polarization within the kagome planes was previously proposed to account for the broadening of ¹⁷O NMR lineshapes in herbertsmithite, but it was attributed to the Zn²⁺ anti-site defects within the kagome planes^12,20. Our findings here indicate that such spin vacancies within the kagome planes are not required to account for the NMR line broadening in these kagome materials; after all, Cu₄ has no Zn²⁺ ions in its composition. Moreover, the upper bound of the Zn²⁺ anti-site defect population is as little as ~1%^14,25. At high temperatures (T > 160 K) dominated by short-range spin fluctuations, it is striking that the behavior of Zn_0.95 (¹⁹K_main and $1/{T}_{1}^{{{{\rm{main}}}}}$ decrease upon cooling) clearly differs from Cu₄ (¹⁹K and 1/T₁ increase upon cooling). This hints at an approach to two different ground states.

Another important point to emphasize is that the fraction of the paramagnetic ¹⁹F sites in Fig. 1d gradually grows below ~100 K to as large as F_para ~60% toward the base temperature of 2 K. In Fig. 1c, we use a circle with orange-to-yellow shading to schematically present the spin-polarized domain surrounding a Cu²⁺ defect spin occupying the interlayer Zn²⁺ site in the middle. Note that 12 ¹⁹F sites in the field of view belong to the domain out of 20 ¹⁹F sites, corresponding to 60%. Assuming that the interlayer defect sites are randomly located without forming a cluster, Fig. 1c suggests that kagome Cu²⁺ sites within the range up to third or fourth nearest-neighbor may be affected by each interlayer Cu²⁺ defect spin at 2 K. As F_para decreases at higher temperatures, the range of the spin-polarized domain smoothly decreases.

In our recent work, we demonstrated, based on ⁷⁹Br NQR measurements, that the crystal structure develops local distortions below ~160 K (~J), and their slow fluctuations freeze below ~50 K³⁰. That is, the crystal structure of Zn_0.95 is locally distorted without breaking the overall ideal kagome rotational symmetry of the lattice. The crystal structure of Cu₄ is also distorted at low temperatures, as evidenced by the development of the superlattice structural peak in neutron scattering experiments²⁵, prior to the drastic growth of χ_bulk below ~100 K. It is possible that these paramagnetic domains in Zn_0.95 are stabilized and grow in size owing to the local lattice distortion.

Turning our attention to the 1/T₁ results in Fig. 8, $1/{T}_{1}^{{{{\rm{para}}}}}\simeq 100$ s⁻¹ observed for the paramagnetic peak of Zn_0.95 is nearly temperature-independent, and comparable to the high-temperature asymptotic value of 1/T₁ observed for Cu₄. In contrast, Cu₄ exhibits a continuous enhancement of 1/T₁ with decreasing temperature due to the growth of short-range order toward the long-range antiferromagnetic order at T_N = 15 K. According to Eq.(1), the local dynamical spin susceptibility in the vicinity of the interlayer Cu²⁺ defects of Zn_0.95, therefore, grows slowly as χ^″ ~1/T. This means that gapless low-energy spin excitations exist locally in the vicinity of interlayer Cu²⁺ defects. We emphasize, however, that these gapless spin excitations in Zn_0.95 are driven by the presence of interlayer defects and their influence on local magnetic and structural environments of the kagome planes, and should not be confused with the intrinsic gapless spin excitations expected for certain theoretical models of pristine kagome spin liquids^49,57.

In this context, it is worth recalling that some gapless spin liquid models predict a power law behavior in the temperature dependence of 1/T₁⁴⁹. The recent ¹⁷O 1/T₁ results²⁰ as well as our initial ⁶³Cu 1/T₁ results¹¹ for herbertsmithite, both deduced from the empirical analysis of M(t) due to extremely large distributions of 1/T₁, appeared to confirm such expectations of a power law behavior. However, our recent ILTT₁ analysis of 1/T₁ observed at ⁶³Cu sites showed that such an apparent power law temperature dependence previously proposed for herbertsmithite is purely fictitious²¹. Cu sites exhibit either gapped spin-singlet behavior with activation type temperature dependence of 1/T₁ or paramagnetic behavior with large, nearly temperature-independent 1/T₁ ~10³ s⁻¹. The latter is qualitatively similar to $1/{T}_{1}^{{{{\rm{para}}}}}$ observed here at paramagnetic ¹⁹F sites. This suggests that the paramagnetic ⁶³Cu sites observed for both Zn_0.95 and herbertsmithite²¹ may be also located in the same spin-polarized domains observed here.

The absence of the signature of critical slowing down or spin freezing at 15 K for $1/{T}_{1}^{{{{\rm{para}}}}}$ implies that the dimensions of the spin-polarized domain remain finite. Assuming that these domains are centered around 5% interlayer defect spins, F_para ~60% at 2 K implies that the domains grow in space but span only a few triangles from defects, as schematically shown in Fig. 1c. The significant fraction F_para ~60% at 2 K with strongly enhanced $1/{T}_{1}^{{{{\rm{para}}}}}$ also implies that, as far as the low-energy spin excitations probed by NMR are concerned, one can expect to find the intrinsic kagome spin liquid behavior associated with spin-singlet formation only in 1 − F_para ~40% or less of the sample volume. In Fig. 10, we show the temperature dependence of 1 − F_para together with the upper bound of the spin-singlet fraction F_singlet as determined by ⁶³Cu NQR measurements²¹. Indeed the upper bound F_singlet in Zn_0.95 is limited by the fraction of non-paramagnetic domains 1 − F_para without spin polarization induced in the vicinity of interlayer Cu²⁺ defects. We recall that F_singlet is an upper bound²¹ and could easily be overestimated by a factor of ~2 due to difficulties in accounting for the transverse relaxation effects, but the magnitude of F_singlet is still significant in comparison to 1 − F_para ~40%.

Fig. 10: Volume fraction of paramagnetic and spin-singlet sites.

1 − F_para (∘, obtained from the results in Fig. 1d) represents the shrinking volume fraction of Zn_0.95, which is not within the spin-polarized domain and hence could exhibit intrinsic kagome spin liquid behavior without the perturbation caused by interlayer Cu²⁺ defects occupying the Zn²⁺ sites. Notice that only ~40% of such a volume remains at 2 K in Zn_0.95. Also shown (♢) is the upper bound of the spin-singlet fraction F_singlet of Cu sites, estimated from the observable ⁶³Cu NQR signals at the fixed pulse separation time τ = 6 μs²¹ (The lower error bars below 6 K represent the lower bound of F_singlet, estimated based on the assumption that the ⁶³Cu signal loss affects only the paramagnetic component. All other error bars correspond to the absolute maxima or minima of F_singlet²¹.) Blue shading schematically shows the temperature range where the freezing of the structural distortion sets in (~160 K) and completes (~50 K).

In conclusion, we have found two distinct magnetic signatures in the quantum spin liquid candidate Zn-barlowite, Zn_0.95. One signature arises from domains that are nearest to interlayer Cu²⁺ impurities, and the local susceptibility and spin-lattice relaxation rate are similar to the paramagnetic state above T_N of the antiferromagnetic parent compound barlowite, Cu₄. The other signature arises from regions far from the impurities and, therefore more closely reflects the intrinsic behavior of the kagome spins, where spin singlets gradually emerge with spatially inhomogeneous gaps²¹. Finally, it may be also interesting to apply the one- and two-dimensional NMR techniques based on ILT to other kagome spin liquid and related materials, which are known to exhibit distributions in their NMR properties^58,59,60. Moreover, our approach might yield fresh insight into the distributions in local magnetic properties in unrelated disordered materials, such as diluted magnets⁵² and heavy Fermions⁵³.

Methods

Samples

We synthesized the deuterated (D = ²H) powder samples of ZnCu₃(OD)₆FBr and Cu₄(OD)FBr based on the procedures described in detail in ref. ²⁵. We confirmed the sample quality based on powder X-ray diffraction measurements.

NMR

We conducted the spin echo ¹⁹F NMR measurements using standard pulsed NMR spectrometers. We conducted most of the 1/T₁ measurements by applying saturation comb pulses prior to the spin echo sequence, and confirmed that the results are the same if we use an inversion pulse instead. The pulse separation time τ between the 90^∘ and 180^∘ pulses were 10 μs to ensure that we do not overlook any signals with faster transverse relaxation time. We carried out the inverse Laplace transform of the recovery curve M(t) based on Tikhonov regularization using a fixed regularization parameter α = 1 rather than optimizing at different temperatures and frequencies. This ensured that the effective resolution of ILT remains unchanged between different sets of measurements.