The proposal that a quasi-2D Mott insulator α-RuCl 3 may provide a realization [1] of Kitaev's celebrated honeycomb compass model [2] has attracted much attention to this material. While α-RuCl 3 orders antiferromagnetically below T N = 7 K [3], it was found that an in-plane magnetic field h ≡ µ 0 H 10 T is sufficient to suppress the magnetic order. While the nature of the resulting phase is still under intense debate, the observation of approximately quantized value of the thermal Hall conductivity κ xy /T in a narrow range of field (6 < h < 9 T) [4,5] was attributed to the presence of the Majorana edge mode, predicted to exist in Kitaev's spin liquid subjected to an external magnetic field [2,6]. This interpretation has been recently challenged by an independent measurement of the thermal Hall effect [7], in which the authors find a non-quantized, temperature-dependent κ xy , which they attribute to a bosonic, rather than fermionic mechanism [8][9][10][11][12]. Its nature remains controversial, with one recent experimental study suggesting the possibility of quantized Hall effect in high fields h > 10 T [13], while another attributing the origin of the thermal Hall effect to phonons [14].

In this Letter, we investigate the possibility of the bosonic origin of thermal conductivity in a widely accepted spin model of α-RuCl 3 . Since there is a considerable debate on the precise values of the model parameters describing α-RuCl 3 , we perform a careful scan over a wide region in the parameter space to determine the largest possible values of κ xy /T . We find that the bulk magnon excitations alone cannot explain the experimentally measured values of thermal conductivity in α-RuCl 3 , even under the most favourable circumstances. Instead, we find that it is crucial to take the magneto-elastic coupling into consideration, whereby acoustic phonons hybridize with the magnon excitations, boosting the value of κ xy /T . Even then, it turns out that in order to explain the experimental measurements, one must consider not only intrinsic but also extrinsic contributions to the thermal Hall effect, such as the skew-scattering of phonons/magnons off of impurities. We deduce the realistic value of this extrinsic contribution from a recent measurement on α-RuCl 3 .

When magnon, phonon and extrinsic contributions are taken into account, we are able to quantitatively reproduce the recent experimental data [7] on thermal conductivity in this material.

Model and Phases. α-RuCl 3 emerged as a candidate material to study Kitaev physics on the honeycomb lattice because of its purported proximity to the spin-liquid state [1]. In addition to the Kitaev's bond-dependent interactions stemming from the interplay of spin-orbit coupling and superexchange between Ru 3+ ions [15], the importance of nearest-neighbor Heisenberg interactions J 1 and the off-diagonal exchanges, so-called Γ and Γ terms [16] has been established. Much theoretical work [9,[17][18][19][20][21][22][23][24][25][26][27] has since focused on deducing the values of these parameters in α-RuCl 3 , leading to the minimal effective spin- 1 2 model of the form [10,18,19,26,28,29]

where the third-neighbour Heisenberg exchange J 3 [22] was also added. The index α = (x, y, z) enumerates the three nearest bonds on the honeycomb lattice and also labels the bond-dependent spin couplings, with the remaining indices β, γ taking values among the cyclic permutations of (x, y, z) indices, for a given α (see SM).

Since the experiments are conducted under the applied magnetic field along the a-axis, its effect is captured by the last term in Eq. (1) with the Landé g-factor g = 2.5 [21,24,30]. It is important to emphasize that the Kitaev axes (x, y, z) are the so-called cubic axes [22] that do not coincide with the crystallographic ones. In particular, the magnetic field along the a-axis has nonzero components along all three Cartesian x, y, z axes, which in the pure Kitaev model is predicted to open a spectral gap proportional to the third power of the field [2].

In all the parameter sets proposed in previous works, obtained either from first-principles calculations or from phenomenological analysis (see e.g. Ref. [27] for review), the leading coupling is believed to be the ferromagnetic Kitaev term K < 0 [31], with the off-diagonal term Γ > 0 large and potentially comparable to |K|. In what follows, we assume the strength of the Kitaev interaction to be K = -7.2 meV, which is close to a recent ab initio derived value of 80 K and is in the middle of the "realistic parameter regime" proposed in Ref. [27]. The subleading Heisenberg interactions J 1 < 0 and J 3 > 0 are also necessary to explain the ordered zigzag phase of α-RuCl 3 . The behaviour of this model in the applied field is illustrated in Fig. 1 for a representative parameter set (K, J 1 , J 3 , Γ, Γ ) = (-7.2, 0, 0.8, 0, -0.2) meV. As the strength of the magnetic field (along the a-axis) increases, the spins tilt along the field direction, resulting in the canted zigzag phases ZZ1 and ZZ2 depicted schematically in Figs. 1(b,c) -what distinguishes these two phases is the plane in which the spins of the two sublattices lie. At a sufficiently large field (whose value depends on the model parameters, and here h sat = 19 T), a fully polarized (PL) phase is reached. We use the standard linear spin wave theory (LSWT) (see Supplementary Materials (SM)) to compute the magnon spectrum of the model, and hence the thermal conductivity, given by the well known formula [32][33][34][35][36]:

