Twist-tunable polaritonic nanoresonators in a van der Waals crystal

Optical nanoresonators are key building blocks in various nanotechnological applications (e.g., spectroscopy) due to their ability to effectively confine light at the nanoscale. Recently, nanoresonators based on phonon polaritons (PhPs)—light coupled to lattice vibrations—in polar crystals (e.g., SiC, or h-BN) have attracted much attention due to their strong field confinement, high quality factors, and their potential to enhance the photonic density of states at mid-infrared (mid-IR) frequencies, where numerous molecular vibrations reside. Here, we introduce a new class of mid-IR nanoresonators that not only exhibit the extraordinary properties previously reported, but also incorporate a new degree of freedom: twist tuning, i.e., the possibility of controlling their spectral response by simply rotating the constituent material. To achieve this result, we place a pristine slab of the van der Waals (vdW) α-MoO3 crystal, which supports in-plane hyperbolic PhPs, on an array of metallic ribbons. This sample design based on electromagnetic engineering, not only allows the definition of α-MoO3 nanoresonators with low losses (quality factors, Q, up to 200), but also enables a broad spectral tuning of the polaritonic resonances (up to 32 cm−1, i.e., up to ~6 times their full width at half maximum, FWHM ~5 cm−1) by a simple in-plane rotation of the same slab (from 0 to 45°). These results open the door to the development of tunable and low-loss IR nanotechnologies, fundamental requirements for their implementation in molecular sensing, emission or photodetection applications.

