Introduction

Magnons and phonons, quanta of spin waves and lattice vibrations respectively, constitute two fundamental collective excitations in ordered magnets. When there is a strong magnon-phonon coupling, they can be hybridized to form a gap at the intersection of the original magnon and phonon bands (Fig. 1a)1,2,3. The hybridized bands feature mixed, interconverted and reversed magnonic and phononic characters, and the associated quasiparticles are defined as magnon polarons1,2,3. Magnon polarons can resonantly enhance the spin-pumping effect4, and provide a phonon-involved way to generate and manipulate spin currents carried by magnons thanks to their hybrid nature5,6, signifying promising potentials in spintronics technology7,8. More recently, it has been predicted that magnon-polaron bands can exhibit nonzero Chern numbers and large Berry curvatures, giving rise to the thermal Hall effect9,10,11,12,13,14,15.

Fig. 1: Schematic of magnon polarons and two-dimensional magnetism in a multiferroic Fe2Mo3O8.
figure 1

a Schematic of the magnon-polaron bands after the hybridization between the original magnon and phonon bands illustrated by the thin dashed lines. Colors indicate the weights of the magnonic and phononic components. b Crystal and magnetic structures of Fe2Mo3O8. The material belongs to the polar P63mc space group (No. 186) and exhibits a collinear antiferromagnetic configuration. Fet and Feo label the ions in the centers of the tetrahedra and octahedra, respectively. Arrows indicate the magnetic moments on Fe atoms. c Top view of one Fe-O layer. Black arrows mark the in-plane DM vectors for the NN sites. d Temperature dependence of the integrated intensities of the magnetic Bragg peak (1, 0, 0). Error bars represent one standard deviation throughout the paper. The solid curve is a fit to the intensity with \(I\propto {(1-T/{T}_{{{{{{\rm{N}}}}}}})}^{2\beta }\), where TN = 59.33 ± 0.07 K and β = 0.129 ± 0.005 are the Néel temperature and critical exponent, respectively. The inset shows elastic scans across (1, 0, L) along [001] direction at 1.5 and 71 K. e Spin excitation spectra around (1, 0, 0) along [100] direction at 1.5 K. The cut-off bright spot near 17 meV is a false signal caused by the small scattering angle. f Energy scans at (1, 0, 0) measured at various temperatures. Solid curves are fits with Gaussian functions. Dashed lines in d–f are guides to the eye.

Full size image

These proposals have motivated sustained experimental efforts in chasing for magnon polarons4,5,6,16,17,18,19,20,21,22,23,24,25,26. However, direct spectroscopic evidence with their delicate band structures being explicitly unveiled by neutron spectroscopy is still rare. This is primarily because: i) materials with strong magnon-phonon coupling that can result in such excitations with prominent features are scarce; ii) magnons and phonons are rarely observed in the same energy-momentum window by neutron spectroscopy due to their different dynamical structure factors, which hinders the exploration of the interaction effects between them. To overcome these difficulties, Fe2Mo3O8 is a prime candidate material10,27,28,29,30,31,32,33.

Fe2Mo3O8 is a multiferroic material with a polar P63mc space group (No. 186). Below the Néel temperature TN ~ 60 K, both the space-inversion and time-reversal symmetries are broken27,28,29,34,35,36, as shown in Fig. 1b. The magnetism in Fe2Mo3O8 arises from Fe2+ ions, while Mo4+ ions form nonmagnetic spin-singlet trimers28,29. The Fe2+ ions on each Fe-O layer form a bipartite honeycomb network with different magnetic moments in corner-shared tetrahedra and octahedra (Fig. 1b, c)34,35,36. The absence of an inversion center between nearest-neighbor (NN) Fe sites allows for a non-zero in-plane Dzyaloshinskii-Moriya (DM) interaction. Fe2Mo3O8 exhibits a long-range collinear antiferromagnetic order below TN, with antiparallel yet uncompensated moments on each Fe-O layer stacking antiferromagnetically along the c axis (Fig. 1b)34,35,36. Furthermore, it is noteworthy that the magnetic configuration of Fe2Mo3O8 can be controlled by either an external magnetic field or chemical doping, leading to a metamagnetic transition into a ferrimagnetic state (Supplementary Fig. 1c, d)28,29. More crucially, there has been accumulating evidence that the spin and lattice degrees of freedom are strongly coupled in Fe2Mo3O810,27,28,29,30,31,32,33, rendering it a promising platform to probe the long-sought magnon polarons4,5,6,16,17,18,19,20,21,22,23,24,25,26.

