Heat diffusion in electronics and optoelectronics is usually disordered, limiting the performance of temperature-sensitive compact macro-/nanodevices. This calls for alternative techniques to control heat flow at the nanoscale. (1−5) One of the emerging strategies, photon-based modulation, and specifically thermal radiation, has attracted extensive interest due to its enhanced heat transfer in the near-field region and the super-Planckian effect. (6−10) Devices such as transistors, (11,12) switches, (13,14) modulators, (15,16) and rectifiers (17−20) have been developed based on near-field energy transport. While most of these works focus on modulating the intensity of heat flux, there is a high demand for thermal routers that can actively configure the heat radiation directions rather than just the intensity.

Generally, a thermal router consists of one heat source and one or more drains (Figure 1a). The direction of heat flow to different drain terminals can be manipulated by configuring the thermal router through various external stimuli, such as magnetic field, temperature, and others. (21−26) For instance, using graphene-coated SiC nanoparticles and via adjusting the Fermi level of the graphene shell, (22,23) thermal routing has been realized since the Fermi level directly modifies the excitation frequency of localized surface mode of nanoparticles, which in turn influences the strength of the localized surface mode coupling between nanoparticles. Tailoring the advanced media with intrinsic direction-dependent optical properties, anisotropic nanoparticles composed of magnetic Weyl semimetals, (24) magneto-optical InSb, (25) or α-MoO₃ (26) are harnessed in the thermal router. Such a scheme, however, remains less experimentally feasible due to difficulties in obtaining such a sample. In general, a simple, configurable, and feasible thermal router is in high demand for thermal applications in integrated high-performance nanodevices.

The thermal router is also related in many-body interactions, (27−29) where more than two terminates (nanoparticles) are included. Radiative thermal heat flow can be manipulated among multiple bodies, thus facilitating the thermal router. Among these, introducing a substrate placed in the vicinity of nanoparticles could enhance near-field energy transports between nanoparticles by orders of magnitude at subwavelength distances, (30−32) which can be further boosted by the coupling with nanophotonic modes such as propagating phonon polaritons. More recent findings revealed richer configurations of thermal routers via harnessing the directional propagation of hyperbolic polaritons (HPs), particularly in in-plane hyperbolic materials. (33−37) Furthermore, with twisted bilayers to modulate the topology of dispersion contours via the twist, the topological transitions at critical angles (i.e., photonic magic angle) can be realized, featuring the collimated group velocity (associated with the canalization effect). (38−44) This raises the potential for long-range directional energy flow between nanoparticles.

We study the concrete design of a radiative thermal router based on three SiC nanoparticles and a rotated monolayer α-MoO₃ slab, as shown in Figure 1c. The material properties of SiC and α-MoO₃ are described by the Lorentz model, (36,37) and the details of permittivity can be seen in the Supporting Information. The centers of the three nanoparticles form an isosceles triangle in the x–y plane with a side length of 1000 nm, i.e., d = 1000 nm, and a vertex angle of 90°. Nanoparticle 1, located at the apex, is a high-temperature source with T₁ = 350 K, and other nanoparticles 2 and 3 are low-temperature drains with T₂ = T₃ = 300 K. The monolayer slab rotates around the z-axis and has an angle of θ₁ with respect to the x-axis (Figure 1c). For the bilayer case, the bottom layer has an interlayer rotation of Δθ with respect to that of the top one. We also define the relative positions of nanoparticles, i.e., the angle between the x-axis and the direction of the source–drain connection line, α_2,3,4···. Without losing the generality, the connecting line between nanoparticles 1 (source) and 2 (drain) is parallel to the x-axis, i.e., α₂ = 0°. For the three-nanoparticle system, α₃ is 90° or varies. For the more nanoparticle system, α₃ and α₄ are fixed to 120 and 240°, respectively. Here, the substrate has the same temperature as the additional drain so that the net power absorbed by drains only comes from the source.

