Observation of two types of anyons in the Kitaev honeycomb magnet

TL;DR

This paper reports the experimental observation of two types of anyons, gauge fluxes and Majorana fermions, in the Kitaev honeycomb magnet $ ext{α}$-RuCl$_3$, confirming fractionalization predicted by the Kitaev model.

Contribution

The study demonstrates fractionalization into gapped anyons in a real material, showing survival of these excitations under various conditions despite non-Kitaev interactions.

Findings

Majorana fermion contribution scales with the cube of magnetic field

Finite zero-field gauge-flux gap matches predictions

Fractionalization persists across a broad temperature and magnetic field range

Abstract

Quantum spin liquid is a disordered magnetic state with fractional spin excitations. Its clearest example is found in an exactly solved Kitaev honeycomb model where a spin flip fractionalizes into two types of anyons, quasiparticles that are neither fermions nor bosons: a pair of gauge fluxes and a Majorana fermion. Here we demonstrate this kind of fractionalization in the Kitaev paramagnetic state of the honeycomb magnet $\alpha$-RuCl$_3$. The spin-excitation gap measured by nuclear magnetic resonance consists of the predicted Majorana fermion contribution following the cube of the applied magnetic field, and a finite zero-field contribution matching the predicted size of the gauge-flux gap. The observed fractionalization into gapped anyons survives in a broad range of temperatures and magnetic fields despite inevitable non-Kitaev interactions between the spins, which are predicted to drive the system towards a gapless ground state. The gapped character of both anyons is crucial for their potential application in topological quantum computing.

Equations

$$T_1^{-1}\propto \frac{1}{T}\exp \left(-\frac{n\Delta }{T}\right),$$

$$\Delta ={\Delta }_0+\frac{\alpha }{3}\frac{{\tilde{B}}^3}{{\Delta }_0^2},$$

$$m(\tau )=1-(1+s)\left[0.1\exp \left(-\frac{\tau }{T_1}\right)+0.9\exp \left(-\frac{6\tau }{T_1}\right)\right],$$

$$\displaystyle\Bigl{<}A$$

$$\displaystyle=A^{2}\Bigl{[}\bigl{<}S_{1}(t)S_{1}\bigr{>}+\bigl{<}S_{2}(t)S_{2}\bigr{>}\,\pm$$

$$\displaystyle\pm\Bigl{(}\bigl{<}S_{1}(t)S_{2}\bigr{>}+\bigl{<}S_{2}(t)S_{1}\bigr{>}\Bigr{)}\Bigr{]},$$

$$T_{1}^{-1}\propto\int g^{2}(E)n(E)\bigl{[}1+n(E)\bigr{]}{\rm d}E,$$

$$T_1^{-1}\propto T^{2D/s-1}\exp \left(-\frac{\Delta }{T}\right){\int }_0^{\infty }\exp (-x)x^{2(D/s-1)}\mathrm{d}x.$$

$$T_1^{-1}\propto T^p\exp \left(-\frac{\Delta }{T}\right)$$

$${\mathcal{H}}_{\mathrm{eff}}^{(3)}=-\alpha \frac{h^3}{k_B^2{\Delta }_0^2}\sum _{jkl}{S}_j^x{S}_k^y{S}_l^z,$$

$$E_{\mathbf{k}}=2\sqrt{k_B^2{J}_K^2|1+e^{i\mathbf{k}\cdot {\mathbf{a}}_1}+e^{i\mathbf{k}\cdot {\mathbf{a}}_2}{|}^2+{\kappa }^2{\sin }^2(\mathbf{k}\cdot {\mathbf{a}}_1)},$$

$${\Delta }_f=\sqrt{3}\frac{\kappa }{k_B}=\frac{\alpha }{3{\Delta }_0^2}{\left(\frac{g{\mu }_{B}B}{k_B}\right)}^3,$$

Full text