to 200) can be tuned in a broad range (up to 32 cm -1 , i.e up to ~ 6 times its full width at half maximum, FWHM ~ 5 cm -1 ).Our results open the door to the development of tunable low-loss nanotechnologies at IR frequencies with application in sensing, emission or photodetection.KEYWORDS: nanoresonators, van der Waals crystals, phonon polaritons, twisted heterostructures The excitation of PhPs in vdW crystals have recently emerged as an attractive strategy for manipulating IR light on deeply subwavelength scales [1][2][3] , enabling fingerprint identification of organic molecules as well as strong coupling phenomena 4,5 .Some of the PhPs, such as those excited in -MoO3 [6][7][8][9][10][11] or -V2O5 12 , present long lifetimes (up to a few ps) and in-plane hyperbolic propagation, being thus very appealing for the development of optical nanoresonators with ultralow losses and potentially other unique properties.However, structuring of vdW crystals -needed for engineering nanoresonators -remains a challenging task, since the fabrication process dramatically degrades their optical properties, and thus the lifetime of PhPs 13 .Moreover, in contrast to plasmonic resonances in 2D materials, which are actively tunable via an external gate 14- 16 , tunability of PhPs nanoresonators has remained elusive.Here, we demonstrate PhPs nanoresonators fabricated in the van der Waals crystal -MoO3 that not only exhibit low losses, but also incorporate a unique twist tuning, i.e. the capability to be spectrally tuned by a simple rotation.The fabrication of these nanoresonators is based on placing a continuous (pristine) thin biaxial vdW crystal slab of -MoO3 on top of a grating formed by metal nanoribbons (see Methods).Remarkably, such design avoids any degradation of the optical properties of the crystal due to fabrication, maintaining the intrinsic low-loss and in-plane hyperbolic propagation of PhPs in -MoO3.The refractive indexes of the PhP modes supported by the -MoO3 slab above the air and metal regions are different, so that propagating PhPs modes can bounce back and forth between the boundaries of these regions (i.e., -MoO3/air and -MoO3/metal), eventually forming resonances.Most importantly, since the in-plane hyperbolicity of -MoO3, implies that different in-plane directions support PhPs with different momenta, the mutual orientation of the crystal axes and the metal grating can serve as a tuning knob, in analogy to recent demonstrations of twistoptics in -MoO3 bilayers [17][18][19][20] .Thus, by twisting the -MoO3 slab with an in-plane angle with respect to the metal ribbons, we demonstrate a broad tunability of the PhPs nanoresonators, as revealed by both far-and near-field techniques: Fourier-transform IR spectroscopy (FTIR) and scattering-type scanning near-field optical microscopy (s-SNOM), respectively.Furthermore, with the help of theoretical analysis, we are able to interpret and disentangle the different resonant PhPs modes observed in the experimental data.Figure 1a illustrates schematically the heterostructure designed to define "dielectric-engineered" PhPs nanoresonators in our work.It consists of an -MoO3 slab placed on top of a periodic array of metal ribbons fabricated on a IR-transparent CaF2 substrate (see the AFM profile in the inset to Figure 1a and the optical image in Figure 1b), similarly to the concept of polaritonic nanoresonators based on in-plane isotropic slabs of h-BN and graphene layers [21][22][23] .The -MoO3 slab is twisted with respect to the direction across the ribbons (x-direction) by an arbitrary angle .The metal arrays (with ribbon widths  = 1480 nm and a separation distance  = 1230 nm) have been designed so that the heterostructure exhibits PhPs resonances at mid-IR frequencies (in the range  = 850 − 1010 cm -1 ).The reflection spectra taken by FTIR upon incident illumination polarized across the ribbons and along the -MoO3 [100] crystal direction ( = 0º) is illustrated in Figure 1d (black curve).To better recognize the resonance spectral features, we normalize the reflection coefficient, , to its moving average,  , obtaining the relative reflection  = ( −  )/ (Figure 1d, red curve), which clearly demonstrates a set of sharp resonance peaks.Particularly, these peaks appear within the RBs of the -MoO3 crystal (Figure 1e), defined in the spectral ranges  = 850 − 956.8 cm -1 , where Re( ) < 0, and  = 963 − 1006 cm -1 , where Re( ) < 0, exhibiting Q as large as 180 at 916 cm -1 and 200 at 993 cm -1 (see Supporting Note V).In the lower (upper) RB PhPs exhibit hyperbolic (elliptic) dispersion 24 , composing a set of in-plane anisotropic electromagnetic modes propagating along the crystal slab.Importantly, the PhPs modes in the region -MoO3/air (Mla) are expected to have a different refractive index (momentum) to those in the region -MoO3/metal (Mlm), even for the same mode number, , the latter representing the quantization of the PhP's electric field in the z-direction 9,25 .This is corroborated by plotting the dispersions of the fundamental modes in both regions within the hyperbolic and elliptic RBs (M0a, M1a and M1m in Figure 1f, with the hyperbolic modes corresponding to the -MoO3 [100] crystal direction).Clearly, the M1m mode exhibits a much shorter wavelength (larger momentum), together with a stronger field confinement to the faces of the -MoO3 slab (Figure 1c).We thus assume that the resonant features seen in the spectra of  To corroborate our assumption on the physical origin of the PhPs resonances and explain the reflection spectra obtained experimentally, we perform a theoretical analysis based on full-wave simulations of a structure mimicking the fabricated nanoresonators.The simulated relative reflection spectrum (grey curve in Figure 2a) shows a good agreement with the experiment (red curve in Figure 2a) for the case of  = 1480 nm and  = 1230 nm.We, therefore, extend the theoretical study for a broad range of  and  in Figure 2b, where , represented by a colorplot, shows a set of maxima that clearly depend on the ribbon width or the distance between the ribbons.These resonances presumably coincide with an anisotropic FPRs condition, which can be written in a simple form according to the PhP mode phase matching as  •  +  = , analogously to that for in-plane isotropic FPRs in graphene ribbons 23,26 .In our case,  coincides with either  or , depending on the region in the -MoO3 slab where the nanoresonators are defined (above the metal or air, respectively),  and  are the momentum and reflection phase of the Mth mode, respectively, and  is the integer numerating the FPRs (quantization of the Mth PhP mode in the x-direction), indicating how many wavelengths of the mode fit across either the distance  or .This numeration includes only the "bright" FPRs, corresponding to the non-zero overlap between the illuminating wave and the PhPs in a FPR cavity 26 .Plugging the  from the dispersion relation of PhP modes in a biaxial crystal slab 9,25 into the FPR condition, we obtain a set of curves for each of the modes (rendered on top of the colorplot in Figure 1b).These curves match the maxima of , thus confirming the Fabry-Pérot origin of the resonances in our polaritonic structures.The FPR can also be recognized in the calculated spatial distribution of the vertical electric field above the unit cell of the periodic structure (see schematics in Figure 2f), represented as a function of frequency  and coordinate  in Figure 2c.Indeed, at the peak frequencies the number of field oscillations across the ribbon matches with the quantization number  in the FPR condition, both for the PhP modes in the -MoO3/air regions (marked with "na") and in the -MoO3/metal regions (marked with "nm").
Importantly, the spectral position of the PhP resonances strongly depends on the twist angle  between the -MoO3 slab and the metallic ribbons, as clearly observed in the reflection spectra of Figure 2d.In particular, the main resonant peak (1a) shifts to lower frequencies (from ~ 909.5 to ~ 877 cm -1 ) with increasing  (the increasing signal-to-noise ratio enables determining the resonant features for  up to 45º), which can also be seen in our theoretical calculations (Figure 2e).Note the much higher Q (see Supporting Note V) of our twist-tunable PhP resonators (by more than an order of magnitude, due to both the absence of the electron scattering and avoiding the direct crystal structuring) compared to their plasmonic counterpart reported at THz frequencies in in-plane anisotropic WTe2 crystals 27 .To further corroborate our results and the structure of the electric field in the anisotropic PhP nanoresonators, we perform near-field measurements by s-SNOM (see Methods).By scanning the -MoO3 slab with the s-SNOM tip (Figure 3a), we record the electromagnetic signal as a function of the tip position, composing the near-field images (Figures 3c-e).At a frequency  = 909 cm -1 (Figure 3c), corresponding to the 1a and 6m resonances in the far-field spectra (Figure 2a), we observe signal fringes across the ribbons, indicating the excitation of PhPs.The much larger distance between the red and blue fringes in the region -MoO3/air (1 oscillation) compared to the region -MoO3/metal (6 oscillations) is consistent with the much longer PhP wavelength (see the simulated PhP field distribution along the black and grey vertical dashed lines, "1a, 6m", in Figure 2c), as expected from their dispersion relation (Figure 1f).For a better visualization, we represent the simulated Re( ) as a function of the in-plane coordinates ,  in Figure 3b, which perfectly matches with the near-field image in Figure 3c.From these observations, we can thus conclude that: (i) the s-SNOM tip-scattered signal represents a measure of Re( ) above the slab 13,28 and (ii) the visualized near-fields at the resonant peaks of the far-field spectra are clearly consistent with the FPR model.Remarkably, changing both the measuring frequency and the twist angle, following the experimental data shown in Figure 2d, the number of oscillations in the near-field images keeps constant (the agreement between the periods of the near-field signal oscillations across the ribbons is shown in Figures 3c-e by the vertical dashed eye-guides).Namely, the calculated PhPs wavelength in the area above the air gap (1a resonance) is equal to  = 1.29 μm, approximately fitting the air gap width (1.23 μm), while in the area above the metal (6m resonance)  = 249 nm, thus being consistent with the ribbon width divided by 6 (247 nm).Therefore, our near-field measurements reconfirm the emergence of FPR of the same order for the chosen parameters.On the other hand, although similar oscillations are observed in the near-field images, the FPR originates from PhPs at different points in the dispersion surface (oscillating at different frequencies and propagating along different directions).In order to illustrate "probing" of the PhP modes in our nanoresonators at different points of the dispersion surface, we mark the positions of the main peaks in the normalized reflection spectra (Figure 2d) directly on the hyperboloids representing the dispersion of the M0a PhP mode (Figure 4a).The experimentally-measured positions of the resonances (green, black, blue and red points in Figure 4a) can be represented in polar coordinates by the frequency, , of the resonant peak (z coordinate), the inverse air gap width, 2/ (radial coordinate), and the twist angle, .The selection of the PhPs momenta for the first-order FPR at any  can be visualized by a cylinder with radius, 2/, passing through the dispersion surface (Figure 4a).Crossings between the latter and the cylinder give the approximate positions of the FPR, which can be more clearly seen in the  ,  plane (Figure 4b), where the isofrequency curves (IFCs) for each resonant frequency are shown.The experimental points fall very close to the crossings between the circle and the IFCs, thus reconfirming our initial assumption on the Fabry-Pérot origin of the resonances.
Alternatively, the PhP dispersion surface can also be probed by cutting it with planes of constant frequencies, .To that end, we performed an additional set of far-field measurements for metal gratings with different air gaps, , and twist angles, , all of them matching the same FPR (see schematics in Figure 4).Namely, we designed our structures to exhibit the "1a" FPR at two fixed frequencies: 908.5 cm -1 and 911 cm -1 (see measurements in the Supporting Note I).The positions of the measured resonant peaks perfectly match with the crossing points between the calculated IFCs and the radial lines from the origin at an angle  (dashed lines in Figure 4d), as seen in Figures 4c,d for both frequencies.Note that a similar analysis can be also easily performed for the FPRs of other orders, and particularly for the FPRs composed by the PhP modes above the metallic region, M1m (see Supporting Note III).This analysis is equally suitable for the upper RB, i.e. for the elliptic regime (see Supporting Note IV).
To summarize, we introduce Fabry-Pérot nanoresonators that exhibit a unique tuning by a simple rotation of the host crystal.This is achieved by fabricating heterostructures composed of metal gratings and twisted vdW crystal slabs supporting in-plane anisotropic PhPs (-MoO3).In contrast to conventional Fabry-Pérot polaritonic resonators, in which the reflecting boundaries are typically fabricated by etching the polaritonic material, the design of our tunable Fabry-Pérot nanoresonators allows preserving the crystalline properties of the slabs and thus obtaining high Q.
Interestingly, the resulting FPRs, visualized in real space by s-SNOM, allow reconstructing the three-dimensional anisotropic dispersion surfaces of the PhPs from data collected at either fixed momentum or fixed frequency.Such reconstruction of the polaritons dispersion opens the door to an alternative characterization of novel vdW materials, particularly for a more reliable determination of their optical properties.From an application point of view, our nanoresonators can be very appealing for tunable mid-IR narrow-resonance sensors, and in particular for molecular bar-coding 29 .

