Polaritons, coupled modes arising from the strong interactions between photons and material excitations, provide a unique way to harness and manipulate light at the nanoscale^1-6. Conventional surface plasmon polaritons (SPPs) propagating along a metal-dielectric interface have been extensively studied in terms of intrinsic electronic scattering and plasmonic loss. The low loss and ultrahigh confinement of emergent phonon polaritons (PhPs) in van der Waals (vdW) materials (e.g., hBN, α-MoO₃, V₂O₅)^7-12 show considerable promise for integrated infrared nanophotonics. In particular, polaritons in so-called hyperbolic media, namely, hyperbolic phonon polaritons (HPhPs), exhibit hyperbolic dispersion, that is, their iso-frequency surface in momentum space is hyperboloidal-like^{13, 14}. Consequently, HPhPs allow an enormous density of optical states at specific wavevectors and ray-like directional polariton propagation, showing remarkable potential for application in subdiffractional focusing^{15, 16}, thermal management¹⁷, in-plane canalization^18-20 and hyper-lensing^{21, 22}.

Dynamical tunability of PhPs is required for integrated nanophotonics, including adjustable propagation characteristics and reconfigurable dispersion. The hyperbolic regime of PhPs is fundamentally controlled by crystal symmetry, making tailoring the polariton topology more challenging. Recently, continuous efforts have been devoted to tuning the PhPs in anisotropic vdW materials^23-40. Artificial nanostructures, such as grating⁴¹, slot⁴², and nanohole²⁰ show the ability in tuning the polaritons. However, constructing these metamaterials or metasurfaces faces difficulties in processing and manufacturing. In comparison, heterostructuring and hybridizing modes with other polaritons can modify the properties of PhPs through simple structures, attracting extensive interest^23-27. Due to the coupling effect of evanescent waves, the PhPs in the twisted bilayer system can undergo a topological transition, i.e., the polariton topology evolves from a hyperbolic shape to an elliptical one, realizing the field canalization and dispersion flatting^28-33. The above measure mainly adopts two or multi-layer structures, restricting its application in bare films. Another approach to manipulating PhPs is to construct the graphene/polar vdW heterostructures^34-36, in which the hybrid modes with both phonon and plasmon characteristics can be actively modified by changing the chemical potential through the gating or doping method. This idea has been widely explored in the dynamic control of polaritons involving hBN and α-MoO₃. Moreover, considering that polaritons are sensitive to the local environment, the effect of substrate plays a significant role in tailoring the polaritons^{37-40, 43}. Image polaritons or acoustic polaritons^{44, 45} are investigated when hBN or α-MoO₃ are placed in proximity to a metal plate with a thin dielectric spacer. However, previous studies mainly focus on materials with large negative permittivity, like Au or positive permittivity like SiO₂, since they are commonly used in practice. The effect of substrate media with moderate negative permittivity on PhPs has so far remained elusive. In addition, a systematic analysis about the effect of the superstrate and substrate has not been explored. Thus, it is critical to quantify how the local dielectric environment of the substrate modifies propagating PhPs.

The natural anisotropic vdW materials provide an ideal platform for exploiting the effect of in-plane anisotropy in the application of planar optical devices. Especially, in the hyperbolic regime, the theoretically infinite wavevectors and extremely large optical density of states open the door to enhance and modulate the thermal radiation and heat transfer. In addition to tuning PhPs through substrate effect, it is equally important to investigate the effect of practical substrate media on the near-field radiative heat transfer (NFRHT) when considering experimental implementations. However, the fundamental studies on the substrate effect on electromagnetic energy transfer are still scarce, hindering its further exploration and applications.

Here, we present a comprehensive study of the effect of substrate with negative permittivity on the hyperbolic phonon polaritons. We theoretically demonstrate that a reorientation of the iso-frequency contour (IFC) of hyperbolic polaritons in a biaxial slab can be engineered to allow propagation along forbidden directions when the slab is laid on a substrate with a moderate negative permittivity. In addition, we reveal that the substrate with negative permittivity can generate the elliptical-like fundamental mode which is suppressed in the configuration with a substrate having positive permittivity. In particular, we demonstrate that the direction of propagation of in-plane HPhPs in an α-MoO₃ slab can be rotated by 90°, namely, along the intrinsically forbidden direction in a bare slab by placing the slab on a SiC substrate. In the heterostructure consisting of α-MoO₃ and Au, the fundamental mode with elliptical IFC is excited. Besides, we systematically study the effect of the polaritonic coupling between heterostructures integrated with different substrates on radiative heat transfer. The effect of the excitation, reorientation, and annihilation of the fundamental mode on photon tunneling is analyzed via dispersion contour and photon transmission coefficient (PTC). The enhancive or suppressive effect depends on the relative magnitude of gap size and slab thickness. Depending on the substrates, the coupling effect can modulate the radiative heat transfer and photon tunneling. Finally, the effect of the crystalline direction on NFRHT between heterostructures integrated with the substrate of SiC and Au is also considered and compared with the suspended configuration.

To elaborate on the theoretical conditions upon which this reorientation of the directional propagation of HPhPs can take place, let us start analyzing the general case of an α-MoO₃ slab embedded between two semi-infinite isotropic media with dielectric functions ε₁ and ε₃ as shown in {L-End} Fig. 1(a). The α-MoO_{Fig. 3} slab has been recently revealed as a biaxial vdW crystal with natural hyperbolicity in the mid-infrared region. The optical response of the α-MoO_{Fig. 3} is dominated by the phonon absorption and its permittivity tensor components can be described by a Lorentz model:

1${ϵ}_j={ϵ}_{\infty }^j\left(1+\frac{{\omega }_{L\mathrm{,}j}^2-{\omega }_{T\mathrm{,}j}^2}{{\omega }_{T\mathrm{,}j}^2-{\omega }^2-\text{i}\omega {\Gamma }_j}\right)\mathrm{,}$

