Dirac magnons in honeycomb ferromagnets

The discovery of the Dirac electron dispersion in graphene led to the question of the Dirac cone stability with respect to interactions. Coulomb interactions between electrons were shown to induce a logarithmic renormalization of the Dirac dispersion. With a rapid expansion of the list of compounds and quasiparticle bands with linear band touching, the concept of bosonic Dirac materials has emerged. We consider a specific case of ferromagnets consisting of the Van der Waals-bonded stacks of honeycomb layers, e.g chromium trihalides CrX3 (X = F, Cl, Br and I), that display two spin wave modes with energy dispersion similar to that for the electrons in graphene. At the single particle level, these materials resemble their fermionic counterparts. However, how different particle statistics and interactions affect the stability of Dirac cones has yet to be determined. To address the role of interacting Dirac magnons, we expand the theory of ferromagnets beyond the standard Dyson theory to a case of non-Bravais honeycomb layers. We demonstrate that magnon-magnon interactions lead to a significant momentum-dependent renormalization of the bare band structure in addition to strongly momentum-dependent magnon lifetimes. We show that our theory qualitatively accounts for hitherto unexplained anomalies in a nearly half century old magnetic neutron scattering data for CrBr3. We also show that honeycomb ferromagnets display dispersive surface and edge states, unlike their electronic analogs.


Introduction
The observation of fermionic quasiparticles with Dirac dispersion was a key finding for graphene [7]. Since then, the list of materials exhibiting the Dirac and Weyl energy spectra for Fermions has been extended further. Materials hosting bosonic excitations with the Dirac cones have opened a new stage in investigation of the Dirac materials such as photonic crystals [8,9], plasmonic systems [10], honeycomb arrays of superconducting grains [11] and magnets [12]. Magnets were studied in the past, especially in the form of transition metal trihalides TMX3 (where X = F, Cl, Br, I and TM=Cr), which consist of weakly coupled honeycomb ferromagnetic planes. These materials have a potential as spin polarizing elements and exhibit a strong Kerr and Faraday effects.
Early spin-wave analysis revealed a Dirac crossing point in the dispersion in the honeycomb layers containing two magnetic atoms per unit cell. Neutron scattering measurements revealed hitherto unexplained anomalies [5,6] in the boson (spin-wave) self-energies near the Dirac points.
The single-particle properties for both bosonic and fermionic Dirac materials derive from the tight-binding model on the honeycomb lattice and are, thus, identical. At the level of quantum statistics, however, there is a difference between Dirac fermions and bosons. For fermions, the excitations occur near the chemical potential, and one can focus on the low-energy Dirac cones shown in Fig. 1(a). In contrast, bosons are not subject to the Pauli exclusion principle. They can explore the entire momentum space, and excitations near zero energy dominate at low temperatures.
The importance of many-body effects was appreciated immediately for fermionic Dirac materials [13]. In particular, the Coulomb interaction between electrons, see Fig. 1(c), leads to a logarithmic renormalization of the Dirac cone velocity. This renormalization was verified by observing an anomalous dependence of the cyclotron frequency on the carrier concentration [14].
Similar analysis for the interacting bosons and their effects on the Dirac node dispersion has not been done systematically to our knowledge. Since the early paper by F. Bloch [15], the physics of magnons has remained a subject of an active experimental and theoretical research [16,17,18]. Despite an immense number of works, the important case of ferromagnets with non-Bravais lattices has received a relatively little attention. In the two milestone papers by Dyson [3,4], magnon thermodynamics is discussed only for crystals with Bravais lattices. Therefore, such theories cannot be directly applied to the honeycomb lattice, a prominent example of the non-Bravais bipartite lattices.
As mentioned we focus on the Dirac bosons, as realized by Cr trihalides with ferromagnetic honeycomb lattices. We start with a 2D honeycomb ferromagnet, where the bosons are spinwaves (magnons) that form Dirac nodes, Fig. 1. We calculate the lowest order self-energy diagrams shown in Fig. 1(d). The Hartree self-energy gives a uniform renormalization of the energy bands consistent with the theory of Bloch [19]. We then consider both the real and imaginary parts of the self-energy, which give the energy renormalization and decay rate of magnon excitations. We find that interactions induce strong temperature-and momentumdependent renormalization of the magnon bands near a Dirac node. Our results allow to explain the outstanding puzzle of spin-wave anomalies seen in CrB3[5, 6].
We also discuss the role of topology and qualitative difference between bosons and fermions for the formation of the surface states. Bosons do not have a Pauli constraint and several bosons can occupy the same state. This has important consequences for single particle properties such as surface states, and for many-body renormalization of the Dirac dispersion. Dispersive states can exist on the surface of a ferromagnetic "graphite", while they do not for fermions [20,21,22].
Our result are applicable to the rapidly growing class of Van der Waals ferromagnets with non-Bravais lattices. We note recent experimental works by Gong et al. [23] and Huang et al. [24] that demonstrate the robustness of honeycomb ferromagnetic phase in a single layer materials including CrI3. Recent theoretical works [25,26] with the focus on antiferromagnets are also relevant for our discussion.
The outline of the paper is as follows. After Introduction in Section I, we present the free spinwave theory in Section II. In Section III, we present the main result of the paperthe evaluation of the self-energies and their effect on magnon spectrum. In Section IV, we discuss properties of topological surface states and conclude in Section V.
We begin with a Heisenberg model on a 2D honeycomb lattice with two sites per unit cell commonly denoted as A and B. The Heisenberg Hamiltonian of the model describes the nearest-neighbor coupling between spins. For convenience, we choose energy units for the coupling J, whereas the magnitude of spin S is dimensionless. Positive coupling J > 0 implies a ferromagnetic ground state at low temperatures ( 32.5 K, = 3/2 for CrBr3). We choose to study the magnetic excitations, the Dirac magnons, above the preexisting ferromagnetic ground state on the honeycomb lattice. For the moment we ignore the interlayer coupling in CrBr3 that will be discussed at the end of the paper.
We follow the standard practice and bosonize the Heisenberg Hamiltonian using the Holstein-Primakoff transformation. Truncated to zeroth order, these transformations relate the spin- (1) Here, the two-by-two single-particle Hamiltonian H0 acts upon the spinor with components corresponding to the two sublattices. The off-diagonal element is ubiquitous for the honeycomb lattice and corresponds to coupling along the three nearestneighbor in-plane bond vectors . The energy spectrum of the single particle consists of two branches: the acoustic "down" and optical "up" branch with the dispersion and wave functions , where for the u (d) state and the momentum-dependent phase is introduced. The magnon energy dispersion is plotted in Fig. 1(b). The down branch touches zero energy quadratically in the center of the Brillouin zone near the point as . The gapless magnons at the point are protected by the Goldstone theorem. The magnon branches have symmetry-protected Dirac crossings at the K and K′ points of the Brillouin zone characteristic of the honeycomb lattice. At the single-particle level, the magnon dispersion law is identical to the energy spectrum of electrons in graphene.

