Active Tuning and Anisotropic Strong Coupling of Terahertz Polaritons in Van der Waals Heterostructures

Electromagnetic field confinement is significant in enhancing light-matter interactions as well as in reducing footprints of photonic devices especially in Terahertz (THz). Polaritons offer a promising platform for the manipulation of light at the deep sub-wavelength scale. However, traditional THz polariton materials lack active tuning and anisotropic propagation simultaneously. In this paper, we design a graphene/α-MoO3 heterostructure and simulate polariton hybridization between isotropic graphene plasmon polaritons and anisotropic α-MoO3 phonon polaritons. The physical fundamentals for polariton hybridizations depend on the evanescent fields coupling originating from the constituent materials as well as the phase match condition, which can be severely affected by the α-MoO3 thickness and actively tuned by the gate voltages. Hybrid polaritons propagate with in-plane anisotropy that exhibit momentum dispersion characterized by elliptical, hyperboloidal and even flattened iso-frequency contours (IFCs) in the THz range. Our results provide a tunable and flexible anisotropic polariton platform for THz sensing, imaging, and modulation.


Introduction
In conventional optics, the ultimate spatial resolution is fundamentally restained by the diffraction limit to roughly half of the incident wavelength [1,2], which prohibits fabricating deeply sub-wavelength scale devices. Especially in the terahertz (THz) spectral range, due to the longer wavelength, the footprints of the THz devices is always about hundreds of micrometers [3,4], which is not compatible with modern integrated optics. Polaritons are hybrid modes formed by light strongly coupled to collective excitations in matter, which offers a promising platform for the manipulation of light at the deeply subwavelength scale [5][6][7][8][9]. In THz range, gold is the commonly used polariton platform for confined THz techonlogies; however, it sustains low confinement factors, high transmission losses and lack of active tunability [10,11]. Graphene, a kind of two-dimensional material, presents THz plasmon polaritons (PPs) with high confinent factors and long lifetimes. Most importantly, the Fermi level can be conveniently adjusted by electrical gating [12,13]. The only existing shortcoming in graphene plasmon is the isotropic propagation. α-MoO 3 is a kind of polar Van der Waals (vdWs) material; it is reported that the anisotropic phonon polaritons (PhPs) in mid-infrared and THz have been verified and can also be tuned by ion intercalation [14][15][16]. However, active tuning can not be executed. Therefore, a new kind of material or platform supporting THz polaritons with active tuning and anisotropic propagation is urgently needed.
VdW heterostructures provide an attractive platform for integrating bulk materials into an atomic interaction level, giving rise to an extraordinary hybridization of different polariton modes [17]. In addition, the physical fundamentals for polariton hybridizations depend on the evanescent fields coupling originating from the collective charge oscillations in constituent materials. However, the prerequisites-phase match conditions-for polariton hybridizations must be satisfied: firstly, the constituent polaritons need to propagate in the same direction with mode volumes coinciding, and secondly, the individual polaritons must possess nearly the same momenta [18][19][20]. Consequently, hybrid polaritons propagate with in-plane anisotropy that exhibits momentum dispersion characterized by elliptical, hyperboloidal and even flattened iso-frequency contours (IFCs) in the THz range.
For that purpose, we assembled graphene/α-MoO 3 heterostructures for hybridizing isotropic graphene PPs with anisotropic α-MoO 3 PhPs, thus permitting substantial in-plane anisotropy to the plasmonic resonance of graphene. It is found that the polariton hybridization can be severely affected by the α-MoO 3 thickness and actively tuned by the gate voltages. Moreover, the anisotropic momentum dispersion can be converted from elliptical IFCs to flattened bands. Finite difference time domain (FDTD) methods were used to simulate the mode field distributions and extract the polariton wavelengths. The anisotropic real-space near-field electromagnetic distributions were visualized by a phenomenological cavity model within a graphene/α-MoO 3 nanostructure, which indicates that the electromagnetic field can be restricted into specific patterns depending on the illuminating frequency. Our results reveal a tunable and scalable anisotropic polariton platform for applications in sensing, imaging, and light modulation.

