Subdiffractional focusing and guiding of polaritonic rays in a natural hyperbolic material

Uniaxial materials whose axial and tangential permittivities have opposite signs are referred to as indefinite or hyperbolic media. In such materials, light propagation is unusual leading to novel and often non-intuitive optical phenomena. Here we report infrared nano-imaging experiments demonstrating that crystals of hexagonal boron nitride, a natural mid-infrared hyperbolic material, can act as a ‘hyper-focusing lens' and as a multi-mode waveguide. The lensing is manifested by subdiffractional focusing of phonon–polaritons launched by metallic disks underneath the hexagonal boron nitride crystal. The waveguiding is revealed through the modal analysis of the periodic patterns observed around such launchers and near the sample edges. Our work opens new opportunities for anisotropic layered insulators in infrared nanophotonics complementing and potentially surpassing concurrent artificial hyperbolic materials with lower losses and higher optical localization.

= 1541 cm -1 . Remarkably, this device gives a /33 focal spot with a focal distance of /6 (1050 nm, hBN thickness). c, Line profile of taken along the red dotted line in (b). The black double-arrow indicates the FWHM measured as the focal spot size. In addition to the Au disks, we have prepared and imaged samples with more complicated shapes, for example, in the form of nano-lettered logos of our research institutions (d, e) at  = 1610 cm -1 (Supplementary Figs. 5d-e). We also see an approximately 1:1 copy of the pattern underneath hBN, in agreement with the theoretical picture of the directional propagation of the polariton rays almost normal to the surface, Fig. 3d (bottom) of the main text. Interestingly, some parts of the image have higher intensity than others. The origin of these intensity variations may be studied in future experiments. Scale bar: 1m.

