THz-frequency magnon-phonon-polaritons in the collective strong-coupling regime

The interaction between an electromagnetic wave (light) and an atomic, a molecular, or a material mode is termed “strong coupling” if the rate of coherent energy transfer between the light and matter is faster than the irreversible decay of the light or the coherent material excitation. In the strong-coupling regime, the electromagnetic and material modes can no longer be treated as separate entities but rather form two hybrid modes called polaritons, with an energy difference between the modes given by the splitting energy $ℏΩ$⁠.^1–4 Polariton systems enable extensive optical control over material behavior and coherent information transfer between light and material degrees of freedom at a rate $Ω$⁠, yielding opportunities for both classical^5,6 and quantum information processing.^7–9 Earlier, the strong coupling of the electric field of light to the electric dipole moments of atoms was explored,¹⁰ and later, it was also demonstrated with superconducting qubits and excitons.^11–15 More recently, coupling to magnetic moments has been studied using magnetically active systems such as nitrogen-vacancy (NV) centers and magnon spin waves.^16–18 Cavity-enhanced versions of this coupling have also been observed, for example, between photons and a diverse range of excitations, including (Rydberg) atom transitions,^19–21 excitons,²² plasmons,²³ NV center spins,²⁴ magnons,²⁵ two-dimensional electron gases,²⁶ and Cooper pairs of a superconductor.¹¹ A subset of demonstrations has even reached the single-photon or few-photon level, which falls under the field of cavity quantum electrodynamics (QED).^8,9 While quantum demonstrations motivate cutting-edge applications such as quantum computation²⁷—for which dispersive coupling between a cavity and a two-level system might be used for repeated measurements of a qubit via a quantum nondemolition method^28,29—strong coupling can also be based on the classical coherence of macroscopic electrodynamic fields. The latter not only serves as important groundwork for its quantum analogs, but also has applications of its own right in fields such as lasing^11,12 and sensing.^13,14

Spintronics, in its quest toward long-range and high-bandwidth operation, would particularly benefit from such strong coupling because light can provide a means for facile transport and interaction with spin information. Moreover, cavity-enhanced versions of this coupling can be used to achieve long-term storage and all-optical addressing of spin states, important elements in realizing spin-based computation. Many of these possibilities already exist in the microwave frequency range. For example, in recent advances, magnons have been used as data carriers,^30,31 because, unlike spin-polarized electrons, they provide transfer of spin information over macroscopic distances without any Joule heating and enable access to wave-based computing concepts. Meanwhile, the first classical demonstration of cavity spintronics was reported by Zhang et al. using microwave cavity photons and magnons in yttrium iron garnet (YIG).¹⁸ They observed a Rabi-like oscillation of the light intensity, which they attributed to the coupling between the magnetic field of the EM wave and the magnetization of YIG.

The first report of polaritons formed by the strong coupling of terahertz (THz) light and magnons was by Sanders et al. in 1978,³² wherein a slab of iron (II) fluoride that exhibits an antiferromagnetic resonance at $1.58THz$ was used. More recently, Kampfrath et al. used advancements in the generation of ultrashort THz pulses to demonstrate control over the coherent superposition of antiferromagnetic magnons in nickel oxide.³³ There have also been demonstrations at THz frequencies of the inverse Faraday effect³⁴ and nonlinear spin control.^35–37 The prospect of harnessing these properties makes it desirable to realize a magnon-based spintronics platform at THz frequencies, but this has been difficult to achieve due to a scarcity of appropriate sources and detectors. Optical light can be used for generation and detection of THz magnons, but the interaction is either indirect^38,39 or nonlinear⁴⁰ and the efficiency is not high. Meanwhile, free-space THz light has been shown to exert direct linear control over magnons^33,41 but is not conducive to realizing miniaturized devices. Moreover, a cavity-enhanced version of this platform has been elusive due to the added constraint of requiring a highly confining THz cavity.

Here, we overcome these limitations using the THz polaritonics on-chip waveguide platform to demonstrate collective strong coupling of both magnons and phonons to THz light, with and without the utilization of a cavity structure. This platform has already been utilized to demonstrate capabilities such as spatiotemporal coherent control of lattice vibrational waves,⁴² time-resolved near-field THz imaging,⁴³ and on-chip all-dielectric devices.⁴⁴ Its great advantages are the high confinement and the efficient all-optical generation and detection of the polaritons on subcycle time scales and subwavelength length scales. Using this platform, optical phonons in a slab of ferroelectric lithium niobate (⁠$LiNbO3$⁠, LN) are strongly coupled to a THz electric field to form phonon-polaritons, which are simultaneously strongly coupled to magnons in an adjacent slab of canted antiferromagnetic erbium orthoferrite (⁠$ErFeO3$⁠, EFO) via the THz magnetic field. In addition, we also confine the THz waves in a photonic crystal cavity fabricated by patterning the hybrid LN–EFO slab with a two-dimensional array of air holes. In both patterned and unpatterned hybrid waveguide samples, the strong coupling leads to the formation of new modes we term magnon-phonon-polaritons, which we experimentally observe in the time domain as beating between the upper and lower polariton branch modes, corresponding in the frequency domain to avoided crossings. Coupling to the magnons in EFO forms the basis for a spintronic platform, while coupling to the phonons in LN enables the efficient generation and detection of light directly inside the cavity. Our simple yet versatile waveguide platform, thus, provides a promising avenue through which to explore ultrafast THz spintronics, quantum electrodynamics, sensing, and spectroscopic applications.

The four different sample structures under investigation are sketched in Fig. 1. For all samples, we started with a $53μm$ thick (100) LN MgO-doped stoichiometric melt single crystal grown by the Czochralski method and diced to dimensions of $11×10mm2$⁠. For the hybrid waveguide/cavity [see Figs. 1(b) and 1(d)], an additional $40μm$ thick (001) EFO slab, grown as a single crystal by the floating-zone method and cut to a size of $5×5mm2$⁠, was used.^45–47 Both LN and EFO crystals were polished using diamond paste (⁠$1μm$ grain size) to reach a surface roughness of less than $100nm$⁠, which is less than 1/10 of the THz wavelength. Large sample planes were prepared with a parallelism accuracy below $1μm$⁠. In case of the hybrid structures, a thin layer of air, approximately $8μm$ thick, separated the two slabs.

