In this work, we study the near-field radiative heat transfer between two suspended sheets of anisotropic 2D materials. It is found that the radiative heat transfer can be enhanced with orders-of-magnitude over the blackbody limit for nanoscale separation. The enhancement is attributed to the excitation of anisotropic and hyperbolic plasmonic modes. Meanwhile, a large thermal modulation effect, depending on the twisted angle of principal axes between the upper and bottom sheets of anisotropic 2D materials, is revealed. The near-field radiative heat transfer for different concentrations of electron is demonstrated and the role of hyperbolic plasmonic modes is analyzed. Our finding of radiative heat transfer between anisotropic 2D materials may find promising applications in thermal nano-devices, such as non-contact thermal modulators, thermal lithography, thermos-photovoltaics, etc.

Thermal radiation is an important physical phenomenon. Any object with temperature *T*>0K emits electromagnetic (EM) waves due to the fluctuating current generated from thermal motion of charge carriers. According to the Stefan-Boltzmann law, the radiative heat flux between two separated black bodies is given as $Sbb=\pi 2kB460\u210f3c2(T14\u2212T24)$, where *k*_{B} is the Boltzmann constant, *ℏ* is the reduced Planck’s constant, *c* is the speed of light in vacuum, *T*_{1} and *T*_{2} is the high and low temperatures, respectively. For the Stefan-Boltzmann law, the separation distance *d* is much larger than the thermal wavelength λ_{th}=*ℏ*c/k_{B}*T* and only the propagation modes are taken into account in the process of radiative heat transfer. However, an additional contribution, i.e., evanescent waves, is dominant and should be considered when the separation distance *d* is much smaller than *λ*_{th} (about 10 μm at room temperature).^{1,2} The evanescent waves near the surface can be surface plasmon polaritions (SPPs),^{3} surface phonon polaritions (SPhPs),^{3,4} or even frustrated modes from hyperbolic materials.^{5,6} Due to large density of states of evanescent waves, near-field radiative heat transfer (NFRHT) can exceed the blackbody limit by several orders of magnitude, which were reported in a number of configurations.^{7–11} The control of NFRHT has many promising applications such as near-field imaging,^{12} thermos-photovoltaics^{13,14} heat-assisted magnetic recording,^{15} etc.

The control of NFRHT based on isotropic two-dimensional (2D) materials such as graphene have been reported.^{16–21} The frequency of SPPs of doped graphene lies in the mid-infrared, which guarantees a high-efficient excitation of SPPs at room temperature or above. Moreover, the property of large electrical tunability enables graphene as an excellent platform for active control of NFRHT.^{16–21} Recently, anisotropic 2D materials, such as black phosphorus (BP),^{22–25} rhenium disulfide (ReS_{2})^{26,27} and trichalcogenides,^{28,29} have received enormous attentions. Particularly, monolayer or few-layer of black phosphorus, showing direct bandgap, high mobility and large current on/off ratios, is a promising 2D material for applications in optoelectronic devices.^{30–32} Compared with the well-known graphene sheets, in-plane anisotropic is one of the most intriguing properties of anisotropic 2D materials. Under properly doping, SPPs can be supported in anisotropic 2D materials in the frequency of THz to infrared regime.^{33} In addition to the features of sub-wavelength and electrical tunability, the SPPs in anisotropic 2D material are highly anisotropic. Moreover, natural hyperbolic SPPs can be supported in anisotropic 2D material under specific conditions.^{33} These unique features of SPPs of anisotropic 2D materials are highly desired for the control of near-field thermal radiation.

In this work, we investigated the NFRHT between two suspended sheets of anisotropic 2D materials. Based on the fluctuation-dissipation theory, the radiative heat flux between two separated sheets of anisotropic 2D materials is calculated. The enhancement can be orders-of-magnitude over the blackbody limit as the separation distance between 100 nm and 10 nm. It is revealed that the enhancement of NFRHT is not only related with the excitation of anisotropic plasmonic modes, but also with hyperbolic plasmonic modes. Moreover, the in-plane anisotropic properties of anisotropic 2D materials can provides a capability for thermal modulation based on rotation.