Supplementary Note 1: Theoretical model for focusing by a slab of a hyperbolic material
The problem we want to solve is computing the distribution of the electric field ⃗ inside a hyperbolic medium (HM) slab when a metallic disk is positioned next to its surface and the entire system is subject to a uniform in-plane field ⃗ 0 . The schematics of the model are shown in Supplementary Fig. 1a.
We assume that the slab with the axial permittivity Re > 0 and the tangential permittivity Re 0 occupies the region 0 and refer to it as medium 1.
The media above and below the slab have isotropic permittivities 0 and 2 , respectively. To simplify the analysis, we assume that the metallic disk is infinitely thin. Of course, disks used in the actual experiments are of a finite thickness. To approximately account for the finite thickness, we choose the position of the disk in the model to be some distance ℎ away from the bottom surface of the slab, i.e., in the = −ℎ plane, which is inside medium 2. We choose the center of the disk to reside on the -axis. The external field ⃗ 0 is taken to be in the -direction.
If the disk radius is much smaller than the free-space photon wavelength, it is sufficient to use the quasi-static approximation for the electric field, where the scalar potentials Φ obey the following conditions. Potentials Φ 0 and Φ 2 satisfy the Laplace equation in media 0 and 2 while potential Φ 1 satisfies the equation inside the slab. The boundary conditions are: Φ 2 = 0 on the disk, Φ ≃ − 0 at infinity for all , and at large distances. Quantity can be recognized as the dipole moment acquired by the disk. At the disk edge potential ϕ has a square-root singularity; hence, the electric field has an inverse square-root divergence. These properties are familiar from classical electrostatics.
In order to make use of Supplementary Equations (4)-(5) our first step is to generalize them for the case of an anisotropic uniaxial medium. As obvious from Supplementary Equation (2), this can be achieved by rescaling of the axial coordinate: → − tan , where (cf. equations (2) and (3) of the main text). If the imaginary parts of the permittivities are negligibly small, the rescaling factor is pure imaginary. This changes the nature of the solution qualitatively. It becomes possible for the arguments of the square roots to vanish not only at the edge but also at a set of points in space whose coordinates satisfy the equation This set of points is a union of two cones of the opening angle | | coaxial with the disk and passing through its edge. The apices of the cones are located at z = ±a cot .
The electric field has the inverse square-root divergence on these cones, which is consistent with the physical picture of hyperbolic phonon polaritons (HP 2 ) launched predominantly at the disk edge. Similar high field intensity cones have been previously studied in plasma physics and dubbed "resonance cones" 2  size. We will discuss numerical estimates of this parameter for hexagonal boron nitride (hBN) and compare them with our experimental results shortly below.
Our next step is to use Supplementary Equation (4) valid for an unbounded medium as a building block for constructing a solution for the case of a finite-thickness slab. We follow the procedure standard in the method of images and consider an approximate solution as follows: is the reflection coefficient at interface of the slab and medium = 0 2. In the absence of damping in the system, these coefficients are complex numbers of unit modulus, i.e., they are phase factors. The relation of these phases to the parameter used in the main text (equation (5)) is The top line in Supplementary Equation (11)  Although the complex amplitudes of the consecutive image terms, e.g., 0 2 do not decay by the absolute value, images with higher are more distant from the slab. As a result, they produce progressively weaker potential inside the slab, which ensures convergence of the series. By construction, the potential given by Supplementary Equation (11) meets the boundary conditions at the slab surfaces. The image terms modify the asymptotical value of the electric field. However, if we add the normalization factor = 1 − 0 2 / 1 − 0 , we bring it back to 0 . This way, the boundary condition at infinity will also be satisfied. Unfortunately, the potential given by Supplementary Equation (11) violates the equipotential boundary condition on the disk and there is no simple way to remedy that. Similar difficulty appears in the classic electrostatic problem of a circular parallel-plate capacitor of finite radius 1 . The potential distribution in such a capacitor is not equal to the sum of the potentials created by each charged plate in isolation. However, it does look qualitatively similar.
In particular, the square-root singularity at the edge is the universal feature, which must also be exhibited by the exact solution. We surmise that in our problem the main inaccuracy of Supplementary Equation (11) is the strength of the square-root edge singularity. Otherwise, our approximate solution for Φ 1 should capture the qualitative aspects of the radiation cones emanating from the edges correctly. Note that the influence of the s-SNOM tip is not considered in our theory since its effect is expected to be only quantitative, and not qualitative. As is well established in literature, the qualitative features of the near-field contrast are adequately described assuming that the near-field signal registered by the tip is proportional to the electric field just above the surface of the sample 3-6 .
To get the desired electric field component , we take the derivative of Φ 1 with respect to , which we can easily do analytically. The results are plotted in the false color in Fig. 3c (The absolute value of δ is taken in this formulas because can and indeed is negative in the cited example.) It is clear from the derivation that the described model has a scaling property: if and are both multiplied by the same factor, the electric field as a function of dimensionless coordinates: / and /| | does not change. Accordingly, the qualitative aspects of the electric field distribution are controlled by the dimensionless ratio /| |. For example, the -index of the smallest ring is the integer closest to 1/2 − /| |. The evolution of the strength and relative arrangement of the hot rings is illustrated in Fig. 3b of the main text for several / |δ| ratios and also in Fig. 3d  Assuming also that the effective disk-slab separation is ℎ = 25 nm, one half of the disk physical thickness, we get 0 nm from Supplementary Equation (10). Hence, the distance between the maxima of Im occurring on the opposite sides of the focal spot is 2 140 nm (Supplementary Fig. 1b). When comparing this estimate with the focal spot size in the experimental images, such as Fig. 2b, one has to keep in mind that the measured quantity is not Im but the demodulated scattering amplitude , which is a certain functional of . Also, in the experiment the incident light has a mixture of different polarizations whereas linear polarization is assumed in the model. These are likely reasons for minor differences between the experimental images and the Im curves in Supplementary Fig. 1b. For example, the calculated Im vanishes along the line = 0 but the scattering amplitude measured in the experiment does not exactly vanish along any direction. In practice, our procedure to extract the focal spot size from the images was to take a linear cut along the direction where the signal looked like a simple peak and then take the full width at half-maximum (FWHM) of this peak (Supplementary Fig. 5b). The above number 2 140 nm gives a rough theoretical estimate of the focal spot size obtained by this procedure. It is in fact in agreement with the FWHM of 185 ~ 210 nm determined for the s-SNOM images (Figs. 2b in the main text and Supplementary  Fig. 5b).
To make connection with the second part of the article we use the fact that the real-space distribution of the scalar potential can be alternatively represented by a two-dimensional (2D) Fourier integral. After some lengthy but straightforward calculation, it is possible to show that Supplementary Equation (11) implies that the electric field just above the top surface of the slab can be written as = sin The product of the last two terms in the integrand of Supplementary Equation (15) represents the amplitude of the Fourier harmonics with momentum ( ). The absolute values of these amplitudes exhibit slow power-law decay −1 .
Accordingly, the Fourier spectrum of the electric field induced in our system is very broad. The inverse square-root divergence of the field at the "hot rings" arises from the constructive interference among the harmonics of this broad spectrum. The dominant contribution comes from the harmonics with momenta near the discrete set of values (same as equation (5) of the main text) at which function exhibits pole singularities. These are the momenta of the guided modes of the slab. This is why we made a statement in the main text that the image formed on the top surface of the slab can be viewed as a coherent superposition of multiple guided waves launched by the disk.
It is worth pointing out that function has the meaning of the transmission coefficient of HP 2 between the metallic disk and the top surface of the slab. (More precisely, a transmission coefficient multiplied by the momentum-independent factor , mentioned above, that enforces the normalization 0 = 1. ) This transmission coefficient includes the factor − due to evanescent decay across the vacuum gap and the Fabry-Pérot-like resonant transmission factor due to the free propagation of HP 2 inside the HM slab. If a slab were made from a non-HM, parameters δ, , and so the guided wave momenta defined by Supplementary Equation (17) would be predominantly imaginary. As a result, instead of the resonant transmission HP 2 going through the slab would suffer yet another exponential decay.
Finally, let us now discuss another interesting effect, which may perhaps be verified by future experiments. The finite thickness of real metallic disks is crudely accounted for in our model by choosing a nonzero disk-slab distance h. This parameter enters the last exponential factor in Supplementary Equation (16) thereby imposing a soft momentum cutoff ℎ −1 in the integral of Supplementary  Equation (15). Accordingly, the characteristic number of the guided modes that can effectively contribute to the image formation on the opposite surface of the slab becomes limited to Therefore, as ℎ increases, the high-guided waves are progressively eliminated. Eventually, at ℎ |δ|, the amplitude of even the lowest-momentum = 0 mode becomes exponentially small. However, since the higher-order modes are suppressed even more, this mode dominates the spatial oscillations of the field. The corresponding period 2 /| 0 | = | / | is in general incommensurate with and several times larger than |δ|, the repeat distance of the "hot rings" at ℎ = 0.
In Supplementary Fig. 1b we illustrate these trends by numerical simulations of the electric field profile in the "surface-focusing" case | / | = 0.5 for several different h, with other parameters kept the same as in Fig. 3c of the main text. We choose to plot the imaginary part of the electric field, i.e., the field component that is /2-out of phase with respect to the external field 0 . (The real part, i.e., the in-phase component, shows a similar behavior.) Two facts that agree with the analytical picture are readily apparent from Supplementary Fig. 1b. First, at nonzero disk-slab separation the magnitude of the field can be much smaller than for disk right next to the slab. Actually, to show ℎ > 0 curves clearly we have to scale them by the factor of ten in Supplementary Fig. 1b. Second, as ℎ grows, the sharply peaked extrema of the electric field producing the equidistant series of "hot rings" separated by | | gradually transform into smooth sinusoidal oscillations of a larger period. This period is approximately four times | |, in agreement with the expected result | / |. (Recall that −0.23 in this example.)