In many-body systems dominated by strong fluctuations, an excitation with an integer quantum number can break up into exotic quasiparticles with fractional quantum numbers. Well known examples include fractionally charged quasiparticles in fractional quantum Hall effect dePicciotto_1997 , spin-charge separation in one-dimensional conductors Jompol_2009 , and magnetic monopoles in spin ice Castelnovo_2008 . A major hunting ground for novel fractional quasiparticles are 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$ Balents_2010 ; Savary_2017 . Even more interesting in this respect is the Kitaev model Kitaev_2006 of $S=1/2$ spins on a two-dimensional (2D) honeycomb lattice with nearest neighbors interacting through an Ising exchange, whose axis depends on the bond direction, as shown in Fig. 1(a). 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 bosonic gauge fluxes and a Majorana fermion Kitaev_2006 ; Baskaran_2007 . As both types of quasiparticles behave as anyons, i.e., neither bosons nor fermions, under exchange, they could be used for decoherence-free topological quantum computation Kitaev_2006 . The experimental detection of such anyons is thus the primary goal of current QSL research.

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 lattices Han_2012 ; Paddison_2016 , instead of sharp magnon modes found in ordered magnets. A Kitaev QSL also exhibits this feature Knolle_2014 ; Nasu_2016 , as well as additional, specific signatures, all related to the fact that fractionalization in this case leads to different quasiparticles. First, the fractionalization proceeds in two steps, with both types of quasiparticles releasing their entropy at different temperatures Nasu_2015 . 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 fluxes Knolle_2014 . And 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 region Kitaev_2006 ; Jiang_2011 ; Nasu_2017 . Currently, $\alpha$-RuCl3 stands as the most promising candidate for the realization of the Kitaev QSL Nasu_2016 ; Sandilands_2015 ; Banerjee_2016 ; Banerjee_2017 ; Do_2017 . Among the listed signatures, a spin-excitation continuum was observed by Raman spectroscopy Sandilands_2015 ; Nasu_2016 and inelastic neutron scattering Banerjee_2016 ; Banerjee_2017 ; Do_2017 , and the two-step thermal fractionalization was confirmed by specific-heat measurements Do_2017 , all in zero field. However, an application of a finite field, which should affect the gaps of both types of quasiparticles differently, is crucial to identify them. Using nuclear magnetic resonance (NMR), we determine the field dependence of the spin-excitation gap $\Delta$ shown in Fig. 1(c), 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 anyons in $\alpha$-RuCl3.

$\alpha$-RuCl3 is structurally related to the other two Kitaev QSL candidates, Na2IrO3 Choi_2012 and $\alpha$-Li2IrO3 Singh_2012 . All three are layered Mott insulators based on the edge-sharing octahedral units, RuCl6 and IrO6 [Fig. 1(a)], respectively, and driven by strong spin-orbit coupling Plumb_2014 , which together lead to a dominant Kitaev exchange coupling between the effective $S=1/2$ spins of Ru3+ and Ir4+ ions, respectively Jackeli_2009 . A monoclinic distortion of the IrO6 octahedra in both iridate compounds results in the presence of 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, these interactions are smaller in $\alpha$-RuCl3 Kubota_2015 ; Sears_2015 ; Majumder_2015 ; Cao_2016 . Signatures of fractional quasiparticles should thus be sought in a region of the phase diagram outside the magnetically ordered phase, at temperatures where the Kitaev physics is not yet destroyed by thermal fluctuations. This is the Kitaev paramagnetic region [Fig. 1(b)] extending to a relatively high temperature around $100$ K where the nearest-neighbor spin correlations vanish Do_2017 .

The boundary of the magnetically ordered phase measured in a large $\alpha$-RuCl3 single crystal Supplemental using 35Cl NMR is displayed in Fig. 1(b). Magnetic properties of $\alpha$-RuCl3 are known to be highly anisotropic Johnson_2015 ; Majumder_2015 , mainly because of the anisotropic Ru3+ $g$-tensor [Fig. 1(a)] with $g_{ab}=2.5$ and $g_{c^{*}}=1.1$ Yadav_2016 . We exploit this anisotropy to scan the phase diagram by varying the direction of the applied fixed field, instead of varying the magnitude of the field applied in the $ab$ plane Johnson_2015 . Namely, as the Zeeman term contains the product $gB$, a magnetic field $B$ applied at an angle $\vartheta$ from the $ab$ plane is equivalent to the effective field $B_{ab}=g(\vartheta)B/g_{ab}$ applied in the $ab$ plane, where $g(\vartheta)=\sqrt{g_{ab}^{2}\cos^{2}{\vartheta}+g_{c^{*}}^{2}\sin^{2}{\vartheta}}$ is the direction-dependent $g$-factor. This is valid if the studied underlying physics is close to isotropic, a condition to be verified at the end. As shown in the inset of Fig. 1(b), we determine the transition temperature $T_{N2}$ as the onset of NMR line broadening Supplemental monitored on the dominant NMR peak [inset of Fig. 2(b)]. The obtained phase boundary extending up to the critical field $B_{c}\approx 8$ T matches the result of a recent reference study Johnson_2015 . The observed transition temperature $T_{N2}$ 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. 1(a)], in addition to the three-layer $ABC$ stacking, which is characterized by a lower transition temperature $T_{N1}$ of around $7$ K in zero field Banerjee_2016 ; Cao_2016 . As our study is focused on the Kitaev paramagnetic region [Fig. 1(b)] governed by the physics of individual layers, it is not affected by the particular stacking type.