The system under study is depicted in Fig. 1. Two parallel sheets of anisotropic 2D materials are suspended in vacuum and the spatial separation between them is *d*. The temperatures for the bottom and upper sheets are respectively *T*_{1} and *T*_{2}. A twisted angle (denoted by ϕ), i.e., relative crystalline orientation between the upper and bottom sheets of anisotropic 2D materials, is tunable through rotation. The optical conductivity for anisotropic 2D materials can be described by a tensor:^{33}

where

*v=x, y,* where *e* is the charge of an electron, *n* is the concentration of electrons, *m*_{v} is the effective mass of electron along the *v* direction, η corresponds to the relaxation time, Θ(*ω*-*ω*_{v}) is a step function defining the absorption due to the interband transition, *s*_{v} represents the strength of interband transitions, and ω_{v} is the frequency of the onset of interband transitions for the *v* component. The optical conductivity of a single sheet of anisotropic 2D is shown in Fig. 2(a) for a set of parameters. As the frequency is low (e.g., ω <0.15 eV/*ℏ* for *n*=5.0×10^{13} cm^{-2}), the imaginary parts of σ_{xx} and σ_{yy} own the same sign, i.e., Im[σ_{xx}] Im[σ_{yy}]>0, which produces elliptic-like anisotropic plasmonic dispersion. At high frequency regime (ω>0.15 eV/*ℏ*), it has Im[σ_{xx}] Im[σ_{yy}]<0, resulting in interesting hyperbolic plasmonic dispersion.^{33}

The radiative heat flux between two sheets of anisotropic 2D materials can be calculated based on the fluctuation-dissipation theory, which is given as follows:^{34}

where Θ(ω, *T*)=*ℏ*ω/(exp(*ℏ*ω/*k*_{B}*T*)-1)is the average energy of a Planck’s oscillator for angular frequency ω at temperature *T*, and *ξ*(ω, *k*_{x}, *k*_{y}) is the energy transmission coefficient, standing for the probability of photon tunneling. The energy transmission coefficient *ξ*(ω, *k*_{x}, *k*_{y}) are contributed from propagating and evanescent modes, expressed as follows:^{34}

where $k||=kx2+ky2$ and $\gamma =k02\u2212k||2$ are respectively in-plan and vertical wavevectors, *k*_{0}=ω/c is the wavevector in vacuum, $D12=(I\u2212R1R2e2i\gamma d)\u22121$ and **R**_{j} (*j*=1, 2) is the 2×2 reflection matrix for the *j*-th sheet, having the form:

where the superscripts *s* and *p* represent the polarizations of transverse electric (**TE**) and transverse magnetic (**TM**), respectively. The matrix element *r*^{α, β} (*α*, *β* = *s*, *p*) represents the reflection coefficient for an incoming α-polarized plane wave turns out to be an outgoing β-polarized wave. Consider an incident plane wave with in-plane wavevector $k||=kxe\u2192x+kye\u2192y$, the reflection coefficients for single sheet of anisotropic 2D materials can be obtained as follows:^{35,36}

where the elements of conductivity tensor are normalized by free-space impedance $\mu 0/\epsilon 0$, and the primes in superscripts stand for the rotation with respect to unit vector of $k||$. For the bottom sheet, we have:

For the upper sheet, the anisotropic conductivity with respect to unit vector of $k||$ is given by $RT\sigma \xaf'R$, due to an additional rotation angle ϕ, where $R=cos\u2061\varphi \u2212sin\u2061\varphi sin\u2061\varphi cos\u2061\varphi $ and $RT=cos\u2061\varphi sin\u2061\varphi \u2212sin\u2061\varphi cos\u2061\varphi $. Interestingly, the off-diagonal elements of reflection matrix $rs,p=rp,s=0$ for isotropic 2D materials because of $\sigma \xaf'xy=\sigma \xaf'yx=0$. However, the off-diagonal terms of reflection matrix may nonzero for anisotropic 2D materials due to the anisotropic conductivity. Noted that the pole of reflection coefficients in Eq. (6) can produce the dispersion of SPPs of anisotropic 2D materials.^{33}

According to above equations, the spectral transmission coefficient *ξ* (ω), i.e., the integral of *ξ* (ω, *k*_{x}, *k*_{y}) over the whole *k*-space, are shown in Fig. 2(b) for twisted angle ϕ=0^{o} and ϕ=90^{o}. Clearly, the spectra, depending on the twisted angle, are consisted of frequency regimes for both anisotropic and hyperbolic SPPs. The spectral transmission coefficient *ξ* (ω) shows a broadband feature for ϕ=0^{o}. However, the magnitude of *ξ* (ω) drop enormously and a strong suppressing of radiative spectrum is found at frequency regime (0.25 *eV*∼0.35 *eV*) for ϕ=90^{o}, which may find application in thermal spectral engineering. Noted that the hyperbolic regime for anisotropic 2D materials appears at high frequency, in contrast to that for patterned graphene ribbons.^{18}

To understand the broadband mechanism of *ξ* (ω) in Fig. 2(b), the contours plots of *ξ* (ω, *k*_{x}, *k*_{y}) at *k*-space are shown in Fig. 3 for three different photon energies. Firstly, we consider the case of ϕ=0^{o} (upper panel of Fig. 3). For low photon energy *ℏ*ω=0.1eV, Im[σ_{xx}]Im[σ_{yy}]>0 and an ellipse-like diagram can be observed in Fig. 3(a), confirming the excitation of anisotropic SPPs. For high photon energies *ℏ*ω=0.17 eV and 0.25 eV, Im[σ_{xx}]Im[σ_{yy}]<0 and hyperbolic configurations can be found in Fig. 3(b) and Fig. 3(c). The broadband spectrum for ϕ=0^{o} is because of the coexistence of strong anisotropic and hyperbolic SPPs. Meanwhile, the large *k* vector (*k*_{x}, *k*_{y} >>*k*_{0}) for both anisotropic and hyperbolic regimes guarantee a giant enhancement of radiative heat flux. For comparison, the corresponding isofrequency curves of SPPs dispersion of a single anisotropic 2D sheet are illustrated by the dashed gray lines. The strong coupling of SPPs for ϕ=0^{o} enables the original dispersion (dashed gray lines) splits into two bright branches in Figs. 3(a–c).

Now we consider the case for ϕ=90^{o} (bottom panel of Fig. 3). The contours of *ξ* (ω, *k*_{x}, *k*_{y}) show quite complicated diagrams, due to the mismatching coupling. For photon energies *ℏ*ω=0.1 eV and 0.17 eV, the most effective coupling appears at *k*_{x}=*k*_{y} as shown in Fig. 3(d) and (e). Nevertheless, the corresponding density of states becomes smaller compared with those of ϕ=0^{o}. As the photon energy increases to *ℏ*ω=0.25 eV, strong suppressing of radiative spectrum is confirmed by the contour plot of *ξ* (ω, *k*_{x}, *k*_{y}) in Fig. 3(f). Clearly, the magnitude of *ξ* (ω, *k*_{x}, *k*_{y}) is almost 100 times smaller than those of ϕ=0^{o}, due to the weak coupling of SPPs. Correspondingly, the dispersion of the two separated sheets is almost the same as those of original dispersion of single sheet (dashed gray lines).

The radiative heat flux between two sheets of anisotropic 2D material, normalized by blackbody limit, are shown in Fig. 4(a). Remarkably, orders-of-magnitude enhancement can be found, especially at small separation for ϕ=0^{o} (e.g., about 2800 for *d*=10 nm in Fig. 4(a)). The enhancement become smaller as the twisted angle become 90^{o}, indicating a thermal modulation effect based on twisted angle. In order to quantify the modulation effect, the ratios of radiative heat flux <*S*(ϕ)>/<*S*(ϕ=0^{o})> are given in Fig. 4(b). It can be seen that the thermal modulation is strong as ϕ increases from 0^{o} and it is almost unchanged as ϕ>45^{o}. The ratio <*S*(ϕ=90^{o})>/<*S*(ϕ=0^{o})> as a function of separation distance *d* is also given in the inset of Fig. 4(b). The ratio can reduce to almost 25% for *d*=10 nm while it increases slightly as *d* increases. It is worth mentioning that similar thermal modulation effect ware reported for two metallic/polar gratings^{37} and two anisotropic hBN plates.^{38} For anisotropic 2D materials, no pattern nanofabrication is needed, which is much simple for practical applications. Besides, the thermal modulation between anisotropic 2D materials can be enlarged further incorporated with tuning the doping level of concentration. The NFRHT as a function of concentration *n* are shown in Fig. 4(c), which decreases greatly as *n* varies from 1×10^{13} cm^{-2} to 10×10^{13} cm^{-2}. Due to the exponentially decaying of Boltzmann factor, low frequency of SPPs (corresponds to low density of concentration) is preferred to enhance NFRHT. The enhancement of NFRHT at *d*=10 nm can reach about 10^{4} times for n∼1×10^{13} cm^{-2}, whereas it can drop to about 10^{3} times for *n*∼10×10^{13} cm^{-2}.

To understand NFRHT contributed from different SPPs modes. The spectral radiative heat flux for *n*∼5×10^{13} cm^{-2} and *n*∼2×10^{13} cm^{-2} are given respectively in Figs. 4(d) and 4(e) with *T*_{1}=310 K, *T*_{2}=*T*_{1}-10 K, and *d*=50 nm. For n∼5×10^{13} cm^{-2}, anisotropic regime is dominant for NFRHT, while the contribution from hyperbolic regime occupies only a small part. As *n* decreases to 2×10^{13} cm^{-2}, the contribution from hyperbolic regime increases, but still smaller than those from anisotropic one. Nevertheless, the radiative spectrum from hyperbolic SPPs can be enlarged as the temperature increases, due to the fact that the hyperbolic regime appears at relative high frequencies. Taking *T*_{1}=610 K as an example in Fig. 4(f), the radiative heat flux generated from hyperbolic SPPs are comparable to those from anisotropic one. Due to the coexistence of anisotropic and hyperbolic SPPs, the spectral radiative heat flux for anisotropic 2D materials is broadband, which can be 1 order broaden over than those of hBN.^{38}

In summary, NFRHT between two suspended sheets of anisotropic 2D materials is investigated. The radiative heat flux can exceed the blackbody limit over 10^{2} to 10^{4} times for separation distance between 100 nm to 10 nm, stemming from the excitation of anisotropic and hyperbolic SPPs. The hyperbolic SPPs play an important role when the concentrations of electron are low and the temperatures are relative high. A large thermal modulation effect is demonstrated based on the twisted angle of principal axes between the upper and bottom sheets. Besides, active modulation of NFRHT can also be realized by tuning the concentrations of electron though external gating. Noted that the effective masses and band gaps of anisotropic 2D materials are tunable through external mechanical strain,^{39,40} or out-of-plane static electric fields,^{41} which may provide alternative modulations of NFRHT without any rotation of the sheets. Our work may pave the way for applications of anisotropic 2D materials in non-contact thermal modulators, thermal lithography, thermos-photovoltaics, etc.

*Noted added.* After submission of this paper, we became aware of a preprint by Y. Zhang et al. (2018),^{42} where the authors consider the NFRHT for suspended monolayer and multilayer black phosphorus. They came to some similar conclusions. It should be pointed out that the reflection coefficients here are calculated by a simple way (the anisotropic 2D material is treated as a conductive surface), which is different from Y. Zhang et al. (2018).^{42} Besides, the role of hyperbolic plasmonic modes in NFRHT is analyzed separately in our work.

This work is supported by the National Natural Science Foundation of China (Grant No. 11747100), and the Innovation Scientists and Technicians Troop Construction Projects of Henan Province. The research of L. X. Ge was further supported by Nanhu Scholars Program for Young Scholars of XYNU.

## REFERENCES

_{2}due to electronic and vibrational decoupling

_{2}probed by Raman spectroscopy and scanning transmission electron microscopy

_{3}nanoribbon transistors