Thank you for visiting nature.com. You are using a browser version with
limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off
compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site
without styles and JavaScript.

Subjects

Abstract

Quantum spin liquid is a disordered but highly entangled magnetic state with fractional spin excitations1. The ground state of an exactly solved Kitaev honeycomb model is perhaps its clearest example2. Under a magnetic field, a spin flip in this model fractionalizes into two types of anyon, a quasiparticle with more complex exchange statistics than standard fermions or bosons: a pair of gauge fluxes and a Majorana fermion2,3. Here, we demonstrate this kind of fractionalization in the Kitaev paramagnetic state of the honeycomb magnet α-RuCl3. The spin excitation gap determined by nuclear magnetic resonance consists of the predicted Majorana fermion contribution following the cube of the applied magnetic field2,4,5, and a finite zero-field contribution matching the predicted size of the gauge flux gap2,6. The observed fractionalization into gapped anyons survives in a broad range of temperatures and magnetic fields, which establishes α-RuCl3 as a unique platform for future investigations of anyons.

Main

In many-body systems dominated by strong fluctuations, an excitation with a well-defined quantum number can break up into exotic quasiparticles with fractional quantum numbers. Well-known examples include fractionally charged quasiparticles in fractional quantum Hall effect7, spin–charge separation in one-dimensional conductors8 and magnetic monopoles in spin ice9. A major hunting ground for novel fractional quasiparticles is disordered magnetic states of interacting spin-1/2 systems governed by strong quantum fluctuations, called quantum spin liquids (QSLs). Most of their models predict that a spin-flip excitation fractionalizes into a pair of spinons, each carrying spin 1/2 (ref. 1). Even more interesting in this respect is the Kitaev model2 of S = 1/2 spins on a two-dimensional (2D) honeycomb lattice with nearest neighbours interacting through an Ising exchange, whose axis depends on the bond direction, as shown in Fig. 1a. This is one of a few exactly solved 2D models supporting a QSL ground state. According to the solution, a spin flip fractionalizes into a pair of gauge fluxes and a Majorana fermion2,3. As both types of quasiparticle behave as anyons under the magnetic field, they could potentially be used for decoherence-free topological quantum computation2. The experimental detection of such anyons is thus the primary goal of current QSL research.

Fig. 1: Structure of α-RuCl3 and the key signature of anyons.

a, The structure of a single layer of α-RuCl3 with the monoclinic C2/m (no. 12) space group and monoclinic axis b (c*⊥a, b). Spin-1/2 Ru3+ ions (red spheres) at the centres of the edge-sharing RuCl6 octahedra (grey) form an almost perfect honeycomb lattice. Ising axes of the Kitaev exchange interactions between nearest-neighbouring spins are perpendicular to the bond directions, pointing along x, y or z for blue, green and orange bonds, respectively. Red arrows show the employed magnetic field directions (described by the angle ϑ from the a–b plane) with respect to the oblate Ru3+g-tensor (yellow ellipsoid) of axial symmetry around c*. The field directions form a fan (red semicircle) perpendicular to the a–b plane, at 15° from the b axis (inset). b, Phase diagram of α-RuCl3 as a function of temperature T and the effective magnetic field Bab = g(ϑ)B/gab (so that Bab = B for ϑ = 0°) for various values of B and ϑ. The boundary of the magnetically ordered phase extending up to Bc ≈ 8 T, obtained from the 35Cl line width δν(T) (inset) and T1-1(T) (Fig. 3), matches the result of ref. 25 (grey line). The error bars represent one standard deviation. c, Δ (obtained from fits in Fig. 3) as a function of Bab follows the theoretically predicted cubic dependence given by equation (2) (grey line) with a finite zero-field offset (corresponding to the two-flux gap Δ0 = 0.065JK)6 leading to JK = 183 ± 10 K. The error bars reflect the change of the determined Δ values upon adding or removing a single edge T1-1(T) data point. Inset shows Δ(ϑ) for B = 9.4 T together with the curve obtained from the theoretically predicted Δ(Bab) (grey line). The only field direction outside the red fan in a is assigned ϑ = 180°.

