Near-field heat transfer between graphene/hBN multilayers

TL;DR

This study investigates how multilayer graphene/hBN structures influence near-field heat transfer, revealing hybrid polaritons' dominant role and the ability to modulate heat transfer via electrical gating and layer number.

Contribution

It provides a detailed analysis of hybrid polaritons in multilayer graphene/hBN structures and highlights the limitations of effective medium theory in predicting heat transfer.

Findings

Hybrid surface polaritons dominate heat transfer in multilayer structures.

Effective hyperbolic behavior enhances photon tunneling and heat transfer.

Heat transfer can be modulated by the number of graphene layers and electrical gating.

Abstract

We study the radiative heat transfer between multilayer structures made by a periodic repetition of a graphene sheet and a hexagonal boron nitride (hBN) slab. Surface plasmons in a monolayer graphene can couple with a hyperbolic phonon polaritons in a single hBN film to form hybrid polaritons that can assist photon tunneling. For periodic multilayer graphene/hBN structures, the stacked metallic/dielectric array can give rise to a further effective hyperbolic behavior, in addition to the intrinsic natural hyperbolic behavior of hBN. The effective hyperbolicity can enable more hyperbolic polaritons that enhance the photon tunneling and hence the near-field heat transfer. However, the hybrid polaritons on the surface, i.e. surface plasmon-phonon polaritons, dominate the near-field heat transfer between multilayer structures when the topmost layer is graphene. The effective hyperbolic regions can be well predicted by the effective medium theory (EMT), thought EMT fails to capture the hybrid surface polaritons and results in a heat transfer rate much lower compared to the exact calculation. The chemical potential of the graphene sheets can be tuned through electrical gating and results in an additional modulation of the heat transfer. We found that the near-field heat transfer between multilayer structure does not increase monotonously with the number of layer in the stack, which provides a way to control the heat transfer rate by the number of graphene layers in the multilayer structure. The results may benefit the applications of near-field energy harvesting and radiative cooling based on hybrid polaritons in two-dimensional materials.

Figures

Equations

$${\sigma }_{\text{D}}=\frac{i}{\omega +i/\tau }\frac{2e^2{k}_{B}T}{\pi {\hslash }^2}\ln \left[2\cosh \left(\frac{\mu }{2k_{B}T}\right)\right]$$

$${\sigma }_{\text{I}}=\frac{e^2}{4\hslash }\left[G\left(\frac{\hslash \omega }{2}\right)+i\frac{4\hslash \omega }{\pi }{\int }_0^{\infty }\frac{G(\xi )-G(\hslash \omega /2)}{(\hslash \omega {)}^2-4{\xi }^2}d\xi \right]$$

$${ϵ}_m={ϵ}_{\infty ,m}\left(1+\frac{{\omega }_{\text{LO},m}^2-{\omega }_{\text{TO},m}^2}{{\omega }_{\text{TO},m}^2-i{\gamma }_m\omega -{\omega }^2}\right)$$

$$q=\frac{1}{4{\pi }^2}{\int }_0^{\infty }\left[\Theta \left(\omega ,T_2\right)-\Theta \left(\omega ,T_1\right)\right]\left[{\int }_0^{\infty }\xi \left(\omega ,\beta \right)\beta d\beta \right]d\omega$$

$$\xi_{j}\left(\omega,\beta\right)=\left\{\begin{array}[]{cr}\frac{\left(1-|r_{1j}|^{2}\right)\left(1-|r_{2j}|^{2}\right)-|t_{1j}|^{2}\left(1-|r_{2j}|^{2}\right)-|t_{2j}|^{2}\left(1-|r_{1j}|^{2}-|t_{1j}|^{2}\right)}{|1-r_{1j}r_{2j}e^{2ik_{z0}d}|^{2}},$$

\;\beta<k_{0}\\ \frac{4\left[\textrm{Im}\left(r_{1j}\right)\textrm{Im}\left(r_{2j}\right)\right]e^{-2|k_{z0}|d}}{|1-r_{1j}r_{2j}e^{2ik_{z0}d}|^{2}},

$$\;\beta>k_{0}\end{array}\right.$$

$$r=r_{12}+\frac{t_{12}{t}_{21}{r}_{23}{e}^{i2k_{z,2}h}}{1-r_{21}{r}_{23}{e}^{i2k_{z,2}h}}$$

$$t=\frac{t_{12}{t}_{23}{e}^{i{k}_{z,2}h}}{1-r_{21}{r}_{23}{e}^{i2k_{z,2}h}}$$

$$r_{ab,s}=\frac{k_{z,a}-{\sigma }_{\text{s}}\omega {\mu }_0-k_{z,b}}{k_{z,a}+{\sigma }_{\text{s}}\omega {\mu }_0+k_{z,b}}$$

$$t_{ab,s}=\frac{2k_{z,a}}{k_{z,a}+{\sigma }_{\text{s}}\omega {\mu }_0+k_{z,b}}$$

$$r_{ab,p}=\frac{k_{z,a}{ϵ}_{\perp ,b}-k_{z,b}{ϵ}_{\perp ,a}+k_{z,a}{k}_{z,b}\frac{{\sigma }_{\text{s}}}{\omega {ϵ}_0}}{k_{z,a}{ϵ}_{\perp ,b}+k_{z,b}{ϵ}_{\perp ,a}+k_{z,a}{k}_{z,b}\frac{{\sigma }_{\text{s}}}{\omega {ϵ}_0}}$$

$$t_{ab,p}=\frac{2k_{z,a}{ϵ}_{\perp ,b}}{k_{z,a}{ϵ}_{\perp ,b}+k_{z,b}{ϵ}_{\perp ,a}+k_{z,a}{k}_{z,b}\frac{{\sigma }_{\text{s}}}{\omega {ϵ}_0}}$$

$${ϵ}_{\perp ,\text{EMT}}=f{ϵ}_{\perp ,\text{G}}+(1-f){ϵ}_{\perp ,\text{hBN}}$$

$${ϵ}_{\parallel ,\text{EMT}}={\left(\frac{f}{{ϵ}_{\parallel ,\text{G}}}+\frac{1-f}{{ϵ}_{\parallel ,\text{hBN}}}\right)}^{-1}$$

$$f=\frac{\Delta }{\Delta +h}$$

Full text

I Introduction

As one of the fundamental modes of heat transfer, radiative heat transfer plays an important role in a wide spectrum of applications from energy harvesting to thermal management [1-5]. In the far field, the maximum radiative heat transfer rate between two objects is restricted by the black-body limit. However, if the two objects are brought very close to a distance comparable to the characteristic wavelength of the thermal radiation, the evanescent waves from each object can couple and assist photons to tunnel through the gap. This is the so-called photon tunneling and the resulted near-field heat transfer rate can be orders of magnitude larger compared to the blackbody limit [6-12]. The enhanced radiative heat transfer finds numerous applications such as thermal energy harvesting, radiative cooling, and thermal imaging [2]. Continuous efforts have been devoted to exploring new materials or structures that can result in large heat transfer rates that can benefit these applications.

Various types of surface polaritons have been extensively studied for their ability to enhance the photon tunneling and greatly boost the near-field heat transfer. Examples include surface phonon polaritons that can exist at the surface of polar dielectric materials such as SiO2 and SiC, or surface plasmon polaritons (SPPs) that can be supported between metallic surfaces or structures [13-17]. Recently, it is demonstrated that surface plasmons in graphene can also achieve a similar role to enhance the photon tunneling between two graphene sheets [18]. Besides the surface polaritons, bulk materials constructed with periodically stacked sub-wavelength metallic and dielectric layers can also enhance near-field heat transfer. The enhancement is originated from the collective response of the multilayers that can be described based on effective medium theory. The effective dielectric function is usually anisotropic, and in some frequency range, the axial and tangential permittivities can even have opposite signs, giving rise to hyperbolic responses [19,20]. In the hyperbolic regions, the isofrequency surfaces become a hyperboloid instead of a sphere or an ellipsoid, and thus such multilayer structures can support resonance modes with unbounded tangential wavevectors. These multilayer structures have found exciting applications in sub-wavelength imaging [21-23] and near-field radiative heat transfer [24,25]. For heat transfers in particular, materials with hyperbolic responses can provide substantial enhancement of heat transfer over a broad frequency region [24].

Recently, it has been shown that the surface plasmons in graphene can couple with the phonon polaritons in hexagonal boron nitride (hBN) films to form hybrid polaritons that greatly enhance the photon tunneling [26]. The structure can yield a larger heat transfer rate than the typical polar materials such as SiO2 and SiC. Since graphene behaves like a thin metallic layer in the structure, one could expect periodically stacked graphene and hBN film would result in a collective hyperbolic response, which may give rise to an enhanced heat transfer rate that can exceed the single-layer heterostructure. Compared to other hyperbolic metamaterials constructed with metal and isotropic dielectrics, such a multilayer structure could enable an actively tunable hyperbolic response by changing the chemical potential of graphene. Note that hBN naturally possesses two mid-infrared Reststrahlen bands that have hyperbolic response [27,28]. Thus, transitions between natural to effective hyperbolic response may also occur in such type of structures.

For these reasons, in this work, we study the near-field heat transfer between multilayers with alternating layers of graphene and hBN film. The coupled modes in hBN/graphene heterostructure are discussed and the evolution of the modes is investigated by exploring the photon tunneling probability between structures with different number of layers. The contributions of the surface polaritons are emphasized by analyzing the spectral heat flux. Predictions based on effective medium theory (EMT) is used to compare with the exact calculations to further illustrate the contributions of the surface polaritons. The heat transfer rates between structures with the same total thickness but different layers are compared.

II CALCULATION OF NEAR-FIELD HEAT TRANSFER

Figure 1 shows the configuration of near-field radiative heat transfer between two periodic multilayer structures separated by a vacuum gap of $d$. Each period of the multilayer stack contains an hBN film with a thickness denoted as $h$ adjacent to a monolayer graphene, and the total number of layers in the structure is denoted by $N$. When $N$ = 1, the structure is simply a heterostructure containing a single layer of hBN and graphene. In Fig. 1, the structure above the vacuum gap is the receiver with a lower temperature $T_{1}$ and the structure below the vacuum gap is the emitter with a higher temperature $T_{2}$. The temperatures are set to be $T_{1}$ = 290 K and $T_{2}$ = 310 K, respectively, and the ambient temperature is assumed to be the same as $T_{1}$. The heat transfer rate is calculated on the receiver [29].

Graphene is modeled with a sheet conductivity, $\sigma_{\textrm{s}}$, that includes the contributions from both the intraband and interband transitions [30], i.e., $\sigma_{\textrm{s}}=\sigma_{\textrm{D}}+\sigma_{\textrm{I}}$, respectively [31]:

[TABLE]

and

[TABLE]

where $G(\xi)=\sinh\left(\xi/k_{B}T\right)/\left[\cosh\left(\mu/k_{B}T\right)+\cosh\left(\xi/k_{B}T\right)\right]$. Here, $e$ is the electron charge, $\hbar$ is the reduced Planck constant, $\mu$ is the chemical potential, $\tau$ is the relaxation time and is chosen to be $10^{-13}$ s for all the calculations, $k_{B}$ is the Boltzmann constant, $\omega$ is the angular frequency, and $T$ is temperature that will be set to either $T_{1}$ or $T_{2}$ in the calculation depending on the location of the graphene sheets [32]. Figure 2 shows the real and imaginary part of the sheet conductivity of graphene at different chemical potentials.

hBN is a uniaxial crystal in the infrared region with two mid-infrared Reststrahlen bands. We assume its optical axis is in the $z$-direction for the structure in Fig. 1. The in-plane and out-of-plane dielectric functions include the contribution from the in-plane phonon vibrations ($\omega_{\textrm{TO},\bot}$ = 1370 $\textrm{cm}^{-1}$ and $\omega_{\textrm{LO},\bot}$ = 1610 $\textrm{cm}^{-1}$) and out-of-plane phonon vibrations ($\omega_{\textrm{TO},\parallel}$ = 780 $\textrm{cm}^{-1}$ and $\omega_{\textrm{LO},\parallel}$ = 830 $\textrm{cm}^{-1}$), respectively, as given by

[TABLE]

where $m=\,\parallel,\bot$[33] denote either the out-of-plane or the in-plane directions, respectively. The other parameters used are $\epsilon_{\infty,\parallel}=2.95$, $\gamma_{\parallel}=4\,\textrm{cm}^{-1}$, $\epsilon_{\infty,\bot}=4.87$, $\gamma_{\bot}=5\,\textrm{cm}^{-1}$. Due to the small damping coefficients as compared to the phonon frequencies, the in-plane and out-of-plane dielectric functions of hBN possess opposite signs in the Reststrahlen bands, making hBN a natural hyperbolic material.

The near-field radiative heat flux $q$ is calculated based on fluctuational electrodynamics [3]

[TABLE]

where $\Theta\left(\omega,T\right)$ is the average energy of a Planck oscillator, $\beta$ designates the magnitude of the wave vector in the $x-y$ plane, and $\xi\left(\omega,\beta\right)$ is the photon tunneling probability (also called energy transmission coefficient). If the integration is done over $\beta$ only, the result is the spectral heat flux. The photon tunneling probability includes contributions of both the transverse electric (TE) waves (or $s$-polarization) and transverse magnetic (TM) waves (or $p$-polarization), that is, $\xi\left(\omega,\beta\right)=\xi_{s}\left(\omega,\beta\right)+\xi_{p}\left(\omega,\beta\right)$. Each polarization contains the contribution from the propagating ($\beta<k_{0}$) and evanescent waves ($\beta>k_{0}$), where $k_{0}=\omega/c_{0}$ is the magnitude of the wave vector in vacuum and $c_{0}$ is the speed of light in vacuum [29]:

[TABLE]

where $j$ is for either $s$ or $p$ polarization, 1 and 2 respectively denote the receiver and emitter, $r$ and $t$ are respectively the corresponding reflection and transmission coefficients, $k_{z0}$ is the $z$-component of the wave vector in vacuum, and Im takes the imaginary part [34].

For structure with $N=1$, the reflection and transmission coefficients for both $s$- and $p$-polarizations take the following forms [3,26]

[TABLE]

and

[TABLE]

where 1, 2, and 3 are the indexes for the vacuum region above hBN film, the hBN film region, and the vacuum region below hBN film, respectively. Also,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, $a$ and $b$ can be 1, 2, or 3, and $\epsilon_{0}$ is the vacuum permittivity, $\mu_{0}$ is the vacuum permeability. Note that the $z$-component of the wavevector in a given region takes different form depends on the polarization. For $s$-polarization, $k_{z,a}=\left(\epsilon_{\bot,a}k^{2}_{0}-\beta^{2}\right)^{1/2}$, and for $p$-polarization, $k_{z,a}=\left(\epsilon_{\bot,a}k^{2}_{0}-\epsilon_{\bot,a}\beta^{2}/\epsilon_{\parallel,a}\right)^{1/2}$. For regions with an isotropic medium like regions 1 and 3, $\epsilon_{1}=\epsilon_{3}=\epsilon_{\bot}=\epsilon_{\parallel}=1$. On the interface without graphene, $\sigma_{\textrm{s}}=0$. For $N>1$, the reflection and transmission coefficients can be obtained by modifying the reflection and transmission coefficient at the interface between the hBN film and the lower vacuum in Eqs. (6) and (7). For example, the structure with $N=2$ has an additional graphene layer and an hBN film added below and above the structure with $N=1$. Therefore, $r_{23}$ is not the reflection coefficient between the interface of two media described by Eqs. (8) and (10), but takes a form that is the same with Eq. (6), excepting that region 1 becomes hBN. The transmission coefficient, $t_{23}$, can be modified in a similar way and takes a form that is the same with Eq. (7) with region 1 being hBN. This process can be repeated to obtain the reflection and transmission coefficients for structures with $N$ layers, and the results are cross-checked using a scattering matrix method [35]. There is an alternative method in which graphene is modeled as a layer of thickness $\Delta$ = 0.3 nm with an effective dielectric function $\epsilon_{\textrm{eff,G}}=1+i\sigma_{\textrm{s}}/\left(\epsilon_{0}\omega\Delta\right)$ [36]. Both methods yield essentially identical results with less than 0.5$\%$ in the predicted total heat flux [26]. The calculations in this work are all based on the above-mentioned analytical expressions. The latter treatment can facilitate the understanding of the physical mechanism, as will be discussed in the following.

III HYBRID POLARITONS IN GRAPHENE/HBN MULTILAYERS

Figure 3(a) demonstrates the photon tunneling probability contours for graphene/hBN heterostructure that contains one unit cell on either side of the vacuum gap ($N$ = 1). The wavevector is normalized using $\beta_{0}=\omega_{0}/c_{0}$ with $\omega_{0}=1\times 10^{14}$ rad/s. Similar to the observations in Ref. [26], when the structure contains only one layer of hBN covered by graphene, hybrid polaritons are formed due to the coupling between hyperbolic phonon polaritons (HPPs) and surface plasmons in graphene. The bright bands indicate the excitation of the hybrid polaritons that enables a high probability of photon tunneling. The polaritons inside the two Reststrahlen bands of hBN are hyperbolic plasmon-phonon polaritons (HPPPs). HPPPs preserve the hyperbolic-waveguide-mode features as in an uncovered hBN film, and they have opposite group velocities in the two Reststrahlen bands. The hybrid polaritons outside the hyperbolic regions are surface plasmon-phonon polaritons (SPPPs), which are surface modes featured with a strong localized field on the interface with graphene. It can be seen from Fig. 3(b) that the SPPPs are the major contribution of the radiative heat transfer, which is 178 kW/m2. There are two branches of SPPPs outside the hyperbolic regions, while the number of branches of HPPPs is affected by the thickness of hBN film. For an infinite thick substrate, the HPPPs merge to form a continuous band.

As the number of layers in the structure increases, more hybrid modes occur as indicated in Figs. 4(a), (b), (c), in which $N$ = 2, 3, and 5, respectively. The heat transfer rate also increases to 196 kW/m2, 208 kW/m2, and 218 kW/m2, respectively. The origin of the increase can be seen clearly from the spectral heat flux for the three cases in Fig. 5. As indicated in Fig. 4, the number of polaritons branches below the lower Reststrahlen band increases as $N$ increases. These additional bands result in a higher spectral heat flux around 1014 rad/s as shown in Fig. 5. The HPPPs extend to the frequencies higher than the upper hyperbolic region of hBN, making the number of the total branches equal to 2$N$. These additional branches do not extend to large frequencies like SPPPs but are bounded within a certain frequency region, and they also lead to a higher spectral heat flux around 3.1 $\times\,10^{14}$ rad/s. The spectral heat flux does not have a noticeable change other than the two mentioned spectral ranges. The observations can be understood by considering the effective behavior of the multilayers based on effective medium theory (EMT), which suggests a uniform property described as [37]

[TABLE]

in which

[TABLE]

is the filling fraction. Note that the optical axis of the multilayer is still in the $z$-direction. Figure 6 shows the dielectric functions of hBN and the effective dielectric functions of the multilayer structure based on Eq. (12) using $\epsilon_{\bot,\textrm{G}}=\epsilon_{\parallel,\textrm{G}}=\epsilon_{\textrm{eff,G}}$. Compared to hBN, EMT predicts a very similar out-of-plane dielectric function, but the in-plane dielectric function is modified significantly due to the metallic behavior of graphene, especially at long wavelengths. Thus, there is a new hyperbolic region formed below the lower Reststrahlen band from 0 to 1.05 $\times$ 1014 rad/s. The high-frequency bond of the upper hyperbolic region of hBN is extended to a slightly higher frequency from 3.03 $\times$ 1014 rad/s to 3.1 $\times$ 1014 rad/s. These changes correspond well with the observations in Fig. 4. The multiple bands in the lower frequency region are bonded by 1.05 $\times$ 1014 rad/s and exhibit a dispersion similar to HPPPs. In the frequency region between 1.05 $\times$ 1014 rad/s to $\omega_{\textrm{TO},\parallel}$, SPPPs are still present due to the lack of hyperbolicity. Meanwhile, due to the extension of the hyperbolic region, HPPPs in the upper Reststrahlen band extend to a slightly higher frequency to 3.1 $\times$ 1014 rad/s. Therefore, EMT provides a qualitative explanation to understand the polariton bands in the multilayer structures. Surprisingly, the effect of the multilayer that creates effective hyperbolic regions can be observed even when $N$ = 2. The thickness of the hBN film can be changed so that the hyperbolic regions can be tuned.

Although EMT gives a qualitative explanation, it cannot capture the details of the polaritons. Figures 7(a) and (b) describe such an effect where the effective properties are used and the thickness of the structure is 50 nm and 250 nm, respectively, corresponding to the cases shown in Fig. 3(a) and Fig. 4(c). The hyperbolic polaritons are well captured by the EMT, though the dispersions are different. The SPPPs, however, do not show up. This can be understood since the mode profile of SPPPs is largely confined on the surface with graphene [33], and this inherent inhomogeneity is not captured by EMT. The effect can be better seen from Fig. 8, which displays the spectral heat flux calculated based on EMT and exact formula for the $N$ = 5 case. The EMT does not capture the peaks in the frequency region where SPPPs exist. Since the heat flux is mainly contributed by the SPPPs, EMT yields a much lower heat flux $q$ = 176 kW/m2 compared to the $q$ = 218 kW/m2. This indicates the importance of the first layer in enhancing the photon tunneling and near-field heat transfer. In fact, exact calculation shows that if the graphene on the top is removed, the SPPPs bands disappear and $q$ drops to 30 kW/m2. The contour plot for $\xi$ looks similar to Fig. 7, though not shown here. Note that surface polaritons also play a critical role in the near-field heat transfer between hyperbolic metamaterials [25,38]. Thus, it can be concluded that EMT is not valid when the surface polaritons dominate the heat transfer.

As shown in the example above, more than half of the contribution to the heat transfer arises from the graphene layers that are immediately adjacent to the vacuum gap, with additional contributions arise from the hyperbolic behavior of the multilayer structures. This observation is in general in agreement with a recent investigation of heat transfer between hyperbolic metamaterials, which showed that a single-layer structure could operate as well or better than hyperbolic metamaterials especially in the limit of small vacuum gap sizes [39]. With the gap sizes that we consider here, there can be significant contributions from both the multi-layers and the top surfaces.

IV EFFECT OF CHEMICAL POTENTIAL AND NUMBER OF LAYERS

In this section, we compare the heat transfer rate for single-layer structure and multilayer structure at different chemical potentials, in order to offer a guidance to choose the optimized chemical potential and number of layers to achieve maximum heat flux. Figure 9 demonstrates the heat transfer rate between two identical structures at different chemical potentials chosen from a set of 0, 0.2, 0.4, 0.6, 0.8 and 1.0 eV, and $N$ is taken from a set of 1, 5, 10, 20, 50, and 100. The total thickness of the structures is fixed at 1 $\mu$m and thus $h$ is one micrometer divided by $N$. At small gap distances, the near-field heat transfer rates are much larger than the black-body limit, i.e., 123 W/m2. Compared to other chemical potentials, $\mu$ = 0.2 eV yields the largest heat transfer rate for all $N$ when $d$ is smaller than 100 nm. The maximum heat transfer rate at $d$ = 20 nm is 303 kW/m2 when $N$ = 50. This rate is much larger compared to that between polar materials (such as SiC and SiO2, yielding 46 kW/m2 and 138 kW/m2, respectively) and plasmonic materials (like heavily doped-Si, yielding 44 kW/m2 at a doping level of 1019 cm*-3* [15]) that are known to yield very large heat transfer rates.

The spectral heat flux corresponding to the maximum heat rate at $d$ = 20 nm (i.e., $N$ = 50 and $\mu$ = 0.2 eV) is shown in Fig. 10. The cases when $N$ = 1 or $\mu$ = 1 eV are also shown in comparison, in which case $q$ = 270 kW/m2 and 68 kW/m2, respectively. The corresponding photon tunneling probability plots are displayed in Fig. 11. The spectral heat flux spectra for the cases with $\mu$ = 0.2 eV are very similar. The majority of the heat flux is contributed by the polaritons below the lower Reststrahlen band and the SPPPs between the two Reststrahlen bands. Increasing the number of graphene layers allows more HPPPs inside the effective hyperbolic regions below the lower Reststrahlen band and above the higher Reststrahlen band as shown in Fig. 11(a), resulting in a higher spectral heat flux in the corresponding frequency regions. The wavevectors of the multiple bands of HPPPs in the hyperbolic regions are very closely spaced and they eventually form a continuous region if $N$ goes to infinity. Note that a larger $N$ does not always result in a larger $q$. For $\mu$ = 0.2 eV, $q$ would decrease to 267 kW/m2 if $N$ = 100. This value is even smaller than the single layer structure with the same total thickness. Thus, one can expect there is a layer number that can maximize the heat transfer, and the number of graphene layer can be used to design a structure that yields a certain heat transfer rate.

The heat transfer rate at $\mu$ = 0.2 eV is relatively large since smaller $\mu$ allows the SPPPs to extend to larger wavevectors [26] as can be seen from Fig. 11(a). It can be seen from Fig. 9 that the heat transfer rate at $\mu$ = 0.2 eV yields the largest $q$ for structures with different $N$. The largest heat transfer rate is achieved around 0.1 eV rather than 0 eV because interband transitions dominate $\sigma_{\textrm{s}}$ in the near-infrared region at $\mu$ = 0 eV and graphene does not support surface plasmons in the wavelength range of interest [40]. When $\mu$ = 1 eV, graphene plasmons expand to cover a wider frequency range as indicated in Fig. 11(c). However, due to the frequency dependence of $\Theta$, the high-frequency SPPPs do not contribute to the spectral heat flux significantly. The polaritons in the lower frequency region do not extend to large wavevectors due to the high $\mu$, and thus the spectral heat flux decrease drastically. The high chemical potential also makes the lower Reststrahlen band become an effective metallic region without a hyperbolicity, which can be seen from the disappearance of the multiple HPPPs bands. For higher emitter temperatures, increasing the chemical potential may result in a larger heat transfer rate since the high-frequency polaritons would be more significant. The strong dependence of the near-field heat transfer on the chemical potential offers another way to actively tune near-field heat transfer besides changing $N$ [41,42]. Note that after the submission of this paper, a paper studying similar system appeared [43].

V CONCLUSIONS

In conclusion, the multilayer structures consisting of graphene and hBN film enable more hybrid polaritons compared to the single-layer structure. HPPPs can be supported in the effectively formed hyperbolic regions that are different from the original hyperbolic regions of hBN. The majority of the near-field heat transfer, however, is still contributed by SPPPs when graphene is the topmost layer. EMT can predict the effective hyperbolic regions but fails to capture the surface polaritons, and thus yields a much lower heat transfer rate compared to the exact calculations. In additional to actively changing the graphene chemical potential, the near-field heat transfer can also be modulated through the number of graphene layers in the structure, which changes the number of the polariton bands. The results demonstrate the possibility to construct hyperbolic metamaterials with two-dimensional materials and may benefit the applications such as near-field energy harvesting and thermal imaging.

Acknowledgement B. Z. and S. F. acknowledges the support from Advanced Research Projects Agency-Energy (ARPA-E), IDEAS program (project title: Demonstration of Near-Field Thermophotovoltaic (TPV) Energy Generation), as well as the support from the Global Climate and Energy Project at Stanford University, and the Department of Energy “Material Interactions in Energy Conversion$\textquotedblright$ Energy Frontier Research Center under Grant No. DE-SC0001293. Z.M.Z. acknowledges the support from the National Science Foundation (CBET-1603761).