In this work, we perform high-resolution neutron spectroscopy measurements and fully map out the magnon and phonon bands in a multiferroic Fe2Mo3O8 with strong magnon-phonon coupling10,27,28,29,30,31,32,33. Due to the acquisition of spin components from magnon conversion, the acoustic phonons show up together with magnons at small momenta. By examining the interaction effects between them, we directly observe the long-sought magnon polarons, which we show to be topologically nontrivial. These results not only unambiguously identify a new type of excitations, but also provide a fresh ground to study topological states.

Results

Magnons at high energies

Figure 1d shows the elastic neutron scattering results for single crystals of Fe2Mo3O8. At 1.5 K, a significant increase in the magnetic scattering intensity is observed at the magnetic Bragg peak (1, 0, 0), indicating the establishment of an antiferromagnetic order. On the other hand, (1, 0, 1) corresponds to the magnetic Bragg peak for the ferrimagnetic order34,36. The fitting of the temperature dependence of the integrated intensities at (1, 0, 0) yields TN = 59.3 K and a critical exponent of 0.129, which closely matches the expected value of 0.125 for a two-dimensional Ising system37. These results are consistent with the magnetic susceptibility measurements (Supplementary Fig. 1c), suggesting that Fe2Mo3O8 is a two-dimensional collinear antiferromagnet with strong magnetocrystalline anisotropy along c axis. Figure 1e shows the inelastic neutron scattering (INS) results along [100] direction at 1.5 K, obtained on a triple-axis spectrometer. No acoustic bands are found to disperse up from the magnetic Bragg peak (1, 0, 0). Instead, only two seemingly flat bands with a sizable gap are observed between 10 and 16 meV. Upon warming, these two modes soften slightly, while the scattering intensities decrease dramatically (Fig. 1f). Above TN, the peaks disappear. This suggests that the two modes correspond to the spin-wave excitations originating from long-range magnetic order, rather than the crystal-electric-field excitations38. Given that there are four sublattices of Fe2+ ions in the magnetic unit cell of Fe2Mo3O8 (Fig. 1b), we consider these two magnon modes to be doubly degenerate in the collinear antiferromagnetic state, as also found in its sister compound Co2Mo3O839. We note that the energy scales of the excitation spectra in the other isostructural compound Ni2Mo3O838 are rather different from those in Fe2Mo3O8 (Fig. 1e) and Co2Mo3O839. In Ni2Mo3O8, the ground state of the crystal-electric-field levels of Ni2+ ions is a nonmagnetic singlet, while the interplay between crystal-electric-field effect and magnetic exchange interactions gives rise to a magnetic order with a relatively low transition temperature of TN = 5.5 K. As a consequence, spin waves are observed below 1.5 meV, and higher-energy excitations are believed to arise from the crystal-electric-field spin excitons with a robust temperature dependence38.

Anomalous phonons at low energies

To examine the excitation spectra of Fe2Mo3O8 in a larger momentum-energy space, we next performed INS measurements on a time-of-flight spectrometer. The two magnon modes between 10 and 16 meV can also be observed at 6 K along both [100] and [110] directions (Fig. 2a, b). These two modes are flat along [001] direction (Fig. 2c), indicating negligible interlayer coupling consistent with the two-dimensional nature of the magnetism deduced from the critical exponent (Fig. 1d). This behavior allows us to integrate the intensities along [001] to improve the statistics when studying the magnon-related excitations, as we did in Fig. 2a, b. Turning to the lower energy window, we observe additional excitations alongside the two intense magnon modes (Fig. 2a–c). These low-energy modes introduce multiple bands that cannot be solely explained by linear spin-wave theory. Even when accounting for the nonlinearity of spin waves, the extended spin-wave theory fails to account for the presence of multiple bands at significantly lower energies40. Therefore, we propose that these additional modes to have a phononic origin, which can also be supported by the following facts. The scattering intensities of these low-energy excitations become stronger as the wavevector Q increases (Supplementary Fig. 2a–c), which is characteristic of phonons. Figure 3a shows a specific mode among these excitations, with an onset energy around 5 meV along the [100] direction. Notably, we find such onset excitations at (1, 0, 1) also correspond to the maximum of the excitations along the [001] direction as shown in Fig. 3b, while no significant signal can be observed at the intense magnetic Bragg peak (1, 0, 0) (Fig. 1d). Furthermore, the scattering intensities of such excitations tend to be stronger at larger Qs as shown in Fig. 3c. These findings indicate that such a specific mode can be interpreted as phonons connected by a saddle point around 5 meV, rather than gapped spin waves. Our theoretical calculations, which reproduce the multiple bands well, provide further support for the interpretation of these low-energy modes as phonons (Supplementary Fig. 4).

Fig. 2: Magnon-polaron excitation spectra.
figure 2