III.
Results-bulk states The interaction vertex is obtained from the next-order terms following from the Holstein-Primakoff transformation (see Supplemental Material for details) where the momentum before and after scattering is conserved . Note that the same function defined above for the single-particle Hamiltonian (1) also occurs in the interaction term (2).
We analyze the effect of interactions to first order and evaluate the Hartree diagram in Fig.  1 where . The matrix elements of the two-by-two Hamiltonian are renormalized by the same temperature-dependent ratio . This is a consequence of a delicate balance between the bare spectrum and interaction term, both containing the honeycomb function .
So, the energy of the renormalized bands becomes . The magnon-magnon interaction leads to an overall temperature-dependent bandwidth renormalization. This is consistent with the observation by Bloch [19], who discussed a similar effect for a single-band case on a cubic lattice.
Next, we analyze the effect of interactions to second order. We consider the "rainbow" diagram in Fig. 1(d) , (4) where is the scattering matrix element, and is the area of the honeycomb lattice unit cell. Let us comment on the labeling of the momenta in the scattering process: two original magnons with momenta are scattered into the two magnons with momenta and . The momentum of the thermally excited magnons is . Note that the numerator of the diagram contains the Bose occupation number of the thermal magnons , which controls the temperature dependence of the self-energy. In the limit of small temperatures, there is a simplification: only the thermal magnons from the down band in the vicinity of the point are excited as is illustrated in Fig. 1(b). The low temperature expansion allows us to integrate out the thermal magnons analytically and reduce Eq. (4) to (5) where , and is an effective matrix element at low temperature. The remaining sum in Eq. (5) is now carried out numerically and the result is plotted in Fig. 2 for momentum along a closed path in momentum space. Note that the magnitude of the diagram is set by the overall temperature-dependent prefactor in Eq. (5). We show the self-energy separately for the up and down magnon branches in the left and right columns of Fig. 2, correspondingly. Generally, we find that the matrix element in the numerator is a well-behaved nonsingular function in momentum space that only slightly modifies the pattern of the self-energy. In contrast, the energy denominator in Eq. (5) is the major factor determining the self-energy (Fig.2) as a function of momentum. The real and imaginary parts of the self-energy have physical significance as the scattering rate and energy renormalization of the magnon excitations, respectively. First, let us comment on the scattering rate, i.e. the imaginary part of the self-energy, shown in the top row of Fig. 2. The scattering process is on-shell, i.e. the energy of the magnons is conserved. Due to the multiband nature of the Hamiltonian, there are multiple scattering channels contributing to each decay process. For example, the scattering rate for the up band shown in Fig. 2(a) is derived from the two scattering channels and . The decay rate displays a complex pattern of peaks and dips. A strong peak near the point corresponds to a mirror symmetry of the magnon bands [27]. The other peaks correspond to saddle points of the energy denominator in the self-energy (5). The decay rate of the down band, shown in Fig. 2(b), is determined by the scattering channel , and in the vicinity of the K point it is sharply suppressed, i.e. the magnon at the Dirac nodal point is a well-defined excitation whose lifetime shrinks rapidly for small excursions away from the K point.
We also plot the real part of the self-energy in the bottom row of the Fig 2: in panel (c) for the up band and in panel (d) for the down band. The real and imaginary parts obey the Kramers-Kronig relation, so a peak in the scattering rate corresponds to a feature in the real part of the self-energy. Results in lower panels in Fig. 2, closely correlate with the upper panels. We pay special attention to the point where the Dirac nodal point is located. Here the values of the real parts of the self-energy for the up and down bands are equal within the numerical precision. The interaction preserves the degeneracy of the Dirac nodal point. To estimate the effect of renormalization on the magnon energy spectrum, we compare the bare spectrum and temperature-dependent renormalized spectrum in the vicinity of the point in momentum space in Fig. 3(a). Interaction leads to renormalization of the energy of the Dirac nodal point together with a decrease of Dirac velocity in the immediate vicinity of the point. The latter behavior is in the stark contrast to that for the fermionic Dirac materials, where the Coulomb interaction, leads to a logarithmic increase of the Dirac velocity [13,14]. There are also sharp kinks, marked by the green and blue arrows, in the self-energies, which correspond to equally sharp features in the magnon decay rate and strongly reshape the original linear Dirac dispersion.
We believe that these calculations explain the old unresolved puzzle of the peculiar magnonband renormalization in CrBr3 observed in [5, 6]. In Fig. 3(b), we plot the calculated self-energy shift together with the experimental points from Ref. . We evaluate from our theory for the coupling constant = 18 K, which best fits the experimental data. We then plot with the red line and overlay it on the experimental data in Fig. 3(b). Surprisingly, the theoretical curve and the experimental data agree reasonably well. In particular, the self-energy correction in both experiment and theory undergo rapid evolution in a small region bounded by the peak, marked by the green arrow, near K.
The strong energy dependence of the self-energy is the result of in-plane interactions and Dirac spectrum. CrBr3 is a three-dimensional (3D) material composed of stacks of weakly-coupled two-dimensional (2D) honeycomb layers. To address three dimension aspects one would need to include the interlayer coupling. We expect that the interlayer coupling may lead to a weak broadening of peaks in the self energy shown in Fig. 2, but the magnitude of the effect is small since the interlayer coupling Jz is much smaller than the intralayer coupling Jz <<J [6].
The results obtained in this work are universal and apply to other types of bosonic Dirac materials. The profile of in momentum space in Fig. 2 is determined predominantly by energy denominators and less so by the matrix elements in Eqs. (4) and (5). So, the momentum dependence of these quantities could be viewed as a unique signature of the honeycomb dispersion. Thus, our findings should be highly relevant to other bosonic honeycomb systems, e.g. photonic crystals [8, 9,30,31,32,33,34].