METHODS
Fabrication of metal ribbons.Arrays of ribbons of both Au and Al were fabricated by electron beam lithography.In the case of Al ribbons, first, a 70 nm-thick Al layer was deposited by electron beam evaporation on a CaF2 substrate.Subsequently, a negative photoresist layer (MA-N2401) was spin coated (3000rpm -> 90 nm) on top of the Al layer for the definition of the ribbon pattern by electron beam lithography (50keV, 200pA, dose 280 uC/cm 2 ).A high-resolution developer AZ726 was used for the lift-off.By reactive-ion etching of Al in BCL3/Cl2 plasma (pressure 40mT, RIE power 100W) and removing the photoresist in oxygen plasma, we got the set of 50 μm (length) × 60 μm (width) × 70 nm (height) gratings with different ribbon widths.In the case of Au ribbons, high-resolution electron beam lithography was employed under 100kV and 100pA.First, the samples were coated with a positive resist layer (PMMA).To dissolve the exposed areas, a conventional high-resolution developer (1:3 MIBK: IPA) was used.Afterwards, 5 nm of Cr and 50 nm of Au were evaporated.Finally, the lift-off was performed to define the 50 µm (length) x 50 µm (width) x 50 nm (height) Au grating with ribbons of 1.48 µm width.We used the set of Al gratings exclusively for measurement presented in the Figures 4c,d.For all other measurements we used the Au grating.
Fabrication and In-plane Twist of -MoO3 Flakes.-MoO3 flakes were mechanically exfoliated using a Nitto tape (Nitto Denko Co., SPV 224P) from commercial -MoO3 bulk crystals (Alfa Aesar).A second exfoliation was performed from the tape to transparent polydimethylsiloxane (PDMS) in order to thin them down.The flakes were examined with an optical microscope in order to select homogeneous pieces with the desired thicknesses (110 nm and 127.5 nm) and large surface areas, which were picked up with the help of polycarbonate (PC).
The dry transfer technique was used, with the help of a micromanipulator, the flakes were precisely aligned and twisted on top of the metal grating.Transferring was carried out by heating up to 250°C in order to liquate the PC.The PC was removed with chloroform at 100°C releasing the flake.
Fourier-Transform Infrared Spectroscopy.The far-field optical response of our -MoO3 nanoresonators was characterized by Fourier-Transform Infrared Spectroscopy (FTIR) using a Varian 620-IR microscope coupled to a Varian 670-IR spectrometer, supplied with a broadband mercury cadmium telluride (MCT) detector (400-6000 cm -1 ).The spectra were collected with a 2 cm -1 spectral resolution and the spatial resolution was adjusted to the size of the metal arrays.The infrared radiation from the thermal source (normal incidence) was linearly polarized employing a wire grid polarizer.CaF2 substrates were used due to their transparency in the spectral range under study.
Scattering-Scanning Near Field Optical Microscopy.Near-field imaging measurements were performed by employing a commercial scattering-type Scanning Near Field Optical Microscope (s-SNOM) from Neaspec GmbH, equipped with a quantum cascade laser from Daylight Solutions (890-1140 cm -1 ).Metal-coated (Pt/Ir) atomic force microscopy (AFM) tips (ARROW-NCPt-50, Nanoworld) at a tapping frequency  ~ 280 kHz and an oscillation amplitude ~ 100 nm were used as source and probe of polaritonic excitations.The light scattered by the tip was focused by a parabolic mirror into an infrared detector (Kolmar Technologies).Demodulation of the detected signals to the 3rd harmonic of the tip frequency (s3) was carried out for background suppression.A pseudo-heterodyne interferometric method was employed to independently extract both amplitude and phase signals.
Full-wave Numerical Simulations.Two types of full-wave numerical simulations were performed using COMSOL software, based on the finite-element method in the frequency domain.
In both cases, the structure was composed of the -MoO3 flake on top of the array of metal ribbons placed on the semi-infinite CaF2 substrate.The first type of simulation was based on the far-field illumination of the structure by a normally-incident plane wave polarized across the ribbons.The ratio of the ribbon to the air gap widths was fixed constant / = 1480/1230 (the experimental measurements for the array of ribbons with such parameters were performed in Figures 1d, 2a and   3b-f).The reflection coefficient and the spatial distribution of the vertical electric field,  , above the -MoO3 slab have been extracted from the simulations (see Fig 2a-e, Fig 3b-e).In the second type of simulation, we have used the quasi-normal eigenmode analysis in order to find the electric field distribution of the modes M0 and M1 (see Fig 1c).
In the analytical analysis of the FPR (Figure 2b and Figure 4), the reflection phase of the PhPs modes,  , in the -MoO3/metal and -MoO3/air regions is taken 0 and π, respectively.