a–c INS results of the excitation spectra measured at 6 K along [100], [110] and [001] directions, respectively. Where applicable, the integration range used to plot the spectra throughout the paper can be found in Supplementary Table 1. In c, raw data on the negative L side have been symmetrised to the positive side to improve the statistics. Due to the limited data coverage, this symmetrisation procedure is applicable only for L ≤ 2 r.l.u. For better display, these data in a–c are plotted in a logarithmic scale of intensity. Bright spots around 2 meV and elongated signal near the elastic line in a–c are spurious peaks caused by saturation of the neutron detectors. d1, e1 Zoom-in view of magnon polarons in a and b, respectively. d2, e2 Spectra obtained by the two-dimensional-curvature method on d1 and e1, respectively. f Three-dimensional Brillouin zones with high-symmetry points and paths. Ticks in d1, e1 mark the energies corresponding to the constant-energy contours plotted in Fig. 4a–f. Dashed lines in d2, e2 are guides to the eye to mark the top (t) and bottom (b) magnon-polaron bands.

Full size image
Fig. 3: Low-energy phonon excitation spectra.
figure 3

INS results of the excitation spectra measured at 6 K along [100] (a) and [001] (b, c) directions, respectively. The spectra in b, c correspond to the out-of-plane variations with different wavevectors Hs. Solid and dashed squares in a–c represent the saddle point of phonon spectra around 5 meV at H = 1 and 2 r.l.u., respectively. d–e Same as in a–c, but measured at 100 K. In b, c, e and f, data on the negative L side have been symmetrised to the positive side to improve the statistics. The raw data obtained with lower Ei = 12 meV have been corrected by Bose factor and are plotted in a logarithmic scale of intensity.

Full size image

On the other hand, these modes exhibit some anomalous behaviors that suggest they are not ordinary phonons. At 100 K, which is above TN, they become invisible at small Qs but persist at large Qs, as shown in Fig. 3d–f and Supplementary Fig. 2d–f. Furthermore, it is noteworthy that the saddle point of the phonon spectra around 5 meV, previously considered as an electromagnon (~1.2 THz) in multiferroics by some optical measurements30,31, will be electric-dipole active in the antiferromagnetic phase30,31. These results indicate these low-energy phonons acquire some spin components17 through the strong magnon-phonon coupling, as discussed later, leading to additional magnetic scattering intensities at small Qs at 6 K. However, at 100 K, as the magnons collapse, phonons recover their original properties and can only be observed at large Qs, where the intrinsic dynamic structure factors are sufficiently large.

Formation of magnon polarons

The appearance of acoustic phonons at small momenta together with magnons enables us to examine the interaction effects between them, and we now focus on the regions where the weak and dispersive acoustic phonons tend to intersect with the intense and relatively flat magnons (Fig. 2a, b). Apparently, there exist spectral discontinuities at the nominal intersections between the phonons and magnons (Fig. 2d1, e1), indicative of the gap opening. To better resolve the gaps, we use a two-dimensional curvature method41 to increase the sharpness of the bands (Fig. 2d2, e2). There are multiple gaps around the two magnon bands at about 11 and 14 meV. Together with the gap opening, the sudden change of the dispersions and the intensities in proximity to the gap unambiguously indicate that magnons and phonons are strongly hybridized and inverted, and consequently two magnon-polaron bands separated by the gap form, as illustrated in Fig. 1a. For the top magnon-polaron band, it changes from dominant magnonic to phononic character as it propagates. As a result, the band’s velocity increases but the intensity decreases dramatically around the original intersection. The trend is opposite for the bottom magnon-polaron band. The observation that the low-energy dispersive phonons convert into magnons with enhanced intensity supports our earlier explanations for the anomalous phonon behaviors, as phonons involved with magnon conversion can carry spins17. Importantly, when the magnons collapse above TN, phonons revert to their original dispersions and scattering intensities (Fig. 3 and Supplementary Fig. 2). The spectroscopic characteristics of hybrid magnon polarons and low-energy phonons acquiring spin components in the resonant and off-resonant regions, respectively, are both manifestations of the strong magnon-phonon coupling in Fe2Mo3O8.