IV.
Results -Surface states in a 3D model.
Dirac nodes have profound consequences for the topological surface states that are analogous to Fermi nodal lines. Because the honeycomb lattice is bipartite, it is a primary candidate to possess topological Shockley edge states (otherwise known as the Su-Schrieffer-Heeger states). Crystal structure of CrBr3 shown in Fig.4(a) shows honeycomb Cr layers stacked in the order that is an analogue of the ABC-stacked graphene layers, which have received a lot of attention because of an unusual spectrum and edge states [20,21,22]. We assume the simplest model in which the honeycomb layers are coupled via the vertical bonds with a strength as illustrated in Fig. 4(a). So, the corresponding free magnon Hamiltonian becomes By comparing Eq. (7) with Eq. (1) one observes that the matrix elements acquire extra terms proportional to . Let us analyze the fate of the Dirac point given by the equation . For small , the solution gives helical contours winding around the corners of the Brillouin zone as shown in Fig. 4(b). This nodal line, i.e. a line where the magnon branches are degenerate, retains the topological properties of the usual Dirac point. In particular, the nodal line has a topological Berry phase which leads to flat band surface states (in the fermionic case) given by the projection of the nodal line onto the 2D momentum space [20,21,22] as illustrated in Fig. 4(b).
Now we evaluate the surface states in the discussed 3D model. Although the surface states are calculated at the single-particle level, it turns out that the surface magnons differ from their fermionic analogues. We compute the magnon dispersion in a slab geometry of N=25 layers and plot the result in Fig. 4(c). The figure contains the spectrum for the closed (in black) and open (in red) boundary conditions. In this representation, the red lines unmatched by the black lines correspond to the surface states.
One can observe the magnon surface states appearing in the vicinity of the K point, which agrees with Fig. 4(b). However, in contrast to what occurs for fermions on stacked honeycomb lattices, the magnon surface states are clearly dispersive. To explain this, we recall the original Heisenberg Hamiltonian and the Holstein-Primakoff transformation, e.g.
for the z-th component of the spin operators . Observe that inter-site coupling in the Heisenberg Hamiltonian generates an on-site potential in the bosonic language, i.e.
( ). For the bulk, this procedure generates the diagonal terms in the Hamiltonians (1) and (7). For the surface, inspection of Fig. 4a shows that one sublatticewe assume it is A -will be coupled (to the B sublattice) in the layer immediately below, but not the other. This generates an extra diagonal surface term for the A but not the B sublattice, in addition to the off-diagonal terms which allow hopping from the A sublattice of the surface to the B sublattice. In the absence of the diagonal surface term, we obtain dispersionless topological surface states that are identical to the fermionic analogues, see e.g. Ref. [22]. Note that the decay length of this surface state in the z direction is inversely proportional to the bulk gap , where is a good quantum number for a surface cleaved parallel to the layers.
Qualitatively, the dispersion of the magnon surface states can then be understood by treating the surface term as a perturbation. Already the first-order perturbation theory produces a kdependent dispersion of the surface states . Indeed, the nonperturbative exact diagonalization of the bosonic Hamiltonian together with the surface term confirms that the surface state acquires a strong dispersion (as shown in the Supplemental Material). We further relax the interplanar Heisenberg coupling Jz as the surface is approached; Fig. 4(c) is plotted for a reduced coupling between the surface and penultimate layers. Our discussion extends the previous analysis of the edge states in 2D photonic materials. We also mention that the described mechanism for generating the dispersion of the surface states is universal and should be of relevance to other types of bosonic Dirac materials, e.g. topological phononic and photonic systems [8,9,30,31,32,33,34,35].

