Enhancing photon tunneling probability is the key to increasing the near-field radiative heat transfer between two objects. It has been shown that hexagonal boron nitride (hBN) and graphene heterostructures can enable plentiful phononic and plasmonic resonance modes. This work demonstrates that heterostructures consisting of a monolayer graphene on an hBN film can support surface plasmon–phonon polaritons that greatly enhance the photon tunneling and outperform individual structures made of either graphene or hBN. Both the thickness of the hBN films and the chemical potential of graphene can affect the tunneling probability, offering potential routes toward passive or active control of near-field heat transfer. The results presented here may facilitate the system design for near-field energy harvesting, thermal imaging, and radiative cooling applications based on two-dimensional materials.

## Introduction

The well-known diffraction limit in optical imaging is caused by the disappearance of evanescent waves in the far field, which contain fine feature information of the sample surface and are critical for constructing images with high resolution [1]. The decay of the evanescent waves at large distances not only reduces the optical imaging resolution but also limits the radiative heat transfer between two objects in the far field. Thermal radiation is initiated by the fluctuational motion of charges inside materials at above absolute 0 K and contains both propagating and evanescent waves [2,3]. If the distance between two media is greater than the characteristic wavelength of thermal radiation, i.e., in the far field, only propagating waves generated by one medium can reach the other and contribute to radiative heat transfer. The disappearance of evanescent waves in the far field limits the radiative heat transfer to a rate that cannot exceed the blackbody limit governed by the well-known Stefan–Boltzmann law. The limited radiative heat transfer greatly hinders the application where thermal radiation plays a major role, such as radiative energy harvesting, thermal management, and local heating and cooling [3–8]. To overcome this limit, the evanescent waves have to be collected, which can be done by enabling photon tunneling.

When two objects are at a distance comparable to or shorter than the characterized thermal radiation wavelength, i.e., in the near-field regime, the forward and backward evanescent waves can couple with each other and open paths for photons to tunnel through. This phenomenon is called photon tunneling and there are more tunneling photons than propagating photons, resulting in a radiative heat transfer rate that can be orders of magnitude higher than the blackbody limit. The huge radiative heat flux in the near field opens the door to various applications like thermal rectification [9], thermophotovoltaics [10,11], noncontact refrigeration [12], and information processing [13]. Since a large heat transfer is of critical importance in these appealing applications, continuous efforts have been devoted to exploring innovative optical materials that can enhance photon tunneling.

Surface modes, such as surface plasmon polaritons (SPPs) or surface phonon polaritons (SPhPs), have been demonstrated to mediate photon tunneling between metallic surfaces and polar materials [14,15]. Nanostructures like photonic crystals, nanowire arrays, and gratings could effectively behave as hyperbolic metamaterials, which support propagating waves with large wave vectors and thus enhance photon tunneling [16–19]. The two-dimensional (2D) materials can enable plentiful resonances [20–26] to facilitate photon tunneling. For example, graphene supports SPPs that can enhance photon tunneling between graphene sheets [27,28]. Graphene SPPs can also couple with the hyperbolic modes in nanowires, resulting in nearly perfect photon tunneling [29]. When graphene is patterned to ribbons, hyperbolic plasmons can be excited that can significantly enhance photon tunneling than with continuous graphene [30,31]. As a natural hyperbolic material, hexagonal boron nitride (hBN) can support multiple orders of phonon–polaritonic waveguide modes in its two infrared Reststrahlen bands [32–34]. Recently, it has been demonstrated both theoretically and experimentally that van der Waals heterostructures assembled by graphene and hBN can support surface plasmon–phonon polaritons (SPPPs), which are resulted from the strong coupling between the phonon polaritons in hBN and the surface plasmons in graphene [35–37]. It is still a question whether SPPPs can enhance photon tunneling between such heterostructures and enable a higher heat flux than other 2D materials like graphene. Meanwhile, since graphene has been proposed for a number of promising nano-electronic applications, and hBN is an ideal substrate supporting high-quality graphene owing to its planar hexagonal lattice structure [38,39], it is imperative to explore the effect of hBN on graphene regarding radiative heat transfer performance.

The present study theoretically investigates the effect of SPPPs in enhancing the photon tunneling between graphene and hBN heterostructures. The near-field radiative heat transfer between the heterostructures calculated based on fluctuational electrodynamics is compared with the scenarios where only graphene monolayers or hBN films are present. The dispersion of the hybrid polaritons in the heterostructure and their contributions to the near-field heat transfer are discussed. The effects of thickness of the hBN film and the chemical potential of graphene on the dispersion are investigated. A second layer of graphene on the backside of the heterostructure is further explored for its impact on the dispersion and near-field heat transfer of the heterostructure.

## Theory

Figure 1(a) shows the configuration of near-field radiative heat transfer between two aligned heterostructures separated by a vacuum gap of *d*. Each heterostructure contains a monolayer graphene covered on an hBN film with a thickness denoted as *h*. The upper one is the emitter with a relatively higher temperature *T*_{1} and the lower one is the receiver with a lower temperature *T*_{2}. Graphene is modeled with a sheet conductivity, $\sigma s$, that includes the contributions from both the interband and intraband transitions [29]. In the mid- and far-infrared region, $\sigma s$ is dominated by the intraband transitions and can be approximately written as

where *e* is the elementary charge, $\u210f$ is the reduced Planck constant, *ω* is the angular frequency, *τ* is the relaxation time, and *μ* is the chemical potential [25]. Here, $\tau =10\u221213\u2009s$ is chosen for all the calculations. Figures 2(a) and 2(b) show the real and imaginary parts of $\sigma s$ with different chemical potentials, where $\sigma 0=e2/4\u210f$. The curves in Figs. 2(a) and 2(b) are almost identical with the prediction from Eq. (1), except for $\mu =0.1\u2009eV$ toward the high-frequency end where interband transitions become important.

*z-*direction, with two mid-infrared Reststrahlen bands. The in-plane (when the electric field is perpendicular to the optical axis) and out-of-plane (when the electric field is parallel to the optical axis) dielectric functions include the contribution from the in-plane phonon vibrations ($\omega TO,\u22a5$ = 1370 cm

^{−1}and $\omega LO,\u22a5$ = 1610 cm

^{−1}) and out-of-plane phonon vibrations ($\omega TO,\u2225$ = 780 cm

^{−1}and $\omega LO,\u2225$ = 830 cm

^{−1}), respectively, as given by

where *m* = $\u2225,\u22a5$ [36]. The other parameters used are $\epsilon \u221e,\u2225=2.95$, $\gamma \u2225=4\u2009cm\u22121$, $\epsilon \u221e,\u22a5=4.87$, and $\gamma \u22a5=5\u2009cm\u22121$. According to Eq. (2), for small damping coefficient *γ*, the dielectric function becomes negative in Reststrahlen band between the TO and LO phonon modes. Therefore, the in-plane and out-of-plane dielectric functions of hBN possess opposite signs in either Reststrahlen band. As a result, the isofrequency contour becomes hyperbolic in these regions, making hBN a natural hyperbolic metamaterial. Figure 2(c) shows the real part of the dielectric function of hBN and the two hyperbolic regions are marked by the shaded areas.

*q*is calculated based on fluctuational electrodynamics using dyadic Green's functions [2]

*ω*,

*T*) is the average energy of a Planck oscillator, $\beta $ designates the magnitude of the wavevector in the

*x–y*plane, and $\xi (\omega ,\beta )$ is the photon tunneling probability (also called energy transmission coefficient). 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 (\omega ,\beta )$ = $\xi s(\omega ,\beta )$ + $\xi p(\omega ,\beta )$. Here

where *j* is for either *s* or *p* polarization, *r _{j}* signifies the corresponding reflection coefficient, and Im takes the imaginary part [18]. The magnitude and

*z*-component of the wavevector in vacuum are denoted as $k0$ and

*k*

_{z}_{0}, respectively.

*a*and

*b*, where

*a*= 1 or 2 and

*b*= 1, 2 or 3. The effect of graphene is included as a current sheet. If there is no graphene in between media

*a*and

*b*, then

*a*) and (6

*b*), $\epsilon 0$ is the vacuum permittivity and $kz,b=(\epsilon \u22a5,bk02\u2212\epsilon \u22a5,b\beta 2/\epsilon \u2225,b)1/2$ with

*b*being 1, 2 or 3 in Eqs. (5) and (6) is the

*z*-component of the wavevector in a given region. For regions with isotropic medium like regions 1 and 3, $\epsilon 1=\epsilon 3=\epsilon \u22a5=\epsilon \u2225=1$. For TE waves, the reflection coefficient can be expressed as

Here, *μ*_{0} is the permeability of vacuum and *μ _{b}* (

*b*= 1, 2 or 3) is the relative permeability for region