where j = x, y, z, and detailed parameters can be found in ref^{9, 46}. Here, the principal axis of α-MoO₃ crystal coincides with the Cartesian coordinate axis. Without rotation, the crystal axes [100], [001], and [010] correspond to the x, y, and z axes in the Cartesian coordinate system, respectively. α-MoO₃ supports phonon polaritons (PhPs) in three Reststrahlen bands (RBs): RB 1 in the range from 545 to 851 cm⁻¹ and RB 2 in the range from 820 to 972 cm⁻¹, in which the in-plane iso-frequency contour is hyperbolic; while an in-plane dispersion of elliptical shape exists for PhPs in RB 3 (from 958 to 1010 cm⁻¹), as shown in {L-End} Fig. 1(b). Without loss of generality, we assume the representative dielectric functions in RB 2, where Re(ε_x) < 0, Re(ε_y) > 0, and Re(ε_z) > 0. The dispersion relation of polaritons in the biaxial slab in the high-momenta approximation reads⁴⁷

2$\begin{array}{ll}k\left(\omega \right)=\hfill & \frac{\rho }{d}\left[\arctan \left(\frac{{ϵ}_1\rho }{{ϵ}_z}\right)+\arctan \left(\frac{{ϵ}_3\rho }{{ϵ}_z}\right)+\pi l\right]\mathrm{,}\hfill \\ \hfill & l=\mathrm{0,1,2,...,}\hfill \end{array}$

where $k=\sqrt{k_x^2+k_y^2}$ is the modulus of the in-plane wavevector component k, $\rho =\text{i}\sqrt{{ϵ}_z/\left({ϵ}_x{\cos }^2\phi +{ϵ}_y{\sin }^2\phi \right)}$ with φ being the angle between the x axis and k, as shown in {L-End} Fig. 1(a), d represents the thickness of the biaxial slab and l denotes the order of different polaritonic modes, where l = 0 represents the fundamental mode, and l = 1, 2, Fig. 3… denotes higher-order mode. Since arctan(x)+ arctan(y)= arctan ((x + y) / (1 − xy)), Eq. (2) can be conveniently rewritten as:

3$\tan \left(\frac{kd}{\rho }\right)=\frac{\frac{\rho \left({ϵ}_1+{ϵ}_3\right)}{{ϵ}_z}}{1-\frac{{\rho }^2{ϵ}_1{ϵ}_3}{{ϵ}_z^2}}\mathrm{.}$

The propagation direction and the shape of the iso-frequency contour are not trivial for the general case of arbitrary complex dielectric functions. Therefore, we assume that the permittivities of the biaxial slab are purely real. This assumption is reasonable under a small loss.

In this section, we consider only the fundamental mode, i.e., l =0. Considering that the dielectric function of the substrate is negative, i.e., ε₃ < 0, and that of the superstrate is positive, i.e., ε₁ > 0, we derive the conditions under which PhPs propagate along the x or y axis in RB 2 where ε_x < 0, ε_y >0, and ε_z > 0, according to Eq. (3). Firstly, let us study the propagation of polaritons along the x axis. In this case, we set φ = 0 and therefore, $\rho =\text{i}\sqrt{{ϵ}_z/{ϵ}_x}$. Since ${ϵ}_x$ < 0 in RB 2, then $\text{i}\sqrt{{ϵ}_z/{ϵ}_x}=-\sqrt{{ϵ}_z/\left|{ϵ}_x\right|}$, so Eq. (3) can be simplified as:

4$\tan \left(-\frac{kd}{\sqrt{{ϵ}_z/\left|{ϵ}_x\right|}}\right)=\frac{-\sqrt{\frac{{ϵ}_z}{\left|{ϵ}_x\right|}}\frac{\left({ϵ}_1+{ϵ}_3\right)}{{ϵ}_z}}{1-\frac{{ϵ}_1{ϵ}_3}{\left|{ϵ}_x\right|{ϵ}_z}}\mathrm{,}$

here $0<\frac{kd}{\sqrt{{ϵ}_z/\left|{ϵ}_x\right|}}<\frac{\pi }{2}$ due to the considered fundamental mode in this section. Consequently, the left-hand side of Eq. (4) is negative. Therefore, for Eq. (4) to have a solution, the right-hand side must also be negative. Note that ε₁ > 0 and ε₃ < 0, thus, for Eq. (4) to have a solution, the numerator on the right-hand side of Eq. (4) must be negative. Hence, we get the necessary conditions for the propagation of polaritons along x axis:

5${ϵ}_1+{ϵ}_3>0.$

Similarly, let us study the propagation of PhPs along the y axis. In this case, we set $\phi =\pi /2$ so that $\rho =\text{i}\sqrt{{ϵ}_z/{ϵ}_y}$. Since ${ϵ}_y$ > 0 in RB 2, $\sqrt{{ϵ}_z/{ϵ}_y}$ is positive and pure real. For the convenience of derivation, we set $\xi =\sqrt{{ϵ}_z/{ϵ}_y}$ ( with ξ also positive and real). Thus Eq. (3) can be simplified as:

6$\tan \left(-\text{i}\frac{kd}{\xi }\right)=\frac{\text{i}\xi \frac{\left({ϵ}_1+{ϵ}_3\right)}{{ϵ}_z}}{1+\frac{{ϵ}_1{ϵ}_3}{{ϵ}_z{ϵ}_y}}\mathrm{.}$

Taking into account that tan(·) is an odd function and tan (iη) = i tanh(η), Eq. (6) can be rewritten as:

7$\tanh \left(\frac{kd}{\xi }\right)=\frac{-\xi \frac{\left({ϵ}_1+{ϵ}_3\right)}{{ϵ}_z}}{1+\frac{{ϵ}_1{ϵ}_3}{{ϵ}_z{ϵ}_y}}\mathrm{.}$

Note that the left-hand side of Eq. (7) is positive, so the right-hand side must also be positive. There are two cases, one is that if ${ϵ}_y{ϵ}_z+{ϵ}_1{ϵ}_3<0$, which holds for a high-permittivity superstrate, the necessary condition for the polaritons to propagate along the y axis reads:

8${ϵ}_1+{ϵ}_3>0.$

The other is that if ${ϵ}_y{ϵ}_z+{ϵ}_1{ϵ}_3>0$ (the superstrate is a low-permittivity medium, such as air), the necessary condition becomes

