### Nanofabrication of the devices

High-quality α-MoO_{3} flakes were mechanically exfoliated from bulk crystals synthesized by the chemical vapour deposition (CVD) method^{19} and then transferred onto either commercial 300 nm SiO_{2}/500 μm Si wafers (SVM) or gold substrates using a deterministic dry transfer process with a polydimethylsiloxane (PDMS) stamp. CVD-grown monolayer graphene on copper foil was transferred onto the α-MoO_{3} samples using the poly(methyl methacrylate) (PMMA)-assisted method following our previous report^{44}.

The launching efficiency of the resonant antenna is mainly determined by its geometry, together with a trade-off between the optimum size and illumination frequency^{45,46}. We designed the gold antenna with a length of 3 μm and a thickness of 50 nm, which provided a high launching efficiency over the spectral range from 890 to 950 cm^{−1} within the α-MoO_{3} reststrahlen band. Alternatively, a thicker antenna with a stronger *z* component of the electric field could be used to launch the polaritons more efficiently in future studies. Note that narrow antennas (50-nm width) were used to prevent their shapes from affecting the polariton wavefronts, especially when their propagation is canalized (such as in Figs. 2 and 3), while wider antennas (250-nm width) were used to obtain a higher launching efficiency and observe polaritons propagating across the SiO_{2}–Au interface in our experiments (for example, Fig. 4).

Gold antenna arrays were patterned on selected α-MoO_{3} flakes using 100 kV electron-beam lithography (Vistec 5000+ES) on an approximately 350 nm PMMA950K lithography resist. Electron-beam evaporation was subsequently used to deposit 5 nm Ti and 50 nm Au in a vacuum chamber at a pressure of <5 × 10^{−6} torr to fabricate the Au antennas. Electron-beam evaporation was also used to deposit a 60-nm-thick gold film onto a low-doped Si substrate. To remove any residual organic material, samples were immersed in a hot acetone bath at 80 °C for 25 min and then subjected to gentle rinsing with isopropyl alcohol (IPA) for 3 min, followed by drying with nitrogen gas and thermal baking (for more details on the fabrication and characterization of the Au–SiO_{2}–Au in-plane sandwich structure, see Supplementary Figs. 26 and 27).

The samples were annealed in a vacuum to remove most of the dopants from the wet transfer process and then transferred to a chamber filled with NO_{2} gas to achieve different doping levels by surface adsorption of gas molecules^{47}. The graphene Fermi energy could be controlled by varying the gas concentration and doping time, achieving values as high as ~0.7 eV (Supplementary Fig. 9). This gas-doping method provides excellent uniformity, reversibility and stability. Indeed, Raman mapping of a gas-doped graphene sample demonstrated the high uniformity of the method (Supplementary Fig. 9). As the deposition of NO_{2} gas molecules on the graphene surface occurs by physical adsorption, the topological transition of hybrid polaritons in graphene/α-MoO_{3} heterostructures can be reversed by controlling the gas doping. For example, after gas doping, the Fermi energy of graphene could be lowered from 0.7 to 0 eV by vacuum annealing at 150 °C for 2 h. The sample could subsequently be re-doped to reach another on-demand Fermi energy (Supplementary Fig. 28). It should be noted that the graphene Fermi energy only decreases from 0.7 to 0.6 eV after being left for 2 weeks under ambient conditions, which demonstrates the high stability of the doping effect (Supplementary Fig. 29). Note that chemical doping has been demonstrated to be an effective way to tune the characteristics of polaritons, such as their strength and in-plane wavelength^{48,49,50,51}.

### Scanning near-field optical microscopy measurements

A scattering scanning near-field optical microscope (Neaspec) equipped with a wavelength-tunable quantum cascade laser (890–2,000 cm^{−1}) was used to image optical near fields. The atomic force microscopy (AFM) tip of the microscope was coated with gold, resulting in an apex radius of ~25 nm (NanoWorld), and the tip-tapping frequency and amplitudes were set to ~270 kHz and ~30–50 nm, respectively. The laser beam was directed towards the AFM tip, with lateral spot sizes of ~25 μm under the tip, which were sufficient to cover the antennas as well as a large area of the graphene/α-MoO_{3} samples. Third-order harmonic demodulation was applied to the near-field amplitude images to strongly suppress background noise.