As fractional quasiparticles are always created in groups, their common signature is a continuous spin excitation spectrum, observed in recent QSL candidates on the kagome and triangular lattices10,11, instead of sharp magnon modes found in ordered magnets. A Kitaev QSL also exhibits this feature6,12, as well as additional, specific signatures, all related to the fact that fractionalization in this case leads to two types of quasiparticle. First, the fractionalization proceeds in two steps, with both types releasing their entropy at different temperatures13. Second, although Majorana fermions themselves are gapless in zero magnetic field, the response of the QSL to a spin flip is gapped due to the inevitable simultaneous creation of a pair of gapped gauge fluxes6. Third, in the presence of an external magnetic field, the Majorana fermions also acquire a gap, which is predicted to grow with the characteristic third power of the field in the low-field region2,4,5. Currently, α-RuCl3 is the most promising candidate for the realization of the Kitaev QSL12,14,15,16,17. Among the listed signatures, a spin excitation continuum was observed by Raman spectroscopy12,14 and inelastic neutron scattering15,16,17, and the two-step thermal fractionalization was confirmed by specific-heat measurements17, all in zero field. However, an application of a finite field, which should affect the gaps of both types of quasiparticle in different ways, is crucial to identify them. Using nuclear magnetic resonance (NMR), we determine the field dependence of the spin excitation gap Δ shown in Fig. 1c, which indeed exhibits a finite zero-field value predicted for gauge fluxes and the cubic growth predicted for Majorana fermions. This result clearly demonstrates the fractionalization of a spin flip into two types of anyon in α-RuCl3.

α-RuCl3 is structurally related to the other two Kitaev QSL candidates, Na2IrO3 (ref. 18) and α-Li2IrO3 (ref. 19). All three are layered Mott insulators based on the edge-sharing octahedral units, RuCl6 and IrO6 (Fig. 1a), respectively, and driven by strong spin–orbit coupling20, which together lead to a dominant Kitaev exchange coupling between the effective S = 1/2 spins of Ru3+ and Ir4+ ions, respectively21. A monoclinic distortion of the IrO6 octahedra in both iridate compounds results in non-Kitaev exchange interactions between the spins, which lead to the low-temperature magnetic ordering and thus prevent the realization of the QSL ground state. Judging by the lower transition temperature with respect to the Kitaev exchange coupling JK, these interactions are smaller in α-RuCl3 (refs 22,23,24). Signatures of fractional quasiparticles should be sought in a region of the phase diagram outside the magnetically ordered phase, at temperatures low enough that the Kitaev physics is not yet destroyed by thermal fluctuations. This is the Kitaev paramagnetic phase (Fig. 1b) extending to a relatively high temperature of around 100 K, roughly half of JK = 190 K (ref. 17).

The boundary of the magnetically ordered phase measured in a large α-RuCl3 single crystal (see Methods) using 35Cl NMR is displayed in Fig. 1b. The magnetic response of α-RuCl3 is known to be highly anisotropic23,25, mainly due to the anisotropic Ru3+g-tensor (Fig. 1a) with the components gab = 2.5 in the a–b plane and gc*=1.1 along the c* axis (ref. 26). Namely, the effect of an applied magnetic field B is described by the Zeeman term, which is proportional to gB with magnitude g(ϑ)B, where g is the g-tensor and g(ϑ)=gab2cos2ϑ+gc*2sin2ϑ is the direction-dependent g-factor. Therefore, if B is applied at an angle ϑ from the a–b plane, it produces the same effect as the field with magnitude Bab = g(ϑ)B/gab applied in the a–b plane. We exploit this to scan efficiently the phase diagram as a function of Bab by varying both the direction and the magnitude of an applied field, instead of varying only the magnitude of a field applied in the a–b plane, as is usually done25. This approach is valid if the g-tensor is the only source of anisotropy, a condition to be verified at the end (Fig. 1c). As shown in the inset of Fig. 1b, we determine the transition temperature TN2 as the onset of NMR line broadening (see Methods) monitored on the dominant NMR peak (inset of Fig. 2d). The obtained phase boundary extending up to the critical field Bc ≈ 8 T matches the result of a recent reference study25. The observed transition temperature TN2 of around 14 K near zero field is consistent with a considerable presence of the two-layer AB stacking in the monoclinic C2/m crystal structure (Fig. 1a), in addition to the three-layer ABC stacking, which is characterized by a lower transition temperature TN1 of around 7 K in zero field15,24. As our study is focused on the Kitaev paramagnetic region (Fig. 1b) governed by the physics of individual layers, it is not affected by the particular stacking type.

Fig. 2: Signature of the Kitaev spin excitations.