9${ϵ}_1+{ϵ}_3<0.$

Our analysis above shows that the HPhPs can propagate along the x axis when ${ϵ}_1+{ϵ}_3>0$. But for PhPs to propagate along the y axis, stricter conditions should be satisfied, as listed in Table 1. Note that in the case when ${ϵ}_1+{ϵ}_3<0$ and ${ϵ}_y{ϵ}_z+{ϵ}_1{ϵ}_3<0$ are simultaneously satisfied, the l = 0 mode is suppressed (i.e., annihilated).

Therefore, we demonstrate that the iso-frequency contour (IFC) of in-plane hyperbolic polaritons can realize a 90° reorientation in a biaxial slab in the neighborhood of Re(ε₁) + Re(ε₃) = 0. Namely, the HPhPs can be engineered to allow propagation along intrinsically forbidden directions by placing the slab on a substrate with a given negative permittivity. For other RBs, the fundamental mode can also be reoriented (see details in the Supplementary information). Note that for RB 3, the l = 0 mode cannot be excited when Re(ε₁) + Re(ε₃) > 0, while it can be excited when Re(ε₁) + Re(ε₃) < 0.

We first demonstrate the reorientation of PhPs in α-MoO₃ /SiC heterostructure through the dispersion properties and electric field distributions and distinguish the type of hyperbolic polaritons based on the real and imaginary parts of the out-of-plane wavevector. Subsequently, we investigate the effect of excitation, annihilation, and reorientation of polaritons on the near-field energy transports between α-MoO₃ crystal with SiC or Au substrate in terms of the effect of vacuum gap width, the effect of thickness of α-MoO₃ slab as well as the effect of the crystalline directions of α-MoO_3.

To address the prediction that the propagating direction of the fundamental mode can be reoriented by changing the permittivity of the substrate, we place the biaxial vdW α-MoO₃ slab on top of the polar crystal SiC, whose relative permittivity can be described by a Lorentz model^{48, 49} (see details in the Supplementary information). Here, air is considered as the superstrate (ε₁ = 1). The permittivity of SiC (shown in Fig. S1) attains a value Re(ε₃) = −1 at ω₀ = 948.6 cm⁻¹ where it approximately fulfills the critical condition Re(ε₁ ) + Re(ε₃) = 0. To explore the possibility, we first investigate the propagation direction of the l = 0 mode of the α-MoO₃ slab placed on top of SiC and compare it with the case of air substrate (ε₃ = 1) where the condition Re(ε₁) + Re (ε₃) > 0 is always fulfilled in the studied spectral range. Note that only if the tangential wavevector is positive, can the PhPs be excited. Therefore, to indicate the excitation of PhPs, we calculate the minimum tangential wavevector of the excited l = 0 mode of the α-MoO₃ slab placed on top of the SiC substrate. Taking RB 2 as an example, we get the minimum wavevector of the excited l = 0 mode along y-axis direction by further simplifying Eq. (7):

10$k=\frac{\xi }{d}\arctan \left(\frac{-\xi \frac{\left({ϵ}_1+{ϵ}_3\right)}{{ϵ}_z}}{1+\frac{{ϵ}_1{ϵ}_3}{{ϵ}_y{ϵ}_z}}\right)\mathrm{.}$

In RB 2, the in-plane HPhPs cannot propagate along y-axis direction for the α-MoO₃ slab embedded in air, while its direction of propagation can be rotated by 90° when the α-MoO₃ slab is placed on the SiC substrate, thus occurring along intrinsically forbidden directions. Namely, k calculated using Eq. (10) takes positive values. Note that the arctanh function parameters must be between −1 and 1 (excluding −1 and 1) due to their pure real number. Therefore, we get the spectral region (labeled as the blue shaded region) where the l = 0 mode propagates along the [001] direction with SiC substrate in {L-End} Fig. 2(a). For other regions like the RB 1, the l = 0 mode propagates along the [001] direction with air substrate while along the [100] direction (forbidden direction) with the SiC substrate as shown in the blue shaded region of {L-End} Fig. 2(b). Another configuration is that the substrate is replaced with Au, described with Drude and Lorentz models^{Fig. 50} (see details in the Supplementary information), which exhibits largely negative permittivity in the considered spectral region. As shown in Fig. S2, in the range from 958.59 cm⁻¹ to Fig. 1004.68 cm⁻¹, the minimum wavevector along the [001] direction is positive, yielding a 90° rotation of the propagation direction of HPhPs in the transition region between RB 2 and RB Fig. 3. In addition, the positive value of the wavevector in RB Fig. 3 exhibits the excitation of l = 0 mode in the α-MoO_{Fig. 3}/Au heterostructure (See details in the Supplementary Information). However, for the other spectral regions, the excitation of l = 0 mode is suppressed.

From Eq. (10), the wavevector of l = 0 mode evidently depends on the permittivity of the substrate. In {L-End} Fig. 3(a), the minimum k along the [001] crystalline direction is positive when the exciting frequency and permittivity of the substrate are located in the regions labeled 1 and 2, which are surrounded by black and blue dashed lines. Although the k value is positive in region Fig. 3 where the ε₂ is more negative, i.e., the absolute value of ε₂ is larger, they do not satisfy that the arctanh function parameters must be between −1 and 1 (excluding −1 and 1) in Eq. (10), thus resulting in no excitation of l = 0 mode. Namely, as the ε_{Fig. 3} decreases, the minimum exciting frequency of l = 0 mode increases until it goes beyond RB 2. This explains that in the case of largely negative permittivity of substrate such as Au, the l =0 mode can not occur in RB 2. Moreover, it is found that the imaginary part of substrate permittivity, i.e., ${{ϵ}^{\prime }}^{\prime }$, also affects the wavevector of the reoriented polaritons as shown in {L-End} Fig. 3(b). Here we fix the real part of substrate permittivity i.e., ${ϵ}^{\prime }$ to −2, and change the ${{ϵ}^{\prime }}^{\prime }$ from 10⁻⁵ to 1. The results show that in RB 2, the k is significantly reduced as ${{ϵ}^{\prime }}^{\prime }$ is increased to 1.