The bare LN slab acts as a waveguide and provides THz field confinement along the $z$-axis (out-of-plane axis), but for the cavity structures [see Figs. 1(a) and 1(c)], we also desire confinement in the $x$- and $y$-directions. To accomplish this, we utilized a photonic crystal cavity^48–50 (PhC) composed of a periodic array of air holes cut into the LN slab using femtosecond laser machining.⁵¹ The cavity we fabricated, often referred to as an L3 cavity,⁵² is composed of a hexagonal lattice of air holes with three holes absent from the array to form the “defect” cavity. To fabricate the hybrid LN-EFO cavity, we first attached the LN and EFO slabs together and then laser machined the cavity into the hybrid slab.

In LN, the material excitation is the polar lattice vibration (i.e., transverse optical phonon mode) with its polarization $P$ along the (ferroelectric) $c$ crystallographic axis of the tetragonal crystal (the $y$-axis as defined in Fig. 1). The electric-field component of a THz-frequency electromagnetic wave, oriented along the $y$-direction, coherently couples to these phonons to produce mixed phonon-polaritons.^42,53 Although the polaritons in our experiments are well detuned from the phonon resonance (the lowest transversal optical phonon resonance is located at $4.6THz$⁸⁸), the phononic nature is still integral to the dynamics. Not only does the lattice carry $∼32%$ of the energy in the $0.1$–$2THz$ range ( Appendix C), the phonon mode is also responsible mechanistically for efficient generation and detection of the THz waves in our experiments. In EFO, the relevant material excitation is the collective magnetic spin wave, i.e., the Brillouin zone center quasiantiferromagnetic (AF) magnon mode wherein the net magnetization $MAF$ along the $c$ axis of the orthorhombic crystal (the $z$-axis as defined in Fig. 1) is modulated in amplitude.^41,54 The magnetic-field component of THz light along the $z$-axis coherently couples to the AF mode at the magnetic resonance transition frequency of $0.67THz$ (with no applied magnetic field) to produce mixed magnon-polaritons. In the hybrid samples, uniform $y$-electric-polarized, $z$-magnetic-polarized THz-frequency electromagnetic waves extend throughout both materials, resulting in new hybridized magnon-phonon-polariton modes.

To study the interaction, we performed pump–probe measurements using a Ti:Sapphire amplifier system (Coherent Inc.) that outputs $90fs$ duration pulses with an 800 nm center wavelength at a 1 kHz repetition rate. The laser power was first split in a ratio of 90/10 for the pump and probe pulses, respectively. The pump pulse was variably delayed by time $t$⁠, frequency doubled to 400 nm in a $β$-barium borate (BBO) crystal, and as shown in Figs. 1(a) and 1(b) focused into the bare $LiNbO3$ portion of the slab to excite a $y$-electric-polarized, $x$-propagating phonon-polariton wavepacket via impulsive stimulated Raman scattering.⁵⁵ This process is also extensively exploited for high-field THz wave generation in LN.⁵⁶ Our phonon-polariton wavepackets have about $0.1$–$2THz$ spectral content. In the LN slab whose thickness (tens of $μm$⁠) is on the order of the generated THz wavelengths (few tens to few hundreds of $μm$⁠), the THz waves propagate in the plane of the slab as dielectric waveguide modes.⁵⁷ Due to the index mismatch between the phonon-polaritons and the optical pulse generating them, the phonon-polaritons also propagate laterally away from the generation region. To detect their progress, we exploit the electro-optic effect wherein the THz electric field modulates the LN refractive index.⁵⁸ More specifically, the lattice displacements of the optical phonon mode change the optical refractive index via electron–lattice interactions.⁵⁹ By passing a time-delayed optical probe pulse (typically at a different wavelength from the $400nm$ pump) through the slab and measuring the THz-induced depolarization, the time-dependent THz electric-field profile can be obtained [see Fig. 2(a)]. The field profile can then be Fourier transformed to retrieve the frequency spectrum. If the probe pulse is spatially expanded and detected on a camera as shown in Fig. 2(b), we can also obtain a spatially resolved image of the THz wave. We thus have a platform for generating and detecting THz waves directly inside the LN slab.

