Interfacing differently oriented biaxial van der Waals crystals for negative refraction

2.1 Proposed structure

Figure 1A and B illustrate the proposed structure and its cross-sectional view (x–z plane), respectively, which comprises air/graphene/α-MoO₃ on an Au substrate. Here, t is the thickness of the α-MoO₃ slab, [100] and [001] are the principal axes of the α-MoO₃ crystal, and θ is the angle between the x-axis and [100] direction. The employment of a metal substrate helps the proposed structure realize significantly higher mode confinement than the conventional HPhP, which is referred to as a hyperbolic image plasmon–phonon polariton (HIPPhP) [23]. The HIPPhP originates from the hybridization of a hyperbolic polariton in the α-MoO₃ slab and their mirror image with an inverted charge distribution in the perfect electric conductor (PEC) [21], which behaves as a mirror. At the mid-IR band, Au can be considered an approximate PEC. Therefore, most of the electric field is confined in the α-MoO₃ slab, which significantly enhances the mode confinement due to the opposing charge distributions of the original and mirror polaritons.

Figure 1:

The structure with single α-MoO₃ slab. (A) Illustration of the air/graphene/α-MoO₃/Au structure where [100], [001], and [010] are the principal axes of the biaxial α-MoO₃ slab and θ is the angle between the x-axis and [100] direction. (B) The x–z plane of the structure, denoted as (A), where t is the thickness of the α-MoO₃ slab.

2.2 Numerical simulations

We specifically investigated the reststrahlen band II (RB II) of α-MoO₃, which encompasses the frequency range of 816–972 cm⁻¹. Within this frequency range, the permittivity (ε) components along the principal axes [100], [001], and [010] exhibit the characteristics of ε _x < 0, ε _y > 0, and ε _z > 0, respectively. The RB is a range of wavelengths where a material exhibits strong absorption due to the resonance between the incident radiation and lattice vibration. Furthermore, an HPhP occurs when the incident radiation couples strongly with the oscillations of the lattice vibrations in the α-MoO₃ slab. Therefore, the RB can give rise to hyperbolic dispersion and the formation of hyperbolic polaritons in the α-MoO₃ slab. Notably, the negative sign of Re (ε _x ) × Re (ε _y ) < 0 implies that the α-MoO₃ has in-plane hyperbolicity in the x–y plane. The complex permittivity of Au can be obtained from the literature [59]. Here, the graphene sheet was modeled as an infinitely thin layer with a surface current density. The surface conductivity of graphene can be calculated by using the Kubo formula [60] (see Supplementary Material, §1). A numerical model was constructed in COMSOL Multiphysics with the mesh resolution set sufficiently fine and the computational domain sufficiently large to ensure simulation accuracy. The permittivity of α-MoO₃ was obtained by employing the conventional Lorentz model [38, 43, 44, 58] (see Supplementary Material, §2.1). The crystallographic directions [100], [001], and [010] of α-MoO₃ coincided with the coordinates x, y, and z, respectively, under the condition θ = 0°. The real and imaginary components of the permittivity of α-MoO₃ were plotted (see Supplementary Material, Figure S1). We placed an electric dipole source with polarization along the z-direction at a position 50 nm above the top graphene layer. The electric fields, represented as Re(E _z ), were calculated at a probing height of 20 nm above the graphene layer. To obtain the corresponding IFCs in wavevector space (k _x , k _y ), we performed a Fourier transform of the Re(E _z ) fields using spatial Fast Fourier Transform (FFT) with spatial sampling resolutions set to N _x = 1000 and N _y = 1000. To minimize undesired reflections, the computational window was enclosed by a perfect-matched-layer boundary condition with a thickness of 1 µm, effectively absorbing the transmitted waves. Simulations were performed to analyze the effects of modifying the α-MoO₃/Au structure with graphene, the effects of interfacing α-MoO₃ slabs at different orientations, and the effects of multiple interfaces.

3.1 Polariton modes with and without graphene

