Analyzing and quantifying symmetry breaking of anisotropic shear polaritons in monoclinic crystal slabs

Abstract

Anisotropic phonon polaritons in materials with high symmetry including hexagonal, trigonal and orthorhombic crystals, enable nanoscale light confinement and manipulation, crucial for nanophotonic on-chip technologies. With reduced material symmetry, monoclinic crystals have been found to endow more intriguing polaritonic phenomena, including axial dispersion and asymmetric propagation patterns, offering new freedoms of nanoscale field manipulations. However, such symmetry-broken anisotropic phonon polaritons, stemming from the intrinsic broken mirror symmetry of dispersion, so far have been observed only in bulk natural materials. Here, we unveil shear polaritons in monoclinic crystal slabs with finite thickness through dispersion analysis and field simulation, and quantify the associated symmetry breaking by introducing a new shear factor based on the k-space integral of the rate of energy dissipation. We demonstrate axial dispersion, mirror symmetry-broken dispersion and asymmetric propagation of shear polaritons. By identifying surface-and volume-confined shear polaritons, we reveal their opposite intensity distribution. In addition, the thickness of monoclinic crystal slab can induce the asymmetry transition, which the surface shear polaritons dominate. Our findings expand the shear polariton platform and provide valuable strategies for manipulating light at the nanoscale through symmetry breaking.

---

1. Introduction

Breaking symmetries is fundamental for the precise control over light propagation.^1,2 A landmark example is the hyperbolic optical response observed in certain polar crystals, where strong phonon-induced anisotropy leads to opposite signs in the permittivity components along different principal axes.^3–7 When interacting with light, the lattice vibrations in such polar materials strongly couple with infrared light, yielding hyperbolic phonon polaritons with hyperbolic iso-frequency contours (IFCs) in momentum space.^8,9 Hyperbolic phonon polaritons exhibit directional, raylike propagation, extreme field confinement and ultralow loss,^10–15 thus enhancing light–matter interactions. In recent years, they have attracted wide attention in subdiffractional focusing and imaging, infrared sensing, and superior light manipulation at the nanoscale.^16–19

A further degree of asymmetry emerges in monoclinic crystals, which not only exhibit three principal axes with distinct lattice constants but also support non-orthogonal microscopic dipolar polarizability in one mirror plane. For example, the monoclinic angle is 103.7° for beta-phase gallium oxide (bGO)^20–22 and 92.1° for cadmium tungstate (CdWO₄).^23,24 Compared with hexagonal, trigonal and orthorhombic crystals, these lower-symmetry crystals (LSC) support a new family of anisotropic phonon polariton, so-called shear polaritons, featured with mirror symmetry breaking. The resulting shear polaritons are characterized by microscopic shear phenomena, such as rotation of their principal axis dependent on frequency, i.e. axial dispersion.^25,26 In addition, the asymmetric loss distribution of shear polaritons leads to enhanced propagation directionality, with energy preferentially guided along specific orientations.²⁷ The exotic light propagation of shear polaritons has brought the light control in nano-optics to a new degree. Recently, the investigation of shear polaritons in natural low symmetry materials has been extended into layered crystals,^28–31 promoting the development of advanced optical devices in subwavelength field. However, the analysis of the symmetry breaking of anisotropic shear polaritons in monoclinic crystal slab has not yet been carried out.

In addition, quantifying the symmetry breaking of shear polaritons is essential for analyzing their intrinsic properties and understanding the microscopic shear phenomena. To date, such effects have mainly been characterized through loss redistribution of shear polaritons, typically using physical quantities such as damping factor²¹ or reflection coefficients.²⁷ However, this quantifying method only evaluates symmetry breaking along the arm of the hyperbolic branches, ignoring the contribution from other regions in momentum space. For the related light–matter interactions such as Purcell factor^32,33 and near-field thermal radiation,^34–37 the in-plane momentum space integrals are significant to understand the exotic physical phenomena induced by the broken symmetry.