To detect and monitor the spin-excitation gap as a function of the magnetic field, we use the NMR spin-lattice relaxation rate $T_{1}^{-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. Fig. 2(a) shows the 35Cl $T_{1}^{-1}(T)$ datasets taken on the dominant NMR peak [inset of Fig. 2(b)] in $9.4$ T for two magnetic field orientations. A noticeable feature of the $T_{1}^{-1}(T)$ dataset for ${\bf B}\perp c^{*}$ (i.e., in the $ab$ plane, $B_{ab}=9.4$ T) is a broad maximum around $30$ K, followed by a steep decrease towards lower temperatures. In the $T_{1}^{-1}(T)$ dataset for ${\bf B}\parallel c^{*}$, such a feature would apparently develop at a lower temperature, if the dataset was not disrupted by the phase transition at $T_{N2}=12$ K [in a field of $B_{ab}=4.1$ T, Fig. 2(b)]. Instead, two $T_{1}^{-1}$ components develop below $T_{N2}$, both exhibiting a steep drop, one below $T_{N2}$ and the other one below $T_{N1}=8$ K. These two phase transitions were observed before and ascribed to the presence of $AB$ and $ABC$ stackings, respectively Banerjee_2016 ; Cao_2016 . The analysis of the data below $T_{N2}$ and $T_{N1}$ using the expression $T_{1}^{-1}\propto T^{2}\exp(-\Delta_{m}/T)$ valid for gapped magnon excitations in the 3D ordered state Supplemental gives comparable values of the magnon gap $\Delta_{m}=32$ K and $35$ K, respectively, implying the same low-energy physics in both cases. The obtained values are compatible with the gap of $29$ K determined by inelastic neutron scattering Banerjee_2016 ; Ran_2017 .

To access the key information held by the $T_{1}^{-1}(T)$ datasets in the Kitaev paramagnetic state, we first observe in Fig. 2(a) that the dataset for ${\bf B}\perp c^{*}$ below $100$ K exhibits the same shape as the theoretical dataset numerically calculated for the ferromagnetic Kitaev model in zero field Supplemental ; Yoshitake_2016 . A characteristic broad maximum of the latter is a sign of thermally excited pairs of gauge fluxes over the two-flux gap Yoshitake_2016 , whose exact value amounts to $\Delta_{0}=0.065J_{K}$ Kitaev_2006 ; Knolle_2014 where $J_{K}$ is the Kitaev exchange coupling. As shown in Fig. 2(a), a large part of the theoretical dataset, up to around $0.2J_{K}$, well above the maximum, can be excellently described by the phenomenological expression

[TABLE]

where $n$ is set to $0.61$ in order for $\Delta$ to match the required value of $\Delta_{0}$. This expression has a useful property: its maximum appears at a temperature $\Delta/n$, which allows for a simple estimate of $\Delta$ directly from the $T_{1}^{-1}(T)$ dataset. The ${\bf B}\perp c^{*}$ ($B_{ab}=9.4$ T) dataset up to $50$ K is excellently reproduced by Eq. (1) using $\Delta=51$ K [Fig. 2(a)]. The validity of Eq. (1) in this case and in the case of the theoretical dataset is clearly demonstrated in the inset of Fig. 2(a), which shows the resulting linear dependence of ${\ln}(T_{1}^{-1}T)$ on $T^{-1}$ below $50$ K. Meanwhile, even in the absence of the characteristic maximum, the ${\bf B}\parallel c^{*}$ ($B_{ab}=4.1$ T) dataset above $17$ K, i.e., slightly above $T_{N2}$, and up to high temperatures matches the theoretical zero-field dataset using the value $J_{K}=190$ K determined by inelastic neutron scattering, also based on the ferromagnetic Kitaev model Do_2017 . This means that the gap for $4.1$ T is already close to the zero-field value $\Delta_{0}=0.065J_{K}=12.4$ K. A large difference between the two determined gaps points to a significant $\Delta(B_{ab})$ variation in the Kitaev paramagnetic state. Finally, the temperature-independent part of both $T_{1}^{-1}(T)$ datasets above $100$ K indicates a crossover into the classical paramagnetic state Moriya_1956 , in line with the result of Ref. Do_2017 .

The expression given by Eq. (1) is not merely phenomenological, but reveals the presence of gapped fractional spin excitations. Similar expressions are obtained for the $T_{1}$ relaxation due to gapped magnons in magnetic insulators at low temperatures $T\ll\Delta$ Supplemental . In this case, the prefactor $T^{-1}$ is replaced by a more general $T^{p}$ originating from the magnon density of states $g(E)$, 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 $g(E)\propto E^{D/2-1}$ and thus $p=D-1\geq 0$, while higher $n$ (multi-magnon scattering) lead to even higher powers $p$ Supplemental . At higher temperatures $T\sim\Delta$, the effective $p$ changes, but always remains positive. As the very unusual $p=-1$ in Eq. (1) valid for $T\lesssim\Delta$ cannot be obtained for magnons, fractional spin excitations should be involved. This is furthermore supported by a fractional $n$ in Eq. (1), implying that fractions of a spin-flip excitation are involved in the relaxation process. In contrast to this unusual gapped $T_{1}^{-1}(T)$ behavior, the temperature dependence of the local susceptibility monitored by the 35Cl NMR shift in Fig. 2(b) is monotonic over the whole covered temperature range, as predicted for the ferromagnetic Kitaev model Yoshitake_2016 . Such a dichotomy between the two observables is a direct sign of spin fractionalization, as different fractional quasiparticles enter the two observables in different ways Yoshitake_2016 .

To obtain the spin-excitation gap $\Delta$ as a function of $B_{ab}$ in Fig. 1(c), the $T_{1}^{-1}(T)$ datasets in Fig. 3 taken in magnetic fields of different directions and magnitudes are fitted to Eq. (1) in the temperature range of the Kitaev paramagnetic phase. As the curve $T_{1}^{-1}\propto T^{-1}$ defined by $\Delta=0$ is steeper than any dataset in this range, the obtained excitation gaps are apparently all nonzero. The inset of Fig. 1(c) showing the symmetric $\Delta(\vartheta)$ dependence around $90^{\circ}$ in $9.4$ T, where $\vartheta$ traverses nonequivalent directions with respect to the Kitaev axes on both sides [inset of Fig. 1(a)], demonstrates that the underlying physics is indeed isotropic as assumed when introducing $B_{ab}$. The obtained $\Delta(B_{ab})$ in Fig. 1(c) can be perfectly reproduced as a sum of two terms: the two-flux gap $\Delta_{0}$ and the gap acquired by Majorana fermions in a weak magnetic field, predicted to be proportional to the cube of the field Kitaev_2006 ; Jiang_2011 ; Nasu_2017 ; Supplemental ,

[TABLE]

using $J_{K}=190$ K Do_2017 to evaluate $\Delta_{0}$, as before, while $\widetilde{B}=g_{ab}\mu_{B}B_{ab}/k_{B}$ is the field in kelvin units, $k_{B}$ is the Boltzmann constant, $\mu_{B}$ the Bohr magneton, and $\alpha=4.5$ (leading to the best fit) accounts for the sum over the excited states in the third-order perturbation theory, which is the origin of the $\widetilde{B}^{3}$ term Kitaev_2006 . This result demonstrates that a spin-flip excitation in $\alpha$-RuCl3 indeed fractionalizes into a gauge-flux pair and a Majorana fermion.

Focusing on the Kitaev paramagnetic region in the phase diagram of $\alpha$-RuCl3 in Fig. 1(b) is essential for our observation of two types of anyons. Instead, other recent experimental studies focused on the low-temperature region above $B_{c}$, observing the spin-excitation continuum Banerjee_2017_2 with either a gapless behavior Leahy_2017 ; Zheng_2017 or the gap opening linearly with $B-B_{c}$ Baek_2017 ; Sears_2017 ; Hentrich_2017 ; Ponomaryov_2017 , but without a definite conclusion about the identity of the involved quasiparticles. Such an ambiguous behavior likely originates in the presence of additional, smaller non-Kitaev interactions between the spins Banerjee_2016 ; Yadav_2016 ; Winter_2016 ; Ran_2017 , whose role should be pronounced particularly at low temperatures and which are indeed predicted to drive the system towards a gapless QSL ground state Song_2016 . Our result shows that spin fractionalization into two types of anyons is robust against these interactions in a broad range of temperatures and magnetic fields. This is the main practical advantage of $\alpha$-RuCl3 with respect to all other anyon realizations, such as the fractional quantum Hall effect in 2D heterostructures dePicciotto_1997 or hybrid nanowire devices Mourik_2012 , where anyons are observed only at extremely low temperatures and at certain field values. Our discovery thus establishes $\alpha$-RuCl3 as a unique platform for future investigations of anyons and braiding operations on them, which form the functional basis of a topological quantum computer Kitaev_2006 .

I Supplemental material

I.1 Crystal growth

Crystals of $\alpha$-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 $\alpha$-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 $\alpha$-RuCl3 was sublimed for the second time in order to obtain bigger crystal plates. The phase and purity of the compounds was 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.

I.2 Magnetic susceptibility

Magnetic susceptibility measurements of $\alpha$-RuCl3 were performed using a Quantum Design MPMS. A powdered sample of the mass $22.5$ mg was placed into a plastic capsule, in a glovebox to avoid contact with air, and then quickly transferred into the MPMS. Fig. 4 shows the measured susceptibility taken with cooling in field and in zero field. The obtained curve with the magnetic transition at $T_{N2}=14$ K (inset of Fig. 4) is almost identical to the corresponding curve in Ref. Johnson_2015 .

I.3 Nuclear magnetic resonance

General. The 35Cl nuclear magnetic resonance (NMR) experiments were performed on a foil-like $\alpha$-RuCl3 single crystal of approximate dimensions $5\times 5\times 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 fitting the sample size 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, in order 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. In order 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, the typical $\pi/2$ pulse length amounted to $2~{}\mu$s. The NMR signals were recorded using the standard spin-echo, $\pi/2-\tau_{d}-\pi$ pulse sequence with a typical delay of $\tau_{d}=70~{}\mu$s (much shorter than the spin-spin relaxation time $T_{2}$) between the $\pi/2$ and $\pi$ pulses.

$T_{1}$** relaxation.** The spin-lattice relaxation (i.e., $T_{1}$) experiment was carried out using an inversion recovery pulse sequence, $\varphi_{i}-\tau-\pi/2-\tau_{d}-\pi$, with an inversion pulse $\varphi_{i}<\pi$ (suitable for broad NMR lines) and a variable delay $\tau$ before the read-out spin-echo sequence. The spin-lattice relaxation datasets were typically taken at $20$ increasing values of $\tau$. The datasets were analyzed using the model of magnetic relaxation for $I=3/2$ spin monitored on the central $-1/2\longleftrightarrow 1/2$ transition:

[TABLE]

where $T_{1}$ 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. 2(b)], this expression reproduces the experimental relaxation curves perfectly. In the magnetically ordered phase, two $T_{1}$ components appear, and the relaxation curves are reproduced as a sum of two terms of the form given by Eq. (3). For instance, the temperature dependence of the corresponding two $T_{1}$ values for $B=9.4$ T with ${\bf B}\parallel c^{*}$ is given in Fig. 2(a). In cases where only a narrow temperature region below the transition was covered, the two components in the relaxation curves were hard to identify, and we used Eq. (3) with a stretched exponent instead.