To gain insights into the mode characteristics, the dispersion properties of the α-MoO₃/Au structure at t = 150 nm were analyzed. The principal crystallographic direction [100] was aligned with the x-direction (i.e., θ = 0°). The results indicated an electromagnetic field with unprecedented confinement (see Supplementary Material, Figure S2A). The relationship between the propagation length (L _p ) and t (see Figure S2B) was found to follow the conventional tradeoff between mode confinement and propagation loss, where an increase in t leads to looser mode confinement and longer L _p . The figure of merit [FOM = Re(k _x )/Im(k _x )] was obtained to assess the benefit–cost ratio between mode confinement and propagation loss [22, 46]. High FOM values were obtained at the wavenumbers ω = 880–920 cm⁻¹ (see Figure S2C), which indicates a favorable tradeoff in this frequency range. The group velocity (V _g ) (see Figure S2D) exhibited a similar trend to L _p . The field distributions of the HIPPhP (see Supplementary Material, Figure S3) clearly demonstrated tighter confinement and a shorter wavelength as the thickness decreased.

Next, the dispersion properties of the graphene/α-MoO₃/Au structure were obtained. Increasing the Fermi energy (E _f ) was found to decrease Re(k _x ) with a more pronounced dependence at higher frequencies (see Figure 2A). In contrast, L _p , FOM, and V _g all increased (see Figure 2B–D). A higher FOM indicates a stronger mode confinement or lower mode loss. Interestingly, when E _f is less than 0.1 eV, the case without graphene exhibits the strongest mode confinement but also the largest energy loss. Conversely, when E _f exceeds 0.2 eV, the situation reverses, resulting in a smaller FOM (see Figure 2C). This phenomenon can be attributed to plasmons reaching the interband transition threshold at very low E _f for specific wavevectors, leading to overdamped outcomes [61, 62]. Notably, the value of E _f reaching the interband transition threshold decreases as the wavevector reduces. In contrast, plasmons transition to the intraband transition as E _f increases. As a result, the HIPPhP modes have larger FOM when E _f > 0.2 eV and smaller when E _f < 0.1 eV while comparing to the case without graphene for the condition given in the proposed system. The Re(E _z ) distributions at different values of E _f (see Figure 2E–H) are also shown that both polariton wavelength and L _p increased with increasing E _f . These results offer valuable insights into the influence of graphene doping and its effects on the field propagation characteristics of the HIPPhP.

Figure 2:

Wavenumber ω of an α-MoO₃ slab at t = 150 nm and covered with monolayer graphene at several values of E _f . against (A) Re(k _x ), (B) L _p , (C) FOM, and (D) V _g . E _z field profiles of an α-MoO₃ slab with t = 150 nm: (E) no graphene and covered by monolayer graphene at (F) E _f = 0.1 eV, (G) E _f = 0.3 eV, and (H) E _f = 0.5 eV.

3.2 Effect of different orientations