V. Concluding remarks
Noninteracting particles in a 2D honeycomb lattice exhibit Dirac excitations regardless of the statistics. We consider the question of difference between bosonic and fermionic excitations on such lattices. Two key features of the bosonic Dirac material were identified as important: i) Interacting bosons lead to a non-divergent velocity renormalization near the Dirac points. In contrast, Coulomb interactions between fermions lead to a logarithmic correction of the Dirac velocity. ii) non-Bravais nature of the honeycomb lattice structure leads to significant modifications of the spin-wave interactions. More generally, non-Bravais lattices, which entail multicomponent wave-functions for quasiparticles, provide a route to Dirac quasiparticles and distinct surface states. iii) we found topological surface states in CrBr3 that are the bosonic analog of Shockley edge and surface states (also known as the Su-Schrieffer-Heeger states).
To prove the point, we consider the case of ferromagnetic Cr trihalides. Magnetic excitations in honeycomb ferromagnetic layers, the essential element of magnetism in Cr trihalides, possess Dirac nodes. Half a century ago, neutron scattering experiments observed an anomalous momentum-dependent renormalization of the magnon spectrum in CrBr3 near Dirac points [4,5]. By evaluating the self-energy due to the magnon-magnon interactions, we obtained a good agreement with data from Refs. [4,5] and, thus, resolved a fifty-year-old puzzle.
Given the large diversity of transition metal trihalides [29,36], metamaterials and cold atom systems which can implement bosonic Hamiltonians for honeycomb and other non-Bravais lattices, our predictions are highly relevant and can be extensively checked using modern spectroscopic methods. For example, neutron scattering has advanced greatly over the last five decades, and so we anticipate new work to measure the lifetimes as well as self-energies of spin waves for honeycomb magnets. There are also potential applications of the surface states, which we have discovered, to the emerging area of "magnon spintronics" or "magnonics" [37], for which distinct surface spin waves would be advantageous because they avoid bulk dissipation.
Recent experimental observation of the single layer honeycomb ferromagnets [23,24] in CrI3 and other materials highlights the important role the Dirac nodes will play in bosonic Dirac materials realized in Van der Waals magnetic structures.