To further illustrate the formation and evolution of the magnon-polaron excitations, we plot a series of constant-energy (E) contours near the two anticrossing regions as shown in Fig. 4a–c (higher energies), and d–f (lower energies). In each region, there are always two sets of excitations around the zone center, representing two magnon-polaron bands separated by the gap. The shapes and intensities of these two excitations change as the energy increases, manifesting the magnon-phonon interconversion within each magnon-polaron band during propagation. To better clarify these, we first plot a set of constant-E scans along [100] direction through the zone center (Fig. 4i), which correspond to the lower anticrossing region (Fig. 4d–f). As the energy increases, it can be found that there is a spectral weight transfer between the two magnon-polaron bands as they propagate. A similar approach can be applied to the constant-E scans at other energies, so that we can extract a series of fitted peak centers and areas. We plot the relative peak position (Δq) and the ratio of the spectral weight for each band to the total spectral weight of the two magnon-polaron bands as a function of energy (Fig. 4j). When the top (inner) and bottom (outer) magnon-polaron bands are approaching the anticrossing point from low energies, they are of primarily magnonic and phononic natures respectively, so the top band dominates the spectral weight as magnons are much stronger. When the two bands reach the anticrossing point, where Δq has its minimum, strongest hybridization happens so that the magnonic and phononic components become comparable. Across this position, the main components of the two bands are reversed, and so are their relative intensity ratios. Eventually, the bottom band which is of primarily magnonic nature dominates the spectral weight. Figure 4g, h shows similar behaviors for the constant-E scans corresponding to the higher anticrossing region (Fig. 4a–c). We also examine the magnon-polaron excitations along the orthogonal [−120] direction, and the results are also similar (Supplementary Fig. 3). These results elaborate the band inversion between the original magnon and phonon bands. Together with the gaps in the dispersions (Fig. 2), these constitute the hallmarks for the magnon polarons.

Fig. 4: Evolution of the magnon-polaron bands.
figure 4

Constant-energy contours at a series of selected energies near the higher anticrossing point around 14.3 meV (a–c) and lower anticrossing point around 11.3 meV (d–f) in the (H, K, 0) plane. The energies are marked in Fig. 2d1, e1. Arrows mark the two magnon-polaron bands. The inner and outer bands correspond to the top and bottom bands respectively. Solid lines indicate the Brillouin zone boundary. g, i Panels labeled 1-3 correspond to the profiles obtained from the cuts in the panels of a–c or d–f. Solid curves are fits with Gaussian functions. The yellow and green fillings illustrate the fitted peaks for the two magnon-polaron bands. The fitted areas are labeled as Sm→p (yellow) or Sp→m (green), based on the direction of the conversion between the magnonic (m) and phononic (p) natures. The asymmetry in the peak shapes on both sides of the origin in g, i is caused by the resolution effect on 4SEASONS. h, j Ratio of the spectral weight (left axis) and the peak interval (Δq, right axis) as a function of energy near the higher and lower anticrossing regions, respectively. The ratio is defined as Sm(p)→p(m)/(Sm→p + Sp→m) counting spectral weights on both sides of (1, 0, 0). Δq is the distance between the peak centers, averaged with the values on both sides. Vertical bars indicate the minimum of Δq. Solid curves are guides to the eye.

Full size image

Dzyaloshinskii–Moriya mechanism and topology

To understand the underlying mechanism of the observed magnon polarons, we develop an effective two-dimensional model that contains both the magnon and phonon terms (Methods). Without the magnon-phonon coupling, the obtained magnon and phonon dispersions are shown in Fig. 5a. Such dispersions capture the main features of the INS spectra (Fig. 2) except for the anticrossings between magnons and phonons. In principle, since magnetic exchange interactions depend on the relative ion positions, magnon-phonon coupling can arise due to the lattice vibrations. In Fe2Mo3O8, the magnon-phonon coupling induced by the DM interaction dominates those by other Heisenberg interactions (See details in the Methods). Generally speaking, only the component of the DM vector parallel to the magnetic moments will enter into the spin waves42. On the other hand, the DM vector here with only the in-plane component allowed, which is perpendicular to the magnetic moments (Fig. 1b, c), typically will not contribute to the magnetic ground state as well as the spin excitation spectra9,10,11. We consider it is the lattice vibrations that make the in-plane DM interaction come into play, which in return closely couples magnons and phonons (Methods). The mechanism that the vibrant DM vectors can disturb the procession of the magnetic moments is illustrated in Fig. 5c. In Fig. 5b, the calculated spectra of the coupled system are shown (Parameters can be found in Supplementary Table 2). The gaps between the magnon and phonon bands and interconversions of their spectral components (Fig. 5b) well reproduce the characteristics of the magnon polarons observed experimentally, indicating that the DM-interaction-induced magnon-phonon coupling can give rise to the magnon-polaron excitations. Importantly, near the gaps, it is found that large Berry curvatures can be induced (Fig. 5d, e). Accordingly, the Chern number for each band is calculated and labeled in Fig. 5b. These results show that the magnon-polaron excitations are topologically nontrivial, in accordance with our experimental observation of the band inversion between magnons and phonons. To account for the three-dimensional nature of the phonons, we have also developed a three-dimensional model, as presented in Supplementary Fig. 4.

Fig. 5: Calculations of the magnon-polaron bands and band topology.
figure 5

a Uncoupled magnons (orange) and phonons (blue) along the high-symmetry directions without the DM interaction. b Topological magnon polarons resulting from the DM-interaction-induced magnon-phonon coupling. The color illustrates the weight of the magnonic and phononic components. Chern numbers from the lowest to highest bands (labeled from 1 to 4) are C1 + C2 = −2, C3 = 4 and C4 = −2. Here, due to the degeneracy of the acoustic phonons at the Γ point, only the total Chern number is well defined for the lowest two bands. The dotted circles labeled α and β correspond to the lower and higher anticrossing regions respectively. c Schematic diagram for the DM-interaction-induced magnon-phonon coupling in Fe2Mo3O8. Due to the in-plane lattice vibrations, the otherwise canceled DM interaction becomes nonzero and comes into play, which induces the magnon-phonon coupling, labeled as mp. d, e Calculated Berry curvatures for bands 4 and 3, respectively.

Full size image

Discussion

By now, combing our experimental spectra and theoretical calculations, we have provided compelling evidence that in Fe2Mo3O8 there exist topological magnon polarons. Remarkably, in addition to the gap opening between the original magnon and phonon bands, we have also observed the band inversion. Such a case is analogous to that in topological insulators induced by the spin-orbit coupling, but involves the interconversion between magnons and phonons induced by the DM interaction. Therefore, our work provides new perspectives in seeking for topological states—that is to go beyond a single type of elementary excitations such as electrons, magnons or phonons only, and to consider the topology in their hybrid form.

We note that a recent magneto-Raman scattering study43 has reported the existence of topological magnon polarons due to the zigzag antiferromagnetic order in the monolayer FePSe3, where the magnon-phonon coupling originates from the anisotropic exchange interactions43. It implies that the topological nature of band-inverted magnon polarons can be inherent and resilient, irrespective of the particular form of magnon-phonon coupling9,10,11,12,13,14,15. In our model calculations, we find that the DM term is larger than the Heisenberg terms (Supplementary Table 2), supporting the strong magnon-phonon coupling induced by the DM interaction in Fe2Mo3O8. Notably, this coupling mechanism involves the interaction between in-plane phonons and in-plane magnons that emerges as the leading order in Fe2Mo3O8, surpassing the general magnetoelastic coupling1,13,14 that arises from perpendicular easy-axis anisotropy (See details in the Methods).

The strong DM-interaction-induced magnon-phonon coupling in Fe2Mo3O8 leads to the formation of topological magnon polarons in the magnon-phonon resonant region, with prominent features that are observed by our neutron spectroscopy measurements. Additionally, it gives rise to other intriguing phenomena even away from the anticrossing region, such as the anomalous phonons carrying spins (Fig. 3 and Supplementary Fig. 2), low-energy excitations with electric dipole activity30,31, and the emergence of thermal Hall effect33. These results showing strong magnon-phonon coupling and hybrid excitations in Fe2Mo3O8 suggest it to be a prime candidate in developing phonon-controllable spintronic devices.

Methods

Single-crystal growth and characterizations

High-quality single crystals of Fe2Mo3O8 were grown by the chemical-vapor-transport method44,45. The well mixed raw powder materials of Fe2O3, MoO2 and Fe in a stoichiometric molar ratio of 2:9:2 were sealed in an evacuated quartz tube with TeCl4 as the transport agent. The tube was placed into a two-zone horizontal tube furnace with the hot and cold end set at 980 °C and 845 °C, respectively. These temperatures were maintained for 10 days for the crystal growth, after which they were cooled naturally in the furnace to room temperature. The magnetization was measured with a 15.1-mg single crystal using the vibrating sample magnetometer option integrated in a Physical Property Measurement System (PPMS-9T) from Quantum Design.

Neutron scattering experiments