To address the theoretical prediction above, we investigate the analytical dispersion of HPhPs in the 100-nm-thick α-MoO₃ slab when it is placed on the top of the SiC and Au substrates, respectively. Meanwhile, the results are compared with the air substrate configuration. {L-End} Fig. 4(a) shows the dispersion of propagating HPhPs along both in-plane crystalline directions ([100] and [001], corresponding to the x and y axes) in the air substrate configuration. The air-embedded α-MoO_{Fig. 3} slab supports the l =0 mode along the [100] direction, while there are no polaritonic modes along the [001] direction in the RB 2, agreeing with our analysis above. In the RB Fig. 3, the lowest order of polaritonic mode is l = 1, hence it suppresses the excitation of the l = 0 mode. Replacing air with SiC, the propagating direction of the l = 0 mode is observed along the [100] direction between the LO phonon frequency of α-MoO_{Fig. 3}, ω_LO = 96Fig. 3 cm⁻¹ and phonon frequency of SiC, ω₀ = 948.6 cm⁻¹ at which Re (${ϵ}_{\mathrm{Fig.\; 3}}$) = −1. On the other hand, the l = 0 mode emerges along the [001] direction, namely, the forbidden direction below ω₀, i.e., when Re (${ϵ}_{\mathrm{Fig.\; 3}}$) < −1 in {L-End} Fig. 4(b). For the case of Au substrate, the lowest mode in RB 1 and RB 2 is l =1 mode due to the influence of the generated image charge^{Fig. 51} while in RB Fig. 3 it is l =0 mode due to largely negative permittivity of Au. Compared with air or SiC substrates, Au can excite the l =0 mode of the α-MoO_{Fig. 3} slab in RB Fig. 3 in {L-End} Fig. 4(c).

To further study the HPhPs for different substrates, the IFC of the HPhPs dispersion relation is plotted in the k_x-k_y momentum space and compared to the false color plot of Im (r_p (k_x, k_y)). Exhibiting a hyperbola-shaped curve, the IFC of air substrate case gives rise to the propagation of hyperbolic polaritons within hyperbolic sectors between asymptotes, which are defined as $k_y=\pm \sqrt{\left|{{{ϵ}^{\prime }}^{\prime }}_{xx}/{{{ϵ}^{\prime }}^{\prime }}_{yy}\right|}{k}_x$, centered along the [100] crystalline direction ({L-End} Fig. 5(a)) at an incident frequency, ω = 9Fig. 36 cm⁻¹, while it leads to a reorientation of the propagation direction of hyperbolic polaritons, rotating from the [100] direction to the [001] direction in the α-MoO_{Fig. 3}/SiC heterostructures ({L-End} Fig. 5(b)). Then, at another incident frequency, ω = 989 cm⁻¹ (in RB Fig. 3), the shape of the IFC of the in-plane polaritonic mode is elliptical as shown in {L-End} Fig. 5(c). However, the lowest order polariton changes from l =1 mode to l = 0 mode when the substrate is placed by Au, as shown in {L-End} Fig. 5(d). To further validate these findings, the full-wave numerical simulations are performed, below which a vertically oriented electric dipole, acting as the polaritonic source, is placed on the top of the heterostructures. For the individual α-MoO_{Fig. 3} slab ({L-End} Fig. 5(e)), HPhPs are propagating with concave wavefronts within hyperbolic sectors centered along the [100] directions. However, in α-MoO_{Fig. 3}/SiC system ({L-End} Fig. 5(f)), we observe HPhPs propagating with concave wavefronts within hyperbolic sectors centered along the [001] direction, which is the intrinsically forbidden polaritonic direction in a single slab of α-MoO_{Fig. 3}. The results of the corresponding Fourier transform also reveal a 90° rotation of the major axis of the polaritonic IFC. These demonstrate that the directional propagation of in-plane HPhPs in α-MoO_{Fig. 3} slab can be effectively engineered along the previously forbidden direction. So far, there are two types of hyperbolic polaritons: volume-confined hyperbolic polaritons (v-HPs) and surface-confined hyperbolic polaritons (s-HPs). The v-HPs propagate directionally inside the crystal with a purely real-value out-of-plane wavevector^{Fig. 52}. By contrast, the s-HPs decay exponentially in the direction perpendicular to the interface^{Fig. 5Fig. 3}. Hence, the value of the out-of-plane wavevector is purely imaginary.

To study reoriented hyperbolic polaritons in the α-MoO₃/SiC and α-MoO₃/Au heterostructures, we first visualize the cross-sectional polariton field distribution. The electric-field map corroborates the existence of the v-HPs ({L-End} Fig. 6(a) and {L-End} Fig. 6(c)), consistent with what has been reported in prior works^{Fig. 5}. For the electric field in {L-End} Fig. 6(b), corresponding to l = 0 mode in the α-MoO_{Fig. 3}/SiC heterostructure, we cannot directly distinguish the type of hyperbolic polaritons due to the lack of exponentially decaying features away from the air/α-MoO_{Fig. 3} interface (z > 0). Therefore, we pay attention to the out-of-plane wavevector, q_ez^{47, Fig. 54}. {L-End} Fig. 7 shows the real and imaginary parts of q_ez, for incident frequency ω = 9Fig. 36 cm⁻¹ and 989 cm⁻¹. It is found that for Point A ({L-End} Fig. 7(a)), the q_ez is the real part dominated and its imaginary part is negligibly small, consistent with the feature of v-HPs. For Point B ({L-End} Fig. 7(b)), the q_ez is imaginary part dominated, yielding the decay away from the interface. Thus, HPhPs along the [001] α-MoO_{Fig. 3} crystalline direction show surface-confined propagation. By contrast, the q_ez with ω = 989 cm⁻¹ ({L-End} Fig. 7(c) and {L-End} Fig. 7(d)) is always the real part dominated in the whole k_x-k_y momentum space, ensuring the propagating feature inside the α-MoO_{Fig. 3} crystal. Hence, the l = 0 mode in the α-MoO_{Fig. 3}/Au heterostructures is also v-HPs.

