Breakdown of magnons in a strongly spin-orbital coupled magnet

The description of quantized collective excitations stands as a landmark in the quantum theory of condensed matter. A prominent example occurs in conventional magnets, which support bosonic magnons—quantized harmonic fluctuations of the ordered spins. In striking contrast is the recent discovery that strongly spin-orbital-coupled magnets, such as α-RuCl3, may display a broad excitation continuum inconsistent with conventional magnons. Due to incomplete knowledge of the underlying interactions unraveling the nature of this continuum remains challenging. The most discussed explanation refers to a coherent continuum of fractional excitations analogous to the celebrated Kitaev spin liquid. Here, we present a more general scenario. We propose that the observed continuum represents incoherent excitations originating from strong magnetic anharmonicity that naturally occurs in such materials. This scenario fully explains the observed inelastic magnetic response of α-RuCl3 and reveals the presence of nontrivial excitations in such materials extending well beyond the Kitaev state.

In this section, we compute the phase diagram for the extended (J 1 , K 1 , Γ 1 , J 3 ) spin Hamiltonian on the honeycomb lattice discussed in the main text, and given by: where {α, β, γ} = {y, z, x} for the X bonds, {z, x, y} for the Y bonds, and {x, y, z} for the Z bonds shown in Fig. 1 of the main text. In order to capture the effect of local quantum fluctuations on the state energies, we computed second-order energy corrections to the classical ordered state energies, extending the method of Supplementary Ref. [1] to include finite Γ 1 and J 3 . As discussed in Supplementary Refs. [1] and [2], an upper bound on the energies per site of the Kitaev spin-liquid states can be estimated from E KIT = ± 3 2 (J 1 + K 1 ) S γ i S γ j , where S γ i S γ j ≈ 0.131 is the analytical result for the first neighbour correlations at the pure K 1 > 0 or K 1 < 0 Kitaev points [3]. Comparison of E KIT with the second-order corrected energies of the ordered states has been shown to reliably predict the position of the phase boundaries, which agree with the results of exact diagonalization (ED) [1,2].
We show in Supplementary Figure 1a,b the phase diagram associated with Supplementary Eq. (1), parameterizing J 1 = cos φ sin θ, K 1 = sin φ sin θ, and Γ 1 = cos θ as in Supplementary Ref. [4], with J 3 = 0 (Supplementary Figure 1a) and J 3 > 0 (Supplementary Figure 1b). The present results may be compared directly with the ED results of Supplementary Ref. [4] for J 3 = 0. The extended model of Supplementary Eq. (1) exhibits six ordered phases, which have been identified in various previous works [1][2][3][4][5][6][7][8][9][10]: FM = collinear ferromagnetic order, AFM = collinear Néel antiferromagnetic order, ST = collinear stripy order, ZZ = collinear zigzag order, 120 = noncollinear 120 • order, and IC = incommensurate spiral order. In the extended model, the classical energies per site are given by: The full expressions with second-order corrections are lengthy, and therefore omitted here for brevity. However, the general trends are already apparent from Supplementary Eq. (2)- (7). Of particular interest are the regions of stability of the zigzag order, observed experimentally in the honeycomb materials α-RuCl 3 and Na 2 IrO 3 , as well as the extent of the K 1 > 0  [4]. (b) With a relatively small J 3 / J 2 1 + K 2 1 + Γ 2 1 = +0.088 consistent with the magnitude in the studied models. Small circles indicate the considered parameters. For the incommensurate spiral state, we did not consider quantum modifications to the ordering wavevector. The models studied in Supplementary Note 5 are highlighted by blue and red points. and K 1 < 0 Kitaev spin-liquids. For J 3 = 0, there are two zigzag regions, appearing at (J 1 < 0, K 1 > 0, Γ 1 ≈ 0), and (J 1 ≈ 0, K 1 < 0, Γ 1 > 0). The addition of a small finite J 3 uniquely stabilizes the zigzag and Néel states, linking the two zigzag regions, and suppressing the spin-liquid and other ordered phases. The existence of such J 3 coupling has been indicated for both α-RuCl 3 [8,10] and Na 2 IrO 3 [8,9,11]. This fact significantly complicates the identification of the magnetic interactions in these materials from investigations of the static properties alone, since the zigzag phase is stable over a very wide region of the phase diagram. The stability region for the K 1 > 0 and K 1 < 0 Kitaev spin-liquid states estimated from the second-order state energies is shown in Supplementary Table 1.
The relevant energy scale for considering the stability of the spin-liquid is the energy gap for spin-excitations at the pure Kitaev points, given by the two-flux gap ∆ ≈ 0.065|K 1 | [12,13]. Indeed, at the pure Kitaev points, the above estimates for the state energies suggest that the Kitaev state is stabilized with respect to adjacent magnetic orders by 0.072|K 1 | ≈ ∆ per site. In this sense, any perturbation on the scale of ∆ has the potential to destabilize the spin-liquid. This explains why the Kitaev spin-liquid occupies a relatively small region of the phase diagram.
In this section, we briefly review previous ab initio studies of the magnetic interactions in α-RuCl 3 . As discussed by some of the present authors in Supplementary Ref. [8], estimation of such interactions from first principles calculations are complicated by several factors: • The layered structure of α-RuCl 3 allows for significant stacking defects in the crystal structure, which have complicated structural solution. Very early studies indicated a trigonal space group P 3 1 12 [14,15], with ∠Ru-Cl-Ru ∼ 88 • . More recent detailed reanalysis of the structure suggested it to be monoclinic C2/m [16,17] at low temperatures, with ∠Ru-Cl-Ru ∼ 94 • . Furthermore, there is now increasing evidence that α-RuCl 3 exhibits a structural phase transition near T ∼ 150 K [18,19], which may be analogous to the C2/m → R3 transition of CrCl 3 [20].
Given that the magnetic interactions are highly sensitive to the local geometry of the RuCl 6 octahedra, accurate estimation of their values has historically been complicated by the uncertainty in the crystal structure. This had led to a variety of reported interactions for α-RuCl 3 , summarized in Supplementary Table 2.
• The low symmetry of the real crystal structures allows many independent terms to appear in the spin Hamiltonian. For example, up to 30 parameters are required to fully define the interactions up to third nearest neighbours [8]. As in the main text, most previous works have treated only a selection of such parameters, therefore assuming the interactions to be of higher symmetry than required by the crystal structure. This corresponds to an effective averaging of the interactions, which is likely to produce variations in the computed magnitudes of each parameter across different computational methods.
• For α-RuCl 3 , the underlying energy scales (Hund's coupling, spin-orbit coupling, and crystal-field splitting, for example) are all of similar magnitude, which makes the computed interactions highly sensitive to fine details (such as structure, and choice Structure Method sgn(K 1 ) of Coulomb parameters). This provides additional variations in computed interaction magnitudes appearing in the literature.
Despite these complications, estimates of the interaction parameters for α-RuCl 3 have broadly agreed across many different methods, when similar crystal structures are taken as input. A summary of previous ab initio calculations for such interactions is shown in Supplementary Table 2.
Ab initio studies based on the early P 3 1 12 structures [14,15] of α-RuCl 3 suggested antiferromagnetic K 1 > 0, with ferromagnetic J 1 < 0 and antiferromagnetic Γ 1 > 0 of similar magnitude [8,21]. Microscopically, the positive K 1 arises from large direct hopping between Ru metal sites [4], which has a dramatic effect for the short Ru-Ru distances in the P 3 1 12 structure. It should be emphasized that this hopping also generates large Γ 1 interactions, such that antiferromagnetic K 1 > 0 must always be accompanied by large Γ 1 in real materials. For this reason, pure (J 1 , K 1 ) nnHK interactions suggested in Supplementary Ref. [25], based on spin-wave fitting, are inconsistent with K 1 > 0, from a microscopic perspective.
Later ab initio studies based on the (more accurate) C2/m structures [16,17] have instead suggested K 1 < 0, J 1 ∼ 0, and Γ 1 > 0 [8,10,[22][23][24]. For reasons suggested above, there has been a relatively large spread of computed values obtained from different computational methods. However, if we do not favour any particular method, we might expect an appropriate starting point for analysis to appear at the average values over all studies: and K 1 < 0. As noted in the main text, excellent agreement between the ED and experimental neutron spectra is obtained for Model 2, with (J 1 /K 1 ) = +0.1, and (Γ 1 /K 1 ) = −0.5, which is completely consistent with the range of ab initio values appearing in the literature. In the main text we write: "a large decay rate is ensured by the following three conditions: large anisotropic interactions, deviation of the ordered moments away from the high-symmetry axes, and strong overlap of the one-magnon states with the multi-magnon continuum." Below we elaborate on these conditions for strong magnon decays.
The Hamiltonian for α-RuCl 3 is given by Supplementary Eq. (1). The diagonalization of the Hamiltonian requires rotation to the local reference frames of spins. Taking the zigzag ordering wavevector Q = Y, the ordered moment is expected to lie in the crystallographic ac-plane provided Γ 1 > 0, as for Models 1 and 2 of the main text, see also Supplementary Figure 2. This plane also contains the cubic z-axis. Thus, defining θ as the angle between the cubic z-axis and the ordered moment directionz, the spin operators can be rotated into the ordered coordinate frame via a rotation matrix: where θ is obtained by a minimization of the classical energy and depends on K 1 /Γ 1 . Then the Hamiltonian is whereJ ij is the bond-dependent exchange matrix in the local spin basis.
The coupling of the one-and two-magnon excitations is generated by the off-diagonal ("odd") terms that containS z iS x j andS z iS y j in the local reference frame. Extracting such terms explicitly from the exchange matrix (10) gives For the zigzag structure with Q = Y, there are five distinct bonds with respect to the values of J odd ij,x , J odd ij,y . Thus, for the X and Y bonds, J odd ij,x(y) = ±A ± , with and for the Z bonds, J odd ij,x(y) = B, where Thus, magnon decay vertices scale as ∼ (A ± , B). It is immediately apparent that the exchange terms do not contribute to the three-magnon coupling because of the collinear spin arrangement in a zigzag structure [26]. For the case of antiferromagnetic Kitaev coupling K 1 > 0, the ordered moments tend to align along the cubic z-axis, thus selecting θ ∼ 0 [27]. In Supplementary Figure 2a, we show the evolution of θ with |Γ 1 /K 1 |. For Γ 1 = 0, as in Model 1 of the main text, the ordered moment is exactly along the high-symmetry cubic z-axis (θ = 0), and the off-diagonal A ± and B terms vanish according to Supplementary Eqs. (12) and (13) above. As such, the decay of one-magnon excitations into the two-magnon continuum is forbidden in the zigzag phase for the pure Heisenberg and Kitaev interactions. In this case, all spin interactions in the local spin basis appear in the form ofS x iS x i ,S y iS y i , orS z iS z i , which contribute only to the terms of even order in magnon operators. As discussed in Supplementary Ref. [27], this ordered moment direction is stable with regard to small Γ 1 > 0, which shifts the moments to θ < 0, with θ ∝ |Γ 1 /K 1 |. In this case, the average magnitude of the two magnon decay vertex is expected to scale linearly with the off-diagonal Γ 1 couplings for small Γ 1 . That is Supplementary  Figure 2b, normalized to the overall magnitude of interactions k,m ∼ K 2 1 + Γ 2 1 , which sets the scale of the one-magnon bandwidth. One can see that the decay terms become significant compared to the one magnon dispersion for large Γ 1 .
For the ferromagnetic Kitaev coupling K 1 < 0, and Γ 1 > 0 (as in Model 2 of the main text), the situation is somewhat different. Finite Γ 1 terms rotate the ordered moment away from the cubic z-axis by an angle θ 90 • for any value of |Γ 1 /K 1 | 0.05 (see Supplementary Ref. [27]). Thus, the moments lie close to the cubicx+ŷ direction. In the rotated coordinate frame, both Γ 1 and Kitaev interactions contribute to the "odd" terms in (11), and thus induce two-magnon decays. Because of that, A ± or B are always large and scale with the magnitude of K 1 and Γ 1 , as shown in Supplementary Figure 2d. For that reason, in this region of the phase diagram, the magnon decay vertex Λ ∼ (A ± , B) ∼ (K 1 , Γ 1 ) is always expected to be of the order of the one-magnon bandwidth.
The Holstein-Primakoff bosonization of (11) yields the three-boson Hamiltonian H 3 where we extracted the coordination number z = 3. Here, Using (12) and (13) for S = 1/2 and the parameters of Model 2, the three-magnon coupling for bonds X, Y and Z are The quantities in (16) can be referred to as the real-space three-magnon vertices. For Model 2, their strength relative to the full magnon bandwidth W ≈ 7 meV is ≈ 0.4 − 0.5. Such a strong anharmonic coupling is a precursor of strong magnon decays.

Decay Kinematics
The consideration above ensures that the amplitude of the anharmonic magnon-coupling terms is significant in case of Model 2 as well as for a wider part of the phase diagram. However, the effect of this coupling on the single-magnon states depends crucially on the availability of the two-magnon states for decays. A strong overlap of one-magnon states with the two-magnon continuum for the high-energy magnon branches is virtually guaranteed by the presence of the low-energy branches, as is exemplified favour of a significant overlap of such kind for all branches of the spectrum, including the lowest one. That argument, together with a strong three-magnon coupling, in turn guarantees substantial magnon line broadenings.
The situation of interest is illustrated in a sketch in Supplementary Figure 3. For a simple magnon branch with a Goldstone mode at k = 0 one can see that the bottom of the two-magnon continuum is precisely degenerate with the one-magnon dispersion. Moreover, the density of the two-magnon states vanishes at the one-magnon energy, rendering decays impossible [26]. If a gap ∆ is introduced in the magnon energy at k = 0, the two-magnon continuum will still have a minimum at the same Γ point, but will be gapped with the energy 2∆.
Next is the case of the gapless modes occurring at finite ±k 0 , keeping an overall shape of the dispersion as simple as possible. It is obvious that now the two magnons with the total momentum q = 0 have the bottom of their continuum at Min( k + −k ) = 0. Thus, the two-magnon continuum must be below the one magnon states in an extended vicinity of q = 0. It is also obvious that this argument is robust against a finite gap at the (pseudo-) Goldstone point as the overlap of the magnon branch with the continuum survives for a finite ∆.
We mention that such a situation is not uncommon and describes almost exactly the case of a spiral antiferromagnet, which has an ordering vector at finite momentum, see [26]. Another generic case are the spin-ladder and 2D valence-bond-like antiferromagnets having band minima at a finite Q. The decay kinematic conditions and magnon decays are well-documented for them, both theoretically and experimentally, see [26] and [28,29] for a recent realization.
While, in general, one can formulate a complete list of checks for decay conditions for an arbitrary form of the magnon spectrum, see [26], the current consideration provides a clear and intuitive picture of why magnon decays must occur in Model 2 and also more broadly.

Decay Rates
Using standard diagrammatic rules with the decay terms, one can straightforwardly obtain magnon self-energies in the one-loop (Born) and on-shell approximations, Σ µk (ω) → Σ µk (ε µk ), both strictly within the 1/S expansion [26]. This neglects such effects as more complicated spectrum renormalizations and spectral weight redistribution away from the quasiparticle pole. One can argue that the real part of the self-energy should be neglected altogether as the renormalization of magnon energies is already built in by the choice of the model parameters. In practice, the LSWT parameters are indeed often chosen to best fit the observed experimental and/or numerical bands.
Altogether, this leaves us with the only remaining and yet the most important and physically distinct effect of the anharmonic terms: magnon decays. Thus, the self-energy of the magnon branch µ is Σ µk (ω) → −iγ µ k , and the decay rate in the Born approximation is We take J odd as an average value from (16), and Φ ηνµ q,k−q;k is a dimensionless function of eigenvectors from the diagonalization of the quadratic Hamiltonian.
The decay rate is related to a much simpler quantity, the on-shell two-magnon density of states (DOS), q = k 1 + k 2 , and the angles 1 , 2 , these are: For ↵  1, as is typical of conventional antiferromagnets, it can be shown that That is, the bottom of the two-magnon continuum is precisely degenerate with the one-magnon dispersion. For this reason, the density of two magnon states g(✏ q , q) is essentially vanishing at the one-magnon energies, which strongly suppresses the decay rate even if the decay vertex ⇤ is large.
We now consider the case where the gapless modes are shifted to ±k 0 . In this case, the low-energy dispersion is given by ✏ k = Min(|k ± k 0 | ↵ ). Consider the case where q = 0.
The choice k 1 = k 2 = k 0 yields Min(E(k 1 , k 2 ; 0)) = 0. That is, a pair of Goldstone modes with opposite momenta always yield a total energy and momentum of zero. However, Thus, the two-magnon continuum must extend below the one magnon states -at least in the vicinity of q = 0. Therefore, the finite density of states of two magnon states g(✏ q⇠0 , q ⇠ 0) implies magnon decay is kinematically allowed near q = 0 whenever the low-energy modes are (symmetrically) shifted away from the -point in the Brillouin zone. This is illustrated below: Finally, we return to the estimation of (k, m) for the models studied in this work. For 5 SUPPLEMENTARY FIGURE 3. Kinematic conditions for two-magnon decay. This sketch shows the kinematic conditions for two-magnon decay emphasizing the importance of the (pseudo-) Goldsone mode at a finite Q-value.
We propose to approximate the square of the dimensionless vertex Φ ηνµ q,k−q;k 2 in the decay rate (17) by a constant, thus eliminating the numerically complex and costly element of the calculation and bypassing its analytical cumbersomeness. As a result, the decay rate (17) is simply proportional to the on-shell two-magnon DOS (18) with f = Φ ηνµ q,k−q;k 2 and brackets implying averaging. We argue that this approximation is well-justified for the gapped systems and for models with lower spin symmetries. This is because Bogolyubov transformations are not singular without the true Goldstone modes and because the lower symmetry implies fewer restrictions on the decay amplitudes beyond the kinematic constraints in the DOS, see [26] and [30].
Another strongà posteriori justification of this procedure comes from the need of a regularization of the imprints of the two-magnon Van Hove singularities in the Born decay rate γ µ k (17) that are inherited from the two-magnon DOS (18). Such singularities are unavoidable, see [26], yet they are unphysical and must be regularized. Regularization procedures result in an averaging of singularities over the momentum space, thus complementing the averaging suggested in our approximation (19). We also note that performing such a regularization in case of the fully numerical calculation of the vertex Φ ηνµ q,k−q;k can be prohibitively costly. The method that we use to address the issue of singularities is referred to as the iDE method, e.g. [31,32]. It is a simplified version of Dyson's equation on the pole with only imaginary part of the equation solved self-consistently, which is in line with the approximations described above, Σ µk (ε µk + iγ µ k ) = −iγ µ k , Allowing the initial magnon to have a finite lifetime relaxes the energy and momentum conservations, thus removing singularities of the Born approximation and mitigating the unphysical largeness of the γ µ k 's. Technically, instead of the integral in (17) with a simplifying assumption of (19), the calculation of γ µ k now requires a recursive solution of which is, typically, a quickly convergent process. To provide a reasonable estimate of the constant f , we use our previous study of the XXZ model on the same (honeycomb) lattice in external field [30], for which analytical expressions for the dimensionless cubic vertices Φ ηνµ q,k−q;k can be obtained. In this model [30], the characteristic values of the three-magnon couplings J odd relative to the magnon bandwidth are close to the ones considered in the current study, single-magnon states significantly overlap with the high-intensity parts of the two-magnon continuum in high fields, also in a close similarity to the kinematics of the present work, and magnon decays were shown to be very significant. By comparing Born-approximation γ µ k for the XXZ model with the corresponding two-magnon density of states [30], we extract the value of f 1/9. Back-of-the-envelope estimates suggest this constant to be f 4/zn 2 , where z is the coordination number and n is the number of magnon branches (sites in the unit cell). This estimate comes from analyzing the structure of the dimensionless cubic vertices Φ ηνµ q,k−q;k in previous studies such as Supplementary Ref. [30]. As is clear from Supplementary Eq. (14), the real-space coupling of spin-flips affects nearest-neighbour Holstein-Primakoff magnons. Hence, the vertex in k-space contains an Energy (meV)  (19) and (20), respectively, for the four magnon modes with the higher values corresponding to the higher-energy modes.
analog of the nearest-neigbour hopping matrix. Averaging of its square yields with ∼ 1/z the inverse coordination number. The number of atoms in the magnetic unit cell n gives the number of independent magnon modes and, therefore, their wavefunctions are normalized by 1/ √ n. Since the vertex couples three magnons, its square is, thus, proportional to ∼ 1/n 3 .
The summation over such modes eliminates one power of n. The factor of 4 comes from the square of the symmetrization factor in the decay term. For the considered problem of n = 4 and z = 3 the value of f 1/12 is in a quantitatively close agreement with the value of f 1/9.

Dynamical Structure Factor
Our results for the dynamical structure factor I(q, ω) for Model 2 are presented in Supplementary Figs. 4 and 5. Supplementary Figure 4 shows I(q, ω) for the zigzag state with the ordering vector at the M -point [ −π, π/ √ 3 ], with that choice motivated by a close similarity of I(q, ω) along the XKΓYΓ MΓ k-path (see Fig. 3 where the broadening is obtained from the approximate Born expression (19) for the middle panel of Supplementary Figure 4 and by the iDE method (20) described above for the rest of the panels. The second panel in Supplementary Figure 5 shows exact diagonalization results from Fig. 3(a) of the main text. The solid lines in the middle panel of Supplementary Figure 4 show Born approximation decay rates γ µ q from (18) with the three-magnon coupling J odd from (16) and f = 1/9 as discussed above. For each k-point, there are four values of γ µ q , one for each branch, with the larger value corresponding to the magnon that is higher in energy. With γ µ q reaching about a half of the total magnon bandwidth, there is hardly anything visible left from the intensity of the highest-energy mode. The second highest mode is also overdamped in most of the k-space. This is in agreement with the high-energy magnons having a significant two-magnon continuum phase space for decays.
The two lower-energy modes are also broadened, with the modulations of the broadening following the k-regions where magnons do or do not overlap with the two-magnon continuum, such as in the vicinity of the Y and M points in Supplementary Figs. 4 and 5 for the latter case.
Last but not the least are the Van Hove singularities in γ µ q from the two-magnon DOS. Their imprints are also visible in I(q, ω) as sharp boundaries between more or less bright regions of intensity, hence more or less well-defined magnons. At k's close to such singularities, the decay rates become unphysically large, violate the perturbative nature of the expansion, and need to be regularized.
Physically, since the Van Hove singularities are affecting magnons that are already within the two-magnon continuum and are, thus, decaying, the most relevant method for such a regularization is the iDE approach described above, which allows to self-consistently account for the broadening of the initial-state magnons.
The iDE broadening γ µ q from (20) is shown in the lower panel of Supplementary Figure 4 by the dashed lines. As one can see, the iDE method yields smooth, completely regular decay rates, with their overall values decreased for the upper-and somewhat increased for the lower-energy modes. We emphasize again that this result is more realistic than that of the Born approximation, because the divergent behaviour violates the perturbative nature of the 1/S expansion and is unphysical.
Supplementary Figure 5 shows I(q, ω) averaged over three zigzag configurations. Our iDE results in the first lower panel of Supplementary Figure 5 clearly capture many of the most notable features seen in the ED data and constitute a clear improvement over the LSWT results. There are still some notable differences. First, in the vicinities of the Mand Y-points, ED branches at lower energies are flatter and asymmetric with a more drastic decrease of intensity in the Y→ Γ direction. There is also only one mode resolved near these points in the ED data, while the SWT predicts two. These may be ascribed to the ignoring of the real part of the self-energies, which neglects the spectrum flattening and the reduction of the quasiparticle peak intensity.
Another major remaining difference is the presence of a significant intensity in the ED data at the energies near and above the single-magnon band maximum, ω 7 meV, which is completely missing in both LSWT and iDE results. This missing feature is beyond standard calculations of the structure factor, which take into account only transverse spin fluctuations (see Supplementary Ref. [33] for an exception). The missing part is the longitudinal component, I zz (q, ω) in the local z-axes, which is directly related to the continuum of the broadened two-magnon excitations. While the full calculation of I zz (q, ω) is beyond the scope of the current work, a simple account of its ω structure is possible.
A naïve approach is to suggest a direct proportionality of I zz (q, ω) to the two-magnon DOS, similarly to Supplementary Eq. (19). However, a better approximation is achieved by modifying the density of states by including the broadenings of the magnon lines inside the continuum This modification is very physical and self-consistent as the continuum is built from the broadened magnons. It also eliminates sharp features in the continuum and because the γ µ q 's are larger for the upper magnon branches, the upper part of the continuum also gets washed out more. The single adjustable parameter f 2 can be argued to be 1/8 using naïve estimates. While the resulting intensity lacks a more involved modulation in the momentum, one can expect an overall better description of the ED results. The results shown in the lowest panel of Supplementary Figure 5 include I zz (q, ω) from (22) with f 2 = 1/8 and the averaged iDE decay rates: γ 1 q = 0.5 meV, γ 2 q = 1.0 meV, γ 3 q = 1.5 meV, and γ 4 q = 2.5 meV, where the numeration is from the lowest to the highest in energy. Thus, the broad features in the full I(q, ω) are a combination of the remnants of the broadened magnon modes from the transverse part of I(q, ω) with the longitudinal I zz (q, ω) part, so the combination has a maximum at the energies between 6 and 7 meV, in a close resemblance of the ED data. While this is only an approximate description, it provides confidence that a complete account should be able to reproduce other features of the ED data in that range of energies as well.
Altogether, we believe we have been able to provide a convincing description of the most significant effects of the three-magnon interaction on the magnon spectrum that agrees with the numerical studies.

SUPPLEMENTARY NOTE 4: FURTHER DETAILS OF EXACT DIAGONALIZATION CALCULATIONS
The exact diagonalization calculations in this work were carried out on the series of 20and 24-site clusters with periodic boundary conditions shown in Supplementary Figure 6. As noted in the Methods section of the main text, ED calculations were performed using the Lanczos algorithm [34], employing the continued fraction method [35] to obtain the desired dynamical correlation functions. While the periodic cluster 24A retains all the symmetries of Supplementary Eq. (1), the remaining clusters are of lower symmetry, resulting in slight anomalies in the symmetries of the computed correlation functions. For 24B, 20A, and 20B the symmetry was partially restored by averaging the results over all symmetryrelated orientations of the clusters, which generates the k-points shown in Supplementary  Figure 6(b).
In the following pages, we show complete results obtained for various models. For each model, comparison of results for each cluster is shown for the high-symmetry Γ, X, M(Y) and Γ points, which live on all four clusters. These results appear in Supplementary Figure 8-10(d,h,l,p). For each case, the lowest energy peak positions are relatively well converged with respect to finite-size effects (compared to the chosen 0.5 meV Gaussian broadening). For excitations representing a continuum, we observe variations in the positions of the higher energy peaks obtained from the various clusters, as might be expected. Averaging over the discrete excitations of the different clusters therefore restores the continuum, improving the validity of the computed intensities at the high-symmetry Γ, X, M(Y) and Γ points. However, we note that away from the high-symmetry points, where averaging is not possible, the intensities are less reliable. This observation does not alter the conclusions drawn from plotting I(k, ω) along various k-paths (as in Supplementary Figure 7), but should be noted.
In particular, in order to avoid spurious features, the plots of the k-dependence of I(k, ω) for various energy intervals (Fig. 2(d)  In this section, we show full results for various additional parameters not appearing in the main text, as well as a comparison of results from the different periodic clusters 20A-24B.

Nearest Neighbour Heisenberg-Kitaev (nnHK) Model
We first show, in Supplementary Figure 8, results obtained for the nnHK model with K 1 > 0 and J 1 < 0. We begin with J 1 = −4.6 meV, K 1 = +7.0 meV (Supplementary Figure 8(a-d)), as suggested from analysis of powder inelastic neutron scattering data in Supplementary Ref. [25], and then consider several models moving towards the spin-liquid, maintaining constant K 2 1 + J 2 1 . The validity of the ED approach can be seen by comparing the results for the pure Kitaev model (Supplementary Figure 8(n)) with the exact results (Supplementary Figure 8(o)). One can see good agreement between I(k, ω) predicted from the two approaches. A similar degree of agreement is seen for J 1 = −4.6 meV, K 1 = +7.0 meV, for which the ED results (Supplementary Figure 8(b)) and the LSWT results ( Supplementary Figure 8(c)) are also in close correspondence. This model exists close to the hidden SU(2) point (at K 1 = −2J 1 ) noted in Supplementary Ref. [7]. For this reason, the dynamics are expected to be well described by conventional spin-waves. Therefore, the LSWT method performs well both close to and far away from the spin liquid.
As discussed in the main text, upon approaching the spin-liquid, the intensive excitations at the 2D Γ-point shift to higher energy, and become increasingly broad as they move deeper into the three-magnon continuum. Away from the Γ-point, the intensive excitations are mostly associated with the lowest magnon band at the level of LSWT. These excitations remain relatively sharp in the ED calculations, and shift to lower energy on approaching the spin-liquid. This effect is clearly seen for the computed intensity at the X-point, shown in Supplementary Figure 8(d,h,l,p), which remains sharply peaked over the entire studied range. As discussed in the main text and Supplementary Note 3, low-energy magnons in the nnHK model are protected from decay due to the absence of three-magnon states at low energies. As shown in Supplementary Figure 8(n,o), these low energy excitations eventually evolve into a flat band in the spin liquid with intensity peaked just above the two-flux gap ∆ ∼ 0.065|K 1 | [12,13].
As in Supplementary Ref. [36], we find that the best agreement with the experimental results within the nnHK model is obtained for |J 1 /K 1 | = 0.3, which shows a star-like pattern at intermediate energy. This ratio of interactions was featured in the main text, for Model 1. However, we also find significant intensity at the X-points, inconsistent with the observed intensities. Moreover, as noted in the main text, the absence of low-energy intensity at the Γ-point in these models is strongly inconsistent with the experimental data of Supplementary Refs. [17,19]. We next consider parameters in the region suggested by ab initio studies of the earlier P 3 1 12 structure of α-RuCl 3 : |J 1 | ∼ |K 1 | ∼ |Γ 1 |, with J 1 < 0, K 1 > 0, Γ 1 > 0 [8,21]. It is worth noting that this region of the phase diagram features competition between ferromagnetic, zigzag, and 120 • order, such that the combination of K 1 > 0 and a finite Γ 1 > 0 tends to destabilize the zigzag order. In order to restore the zigzag ground state, we therefore add a small J 3 coupling. We start with the interactions suggested in Supplementary Ref. [25], namely J 1 = −4.6 meV, K 1 = +7.0 meV. We then add a small J 3 = 0.7 meV, and then vary the ratio of |Γ 1 /K 1 |, holding constant J 1 , J 3 , and J 2 1 + K 2 1 + Γ 2 1 + J 2 3 . Results are shown in Supplementary Figure 9.
Comparing Supplementary Figure 9(a-d) with Supplementary Figure 9(e-h), one can see that the addition of a small J 3 does not significantly influence the spectra. However, similar to the K 1 < 0 region, we observe significant broadening of the ED spectra upon increasing Γ 1 . For (J 1 , K 1 , Γ 1 , J 3 ) = (−4.6, +5.0, +5.0, +0.7) meV (Supplementary Figure 9(m-p)), LSWT predicts a large spin-wave gap, and relatively flat dispersion for the spin-wave bands. Interestingly, the momentum dependence of the predicted intensities from LSWT resemble somewhat those of the K 1 < 0 Kitaev spin-liquid, but shifted to higher energy. That is, there appears a flat band at 4-5 meV, with intensity centered around the 2D Γ-point, and another band at higher energies 10-12 meV, with intensity away from the Γ-point. The vanishing dispersion of these bands at the level of LSWT is likely related to close proximity to the phase boundaries between ferromagnetic, zigzag, and 120 • order, which would typically feature low-energy modes near the Γ, M(Y), and K-points, respectively.
Results for the ED calculations differ significantly from the LSWT intensities at large Γ 1 . In particular, the gap is significantly reduced, such that dispersing modes can be observed near the (M,Y)-points. This may result from shifting of the phase boundaries in the ED calculations compared to the semiclassical LSWT approach. Interestingly, for (J 1 , K 1 , Γ 1 , J 3 ) = (−4.6, +5.0, +5.0, +0.7) meV (Supplementary Figure 9(m-p)) we observe low-energy intensity at the Γ-point in the ED calculations, and a star-like shape at intermediate energies, consistent with the experimental data on α-RuCl 3 . However, in the high energy region, the intensity is mainly located away from the Γ-point, in contradiction with the experiment of Supplementary Ref. [19]. Finally, we show, in Supplementary Figure 10, complete results for a variety of models in the region suggested by various ab initio studies based on recent C2/m structures of α-RuCl 3 , that is J 1 ∼ 0, K 1 < 0, Γ 1 > 0, J 3 > 0 [8,10,[22][23][24]. In each case, we hold J 1 = −0.5 meV, J 3 = +0.5 meV and the overall scale J 2 1 + K 2 1 + Γ 2 1 + J 2 3 constant, and modify the ratio of |K 1 /Γ 1 |. We also show results for the K 1 < 0 Kitaev spin-liquid for comparison.
As discussed in the main text and Supplementary Note 3, a large Γ 1 interaction induces significant deviations between the ED and LSWT results, due to coupling between the oneand two-magnon excitations. For Supplementary Figure 10 For these models, the LSWT predictions for I(k, ω) are very similar despite remarkably different interaction parameters. However, the ED results differ substantially.
As noted in the main text and Supplementary Note 3, the requirements for strong coupling of the one-and two-magnon excitations include a deviation of the ordered moments from the high symmetry cubic axes. While finite Γ 1 interactions generally rotate the ordered moments away from the cubic axes, it is interesting to consider also the case where Γ 1 = 0, J 3 > 0, and K 1 < 0. In this case, the directions of the ordered moments are not completely determined at the classical level. For example, for the ordering wavevector Q parallel to the Z bond, the classical energy is minimized for any orientation of the moments in the xyplane. However, as noted in Supplementary Refs. [27,37] the cubic axes are selected by a quantum order-by-disorder mechanism, such that magnons are expected to remain stable in this limit. Indeed, comparison of the ED and LSWT results in Supplementary Figure 10 Finally, we note that the models studied in this region (such as Model 2 of the main text) display excitation gaps on the order of 0.5 meV (at the M(Y)-points) at the level of both ED and LSWT. This is in contrast with the neutron scattering results, which appear to show a gap on the order of 2 meV [25,38]. However, it is worth noting that the size of the excitation gap is strongly influenced by the relative magnitudes of the K 1 and Γ 1 interactions along each nearest neighbour X,Y and Z bond, which are not constrained to be equal by the symmetry of the real crystals [8,10]. To demonstrate this, we show LSWT results for Model 2 (Supplementary Figure 11), modified with an anisotropic K 1 and Γ 1 , consistent with the results of Supplementary Ref. [8]. Specifically, we show J 1 = −0.5, J 3 = +0.5, with K Z 1 = −5.0 + δ, K XY 1 = −5.0 − δ, Γ Z 1 = +2.5 + δ/2, Γ XY = +2.5 − δ/2. The gap can be reproduced already for small perturbations on the order of δ = 0.1K 1 , while the remainder of the dispersions are not strongly affected. Furthermore, additional small anisotropic interactions [39][40][41] are allowed by symmetry that may also contribute to the gap. For simplicity, we have neglected such additional terms in the main text, but expect that their inclusion would further improve agreement with the experimental observations. (i) [5.5, 8.5] (m)