Here, we introduce and explore the anisotropic phonon polaritons in low symmetry monoclinic crystal slabs with finite thickness. We theoretically derive the analytical dispersion relations of the thin film, based on which the real and imaginary parts of in-plane wavevector are calculated to directly understand the mechanisms of propagation loss redistribution of shear polaritons in monoclinic crystal slab. By reducing the thickness, the polaritonic coupling makes the surface shear polaritons split into symmetric and antisymmetric branches. We find that the surface-confined and volume-confined shear mode branches respectively exhibit clockwise and counterclockwise redistribution. Moreover, by integrating the rate of energy dissipation, we provide a novel global shear factor to quantify the symmetry breaking of anisotropic phonon polaritons in monoclinic crystal slab, offering an in-depth characterization of shear polaritons. Notably, we demonstrate that varying the slab thickness can drive an asymmetry transition, revealing a thickness-dependent control over polaritonic symmetry breaking. Our work offers valuable insights into the light-guiding at deeply sub-wavelength scales.

2. Theoretical aspects

The monoclinic crystal exhibits multiple non-orthogonal dipolar excitations of atomic resonance within lattice planes due to low crystal symmetry. The non-orthogonality of the axes (γ in the Fig. 1a, the angle between a and c axes of the monoclinic plane) results in non-diagonal permittivity tensor,^20,23 with four independent components, written in Cartesian coordinates as:where ε_xy = ε_yx. Here, the a and b axes of the monoclinic plane are along the x and z axes, respectively. The y axis is selected as the direction that forms a right-handed coordinate with x and z axes. Generally speaking, if there is no loss, the real-valued symmetric permittivity matrix can be diagonalized to a diagonal matrix via an orthonormal matrix. The corresponding dispersion is mirror symmetric with respect to the principal axes. However, the complex-valued cannot be diagonalized in real space using unitary rotation matrices. This intrinsic non-diagonalizability leads to asymmetric polariton dispersion and propagation, manifesting a fundamental form of symmetry breaking. To demonstrate this phenomenon, we focus on bGO, a representative monoclinic crystal with a relatively large monoclinic angle (γ = 103.7°) that enhances the shear effect (the permittivity tensor of bGO crystal can be seen in SI1).

Fig. 1 Anisotropic shear polaritons in monoclinic bGO crystal slab with a finite thickness. (a) Schematic of LSC slab with the thickness of d. Two crystalline directions in x–y plane, i.e., a axis and c axis, are non-orthogonal. (b) Axial dispersion, i.e., frequency-dependence of effective optical principal axis orientation. (c)–(f) The distribution of the IFCs in monument space at different frequencies. The thickness of monoclinic bGO crystal slab in (b)–(f) is 100 nm.

To better understand the shear properties in low-symmetry crystal slab structure, we apply Maxwell's equations to analyze their dispersion. The analytic expression of IFCs of hyperbolic shear polaritons, which gives the relation between the in-plane momentum k_ρ and its angle ϕ with respect to the k_x axis, is given by

where ε_d(t) represents the permittivity of superstrate (substrate). l denotes the order of different polaritonic modes, where l = 0 represents the fundamental mode and l = 1, 2, 3… denotes higher-order mode. And , where ε_ϕ is expressed as (see the details in SI2)

ε_ϕ = ε_xx cos² ϕ + ε_yy sin² ϕ + 2ε_xycos ϕ sin ϕ (3)

Using eqn (2), the IFCs at different frequencies are calculated and depicted in Fig. 1c–f. It is evident that the IFCs are asymmetric due to the presence of non-diagonal component, ε_xy. To better analyze the system, we diagonalize the real part of the permittivity tensor of bGO, switching to a frequency-dependent coordinate system [mnz], by rotating the monoclinic plane by the frequency-dependent angle θ(ω):

where positive (negative) value represents counterclockwise (clockwise) rotation, as shown in Fig. 1b. In frequency-dependent mnz coordinate system, a purely imaginary off-diagonal permittivity component is still retained, which is associated with shear phenomena. The rotated permittivity tensor is included in SI. The diagonalization for the real part results in two effective in-plane principal axes, O₁ (m direction) and O₂ (n direction), as indicated by the blue arrows in Fig. 1c–f. These two axes are mutually orthogonal and divide the in-plane momentum space into four quadrants, as highlighted by the shading regions in Fig. 1c. The angle θ represents that between the axis O₁ and the x-axis. The frequency dispersion of the major axes is called axial dispersion.

3. Results and discussion

3.1 The realization of shear polaritons in bGO slab