In our experiments, the p-polarized plane-wave illumination (electric field **E**_{inc}) impinged at an angle of 60° relative to the tip axis^{52}. To avert any effects caused by the light polarization direction relative to the crystallographic orientation of α-MoO_{3}, which is optically anisotropic, the in-plane projection of the polarization vector coincided with the *x* direction ([100] crystal axis) of α-MoO_{3} (Supplementary Fig. 6). Supplementary Fig. 30 shows the method used to extract antenna-launched hybrid polaritons from the complex background signals observed when the polaritons propagate across a Au–SiO_{2}–Au in-plane structure to realize partial focusing.

### Calculation of polariton dispersion and isofrequency dispersion contours (IFCs) of hybrid plasmon–phonon polaritons

The transfer matrix method was adopted to calculate the dispersion and IFCs of hybrid plasmon–phonon polaritons in graphene/α-MoO_{3} heterostructures. Our theoretical model was based on a three-layer structure: layer 1 (*z* > 0, air) is a cover layer, layer 2 (0 > *z* > –*d*_{h}, graphene/α-MoO_{3}) is an intermediate layer and layer 3 (*z* < –*d*_{h}, SiO_{2} or Au) is a substrate where *z* is the value of the vertical axis and *d*_{h} is the thickness of α-MoO_{3} (Supplementary Fig. 31). Each layer was regarded as a homogeneous material represented by the corresponding dielectric tensor. The air and substrate layers were modelled by isotropic tensors diag{*ε*_{a,s}} (ref. ^{53}). The α-MoO_{3} film was modelled by an anisotropic diagonal tensor diag{*ε*_{x}, *ε*_{y}, *ε*_{z}}, where *ε*_{x}, *ε*_{y} and *ε*_{z} are the permittivity components along the *x*, *y* and *z* axes, respectively. Also, monolayer graphene was located on top of α-MoO_{3} at *z* = 0 and described as a zero-thickness current layer characterized by a frequency-dependent surface conductivity taken from the local random-phase approximation model^{54,55}:

$$\begin{array}{rcl} {\sigma \left( \omega \right)}& = &{\frac{{i{{\rm{e}}^2}{k_{\rm{B}}}T}}{{{{\uppi}}{\hbar ^2}\left( {\omega + \frac i \tau} \right)}}\left[ {\frac{{{E_{\rm{F}}}}}{{{k_{\rm{B}}}T}} + 2\ln \left( {{{\rm{e}}^{ – \frac{{{E_{\rm{F}}}}}{{{k_{\rm{B}}}T}}}} + 1} \right)} \right]}\\{}&{}&{ + i\frac{{{{\rm{e}}^2}}}{{4{{\uppi}}\hbar }}\ln \left[ {\frac{{2\left| {{E_{\rm{F}}}} \right| – \hbar \left( {\omega + \frac i \tau} \right)}}{{2\left| {{E_{\rm{F}}}} \right| + \hbar \left( {\omega + \frac i \tau} \right)}}} \right]}\end{array}$$

(1)

which depends on the Fermi energy *E*_{F}, the inelastic relaxation time *τ* and the temperature T; the relaxation time is expressed in terms of the graphene Fermi velocity *v*_{F} = *c*/300 and the carrier mobility *μ*, with \(\tau = \mu E_{\mathrm{F}}/ev_{\mathrm{F}}^2\); *e* is the elementary charge; *k*_{B} is the Boltzmann constant; *ℏ* is the reduced Planck constant; and *ω* is the illumination frequency.

Given the strong field confinement produced by the structure under consideration, we only needed to consider transverse magnetic (TM) modes, because transverse electric (TE) components contribute negligibly. The corresponding p-polarization Fresnel reflection coefficient *r*_{p} of the three-layer system admits the analytical expression

$$\begin{array}{*{20}{c}} {r_{\mathrm{p}} = \frac{{r_{12} + r_{23}\left( {1 – r_{12} – r_{21}} \right){\mathrm{e}}^{i2k_z^{\left( 2 \right)}d_{\mathrm{h}}}}}{{1 + r_{12}r_{23}{\mathrm{e}}^{i2k_z^{\left( 2 \right)}d_{\mathrm{h}}}}},} \end{array}$$

(2)

$$\begin{array}{*{20}{c}} {r_{12} = \frac{{{{Q}}_1 – {{Q}}_2 + SQ_1Q_2}}{{{{Q}}_1 + Q_2 + SQ_1Q_2}},} \end{array}$$

(3)