Optical Conductivity and Dispersion Analysis
In THz, the optical response of α-MoO 3 is dominated by the phonon absorption, and phonon polaritons (PhPs) in α-MoO 3 can be excited within spectral intervals of the reststrahlen bands (RBs), in which the real part of the permittivity is negative. The permittivities of α-MoO 3 in THz can be described by the following Lorentz equation [15]: where j = x, y, z denotes the three principal axes of α-MoO 3 crystal, which corresponds to the crystalline directions [100], [001], and [010], respectively. ω TO,j and ω LO,j refer to the transverse-optical (TO) and longitudinal-optical (LO) phonon frequencies. Γ and ε ∞,j are the damping constant and high-frequency dielectric constant, respectively. The detailed parameters used in calculations are given in reference [15]. As shown in Figure 1a, we observe two spectral bands (marked as shaded blue (low) and yellow (high) ranges) wherein at least one of the principal components of the permittivity is negative. In particular, the low-and high-frequency RB bands are located along the [001] and [100] crystal directions, respectively, indicating their potential to support PhPs with in-plane hyperbolic propagation. Graphene's surface conductivity is modeled using the Kubo formula [21,22]: Graphene's surface conductivity is modeled using the Kubo formula [21,22]: Here, e is the electron charge, E F is the Fermi energy,h is the reduced Planck constants, and τ is the momentum relaxation time. The relaxation time of charge carriers in graphene is set to 0.5 ps, which is typical for encapsulated graphene used in optical experiments at ambient conditions. The permittivity of graphene can be written as ε = −σ/iωt g ε 0 , where t g = 0.34 nm is the thickness of graphene, and ε 0 is the vacuum permittivity. The real parts of the permittivity tensors of graphene with different Fermi levels are negative in the frequency range from 7 THz to 12 THz (Figure 1b), which promises actively tunable PPs.

Polariton Hybridization Analysis
The FDTD method has been performed to investigate the polariton response of graphene and α-MoO 3 in detail. A vertical dipole was placed above the top of the uppermost surface of 70 nm as a polariton launcher, the thickness of α-MoO 3 film was set to 20 nm and the graphene layer was regarded as a surface conductivity sheet. The IFCs The permittivity results suggest the potential formation of hybrid plasmon-phonon polariton modes in two excitation bands. We thereby construct a platform for hybridizing isotropic graphene PPs with anisotropic α-MoO3 PhPs, thus resulting in anisotropic and tunable hybridized polaritons ( Figure 3). However, the preconditions for polariton hybridization must be satisfied. In our analytical calculations and the following numerical simulations, the distance between the graphene layer and the α-MoO3 slab is zero; hence, The permittivity results suggest the potential formation of hybrid plasmon-phonon polariton modes in two excitation bands. We thereby construct a platform for hybridizing isotropic graphene PPs with anisotropic α-MoO 3 PhPs, thus resulting in anisotropic and tunable hybridized polaritons ( Figure 3). However, the preconditions for polariton hybridization must be satisfied. In our analytical calculations and the following numerical simulations, the distance between the graphene layer and the α-MoO 3 slab is zero; hence, the overlap between the graphene PPs and PhPs fields is maximal. Nevertheless, whether the PPs and PhPs couple or not is determined by whether there is a crossing point between the dispersion of PPs (without PhPs) and the dispersion of PhPs (without PPs). If there is a crossing point between the dispersions for a given configuration, there will be coupling. Therefore, the graphene Fermi level, α-MoO 3 thickness, spectra range, and polariton momenta must be optimized. The permittivity results suggest the potential formation of hybrid plasmon-phonon polariton modes in two excitation bands. We thereby construct a platform for hybridizing isotropic graphene PPs with anisotropic α-MoO3 PhPs, thus resulting in anisotropic and tunable hybridized polaritons ( Figure 3). However, the preconditions for polariton hybridization must be satisfied. In our analytical calculations and the following numerical simulations, the distance between the graphene layer and the α-MoO3 slab is zero; hence, the overlap between the graphene PPs and PhPs fields is maximal. Nevertheless, whether the PPs and PhPs couple or not is determined by whether there is a crossing point between the dispersion of PPs (without PhPs) and the dispersion of PhPs (without PPs). If there is a crossing point between the dispersions for a given configuration, there will be coupling. Therefore, the graphene Fermi level, α-MoO3 thickness, spectra range, and polariton momenta must be optimized.  The dispersion of graphene plasmons are governed by a frequency (ω)-momentum (q) relation given by [23]: where κ(ω) = ε o +ε sub (ω)