For general purposes, we consider cluster N nanoparticles made of nonmagnetic isotropic SiC material that exchange heat via radiation with each other and the environment. The radius of all nanoparticles is set as 10 nm, which is much smaller than the thermal wavelength so that particles can be modeled as electric dipoles. In addition, the particle–particle distance and particle–slab distance are three times larger than the radius, meaning the higher multipoles can be neglected. (30−32) Under these conditions, the near-field energy transport to body j with temperature T_j due to thermal noise sources in body i with temperature T_i can be expressed as follows where Q_ij denotes the spectral heat flux and Θ_i(ω) = Θ(ω,T_i) = ?ω/(e^?ω/k_BT_i – 1) represents mean energy of Planck oscillators at the temperatures T_i. ? is the reduced Planck constant and k_B is the Boltzmann’s constant. τ_ij is the transmission coefficient from body i to j, which is given by where χ_j = Im(α_j) – (k₀³|α_j|²)/6π and α_i denotes the electric polarizability. M represents the many-body interactions, including the particle–particle and particle–substrate interactions. More specifically, M = A^–1, and interaction matrix A is given by where δ_ij is the Kronecker delta function, δ_ij = 1 only if i = j, δ_ij = 0 otherwise, and $k_0=\frac{\omega }{c}$ is the vacuum wave vector. I₃ is the unit matrix.

From eq 3, we can see that the total Green’s function is separated into a vacuum part, G_ij⁰, and a reflected part, G_ij^R, the latter of which depends on the interface reflection coefficients and goes to zeros in the absence of the interface. Note that the role of the substrate environment on near-field coupling can be directly visualized in the reflecting Green’s function, G_ij^R, which will be analyzed in detail later on. The vacuum Green’s tensor between two particles located r_i and r_j in a vacuum can be written as in which r_ij = r_i – r_j, r_ij = |r_ij| is the magnitude of the vector and ${\hat{\mathbf{r}}}_{\mathit{ij}}=\frac{{\mathbf{r}}_{\mathit{ij}}}{r_{\mathit{ij}}}$. The reflecting Green’s tensor G_ij^R = G^R (r_i, r_j) relates the electric field at r_i, which is generated from the source particle at r_j through the surface reflection. It can be written as follows Here, r_αβ represents the Fresnel reflection coefficient, which can be obtained using the Transfer Matrix Method. The polarization vectors are defined as ${\mathbf{a}}_s^{\pm }=\frac{1}{k_{\rho }}{(k_y,?k_x,0)}^t$ and ${\mathbf{a}}_p^{\pm }=\frac{k_z}{k_{\rho }{k}_0}{(?k_x,?k_y,k_{\rho }^2/k_z)}^t$ where $k_{\rho }=\sqrt{k_x^2+k_y^2}$and $k_z=\sqrt{k_0^2?k_x^2?k_y^2}$ are the lateral and vertical wave vectors, respectively. The detailed derivation can be seen in the Supporting Information.

In this following, we restrict ourselves to a three-body thermal radiative router, while the theory above applies to more general cases. To describe the performance of the radiative thermal router quantitatively, we define the split ratio as η = P₁₂/(P₁₂ + P₁₃). In this case, the perfect radiative thermal routers feature a split ratio of 1 (or 0) when the heat flux from the source, i.e., nanoparticle 1, is entirely routed to drain 2 (or 3). There is no routing effect if the split ratio is 0.5, a trivial case of which is three isotropic nanoparticles in a vacuum. In our work, three isotropic nanoparticles are placed in the vicinity of the anisotropic slab, enabling new photon tunneling channels that are extremely dependent on the directional HPs.

Such photon tunneling channels can be highly effectively modulated via rotation of the slab. Suppose a single α-MoO₃ slab is rotated around the z-axis with the angle, θ₁, which is defined as the relative angle between the x-axis and α-MoO₃ [100] crystalline direction. When θ₁ < 30° (red-shaded region in Figure 2a), the heat flux P₁₂ is several orders of magnitude larger than P₁₃ so that η approaches the maximum value 1, hinting the high-performance radiative thermal routing. When θ₁ > 60° (blue-shaded region), P₁₃ exceeds P₁₂ by orders of magnitude, suggesting the heat flux from nanoparticle 1 is almost entirely directed to nanoparticle 3. This also indicates the robustness of the radiative thermal router (vividly shown in Figure 2b).