$$\begin{array}{*{20}{c}} {r_{21} = \frac{{{{Q}}_2 – {{Q}}_1 + SQ_1Q_2}}{{{{Q}}_2 + Q_1 + SQ_1Q_2}},} \end{array}$$

(4)

$$\begin{array}{*{20}{c}} {r_{23} = \frac{{Q_2 – Q_3}}{{Q_2 + Q_3}},} \end{array}$$

(5)

$$\begin{array}{*{20}{c}}where {Q_j = \frac{{k_z^{\left( j \right)}}}{{{\it{\epsilon }}_t^{(j)}}},} \end{array}$$

(6)

$$\begin{array}{*{20}{c}} {S = \frac{{\sigma Z_0}}{\omega }.} \end{array}$$

(7)

Here, *r*_{jk} denotes the reflection coefficient of the *j*–*k* interface for illumination from medium *j*, with *j*,*k* = 1–3; \({\it{\epsilon }}_t^{(j)}\) is the tangential component of the in-plane dielectric function of layer *j* for a propagation wave vector *k*_{p}(*θ*) (where *θ* is the angle relative to the *x* axis), which can be expressed as \({\it{\epsilon }}_t^{(j)} = {\it{\epsilon }}_x^{(j)}\mathop {{\cos }}\nolimits^2 \theta + {\it{\epsilon }}_y^{(j)}\mathop {{\sin }}\nolimits^2 \theta\) (where \({\it{\epsilon }}_x^{(j)}\) and \({\it{\epsilon }}_y^{(j)}\) are the diagonal dielectric tensor components of layer *j* along the *x* and *y* axes, respectively); \(k_z^{\left( j \right)} = \sqrt {\varepsilon _t^{\left( j \right)}\frac{{\omega ^2}}{{c^2}} – \frac{{\varepsilon _t^{\left( j \right)}}}{{\varepsilon _z^{\left( j \right)}}}q^2}\) is the out-of-plane wave vector, with \({\it{\epsilon }}_z^{(j)}\) being the dielectric function of layer *j* along the *z* axis; and *Z*_{0} is the vacuum impedance.

We find the polariton dispersion relation *q*(*ω*,*θ*) when the denominator of equation (2) is zero:

$$\begin{array}{*{20}{c}} {1 + r_{12}r_{23}{\mathrm{e}}^{i2k_z^{\left( 2 \right)}d_{\mathrm{h}}} = 0.} \end{array}$$

(8)

For simplicity, we considered a system with small dissipation, so that the maxima of Im{*r*_{p}} (see colour plots in Supplementary Figure 10) approximately solve the condition given by equation (8), and therefore produce the sought-after dispersion relation *q*(*ω*,*θ*) (see additional discussion in Supplementary Note 1).

### Electromagnetic simulations

The electromagnetic fields around the antennas were calculated by a finite-elements method using the COMSOL package. In our experiments, both tip and antenna launching were investigated. For the former, the sharp metallic tip was illuminated by an incident laser beam. The tip acted as a vertical optical antenna, converting the incident light into a strongly confined near field below the tip apex, which can be regarded as a vertically oriented point dipole located at the tip apex. This localized near field provided the necessary momentum to excite polaritons. Consequently, we modelled the tip as a vertical *z*-oriented oscillating dipole in our simulations (Fig. 1c,f), a procedure that has been widely used for tip-launched polaritons in vdW materials^{56}. For the antenna launching, the gold antenna can provide strong near fields of opposite polarity at the two endpoints, thus delivering high-momentum near-field components that match the wave vector of the polaritons and excite propagating modes in the graphene/α-MoO_{3} heterostructure^{45,46}. Our simulations of polariton excitation by means of antennas, such as in Fig. 3b,d, incorporated the same geometrical design as in the experimental structures.

We also used a dipole polarized along the *z* direction to launch polaritons, and the distance between the dipole and the uppermost surface of the sample was set to 100 nm. We obtained the distribution of the real part of the out-of-plane electric field (Re{*E*_{z}}) over a plane 20 nm above the surface of graphene. The boundary conditions were set to perfectly matching layers. Graphene was modelled as a transition interface with a conductivity described by the local random-phase approximation model (see above)^{55,57}. We assumed a graphene carrier mobility of 2,000 cm^{2} V^{–1} s^{–1}. Supplementary Fig. 1c,d shows the permittivity of SiO_{2} and Au, respectively, at the mid-infrared wavelengths used.