To first study the interaction between magnons and phonons in the bare slabs, we performed pump–probe measurements on the unpatterned samples shown in Figs. 1(a) and 1(b). The $400nm$ pump is cylindrically focused into the LN portion of the samples to excite $y$-electric-polarized, $x$-propagating phonon-polariton wavepackets, propagating away from the generation region. In the hybrid case, these waves enter into the hybrid LN–EFO slab. Here, the electric and magnetic fields of the THz wave, polarized along the $y$- and $z$-directions, respectively, drive the polarization of the polar phonon mode in LN and the magnetization of the AF magnon mode in EFO. In Fig. 3(a), we show a finite element method (FEM) simulation of electromagnetic wave transmission from the bare slab into the hybrid slab at $0.67THz$ in the absence of any resonances (detailed simulations can be found in  Appendix A). The transmission can be seen to proceed fairly efficiently, and a uniform electromagnetic mode can be seen extending through the slab(s), evanescently decaying away beyond either face of the structure. These waves are recorded via $800nm$ probe pulses that were shifted spatially (along the $x$-direction) and time-delayed with respect to the pump, enabling spatiotemporal field characterization. A display of the profiles $EyTHz(t)$ measured in successive sample locations appears as a 2D space-time matrix $EyTHz(x,t)$ as shown in Fig. 3(b). Although the pump launches counterpropagating waves, initially we only observe the rightward traveling wave because the pump pulse is to the left of the detection window. The first several oscillation cycles of this wave are due to the lowest-order waveguide mode, and the second oscillation cycles are due to the next waveguide mode. Sometime later, we also observe an additional rightward propagating wave starting at $≈40ps$ and a leftward propagating wave starting at $≈90ps$ that are due to partial reflections at the edges of the sample. Finally, a 2D Fourier transform was performed to obtain the wavevector-frequency dispersion curve $EyTHz(kx,ω)$⁠. In Fig. 3(c), we show the experimentally recorded dispersion curve for the bare LN slab [Fig. 1(a)]. Owing to the almost linear dispersion of the THz phonon-polaritons in our THz pulse bandwidth, two well separated transverse-electric (TE) dielectric phonon-polariton modes appear,⁵⁷ showing almost linear dispersion within our THz pulse bandwidth. In Fig. 3(d), we show the experimentally recorded dispersion curves for the hybrid LN-EFO sample where three well separated modes appear. The existence of three modes in Fig. 3(d) compared to the two modes in Fig. 3(c) is a result of the increased thickness of the hybrid slab. For the proceeding analysis, we make two simplifications. Firstly, the existence of an air gap can be ignored as it does not significantly influence the dispersion curves (see  Appendix A). Secondly, as demonstrated by the simulations, at most frequencies, the hybrid waveguide dispersion curves resemble those in a $93μm$ slab of LN. This is because of the very similar refractive indices of LN (⁠$n=5.0$⁠) and EFO (⁠$n=4.9$⁠) for THz fields with the selected polarization. With these assumptions, we can focus on the notable feature in Fig. 3(d) near the magnon resonance at $0.67THz$⁠: an avoided crossing in both the first and second-order waveguide modes.

The avoided crossing is a hallmark signature of strong coupling and can be explained by the formation of hybrid modes, the magnon-phonon-polaritons. In the following, we discuss this phenomena comprehensively by comparing a classical coupled oscillators’ model with a Lorentz model.

A phenomenological approach of coupled oscillators can be used to describe the interaction, using the Hamiltonian given as

Here, $Hpp$ and $Hm$ are the bare Hamiltonians and $κ=12GHz$ and $γ=8GHz$⁶⁰ are the FWHM linewidths of the uncoupled phonon-polariton and magnon modes, respectively. $Wint$ denotes their interaction, $c$ is the speed of light in vacuum, $ω$ is the eigenfrequency of the coupled system, $ωo/2π=0.67THz$ is the magnon resonance frequency, $β$ is the hybrid waveguide wavevector that acts as the detuning parameter between the two bare modes via the term $βc−ωo$⁠, and $Ω$ is the splitting frequency, which is a measure of the coupling strength. Note that $β$ can be defined in terms of the relation $β=neffω/c$⁠, where $neff$ is termed the effective index and is defined as the ratio of the speed of light to the phase velocity of the waveguide mode. $β$ thereby contains a frequency and a mode-dependent weighted average of the refractive indices of LN, EFO, and air. The solutions to this $2×2$ Hamiltonian correspond to the upper and lower magnon-phonon-polariton branches $ωUP$ and $ωLP$⁠, respectively.

In Fig. 4(a), we show an enlargement of the LN-EFO hybrid slab dispersion curve for the first-order waveguide mode near the magnon resonance [see Fig. 3(d)], in which the avoided crossing feature can be seen clearly. The bare magnon and phonon-polariton modes are only well fitted to the experimental dispersion at large detunings, where the hybrid modes resemble the bare modes [Fig. 4(a), dashed lines]. In contrast, the upper and lower magnon-phonon-polariton branches from Eq. (3) are well fitted at all detunings, reproducing the normal mode splitting [Fig. 4(a), solid lines]. As a complement, in Fig. 4(b), we show the frequency spectrum for the first-order mode at a series of wavevectors spaced by $0.4rad/mm$ and ranging from $61$ to $65rad/mm$⁠. The mode has phonon-polariton character at large detunings but shows a double-peaked structure at and near zero detuning where the magnon and phonon-polariton modes are strongly mixed. Although the second-order mode in Fig. 4(b) is far detuned for this range of wavevectors and thus shows a simple monotonic increase in center frequency, it also demonstrates a double-peaked spectrum at zero detuning. Similarly, Fig. 4(c) shows the frequency spectrum for the second-order mode at a series of wavevectors spaced by $0.4rad/mm$ and ranging from $48.5$ to $52.5rad/mm$⁠. To confirm that the avoided crossing was due to the AF-mode magnon in EFO, we also performed temperature-dependent experiments (data not presented). We clearly observed the avoided crossing shift to $0.75THz$ at $80K$⁠, in agreement with the shift of the AF-mode magnon frequency in EFO.⁶⁰ 

The splitting frequencies for the first and second-order waveguide modes were found to be $Ω1/2π=20GHz$ and $Ω2/2π=18GHz$⁠, respectively. Along with the linewidths for the magnon $γ/2π=8GHz$ and the phonon-polariton mode $κ/2π=12GHz$⁠, these allow us to estimate the cooperativity factor $C=Ω2/κγ$⁠. The cooperativity factor is a dimensionless quantity that compares the coupling strength to the losses in the system, with $C>1$ considered strong coupling. The first and second-order phonon-polariton modes couple to the magnons with cooperativity factors $C=4.0$ and $C=3.4$⁠, respectively, further validating that the system is in the strong-coupling regime.

So far, we have used a classical coupled oscillator model for magnons and phonons because the two can be modelled as damped harmonic oscillators. A Lorentz model (i.e., linear dispersion theory) can also be used to describe the polariton formation. In this case, we assume that the permeability of EFO is modelled as a damped harmonic oscillator near the magnon resonant frequency. As such, the dispersion relation for the electromagnetic wave can be written as