2
, ε o and ε sub are the vacuum and substrate permittivites, respectively, and σ 2D represents the sheet conductance of graphene. In calculation of the plasmon dispersion in a graphene sheet embedded between a semi-infinite air space and a semiinfinite dielectric environment of α-MoO 3 , ε sub = ε ∞,x when calculating the dispersion along the [100] crystal axis and ε sub = ε ∞,y along the [001] crystal axis. The phonons of α-MoO 3 are not considered in this calculation. The dispersion relation of electromagnetic modes in biaxial α-MoO 3 slabs can be described by the following analytical approximations [24]: where ϕ is the angle between the a axis and the in-plane wave vector q (q 2 = q 2 x + q 2 y ), ε 1 and ε 3 are the dielectric constants of the superstrate and substrate of the α-MoO 3 slab, k 0 and d are the free-space wave vector and α-MoO 3 thickness, respectively, and ε x , ε y and ε z are the permittivities of α-MoO 3 along a, b and c axes, respectively. Figure 4a,b shows the dispersion for graphene PPs with the suppressed α-MoO 3 PhPs for different Fermi energies (dashed lines), as well as the dispersions of PhPs in the α-MoO 3 slabs for different crystal thicknesses (solid lines). We can observe that the higher the Fermi energy, the wider the range of thicknesses of the biaxial slab for which there exists coupling between PPs and PhPs. If the slab is too thick, then a higher Fermi energy is required to match the PPs wavevector and achieve the coupling regime. On the contrary, when the slab is thinner, a lower Fermi energy is required. Significantly, for thinner slabs (20 nm), a wider range of graphene PPs will couple with α-MoO 3 PhPs, which permits larger polariton coupling and active tunable range. Therefore, in the following calculations, the thickness of α-MoO 3 is fixed at 20 nm.
are the permittivities of α-MoO3 along a, b and c axes, respectively. Figure 4a,b shows the dispersion for graphene PPs with the suppressed α-MoO3 PhPs for different Fermi energies (dashed lines), as well as the dispersions of PhPs in the α-MoO3 slabs for different crystal thicknesses (solid lines). We can observe that the higher the Fermi energy, the wider the range of thicknesses of the biaxial slab for which there exists coupling between PPs and PhPs. If the slab is too thick, then a higher Fermi energy is required to match the PPs wavevector and achieve the coupling regime. On the contrary, when the slab is thinner, a lower Fermi energy is required. Significantly, for thinner slabs (20 nm), a wider range of graphene PPs will couple with α-MoO3 PhPs, which permits larger polariton coupling and active tunable range. Therefore, in the following calculations, the thickness of α-MoO3 is fixed at 20 nm. To study in more detail the polariton-coupling features as a function of Fermi energy in graphene/α-MoO3 heterostructures, we calculate the Ez distribution at the crossingpoint frequencies in Figure 4. When the graphene Fermi energy is 0.1 eV, the polariton coupling should be enhanced at 11.6 THz; the near-field Ez distribution and the corresponding FFT are shown in Figure 5a. It is found that the near-field distribution To study in more detail the polariton-coupling features as a function of Fermi energy in graphene/α-MoO 3 heterostructures, we calculate the E z distribution at the crossing-point frequencies in Figure 4. When the graphene Fermi energy is 0.1 eV, the polariton coupling should be enhanced at 11.6 THz; the near-field E z distribution and the corresponding FFT are shown in Figure 5a. It is found that the near-field distribution features elliptic wavefronts with long axes centered along the crystallographic directions [100], which is significantly different from that typically observed in bare graphene or α-MoO 3 (Figure 2). By increasing the Fermi energy to 0.2 eV, the polariton coupling is enhanced at a lower frequency of 11.45 THz, which also manifests elliptic electric field wavefronts. However, the hybrid polariton wavelength significantly increased, with the colorplot of FFT centered in a narrower wavevector range. With the Fermi energy further increased to 0.3 eV, coupling frequency is shifted to 11.35 THz with polariton wavelength increased further. The evolutionary process is also the same for 0.4 eV, with the hybrid polaritons enhanced at 11.22 THz with increased wavelength as well as smaller wavevector.
For hybrid polaritons at the low-frequency band, more interesting phenomena present. When the Fermi energy is 0.1 eV, the polaritons are strongly hybridized at 10.1 THz, and the hybrid polaritons propagate in a specific direction (canalization regime). Note that the IFC exhibits a strong flattening along the [001] direction, explaining the canalization of hybrid polaritons along this direction (Figure 5e). When the Fermi energy increases to 0.2 eV, the near-field distribution features elliptic wavefronts with long axes centered along the crystallographic directions of [001], which is opposite to that for the high-frequency band. With the Fermi energy further increased to 0.3 eV and 0.4 eV, the hybrid polaritons also manifest elliptic wavefronts. However, the polariton wavelength increases with Fermi energies, and the corresponding colorplot of FFT centers in a narrower wavevector range. In addition, the hybrid polaritons manifest deeply sub-wavelength confinements of about 60 (λ 0 /λ hybridized ) with respect to the incident THz wave (λ 0 = 25.8 µm), which provides a nanoscale "hot spot" for enhancing light-matter interactions. enhanced at a lower frequency of 11.45 THz, which also manifests elliptic electric field wavefronts. However, the hybrid polariton wavelength significantly increased, with the colorplot of FFT centered in a narrower wavevector range. With the Fermi energy further increased to 0.3 eV, coupling frequency is shifted to 11.35 THz with polariton wavelength increased further. The evolutionary process is also the same for 0.4 eV, with the hybrid polaritons enhanced at 11.22 THz with increased wavelength as well as smaller wavevector. For hybrid polaritons at the low-frequency band, more interesting phenomena present. When the Fermi energy is 0.1 eV, the polaritons are strongly hybridized at 10.1 THz, and the hybrid polaritons propagate in a specific direction (canalization regime). Note that the IFC exhibits a strong flattening along the [001] direction, explaining the canalization of hybrid polaritons along this direction (Figure 5e). When the Fermi energy increases to 0.2 eV, the near-field distribution features elliptic wavefronts with long axes centered along the crystallographic directions of [001], which is opposite to that for the high-frequency band. With the Fermi energy further increased to 0.3 eV and 0.4 eV, the hybrid polaritons also manifest elliptic wavefronts. However, the polariton wavelength increases with Fermi energies, and the corresponding colorplot of FFT centers in a narrower wavevector range. In addition, the hybrid polaritons manifest deeply subwavelength confinements of about 60 (λ0/λhybridized) with respect to the incident THz wave