The substrate media have an obvious influence on the HPhPs excited in the α-MoO₃ crystal. In contrast to far-field thermal radiation, evanescent waves coupling between two bodies with subwavelength separation distance can assist photons to tunnel through the vacuum gap (i.e., photon tunneling), yielding a radiative heat flux that can exceed the blackbody limit by several orders of magnitude^55-61. NFRHT has inspired a lot of research interest for its fundamental scientific relevance and technological application, including thermal rectification, thermal transistors, thermophotovoltaics, and photonic cooling^62-64, etc. To benefit these potential applications, continuous efforts have been devoted to exploring polaritonic modes that could intensify photon tunneling and control thermal radiation. Recently, researchers have found that the hybridization of different kinds of polaritons provides an exciting paradigm to modulate and control polaritons^{65, 66}. Inspired by this concept, many interesting ideas have been put forward in recent years. To date, most proposals for hybridization of polaritons rely strongly on suspended structures, ignoring the possible effects of substrates and not considering experimental implementations. Besides, the effect of the reorientation of hyperbolic polaritons in vdW crystal with substrates remains unexplored. To address this issue, we investigate the influence of polaritonic reorientations induced by substrates on NFRHT between heterostructures by applying the fluctuation-dissipation theorem.

We consider the NFRHT between two aligned heterostructures with the temperatures T₁₍₂₎ = 300 (299) K and vacuum gap d. Each heterostructure consists of a thin α-MoO₃ slab with a thickness of t, and a substrate with infinite thickness as shown in {L-End} Fig. 8(a). In this section, we focus on the analysis of the radiative heat flux (RHF), Q, which is given by^{67, 68}:

11$\begin{array}{ll}Q=\hfill & {\displaystyle {\int }_0^{\infty }}q\left(\omega \right)\text{d}\omega =\frac{1}{8{\pi }^3}{\displaystyle {\int }_0^{\infty }}\left[\Theta \left(\omega \mathrm{,}{T}_1\right)-\Theta \left(\omega \mathrm{,}{T}_2\right)\right]\text{d}\omega \hfill \\ \hfill & \times {\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}\xi \left(\omega \mathrm{,}{k}_x\mathrm{,}{k}_y\right)\text{d}{k}_x\text{d}{k}_y\mathrm{,}\hfill \end{array}$

where $\Theta \left(\omega \mathrm{,}T\right)=\hslash \omega /\left({\text{e}}^{\hslash \omega /k_{\text{B}}T}-1\right)$ is the mean energy of a Planck oscillator at angular frequency ω and temperature T. $\xi \left(\omega \mathrm{,}{k}_x\mathrm{,}{k}_y\right)$ is the PTC, which is expressed as^{65, 66}