where we assume that $Δμk=8×10−4$⁠.⁶⁰ Using this approach, we obtained a dispersion curve as overlaid on the data in Fig. 4(a) (green dash-dotted lines). Outside the anomalous frequency range, the Lorentz model shows excellent agreement to the data and the coupled oscillator model of Eqs. (1)–(3). This is not surprising given that the Lorentz and coupled oscillator models are similar phenomenological approaches to describing strong light-matter coupling. The similarity comes from the fact that the electromagnetic field can be quantized as a harmonic oscillator or defined using dispersion relations. The discrepancy (anomalous dispersion vs an avoided crossing) comes from the fact that the Lorentz model allows for a complex wavevector, i.e., spatial damping. This damping is the reason that the anomalous region was not seen in our experiments. In fact, in the absence of damping (i.e., $γ∼0$⁠) and close to resonance (i.e., $ω∼ωo$⁠), the Lorentz model is equivalent to the coupled oscillator model and its result can be written as

This has the same form as Eq. (3) in the absence of damping and also models an avoided crossing. In this case, the splitting frequency is given by $Δμrωo/2π≃19GHz$⁠, which closely matches the experimentally extracted values for $Ω1/2π$ and $Ω2/2π$⁠.

Despite the convergence of these two models, we should point out that light-matter interactions that allow for light propagation always have an associated complex wavevector, and strong light-matter coupling in these systems always exhibit anomalous dispersion. The Lorentz model is, therefore, often used to describe polaritons in bulk and waveguide geometries, while the coupled oscillator model is more suitable when the coupling occurs in a photonic cavity, for example. Interestingly, while both models are classical, they can also be used to describe the quantum observation of vacuum Rabi splitting.^14,61 The ability of these classical models to describe a quantum interaction stems from the fact that for weak atomic excitations, the quantum theory of interaction between an electromagnetic field and an atom is equivalent to the Lorentz model for a classical oscillator. Thus, to earn the distinction of quantum strong coupling, it is not sufficient to observe a splitting. Instead, a dependence of splitting energy on the electric field or the number of oscillators must be observed.^4,62 Such a dependence was not observed in our experiments but was observed in recent measurements on other samples.⁶³ 

The strong coupling suggests that there is coherent energy exchange between the magnons and phonon-polaritons at a rate $Ω$⁠. Our spatially and temporally resolved measurements permit direct observation of this process. Figure 4(d) shows time-dependent data at a selected wavevector of $β=63rad/mm$⁠, obtained by frequency filtering of the $63rad/mm$ spectrum shown in Fig. 4(b) to exclude all frequencies outside of the $0.60$–$0.73THz$ range, then inverse Fourier transforming to return to the time domain. This allows us to examine the evolution of only the lower-order waveguide mode at frequencies around the avoided crossing in order to illustrate the energy exchange clearly. Note that if there were no strong coupling, then there would be no hybridization and ensuing energy splitting between upper and lower polariton branches as outlined by Eqs. (3) and (6). In this case, there would only ever be a single mode at each wavevector; even on resonance, the magnon and phonon-polariton modes would be degenerate and hence not show any beating. Yet, in Fig. 4(d), we inspect the modes near 63 rad/mm, i.e., on resonance, and the signal shows beating with energy exchange period $2π/Ω=50ps$ and decay time $τ=2/(κ+γ)=100ps$⁠. This is a strong indication of the energy exchange between the two phonon-polariton and magnon components of the mixed magnon-phonon-polariton mode.⁶⁴ 

Although the polariton bandgap in the hybrid waveguide excludes any eigenmodes around the bare magnon frequency, the phonon-polariton wavepacket entering from the LN extension includes frequencies within this forbidden gap. These waves decay exponentially into the hybrid slab. This can be explicitly seen by noting that the imaginary component of $β$ in Eq. (5) can be explicitly written as

The imaginary component of the wavevector corresponds to the inverse of the decay length into the slab. Ignoring other sources of absorption, the decay rate should be the greatest at the center of the gap (i.e., the bare magnon frequency) and the rate should decrease as frequencies at the gap edges are approached.⁵⁰ As shown in Fig. 4(e), the frequency component of $0.67THz$ in the center of the bandgap has a decay length of $0.48mm$⁠, 2–3 times shorter than those either above or below the bandgap, which have decay lengths of $1.27mm$ and $1.32mm$⁠, respectively. The frequency-dependent damping rates outside the forbidden range closely follow those of bare LN.⁶⁵ 

We would also like to note that there exist THz-frequency phonon modes in EFO itself, and it was demonstrated that driving the modes in the infrared (⁠$≈20THz$⁠) can create an effective magnetic field that subsequently drives the AF-mode magnon.³⁶ However, this is a nonresonant Raman-type effect using $≈MV/cm$ 20-THz field strengths. In contrast, we drive THz phonon-polariton waves in LN directly (primarily driving the ions through impulsive stimulated Raman scattering), and we use their magnetic-field component to resonantly drive the magnons in EFO. The driving of the magnons is thus linear in the THz field. In addition, in our system, the magnon response can be resolved through its coupling to the THz pulse at $≈V/cm$ field strengths owing to the efficient electro-optic detection in LN.