Supplementary Note 2: Guided wave dispersion and its thickness dependence
The thickness dependence of the = 0 guided waves in hBN has been studied in our previous work 5 . Here we first summarize the procedure and then apply it to analyze the newly discovered 0 guided wave modes.
The momenta of the guided waves given by Supplementary Equation (17) have been defined above as the poles of the transmission coefficient . These momenta are complex, which simply means that the guided waves exhibit finite damping. As usual, the real part of each determines the wavelength = 2 / Re , while the ratio of the imaginary and real parts specifies the loss factor = Im / Re . Instead of the transmission coefficient , one can examine the surface reflectivity because the poles of the two functions coincide, see Supplementary Equations (13), (16), and (19). However, the reflectivity becomes more useful for the second method of determining the mode dispersions. Namely, the imaginary part of the reflectivity as a function or real momentum is always positive and if the loss factor is small enough, it has sharp maxima at Re . Finding such maxima is easy to implement numerically. For the case = 0 this method 7 becomes advantageous at low momenta where the quasi-static approximation becomes inaccurate and Supplementary Equation (19) has to be replaced by a more complicated expression based on the full Fresnel formulas. However, higher-order modes studied here possess rather high momenta, so either method can be used.
Let us now discuss how this procedure applies specifically to the thickness dependence of the guided wave spectra in hBN. The false color plot of Im at a representative IR frequency  = 1400 cm -1 is shown in Supplementary Fig. 2. It has been calculated using the permittivity functions of hBN and SiO 2 from ref. 5 and 8 as input parameters. The bright lines in this plot give . Whereas Supplementary Equation (17) predicts a strictly linear dependence of on crystal thickness , the bright lines in Supplementary Fig. 2 exhibit a slight curvature at small momenta. Actually, at experimentally relevant momenta the difference between the two methods of determining is negligible. These experimental results are shown by dots and triangles in Supplementary Fig. 2a. They have been found using the Fourier-transform method described in the main text. The measurements were done on several different hBN specimens. The thickness for each of the specimen was measured via the AFM topography simultaneously with the IR images. Matching the symbols with the nearby lines makes it possible to identify the modes as l = 1 for the dots and l = 2 for the triangles. Within the shown range of d, the measured wavelengths of these modes scale linearly with the hBN thickness, in accord with the theory. Supplementary   Figure 2b shows another set of results for the frequency  = 790 cm -1 , which belongs to the lower stop-band (Type I region). Here the polaritons have a negative group velocity but they obey the same linear scaling with as in the upper band, which is again in agreement with the theory.