Magnonic Weyl states in Cu2OSeO3

The multiferroic ferrimagnet Cu$_2$OSeO$_3$ with a chiral crystal structure attracted a lot of recent attention due to the emergence of magnetic skyrmion order in this material. Here, the topological properties of its magnon excitations are systematically investigated by linear spin-wave theory and inelastic neutron scattering. When considering Heisenberg exchange interactions only, two degenerate Weyl magnon nodes with topological charges $\pm$2 are observed at high-symmetry points. Each Weyl point splits into two as the symmetry of the system is further reduced by including into consideration the nearest-neighbor Dzyaloshinsky-Moriya interaction, crucial for obtaining an accurate fit to the experimental spin-wave spectrum. The predicted topological properties are verified by surface state and Chern number analysis. Additionally, we predict that a measurable thermal Hall conductivity can be associated with the emergence of the Weyl points, the position of which can be tuned by changing the crystal symmetry of the material.

Experimentally, topologically nontrivial magnon states were recently identified in a two-dimensional spin- 1 2 kagomelattice ferromagnet [11] and in the three-dimensional (3D) antiferromagnet Cu 3 TeO 6 by two independent groups [12,13] using inelastic neutron scattering (INS). On the theory side, it has been realized that chiral magnets offer a generic route to the realization of topological magnon states, representing a magnon analog of topological insulators. As a result of antisymmetric exchange, known as Dzyaloshinsky-Moriya interaction (DMI) [28,29], the bulk magnon spectrum of a chiral magnet can acquire a topological energy gap that supports a topologically protected gapless Dirac cone in the surface magnon spectrum [17]. A similar mechanism based on DMI was also proposed for the formation of magnonic Weyl crossing points in the spin-wave spectrum of the noncoplanar antiferromagnetic (AFM) state on a breathing-pyrochlore lattice [19,20,22].
The cubic copper(II)-oxoselenite Cu 2 OSeO 3 is a multiferroic ferrimagnet with a chiral crystal structure that came under the focus of recent attention owing to the emergence of skyrmion order in this material [30][31][32][33][34]. Its crystal structure is cubic (space group P2 1 3) with the lattice constant a = 8.925 Å [35]. The magnetic sublattice of Cu 2+ ions can be approximated as a distorted breathing-pyrochlore lattice, consisting of slightly deformed tetrahedral Cu 4 clusters in a face-centered cubic (fcc) arrangement [36]. Magnetic interactions within the tetrahedron lead to a ferrimagnetic ground state, in which one of the Cu 2+ spins is antiparallel to the other three, resulting in the total spin S = 1 of the cluster [37][38][39]. Weaker interactions between the clusters lead to a long-range spin-spiral order that sets in below T C ≈ 57 K. Existing magnetic models [37][38][39][40] consider up to 5 Heisenberg exchange interactions and up to 5 DMI vectors. These models were used to describe the INS spectrum of spin-wave excitations in a broad energy range and in the whole reciprocal space [36], as well as electron spin resonance (ESR) that probes spin-wave excitations at the zone center [41]. However, the DMI was initially neglected in these studies.
This simplified description, that involves only Heisenberg interactions, provides a qualitatively good fit to the experimental spin-wave dispersion over the entire Brillouin zone [36] with the exception of the zone corner (R point), where the magnon bands remain degenerate for any values of the exchange parameters. Tucker et al. [42] recently showed that this degeneracy is removed by DMI, leading to a clearly resolved spin gap of ∼1.6 meV in the magnon spectrum, which they observed by neutron spectroscopy. These observations are a strong indication for the existence of topological magnon states in Cu 2 OSeO 3 , which motivated our present study.
In the following, we present spin-dynamical calculations of the magnon spectrum in the presence of DMI that was adjusted to provide the best fit to the experimental spin-wave dispersion in the vicinity of the R point. Using linear spinwave theory (LSWT) in combination with high-resolution neutron spectroscopy, we show that, in the absence of DMI terms, two pairs of degenerate Weyl nodes with the topo-logical charge +2 and −2 are located at the zone center (Γ point) and at the zone boundary (R point). Consideration of the nearest-neighbor DMI is sufficient to lift the degeneracy of these Weyl nodes, so that they are shifted away from the high-symmetry points into a position that sensitively depends on the direction and magnitude of the DMI vector. A direct observation of the resulting Weyl points would offer a possibility to accurately extract the DMI from INS measurements. We verify the predicted topological properties by the Chern number analysis and give quantitative predictions for the location of magnonic Weyl points in the spin-wave spectrum. We also analyze topologically protected magnon surface states and estimate the magnetic contribution to the thermal Hall conductivity that may serve as robust hallmarks of the emergent topological states in Cu 2 OSeO 3 , awaiting a direct experimental verification.