For higher-symmetry crystal, such as hexagonal (hexagonal boron nitride, h-BN³⁸), trigonal (calcite³⁹), orthorhombic (alpha-phase molybdenum trioxide, α-MoO₃ ^40,41) lattices, the IFCs are mirror symmetry with respect to the principal axes even though the loss is retained. In contrast, the IFCs mirror symmetry with respect to the effective optic axes is broken for the shear polaritons in LSC. To further reveal the shear-induced asymmetry, we calculate the imaginary part of Fresnel coefficient, i.e., Im(r_pp) in monoclinic crystal slab and compare with that in conventional high-symmetry material, such as α-MoO₃ slab whose [100], [001], and [010] crystalline directions are parallel to the x-, y-, and z-axes, respectively, as shown in Fig. 2. The k-space reflection coefficient in monoclinic crystal slab exhibits hyperbolic branches featured with asymmetric loss distribution while that in α-MoO₃ slab is mirror-symmetric. The loss increases in one branch while decreases in the other branch, resulting in asymmetric loss characteristics. This asymmetric distribution makes a significant effect on the polaritonic propagation in the real space, which will be discussed later.

Fig. 2 The analysis of the loss redistribution and asymmetry of shear polaritons in bGO slab. The dispersion of hyperbolic shear polariton in (a) bGO slab and (b) α-MoO₃ slab. The reciprocal space map is characterized by the imaginary part of reflection coefficient, calculated by Transfer Matrix Method, TMM (see details in SI4). The blue dashed curves denote the IFCs based on complex-momentum analysis. (c) and (d) Real part and imaginary part of complex k_ρ of IFCs varying with the azimuthal angle ϕ. The orange vertical lines represent azimuthal angle parallel to the effective optic axis O₂ in (a) and the [001] crystalline direction in (b). The thicknesses of bGO slab and α-MoO₃ slab are 100 nm and 400 nm, respectively. The frequencies are 718 cm⁻¹ and 960 cm⁻¹, respectively.

In addition, to reveal the origination of propagation loss redistribution of shear polaritons, we perform the complex-momentum analysis and plot the real part and imaginary part of k_ρ of IFCs varying with the azimuthal angle, as shown in Fig. 2a and 2c. For hyperbolic polaritons in α-MoO₃, the IFCs is mirror-symmetric (Fig. 2b), also reflecting in Fig. 2d where both the real part and imaginary part of k_ρ is symmetric across the [001] crystalline direction (ϕ = 90°). In contrast, the IFCs of shear polaritons is asymmetric with respect to the effective optic axis (blue curve in Fig. 2a). The corresponding asymmetric real part and imaginary part distribution manifests the mirror symmetry breaking and the shear phenomenon. Further, it can be seen more directly that the imaginary part of wavevector of the hyperbolic branch in quadrant R₁ (the region indicated by red arrows) is larger, making the propagation length less than that of the other branch, thus yielding a more directional propagation.

To directly visualize the effect of this lower symmetry on shear propagation in real space, we perform full-wave simulations to examine dipole-launched polaritonic modes (see details in SI5). The frequency is fixed at 718 cm⁻¹. The shear asymmetry of hyperbolic shear polaritons in real-space propagation is clear in Fig. 3a–d, where the tilted wavefronts propagating at an angle inclined relative to the optic axis. Moreover, the corresponding fast Fourier Transforms (FFT) shown in Fig. 3e–h indicate the redistribution of intensity in momentum space.

Fig. 3 The asymmetric directional propagation of hyperbolic shear polaritons in bGO slab. (a)–(d) Numerical field distribution, that is, real part of E_z, (Re(E_z)), of the dipole-launched polaritons above the surface at ω = 718 cm⁻¹. (e)–(h) The corresponding fast Fourier transform (FFT). The red and while dashed curves represent analytical dispersion of the fundamental and higher-order modes, respectively.

Based on eqn (2), thickness of the slab is a significant parameter for tuning the dispersion contours and the wavevector range of excited hyperbolic shear polaritons. For the case of 50 nm (Fig. 3a and e), the fundamental mode dominates the shear effect which breaks the symmetry of the polariton propagation. With increased thickness, the excited wavevector range of hyperbolic shear polaritons is moved into lower momentum space, so that higher-order mode begins to play a role in inducing the asymmetric propagation of polaritons. Note that the fundamental mode shows a clockwise redistribution of intensity while the l = 1 mode manifests as a counterclockwise intensity redistribution. Correspondingly and in real-space waveforms, the fundamental and high-order modes show counterclockwise and clockwise intensity redistribution, respectively. Yet higher-order mode leads to a slightly less prominent effect when compared to apparent loss redistribution arising from the fundamental mode, which can be seen from the real-space and Fourier space maps. The asymmetric distribution makes the hyperbolic shear polaritons excited exclusively along one dispersion branch, further enhancing their directionality.