Orientation dependence of the NMR spectrum. The 35Cl NMR spectra were recorded point by point in frequency steps of $50$ or $100$ kHz, so that the Fourier transform of the signal was integrated at each step to arrive at the individual spectral point. The covered NMR frequency range was from $34$ MHz, the lower limit of our setup, up to $50$ MHz. The dependence of the corresponding part of the NMR spectrum on the direction of the external magnetic field of $9.4$ T (described by the angle $\vartheta$ from the crystal $ab$ plane) at a temperature of $20$ K is shown in Fig. 5. The spectra are extremely broad because of large 35Cl (with $I=3/2$ spin) quadrupole interaction. As concluded in the following, a large portion of the covered frequency range is associated with the central, $1/2\longleftrightarrow-1/2$ 35Cl NMR transition. As this transition is observed to consist of at least three peaks (Fig. 5), even for the symmetric ${\bf B}\parallel c^{*}$ orientation (with $\vartheta=90^{\circ}$), while there are only two inequivalent Cl sites in the crystal structure, the splitting of the central line is likely a consequence of stacking faults in the layered crystal structure or crystal twinning, or both.

Relation bewteen orientation and field dependence of $T_{1}$. Measuring the $T_{1}$ dependence on the direction of the magnetic field (described by the angle $\vartheta$ from the crystal $ab$ plane) instead of on its magnitude in the $ab$ plane allows us to cover low $B_{ab}$ 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 (Fig. 5): to minimize an already large NMR linewidth 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 $\Delta(B_{ab})$ data points for various angles $\vartheta$ and field values $2.35$, $4.7$ and $9.4$ T in Fig. 1(c) all collapse on a smooth experimental curve. The $\Delta(B_{ab})$ data points taken in lower fields apparently exhibit much larger error bars. Namely, the corresponding $T_{1}^{-1}(T)$ datasets in Fig. 3(b) are more scattered than the datasets taken in $9.4$ T despite a much longer averaging for noise reduction.

