Nanoscale heat flux between nanoporous materials

: By combining stochastic electrodynamics and the Maxwell-Garnett description for effective media we study the radiative heat transfer between two nanoporous materials. We show that the heat ﬂux can be sig-niﬁcantly enhanced by air inclusions, which we explain by:(a) the presence of additional surface waves that give rise to supplementary channels for heat transfer throughout the gap, (b) an increase in the contribution given by the ordinary surface waves at resonance, (c) and the appearance of frustrated modes over a broad spectral range. We generalize the known expression for the nanoscale heat ﬂux for anisotropic metamaterials.

With the modern techniques of nanofabrication it is now possible to explore a whole new level of complexity in material science and to fabricate artificial materials that can exhibit a considerable diversity of optical properties [17][18][19][20][21][22]. In many situations, these composite media possess privileged orientations so that their electromagnetic response depends on the direction of photons propagation. When the photon's wavelength in such a medium is large compared to the size of its representative unit cell, the latter behaves effectively like an anisotropic material and therefore may be described by an effective permittivity tensor (and, when necessary, an effective permeability as well). This naturally points to the question of how anisotropy influences the near-field heat transfer.
In this work, we address this question in the particular case of two semi-infinite uniaxial media characterized by optical axes orthogonally oriented with respect to the surface of interaction. The paper is organized as follows: In Sec. II we derive the expression for the heat flux between two anisotropic media. After a brief description of the relevant composite media to our purposes in Sec. III we investigate in Sec. IV the surface and Brewster modes supported by them and their main features. Next, we compare in Sec. V the near-field heat exchanges between two uniaxial media to the classical ones between two isotropic media. Finally, in order to explain the difference in the behavior of isotropic and anisotropic materials, we discuss in Sec. VI the transmission factor in detail between two uniaxial media and in Sec. VII we present our conclusions.

Radiative heat transfer between anisotropic media
Let B 1 and B 2 be two anisotropic semi-infinite bodies, filling respectively the regions z < 0 and z > d and leaving a vacuum gap between them. In order to ensure a stationary process, we assume that the B i are in local thermal equilibrium at a temperature T i , with T 1 = T 2 . The heat flux P between the two bodies is given by where S = E × H is the Poynting vector and A 12 is any surface that separates the two bodies. By taking such a surface to be a plane defined by z = z 0 (0 < z 0 < d) and using the (transverse) translational invariance of our system, the previous equation simplifies to showing that only the z-component of the Poynting vector is needed. After a straightforward calculation, the latter can be conveniently written as [3] where we identify the mean energy of a harmonic oscillator and also the averaged spectral Poynting vector [3] S ω = 2 Re Tr where r = r + zẑ and G(r, r ) is the electrical Green's dyadic, satisfying Moreover, we have introduced Boltzmann's constant k B , Planck's constant 2πh; the † symbolizes hermitian conjugation and Tr the 3 × 3 trace. In order to evaluate the heat flux in the given geometry we have to determine the Green's dyadic inside the gap region. This can be done by considering the multiple scattering of a plane wave due to a source inside the gap [23]. Details and the final expression for the Green's dyadic can be found in appendix A. When inserting the resulting expression in Eq. (40) into the heat flux formula, we find The integral is carried out over all transverse wave vectors κ = (k x , k y ) t including propagating modes with κ < ω/c and evanescent modes with κ > ω/c, where c is the velocity of light in vacuum. The energy transmission coefficient T (ω, κ; d) is different for propagating and evanescent modes and can be stated as where γ r = ω 2 /c 2 − κ 2 and R 1 , R 2 are the 2 × 2 reflection matrices characterizing interfaces. By writing them a bit more explicitly, we see that their elements r λ ,λ i are the reflection coefficients for the scattering of an incoming λ -polarized plane wave into an outgoing λ -polarized wave. Note that expression (8) is very general, as it in principle applies to any crystallographic anisotropy, both electric and magnetic. In the isotropic limit they reduce to the usual Fresnel coefficients and we see that the matrices become diagonal. In addition, we have introduced the matrix D 12 , defined by which gives rise to a Fabry-Pérot-like denominator for T (ω, κ; d) in the isotropic case. From Eqs. (3), (7) and (8) we see that, once the reflection matrices are known, it is possible to determine the heat flux between two arbitrary anisotropic semi-infinite bodies kept at fixed temperatures T 1 and T 2 . Moreover, in order to have an independent check, we verified that Eq. (7) also can be derived from the general scattering formalism derived on Ref.
[24]. In the following we will use these expressions to discuss the heat flux between two uniaxial anisotropic materials with their optical axes normal to the interface.