a,b, T1-1 and susceptibility χ as a function of temperature T theoretically calculated for the ferromagnetic Kitaev model in zero field in ref. 28 (see Methods). Dotted grey line is a temperature-independent fit in the classical paramagnetic phase. Solid grey line is a fit to equation (1) with Δ = Δ0 = 0.065JK in the Kitaev paramagnetic phase between TL and TH leading to n = 0.67. First inset shows the goodness of the fit χ2 as a function of the power p using the generalized form T1-1∝Tpexp(nΔ∕T), where the minimum is almost exactly at p = −1 used in equation (1). Second inset demonstrates the resulting linear dependence of lnT1-1T on T−1 with the slope −nΔ between TH-1 and TL-1. c, Temperature dependence of 35Cl T1-1 taken on the dominant NMR peak (inset of d) in 9.4 T for two magnetic field orientations, ϑ = 90° (Bab = 4.1 T) and ϑ = 0° (Bab = 9.4 T). The dataset for ϑ = 90° exhibits two transitions into magnetically ordered states, at TN2 = 12 K (AB stacking) and at TN1 = 8 K (ABC stacking). Dashed red lines are fits to T1-1∝T2exp-Δm∕T for gapped magnon excitations in the three-dimensional magnetically ordered state. Dotted blue and red lines are temperature-independent fits in the classical paramagnetic phase (see Methods). Solid blue and red lines are fits to equation (1) with n = 0.67 for Kitaev spin excitations up to TH ≈ 70 K (and down to T* ≈ 17 K for ϑ = 90°, see Methods). Inset demonstrates the linear dependence of lnT1-1T on T−1 (obtained from equation (1)) in the appropriate T−1 range. d, Temperature-dependent 35Cl magnetic shift (NMR shift with subtracted quadrupole shift, see Methods) of the dominant NMR peak (inset shows the whole central line as a function of frequency ν with respect to the Larmor frquency νL = 39.18 MHz) in 9.4 T for two field orientations. The dataset for ϑ = 0° increases on cooling, approaching the maximum at 5 K, in accordance with the theoretical result in b where the maximum is reached at 0.017JK ≈ 3.2 K. The dataset for ϑ = 90° is almost temperature-independent down to TN2, in line with the temperature dependence of susceptibility for ϑ = 90° (ref. 22).

To detect and monitor the spin excitation gap as a function of the magnetic field, we use the NMR spin–lattice relaxation rate T1-1, which directly probes the low-energy limit of the local spin–spin correlation function and thus offers a direct access to the spin excitation gap. Figure 2a shows an exact theoretical T1-1(T) dependence numerically calculated for the ferromagnetic Kitaev model in zero field27,28 adapted to the case of α-RuCl3 (see Methods). It is dominated by a broad maximum, which is a sign of thermally excited pairs of gauge fluxes over the two-flux gap27, whose exact value amounts to Δ0 = 0.065JK (refs 2,6). The Kitaev paramagnetic phase is located between the lower temperatureTL = 0.012JK, where gapped gauge fluxes start to be excited, and the higher temperature TH = 0.375JK, where thermal fluctuations already destroy the short-range spin correlations13,27. As shown in Fig. 2a, in this broad temperature range, the theoretical zero-field T1-1(T) dataset can be nicely reproduced with the empirical expression

T1-1∝1Texp-nΔT

(1)

where Δ is the field-dependent spin excitation gap, and the fitting procedure gives n = 0.67 when setting Δ = Δ0 for zero field. In the following, we use equation (1) with this value of n to extract Δ from the experimental T1-1(T) datasets. Due to the unusual prefactor T−1 on the right-hand side of equation (1), Δ is encoded in the specific concave shape of the curve with the maximum at a temperature nΔ. As shown in the second inset of Fig. 2a, the gap is also directly accessible from the negative slope −nΔ of the linear dependence of lnT1-1T on T−1.

A broad maximum characteristic of the Kitaev spin excitations indeed appears in the 35Cl T1-1(T) dataset recorded in 9.4 T for the magnetic field orientation ϑ = 0° (Bab = 9.4 T), as shown in Fig. 2c. The dataset is excellently reproduced with equation (1) using Δ = 46.8 K in the temperature range up to TH = 0.375JK ≈ 70 K for JK = 190 K (ref. 17). In the T1-1(T) dataset for ϑ = 90° (Bab = 4.1 T), a maximum would apparently develop at a lower temperature, if the dataset was not disrupted by a magnetic ordering transition at TN2 = 12 K (Fig. 2c). Between T* ≈ 17 K, where the critical fluctuations preceding the magnetic ordering vanish (see Methods and Supplementary Information), and TH, the dataset is perfectly reproduced with equation (1) using Δ = 13.8 K. A large difference between the two determined gaps in Fig. 2c points to a considerable Δ(Bab) variation in the Kitaev paramagnetic phase. Below TN2, two T1-1 components develop, both exhibiting a steep drop, one below TN2 and the other one below TN1 = 8 K. These two phase transitions were observed before and ascribed to the presence of AB and ABC stackings, respectively15,24. Fitting the data below TN2 and TN1 with the expression T1-1∝T2exp-Δm∕T valid for gapped magnon excitations in a three-dimensional magnetically ordered state (see Methods) gives comparable values of the magnon gap Δm = 32 K and 35 K, respectively, implying the same low-energy physics in both cases. The obtained values are consistent with the gap of 29 K determined by inelastic neutron scattering15,29. Finally, the temperature-independent part of both T1-1(T) datasets above 120 K (≈2TH) indicates a crossover into a classical paramagnetic state (see Methods), in accordance with the theoretical dataset in Fig. 2a and with the experimental result of ref. 17.