A. Magnetic model, experimental result and magnon spectrum
The crystal structure of Cu 2 OSeO 3 belongs to the chiral space group P2 1 3 and contains 16 magnetic Cu 2+ ions per unit cell with S = 1 2 . They occupy two structurally nonequivalent positions, so that every Cu 4 tetrahedral cluster consists of one Cu(1) ion on the 4a Wyckoff site and three Cu (2) ions on the 12b site [35,43], see Fig. 1. The strong superexchange coupling J FM s between the Cu(2) ions within the cluster is ferromagnetic (FM), whereas the Cu(1) and Cu(2) spins within the same tetrahedron are coupled antiferromagnetically with a coupling constant J AFM s . These exchange constants constitute the dominant magnetic interactions that lead to a ferrimagnetic spin arrangement within the cluster: Three Cu(2) spins align ferromagnetically, and the Cu(1) spin is pointing in the opposite direction, resulting in a total spin of S = 1 [44]. The intercluster interactions are considerably weaker, given by the FM superexchange J FM w between the nearest Cu(2) ions of neighboring clusters, the weak AFM coupling J AFM w between Cu(1) and Cu(2), and a longer-range exchange J AFM O..O that connects Cu(1) and Cu(2) sites across the diagonals of alternating Cu(1)-Cu(2) hexagon loops [38], see Fig. 1.
Numerical values of all five Heisenberg interactions have been calculated from the microscopic electronic structure theory and verified using thermodynamic data [39], terahertz ESR [41], far-infrared [40] and Raman [45] spectroscopy, and INS measurements [36,42] in earlier works. As a result, there are accurate quantitative estimates of all five exchange parameters. On the other hand, antisymmetric DMI is also allowed by crystal symmetry along all mentioned exchange paths. Each DMI channel can add at most three extra parameters, which are the off-diagonal components of theĴ tensor. This results in up to 15 additional parameters in the magnetic Hamiltonian. All previously reported attempts to estimate their magnitude are based on first principles calculations [39]. Hence, the measurable DMI-signatures in the spin-wave spectrum [42] still await experimental verification.
Here we use LSWT to calculate the magnon spectrum of Cu 2 OSeO 3 , starting from the generalized Heisenberg model, where the interaction tensor between the lattice sites i and ĵ includes the symmetric exchange J i j and the antisymmetric off-diagonal DMI terms D i j , caused by the spin-orbit coupling. The DMI vector is defined as . Following earlier works [36], we include five Heisenberg exchange interactions shown in Fig. 1(a), with their numerical values listed in Table II. To deal with the ferrimagnetic system, the rotation matrix O i is introduced, where O i determines the magnetic moment direction at the site i.
To get the magnons excitation spectrum, the LSWT is used [21,[46][47][48], where the Holstein-Primakoff transformation [49] is adopted for the quantum spin operators. A Fourier transformation of the boson operators is given by where N is the number of the unit cells, and k is the vector in the reciprocal k-space of magnons (suppressed below for convenience). The fourier-transformed Hamiltonian part quadratic in the Boson operators, denoted as H 2 , becomes a 2n × 2n matrix, where n is the number of the atoms in the unit cell. From the commutation relation between the bosonic creation (annihilation) operators and H 2 , we arrive at the equation where the dynamical matrix is given by D = gH 2 with g = [( , 0), (0, − )], where is the n × n identity matrix. The positive real eigenvalues of the dynamical matrix D correspond to the magnons excitation spectrum in the system. While the left and right eigenvectors of D, denoted as V L and V R , may differ since D is not necessarily hermitian, their relation is trivially given by V L = gV † R g. The INS experiment has been carried out on the coldneutron triple-axis spectrometer PANDA [50] located at MLZ in Garching, Germany. The sample is a coaligned mosaic of 11 single crystals with a total mass of ∼ 2 grams and a mosaicity of ∼ 2 • . It was mounted in the (H H L) scattering plane, i.e. with the [110] axis vertical, inside the JVM1-5.0T cryomagnet with a base temperature of 1.5 K. The instrument was operated with a fixed final neutron wavelength k f = 1.5 Å −1 . To avoid higher-order contamination from the monochromator, a cold beryllium filter was mounted between the sample and the analyzer.
The magnon spectrum of Cu 2 OSeO 3 has been analyzed in several previous works, which used similar values of the exchange parameters J and neglected the effect of the DMI [36,38,42]. In these works, similar magnon dispersions were obtained, featuring two doubly-degenerate crossing points: one at Γ and one at R high-symmetry points. The results of our calculations, performed without DMI, are very close to previously published data, and they are shown in Fig. 1 (d). On the other hand, the results of our highresolution INS measurements, presented in Fig. 2, clearly mark the formation of a 1.6 meV band gap at R between the bands 2 and 3, according to the enumeration of Fig. 1 (d). To reproduce this band gap in the calculations, we have chosen the value of the nearest-neighbor DMI by fitting it to the experimental data. The values of fitted J and D parameters are listed in Table II, whereas Fig. 2 shows the comparison of the experimental and calculated magnon spectra. The value of the DMI that we use in this work provided by far the best fit to the experiment among other possible DMI choices, which e.g. included more neighbors into considerations, or even with respect to previously published ab-initio results for the DMI in this system [39], see Figs. S2 and S3 of the Supplemental Material [51]. Irrespective of its exact choice, including the DMI into the picture has a drastic effect on the number and position of the degenerate crossings between bands 2 and 3 in the magnonic band structure. The set of DMI parameters we used here (see Table II ously degenerate crossing points at Γ and R, giving rise to overall four crossings: two at R and R (in the vicinity of R), and two at Γ and Γ (in the vicinity of Γ), see Fig. 1 (d,e).
In the next subsection we analyze the topological character of these points.

B. Topological properties
In order to access the topological properties of the system we calculate the Berry curvature of each magnonic band n, defined as [20,48] where ε mk are the magnonic eigenvalues. As we mainly focus on the topological nature of the band crossings arising between the bands 2 and 3, we analyze the cumulative Berry curvature of the bands 1 and 2. In Fig. 1 (b) we present the direction of the normalized projected cumulative Berry curvature vector field and its absolute magnitude in the k y = k z plane, first for the case without DMI. In the mentioned figure, the color scale represents the absolute value of the Berry curvature vector field. As apparent from the figure, the Berry curvature distribution exhibits two monopole-like features at R and Γ, where the band crossings occur, with the crossing at Γ serving as a source, and the crossing at R serving as a sink of the Berry curvature field. The corresponding distribution, obtained after including the DMI, is shown in 2 ) point, where degenerate magnon bands are predicted in the absence of DMI. All measurements were done at 2 K without applying magnetic field. The magnon dispersion calculated with fitted value of the DMI is shown with thin red lines. The straight feature in panel (c), which is not captured by the spin-wave model, is a spurious feature from nonmagnetic multiple scattering. Fig. 1 (c) in the plane which includes Γ and R points and which is perpendicular to the k y -k z plane. In the latter case the distribution of the Berry curvature field, although similar to the previous case, is more complex, owing to the fact that the crossings at Γ and R are very close to the plane so that the overall distribution coming from all four points is plotted.
Next, we compute the monopole charge of the i'th band crossing by evaluating the flux of the cumulative Berry curvature field through an infinitesimal two-dimensional sphere S i surrounding the crossing: where n is the surface normal. According to our calculations, without the DMI, the total topological charge of the two degenerate points at Γ is +2, while it constitutes a value of −2 at R. Upon including the effect of the DMI, each of the double degeneracies splits into two nondegenerate points with charges of +1 at Γ and Γ , and −1 at R and R . The topological analysis is further supported by the Brillouin zone evolution of the first Chern number, defined analogously to the charge as: where P is a two-dimensional slice of the Brillouin zone and n is its normal. By defining the plane P as the k x -k y plane at a given k z with n = (0, 0, 1), we compute the evolution of C(k z ) as a function of k z , presenting the results in Fig. S5 of the Supplemental Material [51]. Without DMI, the Chern number changes by 2 when P passes through the degenerate crossing points, while in the presence of DMI it changes by 1 when P passes through every nondegenerate crossing point. This analysis underlines the main finding of our paper -the emergence of two doubly degenerate type-I Weyl points [52] in the magnonic structure of Cu 2 OSeO 3 , located at R and Γ without DMI, which further split into overall four Weyl points when the symmetry of the system is reduced by including the DMI into consideration.

C. Surface states
As the emergence of the Weyl points in the magnonic band structure of a three-dimensional crystal is expected to give rise to the surface states of a thin film, here, we analyze the magnon band structure of a 75-layer thick twodimensional slab of Cu 2 OSeO 3 cut along the [001] axis, presenting the results in Fig. 3. The spin-wave dispersion is shown along the path which includes the projections of the Weyl points onto (001)-plane, which are further indicated with red and blue small circles in the figure, according to their topological charge. In the magnon band structure, the states are marked with their weight at the surface of the slab (see Supplemental Material [51] for more details). The left plot in Fig. 3 corresponds to the situation without the DMI, with projections of the Weyl points positioned at high symmetry points in the two-dimensional Brillouin zone. We observe that in this case the Weyl points of opposite chirality are connected by the magnon "arc" surface states, which is in accord with our topological analysis from above.
where n enumerates the magnon bands, f B n is the Bose-Einstein distribution function, which can be expressed as f B n = (e ε nk /k B T − 1) −1 , and C 2 is given by with Li 2 denoting the dilogarithm function. The thermal Hall conductivity tensor of the system is then defined as From experiment we know that the Curie temperature of Cu 2 OSeO 3 is around 60 K [39,42]. The computed energydependence and the cumulative components of the thermal Hall conductivity, calculated according to equations above at 60 K, are shown in Figs. 4 (a,b). In these plots we observe that in the energy region between 9 and 10 meV there is a significant enhancement especially in the κ xz component of the thermal Hall conductivity. This enhancement can be attributed to the distribution of the Berry curvature around the Weyl points in that energy region, which correspondingly gives rise to the fingerprint of the Weyl points in the energy is the consequence of the vanishing contribution by the "topologically-trivial" low-lying bands which are basically not affected by the DMI, Fig. 1(d). Respectively, the thermal Hall effect "lifts off" once the region of Weyl points is reached by the distribution of magnons. The overall magnitude of the thermal Hall effect that we predict in the region of higher temperatures is large enough to be observed in experiment.

E. Effect of the DMI on the position of Weyl points
Given the low structural symmetry of Cu 2 OSeO 3 , it seems reasonable to explore the influence of the direction and strength of the DMI vector on the position of the Weyl points in the Brillouin zone. While we envisage that the tuning of the DMI parameters can be realized e.g. by pressure, strain, electric field [55][56][57][58], or doping with defects, knowing the correlation between the Weyl point geometry and the DMI provides a unique tool for accessing the details of the DMI in a given sample, which are challenging to extract with other techniques based e.g. on measuring the properties of domain walls [59,60].
To estimate the influence of the DMI on the Weyl points, we first keep the direction of the DMI along the [110] axis, while scaling its magnitude between 0 and 1.5 meV (color scale in Fig. 5 (d)). The evolution of the Weyl points around R upon increasing the DMI is shown in Fig. 5 (d). Notably, upon starting from a degenerate case at zero DMI, the splitting between the two Weyl point is clearly driven by lowering of symmetry upon including the non-vanishing DMI, while the trajectories of the Weyl points in the Brillouin zone are almost perfectly straight lines.
Further, after fixing the magnitude of the DMI to the value of 1meV, we rotate the direction of the nearest neighbor DMI vector, as specified by angle θ in Fig. 5(b), about the x-axis, and track the position of two Weyl points around Γ and R in Fig. 5(c) and (a), respectively. The results indicate that the Weyl points rotate around the R and Γ points along specific paths when following the rotation of the DMI vector. The corresponding trajectories, while having a relatively complex shape in the three-dimensional Brillouin zone, clearly possess a high degree of symmetry, as apparent from the projections of the trajectories onto the high-symmetry planes, see e.g. Fig. 5(a) and (c). We show further data on the correlation between the DMI and the Weyl point behavior in the Supplemental Material [51].

III. DISCUSSION
In our work, based on the spin-wave theory and experiment, we arrived at several important findings concerning the spin-wave properties of ferrimagnetic Cu 2 OSeO 3 . Firstly, we were able to attribute the origin of the experimentally observed magnon band gap in the spin-wave spectrum at the R point to the effect of the DMI, which was chosen so as to provide the best fit to the high-resolution neutron scattering data. Secondly, after systematically addressing the topological properties of Cu 2 OSeO 3 , we uncovered the emergence of the doubly-degenerate Weyl nodes with topological charge ±2 at high-symmetry points even without the effect of the DMI. We further observed that each Weyl point splits into two as the symmetry of the system is reduced when introducing the DMI into play. Importantly, we find that the position of the Weyl points can be controlled by changing the crystal symmetry of the compound. We further predict that the emergence of the Weyl points in the system goes hand in hand with the formation of topological magnonic surface states, which can be observed for instance at the (001) surface of Cu 2 OSeO 3 .
Our findings open a quest for experimental observation of the Weyl points in this material, and exploring the influence of such points in the spin-wave spectrum on various properties of more complex magnetic phases in Cu 2 OSeO 3 , for example, its skyrmion phase. While we discover that Weyl points play a crucial role in shaping the magnitude and temperature dependence of the thermal Hall effect in its ferrimagnetic phase, we expect that the same holds true also for skyrmions in Cu 2 OSeO 3 . The observation of the exact position of the Weyl points as well as following their dynamics upon structural reconstructions in Cu 2 OSeO 3 can further provide a unique tool for accessing the microscopics of the DMI in this complex compound, which can be of paramount importance for understanding and shaping of chiral dynamics and properties of Cu 2 OSeO 3 . The latter finding also suggests that in special materials of Cu 2 OSeO 3 type one can expect that the topologies in the space of magnons and in the real-space of skyrmions can be closely intertwined.

IV. ACKNOWLEDGEMENTS
We acknowledge fruitful discussions with Flaviano José   Table SII.

II. RESULTS WITH AB-INITIO PARAMETERS
In this paper DMI are introduced to enhance the agreement between the experimental magnon dispersion and the LSWT calculations. While the nearest neighbour DMI of the main text was obtained through fitting to the measured dispersion, this section utilizes 5 distinct ab-initio DMI which were reported by O. Janson et al. [Nat. Commun. 5, 5376 (2014)]. This study also supplied the Heisenberg interaction parameters used in the main text with the spin up moment 0.450 and spin down moment 0.483, which is very similar to the ab-initio results. The resulting magnon dispersion on top of the experimental data is shown in Fig. S2.
Using the ab-initio DMI parameters, four Dirac points were found in the first Brillouin zone (BZ) at the follow-

III. FITTING TO EXPERIMENTAL SPECTRUM
Alternatively, the DMI vectors can be obtained by fitting to the experimental magnon dispersion. This section describes the procedure as applied in the main text and demonstrates analogous results of further interactions. First, 13 representative points of experimental data set were selected to comprise the fitting measure. Subsequently, the Broyden- Fletcher-Goldfarb-Shanno (BFGS) Hessian update strategy is employed to converge the considered DMI vector. This procedure is executed separately for each DMI vector of the first 4 nearest neighbors given in Table 1 (Fig. S4) reveals great accuracy when considering the first nearest neighbor DMI but uncovers poor agreement in all other cases. This justifies the restriction to the first DMI in the main text.

IV. CHERN NUMBER EVOLUTION
The Chern number was calculated in order to analyze the topological character of the Weyl points. Based on Eq. 7 of the main text, the Chern number sum of bands 1 and 2 was computed for Brillouin zone slices perpendicular to k z . The

V. SURFACE LOCALIZATION AND THE SURFACE ARC
The color scale in Fig. 3 of the main text is calculated based on the following equation: where k is the reciprocal space vector, j denotes the band index, i numbers the magnetic atom, and R i z represents the normalized position for atom i along the z-axis. V i L (k, j) and V i R (k, j) are the components of the left and right eigenstates of j at the magnetic atom i.
The above expression provides a reasonable measure to judge the surface character of each state while distinguishing both sides. However, in the left plot of Fig. 3 (main text) surface states at about 9.2 meV along Γ Y and Y R appear to loose and regain their surface character without contact to bulk states, which is highly unusual. This issue is resolved by Fig. 6 which shows the real-space decomposition of one of these apparent bulk states at the Y point. That state exhibits highly localized contributions on both surfaces but none in the bulk, hence the degenerate states are indeed highly localized at the surfaces. Accordingly, Eq. 10 is unable to classify such states equally localized on both surfaces. This justifies the assumption that the other apparent bulk states along that high symmetry line are of surface character as well. Numerically, the origin of the surface character concealment is the exact energetic degeneracy which causes unsuitable eigenstate superpositions. Fig. S7 shows the distribution of the surface states in the surface Brillouin zone at the energy of 9.25 meV. The blue and red represent the states at the left and right surfaces of a 20-layer thick slab, while the grey color indicates bulk states. The plot visualizes clear arcs connecting the projected Weyl point to bulk states which enclose the Weyl points of opposite character.