Introduction

In solid states, band structures provide key information regarding the dynamics of quantized quasiparticles and collective excitations such as electrons or phonons1. For instance, insulators and conductors are classified based on whether the electronic band structures near the Fermi levels are gapped or not. In certain materials, band structures may be gapped (and insulating) only in bulk lattices while remaining gapless (and conducting) along the two-dimensional (2D) edges or three-dimensional (3D) surfaces2. These so-called topological insulators have already been intensively studied as they provide hope for dissipation-free topological edge transports3,4. Investigations in topological insulators have already led to predictions and experimental searches for analogous phenomena in photonics5, phonons6, and magnons7. The necessary ingredients are internal mechanisms to break the time-reversal symmetry of the bulk states, resulting in gapless transports confined to the topological edges8,9. In both electron and magnon band structures, which are charge and spin counterparts of electrons, spin-orbit couplings (SOC) commonly play the most important role in establishing such topological effects10.

Topological magnons are found in magnetic crystals involving hexagonal loops such as 2D kagome lattices11,12,13 or 3D pyrochlores14,15. Among them, the minimal band structures with two magnon modes are available in 2D honeycomb lattices depicted in Fig. 1a16,17,18. Similarly to the electronic band structure of graphene, these honeycomb ferromagnets exhibit crossings between two linearly dispersive magnon bands at the Dirac wave vectors \(\textbf{Q}_{\textrm{K}/\textrm{K}^\prime } = \pm (1/3,1/3)\). The associated magnons propagate effectively as relativistic massless particles and are thus called Dirac magnons in reference to the Dirac electrons in graphene16. These Dirac magnons may also become massive and topological when SOC-driven exchanges break the time-reversal symmetry and open the gap between the otherwise twofold degenerate modes19,20. Recent inelastic neutron scattering (INS) experiments observed finite gap openings larger than 1 meV at the bulk magnon Dirac points of van der Waals honeycomb ferromagnets \(\hbox {CrI}_3\) and \(\hbox {CrGeTe}_3\), suggesting that their edge states should exhibit gapless topological excitations21,22,23,24. Breaking of time-reversal symmetry in these compounds was attributed to antisymmetric Dzyaloshinskii-Moriya (DM) exchanges between the next-nearest neighbor (NNN) \(\hbox {Cr}^{3+}\) spins mediated by heavy ligands such as I or Te17,18. In contrast, \(\hbox {CrBr}_3\) with lighter Br exhibits gapless Dirac magnons, suggesting that the strength of DM exchange depends on the mass of the ligand ions25.

Fig. 1
figure 1

(a) A planar view of the honeycomb lattice with Ti ions omitted for simplicity, and (b) the antiferromagnetic structure of \(\hbox {FeTiO}_3\) showing \(\hbox {Fe}^{2+}\) ions only. The hexagonal unit cell is marked as thick solid lines. Also shown in both figures are the exchange interactions considered in this work. \(J^\parallel\) and \(J^\perp\) are the in-plane bond-directional exchange anisotropies discussed in the main text. (c) Energy level diagram of \(\hbox {Fe}^{2+}\) in the trigonally distorted octahedral crystal field. (d) Two dimensional Brillouin zone boundaries of the honeycomb lattice. Thick lines mark the slices displayed as \(I(\textbf{Q},\hbar \omega )\) plots in Figs. 3 and 4.

Although the symmetry of the honeycomb lattice does not allow the nearest-neighbor (NN) bonds to have finite DM exchanges26,27, symmetric and anisotropic exchanges between them can also break the time reversal symmetry. Bond-dependent Kitaev exchange is a good example that was proposed as an alternative mechanism to open the topological magnon gap and even induce exotic quantum spin liquid ground states in Cr-based honeycomb magnets28,29,30,31,32,33. Although widely successful in explaining the \(J_\textrm{eff}=1/2\) physics of 4d electrons34,35,36, its strength in high-S magnets with quenched orbital moments is estimated to be rather small37,38. It was eventually ruled out in \(\hbox {CrI}_3\) by INS experiments under external magnetic fields22. However, non-Kitaev bond-dependent exchanges can still be relevant in honeycomb magnets with sufficiently strong SOC39,40. For example, bond-directional anisotropic exchanges have recently been used to explain Dirac magnons and nodal line dispersions observed in \(\hbox {CoTiO}_3\)41. Since the \(\hbox {Co}^{2+}\) spins are within the honeycomb planes, however, Dirac magnons in \(\hbox {CoTiO}_3\) remain gapless similarly to \(\hbox {CrCl}_3\)42,43,44. To break the time-reversal symmetry and open the Dirac gaps, it is necessary for the magnetic moments to tilt out of the honeycomb planes.