Temperature dependence of the NMR line. We measured the temperature dependence of the dominant 35Cl NMR peak in a field of $9.4$ T for various sample orientations. From these measurements, we determined the temperature dependence of the frequency width [Fig. 6(a)] and the NMR shift of the peak with respect to the Larmor frequency [Fig. 6(b)]. For $B_{ab}<8$ T, the width exhibits a clear kink as a function of temperature, 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 $B_{ab}$ in the inset of Fig. 6, we obtain the phase boundary of the magnetically ordered state, which perfectly matches the result of the reference study Johnson_2015 . In contrast, the NMR shift does not exhibit any signs of a magnetic transition, except for the $\vartheta=90^{\circ}$ (${\bf B}\parallel c^{*}$) dataset. We find the NMR shift to be a monotonic function of temperature $T$, empirically following a ${\rm log}\,T$ dependence over a broad temperature range.

Contributions to the NMR shift. To separate the magnetic contribution to the NMR frequency shift from the temperature-independent quadrupole contribution, we plot the relative NMR shift (i.e., the NMR shift divided by the 35Cl Larmor frequency $\nu_{L}=39.18$ MHz in a field of $9.4$ T) measured on the dominant NMR peak for ${\bf B}\perp c^{*}$, i.e., in the $ab$ plane [the $\vartheta=0^{\circ}$ dataset in Fig. 6(b)], against the rescaled magnetic susceptibility $\chi_{ab}$ in the inset of Fig. 7. In Ref. Johnson_2015 , an experimental ratio between the susceptibility $\chi$ of the powdered sample and the susceptibility $\chi_{ab}$ of the single crystal with a field applied in the $ab$ plane is obtained as $(2+r)/3$ with $r=0.157$, leading to $\chi_{ab}=3\chi/(2+r)$. We use this empirical relation to evaluate $\chi_{ab}(T)$ from our field-cooled $\chi(T)$ dataset shown in Fig. 4. 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 susceptibility Kubota_2015 . From the observed linear relation between the relative shift and $\chi_{ab}$ up to $20\cdot 10^{-3}$ emu/mol (i.e., down to $35$ K), we obtain the hyperfine coupling constant $A=2.2$ T/$\mu_{B}$ and the zero-temperature relative shift $-0.039$ that, when multiplied by $\nu_{L}$, gives the quadrupole shift $\Delta\nu_{Q}=-1.53$ MHz.