12$\xi \left(\omega \mathrm{,}{k}_x\mathrm{,}{k}_y\right)=\{\begin{array}{c}\text{Tr}\left[\left(I-R_2{}^{*}{R}_2-T_2{}^{*}{T}_2\right)D\left(I-R_1{}^{*}{R}_1-T_1{}^{*}{T}_1\right)D^{*}\right]\mathrm{,}k<k_0\\ \text{Tr}\left[\left(R_2{}^{*}-R_2\right)D\left(R_1-R_1{}^{*}\right)D^{*}\right]{\text{e}}^{-2\left|k_{z0}\right|d}\mathrm{,}k>k_0\end{array}\mathrm{,}$

where $k_{z0}=\sqrt{k_0^2-k^2}$ is the perpendicular wavevector component in vacuum and * denotes the complex conjugate. R and T are matrixes of Fresnel reflection and transmission coefficient for a plane wave incidence from vacuum to medium 1 or 2, respectively, and have the following forms:

13$R_{\mathrm{1,2}}=\left[\begin{array}{cc}{r}_{\text{ss}}^{\mathrm{1,2}}& r_{\text{sp}}^{\mathrm{1,2}}\\ r_{\text{ps}}^{\mathrm{1,2}}& r_{\text{pp}}^{\mathrm{1,2}}\end{array}\right]\mathrm{,}{T}_{\mathrm{1,2}}=\left[\begin{array}{cc}{t}_{\text{ss}}^{\mathrm{1,2}}& t_{\text{sp}}^{\mathrm{1,2}}\\ t_{\text{ps}}^{\mathrm{1,2}}& t_{\text{pp}}^{\mathrm{1,2}}\end{array}\right]\mathrm{.}$

The 2$\times$2 matrix D describes the Fabry-Perot-like denominator, defined as D = (I−R₁R₂${\text{e}}^{2\text{i}{k}_{z0}d}$)⁻¹.

To reveal the effect of the reorientation of HPhPs on NFRHT, we adopt two kinds of practical substrates: SiC and Au and compare them with the individual slab configuration. {L-End} Fig. 8(b) shows the radiative heat flux as a function of the vacuum gap width, d in three kinds of structures: (I) bare α-MoO_{Fig. 3} slab, (II) α-MoO_{Fig. 3} / SiC heterostructure, (III) α-MoO_{Fig. 3} / Au heterostructure. The thickness of the α-MoO_{Fig. 3} slab is fixed to 20 nm. As a benchmark, {L-End} Fig. 8(b) shows the increase of near-field RHF increases by several orders of magnitude over the blackbody limit due to the evanescent contribution. The blackbody limit at room temperature can be given by $h={\sigma }_{\text{B}}{T}^4=6.09{\text{Wm}}^{-2}$, where σ_B is Stefan-Boltzmann constant (gray dashed line in {L-End} Fig. 8(b)). In particular, the results show that the total heat flux of each configuration is enhanced up to 2000-fold above the blackbody limit at d = 10 nm. To directly exhibit the enhancive and annihilative effect of the substrates, we define the relative enhancement coefficient η to describe the ratio of heat fluxes for the case with the substrate to that for suspended α-MoO_{Fig. 3} slabs as:

14$\eta =\frac{Q_{\text{sub}}}{Q_{\text{air}}}\mathrm{,}$

where Q_sub and Q_air refer to the radiative heat flux with and without substrates, respectively.

Further, we observe that the results for SiC and Au substrates exhibit different features when the vacuum gap increases. For the configuration of α-MoO₃ / SiC heterostructure, when d is less than 25 nm (Region I), the radiative energy exchange is basically the same as that of bare α-MoO₃ slab, independent on the type of substrate. This is because the films in this case are in effect bulk materials when the nanoscale gap width becomes comparable to or smaller than the film thickness⁶⁹. In contrast, even when the vacuum gap is smaller than the thickness of the α-MoO₃ film in the α-MoO₃ / Au heterostructure, the α-MoO₃ film can still not be regarded as a bulk material. It can be seen from {L-End} Fig. 8(c) that the Au substrate always suppresses the radiative heat transfer and the suppressive effect strengthens as the vacuum gap width increases. SiC substrate suppresses the radiative energy transport in Region II (25 nm < d < 210 nm), while the effect of SiC substrate on NFRHT transforms from suppression to enhancement as d exceeds 210 nm (Region III). In particular, at d = Fig. 1000 nm, the total heat flux of α-MoO_{Fig. 3} / SiC heterostructure system yields 14.7 W/m², which is 8 times larger than that of α-MoO_{Fig. 3} / air system.

To reveal the physical origin of the enhancement effect at a large vacuum gap width, we calculate the spectral heat flux for three kinds of structures at d =1000 nm as shown in {L-End} Fig. 9(a). Note that at the nanoscale, the main contribution to NFRHT comes from the excitation of surface-confined or volume-confined modes inside the RBs. However, at larger gap width, i.e., far-field region, it is the propagating wave that dominates. As shown in {L-End} Fig. 9(a), there are three spectral peaks for bare α-MoO_{Fig. 3} slab and one spectral peak for α-MoO_{Fig. 3} / Au system and these peaks are all located inside the RBs. The phonon tunneling between two α-MoO_{Fig. 3} heterostructures (i.e., the PTC) can explain the origin of these modes. {L-End} Fig. 10(a) and {L-End} Fig. 10(b) represent the PTC of the first and second peaks of α-MoO_{Fig. 3} / SiC system. One can see that the hyperbolic features come from the coupling of the epsilon near pole modes and hyperbolic volume-confined modes. As shown in {L-End} Fig. 10(c), the PTC at the spectral peak of α-MoO_{Fig. 3} / Au system exhibits elliptical topology, indicating that it results from the coupling the elliptical volume-confined modes. Therefore, combining with {L-End} Fig. 9(b) and {L-End} Fig. 10(a–c), one can see that the contribution to the NFRHT between two heterostructures mainly comes from the inside of RBs for α-MoO_{Fig. 3} / air and α-MoO_{Fig. 3} / Au systems. However, it becomes complicated for α-MoO_{Fig. 3} / SiC system. It is outside of RBs that dominates in energy exchange between two α-MoO_{Fig. 3}/SiC heterostructures at large vacuum gap width, i.e., Fig. 1000 nm. As shown in {L-End} Fig. 10(d) and {L-End} Fig. 10(e), the evanescent wave (k > k₀) and propagating wave (k < k₀) outside of RBs both make contributions. Particularly, in RB 2, the PTC of {L-End} Fig. 10(f) exhibits the closed features, suggesting that the phonon tunneling is mainly affected by the modes excited by SiC substrate. Therefore, we can get a conclusion that the SiC substrate amplifies the enhancive effect on NFRHT compared with the configuration without substrate at a large vacuum gap.

In this section, we fix the vacuum gap d = 20 nm. {L-End} Fig. 11(a) shows the RHF as a function of thickness of α-MoO_{Fig. 3} slab, t. It is obviously seen that in the case of large thickness, the effect of the substrates, i.e., SiC and Au, on radiative heat transfer is negligible, since the slabs can be regarded as bulk materials. As a comparison, the critical thickness of SiC substrate (t = 150 nm), where the effect of the substrate on radiative heat transfer is negligible, is smaller than that of Au substrate (t = Fig. 300 nm). Besides, we obverse that compared with the suspended scenario, the Au substrate always suppresses the radiative heat transfer when t is less than the critical thickness while the SiC substrate slightly enhances RHF as t is larger than 25 nm and smaller than the critical thickness. When t < 25 nm, the SiC substrate plays a negative role in enhancing RHF. The spectral analysis allows us to understand the mechanism for the suppressive effect on NFRHT. {L-End} Fig. 11(b) shows the spectral RHF for two substrate scenarios at d = 20 nm has a prominent feature with new peaks at 9Fig. 39.6 cm⁻¹ and 961.4 cm⁻¹ for SiC substrate and 978.Fig. 3 cm⁻¹ for Au substrate. Moreover, the spectral heat flux of SiC substrate is less than that of the suspended scenario in a spectral region from 766 cm⁻¹ to 927 cm⁻¹, covering RB 1 and RB 2.

The mechanisms by which substrate affects the radiative heat transfer can be explored by phonon tunneling between the two α-MoO₃ heterostructures (i.e., the PTC). Here, we focus on incident frequencies: 780 cm⁻¹ (RB 1), 900 cm⁻¹ (RB 2) and two peaks mentioned above, 939.6 cm⁻¹ and 961.4 cm⁻¹. At the frequency of 780 cm⁻¹, the real part of the permittivity of SiC is positive, resulting in the excitation of the l = 0 mode of HPhPs, as shown in {L-End} Fig. 12(a). The bright branches represent the photon tunneling, originating from the coupling of PhPs between two individual heterostructures. Here, we obtain the dispersion relation from $1-r^{\text{pp}}{r}^{\text{pp}}{\text{e}}^{2\text{i}{k}_{z0}d}=0$⁷⁰, since the contributions to NFRHT in nonmagnetic materials are dominated by transverse magnetic waves, i.e., p-polarized waves^{65-67, 71}. The dispersion curves are hyperbolic and unambiguously located in the bright contour of mode branches. Compared with the suspended case in {L-End} Fig. 12(e), the bright hyperbolic branches are extended to a larger range of wavevectors ~Fig. 300 k₀ for SiC substrate scenario. {L-End} Fig. 12(b) and {L-End} Fig. 12(f) show the PTC at the frequency of 900 cm⁻¹ for suspended and SiC substrate cases, respectively. The bright branches match well with the dispersion curves, yielding hyperbolic features. Combined with {L-End} Fig. 4(a), it can be found that the hyperbolic bright branches in {L-End} Fig. 12(b) represent the l = 1 mode due to Re(${õ}_1$) + Re(${ϵ}_{\mathrm{Fig.\; 3}}$) < 0. Moreover, the wavevector region is also increased to ~250 k₀ compared with ~100 k₀ in the suspended scenario. Whether the increase of wavevector leads to an enhancement or a weakening of spectral RHF can be explained via the coupling of polaritons. Note that the enhancement of NFRHT is dominated by increased photon tunneling channels, which is closely related to the coupling of the surface modes. Note that for small attenuation length, the evanescent wave is easily filtered by the vacuum gap, so the system consisting of two individual surfaces no longer stimulates the obvious polaritonic coupling effect, resulting in the descending of radiative tunnels. Therefore, the strength of coupling depends on the attenuation length and vacuum gap. Previous studies^{68, 71} have suggested that for two bodies with gap distance d, the maximum accessible wavenumber for photon transmission is $\sim \frac{1}{d}$. The attenuation length of HPhPs is $\delta =\frac{1}{\text{Im}\left(k_{z0}\right)}$. We know that $\text{Im}\left(k_{z0}\right)\sim k$ is accessible for large wavevector, so it can be understood that for the maximum accessible wavenumber, the attenuation length is on the order of the gap size⁷². Consequently, only when the attenuation length of the polaritonic mode is larger than the gap distance can the mode make an obvious contribution to the enhancement of total heat flux. As shown in {L-End} Fig. 12(i), the hyperbolic branches move toward the larger-wavevector ones after adopting the SiC substrate. Meanwhile the attenuation length decreases, yielding a weakened coupling effect. Therefore, it explains why the spectral heat flux drops from 0.0276 nW/(m²·rad·s⁻¹) to 0.0215 nW/(m²·rad·s⁻¹) as the α-MoO_{Fig. 3} crystal is placed on the top of SiC substrate. For the case of the 900 cm⁻¹ (shown in {L-End} Fig. 12(j)), the attenuation length of l = 0 mode in the suspended scenario is larger than vacuum gap (20 nm), while that of l = 1 mode in the SiC substrate scenario is smaller than vacuum gap, weakening the radiative heat flux. The new spectral heat flux peaks located at 9Fig. 39.6 cm⁻¹ and 961.4 cm⁻¹ originate from the reorientation of HPs. At 9Fig. 39.6 cm⁻¹ where Re(ε₁) + Re(ε_{Fig. 3})) < 0, the direction of the propagating in-plane l = 0 mode is rotated by 90°(thus occurring along the intrinsically forbidden direction). Additionally, because of the coupling between the evanescent fields of l = 0 mode in the top and bottom heterostructures, the surface mode would split into two resonant branches, i.e., the antisymmetric and symmetric modes, as shown in {L-End} Fig. 12(c) and {L-End} Fig. 12(g). The dispersion relations match well with the bright branches. Different from {L-End} Fig. 12(g) where the antisymmetric and symmetric branches are located in the left and right parts of wavevector regions, these two resonant modes are split into the top and bottom parts, left and right parts of wavevector regions in {L-End} Fig. 12(c). In terms of spectral radiative heat transfer, the reorientation of HPhPs improves the photon tunneling of the system due to stronger coupling in lower k_x-k_y wavevector regions as shown in {L-End} Fig. 12(k). The other spectral heat flux peak is located at ω = 961.4 cm⁻¹, where Re(ε₁) + Re(ε_{Fig. 3}) > 0. For the suspended scenario, the attenuation length of l = 0 mode is so small that the contribution to the radiative heat transfer is negligible ({L-End} Fig. 12(h)). Although no rotation of l = 0 mode occurs, the SiC substrate provides a stronger surface state by reducing the wavevector, improving the photon tunneling of the system effectively as shown in {L-End} Fig. 12(d). This is why the peak occurs.