\(\hbox {FeTiO}_3\) is isostructural to \(\hbox {CoTiO}_3\) while replacing its quantum spins with large spins of \(\hbox {Fe}^{2+}\) (\(3d^4\), \(S=2\)). Although these two oxides also share the same antiferromagnetic (AFM) ordering with stacked ferromagnetic honeycomb layers, \(\hbox {FeTiO}_3\) has spins oriented nearly perpendicular to the honeycomb planes [see Fig. 1b] contrary to \(\hbox {CoTiO}_3\). Therefore, it is considered likely to exhibit a magnon gap opening at the Dirac point. Whereas \(\hbox {Fe}^{2+}\)-based topological electronic band structure has been suggested in \(\hbox {BaFe}_2\hbox {(PO}_4\hbox {)}_2\)45, the possibility of topological magnons was not investigated. Since the \(\hbox {FeO}_6\) octahedra in honeycomb oxides are trigonally distorted and 3d orbitals are only partially filled, a strong SOC is expected [see Fig. 1c]. Using INS on single-crystal \(\hbox {FeTiO}_3\) in this work, we observe that its magnon band indeed exhibits a finite gap of \(\Delta _\textrm{K} = 1.2\) meV at its Dirac points.

Results

Elastic and inelastic neutron scattering

The ordering wave vector of the A-type AFM structure depicted in Fig. 1b is \(\textbf{q}_\textrm{AFM} = \left( 0,0,\frac{3}{2}\right)\). The temperature dependence of the AFM Bragg peaks, plotted in Fig. 2a, shows that \(\hbox {FeTiO}_3\) undergoes an AFM phase transition at \(T_\textrm{N}\) = 59 K. The weak intensity of the Bragg peak at \(\textbf{Q} = \textbf{q}_\textrm{AFM}\) (\(\parallel c\)) demonstrates that its spins are oriented nearly perpendicular to the honeycomb planes. Its rocking scans in Fig. 2b reveal, however, that the spins are not absolutely parallel to the c axis; otherwise, its intensity would be zero. The angle of spins tilting off the c axis has previously been estimated to be as small as \(2-7^\circ\)46,47. Our neutron diffraction data are consistent with these previous reports. Since the tilting angle is sufficiently small, we will ignore it in our linear spin wave calculations. We also note that the single crystal used in this study consists almost entirely of a single domain satisfying the reflection condition \(-h+k+l=3n\).

Fig. 2
figure 2

(a) Temperature dependence of the integrated intensities for three different AFM Bragg peaks where \(\mathbf {q_\textrm{AFM}} = \left( 0, 0, \frac{3}{2}\right)\). (b) Rocking scans of the AFM Bragg peak at \(\textbf{Q} = \mathbf {q_\textrm{AFM}}\). (c) Powder-averaged INS intensity \(S(Q,\hbar \omega )\) at T = 6 K (\(E_i\) = 78.7 meV). (d) The comparison of the integrated \(S(0 < Q\le 2\) \(\hbox {\text{\AA }}^{-1},\hbar \omega )/[n(\hbar \omega )+1]\) between T = 6 and 150 K.