Our neutron scattering measurements were performed on 4SEASONS, a time-of-flight spectrometer located at the MLF of J-PARC in Japan46 and EIGER, a thermal-neutron triple-axis spectrometer located at the SINQ of PSI in Switzerland47. The sample array consisted of ~150 pieces of single crystals weighing about 3.39 g in total. They were coaligned using a backscattering Laue X-ray diffractometer and glued on aluminum plates using a trademarked fluoropolymer CYTOP-M. These plates were assembled together on an aluminum holder, which was then mounted into a closed-cycle refrigerator for measurements, in a manner that the (H, 0, L) plane was the horizontal plane, as shown in Supplementary Fig. 1a. For the measurements on 4SEASONS, we chose a primary incident energy Ei = 30.04 meV and a Fermi chopper frequency of 250 Hz, with an energy resolution of 1.56 meV at the elastic line. Since 4SEASONS was operated in a multiple-Ei mode48, it had other Eis of 11.93 and 17.95 meV, with respective energy resolutions of 0.56 and 0.84 meV at the elastic line. Note that on direct geometry time-of-flight spectrometers such as 4SEASONS, the energy resolution is improved as the energy transfer increases. We set the angle of the incident neutron beam direction parallel to c axis to be zero. Scattering data were collected by rotating the sample about [−120] direction from 60° to 180° in a step of 1° or 2° for 6 and 100 K, respectively. We counted 20 min for each step. In this setup, we found the data with Ei ~ 18 meV were overlapped with and contaminated by the elastic line from the next lower Ei in the region around 15.5 meV. To eliminate this, we used a similar setup as previous measurements, but suppressed the incident neutrons with Ei ~ 12 meV and lower by the disk chopper in the second measurements on 4SEASONS. In the work, the results with Ei ~ 12 meV and Ei ~ 30 meV were all based on the first measurement, while those with Ei ~ 18 meV at 100 and 6 K were based on the first and second measurements, respectively. These data were reduced and analyzed using the software suites Utsusemi49 and Horace50. For the measurements on EIGER, data were collected in the (H, 0, L) scattering plane with a horizontal-focusing analyzer. We fixed the final wavevector kf = 2.662 Å−1 corresponding to an energy of 14.7 meV. We used a hexagonal structure with the refined lattice parameters a = b = 5.773(3) Å and c = 10.054(3) Å28. The wavevector Q was expressed as (H, K, L) in the reciprocal lattice unit (r.l.u.) of \(({a}^{*},{b}^{*},{c}^{*})=(4\pi /\sqrt{3}a,4\pi /\sqrt{3}b,2\pi /c)\). In this paper, the measured neutron scattering intensities S(Q, E) from 4SEASONS were corrected by the magnetic form factor of Fe2+ ions and divided by the Bose factor via \({\chi }^{{\prime\prime} }({{{{{\bf{Q}}}}}},E)=| f({{{{{\bf{Q}}}}}}){| }^{-2}(1-{{{{{{\rm{e}}}}}}}^{-E/{k}_{{{{{{\rm{B}}}}}}}T})S({{{{{\bf{Q}}}}}},E)\), where kB was the Boltzmann constant.

Theoretical calculations of the bands and their topology

We start the calculations with an effective two-dimensional model that contains both the magnon (Hm) and phonon (Hp) terms. The symmetry allowed Hm on a bipartite honeycomb layer with two inequivalent Fe sites in Fig. 1c of the main text can be written as:

$${H}_{{{{{{\rm{m}}}}}}}=\mathop{\sum}\limits_{i < j}{J}_{ij}{{{{{{\bf{S}}}}}}}_{i}\cdot {{{{{{\bf{S}}}}}}}_{j}-\mathop{\sum}\limits_{i}{{{\Delta }}}_{i}{\left({S}_{i}^{z}\right)}^{2}+\mathop{\sum}\limits_{\langle ij\rangle }{{{{{{\bf{D}}}}}}}_{ij}\cdot ({{{{{{\bf{S}}}}}}}_{i}\times {{{{{{\bf{S}}}}}}}_{j}).$$
(1)

Here, Jij is the Heisenberg exchange interaction between the spins Si and Sj, which is considered up to the third-nearest neighbor (TNN) to fit the magnon spectra shown in Fig. 2 of the main text. It is noted that the nearest-neighbor (NN) exchange interaction J1 and TNN exchange interaction J3 are homogeneous over the whole lattice while the next-nearest-neighbor (NNN) exchange interactions can take different values for the bonds between the Fet sites (denoted by \({J}_{2}^{{{{{{\rm{t}}}}}}}\)) and those between the Feo sites (denoted by \({J}_{2}^{{{{{{\rm{o}}}}}}}\)). Δi is the single-ion anisotropy constant of the local spin at the site i, which can be distinct for Fet (denoted by Δt) and Feo (denoted by Δo) sites as well. Dij is the vector of the DM interaction between the NN sites i and j. This term is present because in this case the midpoint between any NN Fe sites is no longer an inversion symmetry center51. Due to the preservation of the mirror symmetry with the mirror plane perpendicular to the honeycomb layer passing through any two NN Fe sites, the NN DM interaction here only has the in-plane component51, as shown in Fig. 1c of the main text. The phonon part Hp considering only Fe2+ ions can be expressed as follows in the harmonic approximation:

$${H}_{{{{{{\rm{p}}}}}}}=\mathop{\sum}\limits_{i}\frac{{{{{{{\bf{P}}}}}}}_{i}^{2}}{2M}+\frac{1}{2}\mathop{\sum}\limits_{i,j}{{{{{{\bf{u}}}}}}}_{i}^{T}{K}_{ij}({{{{{{\bf{R}}}}}}}_{i}^{0}-{{{{{{\bf{R}}}}}}}_{j}^{0}){{{{{{\bf{u}}}}}}}_{j},$$
(2)

where M is the ion mass of Fe2+, Pi is the momentum of the ion at the site i, ui is the displacement of the ion i from its equilibrium position \({{{{{{\bf{R}}}}}}}_{i}^{0}\), and Kij is the dynamical matrix along the bond ij.

Magnon-phonon coupling can naturally arise in Eq. (1) when lattice vibrations are taken into account, because the exchange interaction Jij and the DM vector Dij depend on the positions of the ions i and j. In a collinear antiferromagnet with perpendicular easy-axis anisotropy, to the lowest order of ui and δSi = Si − 〈Si〉, where 〈Si〉 is the ground state expectation value of Si, the isotropic Heisenberg interaction leads to a cubic magnon-phonon coupling term15. This term alone is only capable of causing the softening and broadening of spin waves52. On the other hand, the anisotropic DM interaction gives rise to a quadratic magnon-phonon coupling term9,10,11, which can create a gap between the magnon and phonon bands, leading to the formation of magnon polarons. It is also worth noting that magnetoelastic coupling, which arises from single-ion magnetostriction and is generally present in magnets with crystalline anisotropy1,13,14, primarily couples out-of-plane phonons with in-plane magnons1,13,14. However, this mechanism contradicts the experimental observations where the main hybridizations occur between magnons and in-plane polarized phonons (Fig. 2a, b). Additionally, the calculated acoustic phonons with out-of-plane polarization have lower energy than the magnons, suggesting the negligible hybridization through magnetoelastic coupling (Supplementary Fig. 4). Hence, we consider the DM-interaction-induced magnon-phonon coupling as the dominant term in Fe2Mo3O8, and for simplicity, we neglect the effects of isotropic Heisenberg interactions and magnetoelastic coupling. The final expression for this coupling Hmp is given by,

$${H}_{{{{{{\rm{mp}}}}}}}=\frac{DS}{\left|{{{{{{\bf{R}}}}}}}_{ij}^{0}\right|}\mathop{\sum}\limits_{\langle ij\rangle }\left({{{{{{\bf{u}}}}}}}_{i}-{{{{{{\bf{u}}}}}}}_{j}\right)\left({\hat{I}}_{3}-{\hat{{{{{{\bf{R}}}}}}}}_{ij}^{0}{\hat{{{{{{\bf{R}}}}}}}}_{ij}^{0}\right)\left(\delta {{{{{{\bf{S}}}}}}}_{i}+\delta {{{{{{\bf{S}}}}}}}_{j}\right),$$
(3)

where \(D=\left|{{{{{{\bf{D}}}}}}}_{ij}\right|\) is the magnitude of the DM interaction, \(\left|{{{{{{\bf{R}}}}}}}_{ij}^{0}\right|\) is the bond length, S = 2 is the total electron spin of the Fe2+ ion, \({\hat{I}}_{3}\) is a 3 × 3 identity matrix, \({\hat{{{{{{\bf{R}}}}}}}}_{ij}^{0}\) is the unit vector along bond ij, and \({\hat{{{{{{\bf{R}}}}}}}}_{ij}^{0}{\hat{{{{{{\bf{R}}}}}}}}_{ij}^{0}\) is the Kronecker product between two \({\hat{{{{{{\bf{R}}}}}}}}_{ij}^{0}\) s. Note that for a collinear antiferromagnet with moments aligned along c axis, spins precess about the c axis. In this case, only phonons involving in-plane displacements can couple to magnons according to Eq. (3), because δSi only has the in-plane component. This can be significantly different from the magnetoelastic coupling mentioned earlier1,13,14.

In general, the dynamical matrix Kij at the bond ij has 6 (3) independent parameters for a three- (two-) dimensional system. For our two-dimensional effective model, it can be assumed that Kij only has the longitudinal components along the bond ij to reasonably reduce the number of tuning parameters. Then, the elastic potential energy of the phonons can be expressed as following by the Hooke’s law11,

$$\frac{1}{2}\mathop{\sum}\limits_{i < j}{k}_{ij}{\left[{\hat{{{{{{\bf{R}}}}}}}}_{ij}^{0}\cdot \left({{{{{{\bf{u}}}}}}}_{i}-{{{{{{\bf{u}}}}}}}_{j}\right)\right]}^{2}.$$
(4)