where

is the dilogarithm function, f ( ) is the Bose-Einstein distribution, and the summation is over all the magnon bands. To compare with the prediction of the two-dimensional Kitaev QSL originating from the Majorana edge modes: κ 2D xy /T = πk 2 B /12 , we compute the same quantity in the unit of fermionic quantized value πk 2 B /6 :

where d = 5.72 Å is the interlayer distance of α-RuCl 3 . In these fermionic units (f.u.), the quantized value reported in Ref. [4] would be κ Maj xy = 0.5 f.u. We compute the integral by summing over the Berry flux φ n k in each small plaquette of the Brillouin zone (BZ) [37,38], weighted by the c 2 function (see SM for more details).

xy |/T is plotted, at T = 10 K, as a function of increasing magnetic field in Fig. 1(f), using the same model parameters used to show the different phases in panel (e). Note that |κ xy |/T first increases monotonically in the ZZ1 phase, stable below field h c1 , before changing sign in the ZZ2 phase. In the region h c2 < h < h c3 , the computed magnon band structure becomes unphysical, meaning the failure of the zigzag ansatz to capture the true ground state, which may be a different four-spin order [10], or a magnetic order with an enlarged unit cell [39] that is beyond the scope of the present study -we use UN to represent this unknown phase. Finally, the system enters fully polarized phase for fields h > h c3 , where |κ xy |/T decreases with increasing field. As this plot illustrates, the intensity of the thermal Hall effect thus obtained is always smaller than at most 0.2 in the fermionic units -or about 40% of the value observed in Ref. [4]. Below, we explore what the upper bound is on the thermal conductivity as a function of model parameters.

Upper bound on the thermal Hall effect due to magnons. Since the precise values of the model parameters in Eq. ( 1) corresponding to α-RuCl 3 are still under intense debate (see Table 1 in Ref. [27] for a list of different proposals), we scan a wide range of the physically relevant parameter values with the goal of determining the upper bound on κ 2D xy /T . The parameter ranges we used are 0 < Γ < 7.2 meV, -4.2 < J 1 < 0 meV, 0 < J 3 < 4.0 meV, and -3.6 < Γ < 3.6 meV with step size δ = 0.2 meV (while keeping the magnitude of the Kitaev term K = -7.2 meV fixed as stated earlier).

According to the recent experimental data in Ref. [7], the thermal Hall conductivity tends to be largest in the high field range h ≈ 10 T and at moderately high temperatures 8 T 12 K. Therefore, for concreteness, we investigate the magnitude of κ 2D xy /T under the relevant experimental conditions h = 10 T and T = 10 K. We first start with the fully polarized (PL) phase. Because of the difficulty of representing the plots in the four-dimensional parameter space (Γ, Γ , J 1 , J 3 ), we choose to plot the distribution of the κ 2D xy /T values vs. J 1 in Fig. 2(a), with each data point corresponding to a different choice o the remaining parameters. It is clear from this panel (a) that the largest values of κ 2D xy /T are attained at negative and very small J 1 . This is because a weak ferromagnetic J 1 will destabilize the polarized state, and result in the lower magnon band moving down, which leads to stronger magnon modes contribution to the thermal Hall effect because of the increasing weight of the c 2 (x) function in Eq. ( 2). In Ref. [27] the authors identified the linear combination J 1 + 3J 3 as the relevant parameter, which is used as the x-axis to plot the calculated κ 2D xy /T in Fig. 2(b), with similar conclusions reached.

An alternative way of looking at the data is to plot the maximum value of κ 2D xy /T as a false color on a twodimensional plot with axes given by J 1 and Γ, which we do in Fig. 2(c), or following the strategy proposed in Ref. [27], with the axes formed by effective couplings J = J 1 + 3J 3 and Γ = Γ + 2Γ , shown in Fig. 2(d). Each data point in these panels is taken to be the maximal value of |κ 2D

xy |/T from varying the remaining parameters. We find that the largest |κ 2D

xy |/T in the PL phase (at h = 10 T and T = 10 K) never exceeds about 0.35 f.u. in the fermionic units. To orient the reader, we show with the dashed rectangle what the authors of Ref. [27] call the "realistic parameter regime," and the cross represents the parameter set chosen in Ref. [11]. In both regions, we find |κ 2D xy |/T to be less than 0.2 f.u., far from 0.5 f.u. reported in Ref. [4] and much below the maximum value measured in Ref. [7].

Finally, we select the parameter set (labeled ps1) with the largest value of thermal Hall effect in our data and plot its value as a function of temperature and field, shown in Figs. 2(e,f). While the monotonically increasing temperature dependence observed in Ref. [7] is qualitatively reproduced, the field dependence is opposite -the experiment shows an increasing |κ xy (h)|, while our data invariably decrease monotonically. Its physical reason is clear -in the fully polarized phase, the increase of the magnetic field leads to the (linear in field) growth of the magnon gap min . And since the Berry curvature integrand in Eq. ( 2) is weighted by the function c 2 (f ( )) ∝ exp(-/T ), its value is exponentially suppressed at the experimentally relevant temperatures, leading to the decrease in |κ xy (h)|. We thus conclude that not only is the magnon contribution too low to account for the experimental value of thermal Hall effect, but its field dependence in the fully polarized phase cannot reproduce the experiment neither.

We therefore turn our attention to the zigzag phases, with the distributions of κ 2D xy /T along the J 1 and J = J 1 + 3J 3 shown in Figs. 3 (a xy /T is plotted as a false color map in (J 1 , Γ) and ( J, Γ) coordinates, respectively, in Figs. 3(c) and(d). We find the largest value of κ 2D xy /T to be about 0.2 f.u. When considering the "realistic parameter regime" in Ref.

[27] (dashed region) or the parameters used in Ref. [11] (cross), κ 2D xy /T does not exceed 0.15 f.u., far below the experimentally reported values [4,5,7].

To compare with the experimental data, we choose two sets of parameters (labeled ps2 and ps3 in Fig. 3, the values are listed in the caption) that belong to the ZZ1 and ZZ2 phases, respectively, and which have the largest values of |κ 2D

xy |/T in our studied range. We evaluate the temperature and field dependence of |κ 2D

xy |/T at these parameter sets and plot them against the experimental data from Ref. [7] in Figs. 3(e,f). Under the experimentally relevant conditions 2 < T < 6 K and 6 < h < 7 T, we find that the parameter set ps2 qualitatively matches the trends in the temperature and field dependence of the experimental thermal Hall data. Near T = 10 K and h = 10 T, κ xy /T continues to increase for ps2 (as is the case experimentally); whereas for ps3 it starts to decline. However, as noted already and as seen from Figs. 3(e,f), the computed magnitude of |κ 2D

xy |/T is well below the experimental values. This indicates that the intrinsic magnon contribution alone cannot fully account for the measured thermal Hall conductivity in α-RuCl 3 . We thus turn our attention to additional, bosonic in nature contributions to the thermal Hall effect.

Phonon contribution to the thermal Hall effect. We now turn to investigate other effects that may enhance the thermal Hall effect. Recent experimental data [14] on α-RuCl 3 show that the temperature dependence of κ xy resembles closely that of longitudinal thermal conductivity κ xx , as demonstrated in Fig. 4(a), with the ratio between the two roughly the same (0.03 -0.10%) across different samples, from which the authors of Ref. [14] conclude that phonons must play a key role in the thermal Hall effect. In an insulator, the logitudinal κ xx is dominated by non-chiral acoustic phonons. The thermal Hall effect by contrast is time-reversal odd and chiral in nature. We distinguish two mechanisms of such chiral phonon response: the intrinsic one, due to the Berry curvature induced by magnon-phonon coupling, and the extrinsic one, due to phonon scattering off of defects.

Intrinsic phonon contribution to κ xy . The distance dependence of the superexchange interactions between Ru 3+ ions leads naturally to the magnetoelastic (ME) coupling of the generic form (see SM for further details):

where u i is the displacement of the ion at site i from its equilibrium position. Writing these displacements in terms of the phonon operators u γ i ∼ (a † iγ + a iγ ) (with polarization γ), this results in the hybridization between the magnons and the phonon branches [40][41][42], endowing the phonons with the Berry curvature that contributes to κ xy just as in Eq. ( 2). The magnitude of the additional phonon contribution κ

xy is related to the strength of the ME coupling, whose existence is supported by the softening of the phonon branch at small q-vectors meansured in experiments [43]. Choosing the ME coupling that qualitatively reproduces the phonon softening (see SM for de-tails), the value of κ xy /T is plotted in Fig. 4(b) with the solid blue line using the parameter set ps2 in the ZZ1 phase. As this figure shows, the magnitude of the intrinsic κ I xy /T is enhanced compared to the magnon-only value in the temperature regime of interest T > 3 K -the same conclusion was reached in the recent study Ref. [12]. However, the magnitude κ I xy /T ends up still being smaller than the experimentally measured Hall signal, prompting us to consider additional, phonon contributions. Extrinsic phonon contribution to κ xy . There are multiple sources of phonon contributions to the thermal Hall effect, in analogy to the phonon contribution to the electronic Hall effect (σ xy ) in metals [44]. One candidate mechanism is the intrinsic skew-scattering, which originates from the Lorentz force on ions [45]. Another source is the extrinsic skew-scattering from phonons scattering off of magnetic impurities [45,46]. However, after comparing the experimental data [14] with characteristic features of these effects (see SM for detailed analysis), we came to the conclusion that they are negligible in α-RuCl 3 .

Instead, the very weak temperature dependence of the ratio κ xy /κ xx (see Fig. 4a) and absence of strong sample variability in recent experiments [14] indicate that the dominant extrinsic contribution to κ xy is most likely from the so-called side-jump scattering of phonons off of defects [44], which we demonstrate in the SM using the formalism recently developed in Ref. [47]. Crucially, this effect scales with the phonon mean-free path , just like the longitudinal thermal conductivity, consistent with the ratio of the two being sample independent.

Using the experimental data from Ref. [14] at high temperatures T |K|/k B ≈ 80K, above the magnon bandwidth where the effects of the Kitaev physics and associated Berry curvature are unimportant, we determine the Hall angle due to extrinsic scattering θ E H ≡ κ E xy /κ xx = (0.6 ± 0.2) × 10 -3 . By contrast, at low temperatures T ∼ 10 -15 K where the interpretation of the thermal Hall measurements on α-RuCl 3 is disputed, both the intrinsic (due to the Berry curvature) and extrinsic contributions to κ xy must be taken into account. By comparing the data in this low-T region with the value θ E H from above, we are able to determine the phenomenological ratio η ≡ κ E xy /κ I xy . We obtain (see SM for details) η = 1.2±0.5 with the uncertainty related to the spread of the experimental data among the different samples. Taking the intrinsic κ I xy (solid line in Fig. 4(b)) and multiplying by (1 + η), we thus obtain the total κ xy = κ I xy + κ E xy , marked by the shaded blue region in Fig. 4(b). The experimental data from Ong's group [7] fall inside the yellow shaded region. This indicates that the thermal Hall effect from the bosonic model given our parameter choice can explain the experimental observation.

Furthermore, given the flexibility in the magnitude of the ME coupling, we show (see SM) that choosing the largest physically allowed coupling results in the intrinsic value of κ I xy /T of the order of 0.35 f.u. (for ps2), and further including the extrinsic effects can yield values of thermal Hall effect even in excess of the experimentally measured values. We thus conclude that the magnonphonon mechanism proposed in this work is more than sufficient to explain the experimental thermal Hall data in α-RuCl 3 without overly fine-tuning the model parameters, i.e., in a finite region in the parameter space.

Discussion. Having scanned a broad range of physically motivated parameters of the (K, J 1 , J 3 , Γ, Γ ) generalized Kitaev-Heisenberg model, we conclude that the intrinsic magnon contribution alone is insufficient to explain the large observed magnitude of the thermal Hall effect in α-RuCl 3 . We further found that in order to reconcile the observed magnetic field dependence of the thermal Hall effect, it is necessary to conclude that α-RuCl 3 remains in the canted zigzag phase (as opposed to field polarized) in fields up to ∼ 10 T. This conclusion is supported by the recent study by Li and Okamoto [12] who found that the spin-phonon coupling tends to stabilize the canted zigzag phase for higher applied fields, compared to the pure spin model which would otherwise become fully polarized. Taking into account the spinphonon coupling endows phonons with the chirality, contributing an additional intrinsic term to κ I xy , which however still falls short of explaining the experimental data, necessitating the inclusion of an extrinsic source of Hall effect. With minimal assumptions as to its mechanism, we used the existing experimental data to quantitatively arrive at the measure κ E xy /κ I xy yielding values between 1 and 2, meaning that the extrinsic phonon contribution to κ xy is comparable, or a little larger, than the intrinsic Berry curvature effect. Taking both into account, we are able to explain not only the large magnitude but also the detailed temperature dependence of κ xy , which is bosonic in nature.

Supplementary Materials for "Magnons, Phonons, and Thermal Hall Effect in Candidate Kitaev Magnet α-RuCl 3 "

The magnetic sector of the Hamiltonian is given in Eq. ( 1) in the main text, which we repeat here:

It contains five free parameters J 1 , K, Γ, Γ , J 3 in the spin-spin couplings, and also the external magnetic field h. We note that the components α = (x, y, z) correspond to the Kitaev ("cubic") axes, rather than the crystallographic (a, b, c) axes. The relation between the two coordinate systems is given by:

We note in passing that while additional terms in the spin Hamiltonian have been considered in the literature (see e.g. Ref. [48]), the above J 1 -K-Γ-Γ -J 3 parametrization appears to be widely adopted [10,18,19,[26][27][28][29].

Ordered phases and reference states. For each parameter set, we first determine the magnetically ordered ground state by minimizing the classical energy. Because the spins in α-RuCl 3 form a zigzag order at a low field and transition to a polarized phase at a high field, we make the ground state ansatz that the ground-state spin configurations are captured by the two-sublattice unit cell of the (canted) zigzag-type ordering, with four Ru atoms per unit cell (labelled A, B, C, D in Fig. 1 (a)). There are regions in the phase diagram where the classical calculations have shown that magnetic orders with larger unit cells may exist [11,39], however for the purpose of a manageable computation, we ignore such regions (which turn out to be very small in our parameter scans) and limit our considerations to where the ground states are captured by the four-site ansatz.

Under the application of an external magnetic field (taken to be along the a-axis as in the experiments), these four spins cant along the field direction, which can be parametrized by the polar and azimuthal angles (θ, φ) in the (abc) crystallographic coordinates, as shown in Fig. 1 in the main text, such that:

For fixed model parameters and magnetic field, the classical energy thus becomes a function of the four angles (θ 1 , φ 1 , θ 2 , φ 2 ). After the classical energy minimization, the reference phases that we observe, in the order of increasing magnetic field, are as follows:

• Zigzag phase I (ZZ1): canted zigzag phase with spins in the (ac)-plane, such that φ 1 = φ 2 = 0 and θ 1 = θ 2 ;

• Zigzag phase II (ZZ2): canted zigzag phase with spins in the (ab)-plane, such that θ 1 = θ 2 = π 2 and φ 1 = -φ 2 ; • Polarized phase (PL) with spins along the a field direction:

These spin configurations are shown in Figs. 1 (b,c,d) in the main text. In the small field region 0 < h < h c1 , the system is in the ZZ1 phase, in which the spins are situated in the ac plane, see Fig. 1(b). With increasing field h > h c1 , the system transitions to the ZZ2 phase, where the spins lie in the ab plane, see Fig. 1(c). In the region h c2 < h < h c3 , the magnon band structure becomes unphysical (having imaginary eigenvalues), symptomatic of the failure of the zigzag ansatz to capture the true ground state, which may an indication of other four-spin order, or of the need to consider enlarged magnetic unit cell, such as for instance found in the semiclassical analysis [11,39]. Analyzing such enlarged unit cells is beyond the scope of the present work and we use the abbreviation UN in Fig. 1 of the main text to represent the unknown phases. As mentioned in the previous paragraph, such unknown phases constitute only a very small region of the parameter regime we have surveyed in this work and we do not expect their existence to qualitatively alter our conclusions. Finally, the system enters the fully polarized phase for fields h > h c3 .

A side note is that we sometimes (depending on the model parameters) find a small region of a noncollinear zigzag phase III (ZZ3) where the spins are not confined to either the ac nor ab planes. Occasionally, we also observe a partially polarized (PPL) phase, in which all the spins are collinear but do not point along the magnetic field (this occurs due to strong spin-orbit coupling, for large |Γ| or |Γ |). Given the tiny regime the ZZ3 and PPL phases, and the fact that unlike the previously discussed phases, they do not appear universally in all parameter sets, we focus in what follows on the three main phases ZZ1, ZZ2 and PL.

Linear Spin Wave Theory. In each of the phases above, we perform the linear spin wave theory (LSWT) calculations to obtain the magnon band structure. The quantization axis (local z direction) is chosen such that it coincides with the mean-field spin direction on a given site. In the following discussion, the tilde over the spin operators ( Si ) indicates that these are in the local coordinate frame with the site-dependent quantization axis.

We then perform the standard Holstein-Primakoff transformation in this local basis, expressing the spin operators in terms of the magnon creation and annihilation operators a † and a:

Upon the Fourier transformation to k-space, the Hamiltonian in Eq. (S1) takes a quadratic form in the magnon operators:

where

Here N is the number of sites in the magnetic unit cell. N = 4 for ZZ1 and ZZ2 orders, and N = 2 for PL order. The matrix in the Nambu space is of the form

with H 11 ( k) and H 12 ( k) being N × N matrices. To be concrete, in the case of the PL phase, the matrix acts on the ket vector

) , and the block submatrices H 11 ( k) and H 12 ( k) are given by

where the k-dependent functions are

In the case of ZZ1 and ZZ2 phases,

) is the Nambu ket-vector composed of the magnon operators on sublattices A, B, C, D in Fig. 1(a). The band dispersions n k and eigenvectors |n, k at each k point are obtained by the similarity (Bogoliubov) transformation.