The powder-averaged INS intensities \(I(Q,\hbar \omega )\) at \(T=6\) K in the AFM phase are shown in Fig. 2c. Strong intensities are observed at low Q extending up to \(\hbar \omega \approx\) 16 meV, which are separated into two bands below and above \(\hbar \omega =10\) meV, respectively. The excitation intensities decreasing as a function of momentum transfer suggest that these intensities are due to magnetic scattering. In contrast, the intensities observed at 18 meV and higher, which enhance at higher Q, are ascribed to phonon scattering. The temperature-independent scattering function was obtained by integrating the INS intensities over the low-Q range (\(Q \le 2\) \({\text{\AA }}^{-1}\)) such that \(S(\hbar \omega )=(k_f/k_i)^{-1}[n(\hbar \omega )+1]\int \!I(Q,\hbar \omega )dQ\), where \(n(\hbar \omega ) = [e^{\hbar \omega /k_BT}-1]^{-1}\) is the Bose-Einstein distribution function and \(k_i\) and \(k_f\) are the initial and final neutron wave vectors, respectively. As shown in Fig. 2d, we find that the single-magnon excitations in the AFM phase are limited to below 16 meV.

Shown in the upper row of Fig. 3 are the magnon excitation spectra observed along selected high-symmetry directions on the \(\left[ h\, k\, \frac{3}{2}\right]\) plane crossing the AFM zone center at \(\textbf{q}_\textrm{AFM}\). These directions are marked in Fig. 1d as thick lines. At the 2D Brillouin zone center \(\textbf{Q}_\Gamma =(0, 0)\), a strong magnon mode appears around \(\hbar \omega\) = 4.9 meV and disperses upward in all directions in the honeycomb plane [see Fig. 3a]. This mode is ascribed to the two sublattice spins of the honeycomb ferromagnet fluctuating in-phase with respect to each other and thus called the acoustic magnon by analogy to acoustic phonons. Figure 3a shows that in the [h h] direction the acoustic magnon mode disperses up to \(\approx\)12 meV at the zone boundary \(\textbf{Q}_\textrm{M}=\left( \frac{1}{2},\frac{1}{2}\right)\). The other mode, ascribed to the out-of-phase spin fluctuations and thus called optical magnon, is visible along the \(\left[ h-\frac{1}{2}\; h+\frac{1}{2}\right]\) direction plotted in Fig. 3b. It is clearly shown that the optical magnon continues up to \(\hbar \omega =\) 15.2 meV at \(\textbf{Q}_\Gamma\), which is the maximum energy of the magnon band. Notice that the scattering intensities are not equal between left (\(k<0\)) and right (\(k>0\)) half panels in Figs. 3c and 2d. Such intensity asymmetry is ascribed to the crystal structure of \(\hbox {FeTiO}_3\) belonging to the rhombohedral \(R\bar{3}\) space group, and is also consistent with the fact that the sample used in this work consists of a single domain.

Fig. 3
figure 3

(ad) The upper row shows INS spectra in the selected \(\textbf{Q}\) directions marked in Fig. 1d while (eh) the lower row shows the linear spin wave calculations based on the DM model discussed in the text using the parameters listed in Table 1. The scattering intensities \(I(\textbf{Q},\hbar \omega )\) are integrated over the range \(1.3 \le L \le 1.7\) commonly for the experimental data and calculations. The integration ranges for the transverse in-plane wave vectors are written at the bottom of calculated plots in (eh), which also apply for the corresponding experimental plots. The letters on top of (ad) mark the high-symmetry positions on the 2D Brillouin zone. All calculations are convoluted by the instrumental resolution, the sample mosaic distribution and an additional Gaussian broadening of the 0.8-meV full width half maximum52. The vertical dashed lines indicate the Dirac points marked K.

Upon close inspection, we notice apparent discontinuities in magnon dispersion at the Dirac point \(\textbf{Q}_\textrm{K}=\left( \frac{1}{3}, \frac{1}{3}\right)\) where the two magnon modes are expected to cross each other. These discontinuities are consistently observed at \(h=\pm \frac{1}{3}\), \(h=\pm \frac{1}{6}\), and \(k=\pm \frac{1}{2}\) in Fig. 3a, b and d, respectively, suggesting a finite gap opening via avoided crossing. Since the magnon density of states will vanish at Dirac points, neutron scattering data may show gap-like features with intensity minima even when a true energy gap does not open.25,44 Therefore, to confirm the validity of the suspected gap opening, it is important to compare experimentally observed intensities against model calculations after being integrated over equal momentum widths. First, we plot in Fig. 4a and b the L-dependence of the scattering intensities at \(\textbf{Q}_\textrm{K}\) integrated over a narrow range of \(dh =\pm 0.02\) along [h h] and \(dK =\pm 0.03\) along [\(-k\) k]. The plot clearly shows that two intensity maxima are separated by \(\approx\) 1.2 meV with virtually no L dependence within the range up to \(L = 6\). It is notable that the apparent separation is significantly larger than the energy resolution of 0.65 meV in full-width-half-maximum estimated at the energy transfer of \(\hbar \omega\) = 10.5 meV. We then integrated the data over \(0 \le L \le 6\) to improve the statistics and then plotted a \(\textbf{Q}\)-E slice in Fig. 4c along the [\(-k\) k] direction orthogonal to \(\textbf{Q}_\textbf{K}\). It is apparent that the upper and lower magnon modes clearly deviate from linear dispersion relations near the Dirac point at \(k = 0\), which is more clearly visible in a series of constant-\(\textrm{Q}\) scans plotted in Fig. 4f. We therefore conclude that the Dirac magnons in \(\hbox {FeTiO}_3\) are gapped and massive.

Fig. 4
figure 4

(a,b) L dependence of the INS spectra at the Dirac wave vector displaying the gap opening. Open circles in (a) are the fitted Gaussian peak positions as shown in (b), and the dashed lines mark their averages. In (b), the lower horizontal axis is used to indicate the average L value of each data cut for Gaussian fittings. (c) The INS spectra along the \([-K K]\) direction passing though the Dirac wave vector at \(\textbf{q}_K\). The data are integrated over \(0 \le L \le 6\) to gain statistics. For comparison, model calculations with (d) \(A=\eta\) = 0 and (e) A = 0.04 meV & \(\eta\) = 0, respectively, were performed and integrated over the same momentum range as in (c). The dashed and solid lines in (ce) are the calculated magnon dispersions for \(A=\eta\) = 0 and A = 0.06 meV & \(\eta\) = 0, respectively. The plots in (f) show the constant-\(\textbf{Q}\) cuts taken from (c) revealing magnon peaks. (g) Open circles are the INS intensity at the Dirac wave vector \(\textbf{Q}_\textrm{K} = \left( \frac{1}{3},\frac{1}{3}\right)\). Red solid and black solid lines are the calculations using A = 0.04 meV and A = 0, respectively, with \(\eta =0\). On the other hand, the blue dash-dotted line is calculated using A = 0 and \(\eta\) = 0.71 meV.

Linear spin wave calculations of magnon excitations

To quantitatively evaluate the Dirac magnon gap and also deduce its origin, we setup the spin Hamiltonian including isotropic Heisenberg exchanges and additional SOC terms. We tested two models to account for the SOC on the honeycomb lattice: (1) the NNN antisymmetric DM exchanges (DM model)26,27 and (2) the NN bond-directional anisotropic exchange (BA model) suggested for \(\hbox {CoTiO}_3\)41. In the former DM model, the time-reversal symmetry between two sublattice excitations will be broken when the DM vector \(\textbf{A}\) is parallel to the c axis.18 Examples of such gap opening were experimentally observed for Dirac magnons in \(\hbox {CrI}_3\) and \(\hbox {CrGeTe}_3\), respectively.21,22,23,24 The relevant spin Hamiltonian is written as follows:

$$\begin{aligned} H_\textrm{DM} = -\frac{1}{2}\sum _{ij} [J^{xy}_{ij}(S^x_iS^x_j+S^y_iS^y_j)+J^z_{ij}S^z_iS^z_j] + \frac{1}{2}\sum ^{NNN}_{ij}\textbf{A}\cdot \textbf{S}_i\times \textbf{S}_j. \end{aligned}$$
(1)

Notice that the above Hamiltonian includes exchange anisotropies between in-plane and out-of-plane bonds, which we apply only to the NN exchanges (\(J_1^{xy} < J_1^z\)). Such easy-axis exchange anisotropy may alternatively be replaced by a single-ion anisotropy field, i.e. \(+D(\textbf{S}\cdot \hat{c})^2\), where \(2D=3(J^{xy}_1-J^z_1)<0\) while setting \(J^z_1 = J^{xy}_1\). On the other hand, the spin Hamiltonian for the BA model is written as follows:

$$\begin{aligned} H_\textrm{BA} = -\frac{1}{2}\sum _{ij} [J^{xy}_{ij}(S^\parallel _iS^\parallel _j+S^\perp _iS^\perp _j)+J^z_{ij}S^z_iS^z_j] - \frac{1}{2}\sum ^{NN}_{ij}\eta (S^\parallel _iS^\parallel _j-S^\perp _iS^\perp _j). \end{aligned}$$
(2)

In the above Hamiltonian, an additional bond-directional anisotropy is imposed on the NN in-plane exchange using the term \(\eta\) such that \(\eta \equiv J^\parallel _1-J^{xy}_1 =J^{xy}_1-J^\perp _1\), where \(J^\parallel\) and \(J^\perp\) act on spin components parallel (\(S^\parallel\)) and perpendicular (\(S^\perp\)), respectively, to the NN \(\hbox {Fe}^{3+}\)-\(\hbox {Fe}^{3+}\) bonds [see Fig. 1d]. We stress that this BA model is different from the Kitaev K-\(\Gamma\) model in which three orthogonal directions are common for all NN bonds. We also note that the BA model does not include cross terms between the parallel and perpendicular components.

In all of our calculations, we ignore the additional weak gapless and dispersive mode below \(\hbar \omega \le 4\) meV reported in earlier works46,48. Although its dispersion is reminiscent of an easy-plane anisotropy in \(\hbox {CoTiO}_3\)41, on the contrary \(\hbox {FeTiO}_3\) exhibits strong easy-axis anisotropy in dc magnetization46. Their origin is so far unknown although we tentatively conjecture that it is probably related to the small tilting of ordered spins. Nevertheless, its presence is unlikely to affect the physics of high-energy magnons near the Dirac point. We thus ignore it and adopt the easy-axis anisotropy (\(J_1^{xy} < J_1^z\)).

The best-fit parameters obtained by the least square fitting of observed magnon energies are listed in Table 1, among which those for the DM model were used to simulate the magnon spectra plotted in Fig. 3e–h. Since the experimentally observed spectra were noticeably broader than the instrumental effect, we convoluted the calculated intensities with Gaussian functions of 0.8 meV in full-width-half-maximum in addition to the energy-dependent instrumental resolution. As shown in Fig. 3, the calculated spectra agree excellently with the experimental observations in all directions. In Fig. 3, we opt to show the simulated results only of the DM model because the corresponding results of the BA model are found to be almost indistinguishable. Most importantly, plots in Fig. 3 clearly show finite magnon gap opening at all available Dirac points. To open the energy gap of \(\Delta _\textrm{K}\approx 1.2\) meV, we find that a small DM exchange of \(|\textbf{A}|\) = 0.04 meV is sufficient in the DM model. In the case of BA model, on the other hand, a large bond-directional anisotropy of \(\eta = 0.71\) meV was necessary to open the gap of equivalent magnitude. Finally in Fig. 4c–e, we show the effect of the DM exchange on the gap opening more clearly using the data and simulations integrated over a wide range of \(0 \le L \le 6\). Without the DM exchange or bond-directional anisotropy, \(|\textbf{A}| = \eta = 0\), two magnon modes directly cross each other exhibiting no intensity minimum at the Dirac point [see Fig. 4d]. In contrast, either the finite DM exchange or bond-directional anisotropy introduces a clear two-peak structure with an intensity minimum in the middle.

Table 1 The best-fit parameters for the two model spin Hamiltonians obtained by least square fittings to the INS data.

Discussion

As discussed above, simulations based on the two models, DM and BA, equally well reproduce the experimentally observed magnon excitations. Although the intensity around the magnon gap seems better reproduced by the BA model [see Fig. 4g], the difference between the two simulations is not significant. We tentatively suppose that the apparent asymmetry of the peak intensities is likely to be an experimental artifact. In the case of DM model, the relative strength of the DM exchange, \(|\textbf{A}|/J^{xy}_1 = 0.05\), is comparable to the ratios estimated in van der Waals \(\hbox {CrI}_3\) (0.09/2.11 = 0.04) or \(\hbox {CrGeTe}_3\) (0.20/2.76 = 0.07)22,24. Although the ligand ions are significantly lighter in mass. Since the Fe-O honeycomb planes are significantly more buckled with highly distorted \(\hbox {FeO}_6\) octahedra, it is plausible that \(\hbox {FeTiO}_3\) may have DM exchanges of a strength comparable to those of \(\hbox {CrI}_3\) or \(\hbox {CrGeTe}_3\). In the case of the BA model, we find that the exchange anisotropy, \(J^\parallel _1/J^\perp _1 = 12.7\), is larger than the ratio estimated in \(\hbox {CoTiO}_3\) by an order of magnitude41. Such large anisotropy suggests that the BA model may not be physically plausible49.

To test the plausibility of two considered models, we calculated the strength of exchange interactions using the OpenMX package based on the fully relativistic density functional theory50. Fixing the magnetic moments along the direction parallel to the c axis, we obtained isotropic Heisenberg exchanges \(J_1 = 1.00\) meV, \(J_{c1} = 0.02\) meV, \(J_2 = -0.14\) meV, \(J_{c2} = -0.17\) meV, \(J_3 = -0.04\) meV, and \(J_{c3} = -0.23\) meV, which are broadly consistent with the best-fit exchange parameters. Most importantly, we obtain the out-of-plane DM exchange \(|\textbf{A}| = 0.03\) meV but negligible bond-directional anisotropic exchange anisotropy (\(\eta \approx 0\)). Therefore, we finally conclude that the DM exchange is responsible for the magnon gap opening at the Dirac points.

Finally, we checked the topological characteristics of the massive Dirac magnons in \(\hbox {FeTiO}_3\) by calculating Berry curvatures and Chern numbers. Since the honeycomb ferromagnet consists of two spin sublattices, its magnon Berry curvature can be calculated using the coefficients of the Pauli matrices when its spin Hamiltonian is diagonalized on the basis of spinor operators19. Following the method described in Ref. 18, we first calculate the \(\textbf{h}(\textbf{Q})\) vector as follows:

$$\begin{aligned} \textbf{h}(\textbf{Q}) = \begin{pmatrix} h_x\\ h_y\\ h_z \end{pmatrix} = \Sigma _j^3 \begin{pmatrix} -J_1 S \cos {(\textbf{Q}\cdot \textbf{r}_{1,j})}-J_3 S \cos {(\textbf{Q}\cdot \textbf{r}_{3,j})}\\ J_1 S \sin {(\textbf{Q}\cdot \textbf{r}_{1,j})}+J_3 S \sin {(\textbf{Q}\cdot \textbf{r}_{3,j})}\\ 2|\textbf{A}|S \sin {(\textbf{Q}\cdot \textbf{r}_{2,j})}) \end{pmatrix} \end{aligned}$$
(3)

.

Above, \(\textbf{r}_{1,j}\), \(\textbf{r}_{2,j}\) and \(\textbf{r}_{3,j}\) are the spatial vectors corresponding to the in-plane exchanges \(J_1\), \(J_2\)/\(\textbf{A}\), and \(J_3\), respectively. For simplicity, we treated the system as a two-dimensional lattice and ignored the exchanges between honeycomb planes. Figure 5a shows the calculated distribution of \(\hat{\textbf{h}} = \textbf{h}/|\textbf{h}|\) using the best-fit parameters of the DM model. We note that \(h_z\) is focused around six Dirac wave vectors alternating between \(h_z>0\) and \(h_z<0\) around \(\textbf{Q}_\textrm{K}\) and \(\textbf{Q}_{\textrm{K}^\prime }\), respectively. We subsequently calculated the Berry curvature \(\Omega _\pm = \mp \hat{\textbf{h}}\cdot (\partial _{k_x}\hat{\textbf{h}}\times \partial _{k_y}\hat{\textbf{h}})\) for the magnons mode above and below the Dirac magnon gap, respectively. Figure 5b shows that the Berry curvature is positive (negative) for the upper (lower) mode and mostly focused around the six-fold Dirac wave vectors. The Chern number is calculated to be \(C^\pm = (2\pi )^{-1} \int _{BZ}\Omega _\pm d^2 k = \pm 1\) for the upper and lower modes, respectively, by integrating the Berry curvature within the 1st Brillouin zone boundary. If we set \(\textbf{A} = 0\) instead, then \(h_z\) becomes zero while \(h_x\) and \(h_y\) remain unchanged. The Berry curvature subsequently vanishes throughout the Brillouin zone leaving the Chern number to be zero. We thus conclude that the DM exchanges indeed provide an effective massive to Dirac magnons and thereby establishes topological spin excitations in \(\hbox {FeTiO}_3\).

