Gate voltage and doping effects on near-field radiation heat transfer in plasmonic heterogeneous pairs of graphene and black phosphorene

Plasmon coupling and hybridization in 2D materials plays a significant role for controlling light–matter interaction at the nanoscale. We present a near-field radiation heat transfer (NFRHT) between vertically separated graphene and black phosphorene sheets at different temperatures in nanoscale separations. Radiation exchange from the theory of fluctuation electrodynamics is modulated by the carrier density of graphene and phosphorene. Direct comparison of NFRHT black phosphorene–graphene to symmetric graphene–graphene radiation exchange can be as much as 4 times higher for the selected doping range in both armchair (AC) and zigzag (ZZ) orientations of BP. The strong NFRHT enhancement of the specific optical properties of the heterogenous 2D material is due to the strong coupling of propagating surface plasmon polaritons as demonstrated by the distribution of the heat transfer coefficient. We also demonstrate that the magnitude of the near-field radiation enhancement is found to acutely depend on the vacuum gap of the graphene and BP pair. Interestingly, for separation distances below 200 nm, the total near-field heat transfer between black phosphorene and graphene exceeds that between graphene and graphene by 5 times. The radiation enhancement can be further tuned based on the orientation, AC, and ZZ of black phosphorene. These results prominently enable dynamic control of the total NFRHT relying on tunable anisotropic characteristics of BP irrespective of graphene's optical conductivity. Furthermore, the heterogeneous pairs of 2D materials potentially provide alternative platforms to achieve beyond super-Planckian radiation.


Introduction
Near-eld radiative heat transfer (NFRHT) between two bodies separated by a distance much smaller than or comparable to the thermal wavelength can break the blackbody (BB) limit. 1,2 Indepth understanding of this phenomenon is not only essential for fundamental scientic importance but also offers a wide range of potential applications such thermal nanoimaging, 3 near-eld thermophotovoltaics, 4 and thermal rectiers. 5 The NFRHT enhancement by several orders of magnitude has been widely investigated in materials that support phonon polaritons, 6 resonant surface phonon polaritons (SphPs), 7 surface plasmon polaritons (SPPs) 8 or hyperbolic phonon polariton (HPP) modes 9,10 due to evanescent waves. Control of NFRHT has been covered in several structural arrangements and materials; [11][12][13][14][15][16][17] however, a challenge arises in enhancing NFRHT due to the large number of independent parameters, and more importantly, lack of degree of freedom for active tunability.
In this context, a breakthrough for active tuning mechanics has been effectively achieved in two-dimensional (2D) semiconductors: prominently graphene or MoS 2 and, very recently, black phosphorene. [18][19][20][21] Graphene exhibits highly tunable optical properties via voltage-gating, 18,22 good thermal stability, 23 and high carrier mobility. 18 Moreover, graphene supports strongly conned surface plasmons with a frequency that varies with the wave number. [24][25][26] Theoretically predicted tuning and coupling of surface mode plasmon of a graphene sheet with a nearby body, resulting in a signicant increase of the NFRHT, [27][28][29][30] has been conrmed experimentally. 31 In addition, the near-eld active thermal radiation response of graphene in subwavelength geometries, graphene ribbon, 32,33 periodic nanoribbon, 34 or nanodisks [35][36][37] has demonstrated enhanced radiation. Most recently, Black phosphorus (BP) is also reported to have a layer dependent tunable direct bandgap, 38 thermal stability, 39 high carrier mobility, 20,21 and remarkable in-plane anisotropic electrical and optical properties. 40 In contrast to graphene, BP shows polarization dependency in the anisotropic surface plasmon response due to an order of magnitude difference in the in-plane effective electron masses along the two crystal axes results. 41,42 These plasmon properties allow for the realization of novel 2D polarization dependent and tunable NFRHT. 33,[43][44][45] Another interesting approach to study is the NFRHT between heterogenous 2D material, shown to be possible using a graphene sheet and a MoS 2 sheet. 46 However, both graphene and MoS 2 are nearly isotropic materials. This makes the graphene and MoS 2 orientation indiscriminate in light matter interaction. Recently, the NFRHT of the graphene (GR) heterogenous (GR-BP) in contact with bulk substrate SiC has shown anisotropic effects in near eld radiation heat transfer. 47 However, SiC supports hyperbolic phonon polariton modes in the midinfrared range. In hybrid structures composed of 2D material, graphene/SiC or BP/SiC, SPPs can couple with hyperbolic phonon polariton. 47 The hybrid modes move from the substrate and introduce a change in the propagation, and evanescent waves either greatly improve or reduce the NFRHT by creating a new channel. Therefore, it is important to understand the NFRHT and tuning capability graphene and BP heterogeneous system without dispersive hyperbolic layer.
Through numerical and analytic solutions based on the transfer matrix method combined with uctuation electrodynamics, 48 the near-eld radiation transfers energy between the BP and the graphene, which are compared to two identical suspended graphene sheets by analysing the distribution of the heat transfer coefficient of the emitter and the receiver 2D layers. The NFRHT model shows that BP-GR has a more selective optical modulation compared to the graphene-graphene system. Furthermore, the NFRHT differences between these two systems are explored, considering the inherent BP optical anisotropy nature. The effect of the reorientation of BP is found to be signicant in modulating the NFRHT, perhaps creating possible conditions roughly equivalent to or greater than graphene. The hybrid modes, which evolve from low to high parallel wave vectors in the infrared range, controlled by graphene and BP optical parameters determine the tunnelling probability. We also discuss the total radiation, enhancement in the 2D structure relative to the value of separation distance, offering potential routes toward passive or active control of NFRHT. In particular, we examined the ratio of total heat transfer of BP-GR and GR-GR, and its dependence on the anisotropy, doping, and spacing.