We use the term “strong coupling” in the manner broadly used in the literature to refer to strong light-matter coupling and the formation of polaritons, such as phonon-polaritons, plasmon-polaritons, exciton-polaritons, etc.^1,4 In this context, there is a clear distinction between simple linear absorption, i.e., weak coupling, and the strong-coupling regime, and the transition from one to the other is determined by the relative strength of the coupling compared to the losses in the system. In simple linear absorption, the coupling between the light and matter does not modify the dispersion except to reduce the light intensity at the absorption resonance, which in this case would be at the uncoupled magnon frequency. In the weak-coupling regime, the matter does indeed modify the dispersion but not in a way that produces two distinct coupled modes. In contrast, in the strong-coupling regime, the light and matter strongly interact and become hybridized modes at completely new frequencies, different from any uncoupled resonance frequency. We can use Eqs. (4) and (5) to demonstrate the difference between these regimes. In particular, in Fig. 4(f), we have plotted the dispersion relation of Eq. (5) for different ratios of magnon damping $γ$ to oscillator strength $Δμrωo2$⁠, (the value in experiments was $Δμrωo≃2.4$⁠). As the ratio is increased, linear dispersion is transformed into one approaching an avoided crossing. In reality, a true avoided crossing is only seen in the limit that the damping is small compared to the coupling/oscillator strength, such that Eq. (5) produces two solutions corresponding to the upper and lower hybridized polariton modes as specified by Eq. (6). Thus, strong coupling can be thought to occur when there are two solutions to Eq. (5) at a single wavevector with linewidths narrow enough to appear as distinct peaks, i.e., the birth of two hybrid modes. This transition occurs exactly at $Δμrωo/γ=1$⁠, as evidenced in Fig. 4(f). Any ratio below this value can practically be considered to describe the weak-coupling regime, since there is still a deviation from a linear dispersion, yet there is still only a single peak evident for each wavevector.

To further demonstrate the point that the features in Figs. 4(a)–4(c) are not simply “gaps” in the dispersion curve due to linear absorption, we have performed an experiment that demonstrates the case of linear absorption due to a magnetic resonance by measuring the dispersion curve after transmission through an array of split-ring resonators (SRRs) patterned onto the surface of a $30μm$ LN slab (confer Fig. 5). The SRRs are $200nm$ thick gold on a thin chromium adhesion layer deposited by optical lithography directly on the surface of the LN slab. Due to the orientation of the $E$- and $H$- fields in the experiment, the SRRs only exhibit a magnetic resonance that can be modelled using a Lorentz model for the magnetic permeability.⁶⁶ Thus, they permit close comparison to the AF-mode magnon in EFO. From the experimental dispersion curve of Fig. 5(c), there is no evidence of an avoided crossing. Instead, we observe a reduced amplitude around the resonant frequency. Furthermore, from the experimental wavevector-dependent spectra in Fig. 5(d), there is never a double-peaked structure. Instead, there is always a single peak whose magnitude is reduced when on resonance. This is the behavior expected for linear absorption and is demonstrably different from the avoided crossing shown in Fig. 4(a).

The bare LN waveguide [Fig. 1(a)] provides THz field confinement along the $z$-axis (out-of-plane axis), but to build a cavity, we also desire confinement in the $x$- and $y$-directions. To accomplish this, we utilized a photonic crystal cavity^48–50 composed of a periodic array of air holes cut into the LN slab using femtosecond laser machining.⁵¹ The cavity we fabricated, often referred to as an L3 cavity,⁵² is shown in Fig. 1(c) and is composed of a hexagonal lattice of air holes with three holes removed from the array to form the “defect” cavity. To excite the cavity modes, we focused pump pulses into the defect region of the crystal, which launched two counterpropagating THz phonon-polariton wavepackets. To detect the cavity polaritons, we directed probe pulses into the defect region and recorded the time evolution of the THz E-field as shown in Fig. 6(a) (imaging) and Fig. 6(b) (localized point detection). The wavepackets initially contain a broad range of frequencies as evidenced by a windowed FT [see Fig. 6(b), right]. When they reach the edges of the cavity, those waves with frequencies in the PhC bandgap are strongly reflected, while the rest of the frequency components are partially reflected due to the index mismatch with the surrounding PhC (i.e., Fresnel reflection). The bandwidth in the cavity thus narrows with every successive interaction with the holes, quickly approaching the limit where only the cavity modes are confined to the slab [see Fig. 6(b), right]. We also performed simulations for the E-field profiles within the cavity, which reveal that the different cavity modes have distinct spatial distributions. Since we use small focal spot sizes for the pump (⁠$wx≈10μm$⁠, $wy≈50μm$⁠) and probe (⁠$wx≈wy≈10μm$⁠) for the localized point detection measurements, their locations in the cavity determine the modes with which they predominantly interact. By choosing pump and probe to overlap with the extrema of a particular mode spatially, we can favor its generation and detection, respectively (see  Appendix B).

Thus far, we have discussed the central role that the phonon mode plays in THz wave generation and detection, but the phononic nature of the polaritons is most easily seen by its influence on the cavity modes as a function of temperature. Firstly, the transverse optical phonon in LN is a soft mode associated with the ferroelectric phase transition at $Tc=1418K$⁶⁸ at which the dielectric constant reaches a maximum. As a result, there is a sizable decrease in the THz refractive index from its room temperature values (⁠$neo≈5$ and $no≈6.7$⁠) to those at cryogenic temperatures.⁶⁹ As shown in Fig. 6(c), the decrease in the refractive index increases the cavity mode frequency (⁠$Δω/ω∼Δn/n$⁵⁰). Secondly, the damping of the low-frequency phonon-polaritons in LN is mainly attributed to a coupling of the TO phonon to other low-frequency phonon modes.^70,71 This damping is reduced by a factor of $≈$5 from $290K$ to $79K$ as the low-frequency phonon modes become depopulated.^69,72 The effect on the cavity modes can be deduced from Fig. 6(c), where the mode linewidth is also reduced by a factor of $≈$ 5. To analyze this change more quantitatively, we use a common measure of the cavity lifetime called the quality factor $Q=ωcav/κ$⁠, where $ωcav$ is the resonant frequency of the mode and $1/κ$ is the rate of energy loss from the cavity.⁷³ The quality factor can be interpreted as a dimensionless lifetime that can account for the combined contributions from the intrinsic damping $Qm$ and radiative scattering of light out of the cavity $Qr$ leading to a net $Qnet$ defined as $1/Qnet=1/Qm+1/Qr$⁠. In our case, the increase in $Qnet$ with decreasing temperature is directly proportional to the decrease in phonon damping since this is the dominant mechanism limiting the lifetime, i.e., $Qm≪Qr$⁠, so $Qnet≃Qm$⁠. Hence, in our experiments, the highest $Q=510$ was obtained at our lowest recorded temperature of $79K$⁠. Based on the temperature dependence of the phonon absorption coefficient of LN, the $Q$ can be increased further to $Q≃1000$ at 4K.⁷² Another important measure of the cavity is the effective mode volume $Veff$⁠, defined as