Although equation (1) is empirical, its functional form already bears signs of the involved fractional spin excitations. To show this, we first note that similar expressions are obtained for more conventional gapped magnon excitations in magnetic insulators at low temperatures T≪Δm (see Methods). In that case, the prefactor T−1 is replaced by a more general Tp originating from the magnon density of states g(E) (where E is the energy of magnons), which depends on the dimensionality D, while n is generally the number of magnons involved in the process. For n = 1 (single-magnon scattering) and a quadratic dispersion relation for magnons, one obtains p = D − 1 ≥ 0, while higher n (multi-magnon scattering) lead to even higher powers p (see Methods). At higher temperatures T ≈ Δm, the effective p changes, but always remains positive. Therefore, a very unusual, theoretically confirmed p = −1 in equation (1) valid over a broad temperature range cannot be obtained for magnons. This, together with a fractional number n = 0.67 of the involved spin-flip excitations obtained in Fig. 2a, indicate that equation (1) is actually specific to the presence of fractional spin excitations.

In the ferromagnetic Kitaev model, the magnetic susceptibility χ exhibits a moderate, almost monotonic temperature dependence (Fig. 2b), in contrast to the non-monotonic T1-1(T) behaviour (Fig. 2a)27,28. As shown in Fig. 2c,d, we indeed observe contrasting temperature dependences of T1-1 and the local susceptibility monitored by the 35Cl NMR shift in 9.4 T for ϑ = 0°. While both observables should exhibit a similar gapped behaviour in the presence of gapped magnon excitations, such a contrast between them can arise only in the presence of at least two types of fractional quasiparticle that enter the two observables in different ways27.

To obtain Δ as a function of Bab in Fig. 1c, the T1-1(T) datasets in Fig. 3 taken in magnetic fields of different directions and magnitudes are fitted to equation (1) in the temperature range of the Kitaev paramagnetic phase. As the curve T1-1∝T-1 defined by Δ = 0 is steeper than any dataset in this range or, equivalently, as the datasets in the insets of Fig. 3 all exhibit a negative slope in this range, the obtained excitation gaps are obviously all finite. The inset of Fig. 1c showing the symmetric Δ(ϑ) dependence around 90° in 9.4 T, where ϑ traverses non-equivalent directions with respect to the Ising axes on both sides (inset of Fig. 1a), demonstrates that the g-tensor is indeed the only source of anisotropy as assumed when introducing Bab. The obtained Δ(Bab) in Fig. 1c can be perfectly reproduced as a sum of two terms: the two-flux gap Δ0 = 0.065JK (refs 2,6) and the gap acquired by Majorana fermions in a weak magnetic field, theoretically predicted to be proportional to the cube of the field2,4,5 (see Methods)

Δ=Δ0+αB̃3Δ02

(2)

where B̃=gabμBBab∕kB is the field in kelvin units, kB is the Boltzmann constant, μB is the Bohr magneton and α accounts for the sum over the excited states in the third-order perturbation theory, which is the origin of the B̃3 term2. The fitting procedure leads to α = 1.2 and JK = 183 ± 10 K, in perfect agreement with the value of 190 K determined by inelastic neutron scattering17. This result demonstrates that a spin-flip excitation in α-RuCl3 indeed fractionalizes into a gauge flux pair and a Majorana fermion.

Fig. 3: Determination of the spin excitation gap.

a,b, 35Cl T1-1 as a function of temperature T in 2.35 T, 4.7 T and 9.4 T for various magnetic field orientations given by ϑ. Arrows mark the transition temperatures TN2 into the magnetically ordered state, determined by a weakly pronounced onset of a T1-1 decrease on decreasing T. Solid lines are fits to equation (1) with n = 0.67 for Kitaev spin excitations between T*, where the critical fluctuations related to magnetic ordering vanish, and TH = 70 K (blue background). These allow us to determine Δ(ϑ) and Δ(Bab) dependences shown in Fig. 1c. Dashed lines are the curves T1-1∝T-1 defined by Δ = 0 exhibiting the largest negative slope. Insets demonstrate the linear dependence of lnT1-1T on T−1 (obtained from equation (1)) in the appropriate T−1 range (blue background). Dashed horizontal lines correspond to Δ = 0.