3.2 Distinguishing shear polaritons in bGO slab

Fig. 3a–d exhibit the hyperbolic propagation of shear polaritons in monoclinic crystal slab, however, there are two hyperbolic wavefronts centered along two effective optic axes, respectively. Namely, along the hyperbolic wave asymptotes, the real space is divided into four regions, among which fundamental mode determines the hyperbolic propagation in left and right regions and the higher-order modes dominate in upper and bottom regions. It is crucial to distinguish the type of fundamental and higher-order modes excited in monoclinic crystal slab. Note that there are two types hyperbolic polaritons.^39,42 One is surface-confined hyperbolic polaritons (s-HPs) featured with exponential decaying along the direction normal to interface on both sides. The other is volume-confined hyperbolic polaritons (v-HPs) which directionally propagate inside low-loss crystals. The latter, in thin films, is also in the form of waveguide mode. As marked in Fig. 4a, four specific points-corresponding to l = 0, 1, 2 were chosen to distinguish the fundamental and higher-order modes. The normalized electric field distributions corresponding to the selected points are shown in Fig. 4b–e (the details of mode eigenfield calculation can be seen in SI5). For the red triangle and orange circle, the electric field map of Fig. 4b and c corroborates the existence of s-HPs, which both decay with exponential attenuation inside the vacuum and slab medium. Further, the surface hyperbolic polaritons can be classified into symmetric and antisymmetric modes according to relative value of Re(E_x) in the two vacuum-slab interfaces. In the hyperbolic slab, surface polaritons on the two vacuum-medium interfaces couple with each other and split into two branches, consistent with the isotropic surface modes. Finally, the two blue dashed curves converge in the higher-k space. Fig. 4d and e exhibit the characteristics of v-HPs, which are bound within the layer (that is, waveguide mode).

Fig. 4 Distinguishing the type of hyperbolic shear polaritons in bGO slab: surface-confine or volume-confined. (a) The color map of Im(r_pp) for the thickness of 100 nm at frequency of 718 cm⁻¹. The blue and yellow dashed lines show the analytically calculated polaritons dispersion using eqn (2). The red triangle, orange circle, purple square and green rhombus denotes different order modes, which are further analyzed in b–e, respectively. (b) and (c) Left, cross-sectional normalized electric field distribution Re(E_x) of the fundamental mode (marked by the red triangle (b) and the orange circle (c) in (a)) propagating along the surface. The solid black lines denote the air-bGO interface at z = 0 and d. Right, the corresponding field profile at x = 0. (d and e) Left, normalized Re(E_x) distribution of the (d) l = 1 mode (marked by the purple square in (a)) and (e) l = 2 mode (marked by the green rhombus in (a)) propagating along the surface. Right, the corresponding field profile at x = 0.

3.3 Asymmetry transition driven by thickness of slab

Note that the whole momentum space is divided into four quadrants by two effective optic axes, red and green shadow regions, as shown in Fig. 5a. To quantify the shear effect of hyperbolic shear polaritons, previous strategies are based on the calculation of the line integrals along two arms of the hyperbolic shear polaritons hyperbola in each quadrant of the k-space dispersion, including the damping rate or the reflection coefficient. However, this evaluation method just compares the integrals along two arms of the hyperbolic shear polaritons hyperbola, i.e., the IFCs, ignoring the contribution from the rest of k-space region.

Fig. 5 Analyzing the shear effect. (a) The whole k space is divided into four quadrants according to the effective optic axes, such as R₁ and R₂. The surface and volume confined hyperbolic shear polaritons further divide each quadrant (R₁ as an example) into two sub-regions, denoted by R_s1 and R_v1. (b) Analytical shearing factor dependence on the thickness of bGO slab by integrating the rate of energy dissipation in the k-space viaeqn (5). The integration in adjacent quadrant (R₁ and R₂) of the k-space is compared, denoted by η. Moreover, the contributions from two sub-regions are calculated, respectively, denoted by η_s or η_v. The inset shows the point in vacuum where the rate of energy dissipation is calculated above the bGO slab. h denotes the distance from the vacuum/bGO slab interface, which is fixed to 50 nm.

