Interaction-stabilized topological magnon insulator in ferromagnets

Condensed matter systems admit topological collective excitations above a trivial ground state, an example being Chern insulators formed by Dirac bosons with a gap at finite energies. However, in contrast to electrons, there is no particle-number conservation law for collective excitations. This gives rise to particle number-nonconserving many-body interactions whose influence on single-particle topology is an open issue of fundamental interest in the field of topological quantum materials. Taking magnons in ferromagnets as an example, we uncover topological magnon insulators that are stabilized by interactions through opening Chern-insulating gaps in the magnon spectrum. This can be traced back to the fact that the particle-number nonconserving interactions break the effective time-reversal symmetry of the harmonic theory. Hence, magnon-magnon interactions are a source of topology that can introduce chiral edge states, whose chirality depends on the magnetization direction. Importantly, interactions do not necessarily cause detrimental damping but can give rise to topological magnons with exceptionally long lifetimes. We identify two mechanisms of interaction-induced topological phase transitions---one driven by an external field, the other by temperature---and show that they cause unconventional sign reversals of transverse transport signals, in particular of the thermal Hall conductivity. We identify candidate materials where this many-body mechanism is expected to occur, such as the metal-organic kagome-lattice magnet Cu(1,3-benzenedicarboxylate), the van der Waals honeycomb-lattice magnet CrI$_3$, and the multiferroic kamiokite (Fe$_2$Mo$_3$O$_8$). Our results demonstrate that interactions can play an important role in generating nontrivial topology.