We first analyzed a single α-MoO₃ slab rotated counterclockwise around the z-axis by the angle θ, as shown in Figure 1A. Here, θ corresponds to the orientation of the crystallographic direction [100] with respect to the x-axis. We used this coordinate system to obtain the anisotropic permittivity tensor of the α-MoO₃ slab (see Supplementary Material, §2.1) under the conditions of ω = 910 cm⁻¹, t = 150 nm, and E _f = 0.15 eV. The Re(E _z ) field distributions varied with θ and demonstrated directional canalization along the corresponding rotation angles (see Figure S4). These results provide compelling evidence of the in-plane hyperbolicity of the HIPPhP. To analyze the wavevector properties of the Re(E _z ) fields, we performed a Fourier transform and plotted the resulting isofrequency contours (IFCs) in wavevector space (k _x , k _y ) (see Figure S5). The amplitude of the IFCs reflected the energy distribution in the corresponding wavevector space. Notably, a stronger amplitude indicates a higher concentration of energy. At θ = 0°, a significant portion of the energy was concentrated around k _y /k ₀ = −100 to 100, while k _x /k ₀ = ±20. The symmetric axis of the IFCs (i.e., k _x /k ₀ = 0) exhibited a tilt as θ was varied, which implies that the directional propagation of energy can be flexibly tuned by adjusting θ. To validate the accuracy of the numerically calculated IFCs, we superimposed them with analytically calculated IFCs obtained by solving a four-layer waveguide system (air/graphene/MoO3/Au), as described in reference [58], for the first three cases (refer to Figure S5A–C). The inner elliptical and outer hyperbolic IFCs, indicated by green dashed lines, represent the fundamental and first-order HIPPhP modes. It is evident that there is excellent agreement between the results obtained from the numerical and analytical IFCs. To explore the tunability of the system, we obtained the Re(E _z ) field and corresponding IFC for a bare α-MoO₃ slab and graphene/α-MoO₃ slab with different values of E _f at θ = 20° (see Figure S6). We observed canalization of the field distributions at different doping levels. Notably, the opening angle of the hyperbola asymptote of the IFC increased with E _f , and the bare α-MoO₃ slab exhibited the smallest opening angle. This behavior aligns consistently with the field distribution patterns. The Re(E _z ) field profiles showed the emergence of a hyperboloid-shaped profile along the [100] direction and ellipse-shaped profile along the [001] direction. As E _f was increased to 0.3 eV, the elliptical distribution became dominant. The variations in the field profiles were also reflected in the IFCs, which indicates a doping-driven topological transition from the manipulation of E _f [58]. Notably, increasing E _f reduced the size of the elliptical profile in the IFCs.

Next, we investigated NR when two α-MoO₃ slabs, both covered by monolayer graphene, were interfaced with different orientations. Figure 3A depicts the left and right slabs in light green and light red, respectively. The rotation angles between the x-axis and [100] directions of the individual slabs were denoted as θ ₁ and θ ₂. Under the conditions of ω = 910 cm⁻¹, t = 150 nm, E _f = 0.15 eV, and θ ₁ = 20°, we investigated the Re(E _z ) field distributions and corresponding IFCs for θ ₂ of −10° to −90°. Figure 3B, C, and F–H depict the field distributions of Re(E _z ) with P _i and P _nt . The black dashed line indicates the interface between the two slabs. No significant power scattering or reflection occurred at the interfaces within the range of θ ₂ = −10° to −90°. Figure 3D, E, and I–K depict the corresponding IFCs where they overlapped. Here, we have superimposed the analytically calculated IFCs onto the numerical results, as illustrated in Figure 3D, E, and I–K. In the case of the left α-MoO₃ slab with a fixed orientation angle of θ ₁ = 20°, we have indicated the inner elliptical and outer hyperbolic IFCs using green dashed lines. For the right α-MoO₃ slab, we have indicated the elliptical IFCs with different angles of θ ₂ using cyan dashed lines. Note that the hyperbolic IFCs associated with the first-order HIPPhP modes of the right α-MoO₃ slab are not visible due to their significantly weaker amplitudes. The calculated IFCs for the right slab exhibited inclined orientations that varied with θ ₂. The trajectories of the incident (k _i ) and refracted (k _nt ) wavevectors are depicted with dashed green and red lines, respectively, while the incident (P _i ) and refracted (P _nt ) Poynting vectors (normal to the IFCs) are represented by solid green and red lines, respectively. The dashed white line represents the fixed y-component (k _y ) of the incident wavevector, which was preserved as light passed through the interface. Notably, P _nt was refracted toward the opposite side of the interface normal to P _i . Consequently, NR was achieved by manipulating the crystal orientations of the α-MoO₃ slab, and the NR angle can be arbitrarily determined by θ ₂. These results demonstrate the ability of the proposed structure to achieve all-angle NR in the mid-IR regime.

Figure 3:

The structure with two differently oriented α-MoO₃ slabs. (A) Graphene-covered α-MoO₃ slabs with different crystal orientations on an Au substrate, where θ ₁ and θ ₂ are the angles between the x-axis and [100] directions of the individual slabs. The left and right α-MoO₃ slabs are indicated by light green and light red colors, respectively. Re(E _z ) field distributions under the conditions of ω = 910 cm⁻¹, t = 150 nm, E _f = 0.15 eV, and θ ₁ = 20° at θ ₂ of (B) −10°, (C) −30°, (F) −50°, (G) −70°, and (H) −90° and (D), (E), (I)–(k) the corresponding IFCs, where the green and cyan dashed lines indicate the analytical calculated IFCs of the left and right α-MoO₃ slabs, respectively.