To address this limitation, we introduce a new global shearing factor to quantify the degree of shear through considering the sum of the calculation of surface integrals. We calculate the normalized rate of energy dissipation in the quadrant R of the k-space region, which performs a k-space surface integral, written as:⁴³

where Q is given by

Here, r_αβ represents the Fresnel reflection coefficient which can be obtained using the TMM. The polarization vectors are defined as a_s^± = 1/k_ρ(k_y, −k_x, 0)^t and a_p^± = k_z/(k_ρk₀)(∓k_x, ∓k_x, k_ρ²/k_z)^t where and are the lateral and vertical wave vector components, respectively.

By comparing integrals in the adjacent quadrant of the k-space region, the global shearing factor is defined as

where the indices 1 and 2 represent the R₁ and R₂ quadrant, respectively, illustrated by Fig. 5a. If there is no shear effect, η will be 1, while otherwise, it is asymmetric. We note that eqn (5) is also closely related to enhancement of local density of states or the Purcell factor, which is a direct quantitative measurement on the enhancement of light–matter interactions. As a result, this shearing factor (η) can describe the interplay between the enhancement of light–matter interactions and shear polaritons induced by the symmetry breaking.

To elucidate the shear properties in relation to structural symmetries, we examine the global shearing factors as a function of the thickness of monoclinic crystal slab. The frequency is fixed at 718 cm⁻¹ and the height is 50 nm. It is noteworthy that the thickness induces a transition of asymmetry as shown in Fig. 5b. When d is less than 25 nm, the shearing factor is less than 1, suggesting the rate of energy dissipation in the R₂ quadrant is larger than that in R1 quadrant. As the d further increases, rate of energy dissipation in the R₁ quadrant dominates, leading to the enhanced shear effect. Thus, the dominant quadrant changes from R₂ to R₁ with the increasing of thickness. This undoubtedly enables enhanced control over asymmetry modulation.

Note that hyperbolic shear polaritons can be either volume-or surface-confined. To explore the physics behind the thickness-induced transition, we calculate and compare the contributions from both types of shear polaritons. As shown in Fig. 5a, each quadrant can be further divided into two sub-regions depending on the excitation region of different types of hyperbolic shear polaritons. Correspondingly, the asymmetric factor η_s or η_v is defined as the ratio of the rates of energy dissipation from R₁ and R₂ quadrants, which integrate in surface modes or waveguide mode regions. One can see from Fig. 5b that the surface confined hyperbolic shear polaritons almost dominate the shear effect.

Finally, to further visually understand the asymmetric transition, we plot the imaginary part of reflected Green’ s function Im(G_zz^R) in the momentum space. In Fig. 6a, the bright regions are concentrated in space less than 200k₀, distinct from the fundamental mode regions about 300k₀. This is due to attenuation of evanescent wave in vacuum, which makes the contribution from the fundamental mode of 5 nm monoclinic crystal slab negligible. In contrast, the energy dissipation rate of the non-excited region dominates. For the case of 10 nm slab (Fig. 6b), the brightest regions overlap the fundamental mode regions. If the asymmetry factor is based on the line integral along the arm of hyperbolic shear polaritons hyperbola, the R₁ quadrant dominates the shear factor. However, in term of the rate of energy dissipation adopting the k-space surface integral, the R₂ quadrant dominates. This is because the brightest regions located at R₂ quadrant occupy broader momentum space. As the d further increases, the hyperbolic shear polaritons in R₁ quadrant dominate, resulting in the asymmetry transition.

Fig. 6 The distribution of the imaginary part of reflected Green’ s function in momentum space, Im(G_zz^R(k_x, k_y)) for different thicknesses of bGO slab. The frequency is fixed to 718 cm⁻¹. The white dashed lines are parallel to the effective optic axes, dividing the k space into four quadrants.

Finally, we provide the additional results on the evolution of the shear factor in other polaritonic mode regimes in the Supporting Information (S7).

4. Conclusion