From the obtained quadrupole shift $\Delta\nu_{Q}$, we can estimate the quadrupole splitting $\nu_{Q}$ between the successive 35Cl NMR transitions. For the case of an axially symmetric EFG tensor and the field applied at an angle $\vartheta^{\prime}$ from the principal EFG axis $v_{ZZ}$ with the largest EFG eigenvalue, the second-order quadrupole shift is given by $\Delta\nu_{Q}=-3\nu_{Q}^{2}(1-\cos\vartheta^{\prime 2})(9\cos\vartheta^{\prime 2}-1)/(16\nu_{L})$ for the $I=3/2$ nucleus Abragam_2011 . As the axes of the EFG tensor are not known, we assume a typical tilt $45^{\circ}$ of $v_{ZZ}$ from $c^{*}$, so that $\vartheta^{\prime}\sim 45^{\circ}$. From the previously evaluated $\Delta\nu_{Q}$ we then obtain $\nu_{Q}\sim 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 of Fig. 5 all belong to the central transition.

I.4 Theory

Theoretical $T_{1}^{-1}(T)$ curve. The theoretical temperature dependence of $T_{1}^{-1}$ is numerically calculated for the Kitaev model in zero field Yoshitake_2016 . $T_{1}^{-1}$ contains two contributions, one coming from a single fluctuating spin (i.e., on-site) and the other one coming from fluctuating nearest-neighboring (NN) spins in the Kitaev honeycomb lattice. As the 35Cl nucleus in $\alpha$-RuCl3 is located at equal distances from the closest two Ru3+ $S=1/2$ spins [Fig. 1(a)], $T_{1}^{-1}$ contains both contributions with equal weights. Namely, as the hyperfine coupling constant $A$ of 35Cl to both spins is the same, the relevant spin-spin correlation function can be generally written as