Importantly, we first highlight the effect of PhPs in our thermal routers. We plot the spectral heat flux in Figure 2c for the case with and without the α-MoO₃ slab when θ₁ = 0°, respectively. The results show that all spectral peaks have the same frequency, which coincides with the resonance of localized surface phonon polaritons of SiC particles (ω_NP). It indicates that the coupling between the localized surface mode and the propagating HPs dominates the enhancement of near-field energy transport. The enhancive effect is associated with the coupling strength. The heat flux received by nanoparticle 2 is 1.456 × 10^–15 W, nearly 2113 times larger than that received by nanoparticle 3 and 6931 times larger than the case without the substrate, as shown in Figure 2d. Therefore, the PhPs offer phonon tunneling channels for enhanced thermal radiative power.

Since the PhPs in α-MoO₃ are highly in-plane anisotropic, the photon tunneling channels, as enhanced by polaritons, should also feature strong angular sensitivity, offering the new possibility to configure the thermal routing via the twist of the single anisotropic α-MoO₃ slab. To understand this, we plot the real part of the momentum-space reflected green function Re(G^R(1,1)), which represents the element [G^R]₁₁ of the 3 × 3 matrix G^R, when θ₁ = 0° and θ₁ = 90°, respectively, at the resonant dipolar response of SiC nanoparticles (ω_NP) in Figure 3a–d. Since the integral of Re(G^R) is associated with the local density of state, (45) the distribution of Re(G^R) in k space could reflect the polaritonic energy provided by photon tunneling channels. The direction of group velocity (v?_g = ∇_k?ω) is perpendicular to the iso-frequency contour of phonon polaritons, (37,46) and its direction coincides with the Poynting vector defined as $\mathbf{S}=\frac{1}{2}\mathrm{Re}(\mathbf{E}\times {\mathbf{H}}^{*})$, (47) which represents the direction of energy flow. The blue dashed curves represent the dispersion contours, while the red solid lines denote the direction of the group velocity. We notice that the direction of group velocity coincides with the direction of the strongest energy transport of the nanoparticles, i.e., the direction of green dashed arrows in Figure 3e,3f. As shown in Figure 3a,3b with a rotated angle of 0°, the direction of group velocity, which is perpendicular to the iso-frequency contour, is not consistent in the whole wave vector space. Since v?_g has an angle with the x-axis in the large wave vector region, it is mainly distributed along the connection line between nanoparticles 1 and 2, especially in small wave vector regions. This makes the energy transmission path of nanoparticle 1 tightly bunched in one direction, i.e., the x-axis, resulting in the energy transport between nanoparticles 1 and 2, in this case, being larger than that between nanoparticles 1 and 3. The enhancement of heat transfer between nanoparticles 1 and 2 yields the good performance of the radiative thermal router. The physics behind it is similar to the case of θ = 90°; that is, the energy transmission path of nanoparticle 1 is mainly along the interparticle axis, i.e., the y-axis, as shown in Figure 3c,3d.

Figure 3e,3f displays spatial distributions of the radiated electric field energy density, which can be expressed at a point r₀ above the surface and outside the particles as follows where Vr₀, i = ω²μ₀ ∑_j³G_r0,j⁰M_ji (see details in the Supporting Information). One can see that the higher energy density tends to be distributed along specific directions, i.e., the green arrows as shown in Figure 3e, which coincide with the direction of group velocity in the large wave vector region in Figure 3a. Since the energy emitted by high-temperature nanoparticle 1 is mainly limited along the direction, which is not aligned with the x-axis, more energy is transferred to nanoparticle 2 for θ₁ = 0°. The region between nanoparticles 1 and 2 is brighter than that between nanoparticles 1 and 3, also indicating that the energy transmission between nanoparticles 1 and 2 is dramatically enhanced. This results in the energy received by nanoparticle 2 being orders of magnitude larger than that received by nanoparticle 3. In contrast, when the α-MoO₃ slab is rotated by 90° along the z-axis, the energy distribution is concentrated in the region between nanoparticles 1 and 3, thus enhancing near-field energy transport between nanoparticles 1 and 3. Hence, the heat transfer is significantly influenced by the propagating HPs. The control of directional energy flow provides exciting opportunities for the design of radiative thermal routers.