By exploring the anisotropic phonon polaritons in monoclinic crystal plates, we extend the symmetry breaking to the deep subwavelength scale. We reveal that the thickness of the monoclinic crystal slab is crucial in controlling the asymmetry of phonon polariton propagation and controlling their shear phenomena. The rotation of the effective optic axis with frequency indicates axial dispersion. The asymmetric distribution of the real part and imaginary part of in-plane wavevector of hyperbolic shear polaritons further confirms the in-plane mirror symmetry breaking and loss redistribution, which explains the origination of shear effect and asymmetric directional propagation. In addition, we quantify the symmetry breaking of hyperbolic shear polaritons via the surface integral of the rate of energy dissipation in the momentum space. The results show the thickness of monoclinic crystal slab can induce the asymmetry transition, meaning dominant energy dissipation rate jumps between two adjacent quadrants. We also demonstrate that the broken symmetry induced by surface-confined hyperbolic shear polaritons dominate the shear effect, including the asymmetric polariton propagation and defined asymmetric factor. Our finding with tailorable symmetry-broken hyperbolic polaritons in crystal slabs will be important in sub-diffractional imaging and spectroscopy. Additionally, the exquisite control over the degree of symmetry-breaking at deeply sub-wavelength scales provides a flexible strategy for controlling the asymmetric energy flow at the nanoscale, which are readily valuable for various technologies such as infrared sensing, photonic circuits, near-field radiation management and many others.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplementary information: the permittivity of monoclinic bGO crystal; the derivation of dispersion relations of shear polaritons in monoclinic crystal slab; the reflection matrix of the anisotropic material; the simulation method of electric field; the dielectric function of α-MoO₃; the electric-field Green's function; the symmetry breaking in other polaritonic mode regimes in monoclinic crystal slabs. See DOI: https://doi.org/10.1039/d5nr02806j.

Acknowledgements

This work was supported by the China Scholarship Council (No. 202306010144), the National Natural Science Foundation of China (No. 51576004). G. H. acknowledges the Nanyang Assistant Professorship Start-up Grant, National Research Foundation of Singapore through the Competitive Research Program (CRP22-2019-0064) and Ministry of Education (Singapore) under grant AcRF TIER 1 no. RG61/23 and TIER2 (MOE-T2EP50224-0044).

References