*b*, which is unity for all regions since the materials are all nonmagnetic. Note that in Eqs. (7) and (8), $kz,b=(\epsilon \u22a5,bk02\u2212\beta 2)1/2$ since TE waves are ordinary waves in the hBN film. As mentioned before, the near-field radiative heat transfer is dominated by TM waves. Equations (5) and (7) can also be used for structures with only graphene by setting $r23=0$. They can also be used for the structure without graphene, with one graphene sheet as shown in Fig. 1, or with a graphene monolayer on both sides of the hBN film to be discussed later. An alternative method can also be used is to treat graphene sheet as a layer of thickness Δ = 0.3 nm with an effective dielectric function $\epsilon eff=1+i\sigma s/(\epsilon 0\omega \Delta )$ [40]. Both methods yield essentially identical results with less than 0.5% in the predicted total heat flux. The analytical expressions of the reflection coefficients presented here not only help elucidate the fundamental mechanisms of the coupled plasmonic resonances (to be discussed later) but also can save a lot of simulation time once implanted in the numerical algorithm.

## Results and Discussion

Figure 3 compares the heat fluxes between a pair of graphene sheets, hBN films with *h* = 50 nm, and the heterostructures shown in Fig. 1(a). The chemical potential of graphene is set to 0.37 eV, according to a previous experiment [37]. In the calculations, the temperatures of the emitter and receiver are set as *T*_{1} = 300 K and *T*_{2} = 0 K, respectively. However, the optical properties of graphene and hBN are evaluated at room temperature of 300 K. In terms of the radiative heat flux, the heterostructure outperforms the other configurations, especially at small gap distances. At *d* = 10 nm, the heterostructure yields *q* = 800 kW/m^{2} that is more than twice of that between graphene monolayers or hBN films, which are 305 and 212 kW/m^{2}, respectively. When *d* exceeds about 200 nm, the heat flux for the heterostructures is very close to that between suspended graphene sheets, indicating a negligible effect of the hBN film. Note that the radiative heat flux between blackbodies (in the far field) is 459 W/m^{2}, which is orders of magnitude smaller than the near-field heat flux shown in Fig. 3.

The mechanism for the enhanced heat transfer between the heterostructures can be elucidated by the contours of photon tunneling probability displayed in Fig. 4 for the three scenarios. Note that for the structures considered in this work, $\xi p\u226b\xi s$ since the polaritons discussed here can only be excited for TM waves, which are the dominating contribution to the near-field radiative flux. The bright bands shown in Fig. 4 indicate efficient photon tunneling due to the excitation of different polaritons, corresponding to the dispersion curves where the denominator of $\xi p$ in Eq. (4) approaches zero. The two bands in Fig. 4(a) correspond to the symmetric (lower frequencies) and asymmetric (higher frequencies) branches of the coupled SPPs between two graphene sheets. They are the major contributors to the high near-field heat flux between graphene [30]. For hBN films shown in Fig. 4(b), multiple phononic waveguide modes can be identified in each Reststrahlen band between the horizontal dashed lines. However, there exist strong tunneling branches outside the two Reststrahlen bands for the heterostructure, as shown in Fig. 4(c). These additional polaritons, identified as SPPPs, are one kind of hybrid polaritons resulted from the coupling between surface plasmons in graphene and phonon polaritons in hBN [36,37]. Note that similar to the Fig. 4(a), each order of the polariton bands in Figs. 4(b) and 4(c) splits into two branches. The mechanism of hybrid polaritons and their impact on near-field radiative transfer are elaborated in the following.

When $\epsilon \u22a5,2\epsilon \u2225,2<0$, the isofrequency contour of hBN is hyperbolic and the hybrid polaritons are referred to as hyperbolic plasmon–phonon polaritons (HPPPs) [37] since they preserve the hyperbolic-waveguide-mode features as in an uncovered hBN film. Integer *n* denotes the resonance order. Here, $\delta =\xb1i\epsilon \u2225,2/\epsilon \u22a5,2$ and the plus or minus sign is chosen based on the shape of the HPPPs bands [36]. If $\sigma s$ is set to zero, the first expression in Eq. (9) recovers the dispersion for the waveguide modes in hBN films. When $\epsilon \u22a5,2\epsilon \u2225,2>0$, the isofrequency contour of hBN becomes elliptic, and SPPPs can be supported in the three frequency regions below, between, and above the two Reststrahlen bands, as shown in Fig. 4(c). Unlike HPPPs, SPPPs are surface modes with a strong plasmonic characteristic when they are not so close to the Reststrahlen bands [37]. If $\sigma s$ is set to zero (i.e., without graphene), the second expression in Eq. (9) yields a negative *β* since $\epsilon \u22a5,2>0$, suggesting that hBN films cannot support any resonances outside the two Reststrahlen bands. This provides an explanation why the resonances only exist within the Reststrahlen bands as Fig. 4(b) shows.