For the so-called X1 mode ( Appendix B), $Veff=3.4×10−3λ3≃0.5(λ/n)3$⁠, where $λ$ is the free-space wavelength and $n=5.21$ is the average refractive index of the slab. Although the effective mode volume is not appreciably smaller than the $(λ/n)3$⁠, the large refractive index of the LN slab at THz frequencies results in a mode volume that is nearly three orders of magnitude smaller than the free-space volume.

Having established the particulars of our LN based cavity, we now discuss our experiments regarding the collective strong coupling between magnons, phonons, and cavity photons in a hybrid L3 PhC cavity, depicted in Fig. 1(d). Our system consists again of a thin composite slab of $53μm$ thick (100) LN and $40μm$ thick (001) EFO. As before, the PhC cavity was fabricated using laser machining to cut hexagonal array of air holes through the hybrid slab. The E-field component of the THz light, polarized along the $z$-direction, can couple to the phonons in LN to produce mixed phonon-polariton modes. Simultaneously, the H-field component of light, polarized along the $z$-axis, can couple to the AF mode of EFO to produce mixed magnon-polaritons. In the hybrid cavity, a contiguous $y$-electric-polarized, $z$-magnetic-polarized THz-frequency electromagnetic wave can extend throughout both materials, forming new hybridized cavity magnon-phonon-polariton modes (see  Appendix B for hybrid slab FEM simulations).

To verify the formation of cavity magnon-phonon-polariton modes, we studied the frequency spectrum of the THz waves as the frequencies of the cavity and magnon modes were tuned by varying the temperature. While we have seen that the cavity resonance of the pure LN shifts upon reduction of the temperature from $290K$ to $77K$ by $≈18GHz$ in a weakly quadratic fashion [see the inset of Fig. 6(c)], the magnon resonance shifts by $≈80GHz$ and does so with a higher-order functional form.^60,74 To record the frequency spectrum, we generated and detected THz waves in the LN portion of the hybrid cavity [see Fig. 7(a)]. Despite being excited in LN, the THz waves evolve into the modes that span the entire hybrid structure. Thus, there exists an E-field parallel to the polar lattice vibration in LN and an H-field parallel to the antiferromagnetic magnetization in EFO. However, electro-optic probing is only sensitive to the E-field in LN and can, therefore, only detect the AF-mode magnon response through its interaction with the cavity mode. As shown in Fig. 7(b), at $290K$ we only observed a single Lorentzian peak corresponding to the cavity mode at $0.72THz$ and did not observe the magnon mode at $0.67THz$⁠, which was sufficiently detuned such that any interaction with the cavity mode was negligible. However, as the temperature was reduced and thus the two modes were tuned to be in resonance, a clear double-peaked spectrum appeared. At $210K$⁠, the detuning between the bare modes was only $2GHz$⁠, yet we observed two peaks that were separated by nearly eight times the detuning. This double-peaked structure is characteristic of a normal-mode splitting that can arise from the strong coupling between the cavity and magnon modes.

Like in the unpatterned case, the interaction can be modeled using the Hamiltonian for coupled oscillators given by

Here, $Hcav$ and $Ho$ are the bare Hamiltonians for, respectively, the cavity and magnon modes, and $Wint$ denotes their interaction. Correspondingly, $ωo/2π$ and $γcav$ are the magnon resonance frequency and linewidth, and $ωcav$ and $κcav$ are the cavity resonant frequency and linewidth, respectively. Finally, $Ωcav$ is the splitting frequency and the only undetermined parameter in the equations. The solutions to this Hamiltonian correspond to the upper and lower magnon-phonon-polaritons $ωUP$ and $ωLP$⁠, respectively. As shown in Fig. 7(b) by the blue dashed curve, the coupled oscillator model (SC) shows excellent agreement with our data and indicates that the phonon, magnon, and cavity photons are no longer the eigenmodes of the system. Instead, the THz response leads to collective mode oscillation at the frequency $Ωcav$ between a state wherein the energy is stored purely in the light field and phonons in LN and a state wherein the energy is stored purely in the magnons in EFO. Note that these spectra are distinctly different from the result of simple THz absorption by the magnon resonance in EFO. The magnon linewidth is far narrower than the splitting we observe between the upper and the lower polariton peaks, as shown in Fig. 7(b). Absorption at the magnon frequency would also yield a double-peaked structure around that frequency due to the magnon FID induced in the weak-coupling (WC) limit by the THz field generated in LN. However, in this case the separation between the two peaks would be given by the magnon absorption linewidth [simulation shown in the red dotted curve of Fig. 7(b)] rather than the strength of the coupling. The actual spectrum we measure and the model based on strong coupling (blue dashed curve) show a far larger separation between the two THz spectral peaks.

By identifying the lower and higher-frequency peaks as, respectively, the upper and lower polaritons, we generated a plot for the two branches as a function of the temperature, as shown in Fig. 7(c). The curves for the coupled oscillator modes and bare cavity and magnon modes are also overlaid in Fig. 7(c). The temperature dependence of the bare magnon frequency was recovered from the minima in the frequency spectra [Fig. 7(c), bottom left inset] and agrees well with previous results.^60,74 The temperature dependence of the bare cavity frequencies was determined from a second cavity mode that showed the same dependence [Fig. 7(c), top right inset]. The bare frequency values were used to calculate the upper and lower polariton branches with $Ωcav$ as the fitting parameter. The coupled oscillator model agrees well with the data at all detunings and reveals the existence of an avoided crossing between the two branches, a hallmark signature of strong coupling. From the two branches, we determined the splitting frequency at zero detuning to be $Ωcav/2π=16GHz$⁠. Utilizing the splitting frequency and the linewidths for the magnon $γcav/2π=7GHz$ and cavity mode $κcav/2π=14GHz$ at zero detuning, we also estimated the cooperativity factor $Ccav=2.6$ for our data. This further supports the assertion that we reached the collective strong-coupling regime.