Manipulation of the Field Localization
The near-field simulations review found that the hybrid polaritons in graphene/α-MoO 3 heterostructure can be strongly enhanced, and more importantly, the hybrid polaritons can be actively tuned from elliptic to flattened wavefronts by electrical gating. Moreover, the precise control of light localization is very important in nanophotonics, and the electromagnetic field concentration can be further enhanced by fabricating nanostructures that act as Fabry-Perot resonators. For this purpose, we precisely shaped graphene/α-MoO 3 heterostructures into square nanostructures for tailoring electromagnetic field localizations. We adopted a phenomenological cavity model to simulate the near-field amplitude and phase patterns within graphene/α-MoO 3 nanostructures [25,26]. A point dipole source launches circular surface waves and excites hybrid polariton waves in graphene/α-MoO 3 nanostructures, the excited polariton wave propagates outward and is reflected by the nanostructure boundaries. The dipole-launched polariton waves interfere with the reflected waves perpendicular to the boundaries. The total field detected can be described as: where E 0 and E p,j are the dipole-launched polariton wave and that reflected from four boundaries, respectively. Reflected waves are given by E p,j = R j × E 0 × exp − 4π λ p r j γ p + i , where R j = R 0 exp(i∆ϕ) denotes the complex reflection coefficient, and R 0 and ∆ϕ are the reflectivity and phase shift of the polariton wave, respectively. γ p , r j and λ p represent the polariton damping rate, distance between the nanostructure edge and dipole tip, and the polariton wavelength, respectively. In simulations the amplitude of the dipole-launched polariton wave E 0 was set to 1 and the damping rate γ p to 0.4 eV, whereas the polariton wavelength λ p is set based on Figure 5. Reflectivity was assumed to be 1, and 1.5π was used for the phase shift of the polariton wave [26]. Figure 6 shows the simulated near-field patterns of the hybrid polaritons in graphene/α-MoO 3 micro-flakes at the cross points in Figure 4. Due to the multi-beam interferences caused by additional reflections from the square boundaries, near-field patterns of the hybrid polaritons give rise to a variety of bright and dark features as a function of exciting frequencies and Fermi energies. It is obvious that the bright (dark) spots resulting from constructive (destructive) interferences exhibit a D4-group symmetry. Both the amplitude and phase images manifest standing wave patterns, whereas the adjacent bright and dark fringes close to the boundaries differ by~π/2. In addition, due to the damping of the reflected plasmon waves, more fringes with much weaker intensities can be observed in the interior part with the adjacent fringe separation of 0.5 λ p . In addition, due to the elliptical-type dispersion, hybrid polaritons propagate along x and y axes at different wavelengths. Except for the flattened-band dispersion at 10.1 THz, polariton propagation along y axis is forbidden, which results in interference fringes only appearing in the horizonal direction. As the frequency decreases, the number of bright spots decreases and the field localizations are strongest at the regions close to the corner of the nanostructures. The realization of electromagnetic field localization in two dimensions using graphene/α-MoO3 nanostructures can offer several advantages in terms of efficient focusing and manipulating of light at the nanoscale, which can further be integrated with photonic or optoelectronic devices to enhance light-matter interactions. In addition, due to the elliptical-type dispersion, hybrid polaritons propagate along x and y axes at different wavelengths. Except for the flattened-band dispersion at 10.1 THz, polariton propagation along y axis is forbidden, which results in interference fringes only appearing in the horizonal direction. As the frequency decreases, the number of bright spots decreases and the field localizations are strongest at the regions close to the corner of the nanostructures. The realization of electromagnetic field localization in two dimensions using graphene/α-MoO 3 nanostructures can offer several advantages in terms of efficient focusing and manipulating of light at the nanoscale, which can further be integrated with photonic or optoelectronic devices to enhance light-matter interactions.

Conclusions
In conclusion, we designed and investigated polariton hybridizations in graphene/α-MoO 3 heterostructures. It was found that the polariton hybridization between isotropic graphene PPs and anisotropic α-MoO 3 PhPs strongly depends on the α-MoO 3 thickness and can be actively tuned by the gate voltages. In addition, the hybrid polaritons propagate with in-plane anisotropy that exhibits momentum dispersion characterized by elliptical, hyperboloidal and even flattened IFCs in the THz range. The near-field results suggest that the localized electromagnetic fields can be enhanced in nanostructures, which provide a flexible anisotropic polariton platform and can further be integrated with photonic crystals or quantum well structures to ulteriorly enhance THz-matter interactions. Our findings demonstrate that active tuning of THz polaritons with certain propagations can be used to steer THz light along specific directions at the nanoscale. More significantly, these results will provide guidelines for THz device fabrications and experiments with dynamically controllable properties.