It is worth pointing out that the second expression in Eq. (9) does not require an anisotropic substrate to yield a valid solution. Thus, phononic polaritons in films made of isotropic polar materials like SiO_{2} and SiC could also couple with SPPs in graphene to form a hybridized polariton, whose dispersion can be obtained by setting $|\delta |=1$. Similar phenomena were demonstrated for the coupling between SPPs in graphene or thin metal layer with the SPhPs supported by semi-infinite polar material substrates [41,42]. Moreover, a prominent feature that can be identified in Fig. 4(c) is the mode flatting when the SPPPs approach either Reststrahlen band near $\omega TO$. This is caused by the anticrossing effect or mode repulsion between SPPs in graphene and waveguide modes in hBN, and same effect was found to exist between SPPs in a thin metal layer and SPhPs in isotropic polar substrates [41]. Because of this effect, the two hyperbolic Reststrahlen bands break the otherwise continuous SPPPs into three regions in frequency, allowing high density of state to occur at some band edges and boosting the photon tunneling, as will be discussed in the next.

If Eq. (3) is integrated over $\beta $ only, the result is the spectral heat flux, which is shown in Fig. 5 for the three configurations with *d* = 20 nm. It can be seen that the heat flux for hBN structures is mainly contributed by the waveguide modes, and the spectral heat flux is nearly zero outside the two Reststrahlen bands. Graphene plasmons cover a frequency range up to $6\xd71014\u2009rad/s$ as shown in Fig. 4(a), and thus the spectral heat flux between graphene sheets has a nontrivial value over a broad band but fades away at high frequencies due to the frequency dependence of the Planck oscillator. The heterostructure combines the features of graphene and hBN film, and also shows a high spectral heat flux in the regions outside the Reststrahlen bands due to SPPPs, indicating that the major contribution to the radiative heat transfer between the heterostructures comes from SPPPs instead of HPPPs. The peak around $1.4\xd71014\u2009rad/s$ is caused by the high density of state [43] of the SPPPs near $\omega TO,\u2225$, as indicated by the flat dispersion in Fig. 4(c), which is a result of the above-mentioned anticrossing effect. The flat dispersion allows high *β* modes to exist and makes the term $\u222b0\u221e\xi (\omega ,\beta )\beta d\beta $ in Eq. (3) greater, leading to an increase in the spectral heat flux. However, the anticrossing effect plays the opposite role and suppresses the heat flux in the hyperbolic regions. As shown in Fig. 4(b), the waveguide modes in hBN film are able to be supported when *β* is about 150*k*_{0}. In contrast, the anticrossing causes an early truncation of the HPPPs, as shown in Fig. 4(c), and these modes disappear at about 100*k*_{0}. Similar effect was observed for the anticrossing between graphene SPPs and SPhPs in SiC and SiO_{2}, where SPPs truncate the SPhPs at large wavevectors [42,44]. This early truncation results in a lower spectral heat flux for the heterostructures in the two hyperbolic bands compared to the hBN films, as indicated in Fig. 5. The effect of loss on the hybrid polaritons is also investigated by setting the loss of one of the materials to zero, though not shown here. When the graphene is assumed to have no loss, calculations show that the contribution of SPPPs to the spectral flux becomes much less because the bands of SPPPs are much narrower. However, when the loss of hBN is set to zero, SPPPs are not much affected and the spectral heat flux resembles that shown in Fig. 5. Although HPPPs can slightly extend to larger wavevectors, the spectral heat flux in the two hyperbolic regions actually drops somewhat due to the decrease of the bandwidth of HPPPs.