In contrast to electrons, the absence of a conservation law for the magnon particle number admits numbernonconserving many-body interactions and spontaneous decays [58]. Nonetheless, neglecting magnon-magnon interactions is widely assumed in studies on magnon topology, typically supplemented by the statements that they are anyway frozen out at low temperatures and/or negligible due to their 1/S smallness, where S is the spin length. However, both statements do not generally apply, rendering magnon-magnon interactions one of the most pressing issues in the field of magnon topology.
The few existing studies of how interactions between magnons influence topology have predominantly taken a pessimistic point of view by following the suspicion that interactions may render results obtained within the noninteracting limit obsolete. Chernyshev and Maksimov identified spontaneous magnon decays at zero temperature as the main obstacle [59]. For topological magnon insulators in kagome ferromagnets, they predicted a very large zero-field magnon broadening (damping) of the size of the topological gap. With the gap no longer well defined the notion of chiral in-gap states is jeopardized. Nonetheless, McClarty and Rau pointed out that since topological magnon gaps or band crossings necessarily occur at finite energy, the interaction-induced self-energy is non-Hermitian, opening up possibilities to explore non-Hermitian magnon topology [60]. An anisotropy of the magnon lifetime was identified as a crucial, indirect experimental signature of non-Hermitian topology. A direct detection of typical hallmarks of non-Hermitian topology, e.g., exceptional points, is still impeded by the global magnon linewidth broadening.
In principle, spontaneous magnon decays may be frozen out at high fields, which energetically separate the one-magnon from the two-magnon sector, rendering decays from the for-mer into the latter kinematically impossible [61]. Hence, the damping of magnons can be reduced in large fields, reinstating the notion of gaps and in-gap states. For example, field-polarized Kitaev magnets feature well-defined topological edge magnons despite interactions [62]. Unfortunately, the field polarization trick is not always applicable. If the nontrivial magnon topology relies on noncollinearity as, for example, in kagome antiferromagnets [63][64][65][66], large magnon damping is inevitable [67][68][69]. An exception to the rule are ferromagnetic skyrmion crystals that feature topologically nontrivial magnon bands at very low energy [54,57,70]. Even for ultrasmall skyrmions built from but a few tens of magnetic moments, the energetically lowest nontrivial magnon band exhibits damping much smaller than the band gaps to its adjacent bands [71]. Hence, ferromagnetic skyrmion crystals are a unique platform to study magnonic topology and the validity of the noninteracting theory. Similar results may hold for ferromagnetic bimeron crystals [72] and antiferromagnetic skyrmion crystals [73], a subclass of which [74] was shown to exhibit topological magnons at lowest possible energies [75], thereby minimizing the kinematically allowed phase space for decays.
The above outlined state of the art suggests that magnonmagnon interactions need to be suppressed for Hermitian magnon topology to be appreciated. Herein, we explore particle-number nonconserving interactions beyond their detrimental lifetime broadening effect. Rather than analyzing if magnon topology is present in spite of interactions, we look for nontrivial magnon topology because of interactions. Our main result is that interactions can break symmetries of the harmonic theory, causing topological phase transitions. Taking honeycomb ferromagnets as a contemporarily experimentally relevant example [18,26,31,[76][77][78][79][80][81][82][83][84][85], we account for interactions within many-body perturbation theory to demonstrate how zero-temperature interactions cause a spontaneous mass gap of magnonic Dirac cones. We construct an effective theory to show that the interaction-induced mass gap is topologically nontrivial. Hence, a nontrivial Chern marker and winding number are found, signaling topologically protected chiral edge magnons in finite samples, a prediction that we explicitly verify numerically; the handedness of the chiral edge magnons is linked to the magnetization direction. This finding establishes the existence of nontrivial Hermitian magnon topology of genuinely quantum mechanical origin, i.e., due to spontaneous decays at order 1/S , in sharp contrast to existing studies of Hermitian magnon topology, whose results could equally be arrived at by solving the linearized classical Landau-Lifshitz equation.
Our analysis reveals two unconventional mechanisms for topological magnon phase transitions from left to righthanded chirality (or vice versa): one field induced, the other temperature induced. We relate the field-induced transition to a topological transition in the two-particle decay contours, signaling a saddle point in the two-magnon continuum. The saddle point comes with a sign change in the interactioninduced Dirac mass term, which reverses Chern numbers. Hence, single-particle magnons can exhibit a rich topological phase diagram by virtue of their interplay with the two-magnon manifold. Moreover, we uncover the possibility of temperature-driven topological magnon phase transitions, brought about by thermally activated collision processes that compete with decays in that they cause Dirac masses of opposite sign. This result establishes temperature as an external control of topology, opening up new avenues in topological magnonics. Both types of topological phase transitions come with an experimental signature in transport experiments. Since nontrivial magnon topology has an impact on anomalous magnon-mediated transverse transport of spin and heat at finite temperatures [10,[86][87][88], the field or temperature-induced change in topology causes unconventional sign changes in spin Nernst and thermal Hall conductivities.
Furthermore, we identify two working principles that allow for tailoring interactions in a way that their influence on magnon damping is suppressed, without compromising their effect on the band gap opening: (i) Using magnetic fields to tune Dirac cones close to the lower threshold of the twomagnon continuum. (ii) Externally enhancing interactions because strong interactions-rather than completely wiping out the spectral weight of single-particle bands-expel the quasiparticle peaks from the continuum, reinstating long lifetimes. These two mechanisms, in principle, allow for arbitrarily long-lived topological magnons brought about by magnonmagnon interactions.
By means of exact diagonalization we extend our analysis to the ultimate quantum limit of S = 1/2 and demonstrate the qualitative consistency between the semiclassical spin-wave theory and a full-quantum treatment. We conclude that, in principle, there can be a competition between noninteracting and interacting origins of topology, whose relative influence can be tuned by external fields. Hence, a priori, it is not sufficient to restrict a discussion of magnon topology to the noninteracting limit. We give an account of how interactions could influence the magnon topology and/or magnonmediated anomalous transport of several experimentally available materials, which have previously been discussed in such contexts.
Overall, our results demonstrate that magnon-magnon interactions should not be dismissed as being detrimental to magnon topology. Quite in contrast, they may be utilized to stabilize topology, a finding that gives rise to the paradigm of "interacting topological magnonics" and contributes to the fundamental understanding of interacting topological matter without particle number conservation, in which single and many-particle sectors cannot be considered separately.
The remainder of this paper is organized as follows. In Sec. II, we introduce two models of honeycomb-lattice ferromagnets-one "achiral," the other "chiral." These two models differ in their symmetries, with the achiral model featuring a Dirac cone stabilizing symmetry and the chiral model breaking this symmetry (Sec. II A). Within spin-wave theory, the symmetry-breaking part of the Hamiltonian is shown to cause particle-number nonconserving many-body interactions (Sec. II B), which we account for within lowest-order perturbation theory (Sec. II C). These interactions are shown to gap out the Dirac cone of the chiral magnet (Sec. III B), in- Honeycomb-lattice ferromagnets with Dzyaloshinskii-Moriya interaction (DMI) are minimal models to realize nonconserved interacting Dirac magnons. Red arrows depict fieldpolarized spin moments, small arrows at the bonds indicate DMI vectors, and the large arrow shows the direction of the magnetic field. (a) Achiral in-plane magnetized ferromagnet with second-nearest neighbor DMI. (b) Chiral ferromagnet with interfacial DMI due to structural inversion asymmetry, as indicated by a substrate layer. The texture is out-of-plane field polarized. troducing nontrivial topology (Sec. III B 1) and chiral edge state (Sec. III B 2). The gap is demonstrated to depend on the magnetic field (Sec. III B 3), the strength of interactions (Sec. III B 4), and temperature (Sec. III B 5). A topological phase diagram is mapped out and the topological transitions are predicted to cause unconventional sign changes in transverse transport (Sec. III B 6). We complement the spin-wave analysis by exact diagonalization studies in Sec. IV. After a discussion in Sec. V, we finally conclude in Sec. VI. Appendices A-E provide more detailed information on selected arguments.

II. MODELS: ACHIRAL AND CHIRAL FERROMAGNETS ON THE HONEYCOMB LATTICE
We consider both achiral (A) as well as chiral (C) twodimensional ferromagnets on the honeycomb lattice, as depicted in Fig. 1. Their spin Hamiltonians read where the Zeeman energŷ accounts for a magnetic field B that acts on all spin operatorŝ S r (r points to a lattice site). In both cases, we consider with positive Heisenberg exchange J between nearest neighbors stabilizing a ferromagnetic phase. The two models only differ in their antisymmetric exchangê In Eq. (4a), the out-of-plane relativistic Dzyaloshinskii-Moriya interaction [89,90] (DMI) D z acts between secondnearest neighbors, with ν r,r = ±1 chosen negative (positive) for (counter-)clockwise circulation around a hexagon [cf. Fig. 1(a); clockwise circulation applies]. This type of DMI does not compromise the classical ferromagnetic ground state; hence the name "achiral magnet." To select a particular axis of magnetization, the magnetic field is applied within the honeycomb plane; in Fig. 1(a), we take B −y.
In contrast, in Eq. (4b), the interfacial nearest-neighbor Fig. 1(b)], favors noncollinear ordering, rendering the magnetic ground state a Néel spin spiral with a fixed chirality; hence the name "chiral magnet." Applying B z out of the plane, the classical ground state first changes into a skyrmion crystal [91], and then into the field-polarized phase, which is the case approximately for |B| 0.8D 2 /J [92]. In the remainder of this paper, we exclusively focus on the limit of field polarization.

A. Symmetry analysis
For discussing symmetries, we have to consider both the Hamiltonian and its magnetically ordered classical ground state. We drop the hat symbol to indicate ground state spin directions (Ŝ r → S r ).
Spectral Dirac points on the honeycomb lattice are stable in the simultaneous presence of inversion symmetry (P; also called sublattice or parity symmetry) and time-reversal symmetry (T ). Although ferromagnetism violates actual T symmetry (spin flip), there can be an effective time-reversal symmetry T = R(n, π)T , composed from T and a rotation R(n, π) by π in spin space about an axis n normal to the magnetization M ≡ 1 N N r S r to map the texture back onto itself (N is the total number of spins). If this rotation leaves all terms of the Hamiltonian invariant, T is a good symmetry of the magnet.
Achiral magnet. For M within the honeycomb plane, as in Fig. 1(a), we find that T = R(z, π)T is a symmetry of the system. After the action of T has flipped the spins, a rotation R(z, π) by π about the z axis maps the texture back onto itself and leaves the D z terms invariant. To see so, recall that which has no effect on D z (S x r S y r − S y r S x r ). Moreover, although DMI relies on local inversion asymmetry at a bond's midpoint [90], the achiral magnet is P invariant, with the center of inversion being located at a hexagon's center of mass. (Recall that both the spins as well as the DMI vectors are axial vectors.) Hence, both stabilizing symmetries of the Dirac cone are present.
Chiral magnet. There is no possibility to construct a T symmetry for the chiral magnet in Fig. 1(b) because there is no in-plane axis about which the spins can be rotated back without changing the DMI terms. This lack of a rotation axis traces back to the DMI vectors spanning a plane. Hence, the chiral magnet breaks T symmetry. Surprisingly, in spite of broken interfacial inversion symmetry, there is an effective parity symmetry P = R(z, π)P. After the action of P the substrate is on top of the magnet and, hence, the DMI vectors are opposite; it amounts to the mapping D → −D. By a R(z, π) rotation in spin space, we compensate for the minus sign, because the DMI terms always connect S z r with an in-plane spin component, the latter of which acquires a minus sign upon rotation [cf. Eq. (5)]. Hence, the parity-breaking effect of interfacial DMI can be "rotated away." As summarized in Tab. I, the achiral magnet obeys both symmetries, but the chiral magnet is only parity symmetric. Hence, stable Dirac points are expected for achiral magnets but a mass gap for chiral magnets. Although these symmetry arguments appear to be straightforward, the situation is more intricate, as we elobarate on in the following section.

B. Spin-Wave theory
Relying on a magnetically ordered classical ground state, we proceed with a Holstein-Primakoff (HP) transformation [93] where S is the spin length. The axes of the local reference frame read e ± r = (1, 0, ±i)/ √ 2, and e z r = (0, −1, 0) for the achiral magnet and e ± r = (1, ±i, 0)/ √ 2, and e z r = (0, 0, 1) for the chiral magnet; i 2 = −1. The bosonic operators obey the usual commutation relation [â r ,â † r ] = δ r,r . A Taylor expansion of effectively expands the Hamiltonian asĤ = ∞ p=0Ĥ p , where p denotes the number of bosonic operators. Formally, this is an . Up to quartic order, Hamiltonian (1) readŝ Importantly, the two magnets differ only in their DMI-induced number non-conserving three-particle HamiltonianĤ A/C 3 . The remaining terms, i.e., the ground state energy E 0 , the bilinear HamiltonianĤ 2 , and the number conserving four-particle HamiltonianĤ 4 , derive from the exchange energy, which is identical in both models.
Harmonic Theory. Linear spin-wave theory amounts to diagonalizing the Fourier transformed quadratic Hamiltonian and by a unitary transformation to the normal modesB k = U † kÂ k , whereB † . Theâ k,i 's are the Fourier transformed bosonic annihilators of a spin deviation (or HP boson) at the ith basis site. Similarly, theb k, j 's are the annihilators of a magnon in the jth band ( j = ±). The off-diagonal γ k contains the nearest neighbor bonds The distance between nearest neighbors is set to one. Using we find that E k = U † k H k U k = diag(ε k,− , ε k,+ ) contains the single-particle magnon energies (Note again that B 0.8D 2 /J must hold for the chiral magnet to become field polarized [92], i.e., for the ferromagnetic phase to be the ground state.) Hence, the spectrum of the free theory is identical to that of "spinless" electrons on the honeycomb lattice. The magnonic collective excitations, which come in two flavors, exhibit a graphene-like band structure, as depicted in Fig. 2 by white lines. The two bands touch linearly at the Dirac point energy at the corners (K and K points) of the hexagonal unit cell, where |γ k | = 0. From hereinafter, we measure the excitation energies relative to ε D , such that the harmonic Dirac cone appears at "zero energy." The Dirac cones have zero mass (no gap). Since the harmonic theory of a ferromagnet without quantum fluctuations coincides with the linearized equation of motion of classical spin vectors (Landau-Lifshitz equation without damping), the Dirac magnons are also said to have zero classical mass.
Anharmonic Theory. Note that the DMI did not enter the harmonic theory. This is because DMI vectors orthogonal to the magnetization direction connect longitudinal with transverse spin components (e.g.,Ŝ z iŜ x,y j in a local reference frame, where z is chosen along the magnetization direction), which contribute to theĤ p 's with odd p. The contributions toĤ 1 cancel, such that the lowest-order appearance of the DMI is in H 3 , which comprises magnon decay (one magnon decays into two) and coalescence (two magnons coalesce into one). This Dark blue/yellow color indicates zero/maximal two-magnon density of states. The red rectangle around the Dirac cone at the K point indicates the momentum and energy window shown in Fig. 4. For increasing fields, the single-particle energies are shifted upwards in energy by B, while the two-magnon continuum is shifted by 2B. Hence, for large enough fields, the continuum is shifted past the single-particle energies. In particular, the Dirac cones at the K and K points leave the continuum at B/(JS ) = 1.
particle number nonconserving interaction is in line with the DMI-induced nonconservation of spin, i.e., [ iŜi · m,Ĥ] ∝ D (or ∝ D z ), with m = M/M being the magnetization direction. Thus, the spin-orbit coupling-induced DMI connects the magnons to the lattice bath, into (from) which angular momentum can be dumped (drawn) during a magnon coalescence (decay).
The cubic Hamiltonian readŝ and accounts for an incoming HP bosonâ p,ν with momentum p being destroyed but two bosonsâ † k,λ andâ † q,µ being created. This process has to obey momentum conservation, p = k + q, modulo reciprocal lattice vectors. The Hermitian conjugate (H.c.) part in Eq. (14) describes the opposite process with two incoming and one outgoing boson. The HP interaction vertices V λµ←ν k,q←p A/C are given in Appendix A. Importantly, they are linear in DMI D z (achiral magnet) or D (chiral magnet) and independent of the exchange energy J.
Eventually, we work in the eigenbasis, such that the cubic Hamiltonian must be expressed in terms of normal modes, where the three-magnon vertices are constructed from the HP vertices and the eigenvectors given in Eq. (11). At higher-order spin-wave theory, we encounter the number conserving four-magnon interactions (Ĥ 4 ), which derive from the exchange energy. (The explicit expression ofĤ 4 can be found in Ref. 27.) Apart from being frozen out at zero temperature anyway,Ĥ 4 shares the symmetries ofĤ 2 and, hence, does not play a significant role. Section II C gives more detailed arguments.
We close with noting that we are faced with the theoretically intriguing situation that particular magnetic interactions (here: DMI) between spins appear exclusively as number nonconserving interactions between excitations. Since DMI is the crucial ingredient to break the effective time-reversal symmetry of the chiral magnet (cf. Sec. II A), DMI-derived interactions break accidental symmetries of the harmonic theory and a highly nontrivial renormalization of the Dirac cones is expected. Since it does not take any thermal excitation for a magnon to spontaneously decay, the genuinely quantummechanical effects ofĤ A/C 3 persist down to zero temperature, which is the limit we will explore first. Any gap of the Dirac cone is hence associated with a spontaneous quantum mass.

C. Many-body perturbation theory
To account for the effects ofĤ A/C 3 andĤ 4 , we invoke many-body perturbation theory for the single-particle Green's function, from which we ultimately extract the renormalized magnon spectrum in terms of the spectral function. The latter is related to the dynamical structure factor (i.e., the spin-spin correlation function) probed in neutron scattering experiments [94].
Up to order 1/S , the interacting (Matsubara) one-magnon Green's function is given by [95] where is the noninteracting Green's function, β = (k B T ) −1 is inverse temperature (k B being Boltzmann's constant), T τ orders imaginary times τ i , and averages · con (0) are taken with respect toĤ 2 over connected ('con') diagrams.
After applying Wick's theorem, the first-order perturbation proportional toĤ 4 yields self-energies associated with the Hartree diagram in Fig. 3(a). For the second-order perturbation (Ĥ 3 ), one finds 4 "forward bubble" diagrams [upper diagram in Fig. 3(c)], and 8 each of "upwards tadpole" [left Hartree diagram. The Hartree contribution in Fig. 3(a) is most conveniently evaluated before the diagonalization procedure, as demonstrated by Pershoguba et al. [27]. Extending their analysis to finite fields, we may approximate (given in HP basisÂ k ) at low temperatures, β(S J + B) 1. The constant contains a thermal-activation-like factor (with the dilogarithm Li 2 ), which accounts for field-freezing of four-magnon interactions. 1 Hence, without DMI, the renormalized magnon energies readε There are neither interaction-induced band gaps nor magnon lifetimes. The spectrum is merely compressed which reduces the Dirac velocity. Note also the 1/S 3 dependence of λ in Eq. (19), which reveals that the low-temperature influence of Hartree contributions is even weaker than the nominal order 1/S . Hence, we disregard the Hartree contributions in the remainder of our study. Tadpole diagrams. The self-energy of the two tadpole diagrams in Fig. 3 (given in eigenmode basisB k ). The summation is over momenta q in the Brillouin zone (BZ) and the indices α, β = ± label a self-energy matrix element. Due to the Bose factor, this self-energy is zero at T = 0 and, hence, a priori plays an insignificant role. Moreover, our numerical evaluation of Eq. (21) revealed that Σ tad,αβ k (T ) is zero at all temperatures, rendering tadpole diagrams irrelevant.
Bubble diagrams. We are thus left with "bubbles" [ Fig. 3(c)], whose self-energy (given in eigenmode basisB k ) has two contributions. The first term of Eq. (22) derives from the "forward bubbles" and is called "decay term" and the second one from the "circle bubbles" and is coined "collision term." 2

A. Interaction-induced gap opening at zero temperature: Spontaneous quantum mass
To analyze the many-body renormalized magnon spectrum, we calculate the spectral function from the retarded Green's function matrix At T = 0, only the spontaneous contribution from the +1 summand in Eq. (22) remains and the self-energy reads Hereinafter, we take S = 1 if not stated otherwise. We evaluate A k (ε) in the vicinity of the harmonic Dirac cones (cf. red window in "sharpness" (lifetime broadening), is stable. In contrast, signatures of a spectral gap opening are witnessed for the chiral magnet, suggesting that the Dirac magnons acquired a spontaneous quantum mass. These findings are in accord with the symmetry analysis in Sec. II A. Next, we focus on the K (and K ) point and study the evolution of A K (ε) with D (or D z ) and B. The results in Fig Comparing results at B/J = 0.1 and B/J = 1 (top and middle rows in Fig. 5, respectively), one finds that the quasiparticle line widths are larger in the former case. This is due to the fact that for B/(JS ) ≤ 1 the harmonic magnon energies overlap with the two-magnon continuum, rendering spontaneous decays kinematically possible. The momentum-resolved twomagnon density of states (DOS) reads and is indicated by color in Fig. 2. One can show that at k = K (or k = K ) its lower boundary has energy ε bound = 2JS + 2B, where the 2B arise from the two-magnon excitations being built from two single-particle excitations, each of which grows with B [cf. Eq. (12)]. Hence, by equating ε D = ε bound , we find that the harmonic Dirac cone leaves the continuum at the critical field The imaginary part of the self-energy follows the twomagnon DOS and, thus, exhibits a step at the lower boundary of the two-magnon continuum, see Fig. 6. According to the Kramers-Kronig relation, a step in the imaginary part translates into a logarithmic singularity of the real part. The explicit expressions for these singularities are derived in Appendix B; they read where δB = B − B c and C αβ is a constant. (Note that a tiny numerical linewidth of order 10 −4 J was assumed in the numerical calculations, which is why the divergence of the real part is cut in Fig. 6.) For B > B c , the imaginary part is zero and Sketch of an element of the self-energy matrix at the lower threshold of the two-magnon continuum for D/J = 0.15. At B = B c = JS and S = 1, the lower threshold of the two-magnon continuum is exactly at ε D [see Eq. (13)], where the imaginary part of the self-energy exhibits a step; ImΣ is zero outside of the continuum (ε − ε D < 0) and finite inside (ε − ε D > 0). As enforced by the Kramers-Kronig relations, the real part of the self-energy exhibits a logarithmic singularity which appears to be cut only because of a finite numerical linewidth: 0 + → 10 −4 J. Hence, ReΣ is much larger than ImΣ outside of the continuum. the real part still very large due to the singularity. Hence, considerably renormalized magnon quasiparticles with negligible damping are expected for magnetic fields that are slightly detuned from B c .
Figures 5(e) and (f) explicitly show this behavior for S = 1. For fields B < B c = J, the single-magnon excitations overlap with the two-magnon continuum and the quasiparticle peaks are blurred due to lifetime broadening. However, they become sharp once B > B c . Since the real part of the self-energy exhibits a logarithmic singularity, the magnon renormalization is particularly prominent at B ≈ B c . As B increases further, the degree of energy renormalization reduces only logarithmically. The same is true for the splitting observed in Fig. 5(f), which is also due to a real part of a self-energy as we show below.
To further understand the Dirac cone splitting in the chiral magnet, we continue our analytical considerations. The renormalized magnon energiesε k,± are the solutions of the Dyson equation is the inverse of the interacting matrix Green's function, containing the self-energy matrix as given by Eq. (25). At the Dirac points, k = K or k = K , the harmonic magnons are in resonance, ε K,+ = ε K,− = ε D and we may invoke the on-shell approximation by replacing the argument of the self-energy by ε → ε D . Since , the real part of the renormalized magnon energies becomes from which we read off that the diagonal self-energy Σ diag (ε D ) causes a shift and the off-diagonal self-energy Σ off-diag (ε D ) a splitting of the magnon energy at the Dirac point. Importantly, the achiral magnet has zero off-diagonal selfenergies at the K and K points, such that there is no splitting, i.e., no Dirac mass. In contrast, the chiral magnet has a nonvanishing Σ off-diag , which explains the splitting observed in Figs. 4(b) and 5(b), (d), and (f).
Logarithmic singularity. Within the on-shell solution of the Dyson equation in Eq. (30) the logarithmic singularity encountered in Eq. (28a) directly enters the renormalized spectrum. This indicates the failure of lowest-order (1/S ) perturbative spin-wave expansion, i.e., a violation of the assumption that |Σ k (ε)| ε k . Although the self-energy scales with D 2 and, hence, is hardly larger than 1% of the Dirac energy ε D for D/J = 0.15 [cf. Fig. 6], eventually, its real part grows indefinitely as B → B c . This defines a small field window about B c , within which Re|Σ k (ε)| ε D and the lowest-order-in-1/S perturbation theory is inconsistent. The logarithmic singularity may be cut due to higher-order diagrams which physically account for the two-magnon continuum and, in particular, its thresholds, being renormalized once its building blocks, i.e., the single-particle energies, have been renormalized. Unfortunately, such an analysis is severely hindered by the smallwavelength nature of the problem [67].
However, one notes that the spectral function in Figs. 5(e) and (f) does not show divergent quasiparticle peaks at B = B c = J. This is because it amounts to a non-perturbative solution of the Dyson equation that is no longer on-shell nor strictly of order 1/S . Alternatively, self-consistency may be enforced by iteratively solving the Dyson equation off-shell to cut the singularity [67]. We detail the latter idea in Appendix C. Importantly, the off-shell solution agrees very well with the spectral function. Since we also find independent numerical evidence for the spectrum by means of exact diagonalization in Sec. IV, we conclude that lowest-order perturbation theory is sufficient to qualitatively capture the signatures of magnonmagnon interactions.

B. Interaction-Induced Topology
We restrict our further analysis to the chiral magnet to obtain insight into the nature of the interaction-induced gap.

Effective Hamiltonian
In the limit B B c , the two-magnon continuum (∼ 2B) and the single-particle excitations (∼ B) are well separated, rendering damping zero because decays are kinematically forbidden. Moreover, since the distance between the two magnon bands is negligible compared to their distance to the continuum, the difference between Σ k (ε k,+ ) and Σ k (ε k,− ) is negligible as well. Hence, we set Σ k (ε k, ]/2 is the Hermitian part of the self-energy at ε D . Note the appearance of U k in Eq. (32), necessary for transforming Σ H k (ε D ), which is evaluated in the eigenbasis, back into the HP basis. H eff k describes a two-level system, which may be written as This decomposition in terms of Pauli matrices σ = (σ 1 , σ 2 , σ 3 ) and the unit matrix σ 0 resembles a spin-1/2 Zeeman Hamiltonian with d 0 k and the components of the "effective field" d k given by whereγ k = sgn(γ k ). For notational ease, we suppressed both the dependence on ε D and the "H" label. The eigenvalues of H eff k are given bỹ Importantly, at k = K and k = K , we obtain d 1 k = d 2 k = 0, but d 3 k 0, implying a mass gap. 3 Figure 7 shows a representative 3 At the Dirac point K , the definition of d 3 k in Eq. (34d) agrees with the definition of Σ off-diag discussed in Sec. III A. d 3 k within the entire Brillouin zone, from which it is obvious that the Dirac mass has opposite signs at the K and K points. Hence, the Berry curvature integrates to a nonzero Chern number According to the bulk-boundary correspondence [96,97], the nontrivial winding number [98] suggests a topologically protected chiral edge magnon in finite systems (see Sec. III B 2).

Interaction-Induced Chiral Edge Magnons
To verify that the interacting chiral magnet indeed features chiral edge magnons we proceed with simulating a honeycomb-lattice slab with open (periodic) boundary conditions in the x (y) direction (cf. coordinate system in Fig. 1). Hence, the edges exhibit a zigzag termination. A width of 12 honeycomb-lattice unit cells is chosen. Further details on the numerical implementation are given in Appendix D. At large fields (red and purple lines), the decay contours are centered about the midpoints of the lines from Γ to K . At a critical field, there is a topological transition of the decay contours from the low-field to the high-field case, as indicated by the green lines. This transition corresponds to a saddle point of the two-magnon continuum and, hence, a logarithmic singularity in the imaginary part of the selfenergy.
projections appear at k y = 2π/3 and k y = 4π/3. Here, due to finite-size effects, the Dirac cone projections are shifted in reciprocal space towards each other. Nonetheless, the flat edge states are clearly discernible.
Upon zooming into the relevant energy window [see Fig. 8(b)], we find that the edge states do not coincide with the flat harmonic bands (black lines). The edge modes have acquired both a uniform downwards shift and a dispersion due to many-body interactions. A spatially resolved calculation of the spectral function (details in Appendix D) reveals that they are chiral. States on the left (right) edge have positive (negative) slope as dictated by the winding number w in Eq. (38). This result confirms that interactions open a topologically nontrivial gap, which hosts chiral edge states, and establishes the notion of "interacting magnon Chern insulators." We point out the difference to conventional "harmonic" magnon Chern insulators in Sec. V B.

Field-Induced Topological Phase Transition
Having collected evidence that Σ off-diag (ε D ) = Σ −+ K (ε D )or, analogously, d 3 K in Eq. (34d)-is a useful indicator of topology, we may now explore its behavior in the regime B < B c , in which the Hermitian effective model (32) nominally is a bad approximation both due to lifetime broadening (non-Hermitian contributions) and the self-energy being considerably energy-dependent. Still, crucial information can be extracted, as we explain below.
First, let us study the energy-conserving decay contours of a Dirac magnon at the K point, which are depicted in Fig. 9. At small fields (blue and orange lines), there are two decay channels. The first one is a decay into a magnon close to the Γ point and another one close to the initial K point. The second channel consists of two decay products situated at a ring around the K point. As the field increases, the decay contours grow until they meet for a critical value B c ≈ 0.464JS (green line). Further increase of the field causes a topological transition of the decay contours, which are now centered about the midpoint of the Γ-to-K paths. Eventually, the decay contours shrink to a point and disappear at B = B c (red and purple lines). Similar contours are obtained for Dirac magnons at the K point; only the role of the K and K points is interchanged.
A topological transition of the decay contours at B c signals a saddle point of the two-magnon continuum. Hence, the two-magnon DOS exhibits a logarithmic singularity, which is also found in the imaginary part of the self-energy as shown in Fig. 10 Hence, we demonstrated that the effective model has predictive power even in the limit of overlapping single and twoparticle sectors, thereby linking geometric properties of the two-magnon DOS to single-magnon topology. One might conjecture that there is a more general principle at work, rigorously relating the topology of the decay surface to the band topology of single-particle bands, but this is beyond the scope of the present work and is left for future investigation.

Revival of magnon topology for strong interactions
As already mentioned in the introduction, band gaps ∆ε and, hence, in-gap chiral states are meaningful terms only if the broadening Γ + and Γ − of the gapped bands is smaller than the gap. Hence, ∆ε/ √ Γ + Γ − > 1 must hold. In principle, we may consider the following two sources for damping.
(2) As for the magnon-magnon interactions, we have seen in Sec. III B 3 that for B < B c = JS the single-particle states overlap with the continuum, resulting in damping. However, quite counter-intuitively, this reasoning holds only for weak interactions. As shown by Verresen et al. [99], who built upon the work of Gaveau and Schulman [100], strong interactions do not fully wipe out the quasiparticle but rather expel it from the continuum. The expelled single-particle state retrieves a long lifetime, which comes at the price of a reduced quasiparticle residue.  In principle, this mechanism can "revive" any singleparticle state and, in particular, topological magnon band gaps and chiral in-gap states: although being blurred out for weak interactions, the gapped bands (and the chiral edge states) may get expelled from the continuum by strong interactions, thereby suppressing damping. Let us present this effect for the gapped Dirac magnons of the chiral ferromagnet in Fig. 1(b). Figure 11(a) shows how the Dirac cone is expelled from the two-magnon continuum for increasing interactions, i.e., for increasing DMI D. With S = 1, and B/(JS ) = 0.8, the harmonic Dirac cone is energetically located well within the continuum, whose lower threshold is indicated by the horizontal dotted gray line. At D/J = 0, the Dirac magnon has zero mass and the two quasiparticle peaks are degenerate. Small D/J ratios lead to appreciable lifetime broadening, because decay processes are kinematically allowed. Hence, along a path in reciprocal space through the K point the spectral function resembles the harmonic Dirac cone with additional broadening [ Fig. 11(b)]. For D/J = 0.2, the damping is already so large that only a very blurred quasiparticle feature can be identified. The Dirac cone shape is almost invisible [Fig. 11(c)]. As D/J increases further, spectral weight is transferred from the continuum to the quasiparticle peak clinging to the lower boundary of the continuum and the lower band of the Dirac cone leaves the continuum [ Fig. 11(d)]. Eventually, also the second quasiparticle peak is expelled [ Fig. 11(e)]. Being no longer located within the continuum, both peaks are sharp. Collision density of states in honeycomb ferromagnets. At B = 0, the harmonic single-particle energies [white lines; cf. Eq. (12)] overlap with the collision density of states [color plot; cf. Eq. (40)]. Dark blue/yellow color indicates zero/maximal collision density of states. For increasing fields, the single-particle energies are shifted upwards in energy by B, while the collision continuum is fixed. Hence, for large enough fields, the single-particle bands are shifted past the continuum. In particular, the Dirac cones at the K and K points leave the continuum at B = B c = JS , similar to the two-magnon density of states in Fig. 2.

Thermal Topological Phase Transition
We have already seen in Sec. II C that the most important Feynman diagrams are the bubble diagrams, because the tadpoles integrate to zero and the Hartree term merely compresses the spectrum. Hence, as far as magnon topology at finite temperatures is concerned, we may also restrict to the bubble self-energy in Eq. (22). So far, we considered the decay term at zero temperature. At finite temperatures the spontaneous contribution is accompanied by two Bose factors ρ(ε q, j , T ) + ρ(ε k−q, j , T ), which enhance the efficiency of decays. Still, as long as decays are kinematically forbidden (B > B c ), finite temperatures will not cause damping.
A qualitatively new contribution comes from the thermally activated collision term which is proportional to a difference of Bose functions: ρ(ε q, j , T ) − ρ(ε k+q, j , T ). It may contribute to damping, if the single-particle energies overlap with the collision DOS which is depicted in Fig. 12. However, one verifies that the single-particle energies leave the collision DOS at the same critical field B c at which they also leave the two-magnon DOS. Hence, for B > B c , damping is strongly suppressed at all temperatures (within 1/S perturbation theory). The most relevant influence of temperature is then found in its contribution to the real part of the self-energy, with logarithmic singularities arising both from the (lower) two-magnon as well as the (upper) collision DOS thresholds. Figure 13(a) shows the temperature dependence of the magnonic spectral function A K (ε) for B = 1.1JS > B c and S = 1 at the Dirac cones, taking only bubble diagrams into account. The two quasiparticle peaks stay well defined, in agreement with the aforementioned kinematic reasoning. Just below k B T/J ≈ 2, there is a band gap closing and reopening. Within the effective model (32), this transition is associated with a vanishing of the off-diagonal self-energies and, hence, of the mass term d 3 k in Eq. (34d). One verifies that d 3 k flips sign at the K and K points, as shown in Figs. 13(b)-(d), leading to a temperature-driven topological phase transition from negative chirality (w = −1) to positive chirality (w = +1).
To get some intuition, why such a transition can happen, let us ignore the interaction vertices in the self-energy [Eq. (22)]. Then, decays and collisions are fully characterized by the respective DOS. Note that the two-magnon DOS in Eq. (26), the single-particle energies in Eq. (12), and the collision DOS in Eq. (40) respectively grow with 2B, 1B, and 0B. At the K and K points, for B = B c = JS , the lower threshold of the two-magnon DOS, the Dirac cone in Eq. (13), and the upper threshold of the collision DOS coincide in energy. Thus, for B > B c , the two-magnon (collision) DOS is energetically above (below) the Dirac cone, whence it follows that (for all q) and the two denominators in Eq. (22) have opposite sign. Consequently, decays and collisions cause mass terms of opposite sign. At low temperatures, when collisions are frozen out, decays dominate. However, thermally activated collisions may take eventually over, causing the topological transition.
A comment on the consistency of spin-wave theoy at large temperatures is due. Above, we pushed the lowest-order anharmonic spin-wave theory to temperatures k B T > J. In this limit, temperature-enabled higher-order 1/S corrections can be expected to alter our findings. For example, at order 1/S 2 , sunset diagrams, as obtained within second-order many-body perturbation theory inĤ 4 , are the leading source of damping [27]. The resulting magnon lifetimes will blur the quasiparticle peaks in Fig. 13. Whether or not the topological transition in the spectrum will be visible depends again on the ratio of the spectral gap to the damping (cf. discussion in Sec. III B 4). Moreover, thermal fluctuations eventually destroy the magnetic order. In Sec. III B 6, we estimate the region in parameter space spanned by B and T within which field freezing efficiently counteracts thermal fluctuations.
Thus, the main insight of our analysis of temperature effects is that thermal fluctuations habor, in principle, the potential to cause topological phase transitions. For more quantitative predictions, higher anharmonic spin-wave interactions have to be taken into account.

Interaction-Induced Topological Phase Diagram and Transverse Transport
Armed with the knowledge of both field-induced (cf. Sec. III B 3) and temperature-induced topological phase transitions (cf. Sec. III B 5), we can now map out the full topological phase diagram shown in Fig. 14. To that end, we assume the on-shell approximation and extract the winding number w from the sign of the gap in Eq. (31). We point out that this is an approximation, the accuracy of which may be judged by comparing the gap closing in Fig. 13 (below k B T/J = 2; obtained from the spectral function) with the temperature-induced topological phase transition in Fig. 14 (above k B T/J = 2; obtained within the effective model). Despite this slight quantitative disagreement all qualitative features are captured. We find that there are no topologically trivial phases (w = 0).
Since finite temperatures try to destroy magnetic order, we have to make sure that the relative magnetization (per spin) evaluated within linear spin-wave theory, is well above zero. Dashed lines in Fig. 14 indicate lines of constant M. We find that the temperature-induced topological transition happens at temperatures, at which M is already considerably reduced. Hence, large fields B/(k B T ) > 1 are necessary to freeze the magnetization and to appreciate this transition.
We note that such a topological phase diagram can never be obtained within the noninteraction theory of ferromagnetic topological magnon insulators. This is because (i) the influence of temperature on the spectrum escapes the harmonic theory as it may only be incorporated as an effective scaling and (ii) external magnetic fields cause but a mere uniform energetic shift to the single-particle excitations [cf. Eq. To search for experimental signatures of these transitions we propose to perform transport measurements. We recall that time-reversal breaking is one of the two necessary ingredients for transverse particle transport, the other being compatibility of the magnetic point group with ferromagnetism. With both requirements met by the chiral magnet, we expect anomalous transverse magnon transport at finite temperatures, namely, Hall, Nernst, and Righi-Leduc effects (also called thermal Hall effect). They are respectively quantified by off-diagonal elements of the magnetization conductivity σ, magnetothermal conductivity Υ, and thermal conductivity κ. These conductivities relate magnetization currents j and heat currents q to gradients in magnetic fields, ∇B, and temperature, ∇T . The constitutive equations read and the conductivities are defined as σ = L (0) , Υ = L (1) /T , and Notice that κ is composed from two contributions, with the second one deriving from a particle backflow, establishing a magnonic Wiedemann-Franz law 4 both for longitudinal [101] as well as transverse transport [102]. The transport tensors L (i) with i = 0, 1, 2 are 2 × 2 matrices, whose elements read 4 The magnonic Lorenz number reads L ≡ κ xx /(T σ xx ) = Ck 2 B /(gµ B ) 2 , with g-factor g, Bohr's magneton µ B and a numerical constant C, whose value depends on the dimension [101,102]. By restricting the Boltzmann transport theory laid out in Refs. 20 and 22 to two dimensions, we obtain C = 2. [10,20,22,86,[102][103][104][105][106][107][108][109] where µ = x, y; moreover, note that L (i) yx (T ) = −L (i) xy (T ). We introduced the g-factor g, Bohr magneton µ B , the transport relaxation time τ k, j , the group velocity v µ k, j , and the weights First, we study the off-diagonal elements L (i) xy (T ) in Eq. (45a) and note that we restrict our analysis to intrinsic contributions to transport. Thus, the transverse conductivities are related to the magnon Berry curvature Ω k,± in Eq. (36) of the effective Hamiltonian H eff k . Note that temperature enters both explicitly via the thermal weights c i (x) as well as implicitly because H eff k and Ω k,± inherit the temperature dependence of the self-energy in Eq. (22). Figure 15 shows the dimensionless off-diagonal entries of the transport tensors, in dependence on selected parameters. At low fields,L (i) xy is positive over a large temperature window [cf. Fig. 15(a)], which complies with the positive-winding low-field phase in Fig. 14 (lower yellow phase with w = +1). In contrast, at higher fields,L (i) xy is negative but the temperature-induced transition (from the cyan w = −1 to the yellow w = +1 phase in Fig. 14) eventually flips the sign of the dimensionless conductivities [cf. Fig. 15(b)]. Another sign change is found for a magnetic field sweep, as shown in Fig. 15(c). Large fields eventually freeze out magnonic transport because the singleparticle bands are shifted towards high energies, suppressing thermal occupation. Finally, an increase in DMI D also causes an increase in transverse conductivities [cf. Fig. 15(d)].
The sign changes of L (0) xy and L (1) xy may be detected by transverse spin transport experiments [110]. However, the most frequently measured signature of transverse magnon transport are heat currents [111][112][113]. According to Eq. (44), L (2) xy is not sufficient for a quantitative prediction because the offdiagonal element of L (1) (L (0) ) −1 L (1) has to be subtracted. In contrast to electrons, whose internal energy scale given by the Fermi energy renders this correction negligible at low temperatures, bosons do not come with such a scale [101,102]. Hence, we need the full L (0) and L (1) tensors, in particular, their diagonal elements given by Eq. (45b). We approximate the transport relaxation time by the magnon lifetime due to Gilbert damping α, i.e., τ k, j ≈ /(2αε k, j ). Setting α = 10 −3 results in κ xy /κ xx ∼ 10 −3 , a ratio typically found in experiments [111,113]. Moreover, to obtain more wieldy units, we study the three-dimensional thermal conductivity κ 3D = κ/ , which is given as the ratio of conductivity per layer, κ, and a typical interlayer spacing = 1 nm. Thus, the units of the thermal conductivity are [κ 3D ] = W Km . We decompose the transverse thermal conductivity as Figures 16(a,b) show the two contributions to the thermal conductivity and Fig. 16(c) their sum for a realistic ratio of D/J = 0.15. Importantly, the sign of the correction term [ Fig. 16(b)] is opposite to that of the nominal term [ Fig. 16(a)], which is a result of the backflow current. Due to the missing internal energy scale of magnons (zero chemical potential) both contributions have a similar magnitude. However, since the backflow cannot be larger than what drives the magnonic density imbalance in the first place, |κ 3D, nominal xy | > |κ 3D, correction xy | holds, and the total thermal Hall conductivity κ 3D xy [ Fig. 16(c)] has the same sign as the nominal contribution. Consequently, shows the same features asL (2) xy [cf. Fig. 15], i.e., a sign change with both magnetic field and temperature. (b) The correction term κ 3D, correction xy has a sign opposite to that of the nominal contribution and is of the same magnitude. (c) The total thermal Hall conductivitiy κ 3D xy exhibits features identical to that of κ 3D, nominal xy but its magnitude is reduced. The inset shows κ 3D xy for fields B/J > 1. For each magnetic field value, we define an upper threshold temperature T by M(T ) = 0.5 [see Eq. (42) and Fig. 14]. For T < T , the spinwave theory carried out in this manuscript is expected to be a good approximation. In contrast, temperatures T > T require higherorder spin-wave corrections. Parameters are S = 1 and D/J = 0.15. a measurement of κ 3D xy probes the unconventional topological phase transitions associated with magnon-magnon interactions. Sign changes of κ 3D xy are expected both for the fieldinduced as well as the temperature-induced transition. We find that κ 3D xy is of the order of 10 −4 W Km to 10 −3 W Km . Hence, the interaction-induced thermal Hall conductivity can be as large as that obtained within the free theory [62,87] and, in particular, as that frequently observed in experiments [111][112][113]. It should be considered when interpreting transport experiments.
We reiterate that we concentrated on intrinsic transport and worked with the effective Hamiltonian constructed in Sec. III B 1. If skew scattering and side jump effects due to disorder can be neglected, this approximation is expected to capture the relevant physics. Still, future theoretical stud-ies may develop a full many-body theory of transverse transport effects. In principle, one would not only expect intrinsic contributions due to the magnon Berry curvature but also skew-scattering-like extrinsic contributions due to asymmetric three-magnon scattering [114]. So far, we relied on a large-S approximation, which we pushed to S = 1, a procedure which may be met with scepticism. To confirm the validity of this procedure, we push our analysis even further down to the ultimate quantum limit of S = 1/2 and compare the results to those obtained from exact diagonalization (ED) of the chiral magnet. Simulations are performed on a cluster of 12 honeycomb unit cells (24 spins in total) with tilted periodic boundary conditions. For D = 0, the z-component of spin is conserved and the eigenstates can be labeled by their spin ∆s z = 1, 2, 3, . . ., with ∆s z = 1 being magnons, ∆s z = 2 being two-magnon states, and so on. Large fields energetically separate distinct ∆s z sectors, facilitating a comparison with spin-wave theory. Figure  17(a) shows the eigenspectrum (red crosses) as obtained from ED at B/J = 6. The states lowest in energy are ∆s z = 1 states that excellently agree with the single-particle magnon dispersion obtained within spin-wave theory [Eq. (12) for S = 1/2]. In particular, the two ∆s z = 1 states are degenerate at the K and K points. At these points, one finds that the ∆s z ≥ 2 states cross the single-magnon states for sufficiently low fields, as depicted in Fig. 17(b). At fields slightly above B c = JS = 0.5, the single-magnon states cut into the two-magnon manifold. However, due to spin conservation there is no hybridization between different spin sectors, as can be verified in Fig. 17(c).
Hybridization between spin sectors is brought about by finite DMI (D/J = 0.15) that violates spin conservation. As depicted in Fig. 17(d), the Dirac magnons split and get downwards renormalized as they are approached by two-magnon states. The splitting increases quadratically with growing D/J, as shown in Fig. 17(e) for B/J = 1.
Overall, we find very good qualitative agreement with the anharmonic spin-wave theory, e.g., with Fig. 5(d) and (f). We refrain from a quantitative comparison because ED is prone to finite size effects. By changing the simulated cluster geometry the gap between single-magnon states can vary considerably; the existence of the gap, however, is not subject to finite-size effects. (In Appendix E, we use ED also to show that the Dirac magnons are gapped for S = 1.) We are thus content with concluding that spin-wave theory is a surprisingly good approximation even for S = 1/2 and attribute this finding to the collinear texture and negligible frustration.

A. Working principles for undamped interacting topological magnons
From our results, we identify two working principles to minimize the damping caused by interactions without compromising their band-gap-inducing quality.
(i) Tuning the Dirac point just below the lower threshold of the continuum. To take full advantage of the logarithmic singularity in the real part of the self-energy (cf. Sec. III A) the Dirac cone must be energetically positioned close to but still outside of the continuum. This is easily realized by application of a magnetic field. Real materials may only be quasitwo-dimensional due to some (weak) interlayer coupling. Any degree of three-dimensionality regularizes the DOS of the two-magnon continuum and, hence, also the singular behavior of the self-energy. The logarithmic singularity and step of the real and imaginary part are replaced by a nonsingular squareroot behavior [67]. Still, the gap-inducing real part is strongly enhanced at the threshold of the continuum.
(ii) Artificially enhance interactions to expel the singleparticle states from the continuum. If principle (i) is not an option, e.g., because it would take inaccessibly large magnetic fields, one can try to make the continuum "expel" the single-particle states by artificially increasing interactions (cf. Sec. III B 4). It is known that strong electric fields E cause DMI interaction because they break inversion symmetry [115][116][117]. Hence, by applying an out-of-plane E, which mimicks the electric field due to a structural potential gradient, Rashba-type DMI [as that considered in the chiral magnet in Fig. 1(b)] is enhanced. Moreover, recently, sereval ideas on light-induced effective magnetic interactions were discussed [118][119][120][121][122][123], which, when carefully designed, could cause considerable magnon-magnon interactions.

B. Comparison between harmonic and interacting magnon Chern insulators
We reiterate that it is only due to the particle-number nonconserving interactions that the broken (effective) timereversal symmetry of the chiral magnet in Fig. 1(b) becomes apparent (also cf. Sec. II A). It is worth pointing out the difference to conventional magnon Chern insulators (e.g., Refs. 11 and 124) that exhibit nontrivial topology already at the harmonic level. There, a reversal of the sign of the DMI flips the winding number and, hence, the edge state's chirality. This is because the harmonic theory sees D · M, i.e., the projection of the DMI vectors D onto the magnetization M. Hence, a reversal of D is equivalent to a reversal of M, the latter being the actual crucial ingredient to flip the sign of topological invariants. In contrast, the gap opening due to interactions observed in Sec. III A is quadratic in DMI and, thus, independent of a sign flip of DMI. We argue that there is no necessity for the sign of DMI to be the relevant parameter. Instead, one verifies that it is still the magnetization reversal, which leads to a sign flip of the mass and, hence, also to a sign reversal of the winding number. To see this, we recall that upon magnetization reversal not all spin components acquire a minus sign but that the sign of one in-plane component stays invariant, e.g., for a rotation R(y, π), the sign of S y r stays invariant. Effectively, this is like flipping the sign of the x components of the in-plane DMI vectors. Hence, the DMI-induced phases ϕ δ i = arg(d y δ i − id x δ i ), which enter the interaction vertices via terms of the form e i(ϕ δ i −k·δ i ) in Eqs. (A1a)-(A1d), are reversed, ϕ δ i → −ϕ δ i . For the interaction vertices, such a flip amounts to the mapping (V λ,µ←ν k,q←p ) C → (V λ,µ←ν −k,−q←− p ) C, * , i.e., to complex conjugation and reversal of crystal momentum, the two operations usually associated with time reversal. Hence, upon magnetization reversal, the interaction-induced mass term at the K point is now found at the K point and vice versa. Obviously, the winding number in Eq. (38) has to change sign. One may also discuss what happens under a gradual rotation of M from up-pointing to down-pointing. Somewhere inbetween there has to be a topological phase transition associated with a gap closing. As M is fully rotated into the plane, the chiral magnet in Fig. 1(b) exhibits an effective timereversal symmetry T C 2,z . Taking a hexagon's center of mass as center for the C 2,z rotation, which rotates both spins and space (and, hence, maps the DMI vectors onto each other), actual time-reversal T maps the in-plane texture back onto itself. Hence, Dirac cones are stable, i.e., the band gap has closed.

C. Relevance in real materials
In real materials, magnonic topology is determined by a competition between harmonic and anharmonic mechanisms. For example, if we include a small D z second-nearest neighbor term [as considered in the achiral magnet, Fig. 1(a)] into the model of the chiral magnet [cf. Fig. 1(b)], a Dirac mass and nontrivial topology are already found at the harmonic level [124], with the winding number determined by the sign of D z . Additionally, sublattice-dependent on-site potentials, e.g., in the form of on-site anisotropies, may be considered. They break inversion symmetry and open a topologically trivial gap at the harmonic level. As is known from the electronic Haldane model [5], a rich topological phase diagram with both topologically nontrivial and trivial phases is found at the single-particle level. For magnons, the effects of particlenumber nonconserving interactions due to in-plane DMI (D) come on top. Since D z and D can be independent-the former derives from spin-orbit coupling intrinsic to the material, while the latter is due to structural asymmetry, e.g., in the presence of a substrate-D may be larger than D z . Then, despite the 1/S smallness and the quadratic influence of D on the dispersion, interactions may counteract and even undo or reverse harmonic topology, causing an opposite winding number or a transition from a topologically trivial into a nontrivial phase. With these effects being field as well as temperaturedependent, both magnon topology and transverse transport become highly nontrivial subjects. Since experimental studies of magnon topology and transverse transport have so far relied on the harmonic approximation for the interpretation of results, a careful reexamination may be necessary. Let us quickly go over selected materials for which there is experimental evidence compatible with nontrivial magnon topology and/or thermal Hall effects of magnonic origin. We restrict the discussion to (almost) collinear magnets. 5 5 For collinear magnets, the role of magnetic interactions in either harmonic (S ± i S ± j ) or cubic Hamiltonians (S ± i S z j ) is clear, facilitating the discussion Kagome ferromagnet Cu [1,3-benzenedicarboxylate(bdc)]. The metal-organic kagome ferromagnet Cu(1,3-bdc), which exhibits Dirac gaps [17] as well as transverse thermal transport [113] for an out-of-plane field, orders in-plane with (negligible) antiferromagnetic inter-plane coupling in the absence of an external field [125]. Similar to the achiral model magnet in Fig. 1(a), the in-plane ordered state exhibits D z -induced magnon damping [59]. However, also similar to the achiral magnet, D z does not break time-reversal symmetry, such that magnon-magnon interactions due to D z do not serve as an independent source of topology. However, note that since the kagome layers are no mirror planes in Cu(1,3-bdc), in-plane DMI components that span the entire plane are allowed [126]. An in-plane ordered ferromagnetic state has finite projection onto some of the in-plane DMI vectors but is orthogonal to others. Hence, the time-reversal breaking influence of inplane DMI enters both the harmonic as well as the cubic theory. An intricate interaction and/or counterplay is expected, influencing both magnon topology and the associated thermal Hall effect. 6 Similarly, for the experimentally considered scenario of out-of-plane magnetized Cu(1,3-bdc), D z causes harmonic topology [10,11,15,87] but the in-plane DMI vectors enter the cubic theory. Importantly, similar to the chiral model magnet in Fig. 1(b), the in-plane DMI vectors break the effective time-reversal symmetry even in the absence of D z and, hence, constitute an independent source of topology, which has been neglected so far [87,104]. With the three-magnon interactions being easily "field-frozen," a study of the magnetic-field dependence of topological magnonic gaps could clarify the relative importance of in-plane and out-of-plane DMI. Such results could shed new light on the origin of the experimentally observed sign changes in transverse transport [113].
Bilayer kagome antiferromagnet YMn 6 Sn 6 . The intermetallic room-temperature kagome antiferromagnet YMn 6 Sn 6 exhibits both gapped magnon bands and Dirac magnons [19]. With the kagome lattice admitting of both in-plane as well as out-of-plane DMI vectors (the kagome of their relative influence (S ± i = S x i ± iS y i ). Contrarily, for noncollinear and noncoplanar materials, a single type of magnetic interaction, say exchange interaction, is distributed over all sub-Hamiltonians and an analysis becomes much more complicated; a paradigmatic example is the Heisenberg antiferromanget on the triangular lattice [67]. It may very well be expected that anharmonic interactions in noncollinear magnets have an even larger influence on magnon topology and transverse transport than in collinear magnets. 6 Note that a two-dimensional ferromagnet can exhibit a Chern insulating phase and transverse transport even if the out-of-plane magnetization is zero. The only requirement is that the in-plane magnetization direction breaks all crystalline symmetries that forbid out-of-plane pointing axial vectors. Hence, upon rotating the magnetization direction by 2π within the kagome planes of Cu planes are no mirror planes), any magnetization or Néel vector direction will select a subsection of DMI components entering the harmonic theory, with the remaining components appearing in the cubic Hamiltonian. Hence, with interactions and the free theory sharing symmetries, interactions should be considered for explaining neutron scattering data. Their effect may, however, be overshadowed by the rather large electron-magnon interaction causing Landau damping due to Stoner excitations [19]. Transition metal trihalides CrBr 3 and CrI 3 . A temperature-driven many-body renormalization of Dirac magnons has been studied in the honeycomb-ferromagnet CrBr 3 [26]. The renormalization was successfully explained in terms of number-conserving four-particle interactions (Ĥ 4 ), which derive from the exchange energy [27], and, hence, do not gap out the Dirac magnons. The intrinsically small spin-orbit coupling (SOC) of CrBr 3 renders anisotropic spin-spin interactions irrelevant. However, one may think of growing CrBr 3 on a heavy-metal substrate to induce interfacial DMI for a realization of the chiral magnet model in Fig. 1 The situation is drastically different in CrI 3 , which exhibits mono-layer ferromagnetism due to SOC-induced anisotropies [129]. To explain the experimentally observed gapped Dirac magnons [18] two models-the J-D z -model and the J-K-Γmodel 7 -are debated [131], with the latter one having theoretical support [132,133]. 8 The anisotropic exchange interactions (sometimes called "Γ-terms") between nearest neighbors have a similar time-reversal symmetry breaking influence as the nearest-neighbor DMI in our chiral model magnet in Fig. 1(b). Hence, the associated three-magnon interactions may support or compete with the Kitaev terms that cause nontrivial magnon topology at the harmonic level. Furthermore, one may induce DMI by considering the intrinsically inversion asymmetric Janus monolayers [135]; recent theoretical studies have predicted very large DMI [136][137][138].
Three-dimensional pyrochlore ferromagnet Lu 2 V 2 O 7 and multiferroic ferrimagnet Cu 2 OSeO 3 . Pyrochlore ferromagnets like Lu 2 V 2 O 7 exhibit a thermal Hall effect [111,112] and were predicted to host Weyl magnons at high energies [34,35]. However, neutron scattering data revealed a rather blurred magnon spectrum [139], which may be brought about by spontaneous magnon decays due to the DMI vectors spanning the entire three-dimensional space. Hence, no matter in which direction the magnetization is pointing, DMI is distributed over harmonic as well as anharmonic Hamiltonians, rendering nonlinear spin-wave calculations necessary.
Similar considerations apply to the Weyl magnon-hosting multiferroic ferrimagnet Cu 2 OSeO 3 [39]. However, note that on top of DMI-derived three-magnon interactions there are 7 In the J-D z -model, second-nearest neighbor DMI D z causes the topological gap opening, similar to Ref. [124]. In contrast, the J-K-Γ-model relies on the Kitaev interaction K to open the gap, also topologically nontrivial [62,130]. 8 However, note also that the magnon gap observed in experiments may not be global, if the Dirac cones are shifted in reciprocal space [134].
also number-nonconserving four-magnon interactions (of type b † b † b † b) due to the ferrimagnetic spin alignment. Although they contribute to zero-temperature spontaneous decays only at order 1/S 2 , they are proportional to the exchange energy and do not suffer from the smallness of SOC. Hence, since S = 1/2, the additional 1/S "smallness" does not compensate for the four-magnon vertices being larger by a factor of J/D than the three-magnon vertices. Zero-temperature decays into the three-magnon continuum are expected to be the dominating source of damping. Nonetheless, the Hermitian part of the interaction-induced self-energy due to DMI is expected to renormalize the position of the Weyl points. Moreover, the non-Hermitian part-so far only explored in two dimensions [60]-may introduce a new topological phase of magnons which is the magnonic analog of an "exceptional topological insulator" [140].
Kamiokite  [141], which was argued to derive from phonon skew scattering (extrinsic contribution) and DMI-induced magnon-phonon interactions (intrinsic contribution) [142]; intrinsic contributions due to magnons were dismissed. The geometry between magnetization and DMI vectors is very much similar to the chiral model magnet in Fig. 1(b), and the arguments of Refs. 141 and 142 for dismissing purely magnonic contributions rely on the accidental effective time-reversal symmetry of the harmonic theory. Our results challenge this conclusion; the in-plane DMI vectors of Fe 2 Mo 3 O 8 break the effective time-reversal symmetry via three-magnon interactions and will cause a purely magnonic thermal Hall effect. Hence, a full theory of transverse thermal transport in kamiokite must consider magnons, phonons, the interactions among themselves and among each other.

VI. CONCLUSION
We demonstrated that magnon-magnon interactions may break a symmetry of the harmonic theory, admitting of a different topological phase as that found for non-interacting magnons. In particular, Dirac magnons in out-of-plane polarized honeycomb-lattice ferromagnets obtain a mass gap due to spontaneous quasiparticle decay, giving rise to a topogically nontrivial gapped phase with chiral edge states. Topological phase transitions brought about either by magnetic fields or temperature flip the chirality of the edge magnons, an experimental signature of which is a sign change in the thermal Hall conductivity.
These results show that magnon-magnon interactions do not only harbor detrimental lifetime broadening effects on the magnonic spectrum but constitute an origin of nontrivial topology on their own. Hence, they should not be neglected in studies of magnon topology. Potential material candidates to experimentally identify interaction-induced signatures either in the magnon spectrum or in transverse magnon transport are the kagome ferromagnet Cu(1,3-benzenedicarboxylate), the honeycomb-lattice van der Waals magnet CrI 3 , and the multiferroic Fe 2 Mo 3 O 8 .
Since we relied on a spin-to-boson transformation and extracted our results for the effective bosonic theory, the bosons may retrospectively be interpreted as different collective modes in solids. Hence, our general statements on the effect of interactions on single-particle topology also apply to, for example, phonons [143], triplons [144,145], bosonic spinons (Schwinger bosons) [146], magnon polarons [147,148], and magnon polaritons [149]. Overall, our findings suggest that particle-number nonconserving many-body interactions are a crucial player in the field of topology of collective excitations in quantum condensed matter systems. In Eq. (14), we defined the Holstein-Primakoff three-boson vertices. For the chiral magnet, the only nonzero vertices are k,q←p k,q←p where ϕ δ i = arg(d y δ i −id x δ i ) is the "phase" of the DMI vector belonging to the nearest-neighbor bond δ i [see Eqs. (10a)-(10c)]. Explicitly, the directions of DMI vectors read For the achiral magnet, the nonzero vertices are connect second-nearest neighbors.

Appendix B: Logarithmic singularities
In Sec. III A, we encountered a logarithmic singularity in the real part of the self-energy at the lower threshold of the two-magnon continuum. Here, we derive the analytic expression for this singularity.
We consider the kinematic situation that the Dirac energy ε D coincides with the lower threshold of the two-magnon continuum. This is the case for the critical magnetic field B c = JS . The only decay channel for which the denominator of the self-energy in Eq. (25) can be resonant is that with both decay products in the lower branch: j = j = −. The only possible energy and momentum-conserving decay channels for a Dirac magnon at K are those that involve two identical decay products at K/2 (there are three such momenta, each halfway along the high-symmetry line connecting the Γ point with the K point; see also Sec. III B 3). Hence, we expand the field B about B c , and q about the two-magnon minimum at K/2. In terms of the denominator of the self-energy becomes where a and b characterize the minimum of the two-magnon continuum along the x and y direction, respectively. As far as the singular behavior of the self-energy is concerned, we may replace the interaction vertices by where we made the dependence on DMI explicit. Here, C αβ is a constant. For the achiral magnet C αβ = 0 for α β, indicating the absence of an off-diagonal self-energy and, hence, of the gap. In contrast, for the chiral magnet C αβ 0 for all α and β. The singular part of the self-energy thus reads To cut singularities in the renormalized magnon spectrum encountered within lowest-order perturbation theory, the authors of Ref. [67] proposed a self-consistent off-shell solution of the Dyson equation. The main idea is to iteratively solve the Dyson equation, evaluating the self-energy at the so-obtained complex energiesε k , i.e., where the complex conjugateε * k accounts for causality [67]. ν denotes the iteration step. As an initial value one chooses ε (0) k = ε k − iη, with η ε k being a small numerical linewidth. Hence, after the first iteration one obtains the on-shell spectrumε (1) k . One then feedsε (1) k back into the self-energy, calculates the new energies and repeats this process untilε (ν) k converges. Sinceε (1) k contains a finite imaginary part, i.e., a finite damping Γ k , the on-shell singularities in Eqs. (B5a) and (B5b) are cut. Effectively, the finite damping smears out the threshold of the two-magnon continuum, giving rise to a finite damping also for magnons below the continuum, that is, to those magnons that are nominally stable within the Born approximation. To see so, replace i0 + in Eq. (B4) by iΓ k . Dropping C αβ , and concentrating on the singularity at δB = 0, one obtains where Λ is an artificial ultraviolet integration cutoff. In principle, one would have to solve self-consistently both for the real and the imaginary part, accounting for an energy shift of the point where the single-particle energies cut into the continuum. For the sake of simplicity, we here only solve Eq. (C2b) self-consistently for the damping Γ K = −ImΣ K (ε K + iΓ K ) in the limit Λ Γ K , resulting in Γ K ∼ D 2 /J. Thus, the damping is still perturbatively small but sufficient to cut the singularity of the real part in Eq. (C2a), Hence, the real part exhibits a non-perturbative D 2 log |D| dependence. Both the damping and the real part of the selfenergy vanish as D → 0. The divergence as Λ → ∞ in Eq. (C3) is artificial because Λ is physically bounded from above by the inverse of the lattice constant. In writing Eq. (C1), we ignored to two-band nature of the honeycomb magnet. For the chiral magnet, there will be a spectral gap at K in the renormalized spectrum for ν ≥ 1.
Hence, we decide to evaluate the new self-energy at the average energy, resulting in the equatioñ Note that Σ −− K = Σ ++ K and Σ +− K = Σ −+ K . First, we demonstrate the convergence behavior of Eq. (C4) at the critical field B = B c for S = 1 at D/J = 0.15. Figure  18(a) depicts the renormalized magnon energies as a function of the iteration step ν. Error bars, given by ±Γ K,± , indicate the damping. At ν = 0, we show the harmonic magnon spectrum for which the Dirac cone is closed (zero mass). While the on-shell spectrum (ν = 1) exhibits the logarithmic divergence as η → 0, the converged spectrumε (ν→∞) K,± is independent of the numerical linewidth. We find that ten iteration steps are sufficient to ensure convergence at B = B c . In contrast, more iteration steps (about 30) are necessary to converge results for fields slightly below the critical field, as shown in Fig. 18(b) for B/J = 0.981. This is because the singularity is not only cut but also shifted towards lower energies.
Second, we explore the magnon spectrum at the K point in dependence on B, using η/J = 10 −5 and 30 iteration steps. As shown in Fig. 18(b), the magnon spectrum of the chiral magnet is gapped, in general. For B > B c = J, the two magnon branches are well-separated and have negligible damping. However, for B < B c , decays into the two-magnon continuum cause strong lifetime broadening, as indicated by the colored areas. With the damping being larger than the gap, the notion of the gap ceases to exist. Interestingly, the gap becomes larger than the damping again for small fields. Overall, the off-shell solution for the magnon spectrum agrees very well with the spectral function in Fig. 5(f). In particular, there is no divergence at B = B c but merely a cusp.
Appendix D: Details of slab calculations for the chiral magnet Finite samples (slabs) feature boundary spins that miss neighbors. Hence, the chiral DMI is no longer compensated and there is an edge-located relaxation of the magnetic texture away from the collinear state. This phenomenon is known as "surface twist" and its amplitude decays exponentially towards the bulk [150]. To quantify this effect, we study the tilt angle ϑ of the boundary spins in dependence on DMI and magnetic field, as obtained from a classical energy minimization procedure based on an overdamped Landau-Lifshitz-Gilbert equation. The results in Fig. 19 show that for D/J = 0.2 the tilt decreases from ϑ ≈ 16 • to 6 • as the field increases from B/(JS ) ≈ 0.2 to B/(JS ) ≈ 1. Upon rotating slightly into the plane, the edge spins exhibit a finite projection onto the in-plane DM vectors, hence contributing to the harmonic theory which leads to nonreciprocal features of boundary-located magnons (ε k ε −k ), as expected for a locally broken inversion symmetry. It does not, however, lead to a band gap opening on the harmonic level, because the bulk magnetization is still pointing out of the plane. Hence, we neglect the boundary twist altogether and consider the texture to be strictly collinear. (Doing so, we also neglect the renormalizing effect of many-body interactions on the ground state directions of edge spins.) The reduced coordination of edge spins has another effect: edge excitations are shifted downwards in energy relative to bulk excitations. Colloquially speaking, since edge spins are "more floppy" than bulk spins, it does not cost as much energy to excite them. For our study, this becomes relevant when analyzing edge states. For example, a zigzag-terminated edge of the electronic honeycomb lattice is known to feature a flat state, which connects the projections of the two Dirac cones [151][152][153][154]. For magnons, however, this state is no longer flat but looks like a loose drumhead [155]. Such phenomena have been observed, e.g., for the drumhead surface states associated with magnonic nodal lines [27,40]. Although this is in principle no major road block, it increases the energy window covered by the edge states. With the slab calculations being computationally demanding, any reduction of the relevant energy window is facilitating the analysis. Hence, we apply local magnetic fields only to the edge spins. (Physically, this situation mimicks proximitizing the magnet with another magnet.) The magnitude of the local field is chosen to exactly compensate for the effective field of the missing neighbors. Hence, without interactions, the magnonic edge states of a zigzagterminated honeycomb ferromagnet become flat akin to their electronic analogs. We reiterate that the local field by itself has no influence on topology and just shifts the edges states in energy. Moreover, the edge field also decreases the effect of the surface twist even more, providing an additional reason to neglect any twists.
With the necessary approximations established, we con-tinue with considering a slab with periodic boundary conditions along the y direction but open boundary conditions in the x direction. Hence, the edges are zigzag-terminated. We choose a width of 12 honeycomb-lattice unit cells, resulting in a slab supercell of 24 spins. The harmonic Hamiltonian withÂ † k = (â k,1 , . . . ,â k,24 ) and a 24-by-24 Hamilton kernel H slab k , is readily constructed numerically by a straightforward application of linear spin-wave theory. It is diagonalized by the matrix U slab k . We also use numerics to construct the three-magnon interaction vertices between the 24 types of magnon normal modes. We rely on Refs. 62 and 71 for the general expressions of the vertices. The numerical calculation of the selfenergy Σ slab k (ε) and Green's function G slab k (ε) follow the standard procedure and are done in the eigenbasis. To facilitate numerical efforts, we only evaluate the tridiagonal entries of Σ slab k (ε) (main diagonal and the first diagonals above and below it) because those are sufficient to capture band gap openings between degenerate bands to order 1/S .
To obtain a spatially resolved spectral function where i = 1, . . . , 24 labels a site of the slab unit cell (i = 1 and i = 24 label the left-most and right-most spin, respectively), we transform the Green's function back to the Holstein-Primakoff basis. Then, we define A right k y (ε) = 24 i=13 A slab i,k y (ε) (D3b) and plot A left k y (ε)−A right k y (ε) in Fig. 8(b), choosing blue/red color for positive/negative values.