These results agree well with our above presented results for the hybrid slab without the PhC cavity (⁠$Ω/2π=18GHz$⁠; $γ/2π=8GHz$⁠; $κ/2π=12GHz$⁠; $C=4.0$⁠). The slight differences arise due to the different temperatures and pump intensities, since we slightly reduced the pump intensity to avoid damage to the hybrid PhC structure. At first glance, this agreement might come as a surprise since one expects an increase in coupling strength due to the cavity. However, the collective normal-mode splitting $Ωcav$ of the magnon-phonon-polariton in the hybrid cavity is scaling with the square root of the ratio of the number of oscillators $N$ and the mode volume $V$⁠, i.e., $Ω∝N/V$⁠.^3,75 Indeed, by fabricating the cavity into the hybrid slab of EFO and LN, we decrease $V$⁠, but since we are in the solid phase, we also decrease $N$ by the same factor. Hence, we expect the same coupling strength as in the unpatterned hybrid slab. In addition, note that the quality factor of the hybrid cavity mode ranged from $Q=42$ at room temperature to $Q=56$ at $100K$⁠. This was smaller than in our bare LN cavity and was attributed to both damping from the EFO and significant radiative scattering. The radiative scattering is not entirely surprising given the difficulty of laser machining high-quality air holes through a composite structure. Nonetheless, our Q was sufficient to reach the collective strong-coupling regime, i.e., $C>1$⁠.

We have demonstrated evidence for collective strong coupling among magnons, phonons, and photons at THz frequencies using a unique on-chip platform. For the unpatterned LN-EFO hybrid slab, we have observed an avoided crossing in the dispersion curve for both the first- and the second-order TE dielectric waveguide modes, which exhibit a normal-mode splitting of $Ω1/2π=20GHz$ and $Ω2/2π=18GHz$ with corresponding cooperativity factors $C=4.0$ and $C=3.4$⁠, respectively. Correspondingly, in the time domain, we observed a beating between the upper and the lower magnon-phonon-polaritons with an oscillation period of $2π/Ω=50ps$⁠, which is the period for energy exchange between phonons and magnons. In the spatial domain, we observed evanescent wave decay at frequencies within the bandgap as the wavepacket entered the hybrid waveguide. For the bare LN cavity, we found strongly coupled admixtures of phonons and cavity photons exhibiting a Q-factor of $510$ at $79K$ (the lowest temperature studied) and a mode volume of $V=3.4×10−3λ3≃0.5(λ/n)3μm3$⁠. The Q-factor is limited by the phononic damping and could be increased to $Q≃1000$ at $4K$⁠. Furthermore, we could demonstrate collective strong coupling in a hybrid cavity composed of LN and EFO. Using temperature tuning of the cavity and magnon modes, we observed a normal-mode splitting of $16GHz$⁠, an avoided crossing between the upper and the lower polariton branches, and a cooperativity factor of $C=2.6$⁠, confirming that we reached the collective strong-coupling regime. The data also show excellent agreement with a coupled oscillator model. These observations are understood as evidence for the formation of cavity magnon-phonon-polaritons, wherein phonons in the LN and magnons in the EFO are coupled to the EM wave E- and H-fields, respectively, as they are in the unpatterned hybrid waveguide.

Our demonstration illustrates a system with multifunctional capabilities that should prove promising for realizing ultrafast THz cavity spintronics. For instance, EFO enables long-term storage and all-optical addressing of spin states. As a complement, LN has large electro-optic constants that enable efficient generation and detection of THz waves,^76,77 with prospects for both spatial and temporal shaping of the THz field profile,^42,53 e.g., via chirp-and-delay of the generating femtosecond laser pulse,⁷⁸ and prospects for various forms of electrical control.⁷⁹ The coupling to THz light in a waveguide also allows coherent propagation of the magnetization in a highly confined and efficient manner. Our cavity is able to achieve a Q-factor comparable to that of another THz cavity recently used in strong-coupling work,²⁶ but in a much smaller mode volume that is nearly diffraction limited. The femtosecond laser machining process could be used to fabricate other cavity geometries with considerably smaller model volumes.⁸⁰ Another useful feature is that the Q-factor can be increased without degrading the THz signal injection and readout, since these are achieved via optical-frequency light directly inside the cavity. We believe our approach may prove valuable for experiments pursuing THz cavity QED. These experiments require high $Q$ and low $V$ because the vacuum Rabi frequency must be fast compared with the cavity loss rate and decoherence processes. A high $Q$ minimizes decoherence, while a lower $V$ increases the Rabi frequency (⁠$Ω∼1/V$⁠). There has already been some progress in using electro-optic crystals similar to LN to measure the quantum properties of light at THz frequencies.^81–83 Combining this methodology with our cavity should enable the demonstration of THz cavity QED phenomena in the future.

The on-chip cavity platform might prove useful also for other applications like sensing and spectroscopy, e.g., by depositing liquid or solid samples into unpatterned LN slabs or cavities. The dual-waveguide geometry used here is one example of how a sample could be coupled to a LN waveguide for linear or nonlinear THz spectroscopy with no free-space THz propagation, relying on LN for THz generation and detection. While the hybrid cavity demonstrated here confined the THz polariton field in all three dimensions, cavities with one or two unconfined dimensions can also be created by fabricating two photonic structures with a region between them whose width could be an $(n+1/2)$ multiple of the desired cavity resonance wavelength, with the overall thickness of the structure dictating the extent of confinement in the $z$ dimension through index contrast with the surrounding free space. The dispersion properties of such cavities closely resemble those of the microcavities used in studies of exciton-polaritons.^84,85 It is straightforward to see that the effective mass of magnon-phonon-polaritons in such cavities can be as small as $10−7me$⁠, about two orders of magnitude lower than the typical value for exciton-polaritons. The corresponding de Broglie wavelength $λth=h/2πmkBT$ can be approximately $10μm$ at room temperature and $>50μm$ at low but not ultralow temperatures, $≈10K$⁠. This indicates that magnon-polariton Bose–Einstein condensation should be achievable at quite low excitation densities, with interesting prospects for fundamental study and possible applications in quantum simulation, low-power coherent THz sources, and other areas.