Porous materials
The structures investigated in this paper are depicted in Fig. 1. They are two semi-infinite media composed by a host isotropic material, defined by its complex dielectric function ε h (ω) = ε h (ω)+iε h (ω) (where ε h (ω) > 0), with uniform cylindrical inclusions oriented in the direction orthogonal to the surface as shown in Fig. 1. These inclusions in turn are filled by a medium of dielectric permittivity ε i , that is also assumed to be isotropic. When the size of the representative unit cell is much smaller than all the other characteristic scales involved, a suitable volume average of the material's local electromagnetic response can be made. In our case, the emerging azimuthal symmetry in this long wavelength limit gives rise to effective uniaxial crystals with a permittivity tensor of the form where e x , e y , and e z are orthogonal unit vectors in x, y, and z direction. The parallel and perpendicular components can be derived from the Maxwell-Garnett effective medium theory (EMT) [25,26] where f is the volume fraction of inclusions. For the structure considered in this work the deviation from the exact result of homogenization given in Refs.
[28, 29] is small even for relatively high filling factors such as f = 0.5. Hence, we will discuss the heat flux between porous media with the Maxwell-Garnett expression for f ∈ [0, 0.5] in this work.
The condition of long wavelengths sets a limit to the lattice constant a of the inclusions for which the EMT can be used. In the far-field regime this condition is fullfilled when the thermal wavelength λ th =hc/k B T is much larger than a. In the near-field region the contributing modes at a distance d above the porous material have a lateral wavelength which depends on d. For κ = 2π/a (which corresponds to a lateral wavelength a) the evanescent waves are damped as in the non-retarded near-field region above the porous material. It follows that the contribution to the heat flux is dominated by evanescent waves with lateral wavelength larger than a if d > a/(2π). On the other hand, one can argue that a nonlocal model for the permittivity is necessary if the lateral wave vectors κ are on the order of π/a. Since the exponential in the transmission coefficient in Eq. (8) for κ > ω/c sets a cutoff for κ of the modes contributing to the near-field heat flux which is ≈ 1/d, one finds that a local EMT description is permissible if d a/π. Hence, for a given lattice constant a of the inclusions the validity of the EMT in Eq. (13) in the near field regime is restricted to d a/π. Artificial structures as depicted in Fig. 1   . Furthermore the dark (blue) areas mark the region for which γ p is purely real, whereas the bright (red) areas are the regions for which γ p is purely imaginary.

Surface and Brewster modes in porous media
Let us study the surface waves supported by these media when they are sufficiently far away from each other so that any coupling of evanescent waves can be neglected. By definition, these surface waves are resonant surface modes and therefore are determined by the poles of the reflection coefficients of these media. For out-of-plane uniaxial media the components of the reflection matrix are where γ s,p are given by the solutions of Fresnel equations in the anisotropic material [30] and hence it follows at once that the surface modes are determined by It is straightforward to verify that in this case only the second equation above can be satisfied, meaning that only p-polarized surface waves can exist at the interface of these media. Solving that equation explicitly for κ gives us the sought dispersion relation of surface waves but one must be aware that Eq. (22) has two branches, and only one is connected to surface modes [31]. Since their dispersion relation involves ε and ε ⊥ , these waves are also called extraordinary surface waves [32], and they reflect the material anisotropy. When ε = ε ⊥ = ε, Eq. (22) degenerates into the well-known dispersion relation κ = ω/c ε/(ε + 1) of surface modes supported by a semi-infinite isotropic medium (bounded by vacuum) with a dielectric permittivity ε. In Fig. 2 where ω L = 1.827 · 10 14 s −1 , ω T = 1.495 · 10 14 s −1 , Γ = 0.9 · 10 12 s −1 , and ε ∞ = 6.7 denote respectively the longitudinal and transversal optical phonon pulsation, the damping factor and the high frequency dielectric constant, respectively. In order to avoid the inherent difficulties of multiple possible interpretations of complex dispersion relations [34], we have deliberately neglected the host material losses to represent these curves. The relevance of this approximation can be checked by comparing Fig. 2 with Fig. 3, where we plot the reflection coefficients of dissipating porous material. In order to distinguish between evanescent and propagative waves inside the effective medium, solutions of Eq. (22) are superimposed in Fig. 2 to a two-color background. This background is a binary representation of ζ = sgn(ε ω 2 /c 2 − ε κ 2 /ε ⊥ ). In the blue zones ζ < 0 so that only evanescent modes can exist, and conversely, in the red zones we have ζ > 0 and all modes are propagative. Similarly the light line ω = cκ allows us to distinguish between the radiative (propagative) and the non-radiative (evanescent) modes inside the vacuum. Notice that, in order to satisfy Eq. (22) both ζ and sgn(ω 2 /c 2 − κ 2 ) must be the same. In other words, frustrated modes cannot satisfy the dispersion relation (22). Now let us turn to the description of modes supported by our artificial structures. For low filling factors we note in Fig. 2(a) the existence of two surface modes. The first one (at a lower frequency) corresponds to the classical surface phonon-polariton (SPP) supported by a massive SiC sample [35]. That surface mode is also present in isotropic SiC. The most interesting feature of Fig. 2(a) is, however, the appearance of a second surface mode at higher frequencies, because it is a signature of the anisotropic character of the material and therefore a direct consequence of the vacuum inclusions in the host medium. As the porosity increases, both surface waves split. Beyond a critical filling factor between f = 0.3 and f = 0.5, the upper surface wave disappears as is seen in Fig. 3. Nevertheless the SPP which still exists continue to move toward the smaller frequencies, i.e., to ω T . Above the light line, we see that the anisotropy gives rise to two different types of Brewster modes. At high frequency we recognize the usual modes where the reflection coefficient [ Fig. 3(a)] of the effective medium vanishes. In addition to these modes, different Brewster modes appear depending on the value of filling factor. Also, we see on the reflection curves (Fig. 3) that the Christiansen point [36] for which the reflectivity is zero for all κ does not depend on the porosity. Indeed, an inspection of expressions (13) and (14) shows that the condition for the Christiansen point of the host material ε h = 1 implies that ε = 1 and ε ⊥ = 1 so that, according to Eq. (17), the reflection coefficients vanish.