For the Au substrate scenario, the new spectral heat flux peaks are located at 978.3 cm⁻¹ and 1000 cm⁻¹ as shown in {L-End} Fig. 11(b). From the PTC shown in {L-End} Fig. 13(a), we get that photon tunneling is dominated by the coupling between the evanescent fields of l = 0 mode in the top and bottom α-MoO_{Fig. 3} /Au heterostructures, which is attributed to the spectral heat flux peak. The coupling between l = 1 modes in the suspended scenario is so weak (shown in {L-End} Fig. 13(c)) that the contribution to NFRHT is negligible. At a frequency of 989 cm⁻¹, the coupling effect is obvious in the case with or without Au substrate. However, there is only one resonant branch, matching well with the dispersion contour in {L-End} Fig. 13(b). The blue dashed curve represents the in-plane elliptical IFC, which separates the antisymmetric and symmetric modes. Hence, the only branch is symmetric mode. Namely, the antisymmetric mode is suppressed. This also explains why the spectral heat flux in α-MoO_{Fig. 3} /Au heterostructure is less than that in the suspended scenario.

Considering that α-MoO₃ exhibits different polaritonic properties along the three crystalline directions, it is worthwhile investigating the effect of substrates on NFRHT between heterostructures consisting of α-MoO₃ and substrate along different crystalline directions. We have discussed the heat flux along the [010] crystalline direction above. However, the heat flux along the [001] crystalline direction exhibits distinctive features, as shown in {L-End} Fig. 14 in which the results are obtained with d = 20 nm. Note that in this section the [100], [010] and [001] crystalline directions correspond to the x, y and z axes in the Cartesian coordinate system, respectively. It is clear to see in {L-End} Fig. 14(a) that once the thickness of α-MoO_{Fig. 3} slab exceeds 2.5 nm, the role of Au substrate turns to be enhancive effect and its heat flux is 1.5 folds of that of air scenario when t increases to 7 nm. SiC substrate keeps enhancing the radiative heat transfer until the substrate effect disappears (t > Fig. 300 nm). Particularly, in the case of 1 nm-thickness α-MoO_{Fig. 3}, the heat flux with SiC substrate is more than 2 times as large as that of the air scenario. To explain the physical mechanism behind the enhancing effect, the spectral heat flux of 10 nm-thickness α-MoO_{Fig. 3} is chosen as shown in {L-End} Fig. 14(b) and the inset shows the SiC and Au substrates improve radiative heat flux from Fig. 3.6Fig. 39 kW/m² to Fig. 5.2916 kW/m² and 5.015 kW/m², respectively. Compared {L-End} Fig. 14(b) with {L-End} Fig. 11(b) (heat flux along [010] crystalline direction), the maximum spectral heat flux peak moves toward lower frequency (redshift) and in the case of SiC and Au substrates, the contribution from RB 1 remarkably increases from a qualitative perspective, which attracts our attention. To further clarify the mechanism for the enhancement, the heat flux contribution from every region is calculated and shown in {L-End} Fig. 14(c). The results exhibit that after the substrate is replaced with SiC or Au, the contribution from RB 1 significantly increases from 0.70Fig. 3Fig. 3 kW/m² to 2.2858 kW/m² and 2.7018 kW/m², respectively. Meanwhile, its proportion of total heat flux increases from 19.Fig. 3Fig. 3% to 4Fig. 3.20% and 5Fig. 3.87%, respectively, indicating its dominance. In addition, the heat flux from RB Fig. 3 slightly increases for SiC substrate while dramatically reduces for Au substrate. This phenomenon results from the rotation of the crystal optical axis direction. To be specific, the [001] crystalline direction becomes perpendicular to the surface after rotating 90° around the [100] crystalline direction, i.e., x axis, so that the heat flux is along the [001] direction. This leads to the IFC of PhPs being closed or elliptical in RB 1, while it being open or hyperbolic while that in RB Fig. 3. The coupling of PhPs depends on the polaritons in single slabs. We know that the substrate with negative permittivity can excite the l = 0 mode of elliptical PhPs, which can boost the phonon tunneling. The PTC in {L-End} Fig. 14(d–f) shows the elliptical shape dispersion curves, and only the Au substrate excites the l = 0 mode, leading to the stronger coupling between the two individual heterostructures and the enhancement of heat flux in RB 1. It is noticeable that even though the SiC substrate excites l = 1 mode due to the positive permittivity at a frequency of 760 cm⁻¹, the strong coupling between PhPs in the emitter and the receiver also increases the spectral heat flux, yielding the enhanive effect. Focusing our attention on the RB 2, the reorientation of PhPs occurs for SiC substrate as shown in {L-End} Fig. 14(i). Although the reorientation improves the spectral heat flux in the blue shaded region shown in {L-End} Fig. 14(b), its narrow band limits the enhancive effect on the total heat flux. The largely negative permittivity of Au substrate suppresses the excitation of l = 0 mode, and makes l = 1 mode excited in higher momentum space as shown in {L-End} Fig. 14(h) and {L-End} Fig. 14(k), weakening the coupling of PhPs and decreasing the spectral heat flux. For RB Fig. 3, the positive permittivity of SiC leads to stronger coupling between l = 0 hyperbolic polaritons as shown in {L-End} Fig. 14(l) and increases the photon tunneling.