Here, kij is the longitudinal spring constant of the bond ij. To fit the phonon dispersions of Fe2Mo3O8, it is considered up to the fourth nearest neighbor (FNN), with the NN k1, NNN \({k}_{2}^{{{{{{\rm{t}}}}}}}\) and \({k}_{2}^{{{{{{\rm{o}}}}}}}\), TNN k3, and FNN k4.

By using the standard Holstein–Primakoff transformation, local spins can be mapped to canonical bosons to the leading orders: \({S}_{i}^{+}=\sqrt{2S}{a}_{i}\) and \({S}_{i}^{z}=S-{a}_{i}^{{{\dagger}} }{a}_{i}\) when i belongs to the Fet sites, or \({S}_{i}^{+}=\sqrt{2S}{b}_{i}^{{{\dagger}} }\) and \({S}_{i}^{z}={b}_{i}^{{{\dagger}} }{b}_{i}-S\) when i belongs to the Feo sites. Then, the Hamiltonian of the coupled system can be expressed in a generalized Bogoliubov-de Gennes (BdG) form as \(H=\frac{1}{2}{{{{{{\bf{X}}}}}}}_{{{{{{\bf{k}}}}}}}^{{{\dagger}} }\hat{H}({{{{{\bf{k}}}}}}){{{{{{\bf{X}}}}}}}_{{{{{{\bf{k}}}}}}}\), where \({{{{{{\bf{X}}}}}}}_{{{{{{\bf{k}}}}}}}={\left({a}_{{{{{{\bf{k}}}}}}},{b}_{{{{{{\bf{k}}}}}}},{a}_{-{{{{{\bf{k}}}}}}}^{{{\dagger}} },{b}_{-{{{{{\bf{k}}}}}}}^{{{\dagger}} },{{{{{{\bf{u}}}}}}}_{{{{{{\bf{k}}}}}}},{{{{{{\bf{P}}}}}}}_{{{{{{\bf{k}}}}}}}\right)}^{T}\). Here, uk(Pk) is the four-vector for the two-dimensional displacements (momenta) of the Fet and Feo sites. It is noted that their out-of-plane displacements and momenta are omitted because they are decoupled from the other degrees of freedom in our simple effective model. In this representation, the commutation relation of the Xk vector is

$$g\equiv \left[{{{{{{\bf{X}}}}}}}_{{{{{{\bf{k}}}}}}}^{{{\dagger}} },{{{{{{\bf{X}}}}}}}_{{{{{{\bf{k}}}}}}}\right]=\left[\begin{array}{cccc}{I}_{2}&&&\\ &-{I}_{2}&&\\ &&&-i{I}_{4}\\ &&i{I}_{4}&\end{array}\right].$$
(5)

The eigenvalues and eigenvectors of the coupled system satisfy

$$g\hat{H}({{{{{\bf{k}}}}}})\left|{\phi }_{n{{{{{\bf{k}}}}}}}\right\rangle={\sigma }_{nn}{E}_{n{{{{{\bf{k}}}}}}}\left|{\phi }_{n{{{{{\bf{k}}}}}}}\right\rangle,\langle {\phi }_{n{{{{{\bf{k}}}}}}}| g| {\phi }_{m{{{{{\bf{k}}}}}}}\rangle={\sigma }_{nm},$$
(6)

where σ = σz ⊗ I6 acts on the particle-hole space. In the BdG formulation of the Hamiltonian of the coupled system, the second half of the eigenstates are redundant to the first due to the artificial particle-hole symmetry. In Fig. 5 of the main text, only the lowest four independent energy bands are displayed to make comparison with the experimental results.

The Berry curvature Ωnk of the eigenvector \(\left|{\phi }_{n{{{{{\bf{k}}}}}}}\right\rangle\) can be defined as

$${{{\Omega }}}_{n{{{{{\bf{k}}}}}}}=i\langle {\nabla }_{{{{{{\bf{k}}}}}}}{\phi }_{n{{{{{\bf{k}}}}}}}| {g}^{-1}\times | {\nabla }_{{{{{{\bf{k}}}}}}}{\phi }_{n{{{{{\bf{k}}}}}}}\rangle .$$
(7)

Then, in the spirit of the momentum space discretization method53, the gauge-invariant Chern number and Berry curvatures can be computed from Eq. (7).

We note that our effective two-dimensional model used to obtain the results in Fig. 5 of the main text is a simplified model, but since the magnons of this system are perfectly two dimensional, and the out-of-plane phonons do not couple to the magnons at the leading order (Eq. (3)), it is able to capture the essence of the magnon polarons (Figs. 2 and 4 in the main text). Nevertheless, to account for the three dimensionality of the phonons, we have also developed a three-dimensional model (Supplementary Fig. 4).