Heat flux between porous media
Before we discuss the influence of the inclusions on the heat flux, we show in Fig. 4 the results of the mean Poynting vector S z between two semi-infinite SiC bodies at fixed temperatures T 1 = 300 K and T 2 = 0 K. First of all one can see that the heat flux becomes very large for distances much smaller than the thermal wavelength λ th =hc/k B T (which is about 7.68 μm for T = 300 K). At d = 10 nm the heat flux for the two SiC bodies is about 1000 times larger than the heat flux between two black bodies. This increase is due to the frustrated total internal reflection and to the coupled surface phonon polariton modes [9]. In the propagating regime, i.e., for distances larger than λ th the heat flux is determined by Kirchhoff-Planck's law and is limited by the black-body value. Note, that the heat flux is dominated by the p-polarized modes for distances smaller than 100 nm and larger than 10 μm, whereas for distances in between it is dominated by the s-polarized modes. Now, we introduce the inclusions by using the Maxwell-Garnett expression in Eqs. (13) and (14). We use the same filling factor for both materials, so that we have a symmetric situation. In Fig. 5 we show the resulting heat flux normalized to the values for the two non-pourous SiC plates shown in Fig. 4. We find that for distances smaller than 100 nm and larger than 1 μm the heat flux becomes larger when we add air inclusions, whereas for intermediate distances the heat flux is reduced.
In order to see how the s-and p-mode contribution is changed by the porosity, we show in Fig. 6(a) and 6(b) the plots for the separate contributions of s-and p-polarized modes. It is clear that the p-polarized part of the heat flux gets enhanced for all distances when compared to the isotropic case, regardless of the filling factor. The s-polarized part in turn gives a larger heat flux for distances larger than about 1 μm and a smaller heat flux for distances smaller than 1 μm. Therefore, the smaller heat flux found in Fig. 5 for intermediate distances is associated to the dominance of s-polarized modes in that distance regime.
In summary, by introducing inclusions we find for large and small distances an increase of the heat flux. For the propagating regime (d > λ th ) this can be understood from a simple argument: the vacuum holes simply dilute the material so that, according to Kirchhoff's law, the reflectivity is decreased and hence the emissivity is increased. In fact, for f = 1 one would retrieve the black body result, since in this case the reflectivity is zero. On the other hand, there is no such simple argument for the increased heat flux in the near-field region. Here, it is necessary to study how the coupled surface modes, which give the main contribution to the heat flux for distances smaller than 100 nm, are influenced by the introduction of the inclusions. This will be done by inspection of the transmission coefficient in the next section.

Transmission coefficient
As mentioned before, for the small distance regime (d < 100 nm) the heat flux between two isotropic semi-infinite SiC-bodies is solely dominated by the p-polarized contribution. This remains true for the porous SiC bodies. In fact, the dominance of the p-polarized contribution becomes even greater with increasing filling factors. Hence, to understand the observation that by introducing some porosity the heat flux becomes larger, it suffices to study the p-polarized contribution.
In Fig. 7 we show the transmission coefficient T p (ω, κ; d) in the (ω, κ)-plane for different filling factors and a distance d = 100 nm. In Fig. 7(a) one can see T p (ω, κ; d) for two isotropic SiC plates. Here, T p (ω, κ; d) is one or close to one for the propagating modes, the total internal reflection modes and the coupled surface phonon polaritons. In the plotted region one can mainly see the coupled surface phonon polaritons, which are responsible for the large heat flux at small distances. Now, for f = 0.1 one can see in Fig. 7(b) that a second coupled surface mode appears due to the air inclusions. In addition, the coupled surface mode of the bulk SiC is shifted to smaller frequencies. When increasing the filling factor [ Fig. 7(c) and 7(d)] the upper coupled surface modes shift to higher frequencies and become less important for the transmission coefficient. On the other hand, the low frequency surface modes shift further to lower frequencies.
Between the two coupled surface mode branches a band of frustrated internal reflection modes is formed which gives also a non-negligible contribution to the transmission coefficient.
In order to get further information we now consider the spectral mean Poynting vector S ω defined in Eq. (7) for p-polarization only. We have plotted this quantity in Fig. 8 at the same distance as before, i.e., d = 100 nm, and again for different filling factors. As in Fig. 7 one can see the strong contribution of the two coupled surface mode resonances, which are shifted in frequencies when changing the filling factor. Moreover, the shifting of the primary surface mode to lower frequencies by itself also enhances the flux, as such a shift brings that surface mode closer to the peak wavelength of blackbody radiation as given by the Wien's law. Furthermore, one can now observe, that when increasing the filling factor the low frequency resonance is not only shifted to smaller frequencies, but the resonance is also getting stronger.
The study can now be completed when considering the mean transmission factor for the p-polarized modes, that was introduced in Ref. [9] as with u =hω/k B T and f = u 2 e u /(e u − 1) 2 . It represents the mean transmission coefficient of a mode specified by it's wave vector κ for a given temperature T and a small temperature difference ΔT between the two bodies. By means of this quantity the heat flux can be rewritten in a Landauer-like form [9] S z = π 2 3 Note, that for κd 1 and κ > ω/c the transmission coefficient T p (ω, κ; d) is exponentially damped [see Eq. (8)] and therefore also the mean transmission factor T p (κ). This damping determines the wave vector cutoff and hence the number of states contributing to the heat flux. Now, in Fig. 9 we plot T p (κ) for a given distance of d = 100 nm and different filling factors normalized to the mean transmission factor for two semi-infinite SiC bodies. For f = 0.1 the mean transmission coefficient for the porous SiC increases for intermediate κ but decreases for very large κ. The increased mean transmission factor is due to the second coupled surface mode and the frustrated modes, whereas the lower value for large wave vectors can be attributed to a stronger cutoff in the transmission coefficient, which means that the number of contributing modes is decreased. The enhancement of the transmission factor due to the surface mode prevails and leads to an enhanced heat flux at that distance. The same mechanism is responsible for  the enhanced heat flux for f = 0.3. On the other hand, for larger filling factors the curves change slightly for intermediate κ compared to the curve for f = 0.3. The contribution in that intermediate region is due to the second coupled surface mode branch and the frustrated modes. But for very large κ the mean transmission coefficient increases compared to f = 0.3. This means that by introducing a higher porosity we soften the cutoff of the transmission coefficient. Hence, the number of modes contributing to the heat flux is increased and results for large filling factors in a further enhanced heat flux.
The dependence of the cutoff on the filling factor for large κ can easily be discussed for the where the permittivities have to be evaluated at the surface mode resonance frequency of the semi-infinite anisotropic body (see Appendix B). In Fig. 10 we show a plot of κ uni /κ iso over the filling factor. It is seen that by introducing the air inclusions the cutoff first decreases and then monotonically increases. This is the same qualitative behavior as observed for the mean transmission factor T p (κ) in Fig. 9 for κd 1. This reasoning confirms that the number of contributing modes is the main mechanism for increasing the heat flux at small distances and large filling factors ( f > 0.3).

Conclusion
We have presented a detailed study of near and far field heat transfer between two flat uniaxial media made of polar materials (in our case, SiC) in which cylindrical inclusions drilled orthogonally to surfaces are uniformally distributed.
After applying the classical stochastic electrodynamic theory to anisotropic materials we have shown that, for short distances, the heat flux between such media can be significantly larger than those traditionally measured between two isotropic materials in the same non-equilibrium thermal conditions. For small filling factors we have determined that this enhancement stems from additional surface waves arising at the uniaxial material-vacuum interface, clearly indicating that such increase is intrinsically connected to anisotropy. Indeed, we did calculations for isotropically rarified SiC plates with low filling factors ( f ≤ 0.1) and found that the heat transfer modification for is much smaller. In contrast, for larger filling factors ( f > 0.3) we have shown that, after a thorough analysis of the transmission factor, the enhancement in heat transfer arises mainly from the increased number of modes contributing to the flux.