In order to compute the thermal Hall conductivity, one must compute the integral of the Berry curvature weighted by the c 2 (f ( n k )) function, as explained in Eq. ( 2) in the main text. For numerical purposes, the integral is replaced by a discrete sum (same as Eq. 3(3) in the main text) as follows:

Here φ n k is the Berry flux through a small plaquette formed by δ 1 = (Q a /N p , 0) and δ 2 = (0, Q b /N p ) with Q a = 2π/ √ 3 and Q b = 4π/3 (for N = 2) or 2π/3 (for N = 4), respectively. The elementary flux is given by [37,38]:

where Σ 2N = σ 3 ⊗ I N ×N with the Pauli matrix σ 3 = diag{1, -1} and I N ×N the identity matrix of size N . In our calculations of Eq. (S11), we use the discrete mesh N p × N p , and we show that N p = 71 is sufficiently large for the discretized integration to converge numerically in the pure magnon case (see Fig. S1). FIG. S1. Demonstration of convergence of computed κxy/T for model parameters (J1, J3, K, Γ, Γ ) = (0, 0, -7.2, 2.2, -0.2) meV and field h = 10 T. The k-space mesh Np × Np was chosen to perform numerical integration in Eq. (S11), with varying Np as shown in the legend. The results demonstrate that convergence is reached for Np ≥ 71. In the scan over the parameter space (Figs. 2,3), we used Np = 71.

