• Photonics Research
  • Vol. 10, Issue 10, B14 (2022)
Sicen Tao1、2, Tao Hou1, Yali Zeng1、2, Guangwei Hu2, Zixun Ge1, Junke Liao1, Shan Zhu1, Tan Zhang2, Cheng-Wei Qiu2、3、*, and Huanyang Chen1、4、*
Author Affiliations
  • 1Department of Physics, Institute of Electromagnetics and Acoustics, College of Physical Science and Technology, Xiamen Universityhttps://ror.org/00mcjh785, Xiamen 361005, China
  • 2Department of Electrical and Computer Engineering, National University of Singapore, Singapore 117583, Singapore
  • 3e-mail: chengwei.qiu@nus.edu.sg
  • 4e-mail: kenyon@xmu.edu.cn
  • show less
    DOI: 10.1364/PRJ.463611 Cite this Article
    Sicen Tao, Tao Hou, Yali Zeng, Guangwei Hu, Zixun Ge, Junke Liao, Shan Zhu, Tan Zhang, Cheng-Wei Qiu, Huanyang Chen. Anisotropic Fermat’s principle for controlling hyperbolic van der Waals polaritons[J]. Photonics Research, 2022, 10(10): B14 Copy Citation Text show less


    Transformation optics (TO) facilitates flexible designs of spatial modulation of optical materials via coordinate transformations, thus, enabling on-demand manipulations of electromagnetic waves. However, the application of TO theory in control of hyperbolic waves remains elusive due to the spatial metric signature transition from (+,+) to (-,+) of a two-dimensional hyperbolic geometry. Here, we proposed a distinct Pythagorean theorem, which leads to establishing an anisotropic Fermat’s principle. It helps to construct anisotropic geometries and is a powerful tool for manipulating hyperbolic waves at the nanoscale and polaritons. Making use of absolute instruments, the excellent collimating and focusing behaviors of naturally in-plane hyperbolic polaritons in van der Waals αMoO3 layers are demonstrated, which opens up a new way for polaritons manipulation.


    Hyperbolic media feature extreme optical anisotropy where the component of its permittivity tensor can have different signs, which is widely studied spanning from (known as) metamaterials [15] to natural optical crystals [612]. Those properties will allow the unbounded hyperbolic dispersions with large momentum of waves, and, hence, the high confinement of light at the nanoscale for enhanced light–matter interactions. As a result, hyperbolic metamaterials [4,5], artificial materials made of metal-dielectric multilayers or plasmonic nanowire arrays, are exploited for numerous applications including negative refraction [2,4,5], subwavelength imaging [1,13,14], and enhanced spontaneous emission [15,16]. However, the substantial plasmonic losses are inevitable. More recently, the hyperbolic response is found in bulky crystals and thin van der Waals (vdW) films, such as hexagonal boron nitride (hBN) [68] and α-phase molybdenum trioxide (αMoO3) [912] due to their strong excitation of optical phonons, forming hyperbolic phonon polaritons (PhPs). This kind of hyperbolic response, in principle, is immune from Ohmic loss as a result of the absence of electron–electron scattering and offers tremendous opportunities of nanoscale manipulation of light.

    For further controls, one can make the structure of natural crystals. For instance, the structured array of deeply subwavelength hBN gratings [6,7] can allow the in-plane hyperbolic dispersions, facilitating directional propagations of otherwise isotropic polaritons at the interface. Flexible switching of those directions as well as driving hyperbolic to elliptic topological transition can be realized [17,18] using the twisted bilayer configurations. Moreover, the αMoO3 nanocavities for predescribed field patterns and prolonged lifetime are observed in tailoring the orientation of geometric boundaries [19]. Applications of focusing, wavefront engineering, and electrically dynamical modulation of polaritons are also explored [2023]. In general, those cases show that structuring anisotropic media can supply even more powers for extreme and exotic light manipulations, thus, in this paper we provided a generic and systematic approach of such structuring deriving from transformation optics (TO) theory for the on-demand designer polaritonics.

    TO [2426] applies spatial coordinate transformations to design spatially varying and anisotropic material composites based on the form invariance of Maxwell’s equations, enabling various fascinating metaphotonic applications, such as invisibility cloaks [27], field rotators [28], cosmic string simulators [29], and mimicked Hawking radiation [30]. However, spatial-transformation-based TO theory has not been well constructed for hyperbolic media because the spatial metric signature experiences a transition from (+,+) to (,+), which introduces the “pseudotime coordinate,” causing the effect in similar to a space–time exchange [31,32]. A recent work [33] introduces two-dimensional (2D) Clifford algebra to build a theory of hyperbolic conformal TO, but it requires a coordinate weight ratio of space geometry to be 1, which is isotropic to some extent, limiting the material realizability in experiments.

    Here, we propose a new Pythagorean theorem in an anisotropic system and build up an anisotropic Fermat’s principle which provides guidelines for polaritons and nanoscale extreme light manipulation via reshaping the space geometry to hyperbolic. We apply this theory to the absolute instruments [3437], and specifically, we demonstrate the Luneburg lens (LL) for intrinsic collimation and Maxwell’s fish-eye lens (MFL) for interface focusing of hyperbolic PhPs in αMoO3 layers. Since our theory requires on-site and local modulation of polaritonic characteristics, we show that the simple spatial modulation of thickness of the biaxial αMoO3 would be enough to implement the proposed anisotropic Fermat’s principle within reasonable experimental efforts in the future [38,39].


    We first discuss our approach towards the TO theory for hyperbolic media. In the following and unless otherwise specified, without loss of generality, we focus on the 2D transverse magnetic (TM) polarization, which is usually much more confined. Assume a space filled with the following material: ε=μ=[10001000n2(x,y)],where the z component of the permittivity and permeability follows the profile n(x,y) with respect to Cartesian coordinates x and y depicted as the red square grids as shown in Fig. 1(a). We here consider a mapping that can afford the transformation towards an anisotropic space, defined as x=μyx,y=μxy,where μx and μy are real constants that control the deformation of coordinates. According to the TO recipe [2426], we can easily obtain material parameters in the (x,y) coordinate, shown schematically as the red rectangular grids in Fig. 1(b) from Eq. (1), written as ε=μ=[μxμy000μyμx000μxμyn2(μyx,μxy)].

    Schematic of a transformation relation. (a) The original space, which corresponds to the isotropic space w=x+iy (the red square grids). (b) The physical space, which corresponds to the anisotropic space w′=μyx′+iμxy′ (the red rectangular grids). (c) The transformation relations of Re(x′)−μy with μx=1 and μy≥0 will be changed to Im(x′)−μy when μy<0 with μx=1 and x=−2, −1, 1, and 2, respectively (black stipple lines). The contour curve of r′=1 relevantly changes from hyperbola to parallel straight lines to an ellipse with μy from −3 to 0 to 3.

    Figure 1.Schematic of a transformation relation. (a) The original space, which corresponds to the isotropic space w=x+iy (the red square grids). (b) The physical space, which corresponds to the anisotropic space w=μyx+iμxy (the red rectangular grids). (c) The transformation relations of Re(x)μy with μx=1 and μy0 will be changed to Im(x)μy when μy<0 with μx=1 and x=2, 1, 1, and 2, respectively (black stipple lines). The contour curve of r=1 relevantly changes from hyperbola to parallel straight lines to an ellipse with μy from 3 to 0 to 3.

    From this equation, we know that εz=n2μxμy, μx=μxμy, and μy=μyμx for TM polarization. The physical space after the mapping of Eq. (2) can be expressed by a complex number, w=μyx+iμxy.

    The use of a complex number has geometrical significance: the modulus |w|=w*·w gives a generalized radial coordinate, r=μyx2+μxy2,which obeys a new Pythagorean theorem. Euler’s formula form of the complex number w=reiθ=rcosθ+irsinθ results in cosθ=μyxμyx2+μxy2 and sinθ=μxyμyx2+μxy2.

    Obviously, the constructed anisotropic space depends on the values of μx and μy. The new space can be freely shaped into an elliptic geometry if μx>0,μy>0. However, it is not able to obtain a hyperbolic medium from Eq. (3) as the material parameters are all in the form of square roots, which cannot give negative values. Besides, the complex number of Eq. (4) would lose its meaning if μxμy<0. This weird phenomenon is demonstrated in Fig. 1(c). For instance, setting μx=1 (a mapping where the y-axis coordinate will not change), the relationship of the x-axis coordinate versus μy is plotted in Fig. 1(c). When the value of μy varies from positive to negative, the transformation mappings of x=μyx (illustrated as black stipple lines with x being 2, 1, 1,2, respectively) will jump from the Re(x)μy plane to the Im(x)μy plane. This imaginary coordinate actually changes the original space geometry where the conventional TO is defined. A recent work [33] defined a conformal ultrahyperbolic geometry using 2D Clifford algebra. Nevertheless, the weight ratio of coordinates is restricted as 1, i.e., the conic angle of asymptotes is always 90°, this isotropic geometry is quite confined, and the corresponding materials are not easy to achieve in experiments.

    Next, we show that the mapping of Eq. (2) can be used to define a more general geometry, which accommodates both hyperbolic and elliptic cases, and establish an anisotropic Fermat’s principle theory, which is easier for further implementation.

    From Eq. (2), the line element in the anisotropic space can be defined as ds2=n2dl2=n2(μydx2+μxdy2)=εzμydx2+εzμxdy2,where dl=μydx2+μxdy2 is the infinitesimal increment of the anisotropic geometrical path length, consistent with the generalized radial coordinate of Eq. (5). This generalized system actually conforms to an anisotropic Fermat’s principle. In addition, we can find a perfect analogy of this anisotropic Fermat’s principle in mechanics based on the optical-mechanical analogy theory [37,40] (see Appendix A). The trajectory of the particle with unit mass, total energy E, and potential U in mechanics is the same as the trajectory of light rays in the refractive index n(r), n(r)=2(EU),where it is noted that the radial coordinate r of n(r) should be the generalized one, that is, Eq. (5) in this anisotropic system.

    Thus, from Eq. (6), we can rewrite the material parameters of the anisotropic space as ε=μ=[μx000μy000n2(r)],where εz=n2,μx=μx, and μy=μy for TM polarization. The refractive index n=n(r) is a function of Eq. (5), and the factors μx and μy are exactly the in-plane permeability components μx and μy of the anisotropic space, respectively. In addition, the angular frequency ω and wave-vectors kx and ky satisfy the following dispersion relation: ω2=(cnk)2=c2n2(kx2μy+ky2μx).

    This relationship is reminiscent to the isofrequency surface of uniaxial anisotropic metamaterial for the TM-polarized waves [4], only in our case, here, is a 2D model, where c is the speed of light in free space. Note that a similar process can be extended to 2D transverse electric (TE) polarization, requiring the gradual distribution of permeability, which is relatively difficult in practice, although our theory can be applied well. Besides, we can also take the z coordinate into account to construct a 3D space, but it would be too complex and unnecessary for further study as we focus more on the in-plane effects. Based on the anisotropic Fermat’s principle, we can then construct the space into hyperbolic geometry as known from Eq. (8).

    The hyperbolic media require μxμy<0. This will force that the line element ds2=n2|μy|dx2+n2μxdy2 (if μy<0) has one negative sign of the space metric. It is quite similar to an arbitrary (1+1)-dimensional space–time metric, such as ds2=gxxdx2+gttdt2. Here, we consider the negative space coordinate of the hyperbolic metric as a pseudotime coordinate, which suffers a signature transition [3032] from (+,+) to (,+). It reveals that our theory may provide a new insight for studying various curved space–time of general relativity from studying hyperbolic spaces.

    We draw the curves of r=1 with the value of μy changing from 3 to 3 (μx=1) as shown on grating-like yellow planes in Fig. 1(c). An obvious change in the shape from hyperbola to ellipse can be seen clearly, and thereinto, a critical shape of two straight lines emerges at μy=0. These yellow planes are supposed to be the new anisotropic spaces, and the original space appears at μy=1. This phenomenon is similar to the topological transition occurring in the 2D vdW materials, such as graphene strips [41] and αMoO3 [17]. The dispersion curves of two layers of αMoO3 twisted with different angles vary from hyperbola to ellipse, experiencing a transition stage, which forms a canalization state [17]. We next focus on the hyperbolic part and apply our method to the absolute instruments, such as LL and MFL.


    A. Collimating Effect

    Suppose the original space follows an LL device with refractive index n(r)=2r2 where the radius of the device is 1. Here, without loss of generality, we set μx=1, μy=1 in the new space. From Eqs. (5) and (8), we know that εz=[n(r)]2=[2(x2+y2)2]2. It should be noted that the values of μx and μy can be arbitrarily selected to control the lens to form unconventional shapes. As known from Fig. 1(c), the lens boundary r=x2+y2=1 is an open and infinite hyperbola. As illustrated in Figs. 2(a) and 2(d), the underlying purple profile represents the refractive index distribution n(r), which is uniform outside the lens boundary (black solid curves) and gradient inside. The range of the refractive index can go to infinity but is limited by us to less than 5 in the plot. The LL is a symmetric omnidirectional lens that a point source placed at the rim of the LL forms a collimated light on the other side. The geometrical results show that the upward light rays are emitted from the point source (in yellow) at the upper boundary of the hyperbolic lens at (0,1) or (1,2), and propagate as straight lines in the background. Interestingly, after experiencing the same amount of time, the light rays will form a concave wavefront. The downward light rays entering the lens are bent by the refractive index, which is analogous to the effect that a hyperbolic Hooke’s potential U=r22=x2+y22 applied to particles [see Eq. (7) with total energy E=1], reaching the lower boundary with the same direction right before forming the parallel rays and going out of the lens. The isophase plane in the lower background is parallel to the tangent line of the lens boundary, which passes through the source point. The same phenomena occurred in the wave patterns [Figs. 2(b) and 2(e)]. Obviously, our theory can be used to design an arbitrarily shaped hyperbolic LL with their property maintained and further applied to other devices with different refractive index profiles.

    Hyperbolic Luneburg lens with a collimating effect. (a), (d), (b), (e), (c), and (f), respectively, are the geometrical light behaviors, electromagnetic wave pattern [Re(Ez)], and polaritonic wave pattern [Re(Ez)] from the point source at (0,1) or (1,2).

    Figure 2.Hyperbolic Luneburg lens with a collimating effect. (a), (d), (b), (e), (c), and (f), respectively, are the geometrical light behaviors, electromagnetic wave pattern [Re(Ez)], and polaritonic wave pattern [Re(Ez)] from the point source at (0,1) or (1,2).

    With these numerical examples demonstrated, we now discuss possible experimental demonstrations of our theory in controlling hyperbolic PhPs of αMoO3 flakes for nontrivial integrated devices. Bearing the essential concept of index distributions in mind, importantly, we should seek the local modulation of polariton propagation characteristics. For this purpose, we, here, take advantages of large sensitivity of those highly confined polaritons with the ambient background and the thickness of vdW flakes and, hence, consider the waveguide model composed of three layers: air (zd), αMoO3 slab with finite thickness (0zd), and SiO2 substrate (z0) as shown in Fig. 3(b).

    (a) Schematic of the 2D model; (b) schematic of the 3D waveguide model; (c) relation between the effective refractive index neff and the thickness d obtained from the dispersion relation Eq. (11).

    Figure 3.(a) Schematic of the 2D model; (b) schematic of the 3D waveguide model; (c)  relation between the effective refractive index neff and the thickness d obtained from the dispersion relation Eq. (11).

    The field propagates inside the αMoO3 slab and exponentially decays outside, satisfying Maxwell’s equations, {×E=iωμ0H×H=iωε0ε^E.

    The permittivities of air, αMoO3, and SiO2 are ε1, ε^2, and ε3, respectively. The principal components εi of the permittivity tensor ε^2 of αMoO3 can be calculated by a Lorentz model: εi=εi(1+ωLOi2ωTOi2ωTOi2ω2iωΓi), i=x,y,z. Parameter εi is the high-frequency dielectric constant. Parameters ωLOi and ωTOi are the longitude and transverse optical phonon resonance frequencies, respectively. Γi is the broadening factor of the Lorentzian line shape. The relation between real parts of the permittivity of αMoO3 and frequency is illustrated in Fig. 7 (Appendix B), and the corresponding parameters are listed in Table 1 (see Appendix B, where we refer to the data in Refs. [9,11]). The different in-plane anisotropic polaritonic wave modes can be easily excited by changing the operating frequency.

    In biaxial materials, TE and TM guided modes usually cannot decouple. However, the isofrequency surface of TE modes is a closed surface with finite wave vectors, whereas we focus mainly on opening hyperbolic surface, which requires large wave vectors, so the propagation of PhPs is dominated by the TM modes [9].

    Suppose the in-plane PhPs propagate along the [001] direction, i.e., the y axis; then, the electric and magnetic fields can be expressed as E(y,z,t)=eE(z)exp(iqyiωt) and H(y,z,t)=hH(z)exp(iqyiωt) [11] where the in-plane propagation constant is q=neffk0, with neff being the effective refractive index and k0 being the vacuum wave vector. Considering the TM modes (with only field components of Ey, Hx, and Ez) from Eq. (10) by solving the Helmholtz equations and matching the continuous boundary conditions of Ey and Hx at the boundaries z=d and z=0, we obtain the dispersion relation (see Appendix C), kzd=arctan(α1ε1εykz)+arctan(α3ε3εykz)+mπ,where kz=k02εyεyq2εz and α1,3=q2k02ε1,3>0 are the z components of photon momenta in αMoO3, air, and SiO2, respectively. d denotes the thickness of the αMoO3 slab. m=0,1,2 are the orders of the TM modes. Figure 3(c) shows the effective refractive index neff of the first five TM-polarized guided modes as a function of the thickness d obtained from the dispersion relation of Eq. (11). The effective refractive index has a minimum limit Re(ε3) for the reason of the permittivity of the SiO2 substrate Re(ε3)>ε1=1, but it has no maximum limit. Taking practical experiments into account, the first order (m=0, which is plotted with a black curve in order to distinguish with the other four modes) is only considered because the higher modes are usually inhibited by imperfections of the sample edges and the current signal/noise ratio limitation [8].

    We now calculate the thickness distribution of the αMoO3 flake for demonstration of the collimation effect of the hyperbolic LL in vdW polaritons. We use the refractive index of the 2D model [Fig. 3(a) where the thickness of the biaxial crystal is infinite] as the approximation of the effective refractive index of the biaxial slab in the 3D model [Fig. 3(b) where the thickness of the biaxial crystal is finite]. The effective refractive index neff of q=neffk0 can be obtained approximately from the dispersion relation Eq. (9) of the 2D model through letting kx=0, q2=ky2=μxk02n2giving the results of neff=q/k0=μxn=n(r)=2(x2+y2)2, where μx=1. It turns out that the effective refractive index is proportional to the refractive index in the hyperbolic geometry obtained from our anisotropic Fermat’s principle.

    It should be noted that, we need to limit the thickness of αMoO3 to be relatively small to make sure that the value of neff could be larger than the minimum limitation as shown in Fig. 3(c), which is also consistent with the condition of large momentum. To meet the above requirements, we choose to multiply neff by 25 so that the propagation paths are not affected (only the background is elevated as a whole). Accordingly, the in-plane propagation constant changes to q=25neffk0. Substituting q into Eq. (11), we obtain the thickness distribution d(x,y) at a large momentum approximation as illustrated in Fig. 4(a). The color bar shows the range of d(x,y), which is less than 130 nm. The thickness calculation has a certain degree of flexibility relying on the choice of the constant, but for practical consideration, the smaller the total thickness, the slower the change in thickness, and the fewer high-order effects. Also, as long as the thickness varies slowly enough, the adiabatic theorem will not be disobeyed [38].

    Thickness distribution d(x,y) of (a) hyperbolic collimating LL and (b) hyperbolic focusing MFL with an aerial view (left) and top view (right).

    Figure 4.Thickness distribution d(x,y) of (a) hyperbolic collimating LL and (b) hyperbolic focusing MFL with an aerial view (left) and top view (right).

    The obtained collimating polaritonic waves are shown in Figs. 2(c) and 2(f). Here, the frequency of 674  cm1 (band 1) with εx=8.9943+0.0617i, εy=8.9218+0.2421i, and εz=2.8673+0.0014i is adopted to excite the hyperbolic response along the y axis dominated by TM polarization and to obtain the nearly 90° opening angle with |Re(εx)||Re(εy)| at the same time. A dipole source is positioned at (0 μm, 1 μm) or (1  μm,2  μm) and 220 nm above the bottom of the αMoO3 flake. The real part of the z component of electric fields [Re(Ez)] demonstrates the same responses as that in the above panels of geometrical and wave optical results with obviously long-distance collimating polaritonic waves along corresponding directions. It may be useful for its wavefront shaping in practical applications, such as on-chip integration and directional energy transfer. In addition, the range of operating frequency that our proposed model can be well manipulated is the frequency range in the middle part of RB1 (about 560810  cm1) and RB2 (about 865940  cm1) where the in-plane dispersion of αMoO3 is hyperbolic, and the values of the permittivity components are not too extreme for manipulation.

    B. Focusing Effect

    We further demonstrate the application of our anisotropic Fermat’s principle theory to more functionalities to showcase its enabling power to control the light at the nanoscale. Maxwell’s fish-eye lens, distinguishing from the Luneburg lens, focuses light from one point to another at the boundary. Both collimating and focusing properties are widely studied applications of polaritonics, whereas, our theory provides a more precise and systematic way for implementation. The material parameters of MFL should be εz=[n(r)]2=[21+(x2+y2)2]2, still with μx=1, μy=1. The shape of the lens is the same as that of the hyperbolic LL, but the refractive index distribution is different as shown in Figs. 5(a) and 5(d). The emitted light rays simultaneously walk along curved lines from the source point (0, 1) or (1,2) (the yellow point on the upper boundary of the lens) to the imaging point (0,1) or (1,2) (the green point on the lower boundary of the lens), experiencing the same time. The rays continue to propagate in straight lines to the background space in the same way as they do in the upper half-plane background and eventually form a series of hyperbolic isophase curves. The same phenomena also appear in the corresponding wave patterns in Figs. 5(b) and 5(e) and polaritonic wave patterns in Figs. 5(c) and 5(f). The simulation proceeds under the same frequency at 674  cm1 and the same parameters as that in the LL. Only the thickness distribution d of the αMoO3 flake is different. The thickness distribution d is also less than approximately 130 nm but obtained according to neff=21+r2=21+(x2+y2)2 as shown in Fig. 4(b). Here we note that the value of neff is infinite at areas r2=1 and negative at r2<1. To obtain a valid thickness distribution, we set the value of these areas of neff as a large imaginary number, which leads to a near-zero value of d. This treatment is tricky but reasonable because the lights are supposed to be bounded in the region of neff>0, and the area beyond will not affect the light propagation. As shown in Figs. 5(c) and 5(f), one can see the notable wave converging at the image point and the beautiful hyperbolic wavefront in the background areas. This efficient multidirectional subwavelength focusing effect is applicable to planar nano-optical devices.

    Hyperbolic Maxwell’s fish-eye lens with the focusing effect. (a), (d), (b), (e), (c), and (f), respectively, are the geometrical light behaviors, the electromagnetic wave pattern [Re(Ez)], and the polaritonic wave pattern [Re(Ez)] from the point source at (0, 1) or (1,2).

    Figure 5.Hyperbolic Maxwell’s fish-eye lens with the focusing effect. (a), (d), (b), (e), (c), and (f), respectively, are the geometrical light behaviors, the electromagnetic wave pattern [Re(Ez)], and the polaritonic wave pattern [Re(Ez)] from the point source at (0, 1) or (1,2).

    Finally, we briefly discuss the influence of thickness on the above collimating and focusing effects. The variation of the thickness of the αMoO3 slab does not change the open angle (half the angle between the asymptotic lines of the dispersion curves in momentum space) and the topological nature (i.e., hyperbolic or elliptic) of PhPs, and the behaviors of PhPs show a relatively high tolerance to the thickness variation [17]. Our models of continuously changed thicknesses (Fig. 4) also maintain this robustness against thickness variation as a whole. For example, Fig. 6 displays the collimating behavior of polaritonic waves with different maximum thicknesses of the αMoO3 slab (220 nm, 130 nm, and 65 nm, respectively) working at 634  cm1. We can see that the collimating effect is still well behaved and becomes more local with a smaller total thickness for the reason of larger wave vectors, behaviors more remarkable with larger thickness [few scatterings occur in Fig. 6(a) due to the limited computational domain and computer capacity, although they are tolerable]. The focusing effect will show the similar phenomena, which can be simply verified by simulations. We do not show more details here.

    Polaritonic wave pattern [Re(Ez)] of the hyperbolic Luneburg lens with a collimating effect at a frequency 634 cm−1 where the maximum thicknesses dmax of the α–MoO3 layer are (a) 220 nm, (b) 130 nm, (c) 65 nm, respectively.

    Figure 6.Polaritonic wave pattern [Re(Ez)] of the hyperbolic Luneburg lens with a collimating effect at a frequency 634  cm1 where the maximum thicknesses dmax of the αMoO3 layer are (a) 220 nm, (b) 130 nm, (c) 65 nm, respectively.


    In this paper, we proposed a new Pythagorean theorem that constituted an anisotropic Fermat’s principle, which provided a systematic and novel way to manipulate polaritonic waves. The collimating and focusing properties of LL and MFL were perfectly demonstrated with polaritonic waves by solely adjusting the thickness of the αMoO3 flake. The possible way to fabricate the varied thickness is by the focused ion-beam etching fabrication or the chemical vapor deposition system for the growth of αMoO3. By carefully monitoring the thickness of αMoO3 using an optical microscope and atomic force microscope, the curved surface can be made with the peeling off of some of the αMoO3 using the Omniprobe micromanipulator to leave only the expected shape [9,10,17,23]. It is no doubt easier for experiments compared with many other methods, for example, punching holes or designing structures in substrates or vdW materials. We provided a theory basis for designing nanodevices. In principle, the method can be generalized to other structures and other materials with the in-plane hyperbolic response by choosing appropriate hyperbolic bands, such as the natural materials black phosphorus and Weyl semimetal WTe2 [41]. The study of hyperbolic polaritons is of great significance for us to understand the novel phenomena of electromagnetic waves at the nanoscale and the physical principles behind them. Therefore, we expect that our paper opens up new opportunities for manipulating, reshaping, and focusing electromagnetic waves at a subwavelength scale [42] and have some impacts on in-plane or on-chip optics.

    In addition, curved space–time of general relativity can be studied from a distinct view by studying hyperbolic spaces. What is more, our anisotropic Fermat’s principle can be related to the anisotropic Kepler problem (AKP) [43,44], which is an electronic system with an electron having an anisotropic mass tensor. The AKP possesses chaotic behavior, which makes it a hot topic for the study of classical and quantum correspondences. Thus, further research is necessary on the AKP based on our proposed theory in the future, especially for the hyperbolic anisotropic cases. In a word, from our versatile approach, distinctive insights may be inspired for studying space–time of general relativity or anisotropic systems of other fields, such as acoustics, classical, and quantum mechanics.


    Acknowledgment. C.-W.Q. acknowledges the Advanced Research and Technology Innovation Centre (ARTIC), National University of Singapore for financial support.


    We briefly introduce the optical-mechanical analogy and show that our anisotropic Fermat’s principle leads to anisotropic dynamics, which is still in the framework of Newton’s law, or we will call it the anisotropic Newton’s law.

    The optical-mechanical analogy is the analogy between Fermat’s principle of the shortest/longest time and the Maupertuis principle of least action. These two principles all obey the variational principle and are thought to be the specific formulations of Hamilton’s more general principle [40].

    With the help of variational calculus, Fermat’s principle can be expressed with the Euler–Lagrange equation or Hamilton’s equation. One can derive the Euler–Lagrange equation for the extremal trajectory as [37] ddξLv=Lr,where L is the Lagrangian. r is the position, parameterized by ξ. v is velocity, the derivative of r with respect to ξ with v=drdξ.

    The line element of Eq. (6) can be expressed in terms of the Lagrangian as ds=nμydx2+μxdy2=nμydx2dξ2+μxdy2dξ2dξ=Ldξ.

    The Lagrangian is thereby anisotropic and expressed as L=nμydx2dξ2+μxdy2dξ2=nμyvx2+μxvy2=n2v2.

    The parameter ξ plays the role of time but has the dimension of length, which can be physically interpreted as “optical action.” For ray trajectories, the parameter ξ satisfies dξ=dln.

    It is a vital stepping parameter for optical-mechanical analogy. dl=μydx2+μxdy2 is the infinitesimal increment of the anisotropic geometrical path length. From Eq. (A4), the corresponding “speed” is obtained as v=|drdξ|=n|drdl|=n.

    Thus, v=n=μyvx2+μxvy2. Then, substituting the Lagrangian Eq. (A3) into the Euler–Lagrange Eq. (A1), we have d2rdξ2=n22.

    It is anisotropic Newton’s law for a particle with unit mass moving in time ξ with potential, U=n22+E,where E is the total energy of the particle. When μx=μy=1, all quantities and equations above will degenerate to the regular forms.

    Therefore, it establishes the relation between anisotropic Fermat’s principle and anisotropic Newton’s law, The operation means that Newton’s law can be extended to a more general scope to describe the motion in anisotropic backgrounds. It indicates that our theory can be extended to many other anisotropic systems, such as the AKP (to study the chaotic behavior of particles).


    There are three Reststrahlen bands (RBs) (shaded in different colors in Fig. 7) of αMoO3 in the mid-infrared range of 545to1010  cm1: bands 1–3 originate from the in-plane phonon mode along the [001] (y-axis), [100] (x-axis), and [010] (z-axis) crystalline direction, respectively [9]. In bands 1 and 2, ranging from 545 to 851  cm1 and 820 to 972  cm1, in-plane hyperbolic responses with highly confined ultrahigh wave vectors exist with Re(εy)<0, Re(εx)Re(εz)>0, and Re(εx)<0, Re(εy)Re(εz)>0, respectively. In band 3, frequency ranges from 958 to 1010  cm1, where there are Re(εz)<0, Re(εx)Re(εy)>0, and the in-plane dispersion is elliptical.

    Real-part permittivities of α–MoO3 where three different Reststrahlen bands of α–MoO3 are shaded in different colors.

    Figure 7.Real-part permittivities of αMoO3 where three different Reststrahlen bands of αMoO3 are shaded in different colors.

    Parameters for Calculating the Permittivity of αMoO3 [9]

    Parametersx [100]y [001]z [010]


    Considering the TM modes (with only field components of Ey, Hx, and Ez [11]) and substituting the electric and magnetic fields E(y,z,t)=eE(z)exp(iqyiωt) and H(y,z,t)=hH(z)exp(iqyiωt) into Maxwell’s equation Eq. (10), we derive the following equations: iqEzEyz=iωμ0Hx,Hxz=iωε0εyEy,iqHx=iωε0εzEz,which give the Helmholtz equation, 2Hxz2+(k02εyεyq2εz)Hx=0,0zd.

    The equations in the other two layers are similarly expressed as 2Hxz2+(k02ε1q2)Hx=0,zd,2Hxz2+(k02ε3q2)Hx=0,z0.

    From the above Eqs. (C2)–(C4), we suppose kz=k02εyεyq2εz and α1,3=q2k02ε1,3, which are the z-components of photon momenta in αMoO3, air, and SiO2, respectively. The solutions of the above equations are Hx={[Acos(kzd)+Bsin(kzd)]exp[α1(zd)],zd,Acos(kzz)+Bsin(kzz),0zd,Aexp(α3z),z0,where A and B are the coefficients to be determined. Then, the y-component (Ey) of the electric field can be obtained as Ey={α1iωε0ε1[Acos(kzd)+Bsin(kzd)]exp[α1(zd)],zd,kziωε0εy[Asin(kzz)+Bcos(kzz)],0zd,α3iωε0ε3Aexp(α3z),z0.

    By matching the continuous boundary conditions of Ey and Hx at the boundaries z=d and z=0, one can obtain M(AB)=(α3ε3kzεyα1ε1cos(kzd)+kzεysin(kzd)α1ε1sin(kzd)kzεycos(kzd))(AB)=0.The determinant of M should equal zero, detM=0; then, the dispersion relation of hyperbolic PhPs can be obtained: kzd=arctan(α1ε1εykz)+arctan(α3ε3εykz)+mπ.

    It is Eq. (11). This dispersion relation can be expressed more specifically as k02εyεyq2εzd=arctan(q2k02ε1ε1εzεyk02εzq2)+arctan(q2k02ε3ε3εzεyk02εzq2)+mπ.


    [1] Z. Liu, H. Lee, Y. Xiong, C. Sun, X. Zhang. Far-field optical hyperlens magnifying sub-diffraction-limited objects. Science, 315, 1686(2007).

    [2] J. Yao, Z. Liu, Y. Liu, Y. Wang, C. Sun, G. Bartal, A. M. Stacy, X. Zhang. Optical negative refraction in bulk metamaterials of nanowires. Science, 321, 930(2008).

    [3] Z. Jacob, J.-Y. Kim, G. V. Naik, A. Boltasseva, E. E. Narimanov, V. M. Shalaev. Engineering photonic density of states using metamaterials. Appl. Phys. B, 100, 215-218(2010).

    [4] A. Poddubny, I. Iorsh, P. Belov, Y. Kivshar. Hyperbolic metamaterials. Nat. Photonics, 7, 948-957(2013).

    [5] D. Lee, S. So, G. Hu, M. Kim, T. Badloe, H. Cho, J. Kim, H. Kim, C.-W. Qiu, J. Rho. Hyperbolic metamaterials: fusing artificial structures to natural 2D materials. eLight, 2, 1(2022).

    [6] P. Li, I. Dolado, F. J. Alfaro-Mozaz, F. Casanova, L. E. Hueso, S. Liu, J. H. Edgar, A. Y. Nikitin, S. Vélez, R. Hillenbrand. Infrared hyperbolic metasurface based on nanostructured van der Waals materials. Science, 359, 892-896(2018).

    [7] P. Li, G. Hu, I. Dolado, M. Tymchenko, C.-W. Qiu, F. J. Alfaro-Mozaz, F. Casanova, L. E. Hueso, S. Liu, J. H. Edgar, S. Vélez, A. Alu, R. Hillenbrand. Collective near-field coupling and nonlocal phenomena in infrared-phononic metasurfaces for nano-light canalization. Nat. Commun., 11, 3663(2020).

    [8] S. Dai, Z. Fei, Q. Ma, A. S. Rodin, M. Wagner, A. S. Mcleod, M. K. Liu, W. Gannett, W. Regan, K. Watanabe, T. Taniguchi, M. Thiemens, G. Dominguez, A. H. C. Neto, A. Zettl, F. Keilmann, P. Jarillo-Herrero, M. M. Fogler, D. N. Basov. Tunable phonon polaritons in atomically thin van der Waals crystals of boron nitride. Science, 343, 1125-1129(2014).

    [9] Z. Zheng, N. Xu, S. L. Oscurato, M. Tamagnone, F. Sun, Y. Jiang, Y. Ke, J. Chen, W. Huang, W. L. Wilson, A. Ambrosio, S. Deng, H. Chen. A mid-infrared biaxial hyperbolic van der Waals crystal. Sci. Adv., 5, eaav8690(2019).

    [10] W. Ma, P. Alonso-González, S. Li, A. Y. Nikitin, J. Yuan, J. Martín-Sánchez, J. Taboada-Gutiérrez, I. Amenabar, P. Li, S. Vélez, C. Tollan, Z. Dai, Y. Zhang, S. Sriram, K. Kalantar-Zadeh, S.-T. Lee, R. Hillenbrand, Q. Bao. In-plane anisotropic and ultra-low-loss polaritons in a natural van der Waals crystal. Nature, 562, 557-562(2018).

    [11] F. Sun, W. Huang, Z. Zheng, N. Xu, Y. Ke, R. Zhan, H. Chen, S. Deng. Polariton waveguide modes in two-dimensional van der Waals crystals: an analytical model and correlative nano-imaging. Nanoscale, 13, 4845-4854(2021).

    [12] W. Ma, G. Hu, D. Hu, R. Chen, T. Sun, X. Zhang, Q. Dai, Y. Zeng, A. Alù, C.-W. Qiu, P. Li. Ghost hyperbolic surface polaritons in bulk anisotropic crystals. Nature, 596, 362-366(2021).

    [13] Z. Jacob, L. V. Alekseyev, E. Narimanov. Optical hyperlens: far-field imaging beyond the diffraction limit. Opt. Express, 14, 8247-8256(2006).

    [14] A. Salandrino, N. Engheta. Far-field subdiffraction optical microscopy using metamaterial crystals: theory and simulations. Phys. Rev. B, 74, 075103(2006).

    [15] A. S. Potemkin, A. N. Poddubny, P. A. Belov, Y. S. Kivshar. Green function for hyperbolic media. Phys. Rev. A, 86, 023848(2012).

    [16] A. N. Poddubny, P. A. Belov, P. Ginzburg, A. V. Zayats, Y. S. Kivshar. Microscopic model of Purcell enhancement in hyperbolic metamaterials. Phys. Rev. B, 86, 035148(2012).

    [17] G. Hu, Q. Ou, G. Si, Y. Wu, J. Wu, Z. Dai, A. Krasnok, Y. Mazor, Q. Zhang, Q. Bao, C.-W. Qiu, A. Alù. Topological polaritons and photonic magic angles in twisted α-MoO3 bilayers. Nature, 582, 209-213(2020).

    [18] J. Duan, N. Capote-Robayna, J. Taboada-Gutiérrez, G. Álvarez-Pérez, I. Prieto, J. Martín-Sánchez, A. Y. Nikitin, P. Alonso-González. Twisted nano-optics: manipulating light at the nanoscale with twisted phonon polaritonic slabs. Nano Lett., 20, 5323-5329(2020).

    [19] Z. Dai, G. Hu, G. Si, Q. Ou, Q. Zhang, S. Balendhran, F. Rahman, B. Y. Zhang, J. Z. Ou, G. Li, A. Alù, C.-W. Qiu, Q. Bao. Edge-oriented and steerable hyperbolic polaritons in anisotropic van der Waals nanocavities. Nat. Commun., 11, 6086(2020).

    [20] K. Chaudhary, M. Tamagnone, X. Yin, C. M. Spägele, S. L. Oscurato, J. Li, C. Persch, R. Li, N. A. Rubin, L. A. Jauregui, K. Watanabe, T. Taniguchi, P. Kim, M. Wuttig, J. H. Edgar, A. Ambrosio, F. Capasso. Polariton nanophotonics using phase-change materials. Nat. Commun., 10, 4487(2019).

    [21] Z. Zheng, J. Jiang, N. Xu, X. Wang, W. Huang, Y. Ke, S. Zhang, H. Chen, S. Deng. Controlling and focusing in-plane hyperbolic phonon polaritons in α-MoO3 with a curved plasmonic antenna. Adv. Mater., 34, 2104164(2022).

    [22] J. Duan, G. Álvarez-Pérez, A. I. F. Tresguerres-Mata, J. Taboada-Gutiérrez, K. V. Voronin, A. Bylinkin, B. Chang, S. Xiao, S. Liu, J. H. Edgar, J. I. Martín, V. S. Volkov, R. Hillenbrand, J. Martín-Sánchez, A. Y. Nikitin, P. Alonso-González. Planar refraction and lensing of highly confined polaritons in anisotropic media. Nat. Commun., 12, 4325(2021).

    [23] 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, Z. Dai. Tailoring topological transitions of anisotropic polaritons by interface engineering in biaxial crystals. Nano Lett., 22, 4260-4268(2022).

    [24] U. Leonhardt. Optical conformal mapping. Science, 312, 1777-1780(2006).

    [25] J. B. Pendry, D. Schurig, D. R. Smith. Controlling electromagnetic fields. Science, 312, 1780-1782(2006).

    [26] J. B. Pendry, Y. Luo, R. K. Zhao. Transforming the optical landscape. Science, 348, 521-524(2015).

    [27] D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, D. R. Smith. Metamaterial electromagnetic cloak at microwave frequencies. Science, 314, 977-980(2006).

    [28] H. Chen, B. Hou, S. Chen, X. Ao, W. Wen, C. T. Chan. Design and experimental realization of a broadband transformation media field rotator at microwave frequencies. Phys. Rev. Lett., 102, 183903(2009).

    [29] C. Sheng, H. Liu, H. Chen, S. Zhu. Definite photon deflections of topological defects in metasurfaces and symmetry-breaking phase transitions with material loss. Nat. Commun., 9, 4271(2018).

    [30] I. I. Smolyaninov, E. Hwang, E. Narimanov. Hyperbolic metamaterial interfaces: Hawking radiation from Rindler horizons and spacetime signature transitions. Phys. Rev. B, 85, 235122(2012).

    [31] I. I. Smolyaninov, E. E. Narimanov. Metric signature transitions in optical metamaterials. Phys. Rev. Lett., 105, 067402(2010).

    [32] S. Fumeron, B. Berche, F. Santos, E. Pereira, F. Moraes. Optics near a hyperbolic defect. Phys. Rev. A, 92, 063806(2015).

    [33] S. Dehdashti, A. Shahsafi, B. Zheng, L. Shen, Z. Wang, R. Zhu, H. Chen, H. Chen. Conformal hyperbolic optics. Phys. Rev. Res., 3, 033281(2021).

    [34] J. L. Synge. The absolute optical instrument. Trans. Amer. Math. Soc., 44, 32-46(1938).

    [35] M. Born, E. Wolf. Principles of Optics(2006).

    [36] R. K. Luneburg. Mathematical Theory of Optics(1964).

    [37] U. Leonhardt, T. Philbin. Geometry and Light: The Science of Invisibility(2010).

    [38] A. Patsyk, U. Sivan, M. Segev, M. A. Bandres. Observation of branched flow of light. Nature, 583, 60-65(2020).

    [39] H. Gao, B. Zhang, S. G. Johnson, G. Barbastathis. Design of thin-film photonic metamaterial Lüneburg lens using analytical approach. Opt. Express, 20, 1617-1628(2012).

    [40] C. Joas, C. Lehner. The classical roots of wave mechanics: Schrödinger’s transformations of the optical-mechanical analogy. Stud. Hist. Philos. Mod. Phys., 40, 338-351(2009).

    [41] G. Hu, A. Krasnok, Y. Mazor, C.-W. Qiu, A. Alù. Moiré hyperbolic metasurfaces. Nano Lett., 20, 3217-3224(2020).

    [42] Z. Chen, M. Segev. Highlighting photonics: looking into the next decade. eLight, 1, 2(2021).

    [43] Z. Chen, W. Zhou, B. Zhang, C. H. Yu, J. Zhu, W. Lu, S. C. Shen. Realization of anisotropic diamagnetic Kepler problem in a solid state environment. Phys. Rev. Lett., 102, 244103(2009).

    [44] W. Zhou, Z. Chen, B. Zhang, C. H. Yu, W. Lu, S. C. Shen. Magnetic field control of the quantum chaotic dynamics of hydrogen analogs in an anisotropic crystal field. Phys. Rev. Lett., 105, 024101(2010).

    Sicen Tao, Tao Hou, Yali Zeng, Guangwei Hu, Zixun Ge, Junke Liao, Shan Zhu, Tan Zhang, Cheng-Wei Qiu, Huanyang Chen. Anisotropic Fermat’s principle for controlling hyperbolic van der Waals polaritons[J]. Photonics Research, 2022, 10(10): B14
    Download Citation