Furthermore, we examined the dependence of NR on the thickness of the α-MoO₃ slab. Figure S7 shows the Re(E _z ) fields at different values of t and corresponding IFCs under the conditions of ω = 910 cm⁻¹, E _f = 0.15 eV, θ ₁ = 20°, and θ ₂ = −30°. Please note that the outer hyperbolic IFCs of the right α-MoO₃ slab with a thickness of 200 nm (see Figure S7K) are also represented by cyan dashed lines that is the same color used for its inner elliptical IFCs. The NR angle increased moderately with increasing t. This behavior can be explained by considering the formation of the resultant HIPPhP, which arises from the coupling of a GSPP with an elliptical IFC and HIPhP with a hyperbolic IFC. The contribution of the GSPP to the resultant mode of the α-MoO₃ slab increases with decreasing t. However, decreasing t also lowers the threshold E _f required for the topological transition to occur [56], which weakens the hyperbolicity of the resultant HIPPhP at a fixed E _f . Thus, the largest NR angle was observed at t = 300 nm, as shown in Figure S7L. In addition to E _f , t provides an additional parameter for controlling the occurrence of the topological transition. Therefore, while mode canalization for any given thickness of α-MoO₃ can be adjusted through graphene doping, this adjustment is accompanied by a notable reduction in field strength along the canalization direction, while simultaneously resulting in a stronger field propagating in all directions. To study the dependence of NR on frequency, the Re(E _z ) fields for different operating frequencies within the RB II of α-MoO₃ along are shown in Figure S8A–F with the corresponding IFCs (see Figure S8G–L). Increasing ω from 830 to 930 cm⁻¹ cause the IFC patterns to transform gradually from a closed-form ellipse to an open-form hyperbola. The broadband and the all-angle NRs highlight the potential of this system for realizing NR over a broad range of frequencies and incident angles.

3.3 Effect of multiple interfaces

Figure 4A shows the extension of the proposed structure to multiple interfaces comprising three α-MoO₃ slabs with different crystal orientations. The angles between the x-axis and [100] directions of the left, upper-right, and lower-right slabs are denoted as θ ₁, θ ₂, and θ ₃, respectively. This configuration allowed the crystal orientations of each α-MoO₃ slab to be manipulated and introduced additional degrees of freedom for controlling the optical properties of the system. The angles θ ₁, θ ₂, and θ ₃ can be tuned to engineer the dispersion, wave propagation, and NR characteristics to realize various intriguing optical phenomena and functionalities. Figure 4B shows a three-dimensional view of the Re(E _z ) field distribution under the conditions ω = 910 cm⁻¹, t = 150 nm, E _f = 0.15 eV, θ ₁ = 0°, θ ₂ = −30°, and θ ₃ = 30°. The in-plane coordinate of a z-polarized electric dipole was located at (−0.5, 0), where the crossing point of the two interfaces served as the origin. As the polariton passed through the origin, it was split and refracted into two preset transmission directions corresponding to θ ₂ = −30° and θ ₃ = 30°. Figure 4C visualizes k _i and the two transmitted wavevectors k _pt and k _nt , where k _pt is the positive refraction in the IFCs derived from the three differently oriented α-MoO₃ slabs. The analytical IFCs representing the fundamental modes of the left, upper-right, and lower-right α-MoO3 slabs are denoted by green, red, and cyan, respectively. Additionally, only the hyperbolic IFCs associated with the first-order mode of the left α-MoO3 slab are observable. The corresponding P _i , P _pt , and P _nt are also indicated, where P _pt is the positive refractive Poynting vector. The three overlapping IFCs resulted from the three interfaces between the α-MoO₃ slabs. Figure 4D shows a top view of the Re(E _z ) field distribution and highlights the power flows represented by the Poynting vectors P _i , P _pt , and P _nt . In this structure, positive refraction and NR coexist, which allows it to function as a nanoscale beam splitter in photonic integrated circuits. The proposed design can serve as a functional Mach–Zehnder modulator [63] while interfacing seven α-MoO₃ slabs with various orientations. This demonstrates the versatility and applicability of the proposed structure for different types of photonic devices.