A. Bosonic Hamiltonian
We follow the standard procedure and bosonize the ferromagnetic Heisenberg Hamiltonian on the honeycomb lattice , ( ).
We introduce the bosonic annihilation operators and corresponding to the two sublattices A and B. We relate the spin and boson operators using the Holstein-Primakoff transformation truncated to the first order in 1/ as follows (similarly for the B-sublattice). We substitute the leading order terms from the Holstein-Primakoff Eq. (10) in the Hamiltonian (9) (12) correspond to the optical and acoustic magnon branches that we also refer to as up and down  k k k k k k k k k k k k k k k k k k k k k k k k k k k (14) Observe that each interaction term consists of the two creation and annihilation operators, which manifestly conserve the magnon number. The four magnon interaction acts upon the states with at least two excited magnons. Thus, in the limit of zero temperature 0 T  , where no magnons are present, the interaction has no effect. In contrast for finite temperatures, the interaction effects should become more pronounced. However, in the case of the canted antiferromagnet, such a rotational symmetry is absent, which permits the three-magnon terms in the interaction Hamiltonian. The effect of the three-magnon process is non-zero even at absolute zero 0 T  . Thus, although the single-particle spectra may be nominally identical in the cases of the honeycomb ferromagnet and antiferromagnet, the effects of interaction terms are distinct. In what follows, we will describe this finite temperature interaction effects.

B.1. Hartree term:
Here we explain the Hartree energy renormalization in Eq.
(3) of the main text according to the diagram shown in Fig. 1(a) . For small temperatures T J , only the low-energy down -magnons with momenta close to 0 q  are excited. In order to calculate the self-energy, we replace the product of two operators by their expectation values. There are four ways to accomplish this: a a  a b a a  a b a a  b a a a  b a So, the Hartree contribution can be rewritten as, (17) where for brevity we define the coefficients,  (19) written in the sublattice basis. We compare the free Hamiltonian (11) with the Hartree selfenergy (19) (11). We observe that the Hartree term leads to a temperature-dependent rescaling of both the branches of the energy spectrum uniformly throughout the Brillouin zone.
B.2. Rainbow diagram: Next, we analyze the effect of interaction to the next order in 1 S . The self-energy corresponding to the diagram in Fig. 1(b)  . However, we choose a simplifying route and notice that Eq. (21) can be further reduced in the limit of small temperatures T J . We observe that the first term in the thermodynamic function , ; F k q p corresponds to the direct scattering process, i.e. where the original magnon with momentum k and the thermal magnon with momentum q scatter into magnons with momenta p and  k q p , whereas the second term to the reverse scattering process. For small temperatures T J , the direct scattering process dominates, and ,; F k q p can be approximated as ,; () d