[TABLE]

for the involved components $S_{1}$ and $S_{2}$ of both Ru3+ spins, where the plus (minus) sign is valid for ferromagnetic (antiferromagnetic) fluctuations. The first two terms on the right side of Eq. (I.4) represent the on-site contributions, while the last two represent the NN-sites contributions, both with apparently equal weights. Accordingly, the theoretical curve for the ferromagnetic case, plotted in Fig. 2(a), is the average of the on-site and NN-sites contributions.

$T_{1}$** 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 by Beeman_1968

[TABLE]

where $E$ is the energy of magnons, $g(E)$ is their density of states, $n(E)=[\exp(\beta E)-1]^{-1}$ is the Bose-Einstein distribution function, $\beta=1/(k_{B}T)$, and $k_{B}$ is the Boltzmann constant. Denoting the magnon gap by $\Delta$ (in kelvin units), we define $\varepsilon=E-k_{B}\Delta$ as the energy measured from the bottom of the magnon band. The power-law dispersion relation $\varepsilon\propto k^{s}$ in $D$ dimensions, which includes the standard parabolic dispersion ($s=2$) and the Dirac dispersion ($s=1$) as special cases, leads to $g(E)\propto\varepsilon^{D/s-1}$. For low temperatures $T\ll\Delta$, the distribution function $n(E)$ can be approximated by the Boltzmann distribution, $n(E)\approx\exp(-\beta E)=\exp(-\Delta/T)\exp(-\beta\varepsilon)$. Plugging these expressions for $g(E)$ and $n(E)$ into Eq. (5), we obtain

[TABLE]

The integral on the right side of Eq. (6) converges if $s<2D$ and evaluates to $\Gamma(2D/s-1)$ where $\Gamma$ is the gamma function. We can thus rewrite Eq. (6) as

[TABLE]

with the power of the prefactor $p=2D/s-1$. In case of $D=2$, which is relevant for the Kitaev honeycomb magnet, $p=1$ for $s=2$ and $p=3$ for $s=1$, so that the power $p$ cannot be negative. Even in case of $D=1$, $p$ can only reach the lowest value of [math] precisely for $s=2$ [although care should be taken in this case, as the integral in Eq. (6) then formally diverges]. If more than a single magnon is involved in the $T_{1}$ process, the power $p$ is also positive and becomes even higher Beeman_1968 . Gapped magnons thus cannot lead to the $T_{1}$ relaxation described by Eq. (7) with $p<0$.

Instead, we can use Eq. (7) in the 3D magnetically ordered state, when the elementary excitations are magnons with a gap $\Delta_{m}$. In this case $D=3$ and $s=2$, and this leads to $T_{1}^{-1}\propto T^{2}\exp(-\Delta_{m}/T)$. We use this expression to analyze the $T_{1}^{-1}(T)$ data [Fig. 2(a)] in the low-temperature ordered state of $\alpha$-RuCl3.

As a side observation, all these examples show that a frequently used simple gapped model $T_{1}^{-1}\propto\exp(-\Delta_{s}/T)$ with the gap $\Delta_{s}$, which was used before to analyze the $T_{1}^{-1}(T)$ datasets in $\alpha$-RuCl3 Baek_2017 , is actually not justified in any region of the phase diagram of $\alpha$-RuCl3.

Majorana fermion gap. In the Kitaev model, Majorana fermions acquire a gap in the presence of an external magnetic field Kitaev_2006 . This is shown for a field applied perpendicularly to the honeycomb plane, i.e., in the $(1,1,1)$ direction in the coordinate system defined by the Kitaev axes $x$, $y$ and $z$. The corresponding Zeeman term then reads $\mathcal{H}_{Z}=-h\sum_{j}(S_{j}^{x}+S_{j}^{y}+S_{j}^{z})$, where $h=g\mu_{B}B/\sqrt{3}$ is a single component of the magnetic field $B$ in energy units, $g$ is the $g$-factor and $\mu_{B}$ is the Bohr magneton. When treated as a perturbation to the Kitaev Hamiltonian, the Zeeman term contributes to the Majorana fermion gap only at third order Kitaev_2006 . The corresponding effective Hamiltonian is thus proportional to $h^{3}$ and can be written as Kitaev_2006 ; Jiang_2011 ; Nasu_2017

[TABLE]