Fig. 5
figure 5

Topological characteristics of \(\hbox {FeTiO}_3\) calculated for the DM model including \(J_1\), \(J_2\), \(J_3\), and \(|\textbf{A}|\) listed in Table 1. (a) The arrows show the in-plane directions of \(\hat{\textbf{h}}(\textbf{Q})\) or \((h_x, h_y)/|\textbf{h}|\) in Eq. (3). The colored contours superposed on them show the out-of-plane component \(h_z/|\textbf{h}|\). (b) Topological Berry curvature of the magnon mode above the Dirac gap calculated using \(\hat{\textbf{h}}(\textbf{Q})\) in (a). By integrating the Berry curvature in (b) over the 1st Brillouin zone, we obtain the Chern number \(C^+ = 1\).

In summary, we investigated the magnon excitations of the three-dimensional honeycomb magnet \(\hbox {FeTiO}_3\) with spins oriented nearly parallel to the c axis. Our INS experiment reveals a finite gap opening of \(\Delta _\textrm{K} = 1.2\) meV at the in-plane Dirac crossings, which appears independently of the out-of-plane wave vectors. Using the linear spin-wave and DFT calculations, we find that the antisymmmetric DM exchange between NNN Fe\(^{2+}\) spins is responsible for the Dirac gap opening. Based on these observations and analysis, we conclude that the Dirac magnons in \(\hbox {FeTiO}_3\) are massive and topological. Our work demonstrates that \(\hbox {Fe}^{2+}\) ions of the honeycomb oxide can have a sufficiently strong spin-orbit coupling to establish topological spin excitations on honeycomb lattices.

Methods

A large single crystal of \(\hbox {FeTiO}_3\) was grown using the floating zone method in an image lamp furnace. To control the oxidation state, we used the mixed gas environment of 37.0% \(\hbox {CO}_2\), 3.7% \(\hbox {H}_2\), and 59.3% \(\hbox {N}_2\) in volumes51. For INS, we used a large piece cut out of the rod in which the mosaic spread of and aligned its HHL direction on the horizontal scattering plane. Although the crystal rod included a few growth domains, their relative tilt was less than \(\hbox {1}^\circ\) in angle within honeycomb planes and approximately \(\hbox {1.5}^\circ\) along the perpendicular directions.

To observe magnon excitations, the BL23 POLANO time-of-flight spectrometer at the Materials and Life Science Experimental Facility of the Japan Proton Accelerator Research Complex was used in the multiple incident energy mode with \(E_i =\) 11.57, 30.2 and 78.7 meV. In addition, the Cold-TAS triple-axis spectrometer at the HANARO Research Reactor Facility of the Korea Atomic Energy Research Institute was used to measure the temperature dependence of the antiferromagnetic Bragg peaks.

Momentum-dependent spin wave energies and intensities were calculated using the SpinW software52. The fully relativistic density functional theory calculations were performed using ‘OpenMX’ software package50. Generalized gradient approximation53 was adopted for the exchange-correlation potential. The magnetic moments were aligned along the c axis to obtain the out-of-plane DM components. To describe the Mott insulating phase, the electron interaction \(U_\textrm{eff}\) = 2.0 eV was used within DFT+U scheme54. The noncollinear exchange interactions \(J_{ij}\) were obtained based on magnetic force theorem as implemented in our ‘Jx’ code55,56,57. The AFM ground state was found to be stabilized with the local magnetic moment of Fe site obtained to be 4.34\(\mu _B\) (3.74\(\mu _B\) for spins, 0.60\(\mu _B\) for orbitals). The calculated exchange interactions consistently indicated the layered AFM spin order.