Focusing on the low-field Kitaev paramagnetic region in the phase diagram of α-RuCl3 in Fig. 1b is essential for our identification of two types of anyon. Instead, other recent experimental studies focused on the low-temperature region above Bc, observing the spin excitation continuum30 with either a gapless behaviour31 or the gap opening linearly32,33,34,35 or sublinearly36 with B − Bc, but without definite conclusions about the identity of the involved quasiparticles. These contradictory conclusions probably originate from the presence of additional, non-Kitaev interactions between the spins15,26,29,37, whose role should be pronounced particularly at low temperatures. Our result shows that spin fractionalization into two types of anyon is robust against these interactions in a broad range of temperatures and magnetic fields. This is the main practical advantage of α-RuCl3 with respect to all other anyon realizations such as the fractional quantum Hall effect in 2D heterostructures7 or hybrid nanowire devices38 where anyons are observed only at extremely low temperatures and for certain field values. Our discovery thus establishes α-RuCl3 as a unique platform for future investigations of anyons.

Methods

Crystal growth

Crystals of α-RuCl3 were synthesized from anhydrous RuCl3 (Strem Chemicals). The starting material was heated in vacuum to 200 °C for one day to remove volatile impurities. In the next step, the powder was sealed in a silica ampoule under vacuum and heated to 650 °C in a tubular furnace. The tip of the ampoule was kept at lower temperature and the material sublimed to the colder end during one week. Phase pure α-RuCl3 (with a high-temperature phase of C2/m crystal structure) was obtained as thin crystalline plates. The residual in the hot part of the ampoule was black RuO2 powder. The purified α-RuCl3 was sublimed for the second time to obtain bigger crystal plates. The phase and purity of the compounds were verified by powder X-ray diffraction. All handling of the material was done under strictly anhydrous and oxygen-free conditions in glove boxes or sealed ampoules. Special care has to be taken when the material is heated in sealed-off ampoules. If gas evolves from the material, this may result in the explosion of the ampoule.

NMR

The 35Cl NMR experiments were performed on a foil-like α-RuCl3 single crystal of approximate dimensions 5 × 5 × 0.1 mm3 in a continuous-flow cryostat, allowing us to reach temperatures down to 4.2 K. When handling the sample, we took extreme care to minimize its exposure to air. A thin NMR coil tightly fitting the sample was made from a thin copper wire with 20–40 turns, depending on the required tuning frequency determined by the external magnetic field. The coil was covered with a mixture of epoxy and ZrO2 powder, which was allowed to harden to ensure the rigidity of the coil. The coil was then mounted on a Teflon holder attached to a rotator, which allowed us to vary the orientation of the sample with respect to the external magnetic field. To reduce the noise of an already weak 35Cl NMR signal, a consequence of the extremely broad 35Cl NMR spectrum, we used a bottom-tuning scheme. With the output radio frequency power of around 20 W, a typical π/2 pulse duration was 2 μs. The NMR signals were recorded using the standard spin-echo π/2 − τd − π pulse sequence with a typical delay of τd = 70 μs (much shorter than the spin–spin relaxation time T2) between the π/2 and π pulses.

T1 relaxation

The spin–lattice relaxation (T1) experiment was carried out using an inversion recovery pulse sequence, φi − τ − π/2 − τd − π, with an inversion pulse φi < π (suitable for broad NMR lines) and a variable delay τ before the read-out spin-echo sequence. The spin–lattice relaxation datasets were typically taken at 20 increasing values of τ. The datasets were analysed using the model of magnetic relaxation for nuclear spin I = 3/2, appropriate for 35Cl, monitored on the central −1/2 ↔ 1/2 transition:

m(τ)=1-(1+s)0.1exp-τT1+0.9exp-6τT1

(3)

where T1 is the spin–lattice relaxation time and s is the inversion factor. In the region of the phase diagram outside the magnetically ordered phase (Fig. 1b), this expression reproduces the experimental relaxation curves perfectly. In the magnetically ordered phase, two T1 components appear, and the relaxation curves are reproduced as a sum of two terms of the form given by equation (3). For instance, the temperature dependence of the corresponding two T1 values for B = 9.4 T with ϑ = 90° is given in Fig. 2c. In cases where only a narrow temperature region below the transition was covered, the two components in the relaxation curves were hard to separate, and we used equation (3) with a stretched exponent instead.

Relation between orientation and field dependence of T1