1. E. Galiffi, G. Carini, X. Ni, G. Álvarez-Pérez, S. Yves, E. M. Renzi, R. Nolen, S. Wasserroth, M. Wolf, P. Alonso-Gonzalez, A. Paarmann and A. Alù, Nat. Rev. Mater., 2023, 9, 9–28.
2. M. He, T. G. Folland, J. Duan, P. Alonso-González, S. De Liberato, A. Paarmann and J. D. Caldwell, ACS Photonics, 2022, 9, 1078–1095.
3. G. Hu, J. Shen, C. W. Qiu, A. Alù and S. Dai, Adv. Opt. Mater., 2019, 8, 1901393.
4. J. D. Caldwell, I. Aharonovich, G. Cassabois, J. H. Edgar, B. Gil and D. N. Basov, Nat. Rev. Mater., 2019, 4, 552–567.
5. Q. Zhang, G. Hu, W. Ma, P. Li, A. Krasnok, R. Hillenbrand, A. Alù and C.-W. Qiu, Nature, 2021, 597, 187–195.
6. Y. Wu, J. Duan, W. Ma, Q. Ou, P. Li, P. Alonso-González, J. D. Caldwell and Q. Bao, Nat. Rev. Phys., 2022, 4, 578–594.
7. L. Liu, L. Xiong, C. Wang, Y. Bai, W. Ma, Y. Wang, P. Li, G. Li, Q. J. Wang, F. J. Garcia-Vidal, Z. Dai and G. Hu, Nature, 2025, 644, 76–82.
8. M. N. Gjerding, R. Petersen, T. G. Pedersen, N. A. Mortensen and K. S. Thygesen, Nat. Commun., 2017, 8, 320.
9. Y. Zeng, Q. Ou, L. Liu, C. Zheng, Z. Wang, Y. Gong, X. Liang, Y. Zhang, G. Hu, Z. Yang, C.-W. Qiu, Q. Bao, H. Chen and Z. Dai, Nano Lett., 2022, 22, 4260–4268.
10. G. Hu, Q. Ou, G. Si, Y. Wu, J. Wu, Z. Dai, A. Krasnok, Y. Mazor, Q. Zhang, Q. Bao, C.-W. Qiu and A. Alù, Nature, 2020, 582, 209–213.
11. M. He, L. Hoogendoorn, S. Dixit, Z. Pan, G. Lu, K. Diaz-Granados, D. Li and J. D. Caldwell, Nano Lett., 2023, 23, 5035–5041.
12. J. Duan, N. Capote-Robayna, J. Taboada-Gutiérrez, G. Álvarez-Pérez, I. Prieto, J. Martín-Sánchez, A. Y. Nikitin and P. Alonso-González, Nano Lett., 2020, 20, 5323–5329.
13. J. Duan, G. Álvarez-Pérez, C. Lanza, K. Voronin, A. I. F. Tresguerres-Mata, N. Capote-Robayna, J. Álvarez-Cuervo, A. Tarazaga Martín-Luengo, J. Martín-Sánchez, V. S. Volkov, A. Y. Nikitin and P. Alonso-González, Nat. Mater., 2023, 22, 867–872.
14. P. Li, M. Lewin, A. V. Kretinin, J. D. Caldwell, K. S. Novoselov, T. Taniguchi, K. Watanabe, F. Gaussmann and T. Taubner, Nat. Commun., 2015, 6, 7507.
15. D. Lee, S. So, G. Hu, M. Kim, T. Badloe, H. Cho, J. Kim, H. Kim, C.-W. Qiu and J. Rho, eLight, 2022, 2 DOI:10.1186/s43593-021-00008-6.
16. J. Martín-Sánchez, J. Duan, J. Taboada-Gutiérrez, G. Álvarez-Pérez, K. V. Voronin, I. Prieto, W. Ma, Q. Bao, V. S. Volkov and R. Hillenbrand, Sci. Adv., 2021, 7, eabj0127.
17. Z. Zheng, J. Jiang, N. Xu, X. Wang, W. Huang, Y. Ke, S. Zhang, H. Chen and S. Deng, Adv. Mater., 2021, 34, 2104164.
18. M. Autore, P. Li, I. Dolado, F. J. Alfaro-Mozaz, R. Esteban, A. Atxabal, F. Casanova, L. E. Hueso, P. Alonso-González, J. Aizpurua, A. Y. Nikitin, S. Vélez and R. Hillenbrand, Light: Sci. Appl., 2017, 7, 17172–17172.
19. A. Bylinkin, M. Schnell, M. Autore, F. Calavalle, P. Li, J. Taboada-Gutièrrez, S. Liu, J. H. Edgar, F. Casanova, L. E. Hueso, P. Alonso-Gonzalez, A. Y. Nikitin and R. Hillenbrand, Nat. Photonics, 2020, 15, 197–202.
20. N. C. Passler, X. Ni, G. Hu, J. R. Matson, G. Carini, M. Wolf, M. Schubert, A. Alù, J. D. Caldwell, T. G. Folland and A. Paarmann, Nature, 2022, 602, 595–600.
21. J. Matson, S. Wasserroth, X. Ni, M. Obst, K. Diaz-Granados, G. Carini, E. M. Renzi, E. Galiffi, T. G. Folland, L. M. Eng, J. Michael Klopf, S. Mastel, S. Armster, V. Gambin, M. Wolf, S. C. Kehr, A. Alù, A. Paarmann and J. D. Caldwell, Nat. Commun., 2023, 14, 843.
22. J. Álvarez-Cuervo, M. Obst, S. Dixit, G. Carini, A. I. F. Tresguerres-Mata, C. Lanza, E. Terán-García, G. Álvarez-Pérez, L. F. Álvarez-Tomillo, K. Diaz-Granados, R. Kowalski, A. S. Senerath, N. S. Mueller, L. Herrer, J. M. De Teresa, S. Wasserroth, J. M. Klopf, T. Beechem, M. Wolf, L. M. Eng, T. G. Folland, A. Tarazaga Martín-Luengo, J. Martín-Sánchez, S. C. Kehr, A. Y. Nikitin, J. D. Caldwell, P. Alonso-González and A. Paarmann, Nat. Commun., 2024, 15, 9042.
23. G. Hu, W. Ma, D. Hu, J. Wu, C. Zheng, K. Liu, X. Zhang, X. Ni, J. Chen, X. Zhang, Q. Dai, J. D. Caldwell, A. Paarmann, A. Alù, P. Li and C.-W. Qiu, Nat. Nanotechnol., 2022, 18, 64–70.
24. A. Mock, R. Korlacki, S. Knight and M. Schubert, Phys. Rev. B, 2017, 95, 165202.
25. E. M. Renzi, E. Galiffi, X. Ni and A. Alù, Phys. Rev. Lett., 2024, 132, 263803.
26. S. Yves, E. Galiffi, X. Ni, E. M. Renzi and A. Alù, Phys. Rev. X, 2024, 14, 021031.
27. L. Zhou, X. Ni, Z. Wang, E. M. Renzi, J. Xu, Z. Zhou, Y. Yin, Y. Yin, R. Song, Z. Zhao, K. Yu, D. Huang, Z. Wang, X. Cheng, A. Alù and T. Jiang, Nat. Commun., 2025, 16, 2953.
28. Y. Zheng, Y. Hu, L. Liu, Y. Zeng, T. Hou, Y. Yin, Z. Dai, S. Zhu and H. Chen, J. Opt. Soc. Am. B, 2025, 42, 639–644.
29. Y. Zheng, Y. Zeng, Y. Hu, L. Liu, Z. Dai and H. Chen, npj Nanophotonics, 2025, 2 DOI:10.1038/s44310-025-00053-w.
30. H. Zhang, Y. Liu, H. Tang, J. Xu, B. Liu and Z. Qi, Opt. Commun., 2025, 591, 132112.
31. P. Díaz-Núñez, C. Lanza, Z. Wang, V. G. Kravets, J. Duan, J. Álvarez-Cuervo, A. T. Martín-Luengo, A. N. Grigorenko, Q. Yang, A. Paarmann, J. Caldwell, P. Alonso-González and A. Mishchenko, Sci. Adv., 2025, 11, eadw3452.
32. G. Hu, A. Krasnok, Y. Mazor, C.-W. Qiu and A. Alù, Nano Lett., 2020, 20, 3217–3224.
33. A. N. Poddubny, P. A. Belov, P. Ginzburg, A. V. Zayats and Y. S. Kivshar, Phys. Rev. B:Condens. Matter Mater. Phys., 2012, 86, 035148.
34. S. Chen, C. Fu and G. Hu, ACS Photonics, 2024, 12, 271–281.
35. K. Kim, B. Song, V. Fernández-Hurtado, W. Lee, W. Jeong, L. Cui, D. Thompson, J. Feist, M. T. H. Reid, F. J. García-Vidal, J. C. Cuevas, E. Meyhofer and P. Reddy, Nature, 2015, 528, 387–391.
36. S. Chen, X. Wu and C. Fu, Nanoscale, 2023, 15, 17889–17898.
37. S. Chen, C. Zeng, H. Liu, Y. Hu, X. Wu and C. Fu, Int. J. Heat Mass Transfer, 2025, 249, 127263.
38. A. Kumar, T. Low, K. H. Fung, P. Avouris and N. X. Fang, Nano Lett., 2015, 15, 3172–3180.
39. W. Ma, G. Hu, D. Hu, R. Chen, T. Sun, X. Zhang, Q. Dai, Y. Zeng, A. Alù, C.-W. Qiu and P. Li, Nature, 2021, 596, 362–366.
40. S.-J. Yu, H. Yao, G. Hu, Y. Jiang, X. Zheng, S. Fan, T. F. Heinz and J. A. Fan, ACS Nano, 2023, 17, 23057–23064.
41. G. Álvarez-Pérez, T. G. Folland, I. Errea, J. Taboada-Gutiérrez, J. Duan, J. Martín-Sánchez, A. I. F. Tresguerres-Mata, J. R. Matson, A. Bylinkin, M. He, W. Ma, Q. Bao, J. I. Martín, J. D. Caldwell, A. Y. Nikitin and P. Alonso-González, Adv. Mater., 2020, 32, 1908176.
42. S. Chen, X. Wu and C. Fu, Opto-Electron. Sci., 2024, 3, 240002.
43. L. Novotny and B. Hecht, Principles of nano-optics, 2012, Cambridge university press.

---