The photon tunneling can be controlled by *h* since a thicker film allows higher orders of HPPPs to present at lower *β*. The upper panel of Fig. 6 demonstrates this effect, in which the hBN film is changed to be 200 nm in (*a*) and semi-infinite in (*b*). Compared to the scenario in Fig. 4(c) when *h* = 50 nm, more orders of HPPPs show up in Fig. 6(a) and they eventually merge to a continuous band when the hBN layer becomes semi-infinite. However, the SPPPs do not experience significant changes. This can be understood since SPPPs are surface modes featured with a strong localized field on the interface of heterostructure that has graphene. The field intensity evanescently decays away from the interface and, thus, making the hBN film thicker has little effect on SPPPs. Since the enhanced heat flux is mainly due to the contribution of SPPPs, a thicker hBN film brings a slight increase in the heat flux, yielding 537 and 547 kW/m^{2} for Figs. 6(a) and (b), respectively, compared to 513 kW/m^{2} for the case in Fig. 4(c). This effect may offer a passive way to modify photon tunneling to some extent. Note that calculations show that the penetration depth in hBN is smaller than 300 nm over the frequency range with high spectral heat flux [45]. Subsequently, a thicker film exceeding 300 nm would have little effect to near-field heat transfer. On the other hand, a lower *μ* can push the SPPPs to larger wavevectors as indicated by Eq. (9), and thus can affect the tunneling probability significantly. Figures 6(c) and 6(d) show the tunneling probability when *μ* is changed to 0.2 and 0.6 eV, respectively. The SPPP bands become flatter in (*c*) and steeper in (*d*) compared to Fig. 4(c), and thus, *q* becomes 1131 kW/m^{2} and 249 kW/m^{2}, respectively. The huge increase at *μ* = 0.2 eV is attributed to the flat dispersion that allows the SPPPs to extend to larger $\beta $ values. The same reason holds for the decrease of heat flux at *μ* = 0.6 eV. Note that this does not mean the maximum heat flux can be achieved at *μ* = 0 eV because at this chemical potential, interband transitions dominate *σ*_{s} in the near-infrared region in which graphene does not support SPPs anymore [46,47]. Calculations show that the maximum heat flux occurs at around *μ* = 0.1 eV with *q* = 1700 kW/m^{2}. Since *μ* can be changed by electrical gating [20], this effect may offer a potential way to actively control photon tunneling.

The photon tunneling probability can be further enhanced by placing an extra graphene layer on the other side of the graphene/hBN heterostructure. In this case, both sides of the hBN film are covered with graphene and the graphene on the backside can modify coefficient $r23$ such that additional solutions exist that can zero the denominator of $\xi p$, as shown in Fig. 7. Compared to Fig. 4(c), two new branches of SPPPs occur outside the two hyperbolic Reststrahlen bands. These SPPP branches further enhance photon tunneling as demonstrated by the spectral heat flux shown in Fig. 8, and the total heat flux *q* is increased to 643 kW/m^{2} (a 25.3% increase for the same structure with only one side covered graphene). Note that the chemical potentials of the two graphene layers in this example are the same and they can be actively changed to tune the SPPPs as discussed previously. However, they do not have to be the same and may be controlled independently to allow more active tunability.

## Conclusions

The work demonstrates that graphene/hBN heterostructured optical materials can enhance photon tunneling and outperform graphene or hBN film in terms of achieving high near-field radiative heat transfer, thanks to the hybrid mode SPPPs. The SPPPs enhance the spectral heat flux outside the two Reststrahlen bands and can be controlled by the thickness of the hBN film and the chemical potential of graphene. Placing an additional graphene layer on the other side of the hBN film further enhances the heat transfer between the heterostructures by allowing extra SPPP branches. The findings in this work may facilitate the design of systems utilizing near-field thermal radiation with passively and actively tunable photon tunneling based on graphene and hBN.

## Acknowledgment

This work was supported by the National Science Foundation (CBET-1235975; CBET-1603761).

## Nomenclature

*d*=vacuum gap width (m)

- $e$ =
charge of an electron (1.602 $\xd7$ 10

^{−19}C)*h*=film thickness (m)

- $\u210f$ =
reduced Planck's constant (1.055 $\xd7$ 10

^{−34}J·s)*i*=$\u22121$

*k*_{0}=wavevector in vacuum (m

^{−1})*k*_{z}_{0}=*z*-component of the wavevector in vacuum (m^{−1})*n*=resonance order

*q*=heat flux (W/m

^{2})*r*=reflection coefficient

*T*=temperature (K)

*β*=magnitude of the wavevector in the

*x–y*plane (m^{−1})*γ*=scattering rate (rad/s)

*δ*=variable defined in Eq. (9)

*ε*=dielectric function

*ε*_{0}=vacuum permittivity (8.854 × 10

^{−12}F·m^{−1})- Θ =
mean energy of the Planck oscillator (J)

*μ*=relative permeability or chemical potential (eV)

*μ*_{0}=permeability of vacuum (4

*π*× 10^{−7}H·m^{−1})- $\xi $ =
photon tunneling probability

*σ*_{s}=sheet conductivity (S)

*τ*=relaxation time (s)

- $\omega $ =
angular frequency (rad·s

^{−1})