Figure 4:

The structure with three differently oriented α-MoO₃ slabs. (A) Interfacing three α-MoO₃ slabs with crystal orientations θ ₁, θ ₂, and θ ₃ between the x-axis and [100] directions of the individual slabs. (B) 3D view of the Re(E _z ) field distribution. (C) IFC with the incident wavevector (k _i ), two transmitted wavevectors (k _pt and k _nt ), and the incident (P _i ), positive refractive (P _pt ), and negative refractive (P _nt ) poynting vectors, where the green, red, and cyan dashed lines indicate the analytical calculated IFCs of the left, upper-right, and lower-right α-MoO₃ slabs, respectively. (D) Top view of the Re(E _z ) field distribution.

To further control the polaritons, we interfaced two α-MoO₃ slabs partially covered with graphene, as shown in Figure 5A. The left α-MoO₃ slab had a crystal orientation of θ ₁ = 0°, as shown in Figure 5B, and was completely covered with monolayer graphene, which was electrically gated with E _f1. The right slab had the crystal orientation θ ₂, and it was partially covered with graphene, which was electrically gated with E _f2. The separation between the graphene layers was denoted as S. A z-polarized electric dipole was positioned at a distance of d = 1 µm from the left interface on the x–y plane. The other conditions were set to ω = 910 cm⁻¹, t = 150 nm, S = 1 µm, E _f1 = 0.3 eV, θ ₁ = 0°, and E _f2 = 0.25 eV. Figure 5D and E display the Re(E _z ) field distributions for θ ₂ = 0° and −20°, respectively. As a comparison, Figure 5C illustrates the results when only the left α-MoO₃ slab was covered with graphene with both θ ₁ and θ ₂ set to 0° [58]. The |E| fields and IFCs for the three cases are also shown. Figure 5C corresponds to the structure proposed by Hu et al. [58], which exhibited NR (Figure 5I) and planar focusing (Figure 5F). However, after passing through the focal point, the field experienced conventional diffraction. In contrast, partially covering the graphene on the right α-MoO₃ slab with a suitable E _f2 = 0.25 eV overcame the diffraction of the light beam, and diffractionless propagation was achieved (Figure 5D and G). In addition, the limitation of a fixed focal point [58] was eliminated. This achievement opens up the possibility for diffractionless propagation in imaging systems and steering the diffractionless beam in wide angles to any desired direction by presetting the crystal orientation (Figure 5E and H). The IFCs for Figure 5D and E are also shown in Figure 5J and K, respectively. The parameters d, S, E _f1, and E _f2 can be adjusted to realize different imaging capabilities.

Figure 5:

The structure with two differently oriented α-MoO₃ slabs partially covered with graphene sheet. (A) Interfacing two α-MoO₃ slabs, in which the left slab is entirely covered with graphene (E _f1) but the right slab is partially covered with graphene (E _f2). (B) The same as (A), except peeling off the graphene layer to clearly show the orientations of α-MoO₃ slabs. The two graphene layers have a separation distance of S = 1 µm, where θ ₁ = 0° (left slab) and θ ₂ = 0° or −20° (right slab) indicate the angles between the x-axis and [100] directions of the individual slabs. Re(E _z ) field distributions for various cases: (C) θ ₁ = 0° and θ ₂ = 0° [56] without the graphene layer (E _f2) covering the right α-MoO₃ slab, (D) θ ₁ = 0° and θ ₂ = 0°, and (E) θ ₁ = 0° and θ ₂ = −20°. (F)–(H) Corresponding |E| fields. (I)–(K) Corresponding IFCs, where the green and cyan dashed lines indicate the analytical calculated IFCs of the left and right α-MoO₃ slabs, respectively.