where $\Delta_{0}$ is a two-flux gap (in kelvin units), while $\alpha$ (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 $-\kappa\sum_{jkl}S_{j}^{x}S_{k}^{y}S_{l}^{z}$ with $\kappa=\alpha h^{3}/(k_{B}^{2}\Delta_{0}^{2})$ is still exactly solvable and the dispersion relation of the Majorana fermions is calculated as Jiang_2011

[TABLE]

where $J_{K}$ is the Kitaev coupling (in kelvin units), while ${\bf a}_{1}$ and ${\bf a}_{2}$ are the unit vectors of the honeycomb lattice. The dispersion relation given by Eq. (9) is gapped for $\kappa\neq 0$, and the corresponding gap $\Delta_{f}$ can be calculated numerically as a function of $\kappa$ and thus as a function of the magnetic field. For small magnetic fields, i.e., for $\kappa\ll k_{B}J_{K}$, the Majorana fermion gap (in kelvin units) simplifies to

[TABLE]

while for high magnetic fields it saturates to $\Delta_{f}=2J_{K}$. The total spin-excitation gap $\Delta$ is obtained by adding $\Delta_{f}$ to the two-flux gap $\Delta_{0}$. The field dependence of $\Delta$ is shown in Fig. 8 for $J_{K}=190$ K (taken from Ref. Do_2017 and used in this work), $g=g_{ab}$ and $\alpha=4.5$ [leading to the best fit of our $\Delta(B_{ab})$ data points]. The cubic approximation given by Eq. (10), which is also plotted in Fig. 8, is apparently valid up to $15$ T, well beyond the field range covered in this work.

I.5 Comparison with recent works

Recent NMR works. Very recently, two 35Cl NMR studies of $\alpha$-RuCl3 appeared Baek_2017 ; Zheng_2017 . The analysis of both studies is focused on the low-temperature region below $15$ K. Using a simple exponential model $T_{1}^{-1}\propto\exp(-\Delta_{s}/T)$ in this region, Ref. Baek_2017 finds that the excitation gap $\Delta_{s}$ opens linearly with the field above the critical field around $10$ T. On the other hand, Ref. Zheng_2017 extends the covered temperature range down to $1.5$ K and finds a low-temperature gapless, power-law behavior of $T_{1}^{-1}(T)$ in the covered high-field region above the critical field around $8$ T.

As these results are very different from our results, also because of quite different analysis, we analyze the $T_{1}^{-1}(T)$ datasets obtained in these two works also with our model given by Eq. (1). As in our work, we focus on the Kitaev paramagnetic region, to the temperature range from $50$ K down to slightly above the transition temperature below $8$ T, and down to $4.2$ K above $8$ T (or a bit higher at higher fields), including the characteristic maximum in $T_{1}^{-1}(T)$ as a main feature. The data in Ref. Zheng_2017 were taken with a field applied in the crystal $ab$ plane, while the data in Ref. Baek_2017 were taken with a field applied at $30^{\circ}$ and $-60^{\circ}$ with respect to the $ab$ plane. In this case, we calculate the effective field values $B_{ab}$ in the same way as in our work. The obtained field dependence of the excitation gap $\Delta(B_{ab})$ for both works is shown in Fig. 8 together with our result. Applying our analysis to the data in all three works apparently leads to relatively consistent results. Nevertheless, the results for the data from Refs. Baek_2017 ; Zheng_2017 alone do not allow to conclude on the cubic field dependence of $\Delta$, mostly due to the lack of important low-field data points. Regarding Ref. Baek_2017 , the two data points at the highest fields seem to deviate from the trend set by the other points. The corresponding two $T_{1}^{-1}(T)$ datasets exhibit a suspicious plateau at low temperatures, not observed in any other dataset of the three works, which casts some doubts on their credibility. Regarding Ref. Zheng_2017 , the obtained $\Delta(B_{ab})$ data points, which complement our field region, nicely continue the trend set by our data points. The obtained trend is approximately linear in field, consistent with the theoretical prediction in this intermediate-field region, although with a different slope. However, the theoretical prediction is based on a perturbative treatment, which may not give reliable results outside the low-field region.

In light of our conclusions, the reported findings of these two works should be understood in the following way. While the true signs of fractionalization into Majorana fermions and gauge fluxes can be found in the Kitaev paramagnetic region covering a broad temperature range up to around $100$ K [Fig. 1(a)], the physics at low temperatures of the order of $10$ K and below is apparently obscured, most likely due to the effect of inevitable non-Kitaev interactions, as predicted already in Ref. Zheng_2017 .