The optical measurements planned and conducted by P.S., A.S., B.D., J.L., and K.A.N. at MIT were supported as part of the Center for Excitonics, an Energy Frontier Research Center (EFRC), funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award No. de-sc0001088, and with support from NSF Grant No. CHE-1665383. The crystal growth and preparation by M.X., W.R., and S.C. at Shanghai University and S.K. at the Academy of Sciences of the Czech Republic were, respectively, supported by the National Natural Science Foundation of China (NNSFC) (Grant Nos. 51372149 and 51672171) and the Czech Science Foundation (Project No. 15-08389S). P.S. was supported in part by a Postgraduate Scholarship (PGSD) from the Natural Sciences and Engineering Research Council of Canada (NSERC). A.S. was supported by a scholarship from the German National Academy of Sciences Leopoldina (No. LPDS 2016-02).

In the main text, we show a finite element method (FEM) simulation of electromagnetic wave transmission from the bare LN slab into the hybrid LN-EFO slab at $0.67THz$ in the absence of any resonances [Fig. 3(a)]. Correspondingly, in Figs. 8(a) and 8(b), we show the FEM calculated field profile for the first-order waveguide modes in bare and hybrid structures at $0.67THz$⁠, respectively. In both the single and hybrid waveguides, a uniform electromagnetic mode can be seen extending through the slab(s) and evanescently decaying away on either side of the structure, illustrating that both samples are acting as THz waveguides. Furthermore, the existence of three modes for the hybrid slab in comparison with two modes for the pure LN slab is a result of the increased thickness as demonstrated in Fig. 8(c), which shows that the waveguide modes become lowered in frequency and more closely spaced as the thickness of the waveguide is increased. We also discuss in the main text two simplifications for the analysis of the dispersion curves. Firstly, the existence of an air gap can be ignored as it does not significantly influence the dispersion curves, as can be seen in Fig. 8(d). Secondly, as demonstrated by Fig. 8(e), at most frequencies, the hybrid waveguide dispersion curves resemble those in a $93μm$ slab of LN, due to the very similar refractive indices of LN (⁠$n=5.0$⁠) and EFO (⁠$n=4.9$⁠).

We also performed FEM simulations for the cavity structures. In Fig. 9(a), we show the simulated E-field profiles for the two cavity modes, labeled X1 and X2, in the bare LN PhC cavity. In general, the cavity modes have distinct spatial distributions, with the higher-frequency modes exhibiting a higher number of nodes and extrema.⁸⁶ Correspondingly, the simulated $Ex$-field and $Hz$-field profiles for the X1 cavity mode are shown in Fig. 9(b). The E-field component of the light, polarized along the $z$-direction, can couple to the phonons in LN to produce mixed phonon-polariton modes. Simultaneously, the H-field component of light, polarized along the $z$-axis, can couple to the AF mode of EFO to produce mixed magnon-polaritons.

Due to the small focal spot size of the pump (⁠$wy≈10μm$⁠, $wz≈50μm$⁠) and probe pulses (⁠$wy≈wz≈10μm$⁠) in our experiments, their locations in the cavity determine the modes with which they predominantly interact. By choosing pump and probe pulses to overlap spatially with the extrema of a particular mode, we can favor its generation and detection, respectively. For instance, in Fig. 9(c), we show the frequency spectrum where we selectively excited and detected either the X1 or the X2 mode.

To calculate the energy distribution between the lattice motion and the electromagnetic field, we begin by describing the coupling between the transverse polar optic phonon mode and the electromagnetic field. The coupled equations are given as⁸⁷ 

Here, $P$ is the macroscopic polarization of the material and $Q$ is the normalized ionic displacement of the polar optic phonon mode where an over-dot indicates the time derivative. In addition, $ωTO$ is the transverse optic phonon frequency, $ϵo$ is the vacuum permittivity, and $ε∞$ and $εo$ are the high- and low-frequency dielectric constants of the material, respectively. Poynting’s theorem relates the energy stored in the electromagnetic field to the work done on the electric dipoles in the medium. In differential form, the energy balance can be expressed as

where $S$ is the Poynting vector and $W$ is the energy density. Equation (C3) states that the rate of energy decrease in the medium is equal to the amount of energy flow out of the medium. To determine the energy distribution between the lattice and the electromagnetic field, we must determine an appropriate energy density that satisfies Eq. (C3). Such an energy density was proposed by Huang⁸⁷ in the absence of damping (⁠$Γ=0$⁠), which is given as

where $E$ and $H$ are the electric and magnetic fields, respectively, $Q$ is the normalized displacement of the phonon mode, and $μo$ is the vacuum permeability. The first term in brackets of Eq. (C4) represents the vibrational energy of the TO phonon mode, given as a sum of the harmonic oscillator and potential and kinetic energy terms. The second term in brackets represents the electromagnetic energy, with an inclusion of the electronic response via the term $ε∞$⁠. As such, we can write the time-averaged energies in the lattice and electromagnetic wave as

The time-averaged fraction of mechanical energy in the mode is then given as

where $⟨Wlatt⟩$ and $⟨WEM⟩$ are the time averaged lattice and electromagnetic energies, respectively. Using Eqs. (C5) and (C6) and assuming material parameters for LN as specified in the work of Feurer et al., ⁵³ the fraction of energy in the lattice is 32% at 1 THz.