Using the direction of the applied magnetic field (described by the angle ϑ from the crystal a–b plane) in addition to its magnitude to change the value of Bab allows us to cover low Bab values while keeping the applied magnetic field B high. This is beneficial for two reasons related to the strong quadrupole broadening of the 35Cl NMR spectrum (inset of Fig. 2d and Supplementary Fig. 2): to minimize an already large NMR line width, and to keep the Larmor frequency well above the quadrupole splitting, which is of the order of 10 MHz as concluded in the following. The validity of this approach is supported by the fact that Δ(Bab) data points for various angles ϑ and field values 2.35 T, 4.7 T and 9.4 T in Fig. 1c all collapse onto a smooth experimental curve. The Δ(Bab) data points taken in lower fields obviously exhibit much larger error bars. Namely, the corresponding T1-1(T) datasets in Fig. 3b are more scattered than the datasets taken in 9.4 T despite a much longer averaging for noise reduction.

The boundary of the magnetically ordered phase

We measured the temperature dependence of the dominant 35Cl NMR peak in magnetic fields of various magnitudes and various directions with respect to the crystal a–b plane of the sample, thus covering various Bab values. For Bab < 8 T, the frequency width of the peak exhibits a clear kink as a function of temperature (inset of Fig. 1b and Supplementary Fig. 3), which indicates the onset of NMR line broadening at the phase transition into the magnetically ordered state. Plotting the temperature of the kink as a function of Bab in Fig. 1b (and in the inset of Supplementary Fig. 3), we obtain the boundary of the magnetically ordered phase, which perfectly matches the result of the reference study25.

Contributions to the NMR shift

The temperature dependence of the NMR frequency shift of the dominant NMR peak measured in 9.4 T with the field in the a–b plane (ϑ = 0°) is plotted in Fig. 2d and reproduced in Supplementary Fig. 4. To separate the magnetic contribution to the NMR shift from the temperature-independent quadrupole contribution, we plot the relative NMR shift (the NMR shift divided by the 35Cl Larmor frequency νL = 39.18 MHz) against the rescaled magnetic susceptibility χab in the inset of Supplementary Fig. 4. In ref. 25, an experimental ratio between the susceptibility χ of the powdered sample and the susceptibility χab of the single crystal with a field applied in the a–b plane is obtained as (2 + r)/3 with r = 0.157, leading to χab = 3χ/(2 + r). We use this empirical relation to evaluate χab(T) from our field-cooled χ(T) dataset shown in Supplementary Fig. 1. As we did not measure susceptibility in high magnetic fields, we rely on the dataset taken in 1.0 T. This is valid in a broad temperature range, except at low temperatures where this dataset starts to deviate from the high-field susceptibility22. The relative shift is found to depend linearly on χab up to 20 × 10−3 e.m.u. mol−1 (inset of Supplementary Fig. 4), that is, the NMR shift follows χab down to 35 K (Supplementary Fig. 4). The proportionality constant 1.95 mol e.m.u.−1 between the relative shift and χab translates to 1.09 T/μB. When further multiplied by the relevant g-factor, gab = 2.5, and divided by 2 (as each 35Cl is coupled equally to two Ru3+S = 1/2 spins), we obtain the component A = 1.35 T (in the a–b plane) of the hyperfine coupling tensor between 35Cl and the Ru3+S = 1/2 spin. In addition, a zero-temperature relative shift −0.039, when multiplied by νL, gives the quadrupole shift ΔνQ = −1.53 MHz.

From the obtained quadrupole shift, we can estimate the quadrupole splitting νQ between the successive 35Cl NMR transitions. For the case of an axially symmetric electric field gradient (EFG) tensor and the field applied at an angle ϑ′ from the principal EFG axis vZZ with the largest EFG eigenvalue, the second-order quadrupole shift is given by ΔνQ=-3νQ2(1 − cosϑ′2)(9cosϑ′2 − 1)/(16νL) for the I = 3/2 nucleus39. As the axes of the EFG tensor are not known, we assume a 45° typical tilt of vZZ from c*, so that ϑ′ ≈ 45°. From the previously evaluated ΔνQ, we then obtain νQ ≈ 14.2 MHz. This is an estimate of the quadrupole splitting between the central 35Cl NMR transition and the satellite transitions. We can thus conclude that the NMR peaks in the covered frequency range (inset of Fig. 2d and Supplementary Fig. 2) all belong to the central transition.

Justification for the onset temperature T*

A simple estimate of the onset temperature T* for the critical spin fluctuations preceding the magnetic ordering at TN2 can be obtained using a random phase approximation to express the dynamic susceptibility χ of α-RuCl3. As the pure Kitaev magnet does not magnetically order, the magnetic ordering of α-RuCl3 is a consequence of additional non-Kitaev interactions. If we denote the effective coupling of all these additional interactions by J′, the random phase approximation reads χ = χ1/(1 − J′χ1), where χ1 is the dynamic susceptibility of a single Kitaev plane. As T1-1(T) can be reasonably well approximated by a negative power law above T* (as evident from the ϑ = 90° dataset in Fig. 2c, which is linear in a log–log scale), the same holds for χ1, that is, χ1 = cT−r with r > 0, where c is the proportionality constant. A condition J′χ1(TN) = 1 for the phase transition at TN can then be written as cJ′TN-r=1 This allows us to express J′ with TN, so that the deviation of the dynamic susceptibility from the pure Kitaev case can finally be written in the form