In conclusion, we present a comprehensive study of the effect of negative permittivity substrate on the hyperbolic phonon polaritons. We have demonstrated the reorientation of the directional propagation of low-loss HPhPs in a natural hyperbolic crystal along intrinsically forbidden directions when the crystal is laid on a substrate with a moderate negative permittivity. In addition, we reveal that the substrate with negative permittivity can generate the fundamental mode which is suppressed in the configuration of the substrate with positive permittivity. Our results show that HPhPs propagate along the naturally forbidden [001] crystalline direction for fundamental polaritonic propagation in bare α-MoO₃ when it is placed on the top of SiC at ω = 936 cm⁻¹. The excitation of the fundamental mode with an elliptical iso-frequency dispersion contour is realized at ω = 989 cm⁻¹ when α-MoO₃ is placed on Au with large negative permittivity. The reoriented and excited fundamental modes are theoretically identified as s-HPs and v-HPs, respectively, via the electrical field distribution and the values of q_ez. Changing the local environment provides new opportunities for dynamically tuning PhPs at the nanoscale.

Thanks to the active modulation of the electromagnetic state of PhPs via the blocal environment, we further investigate the effect of the polaritonic coupling between heterostructures integrated with different substrates on the radiative energy transport. The results show that the reorientation and excitation of the fundamental mode can further improve the NFRHT while the annihilation of the fundamental mode suppresses it. For two kinds of practical substrates: SiC and Au, the enhancement or suppression of radiative energy transport depends on the relative magnitude of the slab thickness and the vacuum gap width. When the nanoscale vacuum gap is smaller than the slab thickness, the nanofilms with and without substrates are as effective as bulk materials, yielding no effect of substrate on radiative energy transport. Finally, the effect of substrates on NFRHT along different crystalline directions is considered. In the case of heat flux along the [001] crystalline direction between two heterostructures integrated with SiC or Au, the spectral band of the excitation of the fundamental mode resulting from the negative permittivity substrates is shifted to RB 1, where the heat flux is greatly improved due to wider band. In this configuration, the contribution from RB 1 dominates and improves the enhancement effect of PhPs compared to the suspended configuration.

It should be pointed out that although α-MoO₃ is chosen for demonstration in this paper, the conclusions obtained apply to other vdW materials, such as α-VO₅ and hBN. Therefore, our findings promise new opportunities for actively modulating heat management at the nanoscale via tailoring the phonon polaritons. Furthermore, this study provides fundamentally relevant insights into understanding in-plane anisotropic polaritons of vdW materials.

This work was supported by the National Natural Science Foundation of China (Nos. 52106099 and 51576004), the Natural Science Foundation of Shandong Province (No. ZR2022YQ57), and the Taishan Scholars Program.

All the authors were involved in designing the research, performing the research and writing the paper.

The authors declare no competing financial interests.

Supplementary information for this paper is available athttps://doi.org/10.29026/oes.2024.240002