Monoclinic crystals make up the largest crystal system, with around one-third of the minerals of Earth belonging to one of its three classes 18 . These low-symmetry Bravais lattices exhibit non-orthogonal principal crystal axes (Fig. 1a), in contrast to orthorhombic (for example, biaxial α-MoO 3 (ref. 1 )), tetragonal, hexagonal, trigonal (for example, uniaxial α-quartz, aQ, Fig. 1b) or cubic crystal systems. As a consequence, their dielectric permittivity tensor has major polarizability directions that strongly depend on the frequency, with off-diagonal terms that cannot be completely removed through coordinate rotation 6,7 , and exhibits shear terms analogous to viscous flow 19 . These features arise due to the non-trivial relative orientation (neither parallel nor orthogonal) of several optical transitions that, at a given frequency, contribute to a net polarization that cannot be aligned with the crystal axes. In turn, this property results in exotic light propagation not supported by higher-symmetry crystals 6,7,20 . Here we show exemplary consequences of these material features for nanophotonics, in particular, the emergence of a new form of waves -hyperbolic shear polaritons (HShPs)which have not been previously observed.
In this work, we theoretically and experimentally demonstrate the emergence of HShPs in monoclinic crystals. As an exemplary material to demonstrate this phenomenon, we study beta-phase Ga 2 O 3 (bGO), which has gained a large amount of research and industrial attention for its high breakdown field 21 and applications in photovoltaics 22 , optical displays 23 and gas sensors 24 . In the low-energy range, bGO features several strong infrared-active, non-orthogonal phonon resonances 6 , making the permittivity tensor of bGO naturally non-diagonalizable. Its low symmetry has two consequences on the polariton propagation when compared with more conventional hyperbolic materials with a diagonal permittivity tensor, such as hBN, aQ and α-MoO 3 . First, both the bGO polariton wavelength and the propagation direction strongly disperse with frequency. Second, we demonstrate that the asymmetric nature of optical loss in such crystals gives rise to shear, resulting in polariton propagation with tilted wavefronts. Such tilted wavefronts are a direct consequence of the low symmetry of the material and are one of the most notable and unique features of HShPs. New opportunities for polaritonics arise for HShPs stemming directly from their non-Hermitian and topological nature. Yet, surprisingly, they can be observed in low-loss, naturally occurring materials, without the need for artificial structuring of a material surface 17 .
To highlight the role of the asymmetry of monoclinic crystals in their polariton response, we compare HShPs with HPs supported by higher-symmetry anisotropic crystals, such as aQ 25 . In this vein, we compare the crystal structure of monoclinic bGO in Fig 1a with the trigonal crystal of aQ in Fig. 1b, illustrating the low crystal symmetry present in bGO. In general, the description of the dielectric response of monoclinic crystals requires inclusion of identical off-diagonal elements in the monoclinic plane within the frequency-dependent, complex-valued dielectric tensor.
The coordinate systems used to define the response of bGO and aQ are sketched in Fig. 1a, b, respectively. To analyse the properties of HShPs in monoclinic materials, we first rigorously solve Maxwell's equations (see Methods) to calculate the dispersion relation of the polaritonic modes supported by bGO and -for comparison -aQ, each at two distinct frequencies. Initially, we consider the lossless case, in which the imaginary part of each term in the dielectric tensor is neglected for both bGO and aQ. The solutions for the polariton wavevectors in both materials at two different frequencies are provided in Fig. 1c,d. For aQ, we observe two open hyperboloid surfaces -as expected for uniaxial hyperbolic materials -in which a change in frequency results in a corresponding change in wavevector, while preserving the hyperboloid orientation, that is, the direction of polariton propagation (Fig. 1d). By contrast, as we change the frequency, not only does the bGO polariton wavevector magnitude change but the direction of the hyperboloid also rotates within the monoclinic plane, as can be appreciated by examining the k z = 0 projections (Fig. 1c). This is a direct consequence of the non-trivial relative orientation of the phonon resonances supporting the hyperbolic response 6 , which results in polariton bands that disperse in azimuth angle as a function of frequency. This feature represents a signature of the reduced symmetry associated with HShPs supported in monoclinic crystals (and is also anticipated in triclinic crystals), in contrast to HPs observed in higher-symmetry lattices.
When we also account for natural material losses resulting from inherent phonon-scattering processes, the polariton propagation in bGO shows a reduced symmetry in comparison with hyperbolic polaritons in aQ, even at individual frequencies, as illustrated in Fig. 1e,f. In these panels, we show the results of full-wave calculations of dipole-launched surface polaritons propagating across the surface of a semi-infinite slab of bGO and y-cut aQ, in which -in both cases -natural material losses were explicitly taken into account. For in-plane hyperbolic materials, these surface waves show a hyperbolic dispersion within the surface plane and are referred to as hyperbolic surface or hyperbolic Dyakonov polaritons 26,27 , constituting a subset of HPs supported in these materials similar to volume-confined HPs in thin films. For aQ, HPs spread out along one crystal axis of the surface and are symmetric with respect to the crystal axes, as can be confirmed by a Fourier transform of the real-space electric field profile (Fig. 1h). However, for bGO (Fig. 1e), we observe that HPs are rotated with respect to the coordinate system of the monoclinic plane, as anticipated by the isofrequency contours (Fig. 1c). In addition, the wavefronts are tilted with respect to the major propagation direction, with no apparent mirror symmetry. This feature can also be clearly seen by examining the Fourier transform of the real-space profile (Fig. 1g), exhibiting a stronger intensity along one side of the hyperbola. These observations constitute the discovery of HShPs in low-symmetry crystals.
To experimentally demonstrate the effects of reduced symmetry in polariton propagation in bGO in contrast to higher-symmetry materials, we compare the azimuthal dispersion of HShPs in bGO to the one of HPs in aQ using an Otto-type prism-coupling geometry 28,29 (sketched in Fig. 2a; for details, see Methods). The experimental azimuthal dispersion of HPs on the surface of aQ is shown in Fig. 2b (see also Extended Data Fig. 3), in excellent agreement with the corresponding simulations (Fig. 2c). The dips in the reflectance spectra show the polariton resonances, which are only observable along specific azimuth angles and are symmetric about the crystal axes, ϕ = 0° (180°) and 90°. By contrast, the experimental azimuthal dispersion of HShPs on monoclinic bGO (Fig. 2d) exhibits no mirror symmetry, again in excellent agreement with the simulated dispersion (Fig. 2e).
To experimentally access the in-plane hyperbolic dispersion of the HShP in bGO observed in Fig. 1e, we mapped out the frequency-momentum dispersion in close spectral proximity of that mode (680-720 cm -1 ) at many azimuth angles (see Extended Data Fig. 4). The resulting map of polariton resonance frequencies is shown in Fig. 2f, in excellent agreement with the simulated resonance frequencies shown in Fig. 2g. These data allow extraction of single-frequency in-plane dispersion curves shown in Fig. 2h, i from experiment and simulations, respectively, for several selected frequencies, clearly demonstrating a hyperbolic dispersion, in excellent agreement with Fig. 1e. Notably, the base of the hyperbola shifts continuously with frequency, as marked by the symmetry axes for each curve in Fig. 2i, which directly leads to an asymmetric distribution of the group velocity along the hyperbolic dispersion curve, as shown in Fig. 2j (see also Extended Data Fig. 7).
The reduced symmetry observed in the polaritonic dispersion of bGO (Fig. 2d) is a direct consequence of the lack of symmetry in its vibrational structure 6 . Therefore, the HShPs are not propagating along fixed axes but show a continuous rotation of the HShP propagation direction as the frequency is varied. To describe the nature of this rotation, we diagonalize the real part of the permittivity tensor of bGO ε ω Re[ ( )] individually at each frequency, by rotating the monoclinic plane by the frequency-dependent angle.
The dispersion of γ(ω) is shown in Fig. 2e (white lines), illustrating that the major polarizability directions within the monoclinic plane, denoted as m and n, vary widely across the range. This frequency-dependent coordinate system enables an easier understanding and classification of the polaritonic response (see Methods and Extended Data Fig. 1 for details). The rotated coordinate axes are shown in Fig. 1e, g (see also Extended Data Fig. 2 for further modes), illustrating their alignment with the hyperbolic dispersion.
Although equation ( 2) describes the frequency variation of the polariton propagation direction, it does not capture the tilted wavefronts observed in Fig. 1e. This is because, as we choose the rotated coordinate system [mnz], we still retain a purely imaginary off-diagonal permittivity component (see Extended Data Fig. 1). These terms are associated with the non-orthogonal relative orientation of the material resonances, coupling the two crystal axes in the monoclinic plane. As a result, even in the rotated coordinate system [mnz], the dielectric tensor has off-diagonal terms associated with shear phenomena.
To selectively analyse the role of these shear terms, we simulate the polariton propagation in the rotated coordinate system [mnz] at 718 cm -1 . In particular, we include a scaling factor for the magnitude of the off-diagonal imaginary component, indicated as i × fIm(ε mn ), with f = 0, 0.5 and 1 (shown in Fig. 3a-c), while retaining the diagonal loss terms. When we remove the off-diagonal component (f = 0), bGO essentially becomes a shear-free biaxial material, akin to α-MoO 3 and similar to uniaxial aQ, with polaritons propagating along the optical axes (Fig. 3a). Therefore, whereas polariton propagation in such a fictional form of bGO is anisotropic in specific spectral ranges (similar to polaritons in MoO 3 (refs. 1,3 )), mode propagation without shear phenomena is symmetric about the (frequency-dependent) major polarizability axes (Fig. 3a). As we gradually increase the magnitude of the off-diagonal imaginary terms back to its natural value (f = 1), the wavefronts become increasingly skewed from the major polarizability axis (Fig. 3b,c). The respective reciprocal space maps in Fig. 3d-f show a strong symmetry breaking in the intensity distribution within the hyperbolic isofrequency curves. This observation provides further evidence that the propagation of polaritons is non-trivial within low-symmetry monoclinic -and, by extension, triclinic -systems, and it cannot be expected in higher-symmetry materials in which polariton propagation patterns are symmetric about the principal crystal axes 1,3,4 .
To connect the reduced symmetry of the surface subset of HShPs observed here experimentally (Fig. 2) and through our simulations (Fig. 3a-c) to the more general HShPs in the bulk, we now calculate isofrequency surfaces for polariton modes in bGO explicitly including loss, to account for the effect of shear phenomena. To this end, we solve Maxwell's equations for real momentum values, yielding complex frequency eigenvalues, whose imaginary part accounts for the finite lifetime of the supported modes (see Methods for details). The results of these calculations are shown in Fig. 3g-i. Here the real part of the eigenfrequency is fixed and we find its imaginary part ω i (Fig. 3h), which is proportional to the modal lifetime, and the corresponding value of k z (Fig. 3g) for each pair of k m and k n . The calculations are performed in the rotated coordinate system for both f = 0 and f = 1, showing that both the shape of the isofrequency surfaces as well as their lifetimes change greatly with the inclusion of the off-diagonal imaginary components. Notably, these calculations prove that, at individual frequencies and in the major polarizability frame, mirror symmetry of polariton propagation is lost in monoclinic materials as a direct consequence of shear.
To relate the isofrequency contours of HShPs to the surface mode dispersions in Fig. 3d-f, we plot the k z = 0 solution in Fig. 3i, with the colour scale indicating the relative loss ω i of the mode. Two important observations can be made: first, the mirror symmetry of the isofrequency curves is broken for f > 0 and it requires higher-order terms to account for the asymmetric shape; second, the mode losses are redistributed asymmetrically, with losses decreasing in one arm of the hyperbolae but increasing on the other arm. We note that, also in the experimental data (Extended Data Fig. 4), we observe an indication of asymmetric distribution of polariton quality factors along the hyperbolic dispersion curves (see Extended Data Fig. 5).
These observations naturally link HShPs in monoclinic crystals to the rich, emerging area of non-Hermitian and topological photonics. Although loss in orthogonal systems alone can already have interesting consequences for polariton propagation 30 , the off-diagonal shear terms highlighted here can provide new opportunities for non-Hermitian photonics and for manipulation of topological polaritons in low-symmetry materials. For instance, we foresee asymmetric topological transitions experienced by HShPs, generalizing previous results in orthorhombic systems 2 by exploiting the unique f, Experimental polariton resonance frequency map for bGO in the 680-720 cm -1 frequency range, extracted from Otto reflectance measurements at various incidence angles θ and azimuth angles ϕ. Experiments were performed at constant gap size d gap ≈ 4.0 μm. g, Simulated polariton resonance frequency map. Experimental (h) and simulated (i) in-plane hyperbolic dispersion for bGO at selected frequencies interpolated from f and g, respectively. Dashed lines mark extrapolated values outside the accessed momentum range. j, Radial component of the group velocity extracted from g. Dash-dotted lines in i and j mark the symmetry axis for each dispersion curve, which shift with frequency. The asymmetric distribution of the radial group velocity shows the asymmetry of energy flow for HShPs. non-Hermitian features emerging in low-symmetry materials. In addition, recent studies suggest the connection between Dyakonov surface waves and surface states emerging from one-dimensional band degeneracy (nodal lines) of topological nature of high-symmetry metacrystals 31 . We anticipate that HShPs may generalize these opportunities to asymmetric topological bands in which non-Hermiticity in the natural materials plays a dominant role.
Here we have demonstrated that low-symmetry crystals can support a new class of hyperbolic polariton modes with broken symmetry due to shear phenomena, which we refer to as HShPs. We introduce bGO as an exemplary material to enable the observation of these phenomena and experimentally demonstrate the symmetry-broken dispersion of the supported surface waves. The non-diagonalizable dielectric permittivity plays a key role in the unique properties of low-symmetry crystals, including monoclinic and triclinic lattices. Our findings are generalizable to engineered photonic systems with at least two non-orthogonal oscillators, including new metasurface designs capturing these physics. Beyond the results provided here for intrinsic, compensation-doped bGO, the presence of free charge carriers in bGO 32 may allow for methods for direct steering of the HShP propagation direction ͠ calculated in the rotated frame for f = 0 and 1 (red and blue, respectively). See Methods for details on the approach. h, Imaginary part ω i for f = 0 and 1 (red and blue, respectively). i, Contour lines of the isofrequency surface at k z = 0 for f = 0, 0.5 and 1. The imaginary part ω i at the corresponding point in k-space is colour-coded. FFT, fast Fourier transform.
(see Extended Data Fig. 6). Finally, exfoliation of thin flakes of single-crystal bGO has also been recently reported 33 , which will allow to make use of volume-confined HShPs in such bGO thin films or -potentially -even in monolayers 34 . We anticipate that HShPs may have important implications in the manipulation of phase and directional energy transfer, including radiative heat transport 35 , ultra-fast asymmetric thermal dissipation in the near field 35 and gate-tunability for on-chip all-optical circuitry 36 . Beyond advances in nanophotonics, infrared polariton propagation has been demonstrated as a means for quantifying crystal strain 37 , polytypes 38 , variations in free-carrier density, as well as phononic and electronic properties around defects 39 , thereby also promising a new metrology tool for characterizing low-symmetry ultra-wide-bandgap semiconductors. We highlight that our results are applicable to any material with non-orthogonal optically active transitions and may therefore be extended to other optical phenomena, such as excitons in triclinic ReSe 2 (ref. 40 ).
The insulating (010)-oriented, 5 × 5 × 0.5-mm 3 bGO substrate was produced by means of Fe compensation doping and was purchased from Novel Crystal Technology, Inc., Japan. The aQ sample was purchased from MaTeck GmbH, Germany. The absolute azimuth orientation of the samples was extracted from a global fit for each of the datasets of aQ and bGO (plotted in Fig. 2b,d, respectively), resulting in a rotation with respect to the principal x axis of the laboratory coordinate system of Δϕ bGO = 27.95° and Δϕ aQ = 26.96°. The aQ data have been rotated accordingly to shift the crystal axes (mirror planes) onto multiples of 90°. On the other hand, the bGO data are plotted as measured, as there is no principal azimuth angle for alignment because of the broken mirror symmetry. Here the simulation was rotated accordingly to match the data.
The Otto-type prism-coupling experiment measures the spectral dependence of surface waves through sharp absorption peaks observed as dips in the reflectance spectra, by using a prism placed near the material surface 28,41 . The crystals are oriented such that the monoclinic plane (bGO) and the optical axis (aQ) are parallel to the sample surface. By following the spectral position of the polariton resonances as a function of azimuth angle, we investigate the dispersion of hyperbolic waves at the surface for both bGO and aQ. The Otto geometry effectively selects a specific in-plane momentum component of those surface waves induced by the dipole excitation in Fig. 1e, f, as set by the incidence angle θ and azimuth angle ϕ that define the magnitude 28,29 and direction of the selected momentum, respectively.
As an excitation source for the Otto-type prism-coupled experiments, we use a mildly focused mid-infrared free-electron laser (FEL) (spot size ~1 mm 2 ) with small bandwidth (~0.3%) and wide tunability of 3-50 μm, covering the spectral range 350-800 cm -1 , in which aQ and bGO support polaritonic modes (details on the FEL have been reported elsewhere 42 ). Although the frequency is scanned by tuning the FEL, different in-plane momenta can be accessed by means of changes in the incidence angle θ by rotating the entire Otto geometry (details on the setup have been reported elsewhere 28,43 ). For the experiments shown in Fig. 2b,d, the incident angle was fixed to 28°, resulting in an in-plane momentum of k k / ≈1.10 0 (at 500 cm -1 ). In contrast to alternative approaches, the Otto geometry features experimental control over the excitation efficiency through tunability of the air gap width d gap . Here the gap was adjusted to a separation in which all excited modes could be observed in the spectra simultaneously, that is, d gap ≈ 8.3 μm for bGO and d gap ≈ 14.4 μm for aQ. Direct readout of d gap with a range of 1-50 μm is realized through white-light interferometry, whereas the contrast of the interference range grants parallel alignment of prism and sample 43 .
Mapping of the in-plane hyperbolic dispersion (Fig. 2f,h) was performed analogously to the Otto reflectance measurements shown in Fig. 2b,d. However, here we additionally varied the incidence angle θ to map out the frequency-momentum dispersion at each azimuth angle. Reflectance spectra were taken in a narrow frequency range of 670-730 cm -1 , at θ = 26°, 28°, 30°, 32° and 34°, corresponding to in-plane momenta of k k / ≈1.03, 1.11, 1.18, 1.25 and 1.32 0 (at 700 cm -1 ), at nine azimuth angles. To allow prism coupling to the polaritons also for larger momenta, these data were taken at a constant air gap d gap ≈ 4.0 μm. The reflectance minima marking the polariton resonance were extracted from these data and are shown in Fig. 2f. The theoretical polariton resonance map, Fig. 2g, was calculated using a transfer matrix formalism 44 , by extracting the peak positions of Im(r pp ) of the air-bGO interface. To extract the single-frequency in-plane hyperbolic dispersion curves (Fig. 2h,i), we interpolated the momentum for a given frequency in the frequency-momentum dispersion for each measured/calculated azimuth angle.
Transfer matrix. The calculations of the optical response shown in Fig. 2c, e and the polariton resonance map in Fig. 2g, as well as the dispersion maps in Extended Data Fig. 2c,d, were performed using a generalized 4 × 4 transfer matrix formalism 44 . In short, the formalism enables the calculation of reflection and transmission coefficients in any number of stratified media with arbitrary dielectric tensor, which enables to account for the anisotropy of our samples.
COMSOL simulations. COMSOL 45 version 5.6 was used for simulating point dipole excitation of HShPs on bGO. A point dipole was placed 100 nm above the surface of an infinite slab of bGO, with a dielectric permittivity matching that of ref. 6 . The dielectric function of aQ was taken from ref. 46 . Perfectly matched impedance boundary conditions were used on the sides of the simulation, which -in principle -absorb all radiation. However, to account for the imperfect behaviour of the boundaries, we ensured that the bGO slab was sufficiently large (250 × 250 × 8 μm), such that the wave is sufficiently damped when it reaches the boundary so as not to influence the results. Reciprocal space maps were generated by 2D Fourier transformation of the real-space electric field profiles.
To obtain the isofrequency contour surface of the bulk wave when losses are considered in the materials, we turn to the complex-frequency method and solve the source-free Maxwell equations as follows:
∼ in which ω i is an effective inverse mode lifetime, calculated neglecting the effect of the complex frequency on the material dispersion, and normalize the wavevector as k =
x y z k k , ,
x y z , , 0 ∼ . Note that we choose a negative sign for the timedependent term e -iωt , so ∼ ω i must be real and negative to reflect the decay- ing nature of the wave. The analytic expression ( )
is found by means of the secular equation of the above matrix, and two equations are obtained by separating the real and imaginary components of F, namely, ( ) ( )
x y z x y z
The isofrequency contour of the bulk wave and the imaginary component ω i ∼ are evaluated from those two equations. The numerical examples at 718 cm -1 are given in Fig. 3g, h on the basis of the above method (k x,y → k m,n ). Notice that, in this procedure, we use the permittivity tensor elements calculated at Re(ω). When k = 0 z ∼ , the analytic expression for the isofrequency contour of the bulk wave can be written as
which turns into two equations by separating the real and imaginary parts, and they can be further simplified into the following equation Therefore, the isofrequency curves in the k = 0 z ∼ plane as well as the imaginary component of the bulk complex frequency are obtained from the above expression.
Assignment of the polariton mode nature. First, we briefly outline how conventional polariton materials are classified. A material supports surface polaritons at frequencies for which the real part of the crystal permittivity fulfils Re(ε) < -1 (ref. 47 ). In uniaxial crystals, the diagonal permittivity elements can be of different sign, leading to hyperbolic behaviour in which either the real part of one element is negative and the other two are positive (type I) or two are negative and one is positive (type II) 9 . However, this classification of anisotropic materials relies on the off-diagonal permittivity tensor elements being zero at all frequencies for an appropriate choice of coordinate system. The lower symmetry of bGO requires the emergence of ε xy ≠ 0. The permittivity elements using coordinates as indicated in Fig. 1a are shown in Extended Data Fig. 1a-e, but no coordinate system exists in which Re(ε x,y ) ≠ 0 at all frequencies. This, as we have demonstrated in the Otto-geometry experiments (Fig. 2d), results in a non-trivial polaritonic response with highly directional modes that propagate along frequency-dictated propagation angles in the a-c plane. Furthermore, because ε xy ≠ 0, it is not straightforward to determine whether the modes are elliptical (ε xx , ε yy , ε zz < -1) or hyperbolic in nature (type I or type II), as the propagation angle is typically not aligned with one of the principal axes.
To unambiguously describe the nature of the polariton modes, we switch to a frequency-dispersive coordinate system [mnz], in which the real part of the permittivity tensor is diagonal. This is achieved by rotating the monoclinic plane by the frequency-dependent angle γ(ω) (equation ( 2)). The dispersion of γ(ω) and the resulting diagonal elements of ε [mnz] are plotted in Extended Data Fig. 1f-j. This new frequency-dispersive coordinate system enables the unique assignment of the supported polariton mode nature, which we have colour-coded in Extended Data Fig. 1g-j. For bGO, we observe the full range of possible combinations of positive and negative real parts of ε mm , ε nn and ε zz , leading to dielectric (white), elliptical (grey) and hyperbolic spectral regimes of type I (in-plane in blue, out-of-plane in red) and type II (in-plane in green, out-of-plane in yellow). By performing such a frequency-dependent rotation of the permittivity tensor, we have simplified the system into a pseudo-biaxial crystal at each frequency. However, as the dielectric tensor of a monoclinic crystal is not diagonalizable 20,48,49 , the in-plane, off-diagonal element of ε [mnz] retains a non-vanishing imaginary part at all frequencies, that is, Im(ε mn ) ≠ 0 (plotted in Extended Data Fig. 1i), giving rise to the reduced symmetry of hyperbolic shear polaritons in monoclinic crystals, as discussed in Fig. 3.
The frequency-dependent rotation of the dielectric permittivity tensor is performed in three subsequent steps. First, the in-plane permittivity tensor as shown in Extended Data Fig. 1a-c is rotated about the angle γ(ω) (equation ( 2)). However, γ(ω) has jumps of 90° at arbitrary frequencies, resulting in abrupt discontinuities in the real parts of ε mm and ε nn , for which the two curves switch values. By analysing the derivative of ε mm , we extract the frequency values ω jump , in which the jumps occur and reassign the permittivity curves, respectively. At the eight in-plane TO frequencies ω TO of bGO, the permittivity curves feature a pole, which is also captured in the analysis of the derivative. Therefore, at this step, the resulting curves are smooth between the TO frequencies, but switch assignment at every ω TO . The switching of the curves at every ω TO is performed in the last step. However, near ω TO , the permittivity features a large imaginary part, which is not accounted for in the rotation angle γ(ω). This leads to poles in γ(ω) at ω TO (see Extended Data Fig. 1f), which -in turn -results in a small avoided crossing of ε mm and ε nn at the TO frequencies. Therefore, the reassignment of the last step results in discontinuous solutions near ω TO , which is clearly not physical. To resolve this issue, we cut out ±1 cm -1 in both ε mm and ε nn at all eight ω TO and smooth the curves by interpolation, resulting in the pseudo-biaxial permittivity curves as shown in Extended Data Fig. 1g-h. The rotation about γ(ω) also leads to discontinuities in the off-diagonal imaginary part Im(ε mn ) at the frequencies ω jump . However, because ε mn = ε nm , the abrupt rotation about 90° only leads to sign changes at every ω jump . The curve of Im(ε mn ) shown in Extended Data Fig. 1i is corrected for these sign changes.
Polariton behaviour in bGO in the rotated frame. To verify the polaritonic behaviour of bGO in the rotated frame, that is, the pseudo-biaxial crystal, we subsequently analyse the surface polariton dispersion in the rotated coordinate system [mnz] in Extended Data Fig. 2. For electric fields in the m-z or n-z planes, the analytical expression describing extraordinary surface polaritons in uniaxial crystals can be used 50 in which ε in-plane = ε mm , ε nn . In bGO, the solutions yield four polariton branches for each direction, m and n, respectively, plotted as red dotted lines in Extended Data Fig. 2c,d. These analytical results are in perfect agreement with the numerically obtained surface polariton dispersion using a transfer matrix formalism 44 . To obtain the polariton propagation properties of the system, we calculate the full electric field patterns by placing a point dipole source above the bGO surface at x = y = 0 and simulating the optical response along the bGO-air interface (z = 0) with COMSOL Multiphysics 45 . The real-space field profiles clearly show the rotation of the major polarizability direction as a function of frequency, demonstrated for six different modes M1-M3 and N1-N3 (See Extended Data Fig. 2e-g, k-m). Frequencies are indicated as black dash-dotted lines in Extended Data Fig. 2a-d. Mode N4 is shown in Fig. 1e,g. The field profiles align with the rotated coordinate system, with basis vectors indicated by the 'm' and 'n' crosshair in each figure.
To relate the calculated dispersion of the polariton branches to the field profiles, we calculate the momentum-k maps of these modes, as obtained by a 2D Fourier transformation of the respective electric field patterns of Extended Data Fig. 2e-g, k-m in Extended Data Fig. 2h-j, n-p, respectively. At all selected frequency positions, the electric field patterns contain a directional wave of large amplitude with low spatial frequency, as well as a wave with high spatial frequency. The observed in-plane momenta of the low-k modes follow the modal dispersion predicted in Extended Data Fig. 2c,d, along the m and n axes for modes M1-M3 and N1-N3, respectively, as indicated by the black circles in Extended Data Fig. 2h-j, n-p. According to the mode characterization provided in Extended Data Fig. 1, these modes are hyperbolic, either of type I in-plane (M1-M3 and N1) or of type II in-plane (N2, N3). For all HShP modes, field patterns and k-space maps are characterized by twofold rotational symmetry only, in agreement with the 2D plane group 2 (no mirror plane symmetries). Further, the distinct peaks in the k-space maps verify the principal polariton propagation direction, whereas the corresponding radial and azimuthal spreads are representative of their decay length and degree of directionality, respectively. Analogous to the model at 718 cm -1 discussed in the main text, the maxima of the k-space maps shown in Extended Data Fig. 2 do not lie on the major polarizability axes (most prominently for the cases in Extended Data Fig. 2n,o), owing to shear phenomena in monoclinic bGO. This is discussed further in the main text and in Fig. 3. 44 . The four supported polaritons along each axis m and n are clearly distinguishable, in perfect agreement with the theoretically calculated polariton dispersion (red dotted lines) 50 . Black horizontal dash-dotted lines mark the frequencies M1-M3 and N1-N4 at which the electric field distribution is plotted. N4 is shown in Fig. 1e,g. e-g, Real-space electric fields at the bGO surface at frequencies M1-M3, respectively. h-j, The respective two-dimensional Fourier transformation. The fields were calculated using COMSOL Multiphysics 45 (see Methods for details). k-m, Real-space electric fields at frequencies N1-N3, respectively. n-p, The respective Fourier transforms. All maps (e-p) were calculated using the non-dispersive permittivity tensor (Extended Data Fig. 1a-e), thus showing rotated field patterns with different orientations, depending on the frequency. The thin black and white crosshairs indicate the principal axes of the respective frequency-dispersive coordinate system, its rotation given by γ at the corresponding frequency. Small black circles in h-j and n-p mark the momentum value of the analytical dispersion in c and d, respectively. Extended Data Fig. 3 | Experimental datasets for bGO and aQ at different gap sizes. The gap size d gap in our Otto geometry setup can be tuned and monitored 43 , enabling control over the excitation efficiency of the polariton modes 44 . Datasets measured for bGO (a-c) and for aQ (d-f) at three different gap sizes each. For smaller gaps, some modes are overcoupled and their resonance features broadened (such as the mode at 500 cm -1 in bGO), whereas for larger gaps, some modes are undercoupled and their resonance features too weak to be clearly distinguishable (in particular, the mode at 725 cm -1 in bGO). The centre gap sizes compromise between these effects. Note that the gap sizes indicated here are the values monitored with a white-light interferometry setup 43 . The fits performed for the datasets shown in Fig. 2, that is, the datasets shown in Extended Data Fig. 3b for bGO and Extended Data Fig. 3e for aQ, however, yielded larger gap sizes of 8.3 μm (bGO) and 10.4 μm (aQ). The offset can be attributed to non-perfect parallel alignment between prism and sample and a lateral offset between the polariton excitation site with the FEL and the white-light spot for the gap measurement.