χχ1=1-TTN-r-1

(4)

The ratio χ/χ1 diverges when approaching TN from above and gradually drops to 1 with increasing temperature. The temperature scale for this drop is apparently determined by TN, which leads to the characteristic estimate T* ≈ TN + TN = 2TN consistent with our experimental determination (see Supplementary Information). Precisely such values of T* with respect to TN are also experimentally found in systems of coupled spin chains or ladders, whose χ1 values exactly obey power-law temperature dependences40,41.

Theoretical T1-1(T) for Kitaev spin excitations

The theoretical temperature dependence of T1-1 is numerically calculated for the Kitaev model in zero field in ref. 27. T1-1 contains two contributions, one coming from a single fluctuating spin (on-site) and the other one coming from fluctuating nearest-neighbouring spins in the Kitaev honeycomb lattice. As the 35Cl nucleus in α-RuCl3 is located at equal distances from the closest two Ru3+S = 1/2 spins (Fig. 1a), T1-1 generally contains both contributions. We evaluate their relative weights for the case of a magnetic field applied along c* (ϑ = 90°) and assuming an isotropic hyperfine coupling A. A general expression for T1-1(ref. 42) can then be written as

T1-1=γ22∫-∞∞dtA2S1a(t)+S2a(t)S1a+S2a+A2S1b(t)+S2b(t)S1b+S2b

(5)

where γ is the nuclear gyromagnetic ratio (26.2 MHz T–1 for 35Cl), t is time, while S1 and S2 are the two involved Ru3+S = 1/2 spins with relevant components perpendicular to the field direction, that is, along a and b (Fig. 1a). In the orthogonal system defined by the Ising axes x, y and z (Fig. 1a), unit vectors along a, b and c* are written as

êa=1611-2,êb=12-110,êc*=13111

(6)

This allows us to switch from the spin components along a and b in equation (5) to the spin components along x, y and z, which are appropriate for the Kitaev model. Taking into account also the isotropy of the Kitaev model, that is, S1x(t)S1x = S1y(t)S1y = S1z(t)S1z and similar for the correlation functions S1x(t)S2x and S2x(t)S2x, equation (5) simplifies to

T1-1=12γ2A2∫-∞∞dt4S1x(t)S1x+43S1x(t)S2x

(7)

Two terms under the integral represent precisely the on-site and nearest-neighbouring-sites contributions with the weights 4 and 4/3, respectively. A similar calculation for the magnetic field applied along a or b gives the same result. The theoretical curve plotted in Fig. 2a is based on these weights. As the curves for the on-site and nearest-neighbouring-sites contributions are almost identical up to TH (ref. 27), the values of the weights do not affect the form of equation (1).

T1-1 in the classical paramagnetic state

In this state, nearest-neighbouring spins are not correlated, so that only the on-site spin correlations contribute to T1-1. A contribution of the correlation function S1x(t)S1x amounts to T1x-1 = πγ2A2ℏ∕4Jz (ref. 43), where ħ is Planck’s constant divided by 2π, J is the exchange coupling (in this case, JK = 183 K, as determined in Fig. 1c) and z is the number of neighbouring spins with the relevant component of the coupling (x in this case, z = 1 for the Kitaev honeycomb lattice). Using equation (7), we then obtain T1-1 = 4T1x-1 = πγ2A2ℏ∕JK, which evaluates to 92 s−1 using A = 1.35 T determined above from Supplementary Fig. 4. This value is consistent with the observed values of 121 s−1 and 146 s−1 above 120 K (Fig. 2c) for ϑ = 0° and ϑ = 90°, respectively. In addition, these temperature-independent values are reached at around 2TH (TH = 0.375JK ≈ 70 K), precisely as in the theoretical dataset in Fig. 2a. Both these conclusions demonstrate that the classical paramagnetic state is indeed reached.

T1 relaxation due to gapped magnons

When spin fluctuations in the magnetic lattice are due to excited magnons, the corresponding spin–lattice relaxation rate for a single-magnon process is given by44

T1-1∝∫g2(E)f(E)[1+f(E)]dE

(8)