Based on the discussion above, we find that the polariton-enhanced channel can be modulated by rotating the α-MoO₃ slab. Nevertheless, for the single α-MoO₃ slab, the group velocity of polaritons is angularly dispersive, which is determined by the fixed polaritons response of a natural slab, thus limiting the degree of freedom to configure a radiative thermal router. Inspired by the canalization in the twisted α-MoO₃ bilayer system, (37) we now show a rotated bilayer radiative thermal router, as shown in Figure 4a. The top layer has an integral rotation of θ₁ around the z-axis and the bottom layer has a clockwise interlayer rotation of Δθ with respect to the top one around the z-axis (Figure 4b). By introducing a bilayer, the iso-frequency contours of phonon polaritons can be flexibly modulated by the relative rotation angle between the two slabs, (37,38) importantly, undergoing a topological transition from open to closed. At the photonic magic twist angle at which the topological transition arises (shown in Figure 4c), the dispersion contours flatten, yielding diffractionless and low-loss directional canalization (Figure 4d). At the photonic magic twist angle, the real part of the z-component of the electric field exhibits the extremely directional propagation feature (Figure 4e), which is consistent with results from refs [ (37)]. The highly directional phonon polaritons give rise to the single energy transmission channel; that is to say, the energy flows unidirectionally, which is validated by the Poynting vector (Figure 4f). This elucidates that the canalization supported by the twisted α-MoO₃ bilayer can realize collimation and channelization of energy transmission. The performance of the proposed radiative thermal router is significantly influenced by the relative position of the three particles and the reflected substrate (for discussion, see the Supporting Information).

Next, the effect of rotation on photon tunneling channels between nanoparticles is investigated. Here, the relative rotation angle between two slabs remains unchanged; that is, we rotate the whole bilayer. For simplicity, we define an overall rotation angle, ϑ. In detail, we consider the configuration of Δθ = 79.52°, θ₁ = 50.24° as the ϑ = 0°. Δθ = 79.52°, θ₁ = 140.24° is considered as the configuration of ϑ = 90°. Namely, the former (latter) corresponds to that where the group velocity of the twisted bilayer is parallel to the connecting line between nanoparticles 1 and 2 (3) (see the Supporting Information). The results of the thermal router using bilayer structures are shown in Figure 5. One can see that as the rotated angle ϑ increases from 0 to 90°, the heat transfer between nanoparticles 1 and 2 remains stable until ϑ = 45° and then rises rapidly, while that between nanoparticles 1 and 3 slumps followed by stabilization from Figure 5a. In the neighborhood of ϑ = 0°, the value of heat flux between nanoparticles 1 and 3 is 4 orders of magnitude greater than that between nanoparticles 1 and 2, so that η approaches zero, indicating a good performance of the thermal router. Meanwhile, when 85° ≤ ϑ ≤ 90°, the heat flux from nanoparticle 1 is entirely directed to nanoparticle 2; thus, η is closed to 1. As shown in Figure 5b, the energy emitted by the heat source is highly controlled to flow to nanoparticle 3 while it is shifted to nanoparticle 2 after a clockwise rotation of roughly 90°. The spectral heat flux (Figure 5c) further demonstrates that the bilayer can enhance the near-field coupling, and the rotation has a significant effect on the improved photon tunneling channels. However, the enhancement is strongly dependent on the energy flow direction of the bilayers. When ϑ = 90°, the heat flux between nanoparticles 1 and 2 reaches 2.51 × 10^–16 W while that between nanoparticles 1 and 3 is 4.03 × 10^–19 W, the former is 623 times that of the latter (Figure 5d). Compared to the single α-MoO₃ slab, the bilayer enables a more concentrated transmission of energy so that the radiative thermal router is more sensitive to the rotation angle.

The flat dispersion of PhPs in the twisted α-MoO₃ bilayer allows more highly concentrated energy propagation compared to hyperbolic PhPs in a single layer. In addition, the diffraction-free and unidirectionally propagating PhPs in the bilayer provide more degree of angular manipulation on photon tunneling channels via the rotation of the bilayer. We demonstrate the angle-dependent radiative thermal router through Re (G^R) and the radiated electric field energy density, as shown in Figure 6. Compared to hyperbolic curves in Figure 3, the dispersion contour flattens, thus yielding the parallel group velocity. In the case of ϑ = 90° (Figure 6a and b), the direction of group velocity is parallel to the y-axis, the same as the interparticle axis, thus greatly enhancing the heat transfer between nanoparticles 1 and 3. Due to the orthogonality of group velocity and connecting line between nanoparticles 1 and 3, the enhancement of additional radiative channels formed by the substrate on heat transfer between nanoparticles 1 and 2 is weakened. The difference in energy exchange between the source nanoparticle and the two drain nanoparticles is widened, as we expected. One can see from the electric field energy density distribution (Figure 6e) that the energy transmission is collimated and channelized along the axis from nanoparticles 1 to 3. This results in unidirectional energy transport compared with the double directions in Figure 3. The energy absorbed by nanoparticle 3 is much larger than that by nanoparticle 2, which yields a superior performance of the thermal router. In contrast, when ϑ = 90°, the direction of energy transport between nanoparticles is highly limited along the axis from nanoparticles 1 to 2, as shown in Figure 6f. The unidirectional energy transmission of phonon tunneling channels brings more interesting strategies to control the direction of heat flow.

Via introducing the twisted bilayer with collimated group velocity, another important merit is the robustness of a high-performance thermal router, even in long ranges. The influence of the interparticle distance on the split ratio is examined. Here, the vertex angle of the isosceles triangle remains 90° while the side length varies from 50 to 3000 nm. The results are shown in Figure 7a with a rotation of 90°. We observe that η increases from 0.5 to 1 and then remains stable, even though the heat flux between nanoparticles declines as d increases. In particular, when d < 70 nm (7R), the energy received by nanoparticles 2 and 3 is nearly identical, thus leading to no splitting effect. In contrast, η remains unity when d > 400 nm (40R), which indicates the robustness of the radiative thermal router. Compared with results from ref [ (22−26)], our proposed thermal router can achieve good performance with a split ratio of 0 or 1 at longer distances, such as several microns. Note that the contributions to near-field energy transport come from three parts, which can be seen from the value of Tr(χ₁ G₁^* χ₂G₂) (we take the particles 1 and 2 as an example) where the three terms on the right-hand side denote contributions from direct and reflected waves, as well as interference between two waves, respectively. We illustrate their contributions in Figure 7b,c for P₁₂ and P₁₃. For the radiative heat transfer between nanoparticles 1 and 2 (Figure 7b), the contribution from the direct part dominates when the interparticle distance is less than 100 nm. However, with a further increase in the distance, the contribution from the direct part declines and is inferior to that of the reflected part. In comparison, the contribution from the direct part is always larger than that from the reflected part for the heat exchange between nanoparticles 1 and 3, as shown in Figure 7c. The results in Figure 7b,c indicate that the reflected channel opened by the α-MoO₃ bilayer determines the long-range radiative thermal router.

