Spin excitations of a proximate Kitaev quantum spin liquid realized in Cu$_2$IrO$_3$

Magnetic moments arranged at the corners of a honeycomb lattice are predicted to form a novel state of matter, Kitaev quantum spin liquid, under the influence of frustration effects between bond-dependent Ising interactions. Some layered honeycomb iridates and related materials, such as Na$_{2}$IrO$_{3}$ and $\alpha$-RuCl$_{3}$, are proximate to Kitaev quantum spin liquid, but bosonic spin-wave excitations associated with undesirable antiferromagnetic long-range order mask the inherent properties of Kitaev Hamiltonian. Here, we use $^{63}$Cu nuclear quadrupole resonance to uncover the low energy spin excitations in the nearly ideal honeycomb lattice of effective spin $S = 1/2$ at the Ir$^{4+}$ sites in Cu$_{2}$IrO$_{3}$. We demonstrate that, unlike Na$_{2}$IrO$_{3}$, Ir spin fluctuations exhibit no evidence for critical slowing down toward magnetic long range order in zero external magnetic field. Moreover, the low energy spin excitation spectrum is dominated by a mode that has a large excitation gap comparable to the Ising interactions, a signature expected for Majorana fermions of Kitaev quantum spin liquid.


I. INTRODUCTION
Quantum spin liquids (QSL's) are a liquid-like state of matter, in which entangled spins remain paramagnetic without undergoing a magnetic long-range order. Identifying a QSL is the holy grail of quantum condensed matter research today [1][2][3][4]. A recent milestone in the search of QSL's is the rigorous theoretical solution of Kitaev QSL for effective spin S = 1/2 arranged in a honeycomb lattice [5]. The Kitaev Hamiltonian may be written as where J γ K (γ = x, y, and z) represents the bonddependent, frustrated Ising interaction as schematically shown in Fig. 1A, and the lattice summation < i, j > is taken between nearest-neighbor spins. Kitaev proved based on analytic calculations that his Hamiltonian supports a quantum spin liquid ground state.
Besides possessing the ground state with no magnetic long range order, ideal Kitaev QSL is known to exhibit a variety of exotic properties [5][6][7][8][9]. For example, fractionalization of quantum spins leads to two types of Majorana fermions, itinerant Majorana fermions and localized fluxes. Therefore, elementary excitations in Kitaev QSL are not conventional magnons (spin waves), and may exhibit an inherent gap in its spin excitation spectrum, depending on the anisotropy of the Ising interactions * imai@mcmaster.ca in eq.(1). Thermodynamic properties could show nontrivial behavior, too. For example, the uniform spin susceptibility χ spin is expected to increase or saturate with decreasing temperature, even if the excitation spectrum of Majorana fermions are inherently gapped [7][8][9]; this is due to the non-conservative nature of the z-component of the total spin, S z .
The Ir 4+ , Rh 4+ , and Ru 3+ ions in Na 2 IrO 3 [10], Li 2 RhO 3 [11], α-RuCl 3 [12], and more recently discovered H 3 LiIr 2 O 6 [13] are under the influence of spin-orbit interactions, possess effective spin S i = 1/2, and interact with each other primarily via frustrated Ising interactions [14]. Therefore, the honeycomb planes formed by these spins have been proposed to harbor Kitaev QSL [14][15][16]. In reality, additional interactions, such as Heisenberg's exchange term H H = <i,j> J H S i · S j and the symmetric off-diagonal exchange term H Γ = <i,j> Γ(S α i S β j + S β i S α j ) [3,14,15] perturb Kitaev Hamiltonian in eq.(1), leading to more complicated Kitaev-Heisenberg Hamiltonian, Theoretical calculations established that Kitaev-Heisenberg Hamiltonian has a long-range ordered ground state rather than the spin liquid state in a broad parameter space [3]. Experimentally, Na 2 IrO 3 [10], the first model material proposed for Kitaev-Heisenberg model [15], indeed has a zig-zag antiferromagnetic longrange order below T N ∼ 17 K, as shown schematically in Fig. 1B [17,18]. Such a long-range ordered ground state conceals the inherent properties of Kitaev QSL. In particular, bosonic spin-wave excitations (with total spin 1) associated with the undesirable antiferromagnetic long-range order dominate the low energy sector of spin excitation spectrum in Na 2 IrO 3 [18] and α-RuCl 3 [19][20][21], as shown in the left panel of Figure 2A. This makes it difficult to identify the intrinsic excitations expected for eq.(1), i.e. the fractionalized Majorana fermions (with total spin 1/2) and pairs of gauge fluxes (with total spin 1). Untangling the influence of H H and H Γ on H K has been at the focus of intense research both theoretically and experimentally, but there have been only a handful of proximate Kitaev   Cu + ions into the Na + sites of Na 2 IrO 3 paved a new avenue to investigate Kitaev QSL. We present the crystal structure of Cu 2 IrO 3 in Figure 1D-E [22]. Both the Ir-Ir-Ir (118.7 • to 122.5 • ) and Ir-O-Ir (95 • to 98 • ) bond angles are closer to the ideal Kitaev geometry of 120 • and 90 • [14,15], respectively, than those of Na 2 IrO 3 (114.9 • to 124.2 • and 98 • to 99.4 • , respectively [10]). There are two distinct types of Cu sites in Cu 2 IrO 3 : Cu(H) sites are located within the honeycomb plane in the middle of each hexagon formed by six Ir sites, while inter-planar Cu(I) sites sit halfway between honeycomb planes (i.e. Cu(I) 1.5 Cu(H) 0.5 IrO 3 ). The uniform magnetic susceptibility χ observed for Cu 2 IrO 3 is nearly identical with that of Na 2 IrO 3 except below ∼30 K, where defect spins dominate χ in Cu 2 IrO 3 [22,23]. In the case of Na 2 IrO 3 , The theoretical fit of χ based on Kitaev-Heisenberg Hamiltonian was very good, and led to the initial estimate of J K ∼ 21 meV and J H ∼ −4 meV [15]. In order to account for the zig-zag ordered ground state, however, the additional off-diagonal exchange terms in eq.(2) are required [3]. Subsequent first principles calcula-tions based on exact diagonalization, quantum chemistry method, and perturbation theory led to J K = −16.8 ∼ −29.4 meV, J H = 0.5 ∼ 3.2 meV, and Γ ∼ 1 meV [24]. In view of the nearly identical χ, we expect that spin-spin interactions in Cu 2 IrO 3 are comparable to these estimations. Unlike Na 2 IrO 3 , however, Cu 2 IrO 3 shows no evidence for antiferromagnetic long-range order [22], and recent µSR [23,25] as well as our NQR measurements presented below establish that only ∼50% of the sample volume exhibits spin freezing below the onset temperature of T f ∼ 10 K. Moreover, muons exhibit no spin precession about static hyperfine field, and Ir spins remain dynamic even below 1 K [23,25]. This means that the low energy sector of the spin excitation spectrum in Cu 2 IrO 3 is largely immune from the spin-waves, as schematically shown in the right panel of Figure 2A, and opens a wide temperature window above T f to probe the inherent properties of Kitaev QSL.
On the other hand, Cu 2 IrO 3 has its own complications arising from stacking faults, as is often the case for previously identified Kitaev candidate materials, and from the aforementioned defect spins [25]. The defect-induced enhanced susceptibility below ∼ 30 K is known to obey a peculiar scaling law [23], reminiscent of analogous scaling behavior reported earlier for other spin liquid candidate materials with disorder [26]. The origin and nature of the defect spins in Cu 2 IrO 3 as well as in other related iridates [13,[27][28][29] are a matter of intense research and yet to be clarified. In the case of Cu 2 IrO 3 , XANES measurements showed that up to ∼1/3 of 63 Cu(H) sites are occupied by Cu 2+ ions rather than Cu + ions [25]. To maintain charge neutrality, the charge state of an Ir 4+ ion(s) in their vicinity may change, possibly to Ir 3+ with null spin. If randomly distributed, they might create bond disorder [30] and/or site dilution in the Kitaev Hamiltonian, in addition to the decoration of the Kitaev lattice with additional spin S = 1/2 at the Cu 2+ sites. Curiously, the bulk susceptibility of Cu 2 IrO 3 and Na 2 IrO 3 shares nearly identical Weiss temperature Θ w = −110 ∼ −123 K [22], whereas Mg 2+ substitution into Na 2 IrO 3 is known to suppress Θ w quickly to almost zero [29].
In this paper, we will use 63 Cu nuclear quadrupole resonance (NQR), a variant of NMR, to probe intrinsic spin excitations in the paramagnetic state of Ir honeycomb layers in Cu 2 IrO 3 . A major thrust of our work is that we can probe the Ir honeycomb lattice using NQR in zero magnetic field, owing to the strong nuclear quadrupole interaction at the Cu sites. Moreover, the magnetic hyperfine coupling A hf between Cu nuclear spins and Ir electron spins via the transferred spin polarization into the empty Cu 4s orbital [31] is generally strong. Accordingly, defect spins do not overwhelm the intrinsic 63 Cu NQR and NMR properties of Cu 2 IrO 3 , unlike 1 H NMR data for the kagome lattice in herbertsmithite ZnCu 3 (OH) 6 Cl 2 [32][33][34] and for the honeycomb lattice of H 3 LiIr 2 O 6 [13]. We will demonstrate that Cu 2 IrO 3 remains paramagnetic, and exhibits a large gap ∆ com-parable to the magnitude of Ising interaction, ∆ ∼ 15 meV (= 0.5|J K | ∼ 0.7|J K |) in the primary channel of spin excitations in zero external magnetic field, B ext =0.

II. EXPERIMENTAL METHODS
We synthesized the powder samples of Na 2 IrO 3 , Cu 2 IrO 3 and Cu 1.5 Na 0.5 SnO 3 for NMR measurements via solid state reaction and ion exchange as described in [22,35]. We conducted all NMR and NQR measurements using home-built NMR spectrometers with the standard π/2 -π spin echo pulse sequence. The typical pulse width for the π/2 pulse was 2 to 7 µs, and the pulse separation time τ ∼ 10 µs. We measured 1/T 1 using the inversion recovery method whereby we applied a π pulse, followed by a spin echo sequence for various delay times t. We fitted the inversion recovery of nuclear magnetization M (t) to the exponential function expected for the observed transition [36,37], as discussed in detail in Appendix A. We measured the NQR intensity near and below T f at the peak of the Cu(I 1 ) and Cu(H) sites using a fixed delay time τ = 10 µs. We uni-axially aligned the powder sample of Cu 2 IrO 3 in Stycast 1266 [38] by curing the mixture for 12 hours in a mold made of Teflon in a magnetic field of 9 T.

III. NMR RESULTS
A. 63,65 Cu NQR lineshapes Naturally occurring 63,65 Cu isotopes at Cu(H) and Cu(I) sites have nuclear spin I = 3/2 accompanied by a nuclear quadrupole moment 63,65 Q, and the latter interacts with the electric field gradient (EFG) generated by electrons and ions in the lattice. As shown schematically in Figs. 2B, this nuclear quadrupole interaction lifts the degeneracy of the nuclear spin |±1/2 > and |±3/2 > states even in the absence of an external magnetic field B ext . It is therefore possible to conduct 63,65 Cu NMR based on the NQR techniques in B ext = 0 without perturbing the Ir electron spins with Zeeman effects, a crucial advantage of Cu 2 IrO 3 for NMR investigation into the inherent Kitaev physics.
In Figs. 3A-B, we compare the 63,65 Cu NQR lineshapes observed for Cu 2 IrO 3 and a non-magnetic reference material Cu 1.5 Na 0.5 SnO 3 [35] with Na 0.5 Sn honeycomb layers. The NQR peak frequency 63,65 ν N QR of the 63,65 Cu isotopes is proportional to 63,65 Q, where 65 Q/ 63 Q = 0.927, and hence each 63 Cu NQR peak observed at 63 ν N QR is accompanied by a smaller 65 Cu NQR peak at a lower frequency 65 ν N QR = 0.927 · 63 ν N QR with the intensity ratio set by their natural abundance, 69% and 31%. The two pairs of similar 63,65 Cu NQR peaks observed below 30 MHz for both Cu 2 IrO 3 and Cu 1.5 Na 0.5 SnO 3 indicate that they arise from two distinct types of the inter-layer Cu(I 1 ) and Cu(I 2 ) sites (the O-Cu(I 2 )-O bond has a kink, see Figure 1E). We assign the remaining 63,65 Cu NQR peaks observed at around 52 MHz to the honeycomb Cu(H) sites. We note that the Kitaev candidate materials tend to suffer from stacking faults. In the case of Cu 2 IrO 3 , due to the weak Cu-O-Cu bonds between the layers of Cu 2 IrO 3 , the material is prone to stacking faults in the form of twinning between the adjacent layers [35]. The stacking faults are extended disorders that happen over several unit cells, however, we have to effectively model them in one unit cell when we solve the crystal structure by Rietveld refinements. The result is two inequivalent dumbbell bonds, one straight and one twisted. If both bonds were perfectly straight, there would have been no twisting between adjacent layers and no stacking disorder. Our observation of two distinct types of dumbbell Cu(I 1 ) and Cu(I 2 ) sites supports this.
We summarize the temperature dependence of 63 Cu NQR lineshapes in Figs. 4A-B. We observed no splitting or broadening of the NQR lineshapes due to static hyperfine magnetic fields arising from ordered magnetic moments, a typical NMR signature of antiferromagnetic long-range order (see Appendix C for magnetic NMR line broadening observed for Na 2 IrO 3 below T N ). The NQR frequency 63 ν N QR summarized in Fig. 4C reflects the local lattice and charge environment. 63 ν N QR at 63 Cu(I 1,2 ) decreases smoothly toward higher temperatures due to thermal expansion, without exhibiting evidence for a structural transition or dimerization. We were unable to keep track of the 63 Cu(H) NQR signals above ∼70 K due to extremely weak signal intensity caused by very fast NMR relaxation rates. This is simply because the hyperfine coupling between 63 Cu(H) nuclear spins with six nearest-neighbor Ir spins is stronger than that at the 63 Cu(I 1,2 ) sites, as discussed in detail in Section III E below. Generally, detection of the 63 Cu NQR and NMR spin echo signals becomes difficult when the NMR spinlattice relaxation rate 1/T 1 reaches ∼ 10 4 sec −1 [32].
Although there is no evidence for magnetic long range order, we found that the integrated intensity of the 63 Cu NQR signals begins to decrease gradually below T f ∼ 10 K. This indicates that NMR relaxation rates become very fast in some domains. In Fig. 5A, we sum- as measured by zero field µSR techniques using the MuSR (triangles) and EMU (downward triangles) spectrometers at the Rutherford Appleton Laboratory [25]. (B) A precursor of the spin freezing observed above T f in the spin-lattice relaxation rate 1/T1. Intrinsic 1/T1, obtained as the slower component 1/T 1,slow from the two component fit with eq.(A3) below ∼ 60 K (i.e. 1/T1 = 1/T 1,slow below ∼ 60 K), continues to slow down at low temperatures. In contrast, the sample averaged behavior represented by 1/T1,str, estimated from the stretched fit of M (t) with eq.(A2), levels off below ∼30 K. See Appendix A for details of the fit.
marize the temperature dependence of the lost fraction f N QR of the NQR spin echo intensity, measured at a constant pulse separation time τ = 10 µs, down to 1.5 K. f N QR shows qualitatively the same behavior as the volume fraction f µSR estimated by µSR experiments [25], in which muons exhibit fast depolarization due to slow spin fluctuations in their vicinity [23,25]. Note that f µSR saturates at 50% below 1 K. (f N QR overestimates the volume fraction with slowed spin fluctuations, because the distribution of the spin-spin relaxation time T 2 prevented us from accurately taking into account the contributions of nuclear spins with fast T 2 .) Combined with the fact that we found no evidence for divergent NMR relaxation rates for the observable NQR signals, we conclude that the loss of the 63 Cu NQR signals are primarily due to slowing of defect spin fluctuations below T f . Since the 63 Cu NQR lineshape for observable signals does not broaden even below T f , the defect spins may be spatially confined in some domains, perhaps in the vicinity of stacking faults [25], rather than uniformly distributed throughout the entire Ir 4+ honeycomb planes.
B. Low energy Ir spin excitations in Bext = 0 We now turn attention to our central results on lowfrequency Ir spin dynamics probed by 63 Cu nuclear spinlattice relaxation rate, 1/T 1 . Quite generally, 1/T 1 probes the low frequency component at the NQR frequency ν N QR of the fluctuating hyperfine magnetic field h perp = A hf S arising from the Ir spin S, where γ n /2π = 11.285 MHz/T is 63 Cu nuclear gyromagnetic ratio, and Perp indicates the component orthogonal to the quantization axis of the observed nuclear spins. The quantization axis of 63 Cu nuclear spins is along the caxis for our NQR measurements, and hence our NQR 1/T 1 results probe the slow component at 63 ν N QR of Ir spin fluctuations within the ab-plane.
One can also write 1/T 1 ∝ Σ q |A hf (q)| 2 S(q, E n ), where A hf (q) is the wave-vector q-dependent hyperfine form factor [39]. S(q, E n ) is the dynamic spin structure factor of Ir spins at the very small excitation energy of E n = h · 63 ν N QR ∼ 0.2 µeV, corresponding to the photon energy at the NQR frequency. 1/T 1 therefore probes the low energy sector of the q-integrated spin excitation spectrum S(q, E n ) at energy E n , weighted by |A hf (q)| 2 .
Summarized in Figure 6A is 1/T 1 at the 63 Cu(H), 63 Cu(I 1 ), and 63 Cu(I 2 ) sites measured at the peak frequency of the NQR lineshapes. The hyperfine coupling A hf of Ir spins with the 63 Cu(H) nuclear spins is stronger than with Cu(I 1,2 ) sites due to their spatial proximity, and hence the magnitude of 1/T 1 at the 63 Cu(H) sites is larger by a factor of ∼ 3. As noted above, this is the underlying reason why the 63 Cu(H) NQR signal detection becomes difficult above ∼70 K. Otherwise, the temperature dependence of 1/T 1 is qualitatively the same between the 63 Cu(H), 63 Cu(I 1 ), and 63 Cu(I 2 ) sites.
1/T 1 grows slightly below 300 K, a typical signature of the gradual development of short-range spin-spin correlations, in the present case between Ir spins. It is followed by an onset of a dramatic suppression of 1/T 1 below ∼150 K. A semi-logarithmic plot in Figure 7A indicates that 1/T 1 follows an activation behavior over two decades, from ∼100 K to ∼20 K, due to a complete suppression of the low frequency component of the fluctuating hyperfine magnetic fields h perp . By fitting the result to an activa- We note that 1/T 1 develops a distribution below ∼60 K, and hence the 1/T 1 results plotted below ∼ 60 K with filled symbols in Figs. 5B, 6A and 7 represent the intrinsic slow component (i.e. 1/T 1,slow ) estimated from the two component fit with eq. (A3), as explained in detail in Appendix A. We emphasize, however, that our key finding of the activation behavior in Figs. 7A does not depend on the details of how we estimate 1/T 1 from the nuclear spin recovery curve M (t). The proof comes from the distribution function P (1/T 1 ) of 1/T 1 directly calcu- lated from the inverse Laplace transform of the nuclear magnetization M (t) [40,41], as explained in the next section C. Also plotted in Fig. 5B and 7 using open bullets are 1/T 1,str obtained from the phenomenological stretched fit of M (t) with eq. (A2) [42,43]. As explained in section III-C, 1/T 1,str is a good measure of the center of gravity of the distributed 1/T 1 . As the influence of inhomogeneous spin fluctuations from defects grow near T f , the aforementioned activation behavior of the intrinsic slow component of 1/T 1 is terminated by the precursor of spin freezing, which also results in a shoulder of 1/T 1,str below ∼30 K as shown in Fig. 5B. This is followed by a gradual loss of 63,65 Cu NQR signal intensity below T f ∼ 10 K in Fig. 5A. The fact that a half of the NQR signal remain observable without exhibiting any magnetic anomalies (such as the line broadening, divergence of 1/T 1 , or signal loss) indicates that a majority of Ir spins themselves are not slowing down toward eventual long range order or freezing. Instead, it is the inhomogeneously distributed defects (perhaps accompanied by minority Ir spins adjacent to them) that are freezing. This is consistent with µSR observation that Ir spins remain dynamic even below T f [23,25], and also with a very wide distribution of 1/T 1 near T f as determined by inverse Laplace transform analysis in section III-C. Whether these defects originate from randomly distributed Cu 2+ ions and/or clusters around the stacking faults encompassing several unit cells remains to be seen.
A more precise, model-independent way to investigate the spin dynamics with 1/T 1 under the presence of distributions in the relaxation mechanisms is to take the inverse Laplace transform of M (t), and directly deduce the distribution function P (1/T 1 ) of 1/T 1 [40,41], where we normalize the overall probability to 1, i P (1/T 1,i ) = 1. We numerically inverted M (t) [41] utilizing Tikhonov regularization (i.e. a smoothing factor) [44,45], and deduced the distribution function P (1/T 1 ) as summarized in Fig. 8 at representative temperatures.
We do not presume any phenomenological functional forms for M (t), such as the existence of two separate components in eq.(A3). Nonetheless, our P (1/T 1 ) results indicate that a new small side peak begins to emerge with fast values of 1/T 1 below ∼60 K. The emergence of the two peak structure in P (1/T 1 ) indicates that the intrinsic, slow T 1 process begins to be short circuited by the faster component(s) in the vicinity of defect spins as the intrinsic T 1 slows down. Since the influence of the defectinduced fast relaxation depends on the distance with observed nuclear spins, defect induced fast component in 1/T 1 generally has a broad distribution.
To show the key features of P (1/T 1 ) visually, we plotted P (1/T 1 ) in a color contour map in Fig. 7B. We can identify the existence of the well-defined, slower component of 1/T 1 in the ridge structure (colored in red or green) that continues from ∼60 K down to ∼15 K (i.e.  [41]. For clarity, the origin is shifted vertically except for 4.2 K. P (1/T1) is centered around a single value of 1/T1 above 60 K; the peak location at 77 K and 200 K agrees very well with 1/T1 estimated from the single component fit with eq. (A1), as marked by thick arrows. This one peak distribution of P (1/T1) breaks down at lower temperatures. The filled and open arrows at and below 45 K mark 1/T 1,slow and 1/T 1,f ast obtained from the two component fit with eq.(A3), respectively. The dashed lines mark 1/T1,str estimated from eq.(A2), which keeps track of the center of gravity of P (1/T1) marked by black filled bullets.
1/T = 0.015 to 0.07). On the other hand, the defect induced fast contribution to 1/T 1 develops a wide distribution centered around ∼ 150 s −1 , and remains roughly constant, see the light blue horizontal section starting from 1/T ∼ 0.04 to 0.1 in Fig. 7B. Interestingly, recent theoretical calculations of 1/T 1 for Kitaev QSL predict that the presence of bond-disorder leads to a constant 1/T 1 with a broad distribution [30]. We emphasize that our observation of the activation behavior for the intrinsic slow component and the roughly constant behavior induced by defects, both nicely captured in Fig. 7B, does not depend on fitting procedures of M (t), such as the stretched nature of the relaxation function in eq.(A2) or the presence of two components in eq.(A3).
Also notice that 1/T 1,slow and 1/T 1,f ast estimated from the two component fit (marked by thin filled and thin open arrows, respectively, in Fig. 8) agree fairly well with the two peaks observed for P (1/T 1 ) in a wide temperature range except near and below T f , where the two component fit itself becomes invalid. This justifies the two component fit with eq.(A3) used above T f to generate intrinsic 1/T 1 as 1/T 1,slow in Figs. 5B, 6A and 7. Interestingly, 1/T 1,str estimated using eq.(A2) (marked by vertical dashed lines in Fig. 8) agrees well with the center of gravity of the distribution P (1/T 1 ) (marked by black filled bullets). This means that 1/T 1,str probes only the spatially averaged behavior of the entire sample, and does not always represent the intrinsic behavior of 1/T 1 .

D. Comparison with low energy Ir spin excitations
of Na2IrO3, α-RuCl3, and H3LiIr2O6 It is illustrative to compare the 1/T 1 results of Cu 2 IrO 3 with those observed for other Kitaev candidate materials. In Fig. 6B, we summarize 1/T N a 1 measured for the inter-planar 23 Na(I) sites of Na 2 IrO 3 at the central peak of the powder 23 Na NMR lineshape. We refer readers to Appendix C for the details of 23 Na NMR measurements. Simutis et al. also reported analogous 1/T N a 1 data in ambient and applied pressure [46]. Their results, reported in a limited temperature range up to 25 K, are in good agreement with ours. The overall magnitude of 1/T N a 1 is nearly two orders of magnitude slower than at the Cu(I 1,2 ) sites in Fig. 6A, because the hyperfine coupling A hf is much smaller at the 23 Na sites, as evidenced by the 23 Na NMR Knight shift results (see Fig. 11B in Section III E).
Below 300 K down to ∼25 K, 1/T N a 1 gradually increases as the short-range Ir-Ir spin correlations grow, without exhibiting a downturn observed for Cu 2 IrO 3 . The critical slowing down of Ir spin fluctuations below ∼25 K toward the three dimensional zig-zag antiferromagnetic long range order results in a sharp divergent behavior of 1/T N a 1 at T N = 16.5 ± 0.5 K. This is consistent with the fact that the hyperfine form factor does not cancel at the 23 Na(I) sites for the zig-zag order pattern (schematically shown in Fig. 2B). Analogous divergent behavior is commonly observed for conventional antiferromagnets, such as CuO [47].
is quickly suppressed, accompanied by dramatic broadening of the 23 Na powder NMR lineshape due to the static hyperfine magnetic fields from zig-zag ordered Ir spins (see Fig. 17 in Appendix C as well as [46]). In general, 1/T 1 in the antiferromagnetically ordered state below T N is dominated by two-or three-magnon Raman process, resulting in suppression of 1/T 1 obeying a power law; it is generally followed by an activation behavior 1/T 1 ∼ T 2 ·exp(−∆ magnon /k B T ) due to the anisotropy gap ∆ magnon for magnons at lower temperatures below ∆ magnon /k B [48,49]. While our 1/T N a 1 results were measured for broad and superposed 23 Na NMR lines and do not allow us to conduct detailed analysis below T N , the observed behavior of 1/T N a 1 below T N is similar to the conventional behavior expected for three-dimensional or quasi-two dimensional antiferromagnets, such as CuO and YBa 2 Cu 3 O 6 [49].
Our conventional findings on Na 2 IrO 3 are no surprise, in view of the fact that earlier inelastic neutron scattering measurements established that the low energy sector of the spin excitation spectrum in Na 2 IrO 3 is dominated by magnons that are not native to pure Kitaev Hamiltonian in eq.(1) [18]. In this context, it is also interesting to compare our results with those previously reported for α-RuCl 3 [50][51][52], because magnons also dominate its spin excitation spectrum in the low field regime [19][20][21]. Not surprisingly, 1/T 1 observed at the 35 Cl site of α-RuCl 3 in the low field regime B ext = 2.35 T with antiferromagnetic ground state [52] is similar to the conventional behavior observed here for Na 2 IrO 3 . Our NQR results for Cu 2 IrO 3 in Figs. 6A and 7 are fundamentally different.
In the case of α-RuCl 3 , application of a modest magnetic field above ∼ 7 T drives the antiferromagnetic ground state to a field-induced paramagnetic state with a gap [50][51][52]. This is often attributed to a spin liquid state. The report of fractionalized thermal Hall effect [53] in the intermediate field regime of α-RuCl 3 as well as peculiar magnetic field effect on magnon dispersion [21] certainly raises the hope that chiral spin-liquid with Majorana edge modes may be induced by a magnetic field, but even the classical spin wave theory can semi-quantitatively account for the former [54]. We also caution that application of an external magnetic field could fundamentally alter the nature of spin excitations of quantum spin systems. For example, B ext exceeding 5.3 T applied to the transverse field Ising chains realized in CoNb 2 O 6 suppresses a magnetic long-range order and induces a spin excitation gap that manifests itself in the 59 Co NMR properties [55]. The nature of the fieldinduced gap and its potential relation to Kitaev physics has been at the center of recent intense debate [21,[50][51][52][53], and deserves further careful examinations.
Finally, we compare Cu 2 IrO 3 and H 3 LiIr 2 O 6 [13]; both of these Kitaev candidate materials remain paramagnetic without undergoing magnetic long-range order. Interestingly, 1/T 1 observed at the intra-layer 7 Li sites of H 3 LiIr 2 O 6 shows qualitatively the same behavior as the intra-layer 63 Cu(H) as well as inter-layer 63 Cu(I 1 ) sites from 300 K down to ∼50 K, where the intrinsic behavior of 7 Li 1/T 1 appears to be taken over by large defect induced contributions. On the other hand, 1/T 1 observed at the inter-layer 1 H sites in H 3 LiIr 2 O 6 , which corresponds to the 63 Cu(I) sites in the present case, shows no hint of a gapped behavior at low temperatures. This apparent discrepancy from the 7 Li or our 63 Cu results, however, does not necessarily mean that H 3 LiIr 2 O 6 does not have analogous activation behavior of 1/T 1 . It is well-known that 1 H NMR 1/T 1 results in herbertsmithite kagome antiferromagnet ZnCu 3 (OH) 6 Cl 2 [32] are dominated by defect spins, especially in low magnetic fields, and do not resemble the intrinsic 1/T 1 results observed at 63 Cu [32] or 17 O [34] sites. This is simply because Fermi's contact hyperfine coupling of 1 H nuclear spin through the 1s orbital is generally very weak, and hence 1/T 1 measured at 1 H sites is comparatively more sensitive to defect induced spin fluctuations. We note that the 1 H as well as 7 Li NMR linewidth in H 3 LiIr 2 O 6 is greater than the small NMR frequency shift (i.e. Knight shift is very small). This suggests that a large fraction of nuclear spins have stronger hyperfine couplings with defect spins than with nearby Ir 4+ spins with intrinsic Kitaev properties, which explains why 7 Li and 1 H NMR 1/T 1 results are overwhelmed by defect contributions. In the case of Cu 2 IrO 3 , this does not happen, because the transferred hyperfine couplings of 63 Cu nuclear spins through the 4s orbital [31] with surrounding Ir sites are larger by an order of magnitude, as demonstrated in the next section.

E. NMR measurements in applied magnetic field
We also conducted high-field NMR measurements using a uni-axially aligned powder sample of Cu 2 IrO 3 cured in a glue (stycast 1266) in the presence of B ext = 9 T. This is a standard trick used in powder NMR experiments, and creates a pseudo single crystal by taking advantage of the anisotropy of magnetization in high magnetic fields [38]. We summarize representative 63 Cu NMR lineshapes observed for the nuclear spin I z = +1/2 to −1/2 central transition of the aligned powder sample in B ext = 9 T applied along the aligned c-and ab-axis in Fig. 9B and C, respectively. See Fig. 15 in Appendix B for the field sweep NMR spectra covering the satellite transitions between the nuclear spin I z = ±1/2 to ±3/2 states. From the comparison with the NMR lineshapes for unaligned powder in Fig. 9A, we estimate the fraction of the uni-axially aligned powder as approximately 50%; the rest of the powder remains randomly oriented.
We could assign the 63,65 Cu NMR peaks and shoulders arising from Cu(I 1,2 ) sites in aligned and unaligned portions of the sample, as marked by arrows, + and x symbols in Fig. 9A-C. However, identification of the NMR signals from less abundant 63 Cu(H) sites proved to be a major challenge, because the 63,65 Cu(H) central transition is shifted by the second and higher order effects caused by the extremely large nuclear quadrupole frequency ( 63 ν N QR ∼ 52 MHz), and superposed by much stronger central and satellite transitions of the 63,65 Cu(I 1,2 ) sites. Moreover, since 63 ν N QR at the 63 Cu(H) sites happens to be almost exactly twice larger than 63 ν N QR at the 63 Cu(I 1,2 ) sites, the small I z = ±3/2 to ±1/2 satellite peaks of the 63 Cu(H) sites are accidentally superposed by the satellite structures arising from the 63 Cu(I 1,2 ) sites. To make matters worse, defect spins polarized by B ext broaden all the NMR lines below ∼50 K, where 63 Cu(H) NQR signal was readily observable. Nonetheless, we were able to resolve the 63 Cu(H) NMR central transition peak below ∼120 K in the B ext || c geometry, as marked by blue arrows in Fig. 9B.
The resonance frequency of the 63 Cu(I 1,2 ) and 63 Cu(H) peak for B ext || c in Fig. 9B is shifted from the bare Zeeman frequency 63 ν o marked by the vertical dashed line due to the paramagnetic NMR Knight shift, 63 K = N nn A hf χ spin /N A µ B , where N nn is the number of the nearest neighbor Ir sites that induce transferred hyperfine magnetic fields at Cu nuclear spins, N A is Avogadro's number, µ B is Bohr magneton; χ spin is the spin contribution to the bulk susceptibility χ = χ spin + χ dia + χ vv , where the diamagnetic term χ dia ∼ −0.09 × 10 −3 emu/mol from the standard table, and χ vv ∼ 0.16 × 10 −3 emu/mol [16]. 63 Cu(I 1,2 ) sites are bonded with two O sites above and below (see Fig. 1E), each of which bonds with two Ir sites, resulting in N nn = 4. On the other hand, N nn = 6 for 63 Cu(H), as may be seen in Fig. 1D.
In the case of the 63 Cu(I 1,2 ) peak in B ext || ab, the apparent Knight shift 63 K apparent is also affected by the second order contribution of the nuclear quadrupole interaction ∆ν (2) Q , which is proportional to 1/(γ n B ext ) 2 . We therefore measured 63 K apparent at different magnetic field values, linearly extrapolated the results to the high field limit of 1/(γ n B ext ) 2 = 0 [38], and estimated 63 K as shown in Fig. 10. On the other hand, within experimental uncertainties, we found no dependence of 63 K apparent on 1/(γ n B ext ) 2 for both 63 Cu(I 1,2 ) and 63 Cu(H) sites in the B ext || c geometry. This implies that the nuclear quadrupole effects ∆ν (2) Q play no significant roles in the NMR frequency shift for B ext || c, because the main principal axis of the EFG tensor is along the c-axis and the asymmetry parameter η of the EFG tensor is negligible (η ∼ 0). This is consistent with the local symmetry of these sites.
We summarize the temperature dependence of 63 K in Fig. 11A in comparison to χ spin measured in 5 T (χ spin does not depend on magnetic field above ∼ 30 K [22]). 63 K increases with decreasing temperature, and saturates below ∼30 K. Interestingly, 63 K observed for Cu 2 IrO 3 shows qualitatively the same site and temperature dependences as the 23 Na NMR shift 23 K observed for Na 2 IrO 3 in the paramagnetic state above T N ∼ 16.5 K, as shown in Fig. 11B. Noting that 63 K probes local spin susceptibility rather than the bulk averaged result, our 63 K results indicate that the upturn of the bulk χ data observed below ∼30 K does not reflect the intrinsic behavior of the Ir honeycomb planes, but should be attributed to defect spins. It is the enhanced local spin susceptibility of the defect spins that inhomogeneously broadens our high field NMR peaks at low temperatures.
We plotted 63 K vs. χ spin choosing temperature as the implicit parameter in Fig. 12, and estimated the hyperfine coupling from the slope of the linear fit as A hf = (N A µ B /N nn ) · d( 63 K)/dχ spin = 2.3 and 7.6 kOe/µ B for 63 Cu(I 1,2 ) in the B ext || c and || ab geometry, respectively, and A hf = 4.4 kOe/µ B for 63 Cu(H) in B ext || c. We were unable to determine A hf at the 63 Cu(H) sites along the ab direction, because the small central peak is buried underneath the 65 Cu(I 1,2 ) signals due to extremely large ∆ν (2) Q as explained earlier.

IV. DISCUSSIONS
The activation behavior of 1/T 1 is generically observed when the spin excitation spectrum represented by S(q, E) has a gap in the low energy sector. For example, the collective spin-singlet state realized in the twoleg spin ladders SrCu 2 O 3 [56] and La 6 Ca 8 Cu 24 O 41 [57] has a spin excitation gap as large as ∆/k B ∼ 450 K (∼ J H /2 along the rung formed by a pair of Cu ions), and exhibits an activation behavior of 1/T 1 at both 63 Cu [56,57] and 17 O [57] sites below ∼500 K. Our 1/T 1 results in Fig. 7 is therefore consistent with the presence of a gap ∆/k B = 175±30 K (∆ ∼ 15 meV) in the primary mode of the intrinsic spin excitation spectrum in Cu 2 IrO 3 . Note that this gap is comparable with the Ising interaction energy scale estimated above, |J K | = 17 ∼ 30 meV. In the pure Kitaev model, earlier theoretical analyses showed that the dispersion of the spin 1/2 excitations of fractionalized Majorana fermions may indeed have a gap as large as ∆ M ajorana ∼ 1.2J K [6,8]. It is therefore tempt- ing to associate our finding with the gapped excitations expected for Majorana fermions. Broad high-energy spin excitations centered around q = 0 with a gap ∆ ∼ 1.2|J K | was previously reported for α-RuCl 3 by inelastic neutron scattering, but the low energy sector is filled by additional spin-wave like excitations [19,20]. In the present case, our observation of the faster component in 1/T 1 in the recovery curve M (t) indicates the presence of additional, spatially inhomogeneous relaxation mechanism(s). Theoretically, excitations of a pair of gauge fluxes with total spin 1 could short-circuit the slower NMR relaxation process in the pure Kitaev model, and such gauge flux excitations might morph into the spin-waves (also with net spin 1) under the perturbation caused by the Heisenberg exchange term. It is not clear what roles the gauge fluxes and/or spin waves may be playing in our results. On the other hand, the growing distribution of P (1/T 1 ) toward T f seems to suggest that the faster component of 1/T 1 originates primarily from extrinsic defect spins.
It is important to distinguish our observation of the gapped behavior in zero magnetic field from earlier observations of a small field-induced gap in α-RuCl 3 [50][51][52]. To test the potential influence of the applied magnetic field on the gapped behavior observed for Cu 2 IrO 3 , we measured 1/T 1 at the lower satellite transition of the 63 Cu(I 1,2 ) sites in B ext = 9 T applied along the aligned ab-plane. In this field geometry, 1/T 1 probes the fluctuations of h perp along both the ab-plane and c-axis. We summarize the NMR results with × symbols in Figs. 6A and 7A. Comparison with the NQR results indicates that the gapped behavior of 1/T 1 does not depend significantly on the magnetic field, at least up to 9 T.
A potential caveat in our interpretation of the gapped behavior of S(q, E) is that, in principle, 1/T 1 can be suppressed at low temperatures if the form factor satisfies A hf (q AF ) = 0 at a staggered vector q AF , where S(q AF , E n ) grows due to short-range spin-spin correlations [39]. In this alternate scenario, Ir spins on the entire honeycomb planes continuously grow strong short-range order down to T f without developing long-range order, resulting in diminishing 1/T 1 due to geometrical cancellation of A hf (q AF ) by the high symmetry at Cu sites.
To illustrate the underlying mechanism in this alternate scenario, we present the possible short-range order patterns of Ir spins [18] in Fig. 1B-C. Unlike Na 2 IrO 3 with strong zig-zag short range order, 1/T 1 in Cu 2 IrO 3 does not grow below ∼ 150 K. Accordingly, we can rule out strong zigzag short range order for Cu 2 IrO 3 . However, if Ir spins develop extremely strong Néel-type short- range spin-spin correlations shown in Fig. 1C, the transferred hyperfine field h perp from 6 nearest-neighbor Ir sites of 63 Cu(H) would cancel out, because the sign of h perp alternates between six nearest-neighbors. Although the location of the 63 Cu(I 1,2 ) sites is somewhat shifted from the mid-point between two adjacent Ir sites, h perp should also nearly cancel out for 63 Cu(I 1,2 ) sites. Thus, in principle, Néel-type short-range spin-spin correlations can suppress 1/T 1 in Cu 2 IrO 3 . However, we recall that such perfect cancellation of h perp expected for 1/T 1 ∼ 0 observed here generally requires extremely long spin-spin correlation length ξ. For example, in the case of 17 O NMR measurements in the square-lattice Heisenberg antiferromagnet Sr 2 CuO 2 Cl 2 with no spin excitation gap above T N [58], ξ reaches tens of lattice spacings [59] when 1/T 1 is completely suppressed due to the high symmetry at the 17 O sites. In view of the frustrated geometry of Ir spins in the present case, validity of such a scenario starting from as high as ∼ 150 K (comparable to ∼ |J K |/k B , and much higher than |J H |/k B ) seems somewhat remote. In addition, the intrinsic uniform spin susceptibility deduced from 63 K in Figure 11A does not exhibit any hint of a downturn below ∼150 K, which is expected for such strong quasi-two-dimensional short-range Néel order. It is also worth recalling that 17 O as well as 63 Cu NMR 1/T 1 measurements of two-leg spin-ladders in La 6 Ca 8 Cu 24 O 41 successfully singled out the gapped spin excitation mode at ∼40 meV, even though A hf (q) cancels at the 17 O sites [57]. Next, we wish to address the nature of defect spins and their influence on the NMR properties. As noted in section I, recent XANES measurements for Cu 2 IrO 3 showed that up to ∼1/3 of 63 Cu(H) sites are occupied by Cu 2+ ions rather than Cu + ions. On the other hand, the linear fit of the 63 K vs. χ plot in Fig. 12 extrapolates to the origin for both 63 Cu(H) and 63 Cu(I 1,2 ) sites. This means that the van Vleck contribution to the Knight shift is vanishingly small, as expected for Cu + ions with filled 3d orbitals. If the observed 63 Cu(H) NMR signals arise from Cu 2+ ions, the intercept with the vertical axis of the extrapolated linear fit would be as large as 0.3% to 1.5% due to the temperature independent van Vleck contribution for Cu 2+ ions [31].
The lack of Cu 2+ contributions in our 63 Cu(H) NMR data may be easily understood, if we recall that NMR signals at paramagnetic cation sites are generally not observable due to extremely fast NMR relaxation rates. In fact, 63 Cu NQR and NMR signals in the paramagnetic state of the aforementioned ZnCu 3 (OH) 6 Cl 2 [32], CuO [47], SrCu 2 O 3 [56], La 6 Ca 8 Cu 24 O 41 [57], and Sr 2 CuO 2 Cl 2 [58] etc. are observable only because exchange narrowing effects induced by strong antiferrmagnetic exchange interaction J H /k B of the order of 200 to 1700 K between cation sites suppresses the NMR relaxation rates, which certainly would not be the case for the Cu 2+ defect spins in Cu 2 IrO 3 . We therefore conclude that the primary contributions of Cu 2+ defect spins to our NMR results are through the line broadening effects for the observable 63 Cu(H) and 63 Cu(I 1,2 ) NMR signals below ∼50 K and the enhancement of 1/T 1 in their vicinity. The latter leads to a large distribution of 1/T 1 , as seen in P (1/T 1 ) in Figs. 7B and 8.

V. SUMMARY AND CONCLUSIONS
We have reported comprehensive 63 Cu NQR and NMR measurements in Cu 2 IrO 3 , a proximate Kitaev QSL material that does not undergo magnetic long-range order. We showed that the intrinsic uniform spin susceptibility χ spin of Cu 2 IrO 3 shows nearly identical behavior as in the paramagnetic state of antiferromagnetic Na 2 IrO 3 above T N ∼ 17 K. On the other hand, the low energy sector of Ir spin excitations reflected on the temperature dependence of 1/T 1 is very different in Cu 2 IrO 3 . 1/T 1 shows activation behavior and is qualitatively different from the conventional behavior exhibited by antiferromagnetic Na 2 IrO 3 . Our finding strongly suggests that the low energy sector of the spin excitation spectrum in Cu 2 IrO 3 is dominated by a mode with a gap comparable to the magnitude of Ising interaction |J K |, which may be related to gapped excitations expected for fractionalized Majonara fermions perturbed by Heisenberg exchange term. In such an idealized scenario of the Kitaev lattice, it is not clear what roles the excitations of a pair of gauge fluxes might play. Analogous dilemma was previously encountered in the analysis of the spin excitations of α-RuCl 3 , because the low energy sector is dominated by magnons rather than gauge fluxes [20]. From the measurements of 1/T 1 alone, we cannot entirely rule out an alternate scenario that Ir spin-spin correlation develops strong Néel type short range order with diver-gently long spin-spin correlation length. However, that would normally require suppression of χ spin in quasi twodimensional Heisenberg antiferromagnets, and we did not find such a signature in our 63 K results.
The results of the nuclear magnetization M (t) for the 1/T 1 measurements at 63 Cu(I 1 ), 63 Cu(I 2 ) and 63 Cu(H) sites are qualitatively similar, and we applied the same fitting procedures. In what follows, we will focus on the results for the 63 Cu(I 1 ) sites that have the highest signal to noise ratio. In Figure 13A-B, we summarize T 1 recovery curves M (t) observed for the 63 Cu(I 1 ) site at selected temperatures. For the 63 Cu NQR measurements between the I z = ±3/2 and ±1/2 states in zero magnetic field, magnetic transition would result in the exponential form of the recovery function [36,37], where 1/T 1 , M o and A are the fitting parameters. We found that this single exponential form fits the recovery data nearly perfectly above ∼60 K, as shown by dashed straight lines. Below ∼60 K, 1/T 1 begins to exhibit a distribution and hence the fit with eq.(A1) becomes progressively poorer. As a remedy, one can introduce a phenomenological stretched exponent β introduced for ranfomly diluted antiferromagnets [42,43], If there is no distribution in the relaxation process, the stretched fit would give 1/T 1,str = 1/T 1 and β = 1. As summarized in Fig. 13C, the stretched exponent begins to deviate from β = 1 below ∼ 60 K as the intrinsic relaxation process slows down and becomes susceptible to the faster defect contributions. 1/T 1,str results are approximate, but tend to elucidate the sample averaged behavior, as proven by the P (1/T 1 ) results in Fig. 8 based on inverse Laplace transform of M (t). That is, the center of gravity of the distribution function P (1/T 1 ) agrees fairly well with 1/T 1,str . Our goal, however, is not merely to find the sample averaged behavior. A better approach to capture the intrinsic slow component 1/T 1,slow is the more traditional two-component fit of M (t), Here, M o , A 1 , the intrinsic slow component 1/T 1,slow , A 2 , the defect-induced fast component 1/T 1,f ast are the fitting parameters. We also introduced the phenomenological stretched exponent β in the second term to improve the fit (the estimated values of β from eq.(A3) are similar to those from eq.(A2)). It is worth noting that analogous two component fit was previously applied successfully to diluted antiferromagnets Mn 1−x Zn x F 2 [43]. range below ∼60 K down to T f , as shown in Fig. 8. Fig. 14A  In Fig. 15, we compare the 63,65 Cu NMR lineshapes measured at 90.307 MHz for uniaxially aligned and unaligned powder samples at 10 K while sweeping the external magnetic field. For the purpose of determining the NMR Knight shift accurately, we measured only the I z = −1/2 to +1/2 central peak region in a fixed magnetic field of B ext = 9 T while sweeping the frequency, as shown in Fig. 9A-C.
The lineshape for the unaligned powder sample in Fig.  15C exhibits a typical powder pattern with a double horn structure arising from the I z = −1/2 to +1/2 central transition of the 63 Cu and 65 Cu isotopes, as identified by filled green arrows. The horn at the lower field side originates from nuclear spins in the grains whose primary principal axis of the electric field gradient (EFG) ten-  sor (i.e. the crystal c-axis) forms a right angle θ ∼ 90 • with the external magnetic field, whereas the horn at the higher field side has θ ∼ 41 • . (Their relative locations are reversed in the frequency sweep lineshapes in Fig. 9.) These structures are smeared by magnetic line broadening at 10 K caused by defect spins.
When we apply B ext along the ab-plane of the aligned powder sample, the θ ∼ 90 • peak is enhanced, as seen in both Fig. 15B and Fig. 9C. On the other hand, when we apply B ext along the aligned c-axis, a new peak emerges in the middle of the θ ∼ 90 • and θ ∼ 41 • peaks as seen in Fig. 15A and Fig. 9B; this is because the effects of the nuclear quadrupole interaction on the NMR frequency shift, ∆ν (2) Q , vanishes for both Cu(I) and Cu(H) along θ ∼ 0 • due to their negligibly small asymmetry of the EFG tensor. Since the NMR properties of the Cu(I 1 ) and Cu(I 2 ) sites are very similar, we were unable to clearly resolve their broad NMR lines.
At 10 K, the Cu(H) central peaks are hidden by stronger signals arising from the Cu(I 1 ) and Cu(I 2 ) sites, because the quadrupole effects are very strong at Cu(H) sites and the NMR line is broadened by defect spins. Notice, however, that the 63,65 Cu NMR signal in Fig. 15 extends below the lower bound expected for the I z = ±1/2 to ±3/2 satellite transitions of Cu(I) sites, and diminishes instead at the edge near B ext ∼ 3.4 T as expected from 63 ν Q ∼ 52 MHz at 63 Cu(H) sites. We tested the potential effect of external magnetic field B ext on Ir spin dynamics by measuring 1/T 1 for the aligned powder sample in B ext = 9 T, using the θ ∼ 90 • I z = ±1/2 to ±3/2 lower satellite transition for the 63 Cu(I) sites. The standard recovery function in this case takes the form of [36,37] as shown in Fig. 16. The high field NMR 1/T 1 results presented with × symbols in Fig. 6A and 7A, are nothing but 1/T 1,slow from the first term in eq.(B2), and the temperature dependence is nearly identical with the NQR results. 1/T 1,f ast and β estimated from eq.(B2) are also similar to the NQR results.  I z = +1/2 to -1/2 transition. We assign the main peak centered at ∼ 101.37 MHz to the interlayer 23 Na(I) sites. The smaller side peak at higher frequency ∼ 101. 46 MHz with larger Knight shift should be attributed to the less abundant sodium site within the honeycomb plane, referred to as 23 Na(H). We were also able to resolve shoulders from I z = ±3/2 to ±1/2 satellite transitions at ∼ 102.1 MHz and ∼ 100.8 MHz. An additional small peak/shoulder at ∼ 101.3 MHz, marked with *, is at the non-shifted frequency, and should be attributed to a small amount of non-magnetic impurity phase. The lineshape broadens dramatically below T N ∼ 16.5 K due to the emergence of static hyperfine magnetic fields from the ordered Ir moments. We also measured 1/T 1 at the relatively isolated 23 Na(I) central peak. We summarize the recovery curves M (t) in Fig. 18, and the temperature dependence of 1/T N a 1 in Fig. 6B. Despite the superposition of NMR signals arising from the satellite transitions of the 23 Na(I) sites as well as 23 Na(H) sites, the fit of the recovery curve M (t) with the standard form for the central transition [36,37], was good with fixed β = 1 down to ∼30 K. To improve the fit below T N , we introduced the stretched exponent β as summarized in the inset of Fig. 18. Since the value of β remains close to ∼ 1 above T N , the qualitatively important aspects of 1/T 1 results in Fig. 6B do not depend significantly on the details of the fit, with or without β.