To compare with the experimental data, we choose three sets of parameters where we find the largest values of |κ xy /T | in the scanned parameter range. These parameters sets (ps) belong to the PL, ZZ1 and ZZ2 phases, respectively, in the in-plane field of h = 10 T relevant for comparison with experiments:

These parameter sets are labeled by the colored diamond symbols in Figs. 2(c,d) and Figs. 3(c,d).

In addition, we compare our model parameters with those used by previous studies: our ps1 and ps2 lie close to the parameter set used by Zhang et al. [11] (labeled by a cross in Figs. 2 and3), while ps3 lies on the edge of what Maksimov and Chernyshev call the "realistic parameter regime" in their work [27], shown with a dashed rectangle in Figs. 2 and3. In all cases, we find that the maximum contributions of magnons to |κ 2D

xy |/T do not exceed about 0.3 fermionic units, significantly lower than the experimentally reported values of thermal Hall conductivity around 0.5 fermionic unit [4,7].

The temperature T and magnetic field h dependence of κ xy /T in these three parameter sets are plotted against the experimental data in Figs. 2 and3 in the main text. In the case of the polarized phase ps1, the field dependence of |κ xy |/T is qualitatively different from the experimental data, monotonically decreasing with the increasing magnetic field. The possibility of being in the polarized phase at 10 T is thus ruled out, as explained in the main text.