Last, we discuss the effect of spatial configurations of many particles in radiative routers. We first consider the case when the vertex angle of the three particle centers is arbitrary. We change the angular position of nanoparticle 3. For the sake of simplicity, we use the particle relative position representation defined above, that is, the angle of the line connecting nanoparticle 1 and other particles to the x-axis, α_2,3···. For the three-particle system, the connecting line between nanoparticle 1 and nanoparticle 2 is parallel to the x-axis, i.e., α₂ = 0, and we change α₃. By comparing the monolayer configuration in Figure 2a and the bilayer configuration in Figure 5a, one can find that the angle range where η remains 0 or 1 is about 30° for the former and 5° for the latter. In addition, the value of the heat flux of the former is larger than that of the latter. This is due to the orthogonality of the connecting lines between the source and the two drains. We know that the bilayer structure can provide highly unidirectional photon tunneling channels; hence, it may be sensitive to variations of vertex angle of the isosceles triangle. To demonstrate this, the positions of the source and nanoparticle 2 are fixed, while the vertex angle, α₃, varies from 90 to 0°, as shown in Figure 8a, where the distance between the source and drain particles is still 1000 nm. Note that the near-field energy transport between SiC nanoparticles dominates at the nanoparticle resonance. The trend of the spectral heat flux at nanoparticle resonance of ω_NP is almost the same as that for the total heat flux. Hence, we calculate the spectral heat flux, Q, at ω_NP in only the following. The distributions of spectral heat flux versus the rotated angle of the α-MoO₃ bilayer ϑ and the vertex angle of the isosceles triangle α₃ are presented in Figure 8b,c, respectively. We observe that the rotated angle corresponding to the brightest region for Q₁₃ linearly increases as the vertex angle decreases (Figure 8b). This is because only when the direction of nanoparticle arrangement matches with the direction of energy flow controlled by the rotation does the energy transport between nanoparticles reach the maximum. In contrast, the maximum value of Q₁₂ is located at a rotated angle of 90°, independent of the vertex angle (Figure 8c) due to the invariability of the position of nanoparticles 1 and 2. This indicates that the radiative thermal router is sensitive to the variability of the relative positions of the nanoparticles.

We now consider the four-body system as an example. Previous studies (22,24−26) investigated radiative thermal routers consisting of three particles with one source and two drains. Due to the function of unidirectional energy transmission provided by the twisted bilayer structure, the context of the thermal router can be expanded to systems that include more drains. We take a four-particle system as an example, as shown in Figure 8d. The centers of the three drain particles form an equilateral triangle, and the source is located at the center of the triangle; that is, α₃ and α₄ are fixed to 120 and 240°. The distance between the source and any one of the drains is the same, 1000 nm. The spectral heat fluxes between the source and different drains, varying with the rotation of the bilayer structure, are shown in Figure 8e. We observe that the spectral peaks absorbed by the three nanoparticles are located at ϑ = 0, 60, and 120°, respectively. The coupling between LSPhPs from NPs and hybrid HPs from the α-MoO₃ bilayer can reach its strongest when the direction of interparticle energy transmission coincides with that of directional HPs. This indicates that the thermal router can achieve a good splitting effect via rotating the twisted bilayer structure, where the energy emitted by the source can flow in a controlled manner to any of the drains. The ability to control energy flow opens up a new path to design radiative thermal routers based on many-particle systems.

The radiative thermal router could enable exciting opportunities in energy harvesting and heat transfer applications. We report a radiative thermal router with three isotropic particles that could achieve the control of heat flow via the rotation of an α-MoO₃ monolayer or bilayer structure placed in proximity of the particles. For the α-MoO₃ monolayer, it is numerically demonstrated that the energy absorbed by one drain could be several orders of magnitude greater than that received by the other drain, making it possible to change the heat flow routes. Such a thermal router possesses excellent tunability in modulating the split ratio of two drains and maintains an extreme splitting effect in an angle range of up to 30°. This is due to the directional propagating HPs of the α-MoO₃ slab. For the α-MoO₃ bilayer structure, the energy can be canalized along a specific direction at the magic angle. The thermal router can direct a heat flow to a drain of choice by matching the direction of energy transmission and the connecting line between the source and a drain. Our proposed thermal router can realize the control of heat flow over a long distance due to the additional channel for energy transmission between nanoparticles provided by the coupling of the propagating hybrid HPs and the localized surface mode of SiC nanoparticles. In addition, we find that the rotated angle, where the splitting effect is best, is deeply sensitive to the vertex angle of the isosceles triangle. We expand the thermal router to four- or more-particle systems due to unidirectional energy transmission. Hence, this work can provide guidance to the active control of photon-based energy transport at the nanoscale and holds potential applications in energy harvesting and thermal management.