where f(E) = [exp(βE) − 1]−1 is the Bose–Einstein distribution function and β = 1/(kBT). Denoting the magnon gap by Δm (in kelvin units), we define ε = E − kBΔm as the energy measured from the bottom of the magnon band. The power-law dispersion relation ε∝ks in D dimensions, which includes the standard parabolic dispersion (s = 2) and the Dirac dispersion (s = 1) as special cases, leads to g(E) ∝εD/s − 1. For low temperatures T≪Δm, f(E) can be approximated with the Boltzmann distribution, f(E) ≈ exp(−βE) = exp(−Δm/T)exp(−βε). Plugging these expressions for g(E) and f(E) into equation (8), we obtain

T1-1∝T2D∕s-1exp-ΔmT∫0∞exp(-x)x2(D∕s-1)dx

(9)

The integral on the right-hand side of equation (9) converges if s < 2D and evaluates to Γ(2D/s − 1) where Γ is the gamma function. We can thus rewrite equation (9) as

T1-1∝Tpexp-ΔmT

(10)

with the power of the prefactor p = 2D/s − 1. In the case of D = 2, which is relevant for the Kitaev honeycomb magnet, we obtain p = 1 for s = 2, and p = 3 for s = 1, so that the power p cannot be negative. Even in the case of D = 1, p can only approach the lowest value of 0 precisely for s = 2 (although care should be taken in this case, as the integral in equation (9) then formally diverges). If more than a single magnon is involved in the T1 process, the power p is also positive and becomes even higher44. Gapped magnons thus cannot lead to the T1 relaxation described by equation (10) with p < 0, which is the case in equation (1).

Instead, we can use equation (9) in the three-dimensional magnetically ordered state, when the elementary excitations are indeed magnons with a gap Δm. In this case, D = 3 and s = 2, and this leads to T1-1∝T2exp-Δm∕T. We use this expression to analyse the T1-1(T) data (Fig. 2c) in the low-temperature ordered state of α-RuCl3.

All these examples show that a frequently used simple gapped model T1-1∝exp(-Δs∕T) with the gap Δs, which was used before to analyse the T1-1(T) datasets in α-RuCl3 (ref. 32), is actually not justified in any region of the phase diagram of α-RuCl3.

Majorana fermion gap

In the presence of an external magnetic field, Majorana fermions in the Kitaev model acquire a gap2. This is shown for a field applied perpendicularly to the honeycomb plane, that is, in the (111) direction in the coordinate system defined by the Kitaev axes x, y and z. In this case, the Zeeman term reads HZ=-h∑jSjx+Sjy+Sjz, where h=gμBB∕3 is a single component of the magnetic field B in energy units and g is the g-factor. When treated as a perturbation to the Kitaev Hamiltonian, the Zeeman term contributes to the Majorana fermion gap only at third order2. The corresponding effective Hamiltonian is thus proportional to h3 and can be written as2,4,5

Heff(3)=-3αh3kB2Δ02 ∑jklSjxSkySlz

(11)

where Δ0 is the two-flux gap (in kelvin units), while α (of the order of unity) accounts for the sum over the excited states, and its exact value is not known. The Kitaev model extended with such a three-spin exchange term -κ∑jklSjxSkySlz with κ=3αh3∕kB2Δ02 is still exactly solvable, and the dispersion relation of the Majorana fermions is calculated as4

Ek=2kB2JK21+eik⋅a1+eik⋅a22+κ2sin2k⋅a1

(12)

where i is the square root of −1, k is the Majorana fermion momentum, and a1 and a2 are the base vectors of the honeycomb lattice. The dispersion relation given by equation (12) is gapped for κ ≠ 0, and the corresponding Majorana fermion gap Δf can be calculated numerically as a function of κ and thus as a function of the magnetic field. For small magnetic fields, that is, for κ≪kBJK, the Majorana fermion gap (in kelvin units) simplifies to

Δf=3κkB=αΔ02gμBBkB3

(13)

while for high magnetic fields, it saturates to Δf = 2JK. The total spin excitation gap Δ is obtained by adding Δf to the two-flux gap Δ0 as in equation (2). The whole field dependence of Δ is shown in Supplementary Fig. 8 for JK = 183 K as determined in Fig. 1c, g = gab and α = 1.2 (leading to the best fit of our Δ(Bab) data points). The cubic approximation given by equation (13), which is plotted in Fig. 1c and, for comparison, in Supplementary Fig. 8, is apparently valid up to 15 T, well beyond the field range covered in this work.

Data availability

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Search for Christian Rüegg in:

Search for Martin Klanjšek in:

Contributions

M.K. conceived, designed and led the project. N.J. and M.K. performed the NMR experiments and analysed the data. K.W.K. and D.B. grew the samples. A.B. performed the magnetic susceptibility measurements. All the authors discussed the results. M.K. wrote the paper with feedback from all the authors.