Supporting Information.
The Supporting Information is available free of charge in Supporting_information.pdf .
In this Supporting Information we will consider the Fabry-Pérot resonances (FBRs) with their origin in: i) the M1 PhP mode propagating in the -MoO3 flake above the metal region in the hyperbolic frequency range, and ii) the M1 PhP mode propagating in the flake above both the metal ribbon and the air gap regions in the elliptic frequency range.

I.
Far-field Characterization  Additionally, we fabricated Al gratings (as the type of metal does not affect our results in the mid-IR, in some cases we use Al instead of Au for convenience) with different periods (FigureS1b).In this case, the  flake placed on top had a thickness t = 127.5 nm.By varying the twisting angle of the flake with respect to the ribbons axes, we observe the same resonant frequencies as when using Au in both hyperbolic and elliptic frequency ranges (FigureS2a,b, and FigureS2c, respectively).The FPRs modes 1a and 6m in the hyperbolic range appear at = 908.5 cm -1 (Figure S2a) and= 911 cm -1 (Figure S2b), while the FBRs 2m and 4a in the elliptic range appear at = 991 cm -1 (FigureS2c).The ribbons and air gaps widths were measured by AFM.

II. Analysis of the Anisotropic PhP Resonances in the Elliptic -MoO3 Frequency Range
The FTIR relative reflection spectra (red dots in FigureS3a), hereinafter referred to as δR, of nanoresonators in the elliptic -MoO3 frequency range was measured for a 110nm-thick -MoO3 flake placed on top of an array of Au-ribbons (with a period of 2.71 μm and a ribbon width of 1.48 μm, corresponding to d = 1.23 μm) at φ = 0°.For comparison, full-wave numerical simulations of δR (black curve in FigureS3a) were carried out mimicking the experiment (in both cases using plane wave illumination polarized across the ribbons).Although the agreement between experiment and theory is reasonably good, the spectral resolution in our FTIR measurements (2 cm -1 ) is too low to capture the TO-phonon along the [100] crystal direction at  = 998.7cm - .The colorplot in FigureS3b shows δR as a function of frequency and inversed ribbon width (while the ratio between the ribbon width and the width of the air gap is fixed to 1.48/1.23).
In our system, we consider the FPRs in both regions of the flake (above the metal and above air) as independent (i.e.completely neglecting any coupling between them).In this approximation, we can write the phase matching conditions 1 : where ka and km are the wavevectors in the regions above the air gaps and the metal ribbons, Фa and Фm are the reflection phases from the air/metal and metal/air boundaries, and na and nm represent the number of polariton wavelengths fitting in the areas above the air gaps and metal ribbons, respectively.We take the wavevectors in Eqs.(1,2) from the analytical equation for the dispersion of anisotropic modes Ml in a thin biaxial slab of thickness d, propagating in the plane at an angle φ with respect to the crystallographic direction 2 : =  .
For simplicity, we neglect the reflection phases (Фa = Фm = 0).For each number na and nm and fixed relation w/d, the wavevector ka,m is a function of a single parameter, w, which we represent by the solid and dashed curves in FigureS3b (labeled as "na" and "nm", respectively).
Analogously to the case analyzed in the main text for the hyperbolic -MoO3 spectral range, for the elliptic spectral range we identified the specific FPRs corresponding to the peak positions in the δR spectrum in FigureS3a.To that end we monitor the intersection between a horizontal line representing the inverse ribbon width (1.48μm) and the dispersion curves (Figure S3b).In the experimental spectrum, the resonant peaks at the frequencies matching na = 2nm show the highest intensity.This can be explained by the 2 times larger PhP wavevectors of the M1 mode in the areas above the air gaps compared to those above the metal ribbons (due to the effective doubling of the flake thickness caused by the mirror effect of the metal).In contrast, in the hyperbolic range the modes are intrinsically different: while the regions above the air gaps support the M0 mode, the regions above the metal ribbons support the M1 mode.Consequently, the doubling effect is not observed in the hyperbolic frequency range.Note that the M0 mode in the flake above the air gaps in the elliptic range is not supported since the dielectric permittivities εx and εy are both positive.
FigureS3c shows the real part of the out-of-plane component of the vertical electric field across the ribbons as a function of frequency.The number of field oscillations across the ribbon coincides with the resonance numbers, na and nm, labelling the resonance peak in Fig S3a .In FigureS3d, we illustrate the measured δR spectra for  = 0°, 15°, 30°, and 45°, in which the positions of the resonances 2m and 4a agree well with the simulated spectra for all angles  (color plot in FigureS3e).Thus, we also see a strong dependence of the resonant frequencies on the rotation angles within the elliptic range, both in the experiment and theory.In the main text probing the dispersion surface of the M0a PhPs mode (Figure 4) in the hyperbolic regime is performed for the 1a FBR in the area of the flake above air gaps.At the same time, the resonance 6m (the M1m mode) in the area of the flake above the metal ribbons can also be achieved at the same frequency as the resonance 1a for all twisting angles, .Thus, the same resonant peaks in the reflection spectra (Figure 2d) can be interpreted both as the 1a FBR of the M0a mode and the 6m FBR of the M1m mode.Here, we perform probing the dispersion surface of the M1m PhPs mode (Figure S4) in the hyperbolic regime on the example of the 6m FBR in the area of the flake above the metal ribbons.In Figure S4a the experimentally-measured positions of the FBRs are represented as green, black, blue and red points in cylindrical coordinates (the resonant frequency as a z coordinate, the inverse ribbon width, 2/, n = 6, as a radial coordinate and the twist angle, , as the polar angle).The purple cylinder with radius 2/, n = 6, represent the constantmomentum surface meaning that we perform the measurements for the same grating with fixed ribbon width, w, for all twist angles and frequencies, whereas the analytical dispersion surface of the M1m mode is represented by a hyperboloid.The projections of the points from Figure S4a on the (kx, ky) plane are positioned very closely to the intersection points between the analytical isofrequency curves (IFCs) for the corresponding frequencies and the cylinder projection on the (kx, ky) plane (Figure S4b).
An alternative probing of the PhP dispersion surface can be realized by constructing the crosssections of the hyperboloid by constant-frequency planes (Figure S4c), i.e.IFCs corresponding to the resonant frequencies in Figure S2a,b.The projections of the experimentally-measured positions of the resonances, represented in cylindrical coordinates by red and black points in Figure S4c, fit well the IFCs.
Our results indicate that the dispersion surface of the M1m PhP mode in the hyperbolic frequency range can be experimentally reconstructed by measuring the resonance 6m in the region of the flake above the metal ribbons.
IV. Dispersion Surfaces of the M1a and M1m PhP Modes above both the Air Gap and the Metal Ribbon Regions in the -MoO3 Elliptic Frequency Range In analogy to the results shown in Figure 4 of the main text for probing the dispersion surface of M0a PhP mode in the hyperbolic range (performed for the 1a FBR emerging in the flake area above the air gaps), we represent in Figure S5 the results of the probing of the dispersion surfaces of the M1a and M1m PhP modes in the elliptic range for the 4a and 2m FPRs in the flake area above the air gaps and above the metal ribbons, respectively.The experimental data for Figure S5a-c was taken from Figure S3d, while the one for Figure S5d-f was taken from Figure S2c.
V. Experimental Quality Factors of the Nanoresonators To evaluate the quality factors (Q) of the nanoresonators, we choose resonant peaks corresponding to the 1a and 6m resonances in the α-MoO3 hyperbolic range and the 4a and 2m resonances in the elliptic range.In order to extract Q, we fit the reflection spectra by two Lorentzian curves.The red symbols represent the experimental data measured, while the gray areas represent the fitted Lorentzian line shape to the experimental spectra (Figure S7).The quality factor can be defined as: where Q is the relation between the resonance frequency  and the full width at half maximum Δω.We used the formula: obtaining the following results for the hyperbolic range:

Figure 1 .
Figure 1.PhPs nanoresonators in -MoO3 defined by placing a pristine -MoO3 slab on top of

Figure
Figure 1d correspond to Fabry-Pérot resonances (FPR) arising from multiple reflections of the PhP

Figure 2 .
Figure 2. Analysis of the PhPs resonances and their twist tuning.a) Measured and simulated

Figure 3 .
Figure 3. Near-field measurements of twist-tunable PhPs nanoresonators.a) Schematic of the

Figure 4 .
Figure 4. Probing the dispersion surface of the M0a PhPs mode.a) Analytical dispersion surface The far-field characterization of our samples was carried out by performing FTIR measurements on an -MoO3 flake with thickness t = 110 nm placed on top of Au ribbons with a width w = 1.48 μm and a separation d = 1.23 μm (FigureS1a).Within the hyperbolic frequency range we have identified the FPRs 1a and 6m, matching the same resonance frequency (= 909.5 cm -1 , see Figure1d of the main text).In contrast, in the elliptic range we have identified the resonances 2m and 4a close to = 995 cm -1 (Figure1d in the main text).The analysis of the resonances 2m and 4a for different twisting angles  is presented in FigureS3d of Section II.

Figure S1 .
Figure S1.Optical images of the metal gratings.a) Au grating with a ribbon width w = 1.48 μm and a separation (air gaps) d = 1.23 μm.The  flake placed on top of them has a thickness t = 110 nm.b) Sets of Al gratings with different ribbon widths and air gaps.The substrate is in all cases CaF2.