For the remaining two parameters sets (ps2 and ps3), under the experimentally relevant conditions 2 < T < 6 K and 6 < h < 7 T, we find that the parameter set ps2 qualitatively matches the trends in the experimental data of κ xy /T . Near T = 10 K and h = 10 T, κ xy /T continues to increase for ps2 (as is the case experimentally); whereas for ps3 it starts to decline slowly. The parameter set ps2 (in the ZZ1 phase) thus captures the experimental behavior of κ xy /T qualitatively in the temperature and magnetic regions of the measurement. However, we find the strength of κ xy /T in all cases to be much smaller than the experimental result. This mismatch shows that the intrinsic magnon contribution alone cannot fully account for the experimentally measured thermal Hall conductivity in α-RuCl 3 . This corroborates the experimental suggestions [14] that other sources, in particular phonons, must contribute significantly to the thermal Hall effect in this material.

Given the suggested importance of acoustic phonons in contributing to the thermal Hall in α-RuCl 3 [14], in this section we describe the modeling of phonons and phonon-magnon coupling, in order to explore their intrinsic contributions to the thermal Hall conductivity. The extrinsic (scatterer-dependent) contribution will be discussed separately in section V.

To describe the phonon modes in α-RuCl 3 , we use a nearest neighbor elastic model [41] on the honeycomb lattice:

where u i = (u a i , u b i ) and Π i = (Π a i , Π b i ) are displacement and its momentum at site i, α = a 1 (x), a 2 (y), a 3 (z) are vectors of three types of nearest-neighbor bonds. The matrix Φ α are given by:

(S14) Here, γ 1,2 are free parameters allowed by the hexagonal symmetry of the lattice, whose values we fix to match the experiments [42]. After the second quantization and Fourier transform, the phonon Hamiltonian can be written in a quadratic form:

where H ph ( k) is a 8 × 8 matrix (or 16 × 16 matrix for zigzag order unit cell) and

represents the phonon annihilation operator. The matrix H PH ( k) is given by

written in terms of the block matrices

Here the functions used are

where

To capture the behavior of phonons measured in experiments by Lebert et al. [42], we perform a linear square fit to the experimental phonon dispersions, and obtain E = 13.5 and c = 2.6. The band structure is shown in Fig. S2.

We now introduce the coupling between phonons and magnons, which is essential to endow the phonons with chirality and hence enable their intrinsic contribution to the thermal Hall effect (for extrinsic contribution, see Section V). A distance dependence of the superexchange interactions J(R i -R j ) between Ru 3+ ions leads naturally to the magnetoelastic coupling in a general form

where u i is the displacement of the i th ions from its equilibrium position, and ˆ r ij is the unit vector between site i and j. The bold J is the matrix notation of J αβ . In the second last line, the spin components S α and the matrix J are expressed in the Kitaev coordinate system (xyz). In the last line, the spins ˜ S are in the local coordinate determined by the magnetic order, where Sz is aligned with the direction of ordered spin. Their relation is given by S i = R i ˜ S i , where the rotational matrix R i is

where (θ i , φ i ) are the ordered spin polar angles in Kitaev coordinates.

The above equation (S22) is sufficient to derive all the terms for the magnon-phonon coupling, which will contain many free parameters of the form ∂J(rij ) ∂rij . Here we take the term from Sx A,i

B,i+ a1 as an example. Its magnon-phonon coupling term is:

where the magneto-elastic coupling can be as follows:

We then perform the standard Holstein-Primakoff transformation and write these displacements in terms of the phonon operators u γ i ∼ (A γ † i + A γ i ) (with polarization γ). This results in the hybridization between the magnons and phonons:

Other phonon-magnon coupling terms can be obtained in a similar fashion. Now the total Hamiltonian is given by

In the case of zigzag phase, it takes a quadratic form in the magnon and phonon operators:

where H( k) is a 24 × 24 matrix and

, ...), a k and A k representing the magnon and phonon annihilation operators.

The phonon-magnon coupling contains a few free parameters g 1 , g 2 , g 3 , g 4 defined in Eq. (S25). In order to fit these parameters, we refer to the recent experiment by Li et al in Ref. [43], where the phonon dispersions are measured at both high temperatures (in the paramagnetic phase) and low temperatures (ordered phase) (see Fig. S3(c)). The difference in these two scenarios is that the magnons are non-existent in the high-temperature paramagnetic phase, but interact with phonons at low temperatures. Hence, the downward shift of the phonon dispersion (i.e. phonon softening) seen in the low-temperature data is explicitly the result of the band anti-crossing, due to the phononmagnon coupling, as shown in Fig. S3(a). Hence, the magnitude of the phonon softening can be used to determine the magnitude of the magnon-phonon couplings, for parameter set ps2, given in Eqs. (S29,S31). For the magnon sector, we use the parameter set 2 (ps2, see Section II):

While the experimental data is insufficient to obtain a unique fit, using a simplifying assumption that all g coefficients are of the same order of magnitude. In the main text, we use

to qualitatively reproduce the phonon softening seen in the experiment. Please note that our units (defined Eq. ( S25)) are such that g = 1 corresponds to ∂J (r a1 )/∂r a1 = 4.9 meV/Å. We recognize that in reality the coefficients g i will likely have different magnitude, as shown by Winter et al. [49], however the above simplifying assumption is sufficient to illustrate the effect that the magneto-elastic coupling has on the thermal Hall effect. We also find that the maximal g i 's without inducing a phase transition are

Having thus determined all the coupling parameters in the Hamiltonian (Eq. S27), we compute the dispersions of the coupled magnon and phonon branches, which are shown in Fig. S4 for g i 's taking values in Eq. (S31). Figs. S3(a) shows the zoomed-in view of the phonon and magnon dispersions before (dashed lines) and after (solid lines) turning on the magneto-elastic coupling, which is to be compared with the experimental data in Figs. S3(b).

After having obtained the hybridized magnon and phonon bands, we then compute the κ xy again, which now includes intrinsic phonon contributions. For g i 's taking values in Eq. (S30), the result is shown in Fig. 4 in the main text. For g i 's taking the largest limiting values in Eq. (S31), the result is shown in Fig. S6 below. The convergence of our calculation for κ xy , which involves numerical integration over the Brillouin zone in Eq. (S11), is shown in Fig. S5, for the case of Eq. (S31).

Our final note of this section is that the coupling strength between magnons and phonons cannot be arbitrarily large. A coupling too strong will either drive the system out of the presumed magnetic order and/or induce a lattice distortion. In terms of the original band theory, this is manifested in unphysical band dispersions for values of magneto-elastic constants g in Eq. (S25) that are too large.

Given the complexity of the 12-band Hamiltonian (4 magnons and 8 phonon branches in the magnetic unit cell), we illustrate this point in the simplified picture. Let us consider one phonon and one magnon branch at a particular momentum. Both the phonon and the magnon can be thought of as harmonic oscillators. In the classical picture, they are scalar variables U and S sitting in two quadratic potential wells

where ω U and ω S correspond to the frequencies of the harmonic oscillators. As long as ω U and ω S are positive, the ground state for the system is a stable minimum U = S = 0.

Upon introduction of the interaction between the two harmonic oscillators

The two eigenvalues of the resulting matrix are given by

We can see here that when g 2 > ω S ω U , the lowest eigenvalue ω -will become negative, indicating that the corresponding quantum eigenstate (a linear combination of the U and S modes) can have an arbitrarily large occupation number, resulting in the unbounded negative total energy (unless one considers higher-order terms in the Hamiltonian on physical grounds). This indicates that the hybridized magnon-phonon eigenmode will condense, drive the system into a different phase. In our calculation, such types of phase transitions are not considered, nor are they realized experimentally in α-RuCl 3 [43].

Because of the strong spin-orbit coupling in α-RuCl 3 , the phonon momentum that is proportional to the elastic strain ε ij ∼ ∂ i u j couples to the defect spin as follows:

Following the derivation in Ref. [47], we arrive at the formula for the extrinsic (E) thermal Hall conductivity

where v is the phonon velocity and K is an appropriate skew-symmetric combination of the matrix elements of K ij,α (see Ref. [47]) for more detail. The universal function Φ(x) of the ratio x = ∆/T captures the distribution of two-level splittings on defects: if all defects have identical ∆, then Φ(x) = 1/ sinh(x), and if the splittings are drawn from a distribution, then the function must be averaged over this distribution (which generally results in a power-law in T ). Crucially, it follows that this extrinsic Hall conductivity is proportional to the phonon mean-free path = vτ ph , as is the longitudinal thermal conductivity κ xx , thus explaining the apparent sample-independent Hall angle θ E H = κ E xy /κ xx in the experiments [14].

We proceed to phenomenologically determine the extrinsic Hall angle θ E H from the experimental data in Ref. [14] at high temperatures T |K|/k B ≈ 80K, above the magnon bandwidth where the effects of the Kitaev physics and associated Berry curvature are unimportant:

We arrive at θ E H ≈ (0.6 ± 0.2) × 10 -3 . By contrast, at low temperatures of the order T ∼ 10 -15 K where the interpretation of the thermal Hall measurements on α-RuCl 3 is disputed, both the intrinsic (due to the Berry curvature) and extrinsic contributions to κ xy must be taken into account. In this low-temperature region

where η ≡ κ E xy /κ I xy is the phenomenological ratio of the extrinsic and intrinsic contributions. From the analysis of the data in Ref. [14], we thus obtain η = 1.2 ± 0.5, with the uncertainty related to the spread of the experimental data among the (five) samples. This range of obtained η ratios is used to determine the shaded blue region in Fig. 4(b) in the main text.

In the main text, a value of the magneto-elastic (ME) coupling g=4 (defined in Eqs. (S30) and (S25)) was used as it provided a reasonable match to the experimental data. Here, we would like to remark that if instead one used the largest allowed value of g ≈ 6 (see section III D for the definition), the intrinsic magnon+phonon contribution becomes even larger, resulting in κ I xy as large as 0.36 fermionic units at T = 10 K, as shown with a solid blue line in Fig. S6. This is still below the experimental value (red circles in Fig. S6), indicating the necessity to include extrinsic phonon contribution as discussed above. What this demonstrates however is that, upon including said contributions, the resulting κ 2D xy (T ) would fall inside the shaded blue region in Fig. S6, exceeding the experimental values for the given set of model parameters (ps2, see section II). Hence, even if the model parameters are such that they do not exactly coincide with ps2 chosen to maximize the intrinsic κ I xy , a significant region in the parameter space may host strong enough thermal Hall effect that will match with the experiment. FIG. S6. Total computed κxy/T , contributed from different sources for parameters given in Eqs. (S29,S31), compared with the experimental data (red circles) from Ref. [7]. Here, a larger value of the ME coupling g = 6 is used than in the main text. The magnon and phonon intrinsic component κ I xy (solid line), summed together with the extrinsic phonon contribution κ E xy , is indicated with the blue shaded region (whose width is given by the experimental uncertainty in determining κ E xy = ηκ I xy , see text). For these parameters, the theoretical thermal Hall effect is stronger than measured in the experiment, as a proof of principle that the model is capable of reproducing large κxy/T , provided extrinsic contributions are taken into account.