Optical properties of graphene and black phosphorene
The optical conductivity in graphene, s g (u), is modelled with a 2D sheet conductivity that includes the contribution of an intraband transition term s D (i.e., Drude), and interband transitions term, s I , where s g (u) ¼ s D (u) + s I (u). 49 and GðxÞ À Gðuħ=2Þ ], e is the electron charge, ħ is the reduced Planck constant, k B is the Boltzmann constant, m is the chemical potential, T is temperature, u is the angular frequency, and s is the relaxation scattering time, which we have xed to the value s ¼ 10 À13 s. 50 The anisotropic layer dependent surface conductivity, s jj , of BP is given by a simple semiclassical Drude model in the infrared regime, 40 where j represents AC and ZZ-direction h describes the relaxation rate and is chosen to be 10 meV 40 which is within the range obtained from ab-simulation 51 and D i ¼ pne 2 /m j is the anisotropic Drude weight, n is the electron doping, and m j represents the layer dependent effective mass of electrons along the j, (AC or ZZ) direction. In this work we use monolayer BP with armchair, m x ¼ 0.15m 0 , and zigzag, m y ¼ 0.7m 0 , direction masses, where m 0 is the free electron mass. 40 Near-eld radiation heat transfer Before setting the model of radiation heat transfer, we examined the light-matter interaction for the suspended isotropic 2D materials. As shown in Fig. 1, each two-dimensional layer material is represented by the conductivity tensor s ij , consisting of transverse conductivities, longitudinal conductivities, and crossed conductivity components, which lies in the xy-plane. We obtain the Fresnel coefficients solving Maxwell's electromagnetic equation imposing the appropriate electromagnetic eld boundary condition and representing the 2D conducting sheet surface current (J s ) from Ohm's law. We assume that the electric eld E i in medium i, and E j in medium j, i and j indicated medium above and below the 2D layer, with dielectric permittivity 3 i ¼ 3 j ¼ 1.0. The corresponding magnetic eld components, H i and H j , can be found using Maxwell's equation.
The boundary conditions of the parallel and perpendicular The unit vector normal to the 2D layer isn ¼ Àẑ, and J s ¼ s ij E. Solving the boundary conditions, the Fresnel coefficients of reection for s, p and the cross-polarization term expressions are 52,53 where i and j are the mediums that encapsulate graphene or BP. r pp ij , r ss ij and r sp ij are the Fresnel coefficients between the medium i and j for s, p and the mixture sp polarization. The z-component wavevector for s and p-polarization along each region is , and 3 0 is the vacuum permittivity, m 0 is the vacuum permeability and b is the wave vector parallel to the surfaces of the 2D layers. This work studies BP AC orientation aligned with x-direction and ZZ along y-axis and vice versa. The incident polarization direction is collinear to one of the principal BP directions so that the crossed conductivity effect, s xy ¼ 0, results in the cross-polarization r ps ¼ 0. A similar method has been followed recently by Xiao-Jie Yi et al. 47 A rotated BP from the optical axis does change the momentum vector of along x and y and hence the expression of the reection spectra. 54 Now the NFRHT initiated by the coupling of the propagating surface modes meditated between the thermally hot (T E ) and the thermally cold (T R ) surfaces can be formulated following the framework of uctuational electrodynamics as: 55 is the average energy of Planck's oscillator, and x s(p) (u,b) is the energy transmission coefficient that depends on the frequency and wave vector corresponding to s(p)polarization of light. The energy transmission coefficient considering polarization coupling and contribution of evanescent modes, b > k 0 , and propagating modes, b < k 0 , is given by 56 where j stands for either s or p polarization, 01 and 12 respectively denote the receiver and emitter medium reection, k 0z is the normal wavevector component in vacuum of BP and graphene, Im() takes the imaginary part of a complex quantity, and d is the distance of the vacuum gap. Throughout the paper, we consider the conguration in Fig. 1, where two 2D bodies with conductivity labelled as s jj,E or s E (emitter) and s R (receiver), respectively, are brought into proximity with a vacuum gap separation d. The temperatures of the emitter and the receiver are labeled T E ¼ 300 K and T R ¼ 290 K, (T E > T R ), respectively. We represent gate modulation or applied voltage bias effect in terms of changing the carrier density n, of graphene and BP. This indirect relation of gatemodulation of BP and graphene can be controlled independently leading to more degrees of freedom. The carrier density potential of graphene can be described also in terms of chemical potential, The values of n g used in this work ranges from 0.74 Â 10 12 cm À2 to 4.7 Â 10 13 cm À2 , with an equivalent range of m ¼ 0.1 eV to 0.9 eV, and n b ¼ 2.0 Â 10 12 cm À2 to 9.0 Â 10 13 cm À2 in order to get optimum NFRHT.

Results and discussion
Doping and gate modulation effect First, we consider the scenario in Fig. 2a-c, BP emitter (BPE) at temperature of 300 K and graphene receiver (GRR) at 290 K with the AC direction of BP aligned with x-direction and ZZ is in the y-direction. The separation distance d between graphene and BP is 10 nm. The gure shows the spectral heat ux (SHF) calculated from eqn (7) as a function of wavelength for GR low doped n g ¼ 0.74 Â 10 12 cm À2 , 6.6 Â 10 12 cm À2 and high doped 2.65 Â 10 13 cm À2 , which corresponds to a chemical potential of m ¼ 0.1 eV, 0.3 and 0.6 eV, respectively, and BP varying the dopinglevel, n b ¼ 0.5 Â 10 13 cm À2 , 1.0 Â 10 13 cm À2 , 5.0 Â 10 12 cm À2 .
In all gures of Fig. 2a-c for comparison, we also included spectral heat ux (SHF) results calculated between two graphene sheets, graphene emitter (GRE)-graphene receiver (GRR) (blue broken line labelled as GG), see Fig. 1b. The same temperatures are used in all cases: T E ¼ 300 K and T R ¼ 290 K. We selected identical gate voltage for the emitter and receiver graphene, n gE ¼ n gR ¼ n g since the latter case shows strong spectral heat ux. 38 Therefore, low gate modulated NFRHT between identical graphene shows enhanced radiative energy transport than between low and highly gate modulated graphene. Further support of this conclusion is provided in the ESI Fig. 1 Schematic illustration of the emitter and receiver with temperatures T E (hot) and T R (cold) 2D material. The vacuum gap distance, d, and optical conductivity, s, of the graphene sheet and the black phosphorene sheet are labeled. (a) Graphene sheet emitter (E) and black phosphorene sheet receiver (R), and (b) graphene sheet emitter and graphene sheet receiver. Fig. 1 (a) and (b). † As seen in Fig. 2a, b and c, the SHF of identical graphene has higher n g ¼ 0.74 Â 10 12 cm À2 and decreases more than ten times when the doping is 2.65 Â 10 13 cm À2 . In addition, at higher graphene carrier density, n g ¼ 2.65 Â 10 13 cm À2 , SRH becomes broader due to weakening of interband effect. 27,29,38 In the case of BPE-GRR, increasing number density n b from 0.2 Â 10 13 cm À2 to 5.0 Â 10 13 cm À2 has a dramatic effect on the SHF at low carrier density, n g ¼ 0.74 Â 10 12 cm À2 , Fig. 2a. The SHF peaks decrease gradually with the increase of number density by more than an order of magnitude. In addition, the SHF is broadening gradually increases and the peak value shis slightly to a higher wavelength with increasing n. Compared to BPE-GRR for all n g and n b , GRE-GRR is higher than BPE-GRR at n g ¼ 2.65 Â 10 13 cm À2 . The role of the dependency of n g and n b on the SHF is further evaluated for 6.6 Â 10 12 cm À2 and 2.65 Â 10 13 cm À2 , Fig. 2b and c. On the contrary, for n g ¼ 4 Â 10 12 cm À2 SHF for BPE-GRR is higher than the that of GRE-GRR from moderate to high doping, Fig. 2b, and for all doping level of BP, Fig. 2c. More specically, for n g ¼ 6.6 Â 10 12 cm À2 the SHF of BPE-GRR is smaller than GRE-GRR and larger by a factor of two for n b ¼ 5.0 Â 10 12 cm À2 and 5.0 Â 10 13 cm À2 . This indicates in BPE-GRR there is a critical n b value between 5.0 Â 10 12 cm À2 and 1.0 Â 10 13 cm À2 making the transition of enhanced SHF than GRE-GRR, in Fig. 2b. For highly gated graphene, n g ¼ 2.65 Â 10 13 cm À2 , we observe that BPE-GRR 2D system outperforms identical graphene sheets, within the whole range of n b , having two times higher SHF for n g ¼ 1.0 Â 10 13 cm À2 and 5.0 Â 10 13 cm À2 .
So far, we have discussed the case of doping effect of BPE-GRR in AC orientation. However, besides electrostatic tuning of graphene and BP, effective mass of BP is readily different in the AC and ZZ direction, as shown in the optical conductivity eqn (2). In what follows, we rotate the BP arrangement as such the ZZ direction is along the x-axis and AC direction is along the yaxis. For this purpose, we show in Fig. 3 the dependence of the SHF (all parameters are same as in Fig. 2). The change in the orientation of the BP strongly modulates the height of the NFRHT spectra by the maximum peak energy and slightly changes the peak position. For example in Fig. 3a, the intermediate doping number density, when n b is close to 1.0 Â 10 13 cm À2 for n g ¼ 0.74 Â 10 12 cm À2 , reaches a minimum or a comparable density to GRE-GRR system. In the ZZ direction the SHF can be tuned by a factor of 2 to 4 between n b ¼ 0.5 Â 10 13 cm À2 and n b ¼ 1.0 Â 10 13 cm À2 .
Noticeably at higher doping and lower doping, the SHF falls well below the GRE-GRR system. However, at both lower and higher doping levels, the NFRHT spectra of the BPE-GRR system does not reach the GRE-GRR value at 0.74 Â 10 12 cm À2 , a consistent result of AC orientation (Fig. 2a). Likewise, as in the AC orientation case, the NFRHT of BPE-GRR is greater than the GRE-GRR system when n g ¼ 6.6 Â 10 12 cm À2 and n g ¼ 2.65 Â 10 13 cm À2 (Fig. 3a and b). To be more precise, when the doping level is bellow n b ¼ 0.5 Â 10 13 cm À2 we obtain a result less than the BPE-GRR, not shown in the gure. However, different from Fig. 3a, when n g ¼ 6.6 Â 10 12 cm À2 and n g ¼ 2.65 Â 10 13 cm À2 , the SHF gradual spectra broadening and amplitude increase is also noticeable in the BPE-GRR system with an increase in n. This is due to the weakening of the interband loss in graphene. Both results in Fig. 2 and 3 indicate that irrespective of orientation and doping level of BP, the GRE-GRR system leads to higher radiation transfer than the BPE-GRR system for a small applied bias on the graphene. The gate modulation of graphene, anisotropic mass effect and the carrier density of the BP, manifests itself through the variations in the radiation transfer response in the heterogeneous 2D material.
The contours show two bright bands covering a wide range of frequencies and momentum. These bright bands correspond to the symmetric (lower frequencies) and antisymmetric (higher frequencies) branches of the coupled SPPs between 2D materials. 27,28 Here one mode is originated from the BPE side and a second from the GRR; they are labelled SPB and SPG, respectively. The band are characteristics reections of the SPPs dispersion band and an indication of electromagnetic coupling process of BP and graphene.
The strength of the two bright band colors indicates a high transmission coefficient and the coupling strength of the SPPs of black phosphorene and graphene. Therefore, they are the major contributors to the SHF between black phosphorene and graphene. In Fig. 4a, we observe strong nearly overlapping and small separation of the bright band transition probability for n b ¼ 0.5 Â 10 13 cm À2 and for n g ¼ 0.74 Â 10 12 cm À2 . The upper branch is BP like, and lower branch is graphene like. Short momentum and frequency spanning in the lower branch are due to the interband transition damping, specically at energies ħu < 2m. 26,28,34 Fig. 3 Spectral heat flux between suspended GRE-GRR and BPE-GRR, at T E ¼ 300 K, and black phosphorene sheet zigzag (ZZ) direction along x-axis orientation, T R ¼ 290 K, separation distance, d ¼ 10 nm, electron carrier density n b ¼ 5.0 Â 10 12 cm À2 , 1.0 Â 10 13 cm À2 , 5.0 Â 10 13 cm À2 , and blue dash line in all the figures is between two suspended graphene sheets at identical optical property n g ¼ 0.74 Â 10 12 cm À2 (a), n g ¼ 6.6 Â 10 12 cm À2 (b), and n g ¼ 2.65 Â 10 13 cm À2 (c). Fig. 4 Color map of photon transmission probability (x) between suspended BP emitter and graphene sheet receiver (G), and black phosphorene sheet (B) emitter arm chair (AC) direction along x-axis orientation separation distance, d ¼ 10 nm, electron density of graphene, n g ¼ 0.74 Â 10 12 cm À2 , and BP, n b , (a) 5.0 Â 10 12 cm À2 , (b) 1.0 Â 10 13 cm À2 , (c) 5.0 Â 10 13 cm À2 , respectively and n g ¼ 6.6 Â 10 12 cm À2 , and BP, n b , (d) 5.0 Â 10 12 cm À2 , (e) 1.0 Â 10 13 cm À2 , (f) 5.0 Â 10 13 cm À2 , respectively. However, strongly overlapping bright bands are produced for low doping BP and graphene, resulting in overall a higher NFRHT, Fig. 3a. The two resonances branches show a progressive shi and move further apart, especially for a large n b , indicating mismatch between the hybrid SPPs modes of the BP emitter and the graphene receiver, Fig. 4b and c. The motions toward a lower wave vector result in a decrease in transmission and the NFRHT, Fig. 2b and c, due to weak coupling. In case of n g ¼ 6.6 Â 10 12 cm À2 , Fig. 4d-f, both branches cover a large frequency range as the interband damping effect of graphene weakens. However, progressive increase of n b similarly increases movement of the bright band and shis towards higher u and a lower b. The rise in strong offsets of the branches lead to the overall reduction of the NFRHT observed at a large ng (Fig. 2b-c) and ESI Fig. 3. † In order to understand the signicance of anisotropy effects in the photon transmission and near-eld radiation heat transfer, we look in detail at Fig. 5, where BP ZZ-direction is aligned in x and AC is aligned in the y-direction, having the same parameters as in Fig. 4. Similar to the hybridization effect of the surface plasmons of a BP-GR discussed above, two bright bands are observed. As expected, the effective mass difference has a strong inuence on the heat transmission. On the one hand, the greater effective mass of electrons in ZZ direction results in energies of the bright bands that extends into longer momentum vector, hence extremely conned modes of SPPs. On the other hand, the momentum vector can be highly damped, spanning a smaller energy range. For example, when n g ¼ 0.74 Â 10 12 cm À2 and n b ¼ 5.0 Â 10 12 cm À2 or 1.0 Â 10 13 cm À2 , with increasing b, we observe that the upper bright band and the lower bright band span lower frequency ranges ( Fig. 4a and b). In this case the hybrid modes upper branch is graphene like, and lower branch is BP like. For larger b, this gap narrows because of the bending of the upper band as interband transitions comes into play. The maximum hybridization for n g ¼ 0.74 Â 10 12 cm À2 and n b ¼ 1.0 Â 10 13 cm À2 , Fig. 4b, results in a noticeable illustration of the physical origins of higher SHF peaks observed in Fig. 3a. On the other hand, SHF decreases when n b increases to 5.0 Â 10 13 cm À2 . Since the top branch band that is BP like spans most of the frequency range, the lower bright band that is graphene like stays in a small frequency range with decreased b. In graphene, the carrier density is n g ¼ 6.6 Â 10 12 cm À2 , Fig. 5c-e, top branch extends covering most of the frequency range, so interband transitions of graphene are negligible. The lower branch covers a small frequency range with strong connement for n b ¼ 5.0 Â 10 12 cm À2 and 1.0 Â 10 13 cm À2 moving to close to GR branch. The bright modes supported by BP for n b ¼ 5.0 Â 10 13 cm À2 move closer to GR and the heat transmission increases signicantly. Although the bands do not reach a point of overlap, in contrast to lower number density, there is an indication for the origin of maximum NFRHT in Fig. 4c. This also consistent for n g ¼ 2.65 Â 10 13 cm À2 as shown in ESI Fig. 3. † The large range of u covered in the heat transmission factor of the AC direction in BP clearly demonstrate the broad NFRHT as compared to ZZ oriented BP and preferable enhancement in the heterogeneous 2D, this is also further supported in the transmission probability description for highly doped graphene, see ESI Fig. 3. † Fig. 5 Color map of photon transmission probability (x) between suspended BP emitter and graphene sheet receiver (G), and black phosphorene sheet (B) receiver zigzag (ZZ) direction along x-axis orientation separation distance, d ¼ 10 nm, electron density of graphene, n g ¼ 0.74 Â 10 12 cm À2 , and BP, n b , (a) 5.0 Â 10 13 cm À2 , (b) 1.0 Â 10 13 cm À2 , (c) 5.0 Â 10 13 cm À2 , respectively and n g ¼ 6.6 Â 10 12 cm À2 , and BP, n b , (d) 5.0 Â 10 13 cm À2 , (e) 1.0 Â 10 13 cm À2 , (f) 5.0 Â 10 13 cm À2 , respectively.

Distance effect
The strong modulation of the NFRHT between GR and BP is shown to be quite generic, and it appears to be inuenced by a wide range of parameters. We test this concept by exploring the total radiation, Q, computed from eqn (7) as a function of m, n and the gap size, d. We quantify the various parameters signicance, rst by dening the total heat radiation heat transfer ux ratio, h ¼ Q(n g ,n b ) GB /Q(n g ,n g ) GG , as an efficiency factor to gain better insight into the performance of BP and graphene, where Q(n g ,n b ) GB is the BPE-GRR system total heat ux and Q(n g ,n g ) GG is the GRE-GRR system total heat ux.
The total heat ux, Q, is computed for a range of BP and graphene carrier densities for extremely narrow separation d ¼ 10 nm, intermediate d ¼ 100 nm, and a larger separation d ¼ 750 nm, which are displayed in Fig. 6a-g for AC orientation and ratio, ZZ orientation system, and graphene-graphene system. As the vacuum gap distance increases from d ¼ 10 nm to d ¼ 500 and 750 nm, the monotonic decrease of both maximum and minimum total heat ux is obvious because of the weaker SPPs coupling at a larger gap distance. In addition, this reduction is a result of weak near-eld coupling and short propagation distance of the evanescence eld, which decreases the available of plasmonically-activated modes. The heat ux ratio reaches up to the maximum scale 1.5, at d ¼ 10 nm in both AC and ZZ orientation, Fig. 6a and d, respectively. More importantly, h, diminishes for n b of BP and for all n g in AC orientation, while in ZZ orientation of BP leads to a greater h. From this theoretical analysis, we can conclude that BP-GR system shows superiorities compared to GR-GR in designing super Plank radiation heat transfer. As the separation increases further, efficiency of h > 1 is still sustained for the selected doping range, Fig. 6b, c, e and f. For n g from 0.74 Â 10 12 cm À2 to 4.7 Â 10 13 cm À2 and n b from 2.0 Â 10 12 cm À2 to 9.0 Â 10 13 cm À2 , there is a linear bright region that shows efficiency h > 1. Specically, for n g < 1.18 Â 10 13 cm À2 and n b above 2.0 Â 10 12 cm À2 in ZZ orientation and from 0.1 eV to 0.8 eV and BP doping level up to n b ¼ 5.0 Â 10 12 cm À2 AC orientation. Away from this BP and graphene carrier density range, the ratio is less than 1 irrespective of the BP orientation, which is an undesirable case for NFRHT enhancement of the heterogeneous 2D BP-GR system. However, this theoretical result suggests an alternate solution of realizing enhanced and reduced near-eld radiation heat transfer mechanism. In particular, we can freely control the functionality of this heterogeneous system enhancing and reducing heat transfer by changing the orientation of BP and doping of both BP and graphene. Furthermore, our results demonstrate the superiority of this heterogeneous 2D material pair, which give rise to the anisotropic thermal emission by physical reorientation of BP in addition to gate-modulation of graphene and doping of BP.

Conclusions
Our work explores NFRHT originated from SPPs assisted photon transmission for a pair coupled graphene and BP heterogeneous 2D materials. We found that the radiation beyond the Plank limit and enhanced photon tunnelling is favourable in a homogenous graphene-graphene system for low-gate modulation (<3.0 Â 10 12 cm À2 ). However, in the BPE-GRR system, the NFRHT critical value is signicantly enhanced when a gate (>3.0 Â 10 12 cm À2 ) is applied to graphene and BP for selective doping. By choosing appropriate parameters of doping level of BP and graphene, a NFRHT 4 times higher can be achieved in BPE-GRR than GRE-GRR. This enhanced Fig. 6 Shows the heat flux ratio between graphene sheet-BP sheet Q GB and graphene-graphene sheet Q GG for carrier density of graphene 0.74 Â 10 12 cm À2 to 5.0 Â 10 13 cm À2 versus BP carrier density (n b ) 2.0 Â 10 12 cm À2 to 9.0 Â 10 13 cm À2 vacuum gap distance d 10 nm, 100 nm and 750 nm black phosphorene sheet, receiver AC direction along x-axis orientation for (a, b and c) and ZZ direction along x-axis orientation for d, e and f.
NFRHT effect is also realized in the total heat ux BP-GR pairs leading to more than 5 times higher than GR-GR system at d < 100 nm. In addition, by varying the relative orientation of the BP layer, one can switch from amplied NFRHT to reduced or vice versa. The NFRHT in 2D heterogeneous system resulting from preferential orientation would be dependent on the combination of gate-voltage (carrier density). Thus, for a given carrier density of BP and graphene, the heat tunnelling should be adjusted to obtain the desired radiation transfer either enhanced or reduced. We believe that our work contributes to the active pursuit in the eld of SPPs supporting hybrid nonidentical 2D structures both made from atomically thin (few layers) materials. This study can further motivate future studies of hybrid surface plasmon induced NFRHT between new discovered anisotropic 2D ReS 2 and ReSe 2 . 57,58

Conflicts of interest
The authors declare no conict of interest.