Figure S2 .
Figure S2.Far-field spectra of PhP nanoresonators in the case of using Al gratings.a) Experimental and simulated relative reflection spectra in the hyperbolic range for different air gap widths, d, and flake rotation angles, , showing resonances (FPRs 1a and 6m for the M0 and M1m PhP modes, respectively) at the same frequency, ω = 908.5 cm -1 .b) Analogous spectra as in (a) for gratings with a resonant frequency ω = 911 cm -1 .c) Experimental and simulated relative reflection spectra in the elliptic range for nanoresonators defined on gratings with different ribbon widths w and twisting angles .The nanoresonators show maximum of their relative reflection at the same frequency ω = 991 cm -1 (resonances 2m and 4a for the M1m and M1a PhP modes, respectively).Solid and dash lines in (a-c) represent experimental and simulated data, respectively.The thickness of the flake in was 127.5 nm.

Figure S3 .
Figure S3.Analysis of the PhP FBRs in the elliptical range and their twist-tuning.a) Measured and simulated relative reflection spectra,  (dash red and solid grey curves, respectively), for a ribbon width w = 1480 nm, separation distance d = 1230 nm, and twist angle  = 0°.b) Simulated  as a function of frequency  and the inverse ribbon width, 1/w, for a fixed ratio w/d.c) Simulated field distributions of the M1a and M1m PhP modes in the -MoO3/air and -MoO3/Au regions, respectively, as a function of  and the x coordinate (across the ribbons).d) Measured  spectra for the twist angles  = 0, 15, 30, and 45° (red, blue, black, and green, respectively).e Simulated relative reflection as a function of  and .f) Schematics of the top view of one lattice unit cell.

Figure S5 .
Figure S5.Probing the dispersion surfaces of the M1a and M1m PhP modes in the elliptic range.a) Analytical dispersion surfaces of the M1a and M1m PhP modes crossed by two cylinders representing constant momenta.The color dots mark the positions of the measured FPR peaks, 2m (near the yellow surface of the M1m PhP mode) and 4a (near the blue surface of the M1a PhP mode) in cylindrical coordinates (resonant frequency from Figure S3d; either the inverse ribbon width, 2πn/w, n = 2, or the inverse air gap width, 2πn/d, n = 4, for 2m and 4a resonances, respectively; and the twisting angle, ) for the same Au grating twisted at different angles, .b,e) Isofrequency curves for the M1m mode at different frequencies, . c,f) Isofrequency curves for the M1a mode at different frequencies, . d) Analogous dispersion surface as in (a), but crossed by a plane of a fixed frequency, ω = 991cm -1 .The points mark the 2m and 4a FBR peak positions for Al gratings (resonant frequencies from Figure S2c, the inverse ribbon width of the grating, 2πn/w, n = 2, and the inverse air gap width of the grating, 2πn/d, n = 4, for 2m and 4a resonances, respectively, and the twisting angle, ).The thicknesses of the flakes in (a,b,c) and (d,e,f) are 110 and 127.5 nm, respectively.
angle °    cm -1  ∆cm -1 Analogously, we fitted (in this case by one Lorentzian curve) the